Ermitteln der Entfernung zwischen zwei Punkten basierend auf Längen- und Breitengrad

155

Ich habe versucht, diese Formel zu implementieren: http://andrew.hedges.name/experiments/haversine/ Das Aplet ist gut für die beiden Punkte, die ich teste:

Geben Sie hier die Bildbeschreibung ein

Mein Code funktioniert jedoch nicht.

from math import sin, cos, sqrt, atan2

R = 6373.0

lat1 = 52.2296756
lon1 = 21.0122287
lat2 = 52.406374
lon2 = 16.9251681

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

print "Result", distance
print "Should be", 278.546

Die zurückgegebene Entfernung beträgt 5447.05546147 . Warum?

gwaramadze
quelle

Antworten:

206

Bearbeiten: Nur als Hinweis, wenn Sie nur eine schnelle und einfache Methode zum Ermitteln des Abstands zwischen zwei Punkten benötigen, empfehle ich dringend, den in Kurts Antwort unten beschriebenen Ansatz zu verwenden, anstatt Haversine erneut zu implementieren - siehe seinen Beitrag zur Begründung.

Diese Antwort konzentriert sich nur auf die Beantwortung des spezifischen Fehlers, auf den OP gestoßen ist.


Dies liegt daran, dass in Python alle Triggerfunktionen Bogenmaß und nicht Grad verwenden.

Sie können die Zahlen entweder manuell in Bogenmaß konvertieren oder die radiansFunktion aus dem Mathematikmodul verwenden:

from math import sin, cos, sqrt, atan2, radians

# approximate radius of earth in km
R = 6373.0

lat1 = radians(52.2296756)
lon1 = radians(21.0122287)
lat2 = radians(52.406374)
lon2 = radians(16.9251681)

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

print("Result:", distance)
print("Should be:", 278.546, "km")

Die Entfernung gibt jetzt den korrekten Wert von 278.545589351km zurück.

Michael0x2a
quelle
13
Dies gilt für jede Programmiersprache und auch für die Differentialrechnung. Die Verwendung von Abschlüssen ist die Ausnahme und wird nur in der menschlichen Sprache verwendet.
Bluesmoon
11
Diese Formel erfordert, dass alle Grade positiv sind. radians(abs(52.123))sollte den Trick machen ...
Richard Dunn
1
Sind Sie sicher, dass alle Grade (Winkel?) Positiv sind? Ich denke das ist falsch. Überlegen Sie, ob lat1, lon1 = 10, 10 (Grad) und lat2, lon2 = -10, -10 (Grad). Durch Hinzufügen eines abs () um die Grad wäre der Abstand Null, was falsch ist. Vielleicht wollten Sie den absoluten Wert von dlon und / oder dlat nehmen, aber wenn Sie sich den dlon ansehen, sind dlat-Werte bei der Berechnung von a, Sinus eine gerade Funktion und Cosinusquadrat eine gerade Funktion, also tue ich das nicht Sehen Sie einen Vorteil darin, einen absoluten Wert von dlat oder dlon zu nehmen.
Dave LeCompte
236

Update: 04/2018: Beachten Sie, dass die Vincenty-Entfernung seit GeoPy Version 1.13 veraltet ist - Sie sollten stattdessen geopy.distance.distance () verwenden!


Die obigen Antworten basieren auf der Haversine-Formel , die davon ausgeht, dass die Erde eine Kugel ist, was zu Fehlern von bis zu etwa 0,5% führt (gemäß help(geopy.distance)). Vincenty Distance verwendet genauere Ellipsoidmodelle wie WGS-84 und ist in Geopy implementiert . Beispielsweise,

import geopy.distance

coords_1 = (52.2296756, 21.0122287)
coords_2 = (52.406374, 16.9251681)

print geopy.distance.vincenty(coords_1, coords_2).km

druckt die Entfernung von 279.352901604Kilometern mit dem Standardellipsoid WGS-84. (Sie können auch eine .milesoder mehrere andere Entfernungseinheiten auswählen ).

Kurt Peek
quelle
1
Vielen Dank. Können Sie bitte Ihre Antwort mit den von mir angegebenen Koordinaten anstelle von Newport und Cleveland aktualisieren? Es wird zukünftigen Lesern ein besseres Verständnis geben.
Gwaramadze
1
Die willkürlichen Standorte von Newport und Cleveland stammen aus der Beispiel-Geopy-Dokumentation in der PyPI-Liste: pypi.python.org/pypi/geopy
Jason Parham
Ich musste Kurt Peeks Antwort darauf ändern: Kapitalisierung erforderlich:print geopy.distance.VincentyDistance(coords_1, coords_2).km 279.352901604
Jim
4
Sie sollten wahrscheinlich geopy.distance.distance(…)Code verwenden, der ein Alias ​​der derzeit besten (= genauesten) Entfernungsformel ist. (Vincenty im Moment.)
Geburt
9
Verwenden von geopy.distance.vincenty in Geopy-1.18.1-Ausgaben: Vincenty ist veraltet und wird in Geopy 2.0 entfernt. Verwenden Sie stattdessen geopy.distance.geodesic(oder die Standardeinstellung geopy.distance.distance), was genauer ist und immer konvergiert.
Juanmah
88

Für Leute (wie mich), die über eine Suchmaschine hierher kommen und nur nach einer Lösung suchen, die sofort funktioniert, empfehle ich die Installation mpu. Installieren Sie es über pip install mpu --userund verwenden Sie es wie folgt , um den Haversine-Abstand zu ermitteln :

import mpu

# Point one
lat1 = 52.2296756
lon1 = 21.0122287

# Point two
lat2 = 52.406374
lon2 = 16.9251681

# What you were looking for
dist = mpu.haversine_distance((lat1, lon1), (lat2, lon2))
print(dist)  # gives 278.45817507541943.

Ein alternatives Paket ist gpxpy.

Wenn Sie keine Abhängigkeiten möchten, können Sie Folgendes verwenden:

import math


def distance(origin, destination):
    """
    Calculate the Haversine distance.

    Parameters
    ----------
    origin : tuple of float
        (lat, long)
    destination : tuple of float
        (lat, long)

    Returns
    -------
    distance_in_km : float

    Examples
    --------
    >>> origin = (48.1372, 11.5756)  # Munich
    >>> destination = (52.5186, 13.4083)  # Berlin
    >>> round(distance(origin, destination), 1)
    504.2
    """
    lat1, lon1 = origin
    lat2, lon2 = destination
    radius = 6371  # km

    dlat = math.radians(lat2 - lat1)
    dlon = math.radians(lon2 - lon1)
    a = (math.sin(dlat / 2) * math.sin(dlat / 2) +
         math.cos(math.radians(lat1)) * math.cos(math.radians(lat2)) *
         math.sin(dlon / 2) * math.sin(dlon / 2))
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a))
    d = radius * c

    return d


if __name__ == '__main__':
    import doctest
    doctest.testmod()

Das andere alternative Paket ist [haversine][1]

from haversine import haversine, Unit

lyon = (45.7597, 4.8422) # (lat, lon)
paris = (48.8567, 2.3508)

haversine(lyon, paris)
>> 392.2172595594006  # in kilometers

haversine(lyon, paris, unit=Unit.MILES)
>> 243.71201856934454  # in miles

# you can also use the string abbreviation for units:
haversine(lyon, paris, unit='mi')
>> 243.71201856934454  # in miles

haversine(lyon, paris, unit=Unit.NAUTICAL_MILES)
>> 211.78037755311516  # in nautical miles

Sie behaupten, eine Leistungsoptimierung für Abstände zwischen allen Punkten in zwei Vektoren zu haben

from haversine import haversine_vector, Unit

lyon = (45.7597, 4.8422) # (lat, lon)
paris = (48.8567, 2.3508)
new_york = (40.7033962, -74.2351462)

haversine_vector([lyon, lyon], [paris, new_york], Unit.KILOMETERS)

>> array([ 392.21725956, 6163.43638211])
Martin Thoma
quelle
Gibt es eine Möglichkeit, das gegebene Highet eines der Punkte zu ändern?
Yovel Cohen
Sie können einfach den Höhenunterschied zur Entfernung hinzufügen. Das würde ich aber nicht tun.
Martin Thoma
15

Ich bin zu einer viel einfacheren und robusteren Lösung gekommen, die geodesicaus dem geopyPaket verwendet wird, da Sie sie höchstwahrscheinlich sowieso in Ihrem Projekt verwenden werden, sodass keine zusätzliche Paketinstallation erforderlich ist.

Hier ist meine Lösung:

from geopy.distance import geodesic


origin = (30.172705, 31.526725)  # (latitude, longitude) don't confuse
dist = (30.288281, 31.732326)

print(geodesic(origin, dist).meters)  # 23576.805481751613
print(geodesic(origin, dist).kilometers)  # 23.576805481751613
print(geodesic(origin, dist).miles)  # 14.64994773134371

Geopy

Ramy M. Mousa
quelle
5
import numpy as np


def Haversine(lat1,lon1,lat2,lon2, **kwarg):
    """
    This uses the ‘haversine’ formula to calculate the great-circle distance between two points – that is, 
    the shortest distance over the earth’s surface – giving an ‘as-the-crow-flies’ distance between the points 
    (ignoring any hills they fly over, of course!).
    Haversine
    formula:    a = sin²(Δφ/2) + cos φ1 ⋅ cos φ2 ⋅ sin²(Δλ/2)
    c = 2 ⋅ atan2( √a, √(1−a) )
    d = R ⋅ c
    where   φ is latitude, λ is longitude, R is earth’s radius (mean radius = 6,371km);
    note that angles need to be in radians to pass to trig functions!
    """
    R = 6371.0088
    lat1,lon1,lat2,lon2 = map(np.radians, [lat1,lon1,lat2,lon2])

    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2) **2
    c = 2 * np.arctan2(a**0.5, (1-a)**0.5)
    d = R * c
    return round(d,4)
Sudhirln92
quelle