A few years ago, I was working on a project where I needed to calculate the distance between pairs of points on the Earth’s surface. With a bit of research, I found an implementation of the haversine formula written in Python on John D. Cook’s Standalone Numerical Code site. Given a pair of latitude and longitude values, John’s code produces a unit coefficient that can be multiplied by the radius of the sphere to yield the distance in whatever unit of measurement you might want.
Having the Python code was great but I needed a version of the algorithm in F#. After searching for some time and not finding an F# implementation, I decided to write one based on John’s Python version and put it back into the public domain. John published my F# code into his Standalone Numerical Code library here for everyone to use freely.
The exercise for today is to write the haversine formula in Transact-SQL. I’ll start by proposing a test. Between two points that are well-known and using Google Earth as the standard, I define the following variables.
-- Richmond, VA USA DECLARE @latitudeA FLOAT = 37.542979; DECLARE @longitudeA FLOAT = -77.469092; -- Sao Paulo, Brazil DECLARE @latitudeB FLOAT = -23.548943; DECLARE @longitudeB FLOAT = -46.638818; DECLARE @expectedDistanceMiles FLOAT = 4677.562; DECLARE @expectedDistanceKilometers FLOAT = 7527.806;
Google Earth estimates that these two pairs of coordinates between Richmond, Virginia USA and São Paulo, Brazil are roughly 4,678 miles or 7,528 kilometers apart. Having no better standard handy, those will have to suffice for some testing later on. Next comes the heart of the haversine formula.
DECLARE @phi1 FLOAT = RADIANS(90 - @latitudeA); DECLARE @phi2 FLOAT = RADIANS(90 - @latitudeB); DECLARE @theta1 FLOAT = RADIANS(@longitudeA); DECLARE @theta2 FLOAT = RADIANS(@longitudeB); DECLARE @arcUnitDistance FLOAT = ACOS(SIN(@phi1) * SIN(@phi2) * COS(@theta1 - @theta2) + COS(@phi1) * COS(@phi2));
I recommend saving this as a user-defined function in your SQL repository called something like ArcUnitDistance. The reason for that naming is because the value that this calculation produces is the distance between the two points supplied on a surface that has a radius of one in whatever unit of measure you’re going to apply. Now it’s time to put the code to the test.
DECLARE @earthRadiusMiles FLOAT = 3960; DECLARE @earthRadiusKilometers FLOAT = 6373; DECLARE @calculatedDistanceMiles FLOAT = @arcUnitDistance * @earthRadiusMiles; DECLARE @calculatedDistanceKilometers FLOAT = @arcUnitDistance * @earthRadiusKilometers; SELECT N'Kilometers' AS Units, @expectedDistanceKilometers AS Expected, @calculatedDistanceKilometers AS Calculated, ABS(@expectedDistanceKilometers - @calculatedDistanceKilometers) / @expectedDistanceKilometers AS Skew UNION ALL SELECT N'Miles' AS Units, @expectedDistanceMiles AS Expected, @calculatedDistanceMiles AS Calculated, ABS(@expectedDistanceMiles - @calculatedDistanceMiles) / @expectedDistanceMiles AS Skew
Commonly used values for the radius of the Earth when using the haversine formula are 3,960 miles or 6,373 kilometers. Of course, the Earth isn’t a sphere. It’s a spheroid so you may find other radius values that you trust more. To find the distance between my test points in those units of measurement, I simply need to multiply the @arcUnitDistance by the radius values. Then, it’s easy to calculate the skew from the expectation and print it out. The test yields these results:
If you want much greater accuracy when calculating distances on spheroids like planet Earth, you should check out the Vincenty Formula instead. It’s marginally more complex than haversine but yields typically better results. I hope you find this version of the haversine formula written in Transact-SQL useful for a variety of distance measurement applications.