- Vincenty's formulae
Vincenty's formulae form a fast algorithm to calculate the distance between two latitude/longitude points on an ellipsoid. They have been widely used in
geodesy because they are accurate to within 0.5mm (0.000015″) on the Earth ellipsoid. It has been developed byThaddeus Vincenty in 1975.Algorithm
The following algorithm is an explicit and useful expression of the formulae: a, b = major & minor semiaxes of the ellipsoid φ1, φ2 = geodetic latitude L = difference in longitude f = flattening (a−b)/a U1 = atan((1−f).tanφ1) (U is ‘reduced latitude’) U2 = atan((1−f).tanφ2) λ = L, λ′ = 2π do while abs(λ−λ′) > 10-12 (i.e. 0.06mm error) { sinσ = √ [ (cosU2.sinλ)² + (cosU1.sinU2 − sinU1.cosU2.cosλ)² ] cosσ = sinU1.sinU2 + cosU1.cosU2.cosλ σ = atan2(sinσ, cosσ) sinα = cosU1.cosU2.sinλ / sinσ cos²α = 1 − sin²α (trig identity; §6) cos2σm = cosσ − 2.sinU1.sinU2/cos²α C = f/16.cos²α. [4+f.(4−3.cos²α)] λ′ = λ λ = L + (1−C).f.sinα.{σ+C.sinσ. [cos2σm+C.cosσ.(−1+2.cos²2σm)] } (use cos2σm=0 when over equatorial lines) } u² = cos²α.(a²−b²)/b² A = 1+u²/16384.{4096+u². [−768+u².(320−175.u²)] } B = u²/1024.{256+u². [−128+u².(74−47.u²)] } Δσ = B.sinσ.{cos2σm+B/4. [cosσ.(−1+2.cos²2σm) − B/6.cos2σm.(−3+4.sin²σ).(−3+4.cos²2σm)] } s = b.A.(σ−Δσ) α1 = atan2(cosU2.sinλ, cosU1.sinU2 − sinU1.cosU2.cosλ) α2 = atan2(cosU1.sinλ, −sinU1.cosU2 + cosU1.sinU2.cosλ)
Where: *s is the distance (in the same units as a & b) *α1 is the initial bearing, or forward azimuth *α2 is the final bearing (in direction p1→p2)
The formula may have no solution between two nearly antipodal points; an iteration limit avoids this case.
The most accurate and widely used earth ellipsoid is
WGS-84 , for which: a = 6 378 137 m (±2 m) b = 6 356 752.3142 m f = 1 / 298.257223563Test case
Test case from [http://www.ga.gov.au/geodesy/datums/vincenty_inverse.jsp 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 (spherical distance would yield 54916.730 m) α1 306°52′05.37″ α2 127°10′25.07″ (≡ 307°10′25.07″ p1→p2)
References
[http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf "Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations". Vincenty's original paper] .
External links
[http://www.movable-type.co.uk/scripts/latlong-vincenty.html Computer source code available here (LGPL license) by Chris Veness]
Wikimedia Foundation. 2010.
