This FAQ is copied from www.faqs.org/faqs/geography/infosystemsfaq
since the original at www.census.gov/cgibin/geo/gisfaq?Q5.1
has gone offline. I originally used it in researching formulae for calculating
distance between two Latitude/Longitude points.
Q5.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 flatEarth 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 flatEarth
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 illconditioned (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 flatEarth 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 FlatEarth 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 illconditioned. 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".
5.1a: What value should I use for the radius of the Earth, R?
The historical definition of a "nautical mile" is "one minute
of arc of a great circle of the earth". Since the earth is not a perfect
sphere, that definition is ambiguous. However, the internationally accepted
(SI) value for the length of a nautical mile is (exactly, by definition) 1.852
km or exactly 1.852/1.609344 international miles (that is, approximately 1.15078
miles  either "international" or "U.S. statute"). Thus,
the implied "official" circumference is 360 degrees times 60 minutes/degree
times 1.852 km/minute = 40003.2 km. The implied radius is the circumference
divided by 2 pi:
R = 6367 km = 3956 mi
The shape of the Earth is well approximated by an oblate spheroid with a polar
radius of 6357 km and an equatorial radius of 6378 km. PROVIDED a spherical
approximation is satisfactory, any value in that range will do, such as
R (in km) = 6378  21 * sin(lat) See the WARNING below!
R (in mi) = 3963  13 * sin(lat)
where lat is a latitude near which the bulk of your calculations occur.
WARNING: This formula for R gives but a rough approximation to the radius
of curvature as a function of latitude. The radius of curvature varies with
direction and latitude; according to Snyder ("Map Projections  A Working
Manual", by John P. Snyder, U.S. Geological Survey Professional Paper 1395,
United States Government Printing Office, Washington DC, 1987, p24), in the
plane of the meridian it is given by
R' = a * (1  e^2) / (1  e^2 * sin^2(lat))^(3/2)
where a is the equatorial radius, b is the polar radius, and e is the eccentricity
of the ellipsoid = (1  b^2/a^2)^(1/2).
5.1b: When is it NOT okay to assume the Earth is a sphere?
A quick test is: Compute the values of R produced by the equation with the
WARNING when you use the highest and lowest latitudes that occur in your analysis.
Compare the results produced by using these two values in your analysis. If
the different results are different enough to cause you to change your action
(or your recommendation, or your interpretation of the implication of the results,
etc.), then assuming the Earth is spherical is NOT okay.
For most purposes, it is quite satisfactory to treat the Earth as a sphere.
If necessary, an ellipsoid can provide a better approximation. Some standard
textbooks that may be helpful follow (reviews are by Steve Robertson of Tangent
Survey Systems in Canada: stever@mindlink.bc.ca):
Bomford, Guy 1980 Geodesy Clarendon Press, Oxford
ISBN 019851946X
Review: For geodetic computations, this is pretty well the standard in
English. It's a cookbook and offers no development, however.
Vanicek, Petr, and Krakiwsky, Edward 1986 Geodesy, the Concepts NorthHolland,
Amsterdam
ISBN 0444877754
Review: This offers a great, but quite involved, discussion of the concepts
behind geometrical (and all other) geodesy.
Torge, Wolfgang 1980 Geodesy de Gruyter, Berlin (translated to English
by C. Jekeli)
ISBN 3110072327
Review: This concentrates mostly on gravimetric geodesy, but has some
geometrical stuff, well explained without too much mathematics.
Software for solving distance and azimuth problems on the ellipsoid can be
obtained (as of 10/10/96) by anonymous ftp from several sources, two of which
are listed below:
The URL of the National Geodetic Survey (of the National Oceanic and Atmospheric
Administration in the US Department of Commerce) is:
ftp://www.ngs.noaa.gov/pub/pcsoft/for_inv.3d/
Review (by Ronald C. McConnell of Bellcore: rcmcc@cc.bellcore.com): They
have Fortran source and PC executable versions of both the normal "inverse"
great circle calculations (two lat/long pairs to distance and bearing), and
the less used "forward" calculation (one lat/long pair plus bearing
and distance to the second lat/long pair). They have both 2dimensional and
3dimensional versions of each. The inverse program works to within a few
seconds or a few minutes, depending on the fortran compiler, of the antipodal
points. The forward program seems immune to any and all problem locations
and pairs of locations. You can choose among a couple of dozen ellipsoids.
See the read.me file for explanations. The NGS software directory may contain
other listing of interest. Its URL is:
http://www.ngs.noaa.gov/PC_PROD/pc_prod.html/
Case is relevant in many URLs  eg: this one.
Another anonymous ftp source for ellipsoid software is the US Geological Survey
(of the US Department of the Interior), at:
http://kai.er.usgs.gov/ftp/PROJ.4/proj.html
Again, see the README file for explanations. The URLs for the USGS directory
and home page are:
http://kai.er.usgs.gov/ftp/index.html
http://kai.er.usgs.gov/homepage.html
5.1c: When is it NOT okay to assume the Earth is an ellipsoid?
The shape the Earth would assume if it were all measured at mean sea level
is called the geoid. The geoid varies no more than about a hundred meters above
or below a wellfitting ellipsoid, a variation far less than the ellipsoid varies
from the sphere. Terrain relief is reported relative to the geoid. (Paraphrased
from p. 11 of the book by Snyder cited above.)
Distances on the surface of the geoid are not particularly meaningful. However,
there are applications, such as longterm prediction of orbits of Earth satellites,
that require better approximations than are provided by an ellipsoid. Astrodynamics
texts, such as
Kaula, William M. 1966 Theory of Satellite Geodesy Blaisdell
Publishing Co., Waltham, MA (This book may be out of print.)
Battin, Richard H. 1964 Astronautical Guidance McGrawHill,
New York (There may be later editions.)
may be consulted for further information.
Distances on the surface of the terrain, whether geodesic, on roads, crosscountry,
or straightline, depend on relief (including elevation differences), the status
of engineering projects, and perhaps even route selection. Hence, computation
is idiosyncratic and not well suited to simple approximations.
