Earth Surface Distance in T-SQL

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;

  N'Kilometers' AS Units,
  @expectedDistanceKilometers AS Expected,
  @calculatedDistanceKilometers AS Calculated,
  ABS(@expectedDistanceKilometers - @calculatedDistanceKilometers)
    / @expectedDistanceKilometers AS Skew
  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:

Units Expected Calculated Skew
Kilometers 7527.806 7521.74343697558 0.00080535590641083
Miles 4677.562 4673.79632989539 0.000805049746985879

It seems that for the two test points at least, the skew based on Google Earth as the standard is about 8/100ths of one percent. Over the years, I’ve successfully adapted this haversine code to C#, JavaScript and Java, too. The haversine formula performs better for short distances than using a simple law of cosines-based formula for all sorts of reasons that scientists understand. However, it’s got a margin of error of up to 0.3% in some cases which can be too great for applications that require high precision.

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.

Leave a comment

Your email address will not be published.