GPS Distance Calculations

[mikeytown2]

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

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

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).

expandcollapse popup

#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 
    
    $a = $radius
    
    $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
    $dLon = _Radian($lon2-$lon1)
    $lat1 = _Radian($lat1)
    $lat2 = _Radian($lat2)
    
    $R = $radius

    $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[0]
        Case 1
            $DD = $DMS[1]
        Case 2
            $DD = $DMS[1] + $DMS[2]/60
        Case 3
            $DD = $DMS[1] + $DMS[2]/60 + $DMS[3]/3600
        Case 4
            $temp = "." & $DMS[4]
            $DMS[3] += $temp
            $DD = $DMS[1] + $DMS[2]/60 + $DMS[3]/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
            GUICtrlSetData($Result, distHaversine(GUICtrlRead($lat1), GUICtrlRead($lon1), GUICtrlRead($lat2), GUICtrlRead($lon2))) 
        Case $CalculateV
            GUICtrlSetData($Result, distVincenty(GUICtrlRead($lat1), GUICtrlRead($lon1), GUICtrlRead($lat2), GUICtrlRead($lon2))) 
        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 December 12, 2007 by mikeytown2

[mikeytown2]

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

expandcollapse popup

#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
Func CalcAll($latAdd, $lonAdd, $latBdd, $lonBdd)
    ;Get 3 sides of triangle1, output in meters
    $sideA = distVincenty($latAdd, $lonAdd, $latBdd, $lonAdd)
    $sideB = distVincenty($latAdd, $lonAdd, $latAdd, $lonBdd)
    $sideC = distVincenty($latAdd, $lonAdd, $latBdd, $lonBdd)
    
    ;Get Heights, output in meters
    $heightA = USGS_DDtoElevation($latAdd, $lonAdd)
    $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[15] = [$latAdd, $lonAdd, $latBdd, $lonBdd, $sideA, $sideB, $sideC, $heightA, $heightB, $triangle1[0], $triangle1[1], $triangle1[2], $triangle2[0], $triangle2[1], $triangle2[2]]
    Return $TempArray
EndFunc   ;==>CalcAll

Func DisplayAll($Inarray)
    $Out = "Starting Cords:" & @CRLF & _
            "Latitude: " & DDtoDMS($Inarray[0]) & " Or " & $Inarray[0] & @CRLF & _
            "Longitude: " & DDtoDMS($Inarray[1]) & " Or " & $Inarray[1] & @CRLF & _
            "Elevation: " & $Inarray[7] & " Meters Or " & MetersToFeet($Inarray[7]) & " Feet" & @CRLF & _
            @CRLF & _
            "Ending Cords:" & @CRLF & _
            "Latitude: " & DDtoDMS($Inarray[2]) & " Or " & $Inarray[2] & @CRLF & _
            "Longitude: " & DDtoDMS($Inarray[3]) & " Or " & $Inarray[3] & @CRLF & _
            "Elevation: " & $Inarray[8] & " Meters Or " & MetersToFeet($Inarray[8]) & " Feet" & @CRLF & _
            @CRLF & _
            "Distance North South: " & $Inarray[4] & " Meters" & @CRLF & _
            "Distance East West: " & $Inarray[5] & " Meters" & @CRLF & _
            "Distance: " & $Inarray[6] & " Meters" & @CRLF & _
            "Less Accurate Distance (Vincenty Triangle): " & $Inarray[10] & " Meters" & @CRLF & _
            "Less Accurate Distance (Haversine): " & distHaversine($Inarray[0], $Inarray[1], $Inarray[2], $Inarray[3]) & " Meters" & @CRLF & _
            "Height Adjusted Distance: " & $Inarray[13] & " Meters" & @CRLF & _
            "Elevation Difference: " & $Inarray[7] - $Inarray[8] & " Meters" & @CRLF & _
            @CRLF & _
            "Local Angle: " & $Inarray[9] & @CRLF & _
            "Height Angle From End to Start Cords: " & $Inarray[14] & @CRLF & _
            @CRLF & _
            "Values Not Used:" & @CRLF & _
            "Angle B Triangle 1 (Azimuth Calc): " & $Inarray[11] & @CRLF & _
            "Angle C Triangle 2 (Height Calc): " & $Inarray[12]
    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
    $lat1 = ($lat1) * $degToRad
    $lat2 = ($lat2) * $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[0] > 1 Then
    $url = StringTrimRight(StringTrimLeft($url[5], 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[0]
        Case 1
            $DD = $DMS[1]
        Case 2
            $DD = $DMS[1] + $DMS[2] / 60
        Case 3
            $DD = $DMS[1] + $DMS[2] / 60 + $DMS[3] / 3600
        Case 4
            $temp = "." & $DMS[4]
            $DMS[3] += $temp
            $DD = $DMS[1] + $DMS[2] / 60 + $DMS[3] / 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)
    ;// add leading zeros if required
    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)
    $sideB = $sideA * Sin($angleB * $degToRad) / Sin($angleA * $degToRad)
    $sideC = $sideA * Sin($angleC * $degToRad) / Sin($angleA * $degToRad)
    ;ConsoleWrite("$angleA=" & $angleA & " $angleB=" & $angleB & " $angleC=" & $angleC & " $sideA=" & $sideA & " $sideB=" & $sideB & " $sideC=" & $sideC & @CRLF)
    Dim $TempArray[3] = [$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[3] = [$angleA, $sideB, $angleC]
    Return $TempArray
EndFunc   ;==>SideAngleSide

Func BearingToAzimuth($B)
    $temp = StringSplit($B, " ")
    
    If $temp[0] == 4 Then
        $temp[2] = DMStoDD($temp[2] & " " & $temp[3])
        $temp[3] = $temp[4]
    EndIf
    
    If StringInStr($temp[1], "N") > 0 Then
        If StringInStr($temp[3], "E") > 0 Then
            $a = $temp[2]
        Else
            $a = 360 - $temp[2]
        EndIf
    Else
        If StringInStr($temp[3], "E") > 0 Then
            $a = 180 - $temp[2]
        Else
            $a = 180 + $temp[2]
        EndIf
    EndIf
    Return $a
EndFunc   ;==>BearingToAzimuth

Edited December 12, 2007 by mikeytown2

[GEOSoft]

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

expandcollapse popup

#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
Func CalcAll($latAdd, $lonAdd, $latBdd, $lonBdd)
    ;Get 3 sides of triangle1, output in meters
    $sideA = distVincenty($latAdd, $lonAdd, $latBdd, $lonAdd)
    $sideB = distVincenty($latAdd, $lonAdd, $latAdd, $lonBdd)
    $sideC = distVincenty($latAdd, $lonAdd, $latBdd, $lonBdd)
    
    ;Get Heights, output in meters
    $heightA = USGS_DDtoElevation($latAdd, $lonAdd)
    $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[15] = [$latAdd, $lonAdd, $latBdd, $lonBdd, $sideA, $sideB, $sideC, $heightA, $heightB, $triangle1[0], $triangle1[1], $triangle1[2], $triangle2[0], $triangle2[1], $triangle2[2]]
    Return $TempArray
EndFunc   ;==>CalcAll

Func DisplayAll($Inarray)
    $Out = "Starting Cords:" & @CRLF & _
            "Latitude: " & DDtoDMS($Inarray[0]) & " Or " & $Inarray[0] & @CRLF & _
            "Longitude: " & DDtoDMS($Inarray[1]) & " Or " & $Inarray[1] & @CRLF & _
            "Elevation: " & $Inarray[7] & " Meters Or " & MetersToFeet($Inarray[7]) & " Feet" & @CRLF & _
            @CRLF & _
            "Ending Cords:" & @CRLF & _
            "Latitude: " & DDtoDMS($Inarray[2]) & " Or " & $Inarray[2] & @CRLF & _
            "Longitude: " & DDtoDMS($Inarray[3]) & " Or " & $Inarray[3] & @CRLF & _
            "Elevation: " & $Inarray[8] & " Meters Or " & MetersToFeet($Inarray[8]) & " Feet" & @CRLF & _
            @CRLF & _
            "Distance North South: " & $Inarray[4] & " Meters" & @CRLF & _
            "Distance East West: " & $Inarray[5] & " Meters" & @CRLF & _
            "Distance: " & $Inarray[6] & " Meters" & @CRLF & _
            "Less Accurate Distance (Vincenty Triangle): " & $Inarray[10] & " Meters" & @CRLF & _
            "Less Accurate Distance (Haversine): " & distHaversine($Inarray[0], $Inarray[1], $Inarray[2], $Inarray[3]) & " Meters" & @CRLF & _
            "Height Adjusted Distance: " & $Inarray[13] & " Meters" & @CRLF & _
            "Elevation Difference: " & $Inarray[7] - $Inarray[8] & " Meters" & @CRLF & _
            @CRLF & _
            "Local Angle: " & $Inarray[9] & @CRLF & _
            "Height Angle From End to Start Cords: " & $Inarray[14] & @CRLF & _
            @CRLF & _
            "Values Not Used:" & @CRLF & _
            "Angle B Triangle 1 (Azimuth Calc): " & $Inarray[11] & @CRLF & _
            "Angle C Triangle 2 (Height Calc): " & $Inarray[12]
    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
    $dLon = _Radian($lon2 - $lon1)
    $lat1 = _Radian($lat1)
    $lat2 = _Radian($lat2)

    $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[5], 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[0]
        Case 1
            $DD = $DMS[1]
        Case 2
            $DD = $DMS[1] + $DMS[2] / 60
        Case 3
            $DD = $DMS[1] + $DMS[2] / 60 + $DMS[3] / 3600
        Case 4
            $temp = "." & $DMS[4]
            $DMS[3] += $temp
            $DD = $DMS[1] + $DMS[2] / 60 + $DMS[3] / 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)
    ;// add leading zeros if required
    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)
    $sideB = $sideA * Sin($angleB * $degToRad) / Sin($angleA * $degToRad)
    $sideC = $sideA * Sin($angleC * $degToRad) / Sin($angleA * $degToRad)
    ;ConsoleWrite("$angleA=" & $angleA & " $angleB=" & $angleB & " $angleC=" & $angleC & " $sideA=" & $sideA & " $sideB=" & $sideB & " $sideC=" & $sideC & @CRLF)
    Dim $TempArray[3] = [$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[3] = [$angleA, $sideB, $angleC]
    Return $TempArray
EndFunc   ;==>SideAngleSide

Func BearingToAzimuth($B)
    $temp = StringSplit($B, " ")
    
    If $temp[0] == 4 Then
        $temp[2] = DMStoDD($temp[2] & " " & $temp[3])
        $temp[3] = $temp[4]
    EndIf
    
    If StringInStr($temp[1], "N") > 0 Then
        If StringInStr($temp[3], "E") > 0 Then
            $a = $temp[2]
        Else
            $a = 360 - $temp[2]
        EndIf
    Else
        If StringInStr($temp[3], "E") > 0 Then
            $a = 180 - $temp[2]
        Else
            $a = 180 + $temp[2]
        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?

[mikeytown2]

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.

[GEOSoft]

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

Your Vincenty = 2107.2106102474 km

Which shows a difference of 7 km

[mikeytown2]

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

Google Earth: 724.33 Meters

My Vincenty: 725.42 Meters

My Haversine: 556.60 Meters

Your Haversine: 12082.64 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
    $dLon = _Radian($lon2 - $lon1)
    $lat1 = _Radian($lat1)
    $lat2 = _Radian($lat2)

    $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

[mikeytown2]

Using

Lat 1 = 61

Lon 1 = 149

Lat 2 = 49

Lon 2 = 123

The Results are

My Haversine = 2100.30691449812 Km

Your Vincenty = 2107.2106102474 km

Which shows a difference of 7 km

Google Earth: 2106927.71 Meters

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 December 12, 2007 by mikeytown2

[GEOSoft]

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

Tuned Hav 2100.408

Your Vinc 2107.211

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.
    $Lat1 = _Radian($Lat1);
    $sinl1 = Sin($Lat1);
    $Lat2 = _Radian($Lat2);
    $Lon1 = _Radian($Lon1);
    $Lon2 = _Radian($Lon2);
;; 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 December 12, 2007 by GEOSoft

[mikeytown2]

In my code i have these outputs

"Distance: " & $Inarray[6] & " Meters" & @CRLF & _
"Height Adjusted Distance: " & $Inarray[13] & " 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.

[GEOSoft]

In my code i have these outputs

"Distance: " & $Inarray[6] & " Meters" & @CRLF & _
"Height Adjusted Distance: " & $Inarray[13] & " 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?

[mikeytown2]

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

16,524,957.852 Meters

Using Google Earth

16,497,614.23 Meters

Using Autoit Programs

Your Haversine: 0 M

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

My Haversine: 16,520,009.62 Meters

[GEOSoft]

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

16,524,957.852 Meters

Using Google Earth

16,497,614.23 Meters

Using Autoit Programs

Your Haversine: 0 M

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

Your Modified Vincenty = 16420.2353608579

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 December 12, 2007 by GEOSoft

[Shane0000]

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[0] > 1 Then
;    $url = StringTrimRight(StringTrimLeft($url[5], 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[0] > 1 Then
        Return StringLeft($url[2],  StringInStr($url[2], '</Elevation>') - 1)
    Else
        Return 0
    EndIf
EndFunc   ;==>USGS_DDtoElevation

Edited September 29, 2014 by Shane0000

[iamtheky]

there has to be a better way by now.

[Shane0000]

Im open to suggestions

[iamtheky]

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

[Decibel]

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 April 11, 2017 by Decibel
Citing Google Maps comparison.