This page presents a variety of calculations for latitude/longitude points, with the formulæ and code fragments for implementing them.
All these formulæ are for calculations on the basis of a spherical earth (ignoring ellipsoidal effects) – which is accurate enough* for most purposes… [In fact, the earth is very slightly ellipsoidal; using a spherical model gives errors typically up to 0.3% – see notes for further details].
Enter the co-ordinates into the text boxes to try out the calculations. A variety of formats are accepted, principally:
deg-min-sec suffixed with N/S/E/W (e.g. 40°44′55″N, 73 59 11W), or
signed decimal degrees without compass direction, where negative indicates west/south (e.g. 40.7486, -73.9864):
Point 1: ,
Point 2: ,
Distance: 170.3 km
Initial bearing: 127°12′30″
Final bearing: 128°47′27″
Midpoint: 52°40′48″N, 000°50′37″W
And you can see it on a map (aren’t those Google guys wonderful!)
Distance
This uses the ‘Haversine’ formula to calculate great-circle distances between the two points – that is, the shortest distance over the earth’s surface – giving an ‘as-the-crow-flies’ distance between the points (ignoring any hills!).
Haversine formula:
R = earth’s radius (mean radius = 6,371km)
Δlat = lat2− lat1
Δlong = long2− long1
a = sin²(Δlat/2) + cos(lat1).cos(lat2).sin²(Δlong/2)
c = 2.atan2(√a, √(1−a))
d = R.c
(Note that angles need to be in radians to pass to trig functions).
JavaScript:
var R = 6371; // km
var dLat = (lat2-lat1).toRad();
var dLon = (lon2-lon1).toRad();
var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.cos(lat1.toRad()) * Math.cos(lat2.toRad()) *
Math.sin(dLon/2) * Math.sin(dLon/2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d = R * c;
The Haversine formula ‘remains particularly well-conditioned for numerical computation even at small distances’ – unlike calculations based on the spherical law of cosines. (It was published by R W Sinnott in Sky and Telescope, 1984, though has been known about for much longer; the ‘half-versed-sine’ is (1-cosθ)/2, or sin²(θ/2) – don’t ask, I’m not a mathematician).
Spherical Law of Cosines
In fact, when Sinnott published the Haversine formula, computational precision was limited. Nowadays, JavaScript (and most modern computers & languages) use IEEE 754 64-bit floating-point numbers, which provide 15 significant figures of precision. With this precision, the simple spherical law of cosines formula gives well-conditioned results down to distances as small as around 1 metre. In view of this it is probably worth, in most situations, using either the simpler law of cosines or the more accurate ellipsoidal Vincenty formula in preference to Haversine! (bearing in mind notes below on the limitations in accuracy of the spherical model).
Spherical law
of cosines: d = acos(sin(lat1).sin(lat2)+cos(lat1).cos(lat2).cos(long2−long1)).R
JavaScript:
var R = 6371; // km
var d = Math.acos(Math.sin(lat1)*Math.sin(lat2) +
Math.cos(lat1)*Math.cos(lat2) *
Math.cos(lon2-lon1)) * R;
Excel: =ACOS(SIN(lat1)*SIN(lat2)+COS(lat1)*COS(lat2)*COS(lon2-lon1))*6371
(Note that here and in all subsequent code fragments, for simplicity I do not show conversions from degrees to radians; see below for complete versions).
Bearing
In general, your current bearing will vary as you follow a great circle path (orthodrome); the final bearing will differ from the initial bearing by varying degrees according to distance and latitude (if you were to go from say 35°N,45°E (Baghdad) to 35°N,135°E (Osaka), you would start on a bearing of 60° and end up on a bearing of 120°!).
This is the initial bearing (sometimes referred to as forward azimuth) which if followed in a straight line along a great-circle arc will take you from the start point to the end point:
Formula: θ = atan2( sin(Δlong).cos(lat2),
cos(lat1).sin(lat2) − sin(lat1).cos(lat2).cos(Δlong) )
JavaScript:
var y = Math.sin(dLon) * Math.cos(lat2);
var x = Math.cos(lat1)*Math.sin(lat2) -
Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
var brng = Math.atan2(y, x).toDeg();
Excel: =ATAN2(COS(lat1)*SIN(lat2)-SIN(lat1)*COS(lat2)*COS(lon2-lon1),
SIN(lon2-lon1)*COS(lat2))
* Note that Excel reverses the arguments to ATAN2 – see notes below
Since atan2 returns values in the range -π ... +π (that is, -180° ... +180°), to normalise the result to a compass bearing (in the range 0° ... 360°, with -ve values transformed into the range 180° ... 360°), convert to degrees and then use (θ+360) % 360, where % is modulo.
For final bearing, simply take the initial bearing from the end point to the start point and reverse it (using θ = (θ+180) % 360).
Midpoint
This is the midpoint along a great circle path between the two points.
Formula: Bx = cos(lat2).cos(Δlong)
By = cos(lat2).sin(Δlong)
latm = atan2(sin(lat1) + sin(lat2), √((cos(lat1)+Bx)² + By²))
lonm = lon1 + atan2(By, cos(lat1)+Bx)
JavaScript:
var Bx = Math.cos(lat2) * Math.cos(dLon);
var By = Math.cos(lat2) * Math.sin(dLon);
var lat3 = Math.atan2(Math.sin(lat1)+Math.sin(lat2),
Math.sqrt( (Math.cos(lat1)+Bx)*(Math.cos(lat1)+Bx) + By*By) );
var lon3 = lon1 + Math.atan2(By, Math.cos(lat1) + Bx);
Just as the initial bearing may vary from the final bearing, the midpoint may not be located half-way between latitudes/longitudes; the midpoint between 35°N,45°E and 35°N,135°E is around 45°N,90°E.
Destination point given distance and bearing from start point
Given a start point, initial bearing, and distance, this will calculate the destination point and final bearing travelling along a (shortest distance) great circle arc.
Start point: ,
Bearing:
Distance: km
Destination point: 53°11′18″N, 000°08′00″E
Final bearing: 097°30′52″
view map
Formula: lat2 = asin(sin(lat1)*cos(d/R) + cos(lat1)*sin(d/R)*cos(θ))
lon2 = lon1 + atan2(sin(θ)*sin(d/R)*cos(lat1), cos(d/R)−sin(lat1)*sin(lat2))
d/R is the angular distance (in radians), where d is the distance travelled and R is the earth’s radius
JavaScript:
var lat2 = Math.asin( Math.sin(lat1)*Math.cos(d/R) +
Math.cos(lat1)*Math.sin(d/R)*Math.cos(brng) );
var lon2 = lon1 + Math.atan2(Math.sin(brng)*Math.sin(d/R)*Math.cos(lat1),
Math.cos(d/R)-Math.sin(lat1)*Math.sin(lat2));
Excel: lat2: =ASIN(SIN(lat1)*COS(d/R) + COS(lat1)*SIN(d/R)*COS(brng))
lon2: =lon1 + ATAN2(COS(d/R)-SIN(lat1)*SIN(lat2), SIN(brng)*SIN(d/R)*COS(lat1))
For final bearing, simply take the initial bearing from the end point to the start point and reverse it (using θ = (θ+180) % 360).
Intersection of two paths given start points and bearings
This is a rather more complex calculation than most others on this page, but I've been asked for it a number of times. See below for the JavaScript.
Point 1: , Brng 1:
Point 2: , Brng 2:
Intersection point: 50°54′06″N, 004°29′39″E
Formula:
d12 = 2.asin( √(sin²(Δlat/2) + cos(lat1).cos(lat2).sin²(Δlon/2)) )
φ1 = acos( sin(lat2) − sin(lat1).cos(d12) / sin(d12).cos(lat1) )
φ2 = acos( sin(lat1) − sin(lat2).cos(d12) / sin(d12).cos(lat2) )
if sin(lon2−lon1) > 0
θ12 = φ1, θ21 = 2.π − φ2
else
θ12 = 2.π − φ1, θ21 = φ2
α1 = (θ1 − θ12 + π) % 2.π − π
α2 = (θ21 − θ2 + π) % 2.π − π
α3 = acos( −cos(α1).cos(α2) + sin(α1).sin(α2).cos(d12) )
d13 = atan2( sin(d12).sin(α1).sin(α2), cos(α2)+cos(α1).cos(α3) )
lat3 = asin( sin(lat1).cos(d13) + cos(lat1).sin(d13).cos(θ1) )
Δlon = atan2( sin(θ1).sin(d13).cos(lat1), cos(d13)−sin(lat1).sin(lat3) )
lon3 = (lon1−Δlon+π) % 2.π − π
where
lat1, lon1, θ1 : 1st point & bearing
lat2, lon2, θ2 : 2nd point & bearing
lat3, lon3 : intersection point
% = mod
note – if sin(α1)=0 and sin(α2)=0: infinite solutions
if sin(α1).sin(α2) < 0: ambiguous solution
this formulation is not always well-conditioned for meridional or equatorial lines
Note this can also be solved using vectors rather than trigonometry:
For each point φ,λ (lat=φ, lon=λ), we can define a unit vector pointing to it from the centre of the earth: u{x,y,z} = [ cosφ·cosλ, cosφ·sinλ, sinφ ] (taking x=0º, y=90º, z=north – note that these formulæ depend on convention used for directions and handedness)
And for any great circle defined by two points, we can define a unit vector N normal to the plane of the circle: N(u1, u2) = (u1×u2) / ||u1×u2|| where × is the vector cross product, and ||u|| the norm (length of the vector)
The vector representing the intersection of the two great circles is then ui = ±N( N(u1, u2), N(u3, u4) )
We can then get the latitude and longitude of Pi by φ = atan2(uz, sqrt(ux² + uy²)), λ = atan2(uy, ux)
The antipodal intersection point is (-φ, λ+π)
Cross-track distance
Here’s a new one: I’ve sometimes been asked about distance of a point from a great-circle path (sometimes called cross track error).
Formula: dxt = asin(sin(d13/R)*sin(θ13−θ12)) * R
where d13 is distance from start point to third point
θ13 is (initial) bearing from start point to third point
θ12 is (initial) bearing from start point to end point
R is the earth’s radius
JavaScript:
var dXt = Math.asin(Math.sin(d13/R)*Math.sin(brng13-brng12)) * R;
Here, the great-circle path is identified by a start point and an end point – depending on what initial data you’re working from, you can use the formulæ above to obtain the relevant distance and bearings. The sign of dxt tells you which side of the path the third point is on.
The along-track distance, from the start point to the closest point on the path to the third point, is
Formula: dat = acos(cos(d13/R)/cos(dxt/R)) * R
where d13 is distance from start point to third point
dxt is cross-track distance
R is the earth’s radius
JavaScript:
var dAt = Math.acos(Math.cos(d13/R)/Math.cos(dXt/R)) * R;
Closest point to the poles
And: ‘Clairaut’s formula’ will give you the maximum latitude of a great circle path, given a bearing and latitude on the great circle:
Formula: latmax = acos(abs(sin(θ)*cos(lat)))
JavaScript:
var latMax = Math.acos(Math.abs(Math.sin(brng)*Math.cos(lat)));
Rhumb lines
A ‘rhumb line’ (or loxodrome) is a path of constant bearing, which crosses all meridians at the same angle.
Sailors used to (and sometimes still) navigate along rhumb lines since it is easier to follow a constant compass bearing than to be continually adjusting the bearing, as is needed to follow a great circle. Rhumb lines are straight lines on a Mercator Projection map (also helpful for navigation).
Rhumb lines are generally longer than great-circle (orthodrome) routes. For instance, London to New York is 4% longer along a rhumb line than along a great circle – important for aviation fuel, but not particularly to sailing vessels. New York to Beijing – close to the most extreme example possible (though not sailable!) – is 30% longer along a rhumb line.
Point 1: ,
Point 2: ,
Distance: 5196 km
Bearing: 260°07′38″
view map
Start point: ,
Bearing:
Distance: km
Destination point: 50°57′48″N, 001°51′09″E
view map
Distance/bearing
These formulæ give the distance and (constant) bearing between two points.
Formula: Δφ = ln(tan(lat2/2+π/4)/tan(lat1/2+π/4)) [= the ‘stretched’ latitude difference]
if E:W line, q = cos(lat1)
otherwise, q = Δlat/Δφ
d = √(Δlat² + q².Δlon²).R [pythagoras]
θ = atan2(Δlon, Δφ)
where ln is natural log, Δlon is taking shortest route (<180º), and R is the earth’s radius
JavaScript:
var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
var q = (!isNaN(dLat/dPhi)) ? dLat/dPhi : Math.cos(lat1); // E-W line gives dPhi=0
// if dLon over 180° take shorter rhumb across 180° meridian:
if (Math.abs(dLon) > Math.PI) {
dLon = dLon>0 ? -(2*Math.PI-dLon) : (2*Math.PI+dLon);
}
var d = Math.sqrt(dLat*dLat + q*q*dLon*dLon) * R;
var brng = Math.atan2(dLon, dPhi);
Destination
Given a start point and a distance d along constant bearing θ, this will calculate the destination point. If you maintain a constant bearing along a rhumb line, you will gradually spiral in towards one of the poles.
Formula: α = d/R (angular distance)
lat2 = lat1 + α.cos(θ)
Δφ = ln(tan(lat2/2+π/4)/tan(lat1/2+π/4)) [= the ‘stretched’ latitude difference]
if E:W line q = cos(lat1)
otherwise q = Δlat/Δφ
Δlon = α.sin(θ)/q
lon2 = (lon1+Δlon+π) % 2.π − π
where ln is natural log and % is modulo, Δlon is taking shortest route (<180°), and R is the earth’s radius
JavaScript:
lat2 = lat1 + d*Math.cos(brng);
var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
var q = (!isNaN(dLat/dPhi)) ? dLat/dPhi : Math.cos(lat1); // E-W line gives dPhi=0
var dLon = d*Math.sin(brng)/q;
// check for some daft bugger going past the pole, normalise latitude if so
if (Math.abs(lat2) > Math.PI/2) lat2 = lat2>0 ? Math.PI-lat2 : -(Math.PI-lat2);
lon2 = (lon1+dLon+Math.PI)%(2*Math.PI) - Math.PI;
If you use Ordnance Survey Grid References, I have implemented a script for converting between Lat/Long & OS Grid References.
Convert between degrees-minutes-seconds & decimal degrees
Latitude Longitude 1° ≈ 111 km (110.57 eq’l — 111.70 polar)
1′ ≈ 1.85 km (= 1 nm) 0.01° ≈ 1.11 km
1″ ≈ 30.9 m 0.0001° ≈ 11.1 m
No, I’ve not included decimal minutes: a decimal system is easy, a sexagesimal system has merits, but mixing the two is a complete sow’s ear. Switch off the option on your GPS!
Notes:
Accuracy: since the earth is not quite a sphere, there are small errors in using spherical geometry; the earth is actually roughly ellipsoidal (or more precisely, oblate spheroidal) with a radius varying between about 6,378km (equatorial) and 6,357km (polar), and local radius of curvature varying from 6,336km (equatorial meridian) to 6,399km (polar). 6,371 km is the generally accepted value for the Earth’s mean radius. This means that errors from assuming spherical geometry might be up to 0.55% crossing the equator, though generally below 0.3%, depending on latitude and direction of travel. An accuracy of better than 3m in 1km is mostly good enough for me, but if you want greater accuracy, you could use the Vincenty formula for calculating geodesic distances on ellipsoids, which gives results accurate to within 1mm. (Out of sheer perversity – I’ve never needed such accuracy – I looked up this formula and discovered the JavaScript implementation was simpler than I expected).
Trig functions take arguments in radians, so latitude, longitude, and bearings in degrees (either decimal or degrees/minutes/seconds) need to be converted to radians, rad = π.deg/180. When converting radians back to degrees (deg = 180.rad/π), West is negative if using signed decimal degrees. For bearings, values in the range -π to +π [-180° to +180°] need to be converted to 0 to +2π [0°–360°]; this can be done by (brng+2.π)%2.π [or brng+360)%360] where % is the modulo operator.
The atan2() function widely used here takes two arguments, atan2(y, x), and computes the arc tangent of the ratio y/x. It is more flexible than atan(y/x), since it handles x=0, and it also returns values in all 4 quadrants -π to +π (the atan function returns values in the range -π/2 to +π/2).
If you implement any formula involving atan2 in Microsoft Excel, you will need to reverse the arguments, as Excel has them the opposite way around from JavaScript – conventional order is atan2(y, x), but Excel uses atan2(x, y). To use atan2 in a (VBA) macro, you can use WorksheetFunction.Atan2().
All bearings are with respect to true north, 0°=N, 90°=E, etc; if you are working from a compass, magnetic north varies from true north in a complex way around the earth, and the difference has to be compensated for by variances indicated on local maps.
I learned a lot from the US Census Bureau GIS FAQ which is no longer available, so I’ve made a copy.
Thanks to Ed Williams’ Aviation Formulary for many of the formulæ.
For miles, divide km by 1.609344
For nautical miles, divide km by 1.852
See below for the source code of the JavaScript implementation. These functions should be simple to translate into other languages if required.
Update January 2010: I have revised the scripts to be structured as methods of a LatLon object. Of course, JavaScript is a prototype-based rather than class-based language, so this is only nominally a class, but isolating code into a separate namespace is good JavaScript practice, and this approach may also make it clearer to implement these functions in other languages. If you’re not familiar with JavaScript syntax, LatLon.prototype.distanceTo = function(point) { ... }, for instance, defines a ‘distanceTo’ method of the LatLon object (/class) which takes a LatLon object as a parameter (and returns a number). The Geo namespace acts as a static class for geodesy formatting / parsing / conversion functions. I have extended (polluted, if you like) the base JavaScript object prototypes with trim(), toRad() toDeg(), and toPrecisionFixed() methods. I’ve adopted JSDoc format for the descriptions.
I have also created a page illustrating the use of the spherical law of cosines for selecting points from a database within a specified bounding circle – the example is based on MySQL+PDO, but should be extensible to other DBMS platforms.
Several people have asked about example Excel spreadsheets, so I have implemented the distance & bearing and the destination point formulæ as spreadsheets, in a form which breaks down the all stages involved to illustrate the operation.
I offer these formulæ & scripts for free use and adaptation as my contribution to the open-source info-sphere from which I have received so much. You are welcome to re-use these scripts [under a simple attribution license, without any warranty express or implied] provided solely that you retain my copyright notice and a link to this page.
If you would like to show your appreciation and support continued development of these scripts, I would most gratefully accept donations.
If you need any advice or development work done, I am available for consultancy.
If you have any queries or find any problems, contact me at ku.oc.epyt-elbavom@oeg-stpircs.
© 2002-2010 Chris Veness
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Latitude/longitude spherical geodesy formulae & scripts (c) Chris Veness 2002-2010 */
/* - www.movable-type.co.uk/scripts/latlong.html */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Note that minimal error checking is performed in this example code! */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* Creates a point on the earth's surface at the supplied latitude / longitude
*
* @constructor
* @param {Number} lat: latitude in numeric degrees
* @param {Number} lon: longitude in numeric degrees
* @param {Number} [rad=6371]: radius of earth if different value is required from standard 6,371km
*/
function LatLon(lat, lon, rad) {
if (typeof rad == 'undefined') rad = 6371; // earth's mean radius in km
this._lat = lat;
this._lon = lon;
this._radius = rad;
}
/**
* Returns the distance from this point to the supplied point, in km
* (using Haversine formula)
*
* from: Haversine formula - R. W. Sinnott, "Virtues of the Haversine",
* Sky and Telescope, vol 68, no 2, 1984
*
* @param {LatLon} point: Latitude/longitude of destination point
* @param {Number} [precision=4]: no of significant digits to use for returned value
* @returns {Number} Distance in km between this point and destination point
*/
LatLon.prototype.distanceTo = function(point, precision) {
// default 4 sig figs reflects typical 0.3% accuracy of spherical model
if (typeof precision == 'undefined') precision = 4;
var R = this._radius;
var lat1 = this._lat.toRad(), lon1 = this._lon.toRad();
var lat2 = point._lat.toRad(), lon2 = point._lon.toRad();
var dLat = lat2 - lat1;
var dLon = lon2 - lon1;
var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.cos(lat1) * Math.cos(lat2) *
Math.sin(dLon/2) * Math.sin(dLon/2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d = R * c;
return d.toPrecisionFixed(precision);
}
/**
* Returns the (initial) bearing from this point to the supplied point, in degrees
* see http://williams.best.vwh.net/avform.htm#Crs
*
* @param {LatLon} point: Latitude/longitude of destination point
* @returns {Number} Initial bearing in degrees from North
*/
LatLon.prototype.bearingTo = function(point) {
var lat1 = this._lat.toRad(), lat2 = point._lat.toRad();
var dLon = (point._lon-this._lon).toRad();
var y = Math.sin(dLon) * Math.cos(lat2);
var x = Math.cos(lat1)*Math.sin(lat2) -
Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
var brng = Math.atan2(y, x);
return (brng.toDeg()+360) % 360;
}
/**
* Returns final bearing arriving at supplied destination point from this point; the final bearing
* will differ from the initial bearing by varying degrees according to distance and latitude
*
* @param {LatLon} point: Latitude/longitude of destination point
* @returns {Number} Final bearing in degrees from North
*/
LatLon.prototype.finalBearingTo = function(point) {
// get initial bearing from supplied point back to this point...
var lat1 = point._lat.toRad(), lat2 = this._lat.toRad();
var dLon = (this._lon-point._lon).toRad();
var y = Math.sin(dLon) * Math.cos(lat2);
var x = Math.cos(lat1)*Math.sin(lat2) -
Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
var brng = Math.atan2(y, x);
// ... & reverse it by adding 180°
return (brng.toDeg()+180) % 360;
}
/**
* Returns the midpoint between this point and the supplied point.
* see http://mathforum.org/library/drmath/view/51822.html for derivation
*
* @param {LatLon} point: Latitude/longitude of destination point
* @returns {LatLon} Midpoint between this point and the supplied point
*/
LatLon.prototype.midpointTo = function(point) {
lat1 = this._lat.toRad(), lon1 = this._lon.toRad();
lat2 = point._lat.toRad();
var dLon = (point._lon-this._lon).toRad();
var Bx = Math.cos(lat2) * Math.cos(dLon);
var By = Math.cos(lat2) * Math.sin(dLon);
lat3 = Math.atan2(Math.sin(lat1)+Math.sin(lat2),
Math.sqrt( (Math.cos(lat1)+Bx)*(Math.cos(lat1)+Bx) + By*By) );
lon3 = lon1 + Math.atan2(By, Math.cos(lat1) + Bx);
return new LatLon(lat3.toDeg(), lon3.toDeg());
}
/**
* Returns the destination point from this point having travelled the given distance (in km) on the
* given initial bearing (bearing may vary before destination is reached)
*
* see http://williams.best.vwh.net/avform.htm#LL
*
* @param {Number} brng: Initial bearing in degrees
* @param {Number} dist: Distance in km
* @returns {LatLon} Destination point
*/
LatLon.prototype.destinationPoint = function(brng, dist) {
dist = dist/this._radius; // convert dist to angular distance in radians
brng = brng.toRad(); //
var lat1 = this._lat.toRad(), lon1 = this._lon.toRad();
var lat2 = Math.asin( Math.sin(lat1)*Math.cos(dist) +
Math.cos(lat1)*Math.sin(dist)*Math.cos(brng) );
var lon2 = lon1 + Math.atan2(Math.sin(brng)*Math.sin(dist)*Math.cos(lat1),
Math.cos(dist)-Math.sin(lat1)*Math.sin(lat2));
lon2 = (lon2+3*Math.PI)%(2*Math.PI) - Math.PI; // normalise to -180...+180
if (isNaN(lat2) || isNaN(lon2)) return null;
return new LatLon(lat2.toDeg(), lon2.toDeg());
}
/**
* Returns the point of intersection of two paths defined by point and bearing
*
* see http://williams.best.vwh.net/avform.htm#Intersection
*
* @param {LatLon} p1: First point
* @param {Number} brng1: Initial bearing from first point
* @param {LatLon} p2: Second point
* @param {Number} brng2: Initial bearing from second point
* @returns {LatLon} Destination point (null if no unique intersection defined)
*/
LatLon.intersection = function(p1, brng1, p2, brng2) {
lat1 = p1._lat.toRad(), lon1 = p1._lon.toRad();
lat2 = p2._lat.toRad(), lon2 = p2._lon.toRad();
brng13 = brng1.toRad(), brng23 = brng2.toRad();
dLat = lat2-lat1, dLon = lon2-lon1;
dist12 = 2*Math.asin( Math.sqrt( Math.sin(dLat/2)*Math.sin(dLat/2) +
Math.cos(lat1)*Math.cos(lat2)*Math.sin(dLon/2)*Math.sin(dLon/2) ) );
if (dist12 == 0) return null;
// initial/final bearings between points
brngA = Math.acos( ( Math.sin(lat2) - Math.sin(lat1)*Math.cos(dist12) ) /
( Math.sin(dist12)*Math.cos(lat1) ) );
if (isNaN(brngA)) brngA = 0; // protect against rounding
brngB = Math.acos( ( Math.sin(lat1) - Math.sin(lat2)*Math.cos(dist12) ) /
( Math.sin(dist12)*Math.cos(lat2) ) );
if (Math.sin(lon2-lon1) > 0) {
brng12 = brngA;
brng21 = 2*Math.PI - brngB;
} else {
brng12 = 2*Math.PI - brngA;
brng21 = brngB;
}
alpha1 = (brng13 - brng12 + Math.PI) % (2*Math.PI) - Math.PI; // angle 2-1-3
alpha2 = (brng21 - brng23 + Math.PI) % (2*Math.PI) - Math.PI; // angle 1-2-3
if (Math.sin(alpha1)==0 && Math.sin(alpha2)==0) return null; // infinite intersections
if (Math.sin(alpha1)*Math.sin(alpha2) < 0) return null; // ambiguous intersection
//alpha1 = Math.abs(alpha1);
//alpha2 = Math.abs(alpha2);
// ... Ed Williams takes abs of alpha1/alpha2, but seems to break calculation?
alpha3 = Math.acos( -Math.cos(alpha1)*Math.cos(alpha2) +
Math.sin(alpha1)*Math.sin(alpha2)*Math.cos(dist12) );
dist13 = Math.atan2( Math.sin(dist12)*Math.sin(alpha1)*Math.sin(alpha2),
Math.cos(alpha2)+Math.cos(alpha1)*Math.cos(alpha3) )
lat3 = Math.asin( Math.sin(lat1)*Math.cos(dist13) +
Math.cos(lat1)*Math.sin(dist13)*Math.cos(brng13) );
dLon13 = Math.atan2( Math.sin(brng13)*Math.sin(dist13)*Math.cos(lat1),
Math.cos(dist13)-Math.sin(lat1)*Math.sin(lat3) );
lon3 = lon1+dLon13;
lon3 = (lon3+Math.PI) % (2*Math.PI) - Math.PI; // normalise to -180..180º
return new LatLon(lat3.toDeg(), lon3.toDeg());
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* Returns the distance from this point to the supplied point, in km, travelling along a rhumb line
*
* see http://williams.best.vwh.net/avform.htm#Rhumb
*
* @param {LatLon} point: Latitude/longitude of destination point
* @returns {Number} Distance in km between this point and destination point
*/
LatLon.prototype.rhumbDistanceTo = function(point) {
var R = this._radius;
var lat1 = this._lat.toRad(), lat2 = point._lat.toRad();
var dLat = (point._lat-this._lat).toRad();
var dLon = Math.abs(point._lon-this._lon).toRad();
var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
var q = (!isNaN(dLat/dPhi)) ? dLat/dPhi : Math.cos(lat1); // E-W line gives dPhi=0
// if dLon over 180° take shorter rhumb across 180° meridian:
if (dLon > Math.PI) dLon = 2*Math.PI - dLon;
var dist = Math.sqrt(dLat*dLat + q*q*dLon*dLon) * R;
return dist.toPrecisionFixed(4); // 4 sig figs reflects typical 0.3% accuracy of spherical model
}
/**
* Returns the bearing from this point to the supplied point along a rhumb line, in degrees
*
* @param {LatLon} point: Latitude/longitude of destination point
* @returns {Number} Bearing in degrees from North
*/
LatLon.prototype.rhumbBearingTo = function(point) {
var lat1 = this._lat.toRad(), lat2 = point._lat.toRad();
var dLon = (point._lon-this._lon).toRad();
var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
if (Math.abs(dLon) > Math.PI) dLon = dLon>0 ? -(2*Math.PI-dLon) : (2*Math.PI+dLon);
var brng = Math.atan2(dLon, dPhi);
return (brng.toDeg()+360) % 360;
}
/**
* Returns the destination point from this point having travelled the given distance (in km) on the
* given bearing along a rhumb line
*
* @param {Number} brng: Bearing in degrees from North
* @param {Number} dist: Distance in km
* @returns {LatLon} Destination point
*/
LatLon.prototype.rhumbDestinationPoint = function(brng, dist) {
var R = this._radius;
var d = parseFloat(dist)/R; // d = angular distance covered on earth's surface
var lat1 = this._lat.toRad(), lon1 = this._lon.toRad();
brng = brng.toRad();
var lat2 = lat1 + d*Math.cos(brng);
var dLat = lat2-lat1;
var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
var q = (!isNaN(dLat/dPhi)) ? dLat/dPhi : Math.cos(lat1); // E-W line gives dPhi=0
var dLon = d*Math.sin(brng)/q;
// check for some daft bugger going past the pole
if (Math.abs(lat2) > Math.PI/2) lat2 = lat2>0 ? Math.PI-lat2 : -(Math.PI-lat2);
lon2 = (lon1+dLon+3*Math.PI)%(2*Math.PI) - Math.PI;
return new LatLon(lat2.toDeg(), lon2.toDeg());
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* Returns the latitude of this point; signed numeric degrees if no format, otherwise format & dp
* as per Geo.toLat()
*
* @param {String} [format]: Return value as 'd', 'dm', 'dms'
* @param {Number} [dp=0|2|4]: No of decimal places to display
* @returns {Number|String} Numeric degrees if no format specified, otherwise deg/min/sec
*
* @requires Geo
*/
LatLon.prototype.lat = function(format, dp) {
if (typeof format == 'undefined') return this._lat;
return Geo.toLat(this._lat, format, dp);
}
/**
* Returns the longitude of this point; signed numeric degrees if no format, otherwise format & dp
* as per Geo.toLon()
*
* @param {String} [format]: Return value as 'd', 'dm', 'dms'
* @param {Number} [dp=0|2|4]: No of decimal places to display
* @returns {Number|String} Numeric degrees if no format specified, otherwise deg/min/sec
*
* @requires Geo
*/
LatLon.prototype.lon = function(format, dp) {
if (typeof format == 'undefined') return this._lon;
return Geo.toLon(this._lon, format, dp);
}
/**
* Returns a string representation of this point; format and dp as per lat()/lon()
*
* @param {String} [format]: Return value as 'd', 'dm', 'dms'
* @param {Number} [dp=0|2|4]: No of decimal places to display
* @returns {String} Comma-separated latitude/longitude
*
* @requires Geo
*/
LatLon.prototype.toString = function(format, dp) {
if (typeof format == 'undefined') format = 'dms';
return Geo.toLat(this._lat, format, dp) + ', ' + Geo.toLon(this._lon, format, dp);
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
// extend Number object with methods for converting degrees/radians
/** Convert numeric degrees to radians */
if (typeof(String.prototype.toRad) === "undefined") {
Number.prototype.toRad = function() {
return this * Math.PI / 180;
}
}
/** Convert radians to numeric (signed) degrees */
if (typeof(String.prototype.toDeg) === "undefined") {
Number.prototype.toDeg = function() {
return this * 180 / Math.PI;
}
}
/**
* Format the significant digits of a number, using only fixed-point notation (no exponential)
*
* @param {Number} precision: Number of significant digits to appear in the returned string
* @returns {String} A string representation of number which contains precision significant digits
*/
if (typeof(String.prototype.toDeg) === "undefined") {
Number.prototype.toPrecisionFixed = function(precision) {
var numb = this < 0 ? -this : this; // can't take log of -ve number...
var sign = this < 0 ? '-' : '';
if (numb == 0) { n = '0.'; while (precision--) n += '0'; return n }; // can't take log of zero
var scale = Math.ceil(Math.log(numb)*Math.LOG10E); // no of digits before decimal
var n = String(Math.round(numb * Math.pow(10, precision-scale)));
if (scale > 0) { // add trailing zeros & insert decimal as required
l = scale - n.length;
while (l-- > 0) n = n + '0';
if (scale < n.length) n = n.slice(0,scale) + '.' + n.slice(scale);
} else { // prefix decimal and leading zeros if required
while (scale++ < 0) n = '0' + n;
n = '0.' + n;
}
return sign + n;
}
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
All these formulæ are for calculations on the basis of a spherical earth (ignoring ellipsoidal effects) – which is accurate enough* for most purposes… [In fact, the earth is very slightly ellipsoidal; using a spherical model gives errors typically up to 0.3% – see notes for further details].
Enter the co-ordinates into the text boxes to try out the calculations. A variety of formats are accepted, principally:
deg-min-sec suffixed with N/S/E/W (e.g. 40°44′55″N, 73 59 11W), or
signed decimal degrees without compass direction, where negative indicates west/south (e.g. 40.7486, -73.9864):
Point 1: ,
Point 2: ,
Distance: 170.3 km
Initial bearing: 127°12′30″
Final bearing: 128°47′27″
Midpoint: 52°40′48″N, 000°50′37″W
And you can see it on a map (aren’t those Google guys wonderful!)
Distance
This uses the ‘Haversine’ formula to calculate great-circle distances between the two points – that is, the shortest distance over the earth’s surface – giving an ‘as-the-crow-flies’ distance between the points (ignoring any hills!).
Haversine formula:
R = earth’s radius (mean radius = 6,371km)
Δlat = lat2− lat1
Δlong = long2− long1
a = sin²(Δlat/2) + cos(lat1).cos(lat2).sin²(Δlong/2)
c = 2.atan2(√a, √(1−a))
d = R.c
(Note that angles need to be in radians to pass to trig functions).
JavaScript:
var R = 6371; // km
var dLat = (lat2-lat1).toRad();
var dLon = (lon2-lon1).toRad();
var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.cos(lat1.toRad()) * Math.cos(lat2.toRad()) *
Math.sin(dLon/2) * Math.sin(dLon/2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d = R * c;
The Haversine formula ‘remains particularly well-conditioned for numerical computation even at small distances’ – unlike calculations based on the spherical law of cosines. (It was published by R W Sinnott in Sky and Telescope, 1984, though has been known about for much longer; the ‘half-versed-sine’ is (1-cosθ)/2, or sin²(θ/2) – don’t ask, I’m not a mathematician).
Spherical Law of Cosines
In fact, when Sinnott published the Haversine formula, computational precision was limited. Nowadays, JavaScript (and most modern computers & languages) use IEEE 754 64-bit floating-point numbers, which provide 15 significant figures of precision. With this precision, the simple spherical law of cosines formula gives well-conditioned results down to distances as small as around 1 metre. In view of this it is probably worth, in most situations, using either the simpler law of cosines or the more accurate ellipsoidal Vincenty formula in preference to Haversine! (bearing in mind notes below on the limitations in accuracy of the spherical model).
Spherical law
of cosines: d = acos(sin(lat1).sin(lat2)+cos(lat1).cos(lat2).cos(long2−long1)).R
JavaScript:
var R = 6371; // km
var d = Math.acos(Math.sin(lat1)*Math.sin(lat2) +
Math.cos(lat1)*Math.cos(lat2) *
Math.cos(lon2-lon1)) * R;
Excel: =ACOS(SIN(lat1)*SIN(lat2)+COS(lat1)*COS(lat2)*COS(lon2-lon1))*6371
(Note that here and in all subsequent code fragments, for simplicity I do not show conversions from degrees to radians; see below for complete versions).
Bearing
In general, your current bearing will vary as you follow a great circle path (orthodrome); the final bearing will differ from the initial bearing by varying degrees according to distance and latitude (if you were to go from say 35°N,45°E (Baghdad) to 35°N,135°E (Osaka), you would start on a bearing of 60° and end up on a bearing of 120°!).
This is the initial bearing (sometimes referred to as forward azimuth) which if followed in a straight line along a great-circle arc will take you from the start point to the end point:
Formula: θ = atan2( sin(Δlong).cos(lat2),
cos(lat1).sin(lat2) − sin(lat1).cos(lat2).cos(Δlong) )
JavaScript:
var y = Math.sin(dLon) * Math.cos(lat2);
var x = Math.cos(lat1)*Math.sin(lat2) -
Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
var brng = Math.atan2(y, x).toDeg();
Excel: =ATAN2(COS(lat1)*SIN(lat2)-SIN(lat1)*COS(lat2)*COS(lon2-lon1),
SIN(lon2-lon1)*COS(lat2))
* Note that Excel reverses the arguments to ATAN2 – see notes below
Since atan2 returns values in the range -π ... +π (that is, -180° ... +180°), to normalise the result to a compass bearing (in the range 0° ... 360°, with -ve values transformed into the range 180° ... 360°), convert to degrees and then use (θ+360) % 360, where % is modulo.
For final bearing, simply take the initial bearing from the end point to the start point and reverse it (using θ = (θ+180) % 360).
Midpoint
This is the midpoint along a great circle path between the two points.
Formula: Bx = cos(lat2).cos(Δlong)
By = cos(lat2).sin(Δlong)
latm = atan2(sin(lat1) + sin(lat2), √((cos(lat1)+Bx)² + By²))
lonm = lon1 + atan2(By, cos(lat1)+Bx)
JavaScript:
var Bx = Math.cos(lat2) * Math.cos(dLon);
var By = Math.cos(lat2) * Math.sin(dLon);
var lat3 = Math.atan2(Math.sin(lat1)+Math.sin(lat2),
Math.sqrt( (Math.cos(lat1)+Bx)*(Math.cos(lat1)+Bx) + By*By) );
var lon3 = lon1 + Math.atan2(By, Math.cos(lat1) + Bx);
Just as the initial bearing may vary from the final bearing, the midpoint may not be located half-way between latitudes/longitudes; the midpoint between 35°N,45°E and 35°N,135°E is around 45°N,90°E.
Destination point given distance and bearing from start point
Given a start point, initial bearing, and distance, this will calculate the destination point and final bearing travelling along a (shortest distance) great circle arc.
Start point: ,
Bearing:
Distance: km
Destination point: 53°11′18″N, 000°08′00″E
Final bearing: 097°30′52″
view map
Formula: lat2 = asin(sin(lat1)*cos(d/R) + cos(lat1)*sin(d/R)*cos(θ))
lon2 = lon1 + atan2(sin(θ)*sin(d/R)*cos(lat1), cos(d/R)−sin(lat1)*sin(lat2))
d/R is the angular distance (in radians), where d is the distance travelled and R is the earth’s radius
JavaScript:
var lat2 = Math.asin( Math.sin(lat1)*Math.cos(d/R) +
Math.cos(lat1)*Math.sin(d/R)*Math.cos(brng) );
var lon2 = lon1 + Math.atan2(Math.sin(brng)*Math.sin(d/R)*Math.cos(lat1),
Math.cos(d/R)-Math.sin(lat1)*Math.sin(lat2));
Excel: lat2: =ASIN(SIN(lat1)*COS(d/R) + COS(lat1)*SIN(d/R)*COS(brng))
lon2: =lon1 + ATAN2(COS(d/R)-SIN(lat1)*SIN(lat2), SIN(brng)*SIN(d/R)*COS(lat1))
For final bearing, simply take the initial bearing from the end point to the start point and reverse it (using θ = (θ+180) % 360).
Intersection of two paths given start points and bearings
This is a rather more complex calculation than most others on this page, but I've been asked for it a number of times. See below for the JavaScript.
Point 1: , Brng 1:
Point 2: , Brng 2:
Intersection point: 50°54′06″N, 004°29′39″E
Formula:
d12 = 2.asin( √(sin²(Δlat/2) + cos(lat1).cos(lat2).sin²(Δlon/2)) )
φ1 = acos( sin(lat2) − sin(lat1).cos(d12) / sin(d12).cos(lat1) )
φ2 = acos( sin(lat1) − sin(lat2).cos(d12) / sin(d12).cos(lat2) )
if sin(lon2−lon1) > 0
θ12 = φ1, θ21 = 2.π − φ2
else
θ12 = 2.π − φ1, θ21 = φ2
α1 = (θ1 − θ12 + π) % 2.π − π
α2 = (θ21 − θ2 + π) % 2.π − π
α3 = acos( −cos(α1).cos(α2) + sin(α1).sin(α2).cos(d12) )
d13 = atan2( sin(d12).sin(α1).sin(α2), cos(α2)+cos(α1).cos(α3) )
lat3 = asin( sin(lat1).cos(d13) + cos(lat1).sin(d13).cos(θ1) )
Δlon = atan2( sin(θ1).sin(d13).cos(lat1), cos(d13)−sin(lat1).sin(lat3) )
lon3 = (lon1−Δlon+π) % 2.π − π
where
lat1, lon1, θ1 : 1st point & bearing
lat2, lon2, θ2 : 2nd point & bearing
lat3, lon3 : intersection point
% = mod
note – if sin(α1)=0 and sin(α2)=0: infinite solutions
if sin(α1).sin(α2) < 0: ambiguous solution
this formulation is not always well-conditioned for meridional or equatorial lines
Note this can also be solved using vectors rather than trigonometry:
For each point φ,λ (lat=φ, lon=λ), we can define a unit vector pointing to it from the centre of the earth: u{x,y,z} = [ cosφ·cosλ, cosφ·sinλ, sinφ ] (taking x=0º, y=90º, z=north – note that these formulæ depend on convention used for directions and handedness)
And for any great circle defined by two points, we can define a unit vector N normal to the plane of the circle: N(u1, u2) = (u1×u2) / ||u1×u2|| where × is the vector cross product, and ||u|| the norm (length of the vector)
The vector representing the intersection of the two great circles is then ui = ±N( N(u1, u2), N(u3, u4) )
We can then get the latitude and longitude of Pi by φ = atan2(uz, sqrt(ux² + uy²)), λ = atan2(uy, ux)
The antipodal intersection point is (-φ, λ+π)
Cross-track distance
Here’s a new one: I’ve sometimes been asked about distance of a point from a great-circle path (sometimes called cross track error).
Formula: dxt = asin(sin(d13/R)*sin(θ13−θ12)) * R
where d13 is distance from start point to third point
θ13 is (initial) bearing from start point to third point
θ12 is (initial) bearing from start point to end point
R is the earth’s radius
JavaScript:
var dXt = Math.asin(Math.sin(d13/R)*Math.sin(brng13-brng12)) * R;
Here, the great-circle path is identified by a start point and an end point – depending on what initial data you’re working from, you can use the formulæ above to obtain the relevant distance and bearings. The sign of dxt tells you which side of the path the third point is on.
The along-track distance, from the start point to the closest point on the path to the third point, is
Formula: dat = acos(cos(d13/R)/cos(dxt/R)) * R
where d13 is distance from start point to third point
dxt is cross-track distance
R is the earth’s radius
JavaScript:
var dAt = Math.acos(Math.cos(d13/R)/Math.cos(dXt/R)) * R;
Closest point to the poles
And: ‘Clairaut’s formula’ will give you the maximum latitude of a great circle path, given a bearing and latitude on the great circle:
Formula: latmax = acos(abs(sin(θ)*cos(lat)))
JavaScript:
var latMax = Math.acos(Math.abs(Math.sin(brng)*Math.cos(lat)));
Rhumb lines
A ‘rhumb line’ (or loxodrome) is a path of constant bearing, which crosses all meridians at the same angle.
Sailors used to (and sometimes still) navigate along rhumb lines since it is easier to follow a constant compass bearing than to be continually adjusting the bearing, as is needed to follow a great circle. Rhumb lines are straight lines on a Mercator Projection map (also helpful for navigation).
Rhumb lines are generally longer than great-circle (orthodrome) routes. For instance, London to New York is 4% longer along a rhumb line than along a great circle – important for aviation fuel, but not particularly to sailing vessels. New York to Beijing – close to the most extreme example possible (though not sailable!) – is 30% longer along a rhumb line.
Point 1: ,
Point 2: ,
Distance: 5196 km
Bearing: 260°07′38″
view map
Start point: ,
Bearing:
Distance: km
Destination point: 50°57′48″N, 001°51′09″E
view map
Distance/bearing
These formulæ give the distance and (constant) bearing between two points.
Formula: Δφ = ln(tan(lat2/2+π/4)/tan(lat1/2+π/4)) [= the ‘stretched’ latitude difference]
if E:W line, q = cos(lat1)
otherwise, q = Δlat/Δφ
d = √(Δlat² + q².Δlon²).R [pythagoras]
θ = atan2(Δlon, Δφ)
where ln is natural log, Δlon is taking shortest route (<180º), and R is the earth’s radius
JavaScript:
var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
var q = (!isNaN(dLat/dPhi)) ? dLat/dPhi : Math.cos(lat1); // E-W line gives dPhi=0
// if dLon over 180° take shorter rhumb across 180° meridian:
if (Math.abs(dLon) > Math.PI) {
dLon = dLon>0 ? -(2*Math.PI-dLon) : (2*Math.PI+dLon);
}
var d = Math.sqrt(dLat*dLat + q*q*dLon*dLon) * R;
var brng = Math.atan2(dLon, dPhi);
Destination
Given a start point and a distance d along constant bearing θ, this will calculate the destination point. If you maintain a constant bearing along a rhumb line, you will gradually spiral in towards one of the poles.
Formula: α = d/R (angular distance)
lat2 = lat1 + α.cos(θ)
Δφ = ln(tan(lat2/2+π/4)/tan(lat1/2+π/4)) [= the ‘stretched’ latitude difference]
if E:W line q = cos(lat1)
otherwise q = Δlat/Δφ
Δlon = α.sin(θ)/q
lon2 = (lon1+Δlon+π) % 2.π − π
where ln is natural log and % is modulo, Δlon is taking shortest route (<180°), and R is the earth’s radius
JavaScript:
lat2 = lat1 + d*Math.cos(brng);
var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
var q = (!isNaN(dLat/dPhi)) ? dLat/dPhi : Math.cos(lat1); // E-W line gives dPhi=0
var dLon = d*Math.sin(brng)/q;
// check for some daft bugger going past the pole, normalise latitude if so
if (Math.abs(lat2) > Math.PI/2) lat2 = lat2>0 ? Math.PI-lat2 : -(Math.PI-lat2);
lon2 = (lon1+dLon+Math.PI)%(2*Math.PI) - Math.PI;
If you use Ordnance Survey Grid References, I have implemented a script for converting between Lat/Long & OS Grid References.
Convert between degrees-minutes-seconds & decimal degrees
Latitude Longitude 1° ≈ 111 km (110.57 eq’l — 111.70 polar)
1′ ≈ 1.85 km (= 1 nm) 0.01° ≈ 1.11 km
1″ ≈ 30.9 m 0.0001° ≈ 11.1 m
No, I’ve not included decimal minutes: a decimal system is easy, a sexagesimal system has merits, but mixing the two is a complete sow’s ear. Switch off the option on your GPS!
Notes:
Accuracy: since the earth is not quite a sphere, there are small errors in using spherical geometry; the earth is actually roughly ellipsoidal (or more precisely, oblate spheroidal) with a radius varying between about 6,378km (equatorial) and 6,357km (polar), and local radius of curvature varying from 6,336km (equatorial meridian) to 6,399km (polar). 6,371 km is the generally accepted value for the Earth’s mean radius. This means that errors from assuming spherical geometry might be up to 0.55% crossing the equator, though generally below 0.3%, depending on latitude and direction of travel. An accuracy of better than 3m in 1km is mostly good enough for me, but if you want greater accuracy, you could use the Vincenty formula for calculating geodesic distances on ellipsoids, which gives results accurate to within 1mm. (Out of sheer perversity – I’ve never needed such accuracy – I looked up this formula and discovered the JavaScript implementation was simpler than I expected).
Trig functions take arguments in radians, so latitude, longitude, and bearings in degrees (either decimal or degrees/minutes/seconds) need to be converted to radians, rad = π.deg/180. When converting radians back to degrees (deg = 180.rad/π), West is negative if using signed decimal degrees. For bearings, values in the range -π to +π [-180° to +180°] need to be converted to 0 to +2π [0°–360°]; this can be done by (brng+2.π)%2.π [or brng+360)%360] where % is the modulo operator.
The atan2() function widely used here takes two arguments, atan2(y, x), and computes the arc tangent of the ratio y/x. It is more flexible than atan(y/x), since it handles x=0, and it also returns values in all 4 quadrants -π to +π (the atan function returns values in the range -π/2 to +π/2).
If you implement any formula involving atan2 in Microsoft Excel, you will need to reverse the arguments, as Excel has them the opposite way around from JavaScript – conventional order is atan2(y, x), but Excel uses atan2(x, y). To use atan2 in a (VBA) macro, you can use WorksheetFunction.Atan2().
All bearings are with respect to true north, 0°=N, 90°=E, etc; if you are working from a compass, magnetic north varies from true north in a complex way around the earth, and the difference has to be compensated for by variances indicated on local maps.
I learned a lot from the US Census Bureau GIS FAQ which is no longer available, so I’ve made a copy.
Thanks to Ed Williams’ Aviation Formulary for many of the formulæ.
For miles, divide km by 1.609344
For nautical miles, divide km by 1.852
See below for the source code of the JavaScript implementation. These functions should be simple to translate into other languages if required.
Update January 2010: I have revised the scripts to be structured as methods of a LatLon object. Of course, JavaScript is a prototype-based rather than class-based language, so this is only nominally a class, but isolating code into a separate namespace is good JavaScript practice, and this approach may also make it clearer to implement these functions in other languages. If you’re not familiar with JavaScript syntax, LatLon.prototype.distanceTo = function(point) { ... }, for instance, defines a ‘distanceTo’ method of the LatLon object (/class) which takes a LatLon object as a parameter (and returns a number). The Geo namespace acts as a static class for geodesy formatting / parsing / conversion functions. I have extended (polluted, if you like) the base JavaScript object prototypes with trim(), toRad() toDeg(), and toPrecisionFixed() methods. I’ve adopted JSDoc format for the descriptions.
I have also created a page illustrating the use of the spherical law of cosines for selecting points from a database within a specified bounding circle – the example is based on MySQL+PDO, but should be extensible to other DBMS platforms.
Several people have asked about example Excel spreadsheets, so I have implemented the distance & bearing and the destination point formulæ as spreadsheets, in a form which breaks down the all stages involved to illustrate the operation.
I offer these formulæ & scripts for free use and adaptation as my contribution to the open-source info-sphere from which I have received so much. You are welcome to re-use these scripts [under a simple attribution license, without any warranty express or implied] provided solely that you retain my copyright notice and a link to this page.
If you would like to show your appreciation and support continued development of these scripts, I would most gratefully accept donations.
If you need any advice or development work done, I am available for consultancy.
If you have any queries or find any problems, contact me at ku.oc.epyt-elbavom@oeg-stpircs.
© 2002-2010 Chris Veness
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Latitude/longitude spherical geodesy formulae & scripts (c) Chris Veness 2002-2010 */
/* - www.movable-type.co.uk/scripts/latlong.html */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Note that minimal error checking is performed in this example code! */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* Creates a point on the earth's surface at the supplied latitude / longitude
*
* @constructor
* @param {Number} lat: latitude in numeric degrees
* @param {Number} lon: longitude in numeric degrees
* @param {Number} [rad=6371]: radius of earth if different value is required from standard 6,371km
*/
function LatLon(lat, lon, rad) {
if (typeof rad == 'undefined') rad = 6371; // earth's mean radius in km
this._lat = lat;
this._lon = lon;
this._radius = rad;
}
/**
* Returns the distance from this point to the supplied point, in km
* (using Haversine formula)
*
* from: Haversine formula - R. W. Sinnott, "Virtues of the Haversine",
* Sky and Telescope, vol 68, no 2, 1984
*
* @param {LatLon} point: Latitude/longitude of destination point
* @param {Number} [precision=4]: no of significant digits to use for returned value
* @returns {Number} Distance in km between this point and destination point
*/
LatLon.prototype.distanceTo = function(point, precision) {
// default 4 sig figs reflects typical 0.3% accuracy of spherical model
if (typeof precision == 'undefined') precision = 4;
var R = this._radius;
var lat1 = this._lat.toRad(), lon1 = this._lon.toRad();
var lat2 = point._lat.toRad(), lon2 = point._lon.toRad();
var dLat = lat2 - lat1;
var dLon = lon2 - lon1;
var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.cos(lat1) * Math.cos(lat2) *
Math.sin(dLon/2) * Math.sin(dLon/2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
var d = R * c;
return d.toPrecisionFixed(precision);
}
/**
* Returns the (initial) bearing from this point to the supplied point, in degrees
* see http://williams.best.vwh.net/avform.htm#Crs
*
* @param {LatLon} point: Latitude/longitude of destination point
* @returns {Number} Initial bearing in degrees from North
*/
LatLon.prototype.bearingTo = function(point) {
var lat1 = this._lat.toRad(), lat2 = point._lat.toRad();
var dLon = (point._lon-this._lon).toRad();
var y = Math.sin(dLon) * Math.cos(lat2);
var x = Math.cos(lat1)*Math.sin(lat2) -
Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
var brng = Math.atan2(y, x);
return (brng.toDeg()+360) % 360;
}
/**
* Returns final bearing arriving at supplied destination point from this point; the final bearing
* will differ from the initial bearing by varying degrees according to distance and latitude
*
* @param {LatLon} point: Latitude/longitude of destination point
* @returns {Number} Final bearing in degrees from North
*/
LatLon.prototype.finalBearingTo = function(point) {
// get initial bearing from supplied point back to this point...
var lat1 = point._lat.toRad(), lat2 = this._lat.toRad();
var dLon = (this._lon-point._lon).toRad();
var y = Math.sin(dLon) * Math.cos(lat2);
var x = Math.cos(lat1)*Math.sin(lat2) -
Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
var brng = Math.atan2(y, x);
// ... & reverse it by adding 180°
return (brng.toDeg()+180) % 360;
}
/**
* Returns the midpoint between this point and the supplied point.
* see http://mathforum.org/library/drmath/view/51822.html for derivation
*
* @param {LatLon} point: Latitude/longitude of destination point
* @returns {LatLon} Midpoint between this point and the supplied point
*/
LatLon.prototype.midpointTo = function(point) {
lat1 = this._lat.toRad(), lon1 = this._lon.toRad();
lat2 = point._lat.toRad();
var dLon = (point._lon-this._lon).toRad();
var Bx = Math.cos(lat2) * Math.cos(dLon);
var By = Math.cos(lat2) * Math.sin(dLon);
lat3 = Math.atan2(Math.sin(lat1)+Math.sin(lat2),
Math.sqrt( (Math.cos(lat1)+Bx)*(Math.cos(lat1)+Bx) + By*By) );
lon3 = lon1 + Math.atan2(By, Math.cos(lat1) + Bx);
return new LatLon(lat3.toDeg(), lon3.toDeg());
}
/**
* Returns the destination point from this point having travelled the given distance (in km) on the
* given initial bearing (bearing may vary before destination is reached)
*
* see http://williams.best.vwh.net/avform.htm#LL
*
* @param {Number} brng: Initial bearing in degrees
* @param {Number} dist: Distance in km
* @returns {LatLon} Destination point
*/
LatLon.prototype.destinationPoint = function(brng, dist) {
dist = dist/this._radius; // convert dist to angular distance in radians
brng = brng.toRad(); //
var lat1 = this._lat.toRad(), lon1 = this._lon.toRad();
var lat2 = Math.asin( Math.sin(lat1)*Math.cos(dist) +
Math.cos(lat1)*Math.sin(dist)*Math.cos(brng) );
var lon2 = lon1 + Math.atan2(Math.sin(brng)*Math.sin(dist)*Math.cos(lat1),
Math.cos(dist)-Math.sin(lat1)*Math.sin(lat2));
lon2 = (lon2+3*Math.PI)%(2*Math.PI) - Math.PI; // normalise to -180...+180
if (isNaN(lat2) || isNaN(lon2)) return null;
return new LatLon(lat2.toDeg(), lon2.toDeg());
}
/**
* Returns the point of intersection of two paths defined by point and bearing
*
* see http://williams.best.vwh.net/avform.htm#Intersection
*
* @param {LatLon} p1: First point
* @param {Number} brng1: Initial bearing from first point
* @param {LatLon} p2: Second point
* @param {Number} brng2: Initial bearing from second point
* @returns {LatLon} Destination point (null if no unique intersection defined)
*/
LatLon.intersection = function(p1, brng1, p2, brng2) {
lat1 = p1._lat.toRad(), lon1 = p1._lon.toRad();
lat2 = p2._lat.toRad(), lon2 = p2._lon.toRad();
brng13 = brng1.toRad(), brng23 = brng2.toRad();
dLat = lat2-lat1, dLon = lon2-lon1;
dist12 = 2*Math.asin( Math.sqrt( Math.sin(dLat/2)*Math.sin(dLat/2) +
Math.cos(lat1)*Math.cos(lat2)*Math.sin(dLon/2)*Math.sin(dLon/2) ) );
if (dist12 == 0) return null;
// initial/final bearings between points
brngA = Math.acos( ( Math.sin(lat2) - Math.sin(lat1)*Math.cos(dist12) ) /
( Math.sin(dist12)*Math.cos(lat1) ) );
if (isNaN(brngA)) brngA = 0; // protect against rounding
brngB = Math.acos( ( Math.sin(lat1) - Math.sin(lat2)*Math.cos(dist12) ) /
( Math.sin(dist12)*Math.cos(lat2) ) );
if (Math.sin(lon2-lon1) > 0) {
brng12 = brngA;
brng21 = 2*Math.PI - brngB;
} else {
brng12 = 2*Math.PI - brngA;
brng21 = brngB;
}
alpha1 = (brng13 - brng12 + Math.PI) % (2*Math.PI) - Math.PI; // angle 2-1-3
alpha2 = (brng21 - brng23 + Math.PI) % (2*Math.PI) - Math.PI; // angle 1-2-3
if (Math.sin(alpha1)==0 && Math.sin(alpha2)==0) return null; // infinite intersections
if (Math.sin(alpha1)*Math.sin(alpha2) < 0) return null; // ambiguous intersection
//alpha1 = Math.abs(alpha1);
//alpha2 = Math.abs(alpha2);
// ... Ed Williams takes abs of alpha1/alpha2, but seems to break calculation?
alpha3 = Math.acos( -Math.cos(alpha1)*Math.cos(alpha2) +
Math.sin(alpha1)*Math.sin(alpha2)*Math.cos(dist12) );
dist13 = Math.atan2( Math.sin(dist12)*Math.sin(alpha1)*Math.sin(alpha2),
Math.cos(alpha2)+Math.cos(alpha1)*Math.cos(alpha3) )
lat3 = Math.asin( Math.sin(lat1)*Math.cos(dist13) +
Math.cos(lat1)*Math.sin(dist13)*Math.cos(brng13) );
dLon13 = Math.atan2( Math.sin(brng13)*Math.sin(dist13)*Math.cos(lat1),
Math.cos(dist13)-Math.sin(lat1)*Math.sin(lat3) );
lon3 = lon1+dLon13;
lon3 = (lon3+Math.PI) % (2*Math.PI) - Math.PI; // normalise to -180..180º
return new LatLon(lat3.toDeg(), lon3.toDeg());
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* Returns the distance from this point to the supplied point, in km, travelling along a rhumb line
*
* see http://williams.best.vwh.net/avform.htm#Rhumb
*
* @param {LatLon} point: Latitude/longitude of destination point
* @returns {Number} Distance in km between this point and destination point
*/
LatLon.prototype.rhumbDistanceTo = function(point) {
var R = this._radius;
var lat1 = this._lat.toRad(), lat2 = point._lat.toRad();
var dLat = (point._lat-this._lat).toRad();
var dLon = Math.abs(point._lon-this._lon).toRad();
var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
var q = (!isNaN(dLat/dPhi)) ? dLat/dPhi : Math.cos(lat1); // E-W line gives dPhi=0
// if dLon over 180° take shorter rhumb across 180° meridian:
if (dLon > Math.PI) dLon = 2*Math.PI - dLon;
var dist = Math.sqrt(dLat*dLat + q*q*dLon*dLon) * R;
return dist.toPrecisionFixed(4); // 4 sig figs reflects typical 0.3% accuracy of spherical model
}
/**
* Returns the bearing from this point to the supplied point along a rhumb line, in degrees
*
* @param {LatLon} point: Latitude/longitude of destination point
* @returns {Number} Bearing in degrees from North
*/
LatLon.prototype.rhumbBearingTo = function(point) {
var lat1 = this._lat.toRad(), lat2 = point._lat.toRad();
var dLon = (point._lon-this._lon).toRad();
var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
if (Math.abs(dLon) > Math.PI) dLon = dLon>0 ? -(2*Math.PI-dLon) : (2*Math.PI+dLon);
var brng = Math.atan2(dLon, dPhi);
return (brng.toDeg()+360) % 360;
}
/**
* Returns the destination point from this point having travelled the given distance (in km) on the
* given bearing along a rhumb line
*
* @param {Number} brng: Bearing in degrees from North
* @param {Number} dist: Distance in km
* @returns {LatLon} Destination point
*/
LatLon.prototype.rhumbDestinationPoint = function(brng, dist) {
var R = this._radius;
var d = parseFloat(dist)/R; // d = angular distance covered on earth's surface
var lat1 = this._lat.toRad(), lon1 = this._lon.toRad();
brng = brng.toRad();
var lat2 = lat1 + d*Math.cos(brng);
var dLat = lat2-lat1;
var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
var q = (!isNaN(dLat/dPhi)) ? dLat/dPhi : Math.cos(lat1); // E-W line gives dPhi=0
var dLon = d*Math.sin(brng)/q;
// check for some daft bugger going past the pole
if (Math.abs(lat2) > Math.PI/2) lat2 = lat2>0 ? Math.PI-lat2 : -(Math.PI-lat2);
lon2 = (lon1+dLon+3*Math.PI)%(2*Math.PI) - Math.PI;
return new LatLon(lat2.toDeg(), lon2.toDeg());
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/**
* Returns the latitude of this point; signed numeric degrees if no format, otherwise format & dp
* as per Geo.toLat()
*
* @param {String} [format]: Return value as 'd', 'dm', 'dms'
* @param {Number} [dp=0|2|4]: No of decimal places to display
* @returns {Number|String} Numeric degrees if no format specified, otherwise deg/min/sec
*
* @requires Geo
*/
LatLon.prototype.lat = function(format, dp) {
if (typeof format == 'undefined') return this._lat;
return Geo.toLat(this._lat, format, dp);
}
/**
* Returns the longitude of this point; signed numeric degrees if no format, otherwise format & dp
* as per Geo.toLon()
*
* @param {String} [format]: Return value as 'd', 'dm', 'dms'
* @param {Number} [dp=0|2|4]: No of decimal places to display
* @returns {Number|String} Numeric degrees if no format specified, otherwise deg/min/sec
*
* @requires Geo
*/
LatLon.prototype.lon = function(format, dp) {
if (typeof format == 'undefined') return this._lon;
return Geo.toLon(this._lon, format, dp);
}
/**
* Returns a string representation of this point; format and dp as per lat()/lon()
*
* @param {String} [format]: Return value as 'd', 'dm', 'dms'
* @param {Number} [dp=0|2|4]: No of decimal places to display
* @returns {String} Comma-separated latitude/longitude
*
* @requires Geo
*/
LatLon.prototype.toString = function(format, dp) {
if (typeof format == 'undefined') format = 'dms';
return Geo.toLat(this._lat, format, dp) + ', ' + Geo.toLon(this._lon, format, dp);
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
// extend Number object with methods for converting degrees/radians
/** Convert numeric degrees to radians */
if (typeof(String.prototype.toRad) === "undefined") {
Number.prototype.toRad = function() {
return this * Math.PI / 180;
}
}
/** Convert radians to numeric (signed) degrees */
if (typeof(String.prototype.toDeg) === "undefined") {
Number.prototype.toDeg = function() {
return this * 180 / Math.PI;
}
}
/**
* Format the significant digits of a number, using only fixed-point notation (no exponential)
*
* @param {Number} precision: Number of significant digits to appear in the returned string
* @returns {String} A string representation of number which contains precision significant digits
*/
if (typeof(String.prototype.toDeg) === "undefined") {
Number.prototype.toPrecisionFixed = function(precision) {
var numb = this < 0 ? -this : this; // can't take log of -ve number...
var sign = this < 0 ? '-' : '';
if (numb == 0) { n = '0.'; while (precision--) n += '0'; return n }; // can't take log of zero
var scale = Math.ceil(Math.log(numb)*Math.LOG10E); // no of digits before decimal
var n = String(Math.round(numb * Math.pow(10, precision-scale)));
if (scale > 0) { // add trailing zeros & insert decimal as required
l = scale - n.length;
while (l-- > 0) n = n + '0';
if (scale < n.length) n = n.slice(0,scale) + '.' + n.slice(scale);
} else { // prefix decimal and leading zeros if required
while (scale++ < 0) n = '0' + n;
n = '0.' + n;
}
return sign + n;
}
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
3:34 PM