5.1: What is the best way to calculate the great circle distance (which deliberately ignores elevation differences) between 2 points?
Description
This article is from the Geographic Information Systems FAQ, by Lisa Nyman lnyman@census.gov with numerous contributions by others.
5.1: What is the best way to calculate the great circle distance (which deliberately ignores elevation differences) between 2 points?
(This answer was prepared by Robert G. Chamberlain of Caltech (JPL):
rgc@solstice.jpl.nasa.gov and reviewed on the comp.infosystems.gis
newsgroup in Oct 1996.)
If the distance is less than about 20 km (12 mi) and the locations of the
two points in Cartesian coordinates are X1,Y1 and X2,Y2 then the
Pythagorean Theorem
d = sqrt((X2 - X1)^2 + (Y2 - Y1)^2)
will result in an error of
less than 30 meters (100 ft) for latitudes less than 70 degrees
less than 20 meters ( 66 ft) for latitudes less than 50 degrees
less than 9 meters ( 30 ft) for latitudes less than 30 degrees
(These error statements reflect both the convergence of
the meridians and the curvature of the parallels.)
The flat-Earth distance d will be expressed in the same units as
the coordinates.
If the locations are not already in Cartesian coordinates, the
computational cost of converting from spherical coordinates and
then using the flat-Earth model may exceed that of using the
more accurate spherical model.
Otherwise, presuming a spherical Earth with radius R (see below), and the
locations of the two points in spherical coordinates (longitude and
latitude) are lon1,lat1 and lon2,lat2 then the
Haversine Formula (from R.W. Sinnott, "Virtues of the Haversine",
Sky and Telescope, vol. 68, no. 2, 1984, p. 159):
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin^2(dlat/2) + cos(lat1) * cos(lat2) * sin^2(dlon/2)
c = 2 * arcsin(min(1,sqrt(a)))
d = R * c
will give mathematically and computationally exact results. The
intermediate result c is the great circle distance in radians.
The great circle distance d will be in the same units as R.
The min() function protects against possible roundoff errors that
could sabotage computation of the arcsine if the two points are
very nearly antipodal (that is, on opposide sides of the Earth).
Under these conditions, the Haversine Formula is ill-conditioned
(see the discussion below), but the error, perhaps as large as
2 km (1 mi), is in the context of a distance near 20,000 km
(12,000 mi).
Most computers require the arguments of trignometric functions to
be expressed in radians. To convert lon1,lat1 and lon2,lat2 from
degrees, minutes, and seconds to radians, first convert them to
decimal degrees. To convert decimal degrees to radians, multiply
the number of degrees by pi/180 = 0.017453293 radians/degree.
Inverse trigonometric functions return results expressed in
radians. To express c in decimal degrees, multiply the number of
radians by 180/pi = 57.295780 degrees/radian. (But be sure to
multiply the number of RADIANS by R to get d.)
The problem of determining the great circle distance on a sphere
has been around for hundreds of years, as have both the Law of
Cosines solution (given below but not recommended) and the
Haversine Formula. Sinnott gets the credit here because he was
quoted by Snyder (cited below). Perhaps someone will provide the
truly seminal reference so the proper attribution can be given?
The Pythagorean flat-Earth approximation assumes that meridians are
parallel, that the parallels of latitude are negligibly different from
great circles, and that great circles are negligibly different from
straight lines. Close to the poles, the parallels of latitude are not only
shorter than great circles, but indispensably curved. Taking this into
account leads to the use of polar coordinates and the planar law of cosines
for computing short distances near the poles: The
Polar Coordinate Flat-Earth Formula
a = pi/2 - lat1
b = pi/2 - lat2
c = sqrt(a^2 + b^2 - 2 * a * b * cos(lon2 - lon1)
d = R * c
will give smaller maximum errors than the Pythagorean Theorem for
higher latitudes and greater distances. (The maximum errors, which
depend upon azimuth in addition to separation distance, are equal
at 80 degrees latitude when the separation is 33 km (20 mi),
82 degrees at 18 km (11 mi), 84 degrees at 9 km (5.4 mi).) But
even at 88 degrees the polar error can be as large as 20 meters
(66 ft) when the distance between the points is 20 km (12 mi).
The latitudes lat1 and lat2 must be expressed in radians (see
above); pi/2 = 1.5707963. Again, the intermediate result c is the
distance in radians and the distance d is in the same units as R.
An UNRELIABLE way to calculate distance on a spherical Earth is the
Law of Cosines for Spherical Trigonometry
** NOT RECOMMENDED **
a = sin(lat1) * sin(lat2)
b = cos(lat1) * cos(lat2) * cos(lon2 - lon1)
c = arccos(a + b)
d = R * c
Although this formula is mathematically exact, it is unreliable
for small distances because the inverse cosine is ill-conditioned.
Sinnott (in the article cited above) offers the following table
to illustrate the point:
cos (5 degrees) = 0.996194698
cos (1 degree) = 0.999847695
cos (1 minute) = 0.9999999577
cos (1 second) = 0.9999999999882
cos (0.05 sec) = 0.999999999999971
A computer carrying seven significant figures cannot distinguish
the cosines of any distances smaller than about one minute of arc.
The function min(1,(a + b)) could replace (a + b) as the argument
for the inverse cosine for the same reason as in Sinnott's Formula,
but doing so would "polish a cannonball".
Continue to:
My Books
