两经纬度间的距离——Vincenty公式

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 accurateellipsoidal 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 (a−b)/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:

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

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 scr