If you are interested, you can find my earlier posts on finding the median, the mode, the geometric and harmonic means, ranking every row in a query, selecting random rows out of a table, calculating running sums and averages, calculating running differences, creating histograms, calculating probability masses out of given data, calculating cumulative distributions out of given data, finding percentile scores, percentile values, calculating distinct counts, full outer joins, parameter queries, crosstab queries, working with system objects, listing table fields, finding unmatched rows, calculating statistics with grouping, job-candidate matching, and job-candidate matching with skill levels.

Between any two points on a sphere which are not directly opposite each other, there is a unique great circle. The two points separate the great circle into two arcs. The length of the shorter arc is the great circle distance between the points.

Between two points which are directly opposite each other, called antipodal points, there are infinitely many great circles, but all great circle arcs between antipodal points have the same length, i.e. half the circumference of the circle, or πr, where r is the radius of the sphere.

Great circle distances are very useful for finding the shortest distances between points on a sphere. Since the shape of the earth can be approximated to be a sphere with negligible loss of accuracy, the distance between different points on the earth can be calculated to a great degree of accuracy using great circle distance formulae.

If we have a database that has a table that lists the coordinates of various points on the earth (coordinates of cities, airports, etc.), then we may want to use these formulae to calculate the great circle distance between them at some point. For the purpose of this post, we will assume that our database contains coordinates in the form of degrees, minutes and seconds for both latitude and longitude.

We will also assume that latitude values can go from -90 to +90 degrees, with latitudes of points on the northern hemisphere of the earth being positive and latitudes of points on the southern hemisphere being negative. All points on the equator will have a latitude of 0, while the north pole will have a latitude of +90 and the south pole would have a latitude of -90.

Similarly, we will assume that longitudes range from -180 degrees to +180 degrees. Points on the prime meridian that passes through Greenwich, England, would have a longitude of 0, while points on the opposite side of the earth from that line would have a longitude of either +180 or -180 (the +180 and -180 lines of longitude coincide on top of each other). We will assume that points in the western hemisphere have negative longitudes while points in the eastern hemisphere have positive longitudes.

Now, that we have established the basics of the table structure and conventions regarding signs and quantities, we can look at the formulae used for computing the great circle distance between two points. Let be the coordinates (latitude, longitude) of the start point and final point respectively. Let be the differences in latitude and longitude respectively between the two points. The distance between them is calculated in terms of , which represents the angular distance between the points subtended at the center of the sphere. The actual distance between the points is r*, where r is the radius of the sphere.

In the case of the earth, r is 6371.01 km, 3958.76 statute miles or 3440.07 nautical miles.

The most common formula one encounters is called the ArcCosine formula. It calculates as below:

Unfortunately, this formula has large rounding errors when the distances are very small, so it is not commonly used. One of the alternatives to the ArcCosine formula is the Haversine formula below:

Although this formula is accurate for most distances, it too suffers from rounding errors for the special (and somewhat unusual) case of antipodal points (on opposite ends of the sphere). A more complicated formula that is accurate for all distances is the following special case (a sphere, which is an ellipsoid with equal major and minor axes) of the Vincenty formula (which more generally is a method to compute distances on ellipsoids):

In this post, I will provide code for VBA functions that implement each of these formulae to calculate the great circle distance between any two points, whose coordinates are passed into these functions as arguments. You can compare the performance of these functions against each other and use whichever one suits your purpose best. In most cases, they will come up with the exact same answer, so it is only a matter of preference which one is chosen. But given the caveat about errors in very short distances when using the ArcCosine formula, and the complexity of the Vincenty formula, I prefer the Haversine formula.

Notice that Microsoft Access lacks functions for ArcCos() and ArcSin(), and therefore these quantities have to be derived using other trigonometric functions available in Access. Therefore I will provide these functions first before I go on to the functions for the distance. Those functions for distance will use the ArcCos and ArcSin functions below. The ArcCos and ArcSin functions are written based on the trigonometric identities below:

ArcCos(X) = ArcTan(-X / √(-X * X + 1)) + 2 * ArcTan(1)Below is the function that implements ArcCos:

Arcsin(X) = ArcTan(X / √(-X * X + 1))

Function ArcCos(X As Double) As DoubleNote that the argument to the function must fall in the range between -1 and 1 inclusive (the cosine of any angle can not fall outside this range). When the argument is outside this range I remove the integer portion of the argument using the fix() function and use only the fractional part in my calculations.

If Abs(X) > 1# Then X = X - Fix(X)

ArcCos = Atn(-X / Sqr(-X * X + 1)) + 2 * Atn(1)

End Function

Below is the function that implements ArcSin:

Function ArcSin(X As Double) As DoubleNote again, that the argument to this function must also fall in the range between -1 and 1 inclusive (the sin of any angle can not fall outside this range). When the argument is outside this range, I remove the integer portion of the argument using the fix() function and use only the fractional part in my calculations.

If Abs(X) > 1# Then X = X - Fix(X)

ArcSin = Atn(X / Sqr(-X * X + 1))

End Function

Before we code up the great circle distance calculations, we also need to convert the coordinates of the points into radians. Microsoft Access's trigonometric functions expect angles to be in radians rather than degrees. Thus, we need the function below to convert any given coordinates into radians.

Function DegToRad(Degrees As Double, Minutes As Double, Seconds As Double) As DoubleNotice that we have defined a constant named pi to equal 3.1415926535897932, and we take advantage of the identity that angle in radians = angle in degrees*π/180.0.

Const Pi As Double = 3.14159265358979

Degrees = Degrees + Minutes / 60#

Degrees = Degrees + Seconds / 3600#

DegToRad = Degrees * Pi / 180#

End Function

First we will code up the ArcCosine formula as below:

Function ArcCosineDistance(SourceLatDeg As Double, _Notice that we have constants for the radius of the earth in all three units. One can change the assignment statement to change the variable EarthRadius, and therefore, the great circle distance calculated to any of the three units depending on one's preferences. Also notice how the same variable (in this case, the name of the function, which is the return variable's name) is used to accumulate the calculation step by step all the way to the final step.

SourceLatMin As Double, SourceLatSec As Double, _

SourceLongDeg As Double, SourceLongMin As Double, _

SourceLongSec As Double, DestLatDeg As Double, _

DestLatMin As Double, DestLatSec As Double, _

DestLongDeg As Double, DestLongMin As Double, _

DestLongSec As Double) As Double

Const EarthRadiusKM As Double = 6371.01 'Kilometers

Const EarthRadiusSM As Double = 3958.76 'Statute Miles

Const EarthRadiusNM As Double = 3440.07 'Nautical Miles

Dim EarthRadius As Double

Dim SourceLatRad As Double

Dim SourceLongRad As Double

Dim DestLatRad As Double

Dim DestLongRad As Double

Dim DeltaLongRad As Double

EarthRadius = EarthRadiusKM

SourceLatRad = DegToRad(SourceLatDeg, SourceLatMin, SourceLatSec)

SourceLongRad = DegToRad(SourceLongDeg, SourceLongMin, SourceLongSec)

DestLatRad = DegToRad(DestLatDeg, DestLatMin, DestLatSec)

DestLongRad = DegToRad(DestLongDeg, DestLongMin, DestLongSec)

DeltaLongRad = Abs(SourceLongRad - DestLongRad)

ArcCosineDistance = Sin(SourceLatRad) * Sin(DestLatRad)

ArcCosineDistance = ArcCosineDistance + Cos(SourceLatRad) * Cos(DestLatRad) * Cos(DeltaLongRad)

ArcCosineDistance = ArcCos(ArcCosineDistance)

ArcCosineDistance = ArcCosineDistance * EarthRadius

End Function

Next we code up the Haversine formula for great circle distance as below:

Function HaversineDistance(SourceLatDeg As Double, _And, finally, the Vincenty formula as below:

SourceLatMin As Double, SourceLatSec As Double, _

SourceLongDeg As Double, SourceLongMin As Double, _

SourceLongSec As Double, DestLatDeg As Double, _

DestLatMin As Double, DestLatSec As Double, _

DestLongDeg As Double, DestLongMin As Double, _

DestLongSec As Double) As Double

Const EarthRadiusKM As Double = 6371.01 'Kilometers

Const EarthRadiusSM As Double = 3958.76 'Statute Miles

Const EarthRadiusNM As Double = 3440.07 'Nautical Miles

Dim EarthRadius As Double

Dim SourceLatRad As Double

Dim SourceLongRad As Double

Dim DestLatRad As Double

Dim DestLongRad As Double

Dim DeltaLatRad As Double

Dim DeltaLongRad As Double

EarthRadius = EarthRadiusKM

SourceLatRad = DegToRad(SourceLatDeg, SourceLatMin, SourceLatSec)

SourceLongRad = DegToRad(SourceLongDeg, SourceLongMin, SourceLongSec)

DestLatRad = DegToRad(DestLatDeg, DestLatMin, DestLatSec)

DestLongRad = DegToRad(DestLongDeg, DestLongMin, DestLongSec)

DeltaLatRad = Abs(SourceLatRad - DestLatRad)

DeltaLongRad = Abs(SourceLongRad - DestLongRad)

HaversineDistance = Sin(DeltaLatRad / 2#) * Sin(DeltaLatRad / 2#)

HaversineDistance = HaversineDistance + _

Cos(SourceLatRad) * Cos(DestLatRad) * Sin(DeltaLongRad / 2#) * Sin(DeltaLongRad / 2#)

HaversineDistance = 2 * ArcSin(Sqr(HaversineDistance))

HaversineDistance = HaversineDistance * EarthRadius

End Function

Function VincentyDistance(SourceLatDeg As Double, _Notice that we have calculated the numerator and denominator of the Vincenty formula separately rather than doing the calculations all at once. This is not only a concession to the complexity of the formula itself, but also to the fact that we need to deal with the special case of the denominator becoming zero (this is the only one of the three formulae that has a denominator and involves division). This is accomplished by checking the absolute value of the denominator against a very small constant, Epsilon. When the denominator is less than epsilon, we assume that the denominator is zero, and the expression inside the parentheses becomes infinity. As we know, the ArcTangent of infinity is π/2 radians, so we set the subtended angle to π/2 when the denominator is close to zero rather than attempting a division by a number that close to zero, and risking a run-time error.

SourceLatMin As Double, SourceLatSec As Double, _

SourceLongDeg As Double, SourceLongMin As Double, _

SourceLongSec As Double, DestLatDeg As Double, _

DestLatMin As Double, DestLatSec As Double, _

DestLongDeg As Double, DestLongMin As Double, _

DestLongSec As Double) As Double

Const EarthRadiusKM As Double = 6371.01 'Kilometers

Const EarthRadiusSM As Double = 3958.76 'Statute Miles

Const EarthRadiusNM As Double = 3440.07 'Nautical Miles

Const Pi As Double = 3.14159265358979

Const Epsilon As Double = 0.000000000001

Dim EarthRadius As Double

Dim SourceLatRad As Double

Dim SourceLongRad As Double

Dim DestLatRad As Double

Dim DestLongRad As Double

Dim DeltaLongRad As Double

Dim Numerator1 As Double

Dim Numerator2 As Double

Dim Numerator As Double

Dim Denominator As Double

EarthRadius = EarthRadiusKM

SourceLatRad = DegToRad(SourceLatDeg, SourceLatMin, SourceLatSec)

SourceLongRad = DegToRad(SourceLongDeg, SourceLongMin, SourceLongSec)

DestLatRad = DegToRad(DestLatDeg, DestLatMin, DestLatSec)

DestLongRad = DegToRad(DestLongDeg, DestLongMin, DestLongSec)

DeltaLongRad = Abs(SourceLongRad - DestLongRad)

Numerator1 = Cos(DestLatRad) * Sin(DeltaLongRad)

Numerator1 = Numerator1 * Numerator1

Numerator2 = Cos(SourceLatRad) * Sin(DestLatRad)

Numerator2 = Numerator2 - Sin(SourceLatRad) * Cos(DestLatRad) * Cos(DeltaLongRad)

Numerator2 = Numerator2 * Numerator2

Numerator = Numerator1 + Numerator2

Numerator = Sqr(Numerator)

Denominator = Sin(SourceLatRad) * Sin(DestLatRad)

Denominator = Denominator + Cos(SourceLatRad) * Cos(DestLatRad) * Cos(DeltaLongRad)

If Abs(Denominator) < Epsilon Then

VincentyDistance = Pi / 2#

Else

VincentyDistance = Numerator / Denominator

VincentyDistance = Atn(VincentyDistance)

End If

VincentyDistance = VincentyDistance * EarthRadius

End Function

The code above (all functions) has been checked in Access 2003. It should work in any version of Access from Access 97 on up, but please do let me know through the comments if you have problems using the code.

Hope this post has been helpful in solving any problems you might have had with great circle distance calculations in Access. If you have any problems or concerns with the VBA code in this lesson, please feel free to let me know by posting a comment. If you have other questions on Access that you would like me to address in future lessons, please feel free to let me know through your comments too. Good luck!

## 4 comments:

Great! It helped me a lot

You are welcome!

This can be a really good examine for me, Must admit that you simply are one particular of the best bloggers I ever saw. Thanks for posting this informative article.

Excellent Job.

Thanks for Share it.

Carlos R.

Post a Comment