Great-circle distance

From HandWiki
A diagram illustrating great-circle distance (drawn in red) between two points on a sphere, P and Q. Two antipodal points, u and v are also shown.

The great-circle distance, orthodromic distance, or spherical distance is the distance along a great circle.

It is the shortest distance between two points on the surface of a sphere, measured along the surface of the sphere (as opposed to a straight line through the sphere's interior). The distance between two points in Euclidean space is the length of a straight line between them, i.e., the chord's length, but on the sphere there are no straight lines. In spaces with curvature, straight lines are replaced by geodesics. Geodesics on the sphere are circles on the sphere whose centers coincide with the center of the sphere, and are called 'great circles'.

The determination of the great-circle distance is part of the more general problem of great-circle navigation, which also computes the azimuths at the end points and intermediate way-points.

Through any two points on a sphere that are not antipodal points (directly opposite each other), there is a unique great circle. The two points separate the great circle into two arcs. The length of the shorter arc is the great-circle distance between the points. A great circle endowed with such a distance is called a Riemannian circle in Riemannian geometry.

Between antipodal points, there are infinitely many great circles, and all great circle arcs between antipodal points have a length of half the circumference of the circle, or [math]\displaystyle{ \pi r }[/math], where r is the radius of the sphere.

The Earth is nearly spherical, so great-circle distance formulas give the distance between points on the surface of the Earth correct to within about 0.5%.[1]


An illustration of the central angle, Δσ, between two points, P and Q. λ and φ are the longitudinal and latitudinal angles of P respectively

Let [math]\displaystyle{ \lambda_1, \phi_1 }[/math] and [math]\displaystyle{ \lambda_2, \phi_2 }[/math] be the geographical longitude and latitude of two points 1 and 2, and [math]\displaystyle{ \Delta\lambda, \Delta\phi }[/math] be their absolute differences; then [math]\displaystyle{ \Delta\sigma }[/math], the central angle between them, is given by the spherical law of cosines if one of the poles is used as an auxiliary third point on the sphere:[2]

[math]\displaystyle{ \Delta\sigma = \arccos\bigl(\sin\phi_1\sin\phi_2 + \cos\phi_1\cos\phi_2\cos\Delta\lambda\bigr). }[/math]

The problem is normally expressed in terms of finding the central angle [math]\displaystyle{ \Delta\sigma }[/math]. Given this angle in radians, the actual arc length d on a sphere of radius r can be trivially computed as

[math]\displaystyle{ d = r \, \Delta\sigma. }[/math]

Relation between central angle and chord length

The central angle [math]\displaystyle{ \Delta\sigma }[/math] is related with the chord length of unit sphere [math]\displaystyle{ \Delta\sigma_\text{c}\,\! }[/math]:

[math]\displaystyle{ \begin{align} \Delta\sigma &= 2\arcsin \frac{\Delta\sigma_\text{c}}{2} ,\\ \Delta\sigma_\text{c} &= 2\sin \frac{\Delta\sigma}{2} . \end{align} }[/math]

Computational formulae

On computer systems with low floating point precision, the spherical law of cosines formula can have large rounding errors if the distance is small (if the two points are a kilometer apart on the surface of the Earth, the cosine of the central angle is near 0.99999999). For modern 64-bit floating-point numbers, the spherical law of cosines formula, given above, does not have serious rounding errors for distances larger than a few meters on the surface of the Earth.[3] The haversine formula is numerically better-conditioned for small distances by using the chord-length relation:[4]

[math]\displaystyle{ \begin{align} \Delta\sigma &= \operatorname{archav}\left( \operatorname{hav}\left(\Delta\phi\right) + \left(1 - \operatorname{hav}(\Delta\phi) - \operatorname{hav}(\phi_1 + \phi_2)\right) \operatorname{hav}\left(\Delta\lambda\right)\right) . \end{align} }[/math]

Historically, the use of this formula was simplified by the availability of tables for the haversine function: hav(θ) = sin2(θ/2).

The following shows the equivalent formula expressing the chord length explicitly:

[math]\displaystyle{ \begin{align} \Delta\sigma &= 2\arcsin \frac{\Delta\sigma_\text{c}}{2} ,\\ \Delta\sigma_\text{c} &= 2\sqrt{\sin^2\frac{\Delta\phi}{2} + \left(\cos^2\phi_\text{m} - \sin^2\frac{\Delta\phi}{2}\right) \sin^2\frac{\Delta\lambda}{2}}, \end{align} }[/math]

where [math]\displaystyle{ \phi_\text{m}=\frac{\phi_1+\phi_2}{2} }[/math].

Although this formula is accurate for most distances on a sphere, it too suffers from rounding errors for the special (and somewhat unusual) case of antipodal points. A formula that is accurate for all distances is the following special case of the Vincenty formula for an ellipsoid with equal major and minor axes:[5]

[math]\displaystyle{ \Delta\sigma = \arctan \frac{\sqrt{\left( \cos\phi_2\sin\Delta\lambda\right)^2 + \left(\cos\phi_1\sin\phi_2 - \sin\phi_1\cos\phi_2\cos\Delta\lambda \right)^2}} {\sin\phi_1\sin\phi_2 + \cos\phi_1\cos\phi_2\cos\Delta\lambda} . }[/math]

Here the quadrant for [math]\displaystyle{ \Delta\sigma }[/math] should be governed by the signs of the numerator and denominator of the right hand side, e.g., using the atan2 function.

Vector version

Another representation of similar formulas, but using normal vectors instead of latitude and longitude to describe the positions, is found by means of 3D vector algebra, using the dot product, cross product, or a combination:[6]

[math]\displaystyle{ \begin{align} \Delta\sigma &= \arccos \left(\mathbf n_1 \cdot \mathbf n_2\right) \\ &= \arcsin \left| \mathbf n_1 \times \mathbf n_2 \right| \\ &= \arctan \frac{\left| \mathbf n_1 \times \mathbf n_2 \right|}{\mathbf n_1 \cdot \mathbf n_2} \\ \end{align} }[/math]

where [math]\displaystyle{ \mathbf n_1 }[/math] and [math]\displaystyle{ \mathbf n_2 }[/math] are the normals to the sphere at the two positions 1 and 2. Similarly to the equations above based on latitude and longitude, the expression based on arctan is the only one that is well-conditioned for all angles. The expression based on arctan requires the magnitude of the cross product over the dot product.

From chord length

A line through three-dimensional space between points of interest on a spherical Earth is the chord of the great circle between the points. The central angle between the two points can be determined from the chord length. The great circle distance is proportional to the central angle.

The great circle chord length, [math]\displaystyle{ \Delta\sigma_\text{c}\,\! }[/math], may be calculated as follows for the corresponding unit sphere, by means of Cartesian subtraction:

[math]\displaystyle{ \begin{align} \Delta{X} &= \cos\phi_2\cos\lambda_2 - \cos\phi_1\cos\lambda_1;\\ \Delta{Y} &= \cos\phi_2\sin\lambda_2 - \cos\phi_1\sin\lambda_1;\\ \Delta{Z} &= \sin\phi_2 - \sin\phi_1;\\ \Delta\sigma_\text{c} &= \sqrt{(\Delta{X})^2 + (\Delta{Y})^2 + (\Delta{Z})^2}. \end{align} }[/math]

Radius for spherical Earth

Equatorial (a), polar (b) and mean Earth radii as defined in the 1984 World Geodetic System revision. (Not to scale.)
Main page: Astronomy:Earth radius

The shape of the Earth closely resembles a flattened sphere (a spheroid) with equatorial radius [math]\displaystyle{ a }[/math] of 6378.137 km; distance [math]\displaystyle{ b }[/math] from the center of the spheroid to each pole is 6356.7523142 km. When calculating the length of a short north-south line at the equator, the circle that best approximates that line has a radius of [math]\displaystyle{ \frac{b^2}{a} }[/math] (which equals the meridian's semi-latus rectum), or 6335.439 km, while the spheroid at the poles is best approximated by a sphere of radius [math]\displaystyle{ \frac{a^2}{b} }[/math], or 6399.594 km, a 1% difference. So long as a spherical Earth is assumed, any single formula for distance on the Earth is only guaranteed correct within 0.5% (though better accuracy is possible if the formula is only intended to apply to a limited area). Using the mean Earth radius, [math]\displaystyle{ R_1 = \frac{1}{3}(2a + b) \approx 6371.009\text{ km} }[/math] (for the WGS84 ellipsoid) means that in the limit of small flattening, the mean square relative error in the estimates for distance is minimized.[7]

See also

References and notes

  1. Admiralty Manual of Navigation, Volume 1, The Stationery Office, 1987, p. 10, ISBN 9780117728806, 
  2. Kells, Lyman M.; Kern, Willis F.; Bland, James R. (1940). Plane And Spherical Trigonometry. McGraw Hill Book Company, Inc.. pp. 323-326. Retrieved July 13, 2018. 
  3. "Calculate distance, bearing and more between Latitude/Longitude points". 
  4. Sinnott, Roger W. (August 1984). "Virtues of the Haversine". Sky and Telescope 68 (2): 159. 
  5. Vincenty, Thaddeus (1975-04-01). "Direct and Inverse Solutions of Geodesics on the Ellipsoid with Application of Nested Equations". Survey Review (Kingston Road, Tolworth, Surrey: Directorate of Overseas Surveys) 23 (176): 88–93. doi:10.1179/sre.1975.23.176.88. Retrieved 2008-07-21. 
  6. Gade, Kenneth (2010). "A non-singular horizontal position representation". The Journal of Navigation (Cambridge University Press) 63 (3): 395–417. doi:10.1017/S0373463309990415. 
  7. McCaw, G. T. (1932). "Long lines on the Earth". Empire Survey Review 1 (6): 259–263. doi:10.1179/sre.1932.1.6.259. 

External links