Version 40 (modified by karney, 4 years ago)

Correct title of Krueger's thesis

  1. Geodestic Calculations
  2. Mailing list threads
  3. Terminology
  4. The Math
    1. Spheroidal Approximation
    2. Ellipsoidal Approximation
      1. Thaddeus Vincenty's method, April 1975
      2. PROJ.4 - geod program
      3. GeographicLib::Geodesic
  5. The History

Geodestic Calculations

Geodestic calculations are calculations along lines (great circle) on the surface of the earth. They can answer questions like:

  • What is the distance between these two points?
  • If I travel X meters from point A at bearing phi, where will I be.

They are done in native lat-long coordinates, rather than in projected coordinates.

Mailing list threads

Terminology

The shortest distance on the surface of a solid is generally termed a geodesic, be it an ellipsoid of revolution, aposphere, etc. On a sphere, the geodesic is termed a Great Circle.

HOWEVER, when computing the distance between two points using a projected coordinate system, that is a conformal projection such as Transverse Mercator, Oblique Mercator, Normal Mercator, Stereographic, or Lambert Conformal Conic - that then is a GRID distance which can be converted to an equivalent GEODETIC distance using the function for "Scale Factor at a Point." The conversion is then termed "Grid Distance to Geodetic Distance," even though it will not be as exactly correct as a true ellipsoidal geodesic. Closer to the truth with a TM than with a Lambert or other conformal projection, but still not exactly "on."

So, it can be termed "geodetic distance" or a "geodesic distance," depending on just how you got there ...

The Math

Spheroidal Approximation

The simplest way to compute geodesics is using a sphere as an approximation for the earth. This from Mikael Rittri on the Proj mailing list:

If 1 percent accuracy is enough, I think you can use spherical formulas with a fixed Earth radius. You can find good formulas in the Aviation Formulary of Ed Williams,  http://williams.best.vwh.net/avform.htm.


For the fixed Earth radius, I would choose the average of the:
c = radius of curvature at the poles,
b2 / a = radius of curvature in a meridian plane at the equator,
since these are the extreme values for the local radius of curvature of the earth ellipsoid.


If your coordinates are given in WGS84, then
c = 6 399 593.626 m,
b2 / a = 6 335 439.327 m,
(see  http://home.online.no/~sigurdhu/WGS84_Eng.html) so their average is 6,367,516.477 m. The maximal error for distance calculation should then be less than 0.51 percent.


When computing the azimuth between two points by the spherical formulas, I think the maximal error on WGS84 will be 0.2 degrees, at least if the points are not too far away (less than 1000 km apart, say). The error should be maximal near the equator, for azimuths near northeast etc.


I am not sure about the spherical errors for the forward geodetic problem: point positioning given initial point, distance and azimuth.

Ellipsoidal Approximation

For more accuracy, the earth can be approximated with an ellipsoid, complicating the math somewhat.

Thaddeus Vincenty's method, April 1975

For a very good procedure to calculate inter point distances see:

 http://www.ngs.noaa.gov/PC_PROD/Inv_Fwd/ (Fortan code, DOS executables, and an online app)

and algorithm details published in:

 Vincenty, T. (1975)

Javascript code

Chris Veness has coded Vincenty's formulas as JavaScript.

distance:  http://www.movable-type.co.uk/scripts/latlong-vincenty.html
direct:  http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html

C code

From Gerald Evenden: a library of the converted NGS Vincenty geodesic procedure and an application program, 'geodesic'. In the case of a spherical earth Snyder's preferred equations are used.

Earlier Mr. Evenden had posted to the PROJ.4 mailing list this code for determination of true distance and respective forward and back azimuths between two points on the ellipsoid. Good for any pair of points that are not antipodal. Later he posted that this was not in fact the translation of NGS FORTRAN code, but something else. But, for what it's worth, here is the posted code (source unknown):

PROJ.4 - geod program

The PROJ.4 geod program can be used for great circle distances on an ellipsoid. Currently the underlying geodesic calculation API is not exposed as part of the PROJ.4 library. Gerald writes that geod is based upon a poorer algorithm (than the NGS procedure above) and he no longer supports it. But it is probably good enough for most applications. The method is documented here:

Paul D. Thomas, 1970
"Spheroidal Geodesics, Reference Systems, and Local Geometry"
U.S. Naval Oceanographic Office, p. 162
Engineering Library 526.3 T36s

 http://stinet.dtic.mil/oai/oai?&verb=getRecord&metadataPrefix=html&identifier=AD0703541

GeographicLib::Geodesic

Charles Karney has written a C++ class to do geodesic calculations and a utility Geod to call it. See

This is an attempt to do geodesic calculations "right", i.e.,

  • Accurate to round-off (i.e., about 15 nm);
  • Inverse solution always succeeds (even for near anti-podal points).

In addition, this class computes Christoffel's "reduced length", which gives the azimuthal scale for the corresponding azimuthal equidistant projection. This class is largely based on the work of Bessel (1826) and Helmert (1880) with the series approximations extended with Maxima.

The History

Here is a list of the older mathematical treatments of the geodesic problem for an ellipsoid, together with links to online copies. Unfortunately, the fold-out pages of figures in some books are usually not scanned properly by Google. I (Charles Karney) will add links to scans of the missing figures I can get hold of. Please let me know of errors, omissions, etc.