For the benefit of the terminally obsessive (as well as the genuinely needy, of course), Thaddeus Vincenty devised formulæ 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 solution for the distance between points on an ellipsoidal earth model is accurate to within 0.5 mm distance (!), 0.000015″ bearing, on the ellipsoid being used. Calculations based on a spherical earth model, such as the (much simpler) Haversine, are accurate to around 0.3% – which is still good enough for many (most?) purposes, of course.
Vincenty’s inverse solution can fail on nearly antipodal points. Testing with GeographicLib test data, I’ve found this can happen with distances greater than 19,936 km, or within around 75 km of the antipodal point.
Vincenty wrote this at a time when computing resources were expensive, and made his solutions very computationally efficient – the distance calculation takes around 5 – 10 microseconds using Chrome on a middling Core i5 PC.
Vincenty’s paper presents a set of formulæ rather than a single program; this shows those formulæ pulled together as they are used in the script:
a, b = major & minor semi-axes of the ellipsoid | |
f = flattening (a−b)/a | |
φ_{1}, φ_{2} = geodetic latitude | |
L = difference in longitude_{ } | |
tan U_{1/2} = (1−f) · tan φ_{1/2} | (U is ‘reduced latitude’) |
cos U_{1/2} = 1 / √1 + tan² U_{1/2}, sin U_{1/2} = tan U_{1/2} · cos U_{1/2} | (trig identities; §6) |
λ = L (first approximation) | |
iterate until change in λ is negligible (e.g. 10^{-12} ≈ 0.006mm) { | |
sin σ = √[ (cos U_{2} · sinλ)² + (cos U_{1} · sin U_{2} − sin U_{1} · cos U_{2} · cos λ)² ] | (14) |
cos σ = sin U_{1} · sin U_{2} + cos U_{1} · cos U_{2} · cos λ | (15) |
σ = atan(sin σ / cos σ)_{ } | (16) |
sin α = cos U_{1} · cos U_{2} · sin λ / sin σ | (17) |
cos² α = 1 − sin² α | (trig identity; §6) |
cos 2σ_{m} = cos σ − 2 · sin U_{1} · sin U_{2} / cos² α | (18) |
C = f/16 · cos² α · [4 + f · (4 − 3 · cos² α)]_{ } | (10) |
λ′ = L + (1 − C) · f · sin α · {σ + C · sin σ · [cos 2σ_{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 σ · {cos 2σ_{m} + B/4 · [cos σ · (−1 + 2 · cos² 2σ_{m}) − B / 6 · cos 2σ_{m} · (−3 + 4 · sin² σ) · (−3 + 4 · cos² 2σ_{m})]} |
(6) |
s = b·A·(σ−Δσ) | (19) |
α_{1} = atan(cos U_{2} · sin λ / cos U_{1} · sin U_{2} − sin U_{1} · cos U_{2} · cos λ) | (20) |
α_{2} = atan(cos U_{1} · sin λ / −sin U_{1} · cos U_{2} + cos U_{1} · sin U_{2} · 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 in a small number of cases, the formula will not converge even when λ stays below π).
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. Pythagorean trigonometric identities sin²φ + cos²φ = 1 and 1 + tan²φ = sec²φ are used to reduce the number of trig functions invoked.
const L = λ2 - λ1; const tanU1 = (1-f) * Math.tan(φ1), cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1 * cosU1; const tanU2 = (1-f) * Math.tan(φ2), cosU2 = 1 / Math.sqrt((1 + tanU2*tanU2)), sinU2 = tanU2 * cosU2; const λ = L, λʹ, iterationLimit = 100; do { const sinλ = Math.sin(λ), cosλ = Math.cos(λ); const sinSqσ = (cosU2*sinλ) * (cosU2*sinλ) + (cosU1*sinU2-sinU1*cosU2*cosλ) * (cosU1*sinU2-sinU1*cosU2*cosλ); const sinσ = Math.sqrt(sinSqσ); if (sinσ==0) return 0; // co-incident points const cosσ = sinU1*sinU2 + cosU1*cosU2*cosλ; const σ = Math.atan2(sinσ, cosσ); const sinα = cosU1 * cosU2 * sinλ / sinσ; const cosSqα = 1 - sinα*sinα; const cos2σM = cosσ - 2*sinU1*sinU2/cosSqα; if (isNaN(cos2σM)) cos2σM = 0; // equatorial line: cosSqα=0 (§6) const C = f/16*cosSqα*(4+f*(4-3*cosSqα)); λʹ = λ; λ = L + (1-C) * f * sinα * (σ + C*sinσ*(cos2σM+C*cosσ*(-1+2*cos2σM*cos2σM))); } while (Math.abs(λ-λʹ) > 1e-12 && --iterationLimit>0); if (iterationLimit==0) throw new Error('Formula failed to converge'); const uSq = cosSqα * (a*a - b*b) / (b*b); const A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq))); const B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq))); const Δσ = B*sinσ*(cos2σM+B/4*(cosσ*(-1+2*cos2σM*cos2σM)- B/6*cos2σM*(-3+4*sinσ*sinσ)*(-3+4*cos2σM*cos2σM))); const s = b*A*(σ-Δσ); const fwdAz = Math.atan2(cosU2*sinλ, cosU1*sinU2-sinU1*cosU2*cosλ); const revAz = Math.atan2(cosU1*sinλ, -sinU1*cosU2+cosU1*sinU2*cosλ);
See below for the complete code.
Again, this presents Vincenty’s formulæ arranged to be close to how they are used in the script.
a, b = major & minor semi-axes of the ellipsoid | |
f = flattening (a−b)/a | |
φ_{1}, φ_{2} = geodetic latitude | |
L = difference in longitude | |
s = length of the geodesic along the surface of the ellipsoid (in the same units as a & b) | |
α_{1}, α_{2} = azimuths of the geodesic (initial/final bearing) | |
tan U_{1} = (1−f) · tan φ_{1} | (U is ‘reduced latitude’) |
cos U_{1} = 1 / √1 + tan² U_{1}, sin U_{1} = tan U_{1} · cos U_{1} | (trig identities; §6) |
σ_{1} = atan(tan U_{1} / cos α_{1}) | (1) |
sin α = cos U_{1} · sin α_{1} | (2) |
cos² α = 1 − sin² α | (trig identity; §6) |
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) |
σ = s / (b·A) | (first approximation) |
iterate until change in σ is negligible (e.g. 10^{-12} ≈ 0.006mm) { | |
cos 2σ_{m} = cos(2σ_{1} + σ) | (5) |
Δσ = B · sin σ · {cos 2σ_{m} + B/4 · [cos σ · (−1 + 2 · cos² 2σ_{m}) − B/6 · cos 2σ_{m} · (−3 + 4 · sin² σ) · (−3 + 4 · cos² 2σ_{m})]} |
(6) |
σʹ = s / b·A + Δσ | (7) |
} | |
φ_{2} = atan(sin U_{1} · cos σ + cos U_{1} · sin σ · cos α_{1} / (1−f) · √sin² α + (sin U_{1} · sin σ − cos U_{1} · cos σ · cos α_{1})² ) |
(8) |
λ = atan(sin σ · sin α_{1} / cos U_{1} · cos σ − sin U_{1} · sin σ · cos α_{1}) | (9) |
C = f/16 · cos² α · [4 + f · (4 − 3 · cos² α)]_{ } | (10) |
L = λ − (1−C) · f · sin α · {σ + C · sin σ · [cos 2σ_{m} + C · cos σ · (−1 + 2 · cos² 2σ_{m})]} | (11) |
λ_{2} = λ_{1} + L | |
α_{2} = atan( sin α / −(sin U_{1} · sin σ − cos U_{1} · cos σ · cos α_{1}) ) | (12) |
Where:
const sinα1 = Math.sin(α1); const cosα1 = Math.cos(α1); const tanU1 = (1-f) * Math.tan(φ1), cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1 * cosU1; const σ1 = Math.atan2(tanU1, cosα1); const sinα = cosU1 * sinα1; const cosSqα = 1 - sinα*sinα; const uSq = cosSqα * (a*a - b*b) / (b*b); const A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq))); const B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq))); let σ = s / (b*A), σʹ; do { const cos2σM = Math.cos(2*σ1 + σ); const sinσ = Math.sin(σ); const cosσ = Math.cos(σ); const Δσ = B*sinσ*(cos2σM+B/4*(cosσ*(-1+2*cos2σM*cos2σM)- B/6*cos2σM*(-3+4*sinσ*sinσ)*(-3+4*cos2σM*cos2σM))); σʹ = σ; σ = s / (b*A) + Δσ; } while (Math.abs(σ-σʹ) > 1e-12); const tmp = sinU1*sinσ - cosU1*cosσ*cosα1; const φ2 = Math.atan2(sinU1*cosσ + cosU1*sinσ*cosα1, (1-f)*Math.sqrt(sinα*sinα + tmp*tmp)); const λ = Math.atan2(sinσ*sinα1, cosU1*cosσ - sinU1*sinσ*cosα1); const C = f/16*cosSqα*(4+f*(4-3*cosSqα)); const L = λ - (1-C) * f * sinα * (σ + C*sinσ*(cos2σM+C*cosσ*(-1+2*cos2σM*cos2σM))); const λ2 = (λ1+L+3*Math.PI)%(2*Math.PI) - Math.PI; // normalise to -180...+180 const revAz = Math.atan2(sinα, -tmp);
See below for the complete code.
While is it wonderful to be able to calculate distances to millimetre accuracy, there are certain caveats which it is worth being aware of if you are looking for accuracy better than around ±1 metre.
Millimetre precision is stretching the limits of geodetic survey techniques: very few points will have been surveyed to such precision (it’s easy to provide latitude/longitude references to 6 decimals of degrees; less easy to make that precision meaningful rather than spurious).
Further, the familiar WGS-84 (World Geodetic System 1984) has no ‘physical realisation’ – it is not tied to geodetic groundstations, just to satellites – and is defined to be accurate to no better than around ±1 metre, irrespective of how good your GPS unit is [actually, the ‘G-series’ are a little more accurate, but there are no official transformation parameters]. A metre is good enough for most purposes, but a long way from Vincenty’s millimetre accuracy.
Then, when working to an accuracy of better than ±1 metre, plate tectonic movements become significant. Depending on the datum it’s defined in, any given point will most likely have shifted many millimetres from where it was last year.
Simplifying somewhat (well, actually, a lot!), in order to address requirements for greater precision, the ITRS was developed to give greater accuracy than WGS84, with ‘epoch’-dependant ITRF ‘realisations’, where the latitude/longitude coordinate of a position will vary over time.
Various ‘static’ reference frames are also defined for various continents – NAD83 for North America, ETRS89 for Europe, GDA94 for Australia, etc – within which latitude/longitude coordinates remain essentially fixed (at least to centimetre or so accuracy, major earthquake events excepted). These reference frames have ‘epoch’-dependant mappings to ITRF datums. (Due to plate tectonic continental drift, ETRS89 shifts against ITRF by about 25mm/year; GDA94 by around 80mm/year).
And for any latitude/longitude coordinates dating from earlier datums (such as NAD27, ED50, NTF, AGD66, and any number of others), even metre precision is exceeding their limits.
Also, it must be remembered that these distances are calculated on the surface of the reference ellipsoid, which will differ (to varying degrees) from the earth geoid, and even more from the earth’s physical surface. The geoid is equivalent to the local mean-sea-level (or its equivalent in land-locked locales), and is determined by the local force of gravity; it varies in height from the ellipsoid by around ±100 metres. The earth’s surface varies much more; if you happen to be located in Colorado, 2 kilometres above msl, distances will be 0.03% greater than on the ellipsoid surface; +1 metre over just 30 kilometres.
So for millimetre accuracy to be relevant, you need to know a fair bit about the reference frames, datums, and epochs of your coordinates, and other parameters of your measurements. Caveat emptor!
The ubiquitous WGS-84 is a geocentric datum, based on an ellipsoid with:
Semi-major axis | a | = 6 378 137.0 | metres |
Semi-minor axis | b | ≈ 6 356 752.314 245 | metres |
Inverse flattening | 1/f | = 298.257 223 563 |
Latitude/longitude points can be converted between different datums by converting them to geocentric cartesian coordinates, applying a 7-parameter Helmert transformation, and converting back to a geodetic coordinate in the new datum.
Should even millimetre accuracy not be sufficient for you, Charles Karney has improved on Vincenty’s method with a method which has errors in nanometers (and always converges on antipodal points) – Geodesics on an ellipsoid of revolution, 2011, Algorithms for geodesics, 2012; though Karney’s method is a lot less straightforward (someday I may get around to an implementation).
... to verify the code is correct!
Enter the co-ordinates into the text boxes to try out the calculations. A variety of formats are accepted, principally:
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″ p_{1}→p_{2}) |
Notes:
See below for the JavaScript source code, also available on GitHub. Full documentation is available, as well as a test suite.
Note I use Greek letters in variables representing maths symbols conventionally presented as Greek letters:
I value the great benefit in legibility over the minor inconvenience in typing (if you encounter
any problems, ensure your <head>
includes <meta charset="utf-8">
),
and/or use UTF-8 encoding when saving files.
With its untyped C-style syntax, JavaScript reads remarkably close to pseudo-code: exposing the algorithms with a minimum of syntactic distractions. These functions should be simple to translate into other languages if required, though can also be used as-is in browsers and Node.js.
For convenience & clarity, I have extended the base JavaScript Number
object with
toRadians()
and toDegrees()
methods: I don’t see great likelihood of
conflicts, as these are ubiquitous operations.
February 2019: I have refactored the library to use ES modules, as well as extending it in scope; see the GitHub README and CHANGELOG for details.
Other languages: I cannot support translations in languages I’m not familiar with, but if you have ported the code to another language, I am happy to provide links here.
I offer these scripts for free use and adaptation to balance my debt to the open-source info-verse. You are welcome to re-use these scripts [under an MIT licence, 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-2020 Chris Veness