VBA Tipp: Distanz Luftlinie berechnen

Aus DBWiki
Wechseln zu: Navigation, Suche

Anforderung 1

Berechnung der Entfernung zwischen zwei Orten anhand der Angabe der Längen- und Breitengrade

Lösung 1

'Quelle: http://www.dbwiki.net/
 
 
Option Explicit
 
Public Type Koordinaten
   Laenge As Double
   Breite As Double
End Type
 
 
Private Const EARTH   As Double = 6378137#              'Äquatorialradius in Metern
Private Const FLAT    As Double = 6.69438000426084E-03  'Abflachung
Private Const DPI     As Double = 3.14159265358979      'Kreiszahl PI
Private Const DEG2RAD As Double = DPI / 180#            'Umrechnungsfaktor von Grad in Bogenmaß
 
 
'
' Original kommt von: http://www.dbdev.org/members/schwanke/koordinaten.zip
' © Christa Schwanke, Juni 2000
' http://www.dbdev.org/downloads.html
'
' Das Ergebnis ist in KM
' Längen und Breitengrade können hier ermittelt werden (Stand Januar 2014):
' http://www.linkr.de/2009/03/24/mit-google-maps-laengen-breitengrade-ermitteln-freies-geocoding-skript/
'
 
Private Const EndOfDeclaration As String = "EndOfDeclaration"
 
'
'_______________________________________
'
Private Sub PseudoMain()
 
   Dim Hamburg As Koordinaten
   Dim Berlin  As Koordinaten
 
   Hamburg.Laenge = 9.9936817999: Hamburg.Breite = 53.5510846
   Berlin.Laenge = 13.404954:     Berlin.Breite = 52.5200066
 
   Debug.Print CalcLuftlinie(Hamburg, Berlin)
 
End Sub
 
'
' CalcLuftlinie
'_______________________________________
'
Public Function CalcLuftlinie(Ort1 As Koordinaten, Ort2 As Koordinaten) As Double
 
   On Error GoTo error_lcl
   '------
 
   CalcLuftlinie = Sqr(Abs((fct_dblz(Ort1.Breite) - fct_dblz(Ort2.Breite)) _
                         * (fct_dblz(Ort1.Breite) - fct_dblz(Ort1.Breite)) _
                         + (fct_dblx(Ort1) - fct_dblx(Ort2)) * (fct_dblx(Ort1) - fct_dblx(Ort2)) _
                         + (fct_dbly(Ort1) - fct_dbly(Ort2)) * (fct_dbly(Ort1) - fct_dbly(Ort2)))) / 1000
 
   '------
error_lcl:
 
   Select Case Err
      Case 0
      Case Else
         CalcLuftlinie = -1
   End Select
 
End Function
 
'
'_______________________________________
'
Private Function fct_dblbeta(Breite As Double) As Double
 
   On Error GoTo error_lcl
   '------
 
   fct_dblbeta = EARTH / Sqr(Abs(1 - FLAT * Sin(Breite * DEG2RAD) * Sin(Breite * DEG2RAD)))
 
   '------
error_lcl:
 
   Select Case Err
      Case 0
      Case Else
         fct_dblbeta = -1
   End Select
 
End Function
 
'
'_______________________________________
'
Private Function fct_dblx(Ort As Koordinaten) As Double
 
   On Error GoTo error_lcl
   '------
 
   fct_dblx = fct_dblbeta(Ort.Breite) * Cos(Ort.Breite * DEG2RAD) * Cos(Ort.Laenge * DEG2RAD)
 
   '------
error_lcl:
 
   Select Case Err
      Case 0
      Case Else
         fct_dblx = -1
   End Select
 
End Function
 
'
'_______________________________________
'
Private Function fct_dbly(Ort As Koordinaten) As Double
 
   On Error GoTo error_lcl
   '------
 
   fct_dbly = fct_dblbeta(Ort.Breite) * Cos(Ort.Breite * DEG2RAD) * Sin(Ort.Laenge * DEG2RAD)
 
   '------
error_lcl:
 
   Select Case Err
      Case 0
      Case Else
         fct_dbly = -1
   End Select
 
End Function
 
'
'_______________________________________
'
Private Function fct_dblz(Breite As Double) As Double
 
   On Error GoTo error_lcl
   '------
 
   fct_dblz = fct_dblbeta(Breite) * (1 - FLAT) * Sin(Breite * DEG2RAD)
 
   '------
error_lcl:
 
   Select Case Err
   Case 0:
   Case Else: fct_dblz = -1
   End Select
 
End Function

Anforderung 2

Mir sind die Ergebnisse aus Lösung 1 zu ungenau, weil sie bei größeren Entfernungen Abweichungen von mehreren Hundert Kilometern zu beobachten sind. Ich möchte deshalb Genauigkeiten im Bereich von 50 Metern bei der Berechnung anstreben.

Lösung 2

In der Wikipedia findet sich dazu im Artikel Orthdrome eine Lösung, die nach Jean Meeus im Buch Astronomical Algorithms vom Astronomen Henri Andoyer stammen soll, während die Wikipedia die Herkunft des Algorithmus Thaddeus Vincenty zuschreibt.

Der folgende Code basiert auf dem Wikipedia-Artikel, wobei man bei der Berechnung einen von vier Referenzellipsoid für die Erde vorgeben kann: WGS 72, EARTH 76 (von der IAU lange Zeit verwendet), GRS 80, WGS 84, . Voreingestellt ist WGS84, was z.Z. auch von den gebräuchlichen Navigationssystemen genutzt wird. Interessant ist, dass die Entfernungsberechnung über alle vier Modelle untereinander weniger als 10 Meter Abweichungen aufzeigen.

Die Zusätzlichen Methoden FromSexa bzw. FromSexaSec werden nur vom Beispiel-Code genutzt, um bequem Gradeingaben vornehmen zu können.

Folgender Code ist dazu in einem globalen Modul zu verwenden.

'Quelle: http://www.dbwiki.net/
 
Option Explicit
 
 
Public Enum EllipsoidType
   WGS72
   EARTH76  'IAU 1976 Ellipsoid
   GRS80
   WGS84
End Enum
#If 0 Then
Dim WGS72, EARTH76, GRS80, WGS84
#End If
 
'Coord gibt die geographischen Koordinaten auf der Erde an.
'
'Der Längengrad wird positiv vom Greenwich-Meridian nach Westen gemessen.
Public Type Coord
   Lat As Double  'Breitengrad (phi)
   Lon As Double  'Längengrad (psi oder L)
End Type
 
'Wird von Methode FromSexa verwendet
Public Const DMINUS   As Integer = 45              'ASC-Wert von "-"
 
 
'Distance ist der Abstand zwischen zwei Punkten entlang der Oberfläche eines
'Ellipsoids.
'
'In c1 und in c2 sind die Koordinaten in Dezimalgrad der beiden Orte auf der
'Erde anzugeben.
'
'Ergebniseinheit ist die Einheit von Er (typischerweise Km).
'
'Quelle: https://de.wikipedia.org/wiki/Orthodrome#Genauere_Formel_zur_Abstandsberechnung_auf_der_Erde
Public Function Distance(ByRef c1 As Coord, ByRef c2 As Coord, _
                         Optional ByVal et As EllipsoidType = WGS84) As Double
 
   Const PI  As Double = 3.14159265358979 'Kreiszahl Pi
   Const D2R As Double = PI / 180#        'Umrechnungsfaktor Grad -> Bogenmaß
 
   '===========================================================================
 
   Dim Er As Double  'Erdradius in Km
   Dim Fl As Double  'Erdabflachung
 
   Select Case et
   Case WGS72:   Er = 6378.135: Fl = 1# / 298.26
   Case EARTH76: Er = 6378.14:  Fl = 1# / 298.257
   Case GRS80:   Er = 6378.137: Fl = 1# / 298.257222101
   Case WGS84:   Er = 6378.135: Fl = 1# / 298.257223563
   End Select
 
   '===========================================================================
 
   Dim f          As Double
   Dim sin2F      As Double
   Dim cos2F      As Double
 
   Dim g          As Double
   Dim sin2G      As Double
   Dim cos2G      As Double
 
   Dim lambda     As Double
   Dim sin2Lambda As Double
   Dim cos2Lambda As Double
 
   Dim s          As Double
   Dim c          As Double
   Dim omega      As Double
   Dim d          As Double
   Dim t          As Double
   Dim h1         As Double
   Dim h2         As Double
 
 
   f = (c1.Lat + c2.Lat) / 2#
   g = (c1.Lat - c2.Lat) / 2#
   lambda = (c1.Lon - c2.Lon) / 2#
 
   sincos2 f * D2R, sin2F, cos2F
   sincos2 g * D2R, sin2G, cos2G
   sincos2 lambda * D2R, sin2Lambda, cos2Lambda
 
   s = sin2G * cos2Lambda + cos2F * sin2Lambda
   c = cos2G * cos2Lambda + sin2F * sin2Lambda
 
   omega = Atn(Sqr(s / c))
 
   d = 2# * omega * Er
   t = Sqr(s * c) / omega
 
   h1 = (3 * t - 1) / (2# * c)
   h2 = (3 * t + 1) / (2# * s)
 
   Distance = d * (1# + Fl * (h1 * sin2F * cos2G - h2 * cos2F * sin2G))
 
End Function
 
'Hilfsmethode
Private Sub sincos2(ByVal x As Double, ByRef s2 As Double, ByRef c2 As Double)
   Dim s As Double
   Dim c As Double
 
   s = Sin(x): c = Cos(x)
   s2 = s * s: c2 = c * c
End Sub

Beispiel

Die Ergebnisse von Aufruf DistBerlinTokyo()<tt> werden im VBA-Direktfenster ausgegeben.

'==============================================================================
 
'FromSexa konvertiert von geparsten sexagesimalen Winkelkomponenten in einen
'einzelnen Double-Wert.
'
'Das Ergebnis ist die Einheit d, die erste oder "größte" sexagesimale
'Komponente.
'
'Üblicherweise übergibt man nicht-negative Werte für d, m und s; um einen
'negativen Wert anzugeben, übergibt man "-" für neg.  Jeder andere Wert,
'wie z.B. " ", "+" oder einfach Chr$(0), lässt das Ergebnis nicht negativ
'werden.
'
'Es gibt jedoch keine Grenzen für d, m oder s.  Negative Werte oder
'Werte > 60 für m und s sind z.B. erlaubt.  Die Segmentwerte werden kombi-
'niert und wenn neg "-" ist, wird diese Summe negiert.
'
'Diese Funktion wird allgemein als DMSToDegrees o.ä. bezeichnet, aber die
'Interpretation von d als Grad ist beliebig.  Die Funktion funktioniert
'auch bei Stunden, Minuten und Sekunden.  Unabhängig von den Einheiten
'von d, ist  m ein sexagesimaler Teil von d und s ist ein sexagesimaler
'Teil von m.
Public Function FromSexa(ByVal neg As Byte, ByVal d As Long, _
                         ByVal m As Long, ByVal s As Double) As Double
   FromSexa = FromSexaSec(neg, d, m, s) / 3600
End Function
 
 
'FromSexaSec konvertiert von geparsten sexagesimalen Winkelkomponenten in
'einen einzelnen Double-Wert.
'
'Das Ergebnis ist die Einheit s, die letzte oder "kleinste" sexagesimale
'Komponente.
'
'Ansonsten arbeitet FromSexaSec genauso wie FromSexa.
Public Function FromSexaSec(ByVal neg As Byte, ByVal d As Long, _
                            ByVal m As Long, ByVal s As Double) As Double
   s = ((d * 60 + m) * 60) + s
   If neg = DMINUS Then
      FromSexaSec = -s
   Else
      FromSexaSec = s
   End If
End Function
 
'============================== Beispiel =====================================
 
Sub DistBerlinTokyo()
 
   Dim c1 As Coord  'Berlin
   Dim c2 As Coord  'Tokyo
 
   c1.Lat = FromSexa(0, 52, 31, 12.0288)
   c1.Lon = FromSexa(0, 13, 24, 17.8344)
 
   c2.Lat = FromSexa(0, 35, 38, 10.1952)
   c2.Lon = FromSexa(0, 139, 50, 22.1208)
 
   Debug.Print "WGS72:", Format$(Distance(c1, c2, WGS72), "0.000 km")
   Debug.Print "EARTH76:", Format$(Distance(c1, c2, EARTH76), "0.000 km")
   Debug.Print "GRS80:", Format$(Distance(c1, c2, GRS80), "0.000 km")
   Debug.Print "WGS84:", Format$(Distance(c1, c2, WGS84), "0.000 km")
 
   Dim k1 As Koordinaten
   Dim k2 As Koordinaten
 
   k1.Breite = FromSexa(0, 52, 31, 12.0288)
   k1.Laenge = FromSexa(0, 13, 24, 17.8344)
 
   k2.Breite = FromSexa(0, 35, 38, 10.1952)
   k2.Laenge = FromSexa(0, 139, 50, 22.1208)
 
   'im Vergleich dazu die erste Lösung
   Debug.Print
   Debug.Print "Lös. 1:", Format$(CalcLuftlinie(k1, k2), "0.000 km")
 
End Sub

Wikilinks