Last week we started a series about finding distances on a sphere (which approximates the shape of the earth), using a straightforward formula from spherical geometry. But in practice, that formula turns out not to be ideal, so a different formula is used when accuracy in all circumstances matters. That is this week’s topic: first to state it, and then to prove it.
The haversine formula
This question came to us in 2001:
Deriving the Haversine Formula I'm an SAS programmer at The Coca-Cola Company. It's been a long, long time since I have had to use trigonometry. I need to write a program module to calculate distances given longitude and latitude data. I am trying to find an object within a mile's radius of its location. Would you have the equation for this scenario? Thank you. I hope I am not too old for a response. :)
Some level of precision is needed, maybe even more than we can get assuming the earth is an exact sphere; we certainly want to do the best we can.
Doctor Rick answered, referring two of the pages we saw last week:
Hello, Dena, welcome to Ask Dr. Math. We're allowed to answer "old" people now and then, if the question is interesting enough! :-) You can find several expositions of the "Law of Cosines" for spherical trigonometry in our Dr. Math Archives. For example: Using Longitude and Latitude to Determine Distance http://mathforum.org/library/drmath/view/51711.html Distance using Latitude and Longitude http://mathforum.org/library/drmath/view/54680.html But a correspondent recently directed our attention to a Web site, The Geographic Information Systems FAQ (See Q5.1) http://www.faqs.org/faqs/geography/infosystems-faq/ This site points out that the cosine formula does not give accurate results for small distances, and presents a better formula, the Haversine Formula. The form of the Haversine Formula that uses the atan2() function is probably best for programming, assuming your programming language has an atan2() function. I will copy some excerpts here for your convenience, but I recommend looking at the site for the other information you will find there.
This is an old page from 1997, and not very user-friendly (search for Q5.1)! A somewhat better 2001 version can be found here, and the part with the formulas is here.
Note that in the description that follows, atan2(y, x) = arctan(y/x). However, there does seem to be some confusion about the order of the arguments. C/C++ puts the arguments in the order that I assume, as indicated by the Unix man page as well as the Visual C++ Programmer's Guide. On the other hand, Excel describes the function (not too clearly) as "ATAN2(x_num, y_num) returns the arctangent of the specified x- and y-coordinates, in radians." I tested it, and indeed it takes the denominator x first.
This is an issue we’ve dealt with a number of times in recommending formulas in Excel. This function returns the angle in standard position through the point (x, y); it is understandable that Excel would have you enter the coordinates of the point in that order, but it also makes sense to put the arguments in the order (y, x) to make a direct replacement for “atan(y/x)”, the arctangent of the slope. The difference between the two is that atan2 produces an angle in the appropriate quadrant, as atan can’t.
Just the FAQ
What follows was copied in condensed form from the source (which differed a little from the current version):
Presuming a spherical Earth with radius R (see below), and that 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(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2 c = 2 * atan2(sqrt(a), sqrt(1-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 current version does not include this atan2 version of the formula, but only the other version, using an arcsin with “bulletproofing” to prevent errors under special conditions (“nearly antipodal” points); the atan2 version doesn’t need that protection. We’ll look at that later.
Most computers require the arguments of trigonometric 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.)
Using the wrong units is a very common cause of errors people have written to us about!
A note on earth’s radius
The following comes under the heading “What value should I use for the radius of the Earth, R?” in the GIS FAQ:
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
Last week we used a different value, 3963 miles. Among the possible values, the two main ones from actual earth measurements are, according to Wikipedia:
- Equatorial radius a, or semi-major axis: distance from center to equator = 6,378.1370 km (3,963.1906 mi).
- Polar radius b, or semi-minor axis: distance from center to North and South Poles = 6,356.7523 km (3,949.9028 mi).
Another alternative is the “mean radius” as given by Wikipedia:
- Mean radius: \(\displaystyle\frac {2a+b}{3}\) = 6,371.0088 km (3,958.7613 mi).
The FAQ’s recommendation is equivalent to what is given in Wikipedia as
- Rectifying radius, giving a sphere with circumference equal to the perimeter of the ellipse described by any polar cross section of the ellipsoid: 6,367.4491 km (3,956.5494 mi).
- The mean of the two axes, \(\displaystyle\frac{a+b}{2}\) = about 6,367 km (3,956.547 mi)
What we used last week, 3963 miles, is the equatorial radius, so that those calculations can be expected to be just a little high.
But the FAQ is right not to give as much precision as Wikipedia for the averages, because they are, after all, averages! With about 4 significant digits, it will be reasonable to expect accuracy to within about a mile in 10,000. If we used the average radius to find a length along the equator, we would be about 0.17% low, or 1 mile out of 573, considerably less precise.
What the formula means
Doctor Rick continued:
I'll just add that the intermediate result "a" is the square of half of the straight-line distance (chord length) between the two points.
a = (AB/2)^2
= (sin^2(dlat/2)) + cos(lat1) * cos(lat2) * sin^2(dlon/2)
If you are interested in a derivation of the formula, I could provide it, but the figure would be hard to draw with "character graphics."
In 2001, it was hard to attach (or make) real graphics, so we did as you’ll be seeing below.
Deriving the formula
Dena wrote back:
Thank you very much for your quick response. I'm looking forward to visiting the Web sites you provided to review the expositions of the "Law of Cosines" and the Haversine Formula. I think I will always have a nature of curiosity about me. Therefore, I would LOVE to see a derivation of the formula.
Doctor Rick answered:
Hello again, Dena. You called my bluff! Good - now I have a reason to write up my derivation of the haversine formula. I had derived the cosine formula for great-circle distance, but not the haversine formula. I was reluctant to send the formula to you without having checked it out, so I derived it before I answered you. It turned out to be relatively straightforward.
We often do more work behind the scenes than we show! And this is not the only time we have wanted an occasion to write up the details of a problem.
What’s a haversine?
First, a word about the meaning of "haversine." This is the name of one of several lesser-known trigonometric functions that were once familiar to navigators and the like. The VERSINE (or versed sine) of angle A is 1-cos(A). The HAVERSINE is half the versine, or (1-cos(A))/2. With a little math, you can show that hav(A) = (1-cos(A))/2 = sin^2(A/2) where ^2 means squared. You see the form on the right twice in the haversine formula.
Here is that math: Using the double angle formula, \(\cos(2x) = 1 – 2\sin^2(x)\), we find that $$\text{hav}(A) = \frac{1-\cos(A)}{2}=\frac{1-\cos\left(2\cdot\frac{A}{2}\right)}{2}\\=\frac{1-\left(1 – 2\sin^2\left(\frac{A}{2}\right)\right)}{2}=\frac{2\sin^2\left(\frac{A}{2}\right)}{2}\\=\sin^2\left(\frac{A}{2}\right)$$
Wikipedia has a lot of information about this and a whole family of related functions, in the article on the versine. Jim Wilson has a nice article that also touches on our current topic, here.
Chords and angles
In proving the formula, we will work with a sphere of radius 1. The distance we get should be multiplied by R, the radius of the earth.
Let's start with a formula for the length of a chord subtending an angle theta in the unit circle (circle of radius 1).
O
/|\
/ | \
1 / | \ 1
/ | \
/ | \
/_____|_____\
A C B
O is the center of the circle, and A and B are on the circle, so AB is the chord. If angle AOB = theta, then angle AOC = theta/2 where OC is perpendicular to the chord AB. C bisects AB, since triangle AOB is isosceles. Segment AC = sin(theta/2), so the length of chord AB is 2*sin(theta/2).
$$AB=2\sin\left(\frac{\theta}{2}\right)$$
This formula is used often in working with regular polygons.
Chords of the sphere
Now for the 3-dimensional figure. Here is a sphere of radius 1, with a wedge removed to show the center, labeled O. The points we are interested in are A (lat1,lon1) and B (lat2,lon2). The lines of longitude lon1 and lon2 are shown as curved single lines, meeting at the north pole, N; the equator is shown as a double line (equal signs).
***********
**** N ****
**** / | \ ****
*** / | \ ***
** / | \ **
** / | \ **
* / | \ *
* / | \ *
* A----------------------D *
* /|\ | \ *
* | | \ | | *
* | | \ | | *
* | | \ | | *
* | | \ | | *
* C----------\-----------------B *
* | | - _\ | | *
* | | _O_ | *
* = | | _ _ | = *
* === | G_ _ | === *
* ====E______________________________F ===== *
* ========= ========= *
* ============== *
* *
* *
* *
** **
** **
*** ***
**** ****
**** ****
*************
A and B are opposite vertices of an isosceles trapezoid, ACBD, with additional vertices C (lat2,lon1) and D (lat1,lon2). I haven't tried to draw the straight lines (chords) connecting A to C and B to D, but I show the chords connecting A to D and B to C. I won't bother to prove that the trapezoid ACBD is planar.
Here is a somewhat clearer picture:
Our goal is the length of arc AB; we’ll be first finding the length of segment AB, then using that to find the angle δ = AOB and finally the arc.
If you join A to O and C to O, the angle AOC is the difference between lat1 and lat2, or dlat. The length of chord AC is therefore 2*sin(dlat/2), using the chord-length formula derived above. Chord BD is the same length. Points E and F are the points where the latitude lines lat1 and lat2 meet the equator. The length of EF is 2*sin(dlon/2) by the chord-length formula. Points A and D are on a circle of constant latitude lat1. The radius of this circle is cos(lat1). You can see this by dropping a perpendicular from A to OE, meeting OE at G. The angle EOA is lat1, so OG = cos(lat1). But this is equal to the radius of the latitude circle at A. Therefore the length of chord AD is 2*sin(dlon/2)*cos(lat1). Similarly, the length of chord CB is 2*sin(dlon/2)*cos(lat2).
So far, we have $$AC = BD = 2\sin\left(\frac{dlat}{2}\right)\\ AD = 2\sin\left(\frac{dlon}{2}\right)\cos(lat_1)\\ CB = 2\sin\left(\frac{dlon}{2}\right)\cos(lat_2)$$
Diagonal of a trapezoid
Now we have the lengths of all the sides of (planar) trapezoid ADBC. The straight-line distance AB (passing through the earth) is a diagonal of that isosceles trapezoid.
Now let's go back to 2 dimensions, and find the length of a diagonal of an isosceles trapezoid.
A_________________D
/| * \
/ | * \
/ | * \
/ | * \
/____|_________________*__\
C H B
The length of CH, where AH is perpendicular to CB, is (CB-AD)/2. By the Pythagorean theorem,
AH^2 = AC^2 - CH^2
= AC^2 - (CB-AD)^2/4
The length of HB is (CB+AD)/2. Using Pythagoras again, we find that
AB^2 = AH^2 + HB^2
= AC^2 - (CB-AD)^2/4 + (CB+AD)^2/4
= AC^2 + CB*AD
$$AB^2 = AC^2 + CB\cdot AD$$
The square of the diagonal is the sum of the square of either leg and the product of the bases (which is the square of their geometric mean). Sounds reasonable.
Putting it together
We're ready to plug in the lengths of chords AC, AD, and CB from our work with the sphere:
AB^2 = 4*(sin^2(dlat/2) + cos(lat1)*cos(lat2)*sin^2(dlon/2))
The intermediate result a is the square of half the chord length AB:
a = (AB/2)^2
= (sin^2(dlat/2)) + cos(lat1) * cos(lat2) * sin^2(dlon/2)
So, $$AB^2 = \left[2\sin\left(\frac{dlat}{2}\right)\right]^2 + \left[2\sin\left(\frac{dlon}{2}\right)\cos(lat_1)\right]\left[2\sin\left(\frac{dlon}{2}\right)\cos(lat_2)\right]\\ = 4\sin^2\left(\frac{dlat}{2}\right) + 4\sin^2\left(\frac{dlon}{2}\right)\cos(lat_1)\cos(lat_2)$$
Dividing this by 4, we get the square of the half-chord, $$a = \left(\frac{AB}{2}\right)^2 = \sin^2\left(\frac{dlat}{2}\right) + \cos(lat_1)\cos(lat_2)\sin^2\left(\frac{dlon}{2}\right)$$
This is part of the formula we are trying to prove; it appears that this intermediate value a was chosen partly to eliminate the 4, and partly, as we’ll see in a moment, because half the chord is what we need anyway!
The final step is to find the central angle AOB that corresponds to this chord length. The arctan formula can be found from this figure:
O
/|\
/ | \
/ | \
/ | \
1 / | \ 1
/ | \
/ | \
/ | \
/ | \
/ | \
/ sqrt(1-a)| \
/ | \
/____________|____________\
A sqrt(a) C B
If AC = sqrt(a), then Pythagoras says that
OC = sqrt(OA^2 - AC^2)
= sqrt(1-a)
Therefore tan(<AOC) = AC/OC = sqrt(a)/sqrt(1-a), or
c = 2 * arctan(sqrt(a)/sqrt(1-a))
where c is the angle AOB.
All that’s left for the formula is to multiply this angle (in radians) by the radius R, to get the arc length (great-circle distance). The formula, again, is:
$$dlon = lon_2 – lon_1\\ dlat = lat_2 – lat_1\\ a = \sin^2\left(\frac{dlat}{2}\right) + \cos(lat_1)\cos(lat_2)\sin^2\left(\frac{dlon}{2}\right)\\ \delta = 2 \text{atan2}\left(\sqrt{a}, \sqrt{1-a}\right)\\ d = R\delta$$
The arcsin version
In the FAQ from which we got the formula, the formula is given using a sine instead of the atan2 that we are using. That saves some contortions, as angle AOC is just the arcsin of \(\sqrt{a}\). But then they have to add a check. Here is what they say:
Haversine Formula: 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 opposite 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).
The atan2 function will never have these errors. Instead, in the antipodal case, we will be dividing by a very small number, and might have overflow in a computer.
Comparing the haversine and cosine formulas
The GIS FAQ from which we originally got this formula explicitly warns against using the cosine formula we studied last week. Here is what it says:
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.
And one minute of arc is $$\frac{1}{60}\cdot \frac{1}{360}\cdot 2\pi R = \frac{2\pi\cdot 3956}{60\cdot360} = 1.15\text{ miles}$$
In contrast, when we use the haversine formula for small distances, we just take the sine of small angles, and then the inverse sine (or inverse tangent) of a small number, which behaves nicely, being approximately proportional.
So let’s do one example, at a fairly small distance, using both formulas.
Since last week we found the distance from San Francisco to Paris, let’s suppose we are in Paris and want to find the distance from the Arc de Triomphe, with coordinates 48.8738° N, 2.2950° E, to the Place de la Concorde, with coordinates 48.8656° N, 2.3212° E. I chose this nice straight walk because Google can tell us the distance:
We expect the correct distance to be a little less than 2.2 kilometers.
First, I’ll use the cosine formula, which I’ll do in the format used by the GIS FAQ we saw, which breaks it into steps we can examine:
- a = sin(48.8738°) * sin(48.8656°) = 0.5673338013
- b = cos(48.8738°) * cos(48.8656°) * cos(2.3212° – 2.2950°) = 0.4326661431
- c = arccos(a + b) = arccos(0.9999999445) = 0.00033309707149 rad
- d = 6367 * 0.00033309707149 = 2.1208290542 ≈ 2.12 km
What we see here is that we had to take the inverse cosine of a number very close to 1. When I repeated the calculation on a hand calculator with less precision than the Windows calculator I used initially (which showed many more digits than I have copied), I got the same results (to five decimal places), so the loss of precision was not serious; but in principle, the formula is inadequate – especially if you were using trig tables rather than a modern calculator!
Now I’ll use the haversine formula, again from the FAQ:
- dlon = 2.3212° – 2.2950° = 0.0262°
- dlat = 48.8656° – 48.8738° = -0.0082°
- a = sin^2(-0.0082/2) + cos(48.8738°) * cos(48.8656°) * sin^2(0.0262/2) = 0.00000002773841450
- c = 2 * arcsin(min(1,sqrt(a))) = 0.0003330970714 rad
- d = 6367 * 0.0003330970714 = 2.1208290542 ≈ 2.12 km
The answer is essentially the same; and this, too, worked just as well on a hand calculator (though for a while it looked as if it had lost many digits in calculating a, as it did not display that tiny number in scientific notation, but as 0.000000028).
Pingback: Distances on Earth 3: Planar Approximation – The Math Doctors