_Geo_Diff - Exakte Berechnung der Distanz zweier Punkte auf der Erdoberfläche

**minx**

Titel sagt viel, das Skript den Rest . Beispiel ist mit dabei.

Spoiler anzeigen

[autoit]

; Paris Observatory
$geoParis = 'E-2-20-14-N-48-50-11'

[/autoit] [autoit][/autoit] [autoit]

; US Naval Observatory
$geoUSA = 'W-77-3-55.5-N-38-55-17'

[/autoit] [autoit][/autoit] [autoit]

ConsoleWrite(@LF)
ConsoleWrite("Distance between the two given locations is " & _Geo_Diff($geoParis, $geoUSA) & " ± 0.05 km!" & @LF)
ConsoleWrite(@LF)

[/autoit] [autoit][/autoit] [autoit]

; #FUNCTION# ================================================================================
; Name ..........: _Geo_Diff
; Description ...: Compute the gedesic distance between two Earth surface
; ...............: coordinates to an accuracy of about ±50 meters. To get
; ...............: this degree of accuracy, this function takes into account
; ...............: the spheroidal flattening factor of the Earth rather than
; ...............: assuming that the Earth is a perfect sphere.
; Syntax ........: _Geo_Diff($geoStart, $geoFinish[, $iDec = 2[, $iUnit = 0]])
; Parameters ....: $geoStart - Location Coordinates in GeoData form (see Example)
; $geoFinish - Location Coordinates in GeoData form (see Example)
; $iDec - [optional] Number of decimal places. Default is 2.
; $iUnit - [optional] The unit to use. Default is 0.
; 0 = km (± 0.05)
; 1 = mi (± 0.03)
; 2 = nmi (± 0.03)
; Return values .: Distance in the given unit, rounded to the given number of decimal places.
; Author ........: minx
; Example .......: s. above
; ===========================================================================================
Func _Geo_Diff($geoStart, $geoFinish, $iDec = 2, $iUnit = 0)
Local $aData[4]
Local $fTemp, $aTemp

[/autoit] [autoit][/autoit] [autoit]

$aTemp = StringSplit($geoStart, "-", 3)
$aData[0] = ($aTemp[1] * 3600 + $aTemp[2] * 60 + $aTemp[3]) / 3600
If $aTemp[0] = "E" Then $aData[0] *= -1
$aData[1] = ($aTemp[5] * 3600 + $aTemp[6] * 60 + $aTemp[7]) / 3600
If $aTemp[0] = "S" Then $aData[0] *= -1

[/autoit] [autoit][/autoit] [autoit]

$aTemp = StringSplit($geoFinish, "-", 3)
$aData[2] = ($aTemp[1] * 3600 + $aTemp[2] * 60 + $aTemp[3]) / 3600
If $aTemp[0] = "E" Then $aData[0] *= -1
$aData[3] = ($aTemp[5] * 3600 + $aTemp[6] * 60 + $aTemp[7]) / 3600
If $aTemp[0] = "S" Then $aData[0] *= -1

[/autoit] [autoit][/autoit] [autoit]

Local $ff = (1 / 298.257)

[/autoit] [autoit][/autoit] [autoit]

$UF = 1
If $iUnit = 1 Then $UF = 1.609344
If $iUnit = 2 Then $UF = 1.852

[/autoit] [autoit][/autoit] [autoit]

; auxiliary angles
$F = ($aData[1] + $aData[3]) / 2
$G = ($aData[1] - $aData[3]) / 2
$L = ($aData[0] - $aData[2]) / 2

[/autoit] [autoit][/autoit] [autoit]

; (co)sines of auxliliary angles
$SG = Sin($G * ATan(1) / 45)
$CG = Cos($G * ATan(1) / 45)
$SF = Sin($F * ATan(1) / 45)
$CF = Cos($F * ATan(1) / 45)
$SL = Sin($L * ATan(1) / 45)
$CL = Cos($L * ATan(1) / 45)
$S = ($SG * $CL) ^ 2 + ($CF * $SL) ^ 2
$C = ($CG * $CL) ^ 2 + ($SF * $SL) ^ 2
$O = ATan(Sqrt($S / $C))
$R = Sqrt($S * $C) / $O
$D = 2 * $O * 6378.14
$H1 = (3 * $R - 1) / (2 * $C)
$H2 = (3 * $R + 1) / (2 * $S)

[/autoit] [autoit][/autoit] [autoit]

; distance and angle between start and finish
; on the surface of a geodesic sphereoid
Return Round($D * ((($SF * $CG) ^ 2 * $H1 * $ff + 1) - (($CF * $SG) ^ 2 * $H2 * $ff)) / $UF, $iDec)
EndFunc

[/autoit] [autoit][/autoit] [autoit][/autoit]

**AspirinJunkie**

Interessante Lösung.
Wo stammt dieser Näherungsalgorithmus her und vor allem die Angabe, dass es immer im Bereich ± 0.05 km liegt?
Für dein Beispiel betrug die Differenz zur exakten Lösung auf dem angegeben Ellipsoid übrigens 6m was ja im Bereich liegt.

Vielleicht sollte man aber dazusagen: Bei Koordinatenberechnungen auf dem Ellipsoid müssen die angegebenen Koordinaten in einem Bezugssytem ermittelt wurden sein welchem das verwendete Ellipsoid zu Grunde liegt. Sonst sind die ermittelten Werte verfälscht. In Wikipedia z.B. sind die Angaben meist auf WGS84 bezogen (Referenzsystem für GPS-Koordinaten), also muss man auch dessen Werte für die Halbachse und Abplattung verwenden wenn man mit deren Koordinaten rechnen möchte.

Theoretisch kann man auf dem Ellipsoid bis zu jeder beliebigen Genauigkeit rechnen. Da es sich hier aber um elliptische Integrale handelt welche ein Randwertproblem darstellen, ist die Berechnung ziemlich tricky und nur iterativ zu lösen.
Hier mal eine Version welche auf etwa 1nm genau rechnet:

Strecke zweier Punkte auf einem Rotationsellipsoid mit hoher Genauigkeit

[autoit]

#Region Hilfsgrößen und Strukturdefinitionen
Global Const $PI = ACos(-1)
Global Const $TAG_GEOCOORD = "double B;double L" ; Koordinatenangabe in Dezimalwinkeln (B=Breite/Latitude, L=Länge/Longitude)
Global Const $TAG_GEOCOORD_HOUR = "int bD;int bM;double bS;int lD;int lM;double lS" ; Koordinateenangabe der Form: Winkel:Minute:Sekunde (b=Breite/Latitude, l=Länge/Longitude, D=Degrees, M=Minuten, S=Sekunden)
#EndRegion Hilfsgrößen und Strukturdefinitionen

[/autoit] [autoit][/autoit] [autoit]

#Region Eingabedaten
; Paris Observatory
$geoParis = 'E-2-20-14-N-48-50-11'
ConsoleWrite("Paris: " & @TAB & _Coord_PrintHour(_Coords_Degree2Hour(_Coords_GeoData2Degrees($geoParis))) & @CRLF)
; US Naval Observatory
$geoUSA = 'W-77-3-55.5-N-38-55-17'
ConsoleWrite("USA: " & @TAB & _Coord_PrintHour(_Coords_Degree2Hour(_Coords_GeoData2Degrees($geoUSA))) & @CRLF)
#EndRegion Eingabedaten

[/autoit] [autoit][/autoit] [autoit]

; Berechne die Distanz der beiden Punkte auf dem Ellipsoid
Global $a_Ret = _Ellipsoid_DistanceOf2Points(_Coords_GeoData2Degrees($geoParis), _Coords_GeoData2Degrees($geoUSA), 6378137, 1 / 298.257223563)
ConsoleWrite(StringFormat("Strecke: %15.9f m\n", $a_Ret[0]))
ConsoleWrite(StringFormat("Azimut im Punkt 1: %15.9f°\n", $a_Ret[1] * 180 / $PI))
ConsoleWrite(StringFormat("Azimut im Punkt 2: %15.9f°\n", $a_Ret[2] * 180 / $PI))

[/autoit] [autoit][/autoit] [autoit]

; #FUNCTION# ======================================================================================
; Name ..........: _Ellipsoid_DistanceOf2Points()
; Description ...: Berechnet die Entfernung zweier Punkte auf einem beliebigen Rotationsellipsoid
; Syntax ........: _Ellipsoid_DistanceOf2Points($P1, $P2[, $a = 6378137[, $f = 1 / 298.257222101]])
; Parameters ....: $P1 - Koordinaten von P1 als DezimalGrad-Angaben in der Struktur $TAG_GEOCOORD
; $P2 - Koordinaten von P2 als DezimalGrad-Angaben in der Struktur $TAG_GEOCOORD
; $a - [optional] Große Halbachse des Rotationsellipsoiden (weitere: http://de.wikipedia.org/wiki/Referenze…erenzellipsoide) (default:6378137)
; $f - [optional] Abplattungsfaktor des Rotationsellipsoiden (weitere: http://de.wikipedia.org/wiki/Referenze…erenzellipsoide) (default:1 / 298.257222101)
; $EPS- zu erreichende Genauigkeit
; Return values .: Success Array: [Strecke, Azimut1, Azimut2]
; Failure
; Author ........: AspirinJunkie
; Related .......: __Gn, __Mn, __Simp, __Sn, __Tn, CalcIntegral, CalcIntegralAdaptiv, f, IntMonteCarlo, Sign
; Remarks .......: Lösungsalgorithmus entsprechend H. Schmidt, RWTH Aachen 2000
; =================================================================================================
Func _Ellipsoid_DistanceOf2Points($P1, $P2, $a = 6378137, $f = 1 / 298.257223563, Const $EPS = 1e-15)
Local $s_Funktion ; Variable welche die zu integrierende Funktion f(x) als String enthält
Local $b = $a - $f * $a ; kleine Halbachse des Ellipsoides
Local $e2 = 1 - (($b / $a) ^ 2) ; Quadrat der 1. numerischen Exzentrizität

[/autoit] [autoit][/autoit] [autoit]

Local $B1 = $P1.B * $PI / 180 ; Breite von Punkt 1 in Radiant
Local $B2 = $P2.B * $PI / 180 ; Breite von Punkt 2 in Radiant
Local $L1 = $P1.L * $PI / 180 ; Länge von Punkt 1 in Radiant
Local $L2 = $P2.L * $PI / 180 ; Länge von Punkt 2 in Radiant

[/autoit] [autoit][/autoit] [autoit]

; reduzierte Breiten beta1 und beta2:
Local $beta1 = (Abs($B1 - ($PI / 2)) < 1E-5) ? $B1 : ATan(Sqrt(1 - $e2) * Tan($B1)) ; zusätzlicher Check ob eine Breite dem Pol entspricht
Local $beta2 = (Abs($B2 - ($PI / 2)) < 1E-5) ? $B2 : ATan(Sqrt(1 - $e2) * Tan($B2))

[/autoit] [autoit][/autoit] [autoit]

; Ausgangswerte für Restabweichung und Langendifferenz
Local $deltaL
Local $deltaLsoll = $L2 - $L1
Local $deltaLs = $L2 - $L1
Local $dL = $deltaLsoll
Local $dL1, $dL2 ; Längendifferenzen für beide Fälle

[/autoit] [autoit][/autoit] [autoit]

; Vorzeichen der Längendifferenz
Local $y = Sign($deltaLsoll)

[/autoit] [autoit][/autoit] [autoit]

; iterative Berechnung der Intergrationsgrenzen
Local $tanBetaMax ; tangens der maximalen Breite der geodätischen Linie
Local $h2, $h ; Clairaut-Konstante
Local $w1, $w2, $wu, $wo ; Integrationsgrenzen
Local $deltaL1, $deltaL2 ; Lösung für Fall 1 und Fall 2
Local $FALL ; Angabe welcher Fall vorherrscht
Local $d_ItMax = 50, $it = 0 ; Maximale Iterationsanzahl (um Endlosschleifen zu vermeiden)
Do
; tangens der maximalen Breite der geodätischen Linie
$tanBetaMax = Sqrt(1 - $e2) * (Sqrt(Tan($B1) * Tan($B1) + Tan($B2) * Tan($B2) - 2 * Tan($B1) * Tan($B2) * Cos($deltaLs)) / Sin($deltaLs))

[/autoit] [autoit][/autoit] [autoit]

; Clairaut-Konstante
$h2 = 1 / (1 + $tanBetaMax * $tanBetaMax)

[/autoit] [autoit][/autoit] [autoit]

; Integrationsgrenzen
$w1 = ASin(Tan($beta1) / $tanBetaMax)
$w2 = ASin(Tan($beta2) / $tanBetaMax)

[/autoit] [autoit][/autoit] [autoit]

If Abs($dL) < $EPS Or $it >= $d_ItMax Then ExitLoop ; Integrationsgrenzen hinreichend genau bestimmt

[/autoit] [autoit][/autoit] [autoit]

$wu = (Abs($w1) < Abs($w2) ? Abs($w1) : Abs($w2)) ; = Min(w1,w2)
$wo = (Abs($w1) > Abs($w2) ? Abs($w1) : Abs($w2)) ; = Max(w1,w2)

[/autoit] [autoit][/autoit] [autoit]

; Berechnung von wu und wo für Fall 1 (es gibt zwei mögliche Lösungen - daher beide berechnen und mit $deltaLsoll vergleichen)
$FALL = 1
$s_Funktion = "sqrt( 1-(" & $e2 * $h2 & " / (" & $h2 & " + " & (1 - $h2) & " * sin(x)*sin(x))))"

[/autoit] [autoit][/autoit] [autoit]

$deltaL1 = $y * Abs(_MathInt_Romberg($s_Funktion, $wu, $wo, $EPS)) ; Integralberechnung bis auf 1e-14 Genauigkeit
; Addition der Abschnitte symmetrisch zum Äquator
If ($w1 * $w2) < 0 Then $deltaL1 += $y * 2 * _MathInt_Romberg($s_Funktion, 0, $wu, $EPS)

[/autoit] [autoit][/autoit] [autoit]

; Längendifferenz $dL und Restabweichung DL für Fall I
$deltaL = $deltaL1
$dL1 = $deltaLsoll - $deltaL1
$dL = $dL1

[/autoit] [autoit][/autoit] [autoit]

; Addition der Abschnitte symmetrisch zu Bmax für Fall II
$deltaL2 = $deltaL1 + $y * 2 * _MathInt_Romberg($s_Funktion, $wo, 0.5 * $PI, $EPS)

[/autoit] [autoit][/autoit] [autoit]

; Längendifferenz $dL und Restabweichung DL für Fall II
$dL2 = $deltaLsoll - $deltaL2

[/autoit] [autoit][/autoit] [autoit]

; Entscheidung für Fall 1 oder 2
If (Abs($dL2) < Abs($dL1)) Then ; Fall 2 ist vorhanden - daher dessen Ergebnisse übernehmen
$FALL = 2
$deltaL = $deltaL2
$dL = $dL2
EndIf

[/autoit] [autoit][/autoit] [autoit]

; Korrektur von $deltaLs zur Berechnung von tanBmax
$deltaLs = $deltaLs + (1 + ($dL / ($FALL * $deltaL))) * $dL
$it += 1 ; Iterationszähler inkrementieren
Until 0
; Azimutberechnung
; Vorzeichen der Breitendifferenz
Local $x = Sign($w2 - $w1)
If $FALL = 2 Then $x = Sign(Sign($B1 + $B2) * $PI - $w2 - $w1)

[/autoit] [autoit][/autoit] [autoit]

; Azimute A1 und A2
Local $A1 = (1 - $x) * 0.5 * $PI + $x * $y * ASin(Sqrt($h2) / Cos($beta1))
Local $A2s = (1 - $x) * 0.5 * $PI + $x * $y * ASin(Sqrt($h2) / Cos($beta2))
If $FALL = 2 Then $A2s = $PI - $A2s
If $A1 < 0 Then $A1 += 2 * $PI
If $A2s < 0 Then $A2s += 2 * $PI
Local $A2 = $A2s + $PI
If $A2 > $PI Then $A2 -= 2 * $PI

[/autoit] [autoit][/autoit] [autoit]

; Berechnung der Streckenlänge
; Integrationsgrenzen v1 und v2 bestimmen
Local $v1 = ASin(Sin($beta1) / Sqrt(1 - $h2))
Local $v2 = ASin(Sin($beta2) / Sqrt(1 - $h2))

[/autoit] [autoit][/autoit] [autoit]

If (Abs($deltaLsoll - $PI) < 1E-10) Then $FALL = 2 ; check ob deltaLsoll = Pi
If $FALL = 2 Then $v2 = Sign($B1 + $B2) * $PI - $v2

[/autoit] [autoit][/autoit] [autoit]

; Berechnung der Streckenlänge durch numerische Intergration
$s_Funktion = "sqrt(1-" & $e2 & "*(1-sin(x)*sin(x)*" & 1 - $h2 & "))"
Local $S = $a * _MathInt_Romberg($s_Funktion, $v1, $v2, $EPS) ; Strecke zwischen beiden Punkten in m

[/autoit] [autoit][/autoit] [autoit]

Local $a_Ret = [$S, $A1, $A2]
Return $a_Ret
EndFunc ;==>_Ellipsoid_DistanceOf2Points

[/autoit] [autoit][/autoit] [autoit][/autoit] [autoit]

#Region Hilsfunktionen
; Signum (Vorzeichen) einer Zahl
Func Sign($x)
Return ($x < 0) ? -1 : ($x > 0) ? 1 : 0
EndFunc ;==>Sign

[/autoit] [autoit][/autoit] [autoit]

; #=============================================
; = _Coords_GeoData2Degrees =
; =============================================*
; Wandelt Koordinaten im Geodata-Format (z.B. E-2-20-14-N-48-50-11) - in normale Dezimalwinkel
Func _Coords_GeoData2Degrees($s_GeoDataForm)
; by AspirinJunkie
Local $aTemp = StringSplit($s_GeoDataForm, "-", 3)
Local $t_Ret = DllStructCreate($TAG_GEOCOORD)
$t_Ret.L = ($aTemp[1] * 3600 + $aTemp[2] * 60 + $aTemp[3]) / 3600 * (($aTemp[0] = "W") ? -1 : 1)
$t_Ret.B = ($aTemp[5] * 3600 + $aTemp[6] * 60 + $aTemp[7]) / 3600 * (($aTemp[4] = "S") ? -1 : 1)
Return $t_Ret
EndFunc ;==>_Coords_GeoData2Degrees

[/autoit] [autoit][/autoit] [autoit]

; #=============================================
; = _Coords_Degree2Hour =
; =============================================*
; Wandelt Koordinaten in Dezimalform in Winkel:Minute:Sekunde-Form (Bogenminuten)
Func _Coords_Degree2Hour($t_Degree)
; by AspirinJunkie
Local $t_Ret = DllStructCreate($TAG_GEOCOORD_HOUR)

[/autoit] [autoit][/autoit] [autoit]

$t_Ret.bD = Int($t_Degree.B)
Local $temp = $t_Degree.B - $t_Ret.bD
$t_Ret.bM = Int(60 * $temp)
$t_Ret.bS = Abs(3600 * $temp - 60 * $t_Ret.bM)
$t_Ret.bM = Abs($t_Ret.bM)

[/autoit] [autoit][/autoit] [autoit]

$t_Ret.lD = Int($t_Degree.L)
$temp = $t_Degree.L - $t_Ret.lD
$t_Ret.lM = Int(60 * $temp)
$t_Ret.lS = Abs(3600 * $temp - 60 * $t_Ret.lM)
$t_Ret.lM = Abs($t_Ret.lM)
Return $t_Ret
EndFunc ;==>_Coords_Degree2Hour

[/autoit] [autoit][/autoit] [autoit]

; Gibt die Koordinaten im Bogenminuten-Format als String zurück
Func _Coord_PrintHour($t_Hour)
; by AspirinJunkie
Return StringFormat("Lat: %2d° %2d' %2d"".%05d %1s\tLon: %2d° %2d' %2d"".%05d %1s", _
Abs(Int($t_Hour.bD)), $t_Hour.bM, Int($t_Hour.bS), Int(StringLeft(StringTrimLeft($t_Hour.bS, StringLen(Int($t_Hour.bS)) + 1), 5)), ($t_Hour.bD >= 0) ? "N" : "O", _
Abs(Int($t_Hour.lD)), $t_Hour.lM, Int($t_Hour.lS), Int(StringLeft(StringTrimLeft($t_Hour.lS, StringLen(Int($t_Hour.lS)) + 1), 5)), ($t_Hour.lD >= 0) ? "E" : "W")
EndFunc ;==>_Coord_PrintHour
#EndRegion Hilsfunktionen

[/autoit] [autoit][/autoit] [autoit][/autoit] [autoit][/autoit] [autoit][/autoit] [autoit]

#Region Funktionen zur Berechnung einer numerischen Integration
; #FUNCTION#=================================================================================
; Name...........: _MathInt_Romberg
; Description ...: Berechnet das bestimmte Integral eindimensionalen Funktion mit
; einer festzulegenden relativen Genauigkeit
; Syntax.........: _MathInt_Romberg(Const $a, Const $b, Const $EPS = 0.0001)
; Parameters ....: $cb_F - Funktionsvariable oder Funktionsstring einer 1D-Funktion f(x)
; $a - Untere Integrationsgrenze
; $b - obere Integrationsgrenze
; $EPS - zu erreichender relativer Integrationsfehler
; Hinweis: != absoluter Fehler, rel. Fehler = Abs(1-(berechneter Integralwert / wahrer Integralwert)
; $kMax - Maximale Anzahl an Halbierungsintervallen
; Return values .: Integral in den Grenzen $a-$b, @extended = Anzahl benötigter Iterationen
; @error = 1: Maximale Iterationsanzahl erreicht
; Remarks .......: Romberg-Integrationsalgorithmus
; Author ........: AspirinJunkie
;============================================================================================
Func _MathInt_Romberg(Const $cb_F, Const $a, Const $b, Const $EPS = 0.0001, Const $kMax = 30)
Local $a_R1[$kMax + 1], $a_R2[$kMax + 1]
Local $h = $b - $a, $n = 1, $sumf, $f

[/autoit] [autoit][/autoit] [autoit]

$a_R1[0] = 0.5 * $h * ((IsFunc($cb_F) ? $cb_F($a) : __F($cb_F, StringFormat("%.20f", $a))) + (IsFunc($cb_F) ? $cb_F($b) : __F($cb_F, StringFormat("%.20f", $b)))) ; Erste Approximation

[/autoit] [autoit][/autoit] [autoit]

For $k = 1 To $kMax ; Schrittweise halbieren
$sumf = 0
For $i = 1 To $n
$sumf += (IsFunc($cb_F) ? $cb_F($a + ($i - 0.5) * $h) : __F($cb_F, StringFormat("%.20f", $a + ($i - 0.5) * $h)))
Next
$a_R2[0] = 0.5 * ($a_R1[0] + $h * $sumf) ; Trapezoid-Formel
$f = 1e0
For $j = 1 To $k ; Quadratur-Ordnung erhöhen
$f *= 4
$a_R2[$j] = ($f * $a_R2[$j - 1] - $a_R1[$j - 1]) / ($f - 1) ; Neue Approximation
Next
If $k > 1 Then ; Check auf Konvergenz
If (Abs($a_R2[$k] - $a_R1[$k - 1]) <= ($EPS * Abs($a_R2[$k]))) Then ExitLoop
If (Abs($a_R2[$k]) <= $EPS) And (Abs($a_R2[$k]) <= Abs($a_R2[$k] - $a_R1[$k - 1])) Then ExitLoop
EndIf
$h *= 0.5 ; Halbiere Schrittweite
$n *= 2
For $j = 0 To $k
$a_R1[$j] = $a_R2[$j]
Next
Next
If $k >= $kMax Then Return SetError(1, $k, $a_R2[$k - 1])
Return SetExtended($k, $a_R2[$k])
EndFunc ;==>_MathInt_Romberg

[/autoit] [autoit][/autoit] [autoit]

; Wrapperfunktion für den Fall, dass die zu integrierende Funktion als String übergeben wird
Func __F($s_Func, $x)
$s_Func = StringReplace($s_Func, "exp", "§")
$s_Func = StringReplace($s_Func, "x", $x)
$s_Func = StringReplace($s_Func, "§", "exp")
Return Execute($s_Func)
EndFunc ;==>__F
#EndRegion Funktionen zur Berechnung einer numerischen Integration

[/autoit]

Ein Hinweis vielleicht noch: Bei deinem Koordinatenparsing negierst du die Koordinatenangaben wenn sie Ostkoordinaten sind. Normal ist es aber eigentlich Westkoordinaten zu negieren.
Am Ergebnis ändert das aber nichts - eher ein Schönheitsdetail.

Ein sehr schönes Online-Tool mit dem man sich das ganze berechnen und vor allem anzeigen lassen kann ist hier zu finden: http://geographiclib.sourceforge.net/scripts/geod-google.html (unser Anwendungsfall ist der "Inverse")

Edit: Hab auch mal die Azimutwinkel im Anfangs und Endpunkt mit hinzugefügt, damit man weiß in welche Richtung man loslaufen muss und in welche man ankommt.
Edit2: Auf Romberg-Integration umgestellt.
Edit3: Link zum Online-Tool hinzugefügt

**minx**

Zitat von AspirinJunkie

Wo stammt dieser Näherungsalgorithmus her

Laut meinem Gekritzel ist das eine Methode von Andoyer. Das stammt aus meinen Aufzeichnungen ausm Astronomie-Kurs. Wir hatten das damals in Excel umgesetzt, aber ich wollte es nochmal als AutoIt Funktion umschreiben.

Zitat von AspirinJunkie

Ein Hinweis vielleicht noch: Bei deinem Koordinatenparsing negierst du die Koordinatenangaben wenn sie Ostkoordinaten sind. Normal ist es aber eigentlich Westkoordinaten zu negieren.
Am Ergebnis ändert das aber nichts - eher ein Schönheitsdetail.

Da habe ich genau das Gegenteil gelernt und damit auch in Prüfungen gerechnet

Zitat von AspirinJunkie

Theoretisch kann man auf dem Ellipsoid bis zu jeder beliebigen Genauigkeit rechnen.

Das ist sicher richtig, aber diese Näherung reicht völlig aus, wenn man sie z.B. mit der geläufigen AGPS Triangulation vergleicht.

**AspirinJunkie**

Zitat von minx

Da habe ich genau das Gegenteil gelernt und damit auch in Prüfungen gerechnet

Na dann: Glück gehabt. Wenn du nicht gerade mal derartige Prüfungen schreibst, kannst du es ja von nun an richtig herum machen.

Zitat von minx

Das ist sicher richtig, aber diese Näherung reicht völlig aus, wenn man sie z.B. mit der geläufigen AGPS Triangulation vergleicht.

Das wollte ich Dir hiermit ja zum Ausdruck bringen. Wenn deine Koordinaten aus GPS-Messungen stammen müsstest du auch mit dem WGS84-Ellipsoid rechnen anstatt einem beliebigen.
Für das Beispiel hier hat man schon allein dadurch einen Fehler von 3m der nicht sein müsste.
Herkömmliche GPS-Messungen sind außerdem deutlich genauer als 50m und eine höhere Rechengenauigkeit für die Distanz daher auch hierbei durchaus wünschenswert.

Trubadour

Hallo,
wenn ich das Beispiel (von dir AspirinJunkie) ausprobiere, erhalte ich folgende zwei Richtungen
Azimut im Punkt 1: 291.833555°
Azimut im Punkt 2: 51.793617°

Zusammen müßten sie doch eigentlich 360° ergeben oder denke ich völlig falsch?
Oder ist der WGS84 Geoid so kartoffelig?

**AspirinJunkie**

Ja das ist erstmal bisschen verwirrend, aber richtig so.
Kurz: Die Winkel müssen nicht zusammen 360° ergeben.
Wenn man sich auf einer Kugel ein Dreieck zeichnet und die Kanten dabei die jeweils kürzesten Verbindungen auf der Kugel sind dann ist die Innenwinkelsumme in diesem Dreieck größer als 180°.
Die Gesetzmäßigkeiten der planaren Geometrie gelten nicht mehr vollständig.

Zitat von Trubadour

Oder ist der WGS84 Geoid so kartoffelig?

Nein nein es ist ein ganz glattes normales Rotationsellipsoid.
Kartoffelig ist erst das Geoid.

Trubadour

Vielen Dank, besonders für das Geoid.