Home - Page  
Assn Members Area · Join

Calculate distance, bearing and more between two Latitude/Longitude points

Distance

This script calculates great-circle distances between the two points – that is, the shortest distance over the earth’s surface – using the ‘Haversine’ formula.

It assumes a spherical earth, ignoring ellipsoidal effects – which is accurate enough* for most purposes… – giving an ‘as-the-crow-flies’ distance between the two points (ignoring any hills!).

Enter the co-ordinates into the text boxes to try it out. It accepts a variety of formats:

  • 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):

Lat 1: Long 1:

Lat 2: Long 2:

And you can see it on a map (thanks to the nice guys at Google Maps)

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).

In fact, when Sinnott devised the Haversine formula, computational precision was limited. Nowadays, JavaScript (and most modern computers) 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! (See 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; View Source for complete versions).

Bearing

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).toBrng();
Excel: =ATAN2(COS(lat1)*SIN(lat2)-SIN(lat1)*COS(lat2)*COS(lon2-lon1), SIN(lon2-lon1)*COS(lat2))

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.

This is the initial bearing which if followed in a straight line along a great-circle arc (orthodrome) will take you from the start point to the end point; in general, the bearing you are following will have varied by the time you get to the end point (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°!).

For final bearing, simply take the initial bearing from the end point to the start point and reverse it (using θ = (θ+180) % 360).

Midpoint

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.toRad() + 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

This page is steadily growing! 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 Lat: Start Long:
Bearing (deg): Distance (km):
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/ER) + COS(lat1)*SIN(d/ER)*COS(brng))
lon2: =lon1 + ATAN2(COS(d/ER)-SIN(lat1)*SIN(lat2), SIN(brng)*SIN(d/ER)*COS(lat1))

For final bearing, simply take the initial bearing from the end point to the start point and reverse it (using θ = (θ+180) % 360).

 


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 formulae 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

Also new: ‘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 continually adjust the bearing as is needed to follow a great circle, though they are normally longer than great-circle (orthodrome) routes. Rhumb lines are straight lines on a Mercator Projection map.

Lat 1: Long 1:
Lat 2: Long 2:
Formula: Δφ = ln(tan(lat2/2+π/4)/tan(lat1/2+π/4))
if E:W line q = cos(lat1)
otherwise q = Δlat/Δφ
  d = √[Δlat² + q².Δlon²].R
  θ = 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 = (Math.abs(dLat) > 1e-10) ? dLat/dPhi : Math.cos(lat1);
// 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);

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.

Start Lat: Start Long:
Bearing (deg): Distance (km):

 

Formula: α = d/R (angular distance)
  lat2 = lat1 + α.cos(θ)
  Δφ = ln(tan(lat2/2+π/4)/tan(lat1/2+π/4))
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 = (Math.abs(lat2-lat1) > 1e-10) ? (lat2-lat1)/dPhi : Math.cos(lat1);
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.001° ˜ 111 m
1? ˜ 30.9 m 0.00001° ˜ 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!


Original material published with Thanks to Moveable Type