Search The Web

Today's Headlines

Monday, April 12, 2010

Microsoft Access Tips & Tricks: Great Circle Distance

The great circle distance or orthodromic distance is the shortest distance between any two points on the surface of a sphere measured along a path on the surface of the sphere (as opposed to going through the sphere's interior). Because spherical geometry is rather different from ordinary Euclidean geometry, the equations for distance take on a different form. The distance between two points in Euclidean space is the length of a straight line from one point to the other. On the sphere, however, there are no straight lines. In non-Euclidean geometry, straight lines are replaced with geodesics. Geodesics on the sphere are the great circles (circles on the sphere whose centers are coincident with the center of the sphere).

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 \phi_s,\lambda_s;\ \phi_f,\lambda_f\;\! be the coordinates (latitude, longitude) of the start point and final point respectively. Let \Delta\phi,\Delta\lambda\;\! be the differences in latitude and longitude respectively between the two points. The distance between them is calculated in terms of \Delta\widehat{\sigma}\;\!, which represents the angular distance between the points subtended at the center of the sphere. The actual distance between the points is r*\Delta\widehat{\sigma}\;\!, 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 \Delta\widehat{\sigma}\;\! as below:

arccosine formula for calculating great circle distanceUnfortunately, 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:

haversine formula for calculating great circle distance
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):

Vincenty formula for calculating great circle distance
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)
Arcsin(X) = ArcTan(X / √(-X * X + 1))
Below is the function that implements ArcCos:
Function ArcCos(X As Double) As Double
If Abs(X) > 1# Then X = X - Fix(X)
ArcCos = Atn(-X / Sqr(-X * X + 1)) + 2 * Atn(1)
End Function
Note 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.

Below is the function that implements ArcSin:
Function ArcSin(X As Double) As Double
If Abs(X) > 1# Then X = X - Fix(X)
ArcSin = Atn(X / Sqr(-X * X + 1))
End Function
Note 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.

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 Double
Const Pi As Double = 3.14159265358979

Degrees = Degrees + Minutes / 60#
Degrees = Degrees + Seconds / 3600#
DegToRad = Degrees * Pi / 180#
End Function
Notice 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.

First we will code up the ArcCosine formula as below:
Function ArcCosineDistance(SourceLatDeg As Double, _
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
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.

Next we code up the Haversine formula for great circle distance as below:
Function HaversineDistance(SourceLatDeg As Double, _
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
And, finally, the Vincenty formula as below:
Function VincentyDistance(SourceLatDeg As Double, _
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
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.

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:

Anonymous said...

Great! It helped me a lot

Blogannath said...

You are welcome!

arc length calculator said...

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.

Anonymous said...

Excellent Job.
Thanks for Share it.

Carlos R.

Visitors Country Map

Free counters!

Content From TheFreeDictionary.com

In the News

Article of the Day

This Day in History

Today's Birthday

Quote of the Day

Word of the Day

Match Up
Match each word in the left column with its synonym on the right. When finished, click Answer to see the results. Good luck!

 

Hangman

Spelling Bee
difficulty level:
score: -
please wait...
 
spell the word:

Search The Web