Die Distanz zwischen zwei Breitengrad-Längengrad-Punkten berechnen? (Haversine-Formel)

Question

Wie berechne ich die Entfernung zwischen zwei Punkten, die durch Breiten- und Längengrad angegeben sind?

Zur Klarstellung möchte ich die Entfernung in Kilometern; die Punkte verwenden das WGS84-System und ich möchte die relativen Genauigkeiten der verfügbaren Methoden verstehen.

Answer 1

Wenn Sie die Fahrstrecke/Fahrtroute möchten (hier gepostet, weil dies das erste Ergebnis für die Entfernung zwischen zwei Punkten bei Google ist, aber für die meisten Menschen ist die Fahrstrecke nützlicher), können Sie den Google Maps Distance Matrix Service verwenden:

getDrivingDistanceBetweenTwoLatLong(ursprung, ziel) {

 return new Observable(subscriber => {
  let service = new google.maps.DistanceMatrixService();
  service.getDistanceMatrix(
    {
      origins: [new google.maps.LatLng(ursprung.lat, ursprung.long)],
      destinations: [new google.maps.LatLng(ziel.lat, ziel.long)],
      travelMode: 'DRIVING'
    }, (antwort, status) => {
      if (status !== google.maps.DistanceMatrixStatus.OK) {
        console.log('Fehler:', status);
        subscriber.error({error: status, status: status});
      } else {
        console.log(antwort);
        try {
          let wertInMetern = antwort.rows[0].elements[0].distance.value;
          let wertInKm = wertInMetern / 1000;
          subscriber.next(wertInKm);
          subscriber.complete();
        }
       catch(fehler) {
        subscriber.error({error: fehler, status: status});
       }
      }
    });
});
}

Answer 2

Ich habe die Berechnung vereinfacht, indem ich die Formel vereinfacht habe.

Hier ist es in Ruby:

include Math
earth_radius_mi = 3959
radians = lambda { |deg| deg * PI / 180 }
coord_radians = lambda { |c| { :lat => radians[c[:lat]], :lng => radians[c[:lng]] } }

# from/to = { :lat => (latitude_in_degrees), :lng => (longitude_in_degrees) }
def haversine_distance(from, to)
  from, to = coord_radians[from], coord_radians[to]
  cosines_product = cos(to[:lat]) * cos(from[:lat]) * cos(from[:lng] - to[:lng])
  sines_product = sin(to[:lat]) * sin(from[:lat])
  return earth_radius_mi * acos(cosines_product + sines_product)
end

Answer 3

Hier ist die Implementierung in VB.NET, diese Implementierung gibt Ihnen das Ergebnis in KM oder Meilen basierend auf einem übergebenen Enum-Wert.

Public Enum DistanceType
    Meilen
    Kilometer
End Enum

Public Structure Position
    Public Latitude As Double
    Public Longitude As Double
End Structure

Public Class Haversine

    Public Function Distance(Pos1 As Position,
                             Pos2 As Position,
                             DistType As DistanceType) As Double

        Dim R As Double = If((DistType = DistanceType.Meilen), 3960, 6371)

        Dim dLat As Double = Me.toRadian(Pos2.Latitude - Pos1.Latitude)

        Dim dLon As Double = Me.toRadian(Pos2.Longitude - Pos1.Longitude)

        Dim a As Double = Math.Sin(dLat / 2) * Math.Sin(dLat / 2) + Math.Cos(Me.toRadian(Pos1.Latitude)) * Math.Cos(Me.toRadian(Pos2.Latitude)) * Math.Sin(dLon / 2) * Math.Sin(dLon / 2)

        Dim c As Double = 2 * Math.Asin(Math.Min(1, Math.Sqrt(a)))

        Dim result As Double = R * c

        Return result

    End Function

    Private Function toRadian(val As Double) As Double

        Return (Math.PI / 180) * val

    End Function

End Class

Answer 4

Eine der Hauptprobleme bei der Berechnung von Entfernungen - insbesondere großer Entfernungen - ist die Berücksichtigung der Krümmung der Erde. Wenn nur die Erde flach wäre, wäre die Berechnung der Entfernung zwischen zwei Punkten so einfach wie die einer geraden Linie! Die Haversine-Formel enthält eine Konstante (es ist die Variable R unten), die den Radius der Erde darstellt. Je nachdem, ob in Meilen oder Kilometern gemessen wird, würde es jeweils 3956 mi oder 6367 km betragen.

Die Grundformel lautet:

dlon = lon2 - lon1
dlat = lat2 - lat1
a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
c = 2 * atan2( sqrt(a), sqrt(1-a) )
distance = R * c (wo R der Radius der Erde ist)

R = 6367 km ODER 3956 mi

     lat1, lon1: Die Breite und Länge von Punkt 1 (in Dezimalgraden)
     lat2, lon2: Die Breite und Länge von Punkt 2 (in Dezimalgraden)
     einheit: Die Maßeinheit, in der die Ergebnisse berechnet werden sollen, wobei:
     'M' Meilen sind (Standard)
     'K' Kilometer sind
     'N' nautische Meilen sind

Beispiel

function distance(lat1, lon1, lat2, lon2, unit) {
    try {
        var radlat1 = Math.PI * lat1 / 180
        var radlat2 = Math.PI * lat2 / 180
        var theta = lon1 - lon2
        var radtheta = Math.PI * theta / 180
        var dist = Math.sin(radlat1) * Math.sin(radlat2) + Math.cos(radlat1) * Math.cos(radlat2) * Math.cos(radtheta);
        dist = Math.acos(dist)
        dist = dist * 180 / Math.PI
        dist = dist * 60 * 1.1515
        if (unit == "K") {
            dist = dist * 1.609344
        }
        if (unit == "N") {
            dist = dist * 0.8684
        }
        return dist
    } catch (err) {
        console.log(err);
    }
}

Answer 5

Da dies die beliebteste Diskussion zum Thema ist, werde ich hier meine Erfahrung aus Ende 2019-Anfang 2020 hinzufügen. Um den vorhandenen Antworten hinzuzufügen - mein Fokus lag darauf, eine genaue UND schnelle (d.h. vektorisierte) Lösung zu finden.

Beginnen wir mit dem, was hier hauptsächlich von Antworten verwendet wird - der Haversine-Ansatz. Es ist trivial zu vektorisieren, siehe Beispiel in Python unten:

def haversine(lat1, lon1, lat2, lon2):
    """
    Berechnet den Großkreisabstand zwischen zwei Punkten auf der Erde
    (angegeben in Dezimalgrad)

    Alle Argumente müssen die gleiche Länge haben.
    Die Entfernungen sind in Metern angegeben.

    Ref:
    https://stackoverflow.com/questions/29545704/fast-haversine-approximation-python-pandas
    https://ipython.readthedocs.io/en/stable/interactive/magics.html
    """
    Radius = 6.371e6
    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])

    dlon = lon2 - lon1
    dlat = lat2 - lat1

    a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2

    c = 2 * np.arcsin(np.sqrt(a))
    s12 = Radius * c

    # Anfangs-Azimut in Grad
    y = np.sin(lon2-lon1) * np.cos(lat2)
    x = np.cos(lat1)*np.sin(lat2) - np.sin(lat1)*np.cos(lat2)*np.cos(dlon)
    azi1 = np.arctan2(y, x)*180./math.pi

    return {'s12':s12, 'azi1': azi1}

In Bezug auf die Genauigkeit ist es am wenigsten genau. Wikipedia gibt eine relative Abweichung von durchschnittlich 0,5% ohne Quellen an. Meine Experimente zeigen weniger Abweichung. Hier ist ein Vergleich anhand von 100.000 zufälligen Punkten mit meiner Bibliothek, die auf Millimeterebene genau sein sollte:

np.random.seed(42)
lats1 = np.random.uniform(-90,90,100000)
lons1 = np.random.uniform(-180,180,100000)
lats2 = np.random.uniform(-90,90,100000)
lons2 = np.random.uniform(-180,180,100000)
r1 = inverse(lats1, lons1, lats2, lons2)
r2 = haversine(lats1, lons1, lats2, lons2)
print("Maximaler absoluter Fehler: {:4.2f}m".format(np.max(r1['s12']-r2['s12'])))
print("Durchschnittlicher absoluter Fehler: {:4.2f}m".format(np.mean(r1['s12']-r2['s12'])))
print("Maximaler relativer Fehler: {:4.2f}%".format(np.max((r2['s12']/r1['s12']-1)*100)))
print("Durchschnittlicher relativer Fehler: {:4.2f}%".format(np.mean((r2['s12']/r1['s12']-1)*100)))

Ausgabe:

Maximaler absoluter Fehler: 26671,47m
Durchschnittlicher absoluter Fehler: -2499,84m
Maximaler relativer Fehler: 0,55%
Durchschnittlicher relativer Fehler: -0,02%

Im Durchschnitt also eine Abweichung von 2,5km bei 100.000 zufälligen Koordinatenpaaren, was für die meisten Fälle akzeptabel sein könnte.

Die nächste Option ist die Vincenty-Formeln, die genau bis zu Millimetern sind, abhängig von Konvergenzkriterien und auch vektorisierbar sind. Es gibt jedoch ein Problem mit der Konvergenz in der Nähe von antipodalen Punkten. Sie können die Konvergenz an diesen Punkten erhöhen, indem Sie die Konvergenzkriterien lockern, aber die Genauigkeit sinkt auf 0,25% und mehr. Außerhalb antipodaler Punkte wird Vincenty Ergebnisse liefern, die im Durchschnitt eine relative Fehlerabweichung von weniger als 1.e-6 aufweisen.

Geographiclib, hier erwähnt, ist wirklich der aktuelle Goldstandard. Es hat mehrere Implementierungen und ist ziemlich schnell, insbesondere wenn Sie die C++-Version verwenden.

Wenn Sie also vorhaben, Python für mehr als 10.000 Punkte zu verwenden, würde ich vorschlagen, meine vektorisierte Implementierung in Betracht zu ziehen. Ich habe eine geovectorslib-Bibliothek mit einer vektorisierten Vincenty-Routine für meine eigenen Bedürfnisse erstellt, die Geographiclib als Ausweichlösung für nahe antipodale Punkte verwendet. Hier ist der Vergleich mit Geographiclib für 100k Punkte. Wie Sie sehen können, bietet es eine Verbesserung um bis zu 20x für invers und 100x für direkt Methoden für 100k Punkte und der Unterschied wird mit der Anzahl der Punkte wachsen. In Bezug auf die Genauigkeit wird sie etwa 1.e-5 rtol von Georgraphiclib betragen.

Direkte Methode für 100.000 Punkte
94,9 ms ± 25 ms pro Schleife (Mittelwert ± Std.-Abweichung von 7 Durchläufen, 1 Schleife pro Durchlauf)
9,79 s ± 1,4 s pro Schleife (Mittelwert ± Std.-Abweichung von 7 Durchläufen, 1 Schleife pro Durchlauf)

Inverse Methode für 100.000 Punkte
1,5 s ± 504 ms pro Schleife (Mittelwert ± Std.-Abweichung von 7 Durchläufen, 1 Schleife pro Durchlauf)
24,2 s ± 3,91 s pro Schleife (Mittelwert ± Std.-Abweichung von 7 Durchläufen, 1 Schleife pro Durchlauf)