Movable Type Home Page

Movable Type Scripts


Vincenty solutions of geodesics on the ellipsoid

Distances and bearings between points on an ellipsoidal earth model

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 computa­tionally 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.

Distance/bearing between two points (inverse solution)

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 (ab)/a
φ1, φ2 = geodetic latitude
L = difference in longitude
tan U1/2 = (1−f) · tan φ1/2 (U is ‘reduced latitude’)
cos U1/2 = 1 / √1 + tan² U1/2, sin U1/2 = tan U1/2 · cos U1/2 (trig identities; §6)
λ = L (first approximation)
iterate until change in λ is negligible (e.g. 10-12 ≈ 0.006mm) {
sin σ = √[ (cos U2 · sinλ)² + (cos U1 · sin U2 − sin U1 · cos U2 · cos λ)² ] (14)
cos σ = sin U1 · sin U2 + cos U1 · cos U2 · cos λ (15)
σ = atan(sin σ / cos σ) (16)
sin α = cos U1 · cos U2 · sin λ / sin σ (17)
cos² α = 1 − sin² α (trig identity; §6)
cos 2σm = cos σ − 2 · sin U1 · sin U2 / 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 U2 · sin λ / cos U1 · sin U2 − sin U1 · cos U2 · cos λ) (20)
α2 = atan(cos U1 · sin λ / −sin U1 · cos U2 + cos U1 · sin U2 · 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.

JavaScript

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.

Destination given distance & bearing from start point (direct solution)

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 (ab)/a
φ1, φ2 = geodetic latitude
L = difference in longitude
s = length of the geodesic
α1, α2 = azimuths of the geodesic (initial/final bearing)
 
tan U1 = (1−f) · tan φ1 (U is ‘reduced latitude’)
cos U1 = 1 / √1 + tan² U1, sin U1 = tan U1 · cos U1 (trig identities; §6)
σ1 = atan(tan U1 / cos α1) (1)
sin α = cos U1 · 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 U1 · cos σ + cos U1 · sin σ · cos α1
     / (1−f) · √sin² α + (sin U1 · sin σ − cos U1 · cos σ · cos α1 )
(8)
λ = atan(sin σ · sin α1 / cos U1 · cos σ − sin U1 · 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 U1 · sin σ − cos U1 · cos σ · cos α1) ) (12)

Where:

JavaScript

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.

Ellipsoids

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:


Live examples

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

Inverse solution
Point 1: ,
Point 2: ,
Distance: m
Initial bearing:
Final bearing:
Direct solution
Start point: ,
Bearing:
Distance: m
Destination point:
Final bearing:

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

OSI MIT License 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 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