Differenz zwischen Vincenty- und Großkreisentfernungsberechnung?

Das Geopy-Paket von Python enthält zwei Methoden zur Entfernungsmessung: Great Circle und Vincenty's Formeln .

>>> from geopy.distance import great_circle
>>> from geopy.distance import vincenty
>>> p1 = (31.8300167,35.0662833) # (lat, lon) - https://goo.gl/maps/TQwDd
>>> p2 = (31.8300000,35.0708167) # (lat, lon) - https://goo.gl/maps/lHrrg
>>> vincenty(p1, p2).meters
429.16765838976664
>>> great_circle(p3, p4).meters
428.4088367903001

Was ist der Unterschied? Welche Distanzmessung wird bevorzugt?

Laut Wikipedia ist die Formel von Vincenty langsamer, aber genauer :

Vincentys Formeln sind zwei verwandte iterative Methoden, die in der Geodäsie verwendet werden, um die Entfernung zwischen zwei Punkten auf der Oberfläche eines Sphäroids zu berechnen, die von Thaddeus Vincenty (1975a) entwickelt wurden sind genauer als Methoden wie die Großkreisentfernung, die eine kugelförmige Erde voraussetzen.

Der Genauigkeitsunterschied beträgt ~0.17%in Israel 428 Meter. Ich habe einen schnellen und schmutzigen Geschwindigkeitstest durchgeführt:

<class 'geopy.distance.vincenty'>       : Total 0:00:04.125913, (0:00:00.000041 per calculation)
<class 'geopy.distance.great_circle'>   : Total 0:00:02.467479, (0:00:00.000024 per calculation)

Code:

import datetime
from geopy.distance import great_circle
from geopy.distance import vincenty
p1 = (31.8300167,35.0662833)
p2 = (31.83,35.0708167)

NUM_TESTS = 100000
for strategy in vincenty, great_circle:
    before = datetime.datetime.now()
    for i in range(NUM_TESTS):
        d=strategy(p1, p2).meters
    after = datetime.datetime.now()
    duration = after-before
    print "%-40s: Total %s, (%s per calculation)" % (strategy, duration, duration/NUM_TESTS)

Fazit: Vincentys Formel verdoppelt die Berechnungszeit im Vergleich zum Großkreis und der Genauigkeitsgewinn am getesteten Punkt beträgt ~ 0,17%.

Da die Berechnungszeit vernachlässigbar ist, wird die Formel von Vincenty für jeden praktischen Bedarf bevorzugt.

Update : Nach den aufschlussreichen Kommentaren von whuber und cffk und cffk stimme ich zu, dass der Genauigkeitsgewinn mit dem Fehler verglichen werden sollte, nicht mit der Messung. Daher ist die Formel von Vincenty um einige Größenordnungen genauer, nicht ~ 0,17%.

Wenn Sie Geopy verwenden, sind die Entfernungen great_circle und vincenty gleichermaßen bequem zu ermitteln. In diesem Fall sollten Sie fast immer dasjenige verwenden, das Ihnen das genauere Ergebnis liefert, dh vincenty. Die beiden Überlegungen (wie Sie hervorheben) sind Geschwindigkeit und Genauigkeit.

Vincenty ist zweimal langsamer. Aber wahrscheinlich ist in einer realen Anwendung die erhöhte Laufzeit vernachlässigbar. Selbst wenn Ihre Anwendung eine Million Entfernungsberechnungen erforderte, sprechen wir nur über einen Zeitunterschied von ein paar Sekunden.

Für die von Ihnen verwendeten Punkte beträgt der Fehler in vincenty 6 μm und der Fehler in der Großkreisentfernung 0,75 m. Ich würde dann sagen, dass vincenty 120000-mal genauer ist (anstatt 0,17% genauer). Für allgemeine Punkte kann der Fehler im Großkreisabstand bis zu 0,5% betragen. Können Sie also mit einem Entfernungsfehler von 0,5% leben? Für den gelegentlichen Gebrauch (wie groß ist der Abstand zwischen Kapstadt und Kairo?) Können Sie das wahrscheinlich. Viele GIS-Anwendungen haben jedoch viel strengere Genauigkeitsanforderungen. (0,5% sind 5 m über 1 km. Das macht wirklich einen Unterschied.)

Nahezu alle ernsthaften Kartierungsarbeiten werden am Referenzellipsoid durchgeführt. Daher ist es sinnvoll, Entfernungen auch am Ellipsoid zu messen. Vielleicht können Sie heute mit großen Entfernungen durchkommen. Sie müssen jedoch für jede neue Anwendung prüfen, ob dies noch akzeptabel ist. Besser ist es, nur den Ellipsoidabstand vom Start zu verwenden. Du wirst nachts besser schlafen.

ADDENDUM (Mai 2017)

Als Antwort auf die Antwort von @ craig-hicks. Die vincenty () -Methode in der Geopy weist einen potenziell schwerwiegenden Fehler auf: Sie löst einen Fehler für nahezu antipodale Punkte aus. Die Dokumentation im Code schlägt vor, die Anzahl der Iterationen zu erhöhen. Dies ist jedoch keine allgemeine Lösung, da die von vincenty () verwendete iterative Methode instabil ist für solche Punkte (jede Iteration führt Sie weiter von der richtigen Lösung ).

Warum charakterisiere ich das Problem als "potenziell tödlich"? Denn jede Verwendung der Distanzfunktion in einer anderen Softwarebibliothek muss in der Lage sein, die Ausnahme zu behandeln. Die Behandlung durch Zurückgeben eines NaN oder der Großkreisentfernung ist möglicherweise nicht zufriedenstellend, da die resultierende Entfernungsfunktion nicht der Dreieckungleichung entspricht, die ihre Verwendung beispielsweise in Aussichtspunktbäumen ausschließt.

Die Situation ist nicht ganz trostlos. Mein Python-Paket geographiclib berechnet die geodätische Entfernung fehlerfrei. Die Geopy-Pull-Anforderung Nr. 144 ändert die Entfernungsfunktion des Geopys so, dass das Paket geographiclib verwendet wird, sofern es verfügbar ist. Leider ist diese Pull-Anfrage seit August 2016 in der Schwebe.

ADDENDUM (Mai 2018)

geopy 1.13.0 verwendet jetzt das Paket geographiclib zum Berechnen von Entfernungen. Hier ist ein Beispielaufruf (basierend auf dem Beispiel in der ursprünglichen Frage):

>>> from geopy.distance import great_circle
>>> from geopy.distance import geodesic
>>> p1 = (31.8300167,35.0662833) # (lat, lon) - https://goo.gl/maps/TQwDd
>>> p2 = (31.8300000,35.0708167) # (lat, lon) - https://goo.gl/maps/lHrrg
>>> geodesic(p1, p2).meters
429.1676644986777
>>> great_circle(p1, p2).meters
428.28877358686776

Ich entschuldige mich dafür, dass ich hier eine zweite Antwort gepostet habe, aber ich nutze die Gelegenheit, um auf die Anfrage von @ craig-hicks zu antworten und Genauigkeits- und Zeitvergleiche für verschiedene Algorithmen zur Berechnung der geodätischen Entfernung bereitzustellen. Dies umschreibt einen Kommentar ich zu meinem mache Pull - Anforderung # 144 für geopy die für Geodäten der Verwendung eines von zwei Implementierungen von meinem Algorithmus ermöglicht innerhalb geopy verwendet werden soll, ist man eine nativen Python - Implementierung, geodätische (GeographicLib) , und die anderen Verwendungen eine Implementierung in C, Geodäsie (Pyproj) .

Hier sind einige Zeitdaten. Die Zeiten sind in Mikrosekunden pro Anruf angegeben

method                          dist    dest
geopy great_circle              20.4    17.1
geopy vincenty                  40.3    30.4
geopy geodesic(pyproj)          37.1    31.1
geopy geodesic(geographiclib)  302.9   124.1

Hier ist die Genauigkeit der geodätischen Berechnungen basierend auf meinem geodätischen Test-Set . Die Fehler sind in Einheiten von Mikrometern (1e-6 m) angegeben.

method                        distance destination
geopy vincenty                 205.629  141.945
geopy geodesic(pyproj)           0.007    0.013
geopy geodesic(geographiclib)    0.011    0.010

Ich habe Hannosches Pull-Request # 194 eingefügt, der einen fehlerhaften Fehler in der Zielfunktion behebt. Ohne diesen Fix beträgt der Fehler in der Zielberechnung für vincenty 8,98 Meter.

19,2% der Testfälle schlugen mit vincenty.distance fehl (Iterationen = 20). Der Testsatz ist jedoch zu Fällen verschoben, in denen dieser Fehler auftreten würde.

Mit zufälligen Punkten auf dem WGS84-Ellipsoid ist garantiert, dass der Vincenty-Algorithmus bei 16,6 von 10.000 fehlschlägt (die richtige Lösung ist ein instabiler Fixpunkt der Vincenty-Methode).

Bei der Geopy-Implementierung von Vincenty und Iterationen = 20 beträgt die Fehlerrate 82,8 pro 1000000. Bei Iterationen = 200 beträgt die Fehlerrate 21,2 pro 1000000.

Auch wenn diese Raten gering sind, können Ausfälle durchaus üblich sein. Zum Beispiel würde in einem Datensatz mit 1000 zufälligen Punkten (denken Sie vielleicht an die Flughäfen der Welt) die Berechnung der Matrix für die vollständige Entfernung durchschnittlich 16-mal fehlschlagen (mit Iterationen = 20).

Es scheint, dass das Paket geopy.distance eine Funktion "distance ()" bietet, die standardmäßig vincenty () ist. Ich würde grundsätzlich die Verwendung von distance () empfehlen, da dies die Paketempfehlung ist, falls in Zukunft jemals von vincenty () abgewichen wird (so unwahrscheinlich ist das). Weiterlesen:

Dieser Dokumentationshinweis ist im Quellcode für die von Ihnen angegebene Funktion vincenty () enthalten:

Hinweis: Diese Implementierung von Vincenty distance konvergiert nicht für einige gültige Punkte. In einigen Fällen kann ein Ergebnis erzielt werden, indem die Anzahl der Iterationen erhöht wird ( iterationsSchlüsselwortargument in der Klasse __init__mit einem Standardwert von 20). Es kann vorzuziehen sein: class: zu verwenden .great_circle, was geringfügig ungenauer ist, aber immer ein Ergebnis liefert.

Den Quellcode mit dem obigen Kommentar / Hinweis finden Sie unter https://github.com/geopy/geopy/blob/master/geopy/distance.py Scrollen Sie nach unten zur Definition für vincenty ().

Die voreingestellte Distanzfunktion, die von diesem Paket beim Abrufen von distance () verwendet wird, ist die Funktion vincenty (), die impliziert, dass das Nichtkonvergieren nicht katastrophal ist und eine angemessene Antwort zurückgegeben wird - am wichtigsten ist, dass keine Ausnahme generiert wird.

Update: Wie von "cffk" bemerkt, löst die Funktion vincenty () explizit eine ValueError-Ausnahme aus, wenn der Algorithmus nicht konvergiert - obwohl dies nicht in der Funktionsbeschreibung dokumentiert ist. Daher ist die Dokumentation fehlerhaft.

Unabhängig davon, ob Sie vincenty oder haversine oder das sphärische Gesetz des Cosinus verwenden, ist es ratsam, sich potenzieller Probleme mit dem Code, den Sie verwenden möchten, bewusst zu werden, Dinge zu beachten und abzumildern und wie man mit vincenty vs. haversine vs. sloc-Problemen umgeht werden sich unterscheiden, wenn man sich der lauernden Probleme / Edgecases eines jeden bewusst wird, die im Allgemeinen bekannt sein können oder nicht. Der erfahrene Programmierer weiß das. Neulinge dürfen nicht. Ich hoffe, einige von ihnen vor Frustrationen zu bewahren, wenn ein Ausschnitt aus einem Forum in bestimmten Fällen etwas Unerwartetes bewirkt. Wenn man ernsthaft eine dieser Versionen verwenden möchte, wie zum Beispiel Vincent, Haversine, Sloc, dann haben SE, SO, Reddit, Quora usw. möglicherweise nur eingeschränkte Hilfe bei der anfänglichen Codierung einer Lösung bereitgestellt, aber das bedeutet nicht, dass dies der Fall ist Ihre Lösung oder akzeptierte "Antwort" ist frei von Problemen. Wenn ein Projekt wichtig genug ist, verdient es einen angemessenen Forschungsaufwand. Lesen Sie das Handbuch, lesen Sie die Dokumentation, und falls eine Codeüberprüfung dieses Codes vorliegt, lesen Sie diese. Das Kopieren und Einfügen eines Snippets oder einer Zusammenfassung, das oder die hundertmal oder öfter aktualisiert wurde, bedeutet nicht, dass die Sicherheit umfassend und gewährleistet ist.

Die faszinierende Antwort von cffk macht deutlich, dass es sich bei verpackten Lösungen um lauernde Edgecases handelt, die Ausnahmen oder andere Schwierigkeiten hervorrufen können . Die konkreten Behauptungen, die in diesem Beitrag gemacht werden, übersteigen derzeit mein Zeitbudget, aber ich gehe davon aus, dass in bestimmten Paketen tatsächlich Probleme lauern, einschließlich mindestens einer einzigen Implementierung, für die mindestens eine Person Verbesserungsvorschläge gemacht hat auf die eine oder andere Weise, um das Risiko, auf diese Schwierigkeiten zu stoßen, zu minimieren oder zu beseitigen. Ich werde nicht weiter auf das Thema von Vincent eingehen (da ich es viel zu ignoriert habe), sondern mich stattdessen dem Thema Haversine zuwenden, zumindest teilweise dem Thema mit dem OP.

Die populär veröffentlichte Haversine-Formel, ob in Python oder einer anderen Sprache, da sie wahrscheinlich die IEEE 754-Gleitkomma-Spezifikation für die meisten aktuellen Intel- und Intel-ähnlichen Systeme und ARM-Prozessoren, PowerPCs usw. verwenden wird auch anfällig für seltene, aber reale und wiederholbare Ausnahmefehler in der Nähe oder in einem Bogenabstand von 180 Grad, Antipodenpunkte aufgrund von Gleitkommanäherungen und Abrundungen. Einige Neulinge sind möglicherweise noch nicht von dieser Situation gebissen worden. Da sich diese fp-Spezifikation annähert und rundet, bedeutet dies nicht, dass jeder Code, der fp64 aufruft, Ausnahmefehler verursachen kann, nein. Aber etwas Code, Einige Formeln haben möglicherweise nicht so offensichtliche Randbedingungen, bei denen die Näherungen und Rundungen von IEEE 754 fp64 dazu führen können, dass ein Wert leicht aus dem Bereich einer mathematischen Methode herausfällt, von der erwartet wird, dass sie einen solchen Wert fehlerfrei auswertet. Ein Beispiel ... sqrt (). Wenn ein negativer Wert in ein sqrt () eingeht, z. B. sqrt (-0.00000000000000000122739), liegt ein Ausnahmefehler vor. In der Haversin-Formel gibt es in atan2 () zwei sqrt () -Methoden, um eine Lösung zu finden. Dasa , das in sqrt () berechnet und dann verwendet wird, kann an den Antipodenpunkten auf dem Globus leicht unter 0,0 oder über 1,0 abweichen, sehr leicht aufgrund von fp64-Approximationen und -Rundungen, selten, aber wiederholbar. Die konsistente, zuverlässige Wiederholbarkeit macht dies in diesem Zusammenhang zu einem Ausnahmerisiko, zu einem Edgecase zum Schutz, zur Minderung und nicht zu einem isolierten, zufälligen Zufall. Hier ist ein Beispiel für ein kurzes Python3-Snippet von Haversine ohne den erforderlichen Schutz:

import math as m

a = m.sin(dlat / 2)**2 + m.cos(lat1) * m.cos(lat2) * m.sin(dlon / 2)**2
c = 2 * m.atan2(m.sqrt(a), m.sqrt(1 - a))
distance = Radius * c

Sehr nahe oder an antipodalen Punkten, a in der ersten Zeile der Formel berechnetes negativ sein, selten, aber wiederholt mit denselben Lat-Lon-Koordinaten. Zum Schutz / diese seltenen Ereignisse zu korrigieren, kann man einfach hinzufügen, nach dem einer Berechnung, wie unten zu sehen:

import math as m

note = ''

a = m.sin(dlat / 2)**2 + m.cos(lat1) * m.cos(lat2) * m.sin(dlon / 2)**2
if a < 0.0: a = 0.0 ; note = '*'
if a > 1.0: a = 1.0 ; note = '**'
c = 2 * m.atan2(m.sqrt(a), m.sqrt(1 - a))
distance = Radius * c

# note = '*'  # a went below 0.0 and was normalized back to 0.0
# note = '**' # a went above 1.0 and was normalized back to max of 1.0

Natürlich habe ich hier nicht die gesamte Funktion gezeigt, sondern einen kurzen Ausschnitt, wie er so oft gepostet wird. Aber dieses Beispiel zeigt den Schutz für sqrt (), indem es das a testet und es gegebenenfalls normalisiert, wodurch auch die Notwendigkeit vermieden wird, das Ganze mit Ausnahme von "try" zu testen . Das note = '' up top soll verhindern, dass die Bytecode-Stufe gegen die Verwendung dieser Note protestiert, bevor ihr ein Wert zugewiesen wird, wenn sie mit dem Ergebnis der Funktion zurückgegeben wird.

Mit dieser einfachen Änderung des Hinzufügens der beiden a- Tests sind die Funktionen von sqrt () zufriedenstellend, und der Code enthält jetzt eine zusätzliche Notiz , die an den aufrufenden Code zurückgegeben werden kann, um darauf hinzuweisen , dass ein Ergebnis leicht normalisiert wurde, und warum. Einige kümmern sich vielleicht darum, andere kümmern sich vielleicht nicht darum, aber es verhindert einen Ausnahmefehler, der andernfalls auftreten kann. Ein try except-Block fängt die Ausnahme möglicherweise ab, behebt sie jedoch nicht, es sei denn, er wurde ausdrücklich dazu geschrieben. Es scheint einfacher zu sein, die Korrekturzeile (n) unmittelbar nach der a- Berechnungszeile zu codieren . Eine gründlich bereinigte Eingabe sollte dann keinen Versuch erfordern, außer hier überhaupt zu blockieren.

Zusammenfassung: Wenn Sie Haversine verwenden, das explizit codiert ist, anstatt ein Paket oder eine Bibliothek zu verwenden, ist es unabhängig von der Sprache Ihrer Wahl eine gute Idee, es zu testen und zu normalisieren , einen zurück in den nötigen Bereich von 0,0 <= a <= 1,0 um um die nächste Zeile mit ihren c- Berechnungenzu schützen. Aber die Mehrheit der Haversine-Code-Schnipsel zeigt es nicht an und erwähnt das Risiko nicht.

Erfahrung: Während gründlicher Tests rund um den Globus habe ich in Schritten von 0,001 Grad eine Festplatte mit Lat-Lon-Kombinationen gefüllt, die eine Ausnahme verursacht haben, eine zuverlässig konsistente, wiederholbare Ausnahme, während eines Monats, in dem auch die Zuverlässigkeit der CPU-Kühlung überprüft wurde Fan und meine Geduld. Ja, ich habe seitdem die meisten dieser Protokolle gelöscht, da ihr Zweck hauptsächlich darin bestand, den Sinn zu beweisen (wenn das Wortspiel erlaubt ist). Aber ich habe einige kürzere Protokolle von 'Problem-Lat-Lon-Werten', die zu Testzwecken aufbewahrt werden.

Genauigkeit: Will Verliert a und das gesamte Haversine-Ergebnis an Genauigkeit, wenn es wieder ein wenig in die Domäne normalisiert wird? Nicht viel, vielleicht nicht mehr als die bereits eingeführten fp64-Annäherungen und -Rundungen, die zu dieser leichten Abweichung aus dem Bereich führten. Wenn Sie festgestellt haben, dass Haversine bereits über hundert Mal akzeptabel ist - einfacher, schneller, einfacher anzupassen, zu beheben und zu warten -, ist Haversine möglicherweise eine gute Lösung für Ihr Projekt.

Ich habe Haversine auf einer über Kopf projizierten Himmelskugel verwendet, um Winkelabstände zwischen Objekten am Himmel zu messen, wie von einer Position auf der Erde aus gesehen, Azimut und Alt auf Lat-Lon-Äquivalent-Koordinaten des Himmels abgebildet Die projizierte theoretische Himmelssphäre ist eine perfekte Kugel, wenn es darum geht, Blickwinkel zwischen zwei Objekten von einer Position auf der Erdoberfläche aus zu messen. Es passt perfekt zu meinen Bedürfnissen. Also, Haversine ist immer noch sehr nützlich und in bestimmten Anwendungen sehr genau (genau für meine Zwecke) ... aber wenn Sie es verwenden, ob auf der Erde für GIS oder Navigation oder bei Beobachtungen und Messungen von Himmelsobjekten, schützen Sie es es im Fall von Antipodenpunkten oder sehr nahen Antipodenpunkten durch Testen von aund stupsen Sie es zurück in seine benötigte Domäne, wenn nötig.

Das ungeschützte Haversine ist im gesamten Internet zu finden, und ich habe nur einen alten Usenet-Beitrag gesehen, der, wie ich glaube, vor jemandem von JPL geschützt war und der möglicherweise vor 1985 laut Gleitkomma-Spezifikation IEEE 754 erstellt wurde. Zwei andere Seiten erwähnten mögliche Probleme in der Nähe von Antipodenpunkten, beschrieben diese jedoch nicht oder wie man sie mildern könnte. Aus diesem Grund gibt es Bedenken für Neulinge (wie mich), die die bewährten Methoden möglicherweise nicht immer gut genug verstehen, um einige Codes, die sie kopiert und in ein vertrauenswürdiges Projekt eingefügt haben, weiter zu untersuchen und zu testen. Der faszinierende Beitrag von cffk war insofern erfrischend, als er öffentlich mit diesen Arten von Problemen bekannt war, die nicht oft erwähnt, selten öffentlich zum Schutz in Snippets codiert und selten auf diese Weise diskutiert wurden, verglichen mit der Menge ungeschützter und nicht diskutierter Versionen, die veröffentlicht wurden.

Ab 20190923 wird auf der Wiki-Seite für Haversine-Formeln in der Tat das Problem erwähnt, das aufgrund von Gleitkommaproblemen in Computergeräten an Antipodenpunkten auftreten kann ... ermutigend ...

https://en.wikipedia.org/wiki/Haversine_formula

(Da diese Wiki-Seite zu diesem Zeitpunkt keinen HTML-Anker für den Abschnitt hat, auf den ich direkt verweisen würde, führen Sie nach dem Laden der Seite eine Suche auf dieser Browserseite nach "Wenn Sie diese Formeln verwenden" durch, und Sie werden es tun siehe das Problem der Haversine mit den erwähnten Antipodenpunkten, genauer gesagt.)

Und diese andere Seite hat auch eine sehr kurze Erwähnung davon:

https://www.movable-type.de/scripts/latlong.html

Wenn man auf dieser Seite nach "Schutz vor Rundungsfehlern" sucht, gibt es diese ...

Wenn atan2 nicht verfügbar ist, kann c aus 2 Asin (min (1, √a)) berechnet werden (einschließlich des Schutzes vor Rundungsfehlern).

Jetzt gibt es einen seltenen Fall, in dem Rundungsfehler erwähnt und der Schutz für die Version asin () angezeigt wird, jedoch für die Version atan2 () nicht erwähnt oder angezeigt wird. Es wird aber zumindest auf die Gefahr von Rundungsfehlern hingewiesen.

imho, jede 24/7/365 Anwendung mit Haversine, benötigt diesen Schutz in der Nähe der Antipodenpunkte als wichtiges und einfaches Detail.

Ich weiß nicht, welche haversine-Pakete diesen Schutz enthalten oder nicht, aber wenn Sie mit all dem noch nicht vertraut sind und die im Fachjargon veröffentlichten "Snippet" -Versionen verwenden möchten, wissen Sie jetzt, dass sie geschützt werden müssen, und Dieser Schutz ist sehr einfach zu implementieren, d. h., wenn Sie nicht vincenty und keinen verpackten Haversine ohne einfachen Zugriff zum Ändern des Code des Pakets verwenden.

IOW, egal ob mit vincenty oder haversine oder sloc, man sollte sich über Probleme mit dem Code im Klaren sein, Dinge, auf die man achten und die man mindern muss, und wie man mit vincenty gegen haversine gegen sloc umgeht, wird sich unterscheiden, wenn man sich seiner bewusst wird Lauernde Probleme / Edgecases, die allgemein bekannt sein können oder nicht.