Movable Type Home Page

Movable Type Scripts


Calculate distance, bearing and more between Latitude/Longitude points

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

Great-circle distance between two points

Enter the co-ordinates into the text boxes to try out the calculations. A variety of formats are accepted, principally:

Point 1: ,
Point 2: ,
Distance: km (to 4 SF*)
Initial bearing:
Final bearing:
Midpoint:

And you can see it on a map (aren’t those Google guys wonderful!)

Distance

This uses the ‘haversine’ formula to calculate the great-circle distance between 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, of course!).

Haversine
formula:
a = sin²(Δφ/2) + cos(φ1).cos(φ2).sin²(Δλ/2)
c = 2.atan2(√a, √(1−a))
d = R.c
where φ is latitude, λ is longitude, R is earth’s radius (mean radius = 6,371km)
  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 lat1 = lat1.toRad();
var lat2 = lat2.toRad();

var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
        Math.sin(dLon/2) * Math.sin(dLon/2) * Math.cos(lat1) * Math.cos(lat2); 
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); 
var d = R * c;

The haversine formula1 ‘remains particularly well-conditioned for numerical computation even at small distances’ – unlike calculations based on the spherical law of cosines. The ‘versed sine’ is 1-cosθ, and the ‘half-versed-sine’ (1-cosθ)/2 = sin²(θ/2) as used above. It was published by Roger Sinnott in Sky & Telescope magazine in 1984 (“Virtues of the Haversine”), though known about for much longer by navigators. (For the curious, c is the angular distance in radians, and a is the square of half the chord length between the points). A (surprisingly marginal) performance improvement can be obtained, of course, by factoring out the terms which get squared.

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 (cos c = cos a cos b + sin a sin b cos C) gives well-conditioned results down to distances as small as around 1 metre. (Note that the geodetic form of the law of cosines is rearranged from the canonical one so that the latitude can be used directly, rather than the colatitude).

This makes the simpler law of cosines a reasonable 1-line alternative to the haversine formula for many purposes. The choice may be driven by coding context, available trig functions (in different languages), etc.

Spherical
law of cosines:
d = acos( sin(φ1).sin(φ2) + cos(φ1).cos(φ2).cos(Δλ) ).R
JavaScript:
var φ1 = lat1.toRad(), φ2 = lat2.toRad(), Δλ = (lon2-lon1).toRad(), R = 6371; // gives d in km
var d = Math.acos( Math.sin(φ1)*Math.sin(φ2) +
                   Math.cos(φ1)*Math.cos(φ2) * Math.cos(Δλ) ) * 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 always show conversions from degrees to radians; see code below for complete versions).

Equirectangular approximation

If performance is an issue and accuracy less important, for small distances Pythagoras’ theorem can be used on an equirectangular projection:*

Formula x = Δλ.cos(φ)
y = Δφ
d = R.√x² + y²
JavaScript:
var x = (lon2-lon1) * Math.cos((lat1+lat2)/2);
var y = (lat2-lat1);
var d = Math.sqrt(x*x + y*y) * R;
(lat/lon in radians!)

This uses just one trig and one sqrt function – as against half-a-dozen trig functions for cos law, and 7 trigs + 2 sqrts for haversine. Accuracy is somewhat complex: along meridians there are no errors, otherwise they depend on distance, bearing, and latitude, but are small enough for many purposes* (and often trivial compared with the spherical approximation itself).

Bearing

Baghdad to Osaka
Baghdad to Osaka –
not a constant bearing!

In general, your current heading will vary as you follow a great circle path (orthodrome); the final heading will differ from the initial heading 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 heading of 60° and end up on a heading of 120°!).

This formula is for 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:1

Formula: θ = atan2( sin(Δλ).cos(φ2), cos(φ1).sin(φ2) − sin(φ1).cos(φ2).cos(Δλ) )
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 half-way point along a great circle path between the two points.1

Formula: Bx = cos(φ2).cos(Δλ)
By = cos(φ2).sin(Δλ)
φm = atan2( sin(φ1) + sin(φ2), √((cos(φ1)+Bx)² + By²) )
λm = λ1 + atan2(By, cos(φ1)+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.

Destination point along great-circle given distance and bearing from start point
Start point: ,
Bearing:
Distance: km
Destination point:
Final bearing:

view map

 
Formula: φ2 = asin( sin(φ1)*cos(d/R) + cos(φ1)*sin(d/R)*cos(θ) )
  λ2 = λ1 + atan2( sin(θ)*sin(d/R)*cos(φ1), cos(d/R)−sin(φ1)*sin(φ2) )
where φ is latitude, λ is longitude, θ is the bearing (in radians, clockwise from north), d is the distance travelled, R is the earth’s radius (d/R is the angular distance, in radians)
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))
* Remember that Excel reverses the arguments to ATAN2 – see notes below

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.

Intersection of two great-circle paths
Point 1: , Brng 1:
Point 2: , Brng 2:
Intersection point:

 

Formula:

d12 = 2.asin( √(sin²(Δφ/2) + cos(φ1).cos(φ2).sin²(Δλ/2)) )
φ1 = acos( sin(φ2) − sin(φ1).cos(d12) / sin(d12).cos(φ1) )
φ2 = acos( sin(φ1) − sin(φ2).cos(d12) / sin(d12).cos(φ2) )

if sin(λ2−λ1) > 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) )
φ3 = asin( sin(φ1).cos(d13) + cos(φ1).sin(d13).cos(θ1) )
Δλ13 = atan2( sin(θ1).sin(d13).cos(φ1), cos(d13)−sin(φ1).sin(φ3) )
λ3 = (λ1+Δλ13+π) % 2.π − π

where

φ1, λ1, θ1 : 1st point & bearing
φ2, λ2, θ2 : 2nd point & bearing
φ3, λ3 : 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:


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: φmax = acos(abs(sin(θ)*cos(φ)))
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.

Rhumb-line distance between two points
Point 1: ,
Point 2: ,
Distance: km
Bearing:
Midpoint:

view map

 
Destination point along rhumb line given distance and bearing from start point
Start point: ,
Bearing:
Distance: km
Destination point:

view map

 

 

Distance/bearing

These formulæ give the distance and (constant) bearing between two points.

Formula: Δφ′ = ln( tan(π/4+φ2/2)/tan(π/4+φ1/2) ) (the ‘stretched’ latitude difference)
if E:W line, q = cos(φ1)  
otherwise, q = Δφ/Δφ′  
  d = √(Δφ² + q².Δλ²).R (pythagoras)
  θ = atan2(Δλ, Δφ′)  
where φ is latitude, λ is longitude, Δλ is taking shortest route (<180º), R is the earth’s radius, ln is natural log
JavaScript:
var dPhi = Math.log(Math.tan(Math.PI/4+lat2/2)/Math.tan(Math.PI/4+lat1/2));
var q = (isFinite(dLat/dPhi)) ? dLat/dPhi : Math.cos(lat1);  // E-W line gives dPhi=0

// if dLon over 180° take shorter rhumb across anti-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)
  φ2 = φ1 + α.cos(θ)  
  Δφ′ = ln( tan(π/4+φ2/2) / tan(π/4+φ1/2) ) (the ‘stretched’ latitude difference)
if E:W line, q = cos(φ1)  
otherwise, q = Δφ/Δφ′  
  Δλ = α.sin(θ)/q  
  λ2 = (λ1+Δλ+π) % 2.π − π  
where φ is latitude, λ is longitude, Δλ is taking shortest route (<180°), ln is natural log and % is modulo, R is the earth’s radius
JavaScript:
var dLat = d*Math.cos(brng);
var lat2 = lat1 + dLat;
var dPhi = Math.log(Math.tan(lat2/2+Math.PI/4)/Math.tan(lat1/2+Math.PI/4));
var q = (isFinite(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;

Mid-point

This formula for calculating the ‘loxodromic midpoint’, the point half-way along a rhumb line between two points, is due to Robert Hill and Clive Tooth1 (thx Axel!).

Formula: φm = (φ12)/2  
  f1 = tan(π/4+φ1/2)  
  f2 = tan(π/4+φ2/2)  
  fm = tan(π/4+φm/2)  
  λm = [ (λ2−λ1).ln(fm) + λ1.ln(f2) − λ2.ln(f1) ] / ln(f2/f1)  
where φ is latitude, λ is longitude, ln is natural log
JavaScript:
if (Math.abs(lon2-lon1) > Math.PI) lon1 += 2*Math.PI; // crossing anti-meridian

var lat3 = (lat1+lat2)/2;
var f1 = Math.tan(Math.PI/4 + lat1/2);
var f2 = Math.tan(Math.PI/4 + lat2/2);
var f3 = Math.tan(Math.PI/4 + lat3/2);
var lon3 = ( (lon2-lon1)*Math.log(f3) + lon1*Math.log(f2) - lon2*Math.log(f1) ) / 
   Math.log(f2/f1);

if (!isFinite(lon3)) lon3 = (lon1+lon2)/2; // parallel of latitude

lon3 = (lon3+3*Math.PI) % (2*Math.PI) - Math.PI;  // normalise to -180..+180º
    

Using the scripts in web pages

Using these scripts in web pages would be something like the following:

<script src="latlon.js">/* Latitude/Longitude formulae */</script>
<script src="geo.js">/* Geodesy representation conversions */</script>
...
<form>
  Lat1: <input type="text" name="lat1" id="lat1"> Lon1: <input type="text" name="lon1" id="lon1">
  Lat2: <input type="text" name="lat2" id="lat2"> Lon2: <input type="text" name="lon2" id="lon2">
  <button onClick="var p1 = new LatLon(Geo.parseDMS(f.lat1.value), Geo.parseDMS(f.lon1.value));
                   var p2 = new LatLon(Geo.parseDMS(f.lat2.value), Geo.parseDMS(f.lon2.value));
                   alert(p1.distanceTo(p2));">Calculate distance</button>
</form>

If you use jQuery, the code can be separated from the HTML:

<script src="http://ajax.googleapis.com/ajax/libs/jquery/1.6.0/jquery.min.js"></script>
<script src="latlon.js">/* Latitude/Longitude formulae */</script>
<script src="geo.js">/* Geodesy representation conversions */</script>
<script>
  $(document).ready(function() {
    $('#calc-dist').click(function() {
      var p1 = new LatLon(Geo.parseDMS($('#lat1').val()), Geo.parseDMS($('#lon1').val()));
      var p2 = new LatLon(Geo.parseDMS($('#lat2').val()), Geo.parseDMS($('#lon2').val()));
      $('#result-distance').html(p1.distanceTo(p2)+' km');
    });
  });
</script>
...
<form>
  Lat1: <input type="text" name="lat1" id="lat1"> Lon1: <input type="text" name="lon1" id="lon1">
  Lat2: <input type="text" name="lat2" id="lat2"> Lon2: <input type="text" name="lon2" id="lon2">
  <button id="calc-dist">Calculate distance</button>
  <output id="result-distance"></output>
</form>

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!

Display calculation results as: deg/min/sec decimal degrees

Notes:


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.

Creative Commons License 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 reference to this page.

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-2014 Chris Veness