Movable Type Home Page

Movable Type Scripts


Vincenty formula for distance between two Latitude/Longitude points

 

For the benefit of the terminally obsessive (as well as the genuinely needy), Thaddeus Vincenty (‘TV’) devised formulae for calculating geodesic distances between a pair of latitude/longitude points on the earth’s surface, using an accurate ellipsoidal model of the earth.

When I looked up the references (‘Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations’), I discovered to my surprise that while the mathematics is utterly beyond me, it is actually quite simple to program.

Vincenty’s formula is accurate to within 0.5mm, or 0.000015″ (!), on the ellipsoid being used. Calculations based on a spherical model, such as the (much simpler) Haversine, are accurate to around 0.3% (which is still good enough for most purposes, of course).

Note: the accuracy quoted by Vincenty applies to the theoretical ellipsoid being used, which will differ (to varying degree) from the real earth geoid. If you happen to be located in Colorado, 2km above msl, distances will be 0.03% greater. In the UK, if you measure the distance from Land’s End to John O’ Groats using WGS-84, it will be 28m – 0.003% – less than using the Airy ellipsoid, which provides a better fit for the UK. (Since of course no one can travel on the surface of the theoretical ellipsoid, these differences are generally of no relevance at all).

Functional demo

Enter the co-ordinates into the text boxes to try it out (using deg-min-sec suffixed with N/S/E/W, or signed decimal degrees):

Lat 1: Long 1:

Lat 2: Long 2:

Note: nearly-antipodal points may fail to give a solution, in which case NaN is returned.

Vincenty wrote this at a time when computing resources were very expensive, and made it very computationally efficient. The JavaScript should be quite straightforward to translate into other languages, if required.

See the ‘Direct’ formula to calculate the destination point given the start point, bearing and distance.


Vincenty’s formula as it is used in the script:

a, b = major & minor semiaxes of the ellipsoid  
f = flattening (ab)/a  
φ1, φ2 = geodetic latitude  
L = difference in longitude  
U1 = atan((1−f).tanφ1) (U is ‘reduced latitude’)  
U2 = atan((1−f).tanφ2)  
λ = L (first approximation)  
iterate until change in λ is negligible (e.g. 10-12 ≈ 0.006mm) {  
     sinσ = √[ (cosU2.sinλ)² + (cosU1.sinU2 − sinU1.cosU2.cosλ)² ] (14)
  cosσ = sinU1.sinU2 + cosU1.cosU2.cosλ (15)
  σ = atan2(sinσ, cosσ) (16)
  sinα = cosU1.cosU2.sinλ / sinσ (17)
  cos²α = 1 − sin²α (trig identity; §6)  
  cos2σm = cosσ − 2.sinU1.sinU2/cos²α (18)
  C = f/16.cos²α.[4+f.(4−3.cos²α)] (10)
  λ′ = L + (1−C).f.sinα.{σ+C.sinσ.[cos2σm+C.cosσ.(−1+2.cos²2σm)]} (11)
}    
u² = cos²α.(a²−b²)/b²  
A = 1+u²/16384.{4096+u².[−768+u².(320−175.u²)]} (3)
B = u²/1024.{256+u².[−128+u².(74−47.u²)]} (4)
Δσ = B.sinσ.{cos2σm+B/4.[cosσ.(−1+2.cos²2σm) − B/6.cos2σm.(−3+4.sin²σ).(−3+4.cos²2σm)]} (6)
s = b.A.(σ−Δσ) (19)
α1 = atan2(cosU2.sinλ, cosU1.sinU2 − sinU1.cosU2.cosλ) (20)
α2 = atan2(cosU1.sinλ, −sinU1.cosU2 + cosU1.sinU2.cosλ) (21)

Where:

Note: Vincenty observes that eqn. (18) becomes indeterminate over equatorial lines (since cos²α → 0); in this case, set cos2σm to 0 and the result is computed correctly. He also points out that the formula may have no solution between two nearly antipodal points; an iteration limit traps this case (Vincenty says “this will occur when λ, as computed by eqn. (11), is greater than π in absolute value”, but this is not always a reliable test).
Note: some implementations of Vincenty’s formula inefficiently use a large number of trig functions; Vincenty devised this solution with an eye for efficiency in implementation, and this one uses just one each of sin, cos, sqrt, and atan2 for each iteration – only 3 or 4 iterations are generally required. [Formulation updated Dec 05 to make it closer to Vincenty’s original and computationally more efficient.]


The most accurate and widely used globally-applicable model for the earth ellipsoid is WGS-84, used in this script. Other ellipsoids offering a better fit to the local geoid include Airy (1830) in the UK, International 1924 in much of Europe, Clarke (1880) in Africa, GRS-67 in South America, and many others. America (NAD83) and Australia (GDA) use GRS-80, functionally equivalent to the WGS-84 ellipsoid.

  WGS-84 a = 6 378 137 m (±2 m) b ≈ 6 356 752.314245 m f ≈ 1 / 298.257223563
  GRS-80 a = 6 378 137 m b ≈ 6 356 752.314140 m f = 1 / 298.257222101
  Airy 1830 a = 6 377 563.396 m b = 6 356 256.909 m f ≈ 1 / 299.3249646
  Internat’l 1924 a = 6 378 388 m b ≈ 6 356 911.946 m f = 1 / 297
  Clarke mod.1880 a = 6 378 249.145 m b ≈ 6 356 514.86955 m f = 1 / 293.465
  GRS-67 a = 6 378 160 m b ≈ 6 356 774.719 m f = 1 / 298.247167

Note that to locate latitude/longitude points on these ellipses, they are associated with specific datums: for instance, OSGB36 for Airy in the UK, ED50 for Int’l 1924 in Europe; WGS-84 defines a datum as well as an ellipse. See my convert coordinates page for converting points between different datums.

Aside for fellow geeks: it is often said that WGS-84 and GRS-80 are ‘functionally equivalent’. Well, I would concur that a 0.1 mm difference in Earth radius is pretty equivalent, but why any difference? Apparently, when defining WGS-84, the normalized second degree zonal harmonic gravitational coefficient C2,0, derived from the GRS-80 value for the dynamic form factor J2, was truncated to 8 significant digits in the normalization process! It’s a question of which parameters are defining, and which are derived, in the reference frame.

Some of the terms involved are explained in Ed Williams’ notes on Spheroid Geometry.

Test case (from Geoscience Australia), using WGS-84:
  Flinders Peak 37°57′03.72030″S, 144°25′29.52440″E
  Buninyong 37°39′10.15610″S, 143°55′35.38390″E
  s 54 972.271 m
  α1 306°52′05.37″
  α2 127°10′25.07″ (≡ 307°10′25.07″ p1→p2)

Notes:


See below for the source code of the JavaScript implementation. These functions should be simple to translate into other languages if required.

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 link to this page.

Paypal donation If you would like to show your appreciation and support continued development of these scripts, I would most gratefully accept donations.

If you have any queries or find any problems, contact me at ku.oc.epyt-elbavom@oeg-stpircs.

© 2002-2012 Chris Veness



/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
/* Vincenty Inverse Solution of Geodesics on the Ellipsoid (c) Chris Veness 2002-2012             */
/*                                                                                                */
/* from: Vincenty inverse formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the */
/*       Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975    */
/*       http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf                                             */
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

/**
 * Calculates geodetic distance between two points specified by latitude/longitude using 
 * Vincenty inverse formula for ellipsoids
 *
 * @param   {Number} lat1, lon1: first point in decimal degrees
 * @param   {Number} lat2, lon2: second point in decimal degrees
 * @returns (Number} distance in metres between points
 */
function distVincenty(lat1, lon1, lat2, lon2) {
  var a = 6378137, b = 6356752.314245,  f = 1/298.257223563;  // WGS-84 ellipsoid params
  var L = (lon2-lon1).toRad();
  var U1 = Math.atan((1-f) * Math.tan(lat1.toRad()));
  var U2 = Math.atan((1-f) * Math.tan(lat2.toRad()));
  var sinU1 = Math.sin(U1), cosU1 = Math.cos(U1);
  var sinU2 = Math.sin(U2), cosU2 = Math.cos(U2);
  
  var lambda = L, lambdaP, iterLimit = 100;
  do {
    var sinLambda = Math.sin(lambda), cosLambda = Math.cos(lambda);
    var sinSigma = Math.sqrt((cosU2*sinLambda) * (cosU2*sinLambda) + 
      (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda));
    if (sinSigma==0) return 0;  // co-incident points
    var cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda;
    var sigma = Math.atan2(sinSigma, cosSigma);
    var sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
    var cosSqAlpha = 1 - sinAlpha*sinAlpha;
    var cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha;
    if (isNaN(cos2SigmaM)) cos2SigmaM = 0;  // equatorial line: cosSqAlpha=0 (§6)
    var C = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha));
    lambdaP = lambda;
    lambda = L + (1-C) * f * sinAlpha *
      (sigma + C*sinSigma*(cos2SigmaM+C*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)));
  } while (Math.abs(lambda-lambdaP) > 1e-12 && --iterLimit>0);

  if (iterLimit==0) return NaN  // formula failed to converge

  var uSq = cosSqAlpha * (a*a - b*b) / (b*b);
  var A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
  var B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
  var deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-
    B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
  var s = b*A*(sigma-deltaSigma);
  
  s = s.toFixed(3); // round to 1mm precision
  return s;
  
  // note: to return initial/final bearings in addition to distance, use something like:
  var fwdAz = Math.atan2(cosU2*sinLambda,  cosU1*sinU2-sinU1*cosU2*cosLambda);
  var revAz = Math.atan2(cosU1*sinLambda, -sinU1*cosU2+cosU1*sinU2*cosLambda);
  return { distance: s, initialBearing: fwdAz.toDeg(), finalBearing: revAz.toDeg() };
}

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */
  

Usage:


<p>Lat 1: <input name="lat1" value="53 09 02N"> Long 1: <input name="long1" value="001 50 40W"></p>
<p>Lat 2: <input name="lat2" value="52 12 19N"> Long 2: <input name="long2" value="000 08 33W"></p>

<input type="button" value="calculate distance" 
    onClick="result.value = distVincenty(Geo.parseDMS(lat1.value), Geo.parseDMS(long1.value), 
                                         Geo.parseDMS(lat2.value), Geo.parseDMS(long2.value)) + ' m'">
<input name="result" value="">
  

See the Lat/Long page for the deg/min/sec parsing method Geo.parseDMS()