# GPS Distance Calculations

## Recommended Posts

I decided to make this program to calculate distances based off of GPS coordinates. My next step is to calculate it with differences in elevation between the GPS points. I got the code from here

http://www.movable-type.co.uk/scripts/latlong-vincenty.html

http://www.movable-type.co.uk/scripts/latlong.html

The Autoit _Radian() function doesn't output all the time because sometimes the IsNumber() function thinks the float is not a number... don't ask me why... so the /57.2957795130823 in my code used to be _Radian(). I got the Haversine function to work correctly but the Vincenty one is still broken. If you see anything in my code that doesn't add up, that would be helpful...

The atan2() function missing might be the problem but i kinda doubt it.

http://www.autoitscript.com/forum/index.php?showtopic=17501

The atan2() function widely used here takes two arguments, atan2(y, x), and computes the arc tangent of the ratio y/x. It is more flexible than atan(y/x), since it handles x=0, and it also returns values in all 4 quadrants -π to +π (the atan function returns values in the range -π/2 to +π/2).

```#include <Math.au3>
#include <GUIConstants.au3>
Global \$PI = 4 * ATan(1)
;Global \$PI = 3.14159265358979323846264338327950288419716939937510

; Calculate geodesic distance (in m) between two points specified by latitude/longitude (in DMS)
; using Vincenty inverse formula for ellipsoids
func distVincenty(\$lat1, \$lon1, \$lat2, \$lon2, \$radius = 6378137)
\$lat1 = DMStoDD(\$lat1)
\$lat2 = DMStoDD(\$lat2)
\$lon1 = DMStoDD(\$lon1)
\$lon2 = DMStoDD(\$lon2)

; WGS-84 ellipsiod
\$a = 6378137
\$b_new = 6356752.3142
\$f = 1/298.257223563

\$L = (\$lon2-\$lon1)/57.2957795130823
\$U1 = ATan ((1-\$f) * Tan (\$lat1/57.2957795130823))
\$U2 = ATan ((1-\$f) * Tan (\$lat2/57.2957795130823))
\$sinU1 = Sin (\$U1)
\$cosU1 = Cos (\$U1)
\$sinU2 = Sin (\$U2)
\$cosU2 = Cos (\$U2)

\$lambda = \$L
\$lambdaP = 2 * \$PI

\$iterLimit = 256
While (Abs(\$lambda-\$lambdaP) > 10^-12) And (\$iterLimit > 0)
\$iterLimit -= 1
\$sinLambda = Sin(\$lambda)
\$cosLambda = Cos(\$lambda)
\$sinSigma = Sqrt ((\$cosU2*\$sinLambda) * (\$cosU2*\$sinLambda) + (\$cosU1*\$sinU2-\$sinU1*\$cosU2*\$cosLambda) * (\$cosU1*\$sinU2-\$sinU1*\$cosU2*\$cosLambda))
If \$sinSigma==0 Then Return 0; co-incident points
\$cosSigma = \$sinU1*\$sinU2 + \$cosU1*\$cosU2*\$cosLambda
\$sigma = atan(\$sinSigma/\$cosSigma)
\$sinAlpha = \$cosU1 * \$cosU2 * \$sinLambda / \$sinSigma
\$cosSqAlpha = 1 - \$sinAlpha*\$sinAlpha
\$cos2SigmaM = \$cosSigma - 2*\$sinU1*\$sinU2/\$cosSqAlpha
;If (IsNumber(\$cos2SigmaM)) Then \$cos2SigmaM = 0; equatorial line: cosSqAlpha=0 (§6)
\$C = \$f/16*\$cosSqAlpha*(4+\$f*(4-3*\$cosSqAlpha))
\$lambdaP = \$lambda;
\$lambda = \$L + (1-\$C) * \$f * \$sinAlpha * (\$sigma + \$C*\$sinSigma*(\$cos2SigmaM+\$C*\$cosSigma*(-1+2*\$cos2SigmaM*\$cos2SigmaM)))
WEnd
If (\$iterLimit==0) Then Return "NaN"; formula failed to converge

\$uSq = \$cosSqAlpha * (\$a*\$a - \$b_new*\$b_new) / (\$b_new*\$b_new)
\$A = 1 + \$uSq/16384*(4096+\$uSq*(-768+\$uSq*(320-175*\$uSq)))
\$B = \$uSq/1024 * (256+\$uSq*(-128+\$uSq*(74-47*\$uSq)))

\$innerb = -1 + 2 * \$cos2SigmaM * \$cos2SigmaM
\$innerc = -3 + 4 * \$sinSigma * \$sinSigma
\$innerd = -3 + 4 * \$cos2SigmaM * \$cos2SigmaM
\$innere = \$B/6 * \$cos2SigmaM
\$innera = (\$cosSigma * (\$innerb) - \$innere * (\$innerc) * (\$innerd))

\$innerz = \$cos2SigmaM + \$B/4
\$innery = \$B * \$sinSigma
\$deltaSigma = \$innery * (\$innerz * \$innera)

\$s = \$b_new * \$A * (\$sigma-\$deltaSigma)
Return \$s;
EndFunc

func distHaversine(\$lat1, \$lon1, \$lat2, \$lon2, \$radius = 6378135)
\$lat1 = DMStoDD(\$lat1)
\$lat2 = DMStoDD(\$lat2)
\$lon1 = DMStoDD(\$lon1)
\$lon2 = DMStoDD(\$lon2)

\$R = 6378135; // earth's mean radius in M (wikipedia)
\$dLat = (\$lat2-\$lat1)/57.2957795130823

\$a = sin(\$dLat/2) * sin(\$dLat/2) + cos(\$lat1) * cos(\$lat2) * sin(\$dLon/2) * sin(\$dLon/2)
\$c = 2 * atan(sqrt(\$a)/sqrt(1-\$a));
\$d = \$R * \$c;
return \$d;
EndFunc

; convert DMS to DD // convert from DMS to DD, decimal degrees = whole number of degrees, plus minutes divided by 60, plus seconds divided by 3600.
Func DMStoDD(\$DMS)
\$neg = False
If StringInStr(\$DMS, "S") Or StringInStr(\$DMS, "E") Or StringInStr(\$DMS, "-")Then
\$neg = True
EndIf
\$DMS = StringSplit(StringStripWS(StringRegExpReplace(\$DMS, "[^\d\s]", " "), 7), " ")

Switch \$DMS
Case 1
\$DD = \$DMS
Case 2
\$DD = \$DMS + \$DMS/60
Case 3
\$DD = \$DMS + \$DMS/60 + \$DMS/3600
Case 4
\$temp = "." & \$DMS
\$DMS += \$temp
\$DD = \$DMS + \$DMS/60 + \$DMS/3600
EndSwitch

;MsgBox (0, "", \$DD)
If \$neg Then
Return -\$DD
Else
Return \$DD
EndIf
EndFunc

#Region ### START Koda GUI section ### Form=
\$Form1 = GUICreate("DMS Points to Distance", 633, 150, 199, 123)
\$lat1 = GUICtrlCreateInput("53 09 02N", 72, 24, 217, 21)
\$lon1 = GUICtrlCreateInput("001 50 40W", 360, 24, 217, 21)
\$lat2 = GUICtrlCreateInput("52 12 19N", 72, 56, 217, 21)
\$lon2 = GUICtrlCreateInput("000 08 33W", 360, 56, 217, 21)
\$Label1 = GUICtrlCreateLabel("lat 1", 16, 24, 43, 17)
\$Label2 = GUICtrlCreateLabel("long 1", 304, 24, 43, 17)
\$Label3 = GUICtrlCreateLabel("lat 2", 16, 56, 43, 17)
\$Label4 = GUICtrlCreateLabel("long 2", 304, 56, 43, 17)
\$CalculateH = GUICtrlCreateButton("Calculate Distance With Haversine", 56, 96, 201, 20, 0)
\$CalculateV = GUICtrlCreateButton("Calculate Distance With Vincenty", 56, 120, 201, 20, 0)
\$Result = GUICtrlCreateInput("Results", 312, 104, 241, 21)
GUISetState(@SW_SHOW)
#EndRegion ### END Koda GUI section ###

While 1
\$nMsg = GUIGetMsg()
Switch \$nMsg
Case \$CalculateH
Case \$CalculateV
Case \$GUI_EVENT_CLOSE
Exit
EndSwitch
WEnd```

Let me know what you think!

ps it outputs in KM (Haversine) or M (Vincenty)

Edit Fixed: Vincenty function works kinda... the isNumber() was acting up as well as the 2 different \$b \$B variables. I split up some of the calculations now as well. It still acts up on long distances... any ideas??? Both output in meters now.

Edited by mikeytown2
##### Share on other sites

• 5 months later...

Decided to post my updated version since there seems to be some interest. Code below has no gui, also has some bonus functions.

```#include <INet.au3>

Global Const \$PI = 4 * ATan(1)
Global Const \$degToRad = \$PI / 180

\$lon1dd = -116
\$lat1dd = 33
\$lon2dd = -116.05
\$lat2dd = 33.05

MsgBox(0, "Test", DisplayAll(CalcAll(\$lat1dd, \$lon1dd, \$lat2dd, \$lon2dd)))

; Calculate geodesic distance (in m) between two points specified by latitude/longitude (in DD)
; using Vincenty inverse formula for ellipsoids
;Get 3 sides of triangle1, output in meters

;Get Heights, output in meters
\$heightB = USGS_DDtoElevation(\$latBdd, \$lonBdd)

;Get Triangles
\$triangle1 = SideAngleSide(\$sideA, 90, \$sideB) ; = [\$angleA, \$sideC, \$angleB]
\$triangle2 = SideAngleSide(\$sideC, 90, \$heightA - \$heightB) ; = [\$angleC, \$sideD, \$angleA], SideD is height adjusted distance

;Output Info
Dim \$TempArray = [\$latAdd, \$lonAdd, \$latBdd, \$lonBdd, \$sideA, \$sideB, \$sideC, \$heightA, \$heightB, \$triangle1, \$triangle1, \$triangle1, \$triangle2, \$triangle2, \$triangle2]
Return \$TempArray
EndFunc   ;==>CalcAll

Func DisplayAll(\$Inarray)
\$Out = "Starting Cords:" & @CRLF & _
"Latitude: " & DDtoDMS(\$Inarray) & " Or " & \$Inarray & @CRLF & _
"Longitude: " & DDtoDMS(\$Inarray) & " Or " & \$Inarray & @CRLF & _
"Elevation: " & \$Inarray & " Meters Or " & MetersToFeet(\$Inarray) & " Feet" & @CRLF & _
@CRLF & _
"Ending Cords:" & @CRLF & _
"Latitude: " & DDtoDMS(\$Inarray) & " Or " & \$Inarray & @CRLF & _
"Longitude: " & DDtoDMS(\$Inarray) & " Or " & \$Inarray & @CRLF & _
"Elevation: " & \$Inarray & " Meters Or " & MetersToFeet(\$Inarray) & " Feet" & @CRLF & _
@CRLF & _
"Distance North South: " & \$Inarray & " Meters" & @CRLF & _
"Distance East West: " & \$Inarray & " Meters" & @CRLF & _
"Distance: " & \$Inarray & " Meters" & @CRLF & _
"Less Accurate Distance (Vincenty Triangle): " & \$Inarray & " Meters" & @CRLF & _
"Less Accurate Distance (Haversine): " & distHaversine(\$Inarray, \$Inarray, \$Inarray, \$Inarray) & " Meters" & @CRLF & _
"Height Adjusted Distance: " & \$Inarray & " Meters" & @CRLF & _
"Elevation Difference: " & \$Inarray - \$Inarray & " Meters" & @CRLF & _
@CRLF & _
"Local Angle: " & \$Inarray & @CRLF & _
"Height Angle From End to Start Cords: " & \$Inarray & @CRLF & _
@CRLF & _
"Values Not Used:" & @CRLF & _
"Angle B Triangle 1 (Azimuth Calc): " & \$Inarray & @CRLF & _
"Angle C Triangle 2 (Height Calc): " & \$Inarray
Return \$Out
EndFunc   ;==>DisplayAll

Func distVincenty(\$lat1, \$lon1, \$lat2, \$lon2)
; WGS-84 ellipsiod
\$a = 6378137
\$b_new = 6356752.3142
\$f = 1 / 298.257223563

\$L = (\$lon2 - \$lon1) * \$degToRad
\$U1 = ATan((1 - \$f) * Tan(\$lat1 * \$degToRad))
\$U2 = ATan((1 - \$f) * Tan(\$lat2 * \$degToRad))
\$sinU1 = Sin(\$U1)
\$cosU1 = Cos(\$U1)
\$sinU2 = Sin(\$U2)
\$cosU2 = Cos(\$U2)

\$lambda = \$L
\$lambdaP = 2 * \$PI

\$iterLimit = 256
While (Abs(\$lambda - \$lambdaP) > 10 ^ - 12) And (\$iterLimit > 0)
\$iterLimit -= 1
\$sinLambda = Sin(\$lambda)
\$cosLambda = Cos(\$lambda)
\$sinSigma = Sqrt((\$cosU2 * \$sinLambda) * (\$cosU2 * \$sinLambda) + (\$cosU1 * \$sinU2 - \$sinU1 * \$cosU2 * \$cosLambda) * (\$cosU1 * \$sinU2 - \$sinU1 * \$cosU2 * \$cosLambda))
If \$sinSigma == 0 Then Return 0 ; co-incident points
\$cosSigma = \$sinU1 * \$sinU2 + \$cosU1 * \$cosU2 * \$cosLambda
\$sigma = ATan(\$sinSigma / \$cosSigma)
\$sinAlpha = \$cosU1 * \$cosU2 * \$sinLambda / \$sinSigma
\$cosSqAlpha = 1 - \$sinAlpha * \$sinAlpha
\$cos2SigmaM = \$cosSigma - 2 * \$sinU1 * \$sinU2 / \$cosSqAlpha
;If (IsNumber(\$cos2SigmaM)) Then \$cos2SigmaM = 0 ; equatorial line: cosSqAlpha=0 (§6)
\$C = \$f / 16 * \$cosSqAlpha * (4 + \$f * (4 - 3 * \$cosSqAlpha))
\$lambdaP = \$lambda;
\$lambda = \$L + (1 - \$C) * \$f * \$sinAlpha * (\$sigma + \$C * \$sinSigma * (\$cos2SigmaM + \$C * \$cosSigma * (-1 + 2 * \$cos2SigmaM * \$cos2SigmaM)))
WEnd
If (\$iterLimit == 0) Then Return "NaN" ; formula failed to converge

\$uSq = \$cosSqAlpha * (\$a * \$a - \$b_new * \$b_new) / (\$b_new * \$b_new)
\$a = 1 + \$uSq / 16384 * (4096 + \$uSq * (-768 + \$uSq * (320 - 175 * \$uSq)))
\$B = \$uSq / 1024 * (256 + \$uSq * (-128 + \$uSq * (74 - 47 * \$uSq)))

\$innerb = -1 + 2 * \$cos2SigmaM * \$cos2SigmaM
\$innerc = -3 + 4 * \$sinSigma * \$sinSigma
\$innerd = -3 + 4 * \$cos2SigmaM * \$cos2SigmaM
\$innere = \$B / 6 * \$cos2SigmaM
\$innera = (\$cosSigma * (\$innerb) - \$innere * (\$innerc) * (\$innerd))

\$innerz = \$cos2SigmaM + \$B / 4
\$innery = \$B * \$sinSigma
\$deltaSigma = \$innery * (\$innerz * \$innera)

\$s = \$b_new * \$a * (\$sigma - \$deltaSigma)
Return \$s
EndFunc   ;==>distVincenty

Func distHaversine(\$lat1, \$lon1, \$lat2, \$lon2)
\$R = 6371000
\$dLat = (\$lat2 - \$lat1) * \$degToRad
\$dLon = (\$lon2 - \$lon1) * \$degToRad

\$a = Sin(\$dLat / 2) * Sin(\$dLat / 2) + Cos(\$lat1) * Cos(\$lat2) * Sin(\$dLon / 2) * Sin(\$dLon / 2)
\$C = 2 * ATan(Sqrt(\$a) / Sqrt(1 - \$a));
\$d = \$R * \$C;
Return \$d;
EndFunc   ;==>distHaversine

Func USGS_DDtoElevation(\$lat, \$long)
\$url = "http://gisdata.usgs.gov/XMLWebServices/TNM_Elevation_Service.asmx/getElevation?X_Value=" & \$long & "&Y_Value=" & \$lat & "&Elevation_Units=FEET&Source_Layer=NED.CONUS_NED_13W&Elevation_Only=1"
\$url = _INetGetSource(\$url)
\$url = StringSplit(\$url, "Elevation", 1)
If \$url > 1 Then
\$url = StringTrimRight(StringTrimLeft(\$url, 4), 5)
Else
Return 0
EndIf
Return FeetToMeters(\$url)
EndFunc   ;==>USGS_DDtoElevation

;converts meters to feet
Func MetersToFeet(\$meters)
Return \$meters * 3.2808399
EndFunc   ;==>MetersToFeet

;converts feet to meters
Func FeetToMeters(\$feet)
Return \$feet * (1 / 3.2808399)
EndFunc   ;==>FeetToMeters

; convert DMS to DD // convert from DMS to DD, decimal degrees = whole number of degrees, plus minutes divided by 60, plus seconds divided by 3600.
Func DMStoDD(\$DMS)
\$neg = False
If StringInStr(\$DMS, "S") Or StringInStr(\$DMS, "E") Or StringInStr(\$DMS, "-") Then
\$neg = True
EndIf
\$DMS = StringSplit(StringStripWS(StringRegExpReplace(\$DMS, "[^\d\s]", " "), 7), " ")

Switch \$DMS
Case 1
\$DD = \$DMS
Case 2
\$DD = \$DMS + \$DMS / 60
Case 3
\$DD = \$DMS + \$DMS / 60 + \$DMS / 3600
Case 4
\$temp = "." & \$DMS
\$DMS += \$temp
\$DD = \$DMS + \$DMS / 60 + \$DMS / 3600
EndSwitch

If \$neg Then
Return -\$DD
Else
Return \$DD
EndIf
EndFunc   ;==>DMStoDD

; convert DD to DMS // convert from DD to DMS, decimal degrees = whole number of degrees, plus minutes divided by 60, plus seconds divided by 3600.
Func DDtoDMS(\$DD)
\$d = Abs(\$DD);  // (unsigned result ready for appending compass dir'n)
\$deg = Floor(\$d)
\$min = Floor((\$d - \$deg) * 60)
\$sec = Round((\$d - \$deg - \$min / 60) * 3600, 4)
If (\$deg < 100) Then
\$deg = '0' & \$deg ; if (deg<10) deg = '0' + deg;
EndIf
If (\$min < 10) Then
\$min = '0' & \$min;
EndIf
If (\$sec < 10) Then
\$sec = '0' & \$sec;
EndIf

Return \$deg & "° " & \$min & "' " & \$sec & '"'
EndFunc   ;==>DDtoDMS

;gives the other 2 sides and angle of the triangle
Func AngleSideAngle(\$angleA, \$sideA, \$angleB)
\$angleC = 180 - (\$angleA + \$angleB)
;ConsoleWrite("\$angleA=" & \$angleA & " \$angleB=" & \$angleB & " \$angleC=" & \$angleC & " \$sideA=" & \$sideA & " \$sideB=" & \$sideB & " \$sideC=" & \$sideC & @CRLF)
Dim \$TempArray = [\$sideB, \$angleC, \$sideC]
Return \$TempArray
EndFunc   ;==>AngleSideAngle

Func SideAngleSide(\$sideA, \$angleB, \$sideC)
\$sideB = Sqrt(\$sideA ^ 2 + \$sideC ^ 2 - 2 * \$sideA * \$sideC * Cos(\$angleB * \$degToRad))
\$angleA = ACos((\$sideA ^ 2 - \$sideC ^ 2 - \$sideB ^ 2) / - (2 * \$sideC * \$sideB)) / \$degToRad
If \$angleA == 180 Then
\$angleB = 0
\$sideB = Sqrt(\$sideA ^ 2 + \$sideC ^ 2 - 2 * \$sideA * \$sideC * Cos(\$angleB * \$degToRad))
EndIf
\$angleC = 180 - (\$angleB + \$angleA)
;ConsoleWrite("\$angleA=" & \$angleA & " \$angleB=" & \$angleB & " \$angleC=" & \$angleC & " \$sideA=" & \$sideA & " \$sideB=" & \$sideB & " \$sideC=" & \$sideC & @CRLF)
Dim \$TempArray = [\$angleA, \$sideB, \$angleC]
Return \$TempArray
EndFunc   ;==>SideAngleSide

Func BearingToAzimuth(\$B)
\$temp = StringSplit(\$B, " ")

If \$temp == 4 Then
\$temp = DMStoDD(\$temp & " " & \$temp)
\$temp = \$temp
EndIf

If StringInStr(\$temp, "N") > 0 Then
If StringInStr(\$temp, "E") > 0 Then
\$a = \$temp
Else
\$a = 360 - \$temp
EndIf
Else
If StringInStr(\$temp, "E") > 0 Then
\$a = 180 - \$temp
Else
\$a = 180 + \$temp
EndIf
EndIf
Return \$a
EndFunc   ;==>BearingToAzimuth```
Edited by mikeytown2
##### Share on other sites

Decided to post my updated version since there seems to be some interest. Code below has no gui, also has some bonus functions.

```#include <Math.au3>
#include <INet.au3>

Global Const \$PI = 4 * ATan(1)
Global Const \$degToRad = \$PI / 180

\$lon1dd = -116
\$lat1dd = 33
\$lon2dd = -116.05
\$lat2dd = 33.05

MsgBox(0, "Test", DisplayAll(CalcAll(\$lat1dd, \$lon1dd, \$lat2dd, \$lon2dd)))

; Calculate geodesic distance (in m) between two points specified by latitude/longitude (in DD)
; using Vincenty inverse formula for ellipsoids
;Get 3 sides of triangle1, output in meters

;Get Heights, output in meters
\$heightB = USGS_DDtoElevation(\$latBdd, \$lonBdd)

;Get Triangles
\$triangle1 = SideAngleSide(\$sideA, 90, \$sideB) ; = [\$angleA, \$sideC, \$angleB]
\$triangle2 = SideAngleSide(\$sideC, 90, \$heightA - \$heightB) ; = [\$angleC, \$sideD, \$angleA], SideD is height adjusted distance

;Output Info
Dim \$TempArray = [\$latAdd, \$lonAdd, \$latBdd, \$lonBdd, \$sideA, \$sideB, \$sideC, \$heightA, \$heightB, \$triangle1, \$triangle1, \$triangle1, \$triangle2, \$triangle2, \$triangle2]
Return \$TempArray
EndFunc   ;==>CalcAll

Func DisplayAll(\$Inarray)
\$Out = "Starting Cords:" & @CRLF & _
"Latitude: " & DDtoDMS(\$Inarray) & " Or " & \$Inarray & @CRLF & _
"Longitude: " & DDtoDMS(\$Inarray) & " Or " & \$Inarray & @CRLF & _
"Elevation: " & \$Inarray & " Meters Or " & MetersToFeet(\$Inarray) & " Feet" & @CRLF & _
@CRLF & _
"Ending Cords:" & @CRLF & _
"Latitude: " & DDtoDMS(\$Inarray) & " Or " & \$Inarray & @CRLF & _
"Longitude: " & DDtoDMS(\$Inarray) & " Or " & \$Inarray & @CRLF & _
"Elevation: " & \$Inarray & " Meters Or " & MetersToFeet(\$Inarray) & " Feet" & @CRLF & _
@CRLF & _
"Distance North South: " & \$Inarray & " Meters" & @CRLF & _
"Distance East West: " & \$Inarray & " Meters" & @CRLF & _
"Distance: " & \$Inarray & " Meters" & @CRLF & _
"Less Accurate Distance (Vincenty Triangle): " & \$Inarray & " Meters" & @CRLF & _
"Less Accurate Distance (Haversine): " & distHaversine(\$Inarray, \$Inarray, \$Inarray, \$Inarray) & " Meters" & @CRLF & _
"Height Adjusted Distance: " & \$Inarray & " Meters" & @CRLF & _
"Elevation Difference: " & \$Inarray - \$Inarray & " Meters" & @CRLF & _
@CRLF & _
"Local Angle: " & \$Inarray & @CRLF & _
"Height Angle From End to Start Cords: " & \$Inarray & @CRLF & _
@CRLF & _
"Values Not Used:" & @CRLF & _
"Angle B Triangle 1 (Azimuth Calc): " & \$Inarray & @CRLF & _
"Angle C Triangle 2 (Height Calc): " & \$Inarray
Return \$Out
EndFunc   ;==>DisplayAll

Func distVincenty(\$lat1, \$lon1, \$lat2, \$lon2)
; WGS-84 ellipsiod
\$a = 6378137
\$b_new = 6356752.3142
\$f = 1 / 298.257223563

\$L = (\$lon2 - \$lon1) * \$degToRad
\$U1 = ATan((1 - \$f) * Tan(\$lat1 * \$degToRad))
\$U2 = ATan((1 - \$f) * Tan(\$lat2 * \$degToRad))
\$sinU1 = Sin(\$U1)
\$cosU1 = Cos(\$U1)
\$sinU2 = Sin(\$U2)
\$cosU2 = Cos(\$U2)

\$lambda = \$L
\$lambdaP = 2 * \$PI

\$iterLimit = 256
While (Abs(\$lambda - \$lambdaP) > 10 ^ - 12) And (\$iterLimit > 0)
\$iterLimit -= 1
\$sinLambda = Sin(\$lambda)
\$cosLambda = Cos(\$lambda)
\$sinSigma = Sqrt((\$cosU2 * \$sinLambda) * (\$cosU2 * \$sinLambda) + (\$cosU1 * \$sinU2 - \$sinU1 * \$cosU2 * \$cosLambda) * (\$cosU1 * \$sinU2 - \$sinU1 * \$cosU2 * \$cosLambda))
If \$sinSigma == 0 Then Return 0 ; co-incident points
\$cosSigma = \$sinU1 * \$sinU2 + \$cosU1 * \$cosU2 * \$cosLambda
\$sigma = ATan(\$sinSigma / \$cosSigma)
\$sinAlpha = \$cosU1 * \$cosU2 * \$sinLambda / \$sinSigma
\$cosSqAlpha = 1 - \$sinAlpha * \$sinAlpha
\$cos2SigmaM = \$cosSigma - 2 * \$sinU1 * \$sinU2 / \$cosSqAlpha
;If (IsNumber(\$cos2SigmaM)) Then \$cos2SigmaM = 0 ; equatorial line: cosSqAlpha=0 (§6)
\$C = \$f / 16 * \$cosSqAlpha * (4 + \$f * (4 - 3 * \$cosSqAlpha))
\$lambdaP = \$lambda;
\$lambda = \$L + (1 - \$C) * \$f * \$sinAlpha * (\$sigma + \$C * \$sinSigma * (\$cos2SigmaM + \$C * \$cosSigma * (-1 + 2 * \$cos2SigmaM * \$cos2SigmaM)))
WEnd
If (\$iterLimit == 0) Then Return "NaN" ; formula failed to converge

\$uSq = \$cosSqAlpha * (\$a * \$a - \$b_new * \$b_new) / (\$b_new * \$b_new)
\$a = 1 + \$uSq / 16384 * (4096 + \$uSq * (-768 + \$uSq * (320 - 175 * \$uSq)))
\$B = \$uSq / 1024 * (256 + \$uSq * (-128 + \$uSq * (74 - 47 * \$uSq)))

\$innerb = -1 + 2 * \$cos2SigmaM * \$cos2SigmaM
\$innerc = -3 + 4 * \$sinSigma * \$sinSigma
\$innerd = -3 + 4 * \$cos2SigmaM * \$cos2SigmaM
\$innere = \$B / 6 * \$cos2SigmaM
\$innera = (\$cosSigma * (\$innerb) - \$innere * (\$innerc) * (\$innerd))

\$innerz = \$cos2SigmaM + \$B / 4
\$innery = \$B * \$sinSigma
\$deltaSigma = \$innery * (\$innerz * \$innera)

\$s = \$b_new * \$a * (\$sigma - \$deltaSigma)
Return \$s
EndFunc   ;==>distVincenty

Func distHaversine(\$lat1, \$lon1, \$lat2, \$lon2)
\$R = 6378135 ; // earth's mean radius in Meters (wikipedia)
\$dLat = (\$lat2 - \$lat1) * \$degToRad

\$a = Sin(\$dLat / 2) * Sin(\$dLat / 2) + Cos(\$lat1) * Cos(\$lat2) * Sin(\$dLon / 2) * Sin(\$dLon / 2)
\$C = 2 * ATan(Sqrt(\$a) / Sqrt(1 - \$a));
\$d = \$R * \$C;
Return \$d;
EndFunc   ;==>distHaversine

Func USGS_DDtoElevation(\$lat, \$long)
\$url = "http://gisdata.usgs.gov/XMLWebServices/TNM_Elevation_Service.asmx/getElevation?X_Value=" & \$long & "&Y_Value=" & \$lat & "&Elevation_Units=FEET&Source_Layer=NED.CONUS_NED_13W&Elevation_Only=1"
\$url = _INetGetSource(\$url)
\$url = StringSplit(\$url, "Elevation", 1)
\$url = StringTrimRight(StringTrimLeft(\$url, 4), 5)
Return FeetToMeters(\$url)
EndFunc   ;==>USGS_DDtoElevation

;converts meters to feet
Func MetersToFeet(\$meters)
Return \$meters * 3.2808399
EndFunc   ;==>MetersToFeet

;converts feet to meters
Func FeetToMeters(\$feet)
Return \$feet * (1 / 3.2808399)
EndFunc   ;==>FeetToMeters

; convert DMS to DD // convert from DMS to DD, decimal degrees = whole number of degrees, plus minutes divided by 60, plus seconds divided by 3600.
Func DMStoDD(\$DMS)
\$neg = False
If StringInStr(\$DMS, "S") Or StringInStr(\$DMS, "E") Or StringInStr(\$DMS, "-") Then
\$neg = True
EndIf
\$DMS = StringSplit(StringStripWS(StringRegExpReplace(\$DMS, "[^\d\s]", " "), 7), " ")

Switch \$DMS
Case 1
\$DD = \$DMS
Case 2
\$DD = \$DMS + \$DMS / 60
Case 3
\$DD = \$DMS + \$DMS / 60 + \$DMS / 3600
Case 4
\$temp = "." & \$DMS
\$DMS += \$temp
\$DD = \$DMS + \$DMS / 60 + \$DMS / 3600
EndSwitch

If \$neg Then
Return -\$DD
Else
Return \$DD
EndIf
EndFunc   ;==>DMStoDD

; convert DD to DMS // convert from DD to DMS, decimal degrees = whole number of degrees, plus minutes divided by 60, plus seconds divided by 3600.
Func DDtoDMS(\$DD)
\$d = Abs(\$DD);  // (unsigned result ready for appending compass dir'n)
\$deg = Floor(\$d)
\$min = Floor((\$d - \$deg) * 60)
\$sec = Round((\$d - \$deg - \$min / 60) * 3600, 4)
If (\$deg < 100) Then
\$deg = '0' & \$deg ; if (deg<10) deg = '0' + deg;
EndIf
If (\$min < 10) Then
\$min = '0' & \$min;
EndIf
If (\$sec < 10) Then
\$sec = '0' & \$sec;
EndIf

Return \$deg & "° " & \$min & "' " & \$sec & '"'
EndFunc   ;==>DDtoDMS

;gives the other 2 sides and angle of the triangle
Func AngleSideAngle(\$angleA, \$sideA, \$angleB)
\$angleC = 180 - (\$angleA + \$angleB)
;ConsoleWrite("\$angleA=" & \$angleA & " \$angleB=" & \$angleB & " \$angleC=" & \$angleC & " \$sideA=" & \$sideA & " \$sideB=" & \$sideB & " \$sideC=" & \$sideC & @CRLF)
Dim \$TempArray = [\$sideB, \$angleC, \$sideC]
Return \$TempArray
EndFunc   ;==>AngleSideAngle

Func SideAngleSide(\$sideA, \$angleB, \$sideC)
\$sideB = Sqrt(\$sideA ^ 2 + \$sideC ^ 2 - 2 * \$sideA * \$sideC * Cos(\$angleB * \$degToRad))
\$angleA = ACos((\$sideA ^ 2 - \$sideC ^ 2 - \$sideB ^ 2) / - (2 * \$sideC * \$sideB)) / \$degToRad
If \$angleA == 180 Then
\$angleB = 0
\$sideB = Sqrt(\$sideA ^ 2 + \$sideC ^ 2 - 2 * \$sideA * \$sideC * Cos(\$angleB * \$degToRad))
EndIf
\$angleC = 180 - (\$angleB + \$angleA)
;ConsoleWrite("\$angleA=" & \$angleA & " \$angleB=" & \$angleB & " \$angleC=" & \$angleC & " \$sideA=" & \$sideA & " \$sideB=" & \$sideB & " \$sideC=" & \$sideC & @CRLF)
Dim \$TempArray = [\$angleA, \$sideB, \$angleC]
Return \$TempArray
EndFunc   ;==>SideAngleSide

Func BearingToAzimuth(\$B)
\$temp = StringSplit(\$B, " ")

If \$temp == 4 Then
\$temp = DMStoDD(\$temp & " " & \$temp)
\$temp = \$temp
EndIf

If StringInStr(\$temp, "N") > 0 Then
If StringInStr(\$temp, "E") > 0 Then
\$a = \$temp
Else
\$a = 360 - \$temp
EndIf
Else
If StringInStr(\$temp, "E") > 0 Then
\$a = 180 - \$temp
Else
\$a = 180 + \$temp
EndIf
EndIf
Return \$a
EndFunc   ;==>BearingToAzimuth```
Looks good. Vincenty has been giving me a problem so this is a definite improvement. What have your tests shown when working with different hemispheres?

George

Question about decompiling code? Read the decompiling FAQ and don't bother posting the question in the forums.

Be sure to read and follow the forum rules. -AKA the AutoIt Reading and Comprehension Skills test.***

The PCRE (Regular Expression) ToolKit for AutoIT - (Updated Oct 20, 2011 ver:3.0.1.13) - Please update your current version before filing any bug reports. The installer now includes both 32 and 64 bit versions. No change in version number.

Visit my Blog .. currently not active but it will soon be resplendent with news and views. Also please remove any links you may have to my website. it is soon to be closed and replaced with something else.

"Old age and treachery will always overcome youth and skill!"

##### Share on other sites

Looks good. Vincenty has been giving me a problem so this is a definite improvement. What have your tests shown when working with different hemispheres?

The missing atan2() function seems to be the problem when working in different hemispheres. I didn't need to calculate vast distances so i didn't worry too much.

The National Elevation Dataset (NED), in some areas, has a resolution up to 1/9 arc second. I use 1/3 because of full coverage.

##### Share on other sites

There is definitly a problem in your haversine formula. I took your Vincenty and added it to my functions (without taking elevation into account)

Using

Lat 1 = 61

Lon 1 = 149

Lat 2 = 49

Lon 2 = 123

The Results are

My Haversine = 2100.30691449812 Km

Which shows a difference of 7 km

George

Question about decompiling code? Read the decompiling FAQ and don't bother posting the question in the forums.

Be sure to read and follow the forum rules. -AKA the AutoIt Reading and Comprehension Skills test.***

The PCRE (Regular Expression) ToolKit for AutoIT - (Updated Oct 20, 2011 ver:3.0.1.13) - Please update your current version before filing any bug reports. The installer now includes both 32 and 64 bit versions. No change in version number.

Visit my Blog .. currently not active but it will soon be resplendent with news and views. Also please remove any links you may have to my website. it is soon to be closed and replaced with something else.

"Old age and treachery will always overcome youth and skill!"

##### Share on other sites

When comparing both haversine's to google earth, we are off!

Point A

Lat 33

Lon -116

Point B

Lat 33.005 Or 33° 0'18.00"N

Lon -116.005 Or 116° 0'18.00"W

My Vincenty: 725.42 Meters

My Haversine: 556.60 Meters

My Error has to do with the radius. if i change this to this

```Func distHaversine(\$lat1, \$lon1, \$lat2, \$lon2)
\$R = 6378135 ; // earth's mean radius in Meters (wikipedia)
\$dLat = (\$lat2 - \$lat1) * \$degToRad

\$a = Sin(\$dLat / 2) * Sin(\$dLat / 2) + Cos(\$lat1) * Cos(\$lat2) * Sin(\$dLon / 2) * Sin(\$dLon / 2)
\$C = 2 * ATan(Sqrt(\$a) / Sqrt(1 - \$a));
\$d = \$R * \$C;
Return \$d;
EndFunc   ;==>distHaversineoÝ÷ Ù:+ºÚ"µÍ[ÈÝ]Ú[J    ÌÍÛ]K    ÌÍÛÛK   ÌÍÛ] ÌÍÛÛB   ÌÍÔH
ÍÌL   ÌÍÙ]H
ÌÍÛ]H    ÌÍÛ]JH
ÌÍÙYÕÔY    ÌÍÙÛHÔYX[  ÌÍÛÛH   ÌÍÛÛJB  ÌÍÛ]HHÔYX[  ÌÍÛ]JB   ÌÍÛ]HÔYX[   ÌÍÛ]B    ÌÍØHHÚ[ ÌÍÙ]ÈH
Ú[ ÌÍÙ]ÈH
ÈÛÜÊ    ÌÍÛ]JH
ÛÜÊ  ÌÍÛ]H
Ú[ ÌÍÙÛÈH
Ú[ ÌÍÙÛÈB ÌÍÐÈH
U[Ü
ÌÍØJHÈÜ
HH  ÌÍØJJNÂ ÌÍÙH ÌÍÔ
ÌÍÐÎÂ] ÌÍÙÂ[[ÈÏOIÝÙÝ]Ú[```

I get a distance of 725.61 Meters

##### Share on other sites

Using

Lat 1 = 61

Lon 1 = 149

Lat 2 = 49

Lon 2 = 123

The Results are

My Haversine = 2100.30691449812 Km

Which shows a difference of 7 km

My New Haversine: 2104094.37 Meters

According to this website

http://www.movable-type.co.uk/scripts/latlong-vincenty.html

Vincenty: 2109284.63 Meters

Edited by mikeytown2
##### Share on other sites

New results after tuning my Haversine for accurracy rounded to 3 dp in Km

Tuned Hav 2100.408

Edit; Assuming you are still including Elevation in your calculation, what do you get when you use this for the haversine

```Func _Haversine(\$Lat1, \$Lon1, \$Lat2, \$Lon2);  <<== Get the distance between two points on the globe.
\$sinl1 = Sin(\$Lat1);
;; Diameter of Earth at the equator (referenced at mean sea level) is 7926.3811838861 miles so we do this.....
\$Diam = 7926.3811838861
\$D = (\$Diam - (26 * \$sinl1)) * ASin(_Min(1, 0.707106781186548 * Sqrt((1 - (Sin(\$Lat2) * \$sinl1) - Cos(\$Lat1) * Cos(\$Lat2) * Cos(\$Lon2 - \$Lon1)))));
Return \$D;
EndFunc  ;==>_Haversine```
Edited by GEOSoft

George

Question about decompiling code? Read the decompiling FAQ and don't bother posting the question in the forums.

Be sure to read and follow the forum rules. -AKA the AutoIt Reading and Comprehension Skills test.***

The PCRE (Regular Expression) ToolKit for AutoIT - (Updated Oct 20, 2011 ver:3.0.1.13) - Please update your current version before filing any bug reports. The installer now includes both 32 and 64 bit versions. No change in version number.

Visit my Blog .. currently not active but it will soon be resplendent with news and views. Also please remove any links you may have to my website. it is soon to be closed and replaced with something else.

"Old age and treachery will always overcome youth and skill!"

##### Share on other sites

In my code i have these outputs

```"Distance: " & \$Inarray & " Meters" & @CRLF & _
"Height Adjusted Distance: " & \$Inarray & " Meters" & @CRLF & _```

The output from my haversine/vincenty equations do not take elevation into account. I do that after the fact using trig. Because the distance is so great, the angle is so small that the height adjusted distance is less then 0.01mm with your given coordinates.

In short your haversine is way off when you look at google earth's output. I'm guessing that google earth uses something like PGM2000A, where I use WGS84.

##### Share on other sites

In my code i have these outputs

```"Distance: " & \$Inarray & " Meters" & @CRLF & _
"Height Adjusted Distance: " & \$Inarray & " Meters" & @CRLF & _```

The output from my haversine/vincenty equations do not take elevation into account. I do that after the fact using trig. Because the distance is so great, the angle is so small that the height adjusted distance is less then 0.01mm with your given coordinates.

In short your haversine is way off when you look at google earth's output. I'm guessing that google earth uses something like PGM2000A, where I use WGS84.

I'm having a problem with the vincenty when getting the distance from London to Brisbane with

Lat 1 = 51 32

Lon 1 = 0 5

Lat 2 = 27 29 S

Lon 2 = 153 8 E

What do you get with those coordinates?

George

Question about decompiling code? Read the decompiling FAQ and don't bother posting the question in the forums.

Be sure to read and follow the forum rules. -AKA the AutoIt Reading and Comprehension Skills test.***

The PCRE (Regular Expression) ToolKit for AutoIT - (Updated Oct 20, 2011 ver:3.0.1.13) - Please update your current version before filing any bug reports. The installer now includes both 32 and 64 bit versions. No change in version number.

Visit my Blog .. currently not active but it will soon be resplendent with news and views. Also please remove any links you may have to my website. it is soon to be closed and replaced with something else.

"Old age and treachery will always overcome youth and skill!"

##### Share on other sites

Using http://www.movable-type.co.uk/scripts/latlong-vincenty.html (has the atan2() function)

16,524,957.852 Meters

16,497,614.23 Meters

Using Autoit Programs

My Vincenty: -3,495,508.86 Meters (Error because of missing atan2() function)

My Haversine: 16,520,009.62 Meters

##### Share on other sites

Using http://www.movable-type.co.uk/scripts/latlong-vincenty.html (has the atan2() function)

16,524,957.852 Meters

16,497,614.23 Meters

Using Autoit Programs

My Vincenty: -3,495,508.86 Meters (Error because of missing atan2() function)

My Haversine: 16,520,009.62 Meters

Well your vincenty is closer now anyway I used the following

\$sigma = _Atan2(\$sinSigma, \$cosSigma)

in place of \$sigma = ATan(\$sinSigma / \$cosSigma)

Results in kilometers

My Haversine = 16395.6805508484

That will still land an airplane approx 77 kilometers shy of the runway. I still want to test it some more. I have not gone back and checked the results when points A and B are in the same hemisphere.

Edit: Readings are kilometers NOT meters as first stated

I think we may have to use Abs() in there someplace.

Edited by GEOSoft

George

Question about decompiling code? Read the decompiling FAQ and don't bother posting the question in the forums.

Be sure to read and follow the forum rules. -AKA the AutoIt Reading and Comprehension Skills test.***

The PCRE (Regular Expression) ToolKit for AutoIT - (Updated Oct 20, 2011 ver:3.0.1.13) - Please update your current version before filing any bug reports. The installer now includes both 32 and 64 bit versions. No change in version number.

Visit my Blog .. currently not active but it will soon be resplendent with news and views. Also please remove any links you may have to my website. it is soon to be closed and replaced with something else.

"Old age and treachery will always overcome youth and skill!"

##### Share on other sites

• 6 years later...

MikeyTown,

Thank you for this bit of code. I have recently started using this code for computing short distances (<3km).

The code for retrieving the elevation data was out of date and I have updated it to it's new source.

I thought I would share it so that any one else wanting to use this code could replace the below function and have it work properly. [uSGS_DDtoElevation(\$lat, \$long)]

```Func USGS_DDtoElevation(\$lat, \$long)
;   ConsoleWrite('\$lat: ' & \$lat & ' \$long: ' & \$long & '     ' & @ScriptLineNumber & @CRLF)
;    \$url = "http://gisdata.usgs.gov/XMLWebServices/TNM_Elevation_Service.asmx/getElevation?X_Value=" & \$long & "&Y_Value=" & \$lat & "&Elevation_Units=FEET&Source_Layer=NED.CONUS_NED_13W&Elevation_Only=1"
;    \$url = _INetGetSource(\$url)
;    \$url = StringSplit(\$url, "Elevation", 1)
;    If \$url > 1 Then
;    \$url = StringTrimRight(StringTrimLeft(\$url, 4), 5)
;    Else
;        Return 0
;    EndIf
;    Return FeetToMeters(\$url)
\$url = "http://ned.usgs.gov/epqs/pqs.php?x=" & \$long & "&y=" & \$lat & "units=Meters&output=xml"
\$url = _INetGetSource(\$url, TRUE)
\$url = StringSplit(\$url, "<Elevation>", 1)
If \$url > 1 Then
Return StringLeft(\$url,  StringInStr(\$url, '</Elevation>') - 1)
Else
Return 0
EndIf
EndFunc   ;==>USGS_DDtoElevation```
Edited by Shane0000
##### Share on other sites

there has to be a better way by now.

```,-. .--. ________ .-. .-. ,---. ,-. .-. .-. .-.
|(| / /\ \ |\ /| |__ __||| | | || .-' | |/ / \ \_/ )/
(_) / /__\ \ |(\ / | )| | | `-' | | `-. | | / __ \ (_)
| | | __ | (_)\/ | (_) | | .-. | | .-' | | \ |__| ) (
| | | | |)| | \ / | | | | | |)| | `--. | |) \ | |
`-' |_| (_) | |\/| | `-' /( (_)/( __.' |((_)-' /(_|
'-' '-' (__) (__) (_) (__)```

##### Share on other sites

i looked, and am still looking.  We might have to make a new thread for surprisingly useful necros.

```,-. .--. ________ .-. .-. ,---. ,-. .-. .-. .-.
|(| / /\ \ |\ /| |__ __||| | | || .-' | |/ / \ \_/ )/
(_) / /__\ \ |(\ / | )| | | `-' | | `-. | | / __ \ (_)
| | | __ | (_)\/ | (_) | | .-. | | .-' | | \ |__| ) (
| | | | |)| | \ / | | | | | |)| | `--. | |) \ | |
`-' |_| (_) | |\/| | `-' /( (_)/( __.' |((_)-' /(_|
'-' '-' (__) (__) (_) (__)```

##### Share on other sites

• 2 years later...

This was great. Thanks. All I needed was the distVincentry function with no other Includes. Very well written, and it was spot on when compared to a couple other sources such as an onilne version at  http://boulter.com/gps/distance/ and an Excel version at http://bluemm.blogspot.com/2007/01/excel-formula-to-calculate-distance.html, and of course Google Maps. Thanks!

Edited by Decibel

## Create an account

Register a new account

• ### Recently Browsing   0 members

×

• Wiki

• Back

• #### Beta

• Git
• FAQ
×
• Create New...