For the benefit of the terminally obsessive (as well as the genuinely needy, of course), Thaddeus Vincenty (‘TV’) 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.5mm 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 most purposes, of course).
Vincenty wrote this at a time when computing resources were expensive, and made his solutions very computationally efficient – typically just a dozen or so trig functions for the distance calculation.
Since Vincenty published this method, it has been improved on, particularly by Charles Karney (Algorithms for geodesics, 2012). I will add Karney’s method when I have a chance.
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 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. Pythagorean trigonometric identities sin²φ + cos²φ = 1 and 1 + tan²φ = sec²φ are used to reduce the number of trig functions invoked.
var L = λ2 - λ1; var tanU1 = (1-f) * Math.tan(φ1), cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1 * cosU1; var tanU2 = (1-f) * Math.tan(φ2), cosU2 = 1 / Math.sqrt((1 + tanU2*tanU2)), sinU2 = tanU2 * cosU2; var λ = L, λʹ, iterationLimit = 100; do { var sinλ = Math.sin(λ), cosλ = Math.cos(λ); var sinSqσ = (cosU2*sinλ) * (cosU2*sinλ) + (cosU1*sinU2-sinU1*cosU2*cosλ) * (cosU1*sinU2-sinU1*cosU2*cosλ); var sinσ = Math.sqrt(sinSqσ); if (sinσ==0) return 0; // co-incident points var cosσ = sinU1*sinU2 + cosU1*cosU2*cosλ; var σ = Math.atan2(sinσ, cosσ); var sinα = cosU1 * cosU2 * sinλ / sinσ; var cosSqα = 1 - sinα*sinα; var cos2σM = cosσ - 2*sinU1*sinU2/cosSqα; if (isNaN(cos2σM)) cos2σM = 0; // equatorial line: cosSqα=0 (§6) var 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'); var uSq = cosSqα * (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 Δσ = 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))); var s = b*A*(σ-Δσ); var fwdAz = Math.atan2(cosU2*sinλ, cosU1*sinU2-sinU1*cosU2*cosλ); var 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 | |
α_{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:
var sinα1 = Math.sin(α1); var cosα1 = Math.cos(α1); var tanU1 = (1-f) * Math.tan(φ1), cosU1 = 1 / Math.sqrt((1 + tanU1*tanU1)), sinU1 = tanU1 * cosU1; var σ1 = Math.atan2(tanU1, cosα1); var sinα = cosU1 * sinα1; var cosSqα = 1 - sinα*sinα; var uSq = cosSqα * (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 σ = s / (b*A), σʹ; do { var cos2σM = Math.cos(2*σ1 + σ); var sinσ = Math.sin(σ); var cosσ = Math.cos(σ); var Δσ = 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); var tmp = sinU1*sinσ - cosU1*cosσ*cosα1; var φ2 = Math.atan2(sinU1*cosσ + cosU1*sinσ*cosα1, (1-f)*Math.sqrt(sinα*sinα + tmp*tmp)); var λ = Math.atan2(sinσ*sinα1, cosU1*cosσ - sinU1*sinσ*cosα1); var C = f/16*cosSqα*(4+f*(4-3*cosSqα)); var L = λ - (1-C) * f * sinα * (σ + C*sinσ*(cos2σM+C*cosσ*(-1+2*cos2σM*cos2σM))); var λ2 = (λ1+L+3*Math.PI)%(2*Math.PI) - Math.PI; // normalise to -180...+180 var revAz = Math.atan2(sinα, -tmp);
See below for the complete code.
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 C_{2,0}, derived from the GRS-80 value for the dynamic form factor J_{2}, 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″ p_{1}→p_{2}) |
Notes:
... 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:
See below for the JavaScript source code, also available on GitHub. 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.
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.
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.
May 2014: I have organised the geodesy scripts into a collection of JavaScript objects, some complementary, some inter-dependent. If you are interested only in a certain aspect, I hope you should find it simple to extract elements of interest (eg if you want Vincenty’s formulae only on WGS84, you can add a LatLonE constructor to remove the dependency on latlon-ellipsoid.js, and hardwire the WGS84 ellipsoid parameters).
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-2014 Chris Veness