This FAQ is copied from www.faqs.org/faqs/geography/infosystems-faq since the original at www.census.gov/cgi-bin/geo/gisfaq?Q5.1 has gone off-line. I originally used it in researching formulae for calculating distance between two Latitude/Longitude 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 * cwill 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 * cwill 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 * cAlthough 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".
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).
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 0-19-851946-XReview: 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 North-Holland, Amsterdam
ISBN 0-444-87775-4Review: 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 3-11-007232-7Review: 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 2-dimensional and 3-dimensional 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
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 well-fitting 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 long-term 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 McGraw-Hill,
New York (There may be later editions.)
may be consulted for further information.
Distances on the surface of the terrain, whether geodesic, on roads, cross-country, or straight-line, 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.