Geodaten - Richtungswinkel und Entfernung berechnen

qwert23

Hallo,

ich bitte mal um Unterstützung, da meine Kenntnisse in der Richtungswinkelberechnung nicht ausreichen.

Mein Ziel ist es, aus GPX-Tracks die Entfernung und die Richtung bis zum nächsten Trackpunkt zu berechnen.

Die Trackpunkte liegen nicht weit auseinander. Ich gehe von maximal 1 km aus, eher im Meterbereich. Also ohne Erdkrümmung-berücksichtigung etc..

Im Anhang habe ich eine Testdatei mit extrahierten GPX Koordinaten beigefügt. In diesem Zahlenformat werden GPX Tracks aufgezeichnet. Diese Werte stammen aus einer Aufzeichnung vom Garmin Oregon.

Ich möchte ermitteln, von der Zeile 1 zur Zeile 2 sind es wie viele Meter und in welcher Richtung liegt der nächste Punkt.

Dann von Zeile 2 zur Zeile 3 sind es wie viele Meter und in welcher Richtung liegt der nächste Punkt und so weiter....

Bei der Richtungsangabe, so haben ich gelesen, wird 0 Grad für Norden, 90 Grad für Osten, 180 Grad für Süden und 270 Grad für Westen verwendet. Das würde ich als Ausgabe so übernehmen mit allen Zwischenwerten.

So, hier mein kleiner Versuch, doch auf Grundlage der beiden verschiedenen "kopierten" Ansätze kommen völlig verschieden Werte heraus, die ich nicht nachvollziehen kann. Jetzt ist ein Sehender gefordert

Was ist falsch an der Berechnung? Oder besser, wie lautet die korrekte Berechnung?

Vielen vielen Dank schon mal im Voraus.

**AspirinJunkie**

Das Thema ist deutlich verzwickter als du eventuell momentan vermutest.

Problem ist, dass du keine ebenen Koordinaten hast (dafür würde die atan2-Formel passen) sondern Koordinaten auf einer Kugel bzw. Ellipsoid.

Und darüber hinaus auch keine kartesischen Koordinaten sondern Breiten und Längenangaben.

Die Formeln die du hier hast, kannst du eigentlich gleich wegschmeißen.

Folgende Möglichkeiten hast du:

Überführung der Kugel/Ellipsoidkoordinaten in ein ebenes Koordinatensystem durch Abbildung. Z.B. Gauß-Krüger Abbildung oder UTM-Abbildung.
- Sobald du dann ebene Koordinaten vorliegen hast, kannst du wieder ganz einfach mit atan2 und dem Pythagoras rechnen.
- Nicht wirklich trivial und easy umsetzbar. Am einfachsten wäre noch QGIS zu nehmen, die ellipsoidischen Koordinaten reinzuhauen und in ein entsprechendes ebenes Koordinatensystem zu konvertieren. Dann hast du ebene Koordinaten und kannst diese wie in deinem Ansatz mit AutoIt behandeln.
Berechnung auf einer Kugel (einfacher als Ellipsoid). Das Fachgebiet nennt sich "Sphärische Trigonometrie" - viel Spaß beim einstudieren. Da gibt es so lustiges Zeug wie "Zweiecke"...
DIrekte Berechnung der sogenannten "2. Geodätischen Hauptaufgabe" auf dem Ellipsoid:
- Eigentlich auch alles andere als trivial aber witzigerweise hatte ich sowas vor zig Jahren schonmal in AutoIt umgesetz: [Post]
- Angepasst auf deinen Fall dann ungefähr so:

AutoIt

#include <File.au3>

#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)
#EndRegion Hilfsgrößen und Strukturdefinitionen

#Region Eingabedaten
Global $a_Coords
#_FileReadToArray("testkoord.txt", $a_Coords, $FRTA_INTARRAYS, ";")
_FileReadToArray("testkoord_bn4.txt", $a_Coords, $FRTA_INTARRAYS, ";")

; Koordinaten in die hier verwendete Strukturform übertragen:
Global $a_Coords_Struct[UBound($a_Coords)]
For $i = 0 To UBound($a_Coords) - 1
    $a_Coords_Struct[$i] = _Coords2Struct( Number(($a_Coords[$i])[0], 3), Number(($a_Coords[$i])[1], 3) )
Next
#EndRegion Eingabedaten

#Region Berechnung
; Distanzen und Richtungen zwischen aufeinander folgenden Punkten berechnen:
Global $a_Ergebnisse
ConsoleWrite(" Von   Nach   Strecke [m]   Azimut [°]" & @CRLF & "-------------------------------------" & @CRLF)
For $i = 0 To UBound($a_Coords_Struct) - 2
    $a_Ergebnisse = _Ellipsoid_DistanceOf2Points($a_Coords_Struct[$i], $a_Coords_Struct[$i+1])
    ConsoleWrite(StringFormat("%3d    %3d     %7.1f      %6.2f \n", $i, $i+1, $a_Ergebnisse[0], $a_Ergebnisse[1] * 180 / $PI))
Next
#EndRegion Berechnung


; #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/R…chtige_Referenzellipsoide) (default:6378137)
; $f - [optional] Abplattungsfaktor des Rotationsellipsoiden (weitere: http://de.wikipedia.org/wiki/R…chtige_Referenzellipsoide) (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-5)
    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

    If $P1.B = $P2.B Then $P1.B += 0.000001

    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

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

; 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

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

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

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

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

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

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

        ; 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))))"

        $deltaL1 = $y * Abs(_Math_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 * _Math_Romberg($s_Funktion, 0, $wu, $EPS)

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

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

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

        ; 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

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

    ; 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

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

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

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

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


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

; Bastelt ellipsoidische Koordinaten in die hier verwendeten Strukturen
Func _Coords2Struct($B, $L)
    Local $t_Ret = DllStructCreate($TAG_GEOCOORD)
    $t_Ret.L = $L
    $t_Ret.B = $B
    Return $t_Ret
EndFunc
#EndRegion Hilsfunktionen


#Region Funktionen zur Berechnung einer numerischen Integration
; #FUNCTION#=================================================================================
; Name...........: _Math_Romberg
; Description ...: Berechnet das bestimmte Integral eindimensionalen Funktion mit
; einer festzulegenden relativen Genauigkeit
; Syntax.........: _Math_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 _Math_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

    $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

    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   ;==>_Math_Romberg

; 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

Alles anzeigen

Ist aber zig Jahre her und ich bin auch schon bisschen raus.
Vielleicht ist es aber ein Ansatz für dich.

qwert23

Hallo AspirinJunkie,

oh je, das habe ich befürchtet.

Deshalb vielen vielen Dank für Deine Erläuterungen und dem Script.

Herzlicher Gruß

Bernd