This article is from the Geographic Information Systems FAQ, by Lisa Nyman lnyman@census.gov with numerous contributions by others.

(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: