Search The Web

Today's Headlines

Monday, April 19, 2010

Microsoft Access Tips & Tricks: Great Circle Initial Heading

In the previous post, we saw how to calculate the great circle distance between two points on the earth based on their latitudes and longitudes. Another important property of the great circle route between two points is the initial heading of such a route. In this post, we will explore a couple of different formulae for calculating the great circle initial heading, and how to implement them in Access VBA.

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, job-candidate matching with skill levels, and great circle distances.

As in the previous post, let us establish some conventions first. The great circle initial heading can be some angle between 0 and 360 degrees (0 and 2π radians). A heading of 0 degrees signifies a heading due north. We then move clockwise, getting a heading of 90 degrees for due east, 180 degrees for due south, 270 degrees for due west and 360 degrees for due north. Thus due north can be represented by both 0 degrees and 360 degrees. The corresponding angles in radians are 0 and 2π radians for north, π/2 radians for east, π radians for south, and 3π/2 radians for west.

Also note that the initial heading from the north pole is always towards the south, and hence, 180 degrees (or π radians), and the initial heading from the south pole is always towards the north, and hence, 0 or 360 degrees (0 or 2π radians), regardless of the destination. That is because the great circle path from either pole always follows a line of longitude to the destination (if the destination is the other pole, then the great circle path can be any line of longitude, otherwise, the line of longitude for the great circle path is the longitude of the destination).

As in the previous post, assume that φs, λs; φf, λf are the coordinates (latitude, longitude) of the start point and final point respectively. Let Δλ be the difference in longitude between the start point and final point. Similarly, let Δσ be the angle subtended by the two points at the center of the sphere.

Let θ represent the great circle initial heading of the great circle path between the two points. There are two formulae that are commonly used to find θ. The first of these formulae can be used when Δσ is known, or has already been calculated (note that this is often the case because the great circle initial heading is often calculated in conjunction with the great circle distance, and the calculation of the great circle distance requires the calculation of the central angle between the two points). This formula is given below:

If sin(Δλ) < 0
θ = arccos((sin(φf) - sin(φs)*cos(Δσ))/(sin(Δσ)*cos(φs)))

Otherwise,
θ = 2π - arccos((sin(φf) - sin(φs)*cos(Δσ))/(sin(Δσ)*cos(φs)))

We will refer to the above formula as the ArcCosine formula for great circle initial heading.

The second formula assumes that the central angle, Δσ, is not known. So, this formula can be used when we want to calculate the great circle initial heading without going to the trouble of calculating the great circle distance at the same time or beforehand. This formula is presented below:

θ = mod(arctan((sin(λf - λs)*cos(φf))/(cos(φs)*sin(φf) - sin(φs)*cos(φf)*cos(λf - λs))), 2π)

We will refer to the above formula as the ArcTangent formula for great circle initial heading. Notice that we use the mod() function in the above formula. Mod(y,x) is the remainder obtained by dividing y by x, and always lies in the range 0 <= mod(y,x) < x. For instance: mod(2.3,2) = 0.3, and mod(-2.3,2) = 1.7.

We will assume that our database already contains the code from the previous post for ArcCos(), HaversineDistance(), etc. We will call those functions in our code below, but I will not repeat the code for them in this post. Note that in the VBA that follows, I have used the Haversine formula to calculate the central angle between the points on the sphere for the ArcCosine formula, but you can use any of the three distance formulae in the previous post depending on your personal preference.

Also, in the VBA presented in this post, the great circle initial heading will be calculated in radians. To convert it into degrees if needed, you have to multiply the calculated number by 180/π.

The function implementing the ArcCosine formula is presented below:

Note that the great circle distance between the two points is calculated by using the Haversine formula for great circle distances, and the answer is then divided by the radius of the earth to get the central angle for use with the ArcCosine formula. The cases of the source being the north pole or south pole are dealt with separately since they require no calculations.
Function ArcCosineHeading(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 Pi As Double = 3.14159265358979
Const EarthRadiusKM As Double = 6371.01 'Kilometers
Const Epsilon As Double = 0.000000000000001

Dim SourceLatRad As Double
Dim SourceLongRad As Double
Dim DestLatRad As Double
Dim DestLongRad As Double
Dim DeltaLongRad As Double
Dim CentralAngle As Double

SourceLatRad = DegToRad(SourceLatDeg, SourceLatMin, SourceLatSec)
SourceLongRad = DegToRad(SourceLongDeg, SourceLongMin, SourceLongSec)
DestLatRad = DegToRad(DestLatDeg, DestLatMin, DestLatSec)
DestLongRad = DegToRad(DestLongDeg, DestLongMin, DestLongSec)

DeltaLongRad = Abs(SourceLongRad - DestLongRad)

If SourceLatRad = Pi / 2# Then 'If the source is the north pole
ArcCosineHeading = Pi
ElseIf SourceLatRad = -1 * Pi / 2# Then 'If the source is the south pole
ArcCosineHeading = 0
Else
DeltaLongRad = Abs(SourceLongRad - DestLongRad)
CentralAngle = HaversineDistance(SourceLatDeg, SourceLatMin, SourceLatSec,_
SourceLongDeg, SourceLongMin, SourceLongSec, _
DestLatDeg, DestLatMin, DestLatSec, _
DestLongDeg, DestLongMin, DestLongSec)
CentralAngle = CentralAngle / EarthRadiusKM
If CentralAngle <= Epsilon Then 'the two points are the same
ArcCosineHeading = 0
Else
ArcCosineHeading = Sin(DestLatRad) - Sin(SourceLatRad) * Cos(CentralAngle)
ArcCosineHeading = ArcCosineHeading / (Sin(CentralAngle) * Cos(SourceLatRad))
ArcCosineHeading = ArcCos(ArcCosineHeading)

If Sin(DeltaLongRad) >= 0 Then
ArcCosineHeading = 2 * Pi - ArcCosineHeading
End If
End If
End If
End Function
Note that the calculation of the central angle relies on the great circle distance being calculated in kilometers. If a different unit is used, the radius of the earth in that unit should be used in the division to derive the central angle, otherwise large errors could result!

Since the calculation involves a division, we have to make sure that we don't divide by zero accidentally during run-time. The denominator consists of cos(φs) which is already taken care of, because this quantity becomes zero only when the source is one of the poles. The other quantity in the denominator is sin(Δσ), which can become zero only when the central angle between the two points is zero. This can happen only when the two points are one and the same, so we detect this occurence (by comparing the central angle against a very small quantity) and assign an initial heading of 0 if the central angle is smaller than that infinitesimal quantity.

Now, the implementation of the ArcTangent formula is given below:
Function ArcTangentHeading(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 Pi As Double = 3.14159265358979
Const Epsilon As Double = 0.000000000000001

Dim SourceLatRad As Double
Dim SourceLongRad As Double
Dim DestLatRad As Double
Dim DestLongRad As Double
Dim DeltaLongRad As Double
Dim CentralAngle As Double
Dim Denominator As Double
Dim Numerator As Double

SourceLatRad = DegToRad(SourceLatDeg, SourceLatMin, SourceLatSec)
SourceLongRad = DegToRad(SourceLongDeg, SourceLongMin, SourceLongSec)
DestLatRad = DegToRad(DestLatDeg, DestLatMin, DestLatSec)
DestLongRad = DegToRad(DestLongDeg, DestLongMin, DestLongSec)

If SourceLatRad = Pi / 2# Then 'If the source is the north pole
ArcTangentHeading = Pi
ElseIf SourceLatRad = -1 * Pi / 2# Then 'If the source is the south pole
ArcTangentHeading = 0
Else
Denominator = Cos(DestLatRad) * Sin(DestLatRad)
Denominator = Denominator - _
Sin(SourceLatRad) * Cos(DestLatRad) * Cos(DestLongRad - SourceLongRad)
If Denominator <= Epsilon Then
ArcTangentHeading = Pi / 2#
Else
Numerator = Sin(DestLongRad - SourceLongRad) * Cos(DestLatRad)
ArcTangentHeading = Numerator / Denominator
ArcTangentHeading = Atn(ArcTangentHeading)
If ArcTangentHeading < 0 Then
ArcTangentHeading = ArcTangentHeading + 2# * Pi
End If
End If
End If
End Function
Note that Microsoft Access's implementation of the Mod operator rounds floating point numbers to integers. So, instead of using the mod operator, we take advantage of the fact that the Atn() function always returns a value between -π/2 and +π/2. Thus, if the return value from the Atn() function is less than zero, we simply add 2π to the result to simulate the modulus operation with 2π!

Also note that we have calculated the denominator in the formula separately. If it turns out to be zero, we know that the arctangent of the final resulting quantity will be π/2 (since the arctangent of ∞ is π/2). So, we set the initial heading to that value instead of going through with the calculation and risking a division by zero during run-time.

The code above (both 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 initial heading 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!

1 comment:

Jeff said...

One bug in ArcCosineHeading: after determining source is not one of the poles, calculation of DeltaLongRad should not use absolute value.

Thanks for the code!

ps - anonymous comments not working

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