Obtener distancia entre dos puntos según la latitud / longitud

158

Intenté implementar esta fórmula: http://andrew.hedges.name/experiments/haversine/ El aplet funciona bien para los dos puntos que estoy probando:

ingrese la descripción de la imagen aquí

Sin embargo, mi código no funciona.

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

La distancia que regresa es 5447.05546147 . ¿Por qué?

gwaramadze
fuente

Respuestas:

207

Editar: solo como una nota, si solo necesita una forma rápida y fácil de encontrar la distancia entre dos puntos, recomiendo encarecidamente utilizar el enfoque descrito en la respuesta de Kurt a continuación en lugar de volver a implementar Haversine: consulte su publicación para obtener una justificación.

Esta respuesta se centra solo en responder al error específico con el que se encontró OP.


Es porque en Python, todas las funciones trigonométricas usan radianes , no grados.

Puede convertir los números manualmente a radianes o usar la radiansfunción del módulo matemático:

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

La distancia ahora devuelve el valor correcto de 278.545589351km.

Michael0x2a
fuente
13
Esto es cierto en cualquier lenguaje de programación, y también en cálculo diferencial. usar grados es la excepción, y solo se usa en el habla humana.
bluesmoon
11
Palabra para los sabios, esta fórmula requiere que todos los grados sean positivos. radians(abs(52.123))debería hacer el truco ...
Richard Dunn
1
¿Estás seguro de que todos los grados (ángulos) son positivos? Creo que esto está mal. Considere si lat1, lon1 = 10, 10 (grados) y lat2, lon2 = -10, -10 (grados). Al agregar un abs () alrededor de los grados, la distancia sería cero, lo cual es incorrecto. Quizás quisiste tomar el valor absoluto de dlon y / o dlat, pero si miras los valores de dlon, dlat en el cálculo de a, el seno es una función par, y el coseno al cuadrado es una función par, así que no vea cualquier beneficio al tomar un valor absoluto de dlat o dlon, tampoco.
Dave LeCompte
238

Actualización: 04/2018: tenga en cuenta que la distancia de Vincenty está en desuso desde la versión 1.13 de GeoPy ; en su lugar, debe usar geopy.distance.distance ().


Las respuestas anteriores se basan en la fórmula de Haversine , que supone que la tierra es una esfera, lo que da como resultado errores de hasta aproximadamente 0.5% (según help(geopy.distance)). La distancia de Vincenty utiliza modelos elipsoidales más precisos, como WGS-84 , y se implementa en geopy . Por ejemplo,

import geopy.distance

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

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

imprimirá la distancia de 279.352901604kilómetros utilizando el elipsoide predeterminado WGS-84. (También puede elegir .mileso una de varias otras unidades de distancia).

Kurt Peek
fuente
1
Gracias. ¿Puede actualizar su respuesta con las coordenadas que proporcioné en cuestión en lugar de Newport y Cleveland? Le dará una mejor comprensión a los futuros lectores.
gwaramadze 05 de
1
Las ubicaciones arbitrarias de Newport y Cleveland provienen de la documentación de geopy de ejemplo en la lista de PyPI: pypi.python.org/pypi/geopy
Jason Parham
Tuve que modificar la respuesta de Kurt Peek a esto: Capitalización requerida:print geopy.distance.VincentyDistance(coords_1, coords_2).km 279.352901604
Jim
44
Probablemente deberías usar un geopy.distance.distance(…)código que sea un alias de la fórmula de distancia mejor (= más precisa) actualmente. (Vincenty en este momento.)
mbirth
10
Uso de geopy.distance.vincenty en salidas geopy-1.18.1: Vincenty está en desuso y se eliminará en geopy 2.0. Use geopy.distance.geodesic(o el valor predeterminado geopy.distance.distance) en su lugar, que es más preciso y siempre converge.
juanmah
88

Para las personas (como yo) que vienen aquí a través del motor de búsqueda y solo buscan una solución que funcione de inmediato, recomiendo la instalación mpu. Instálelo a través de pip install mpu --usery utilícelo de esta manera para obtener la distancia haversine :

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.

Un paquete alternativo es gpxpy.

Si no quieres dependencias, puedes usar:

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

El otro paquete alternativo es [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

Afirman tener una optimización del rendimiento para distancias entre todos los puntos en dos vectores

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
fuente
¿Hay alguna manera de cambiar el Highet dado de uno de los puntos?
Yovel Cohen
Simplemente puede agregar la diferencia de altura a la distancia. Sin embargo, no haría eso.
Martin Thoma
16

Llegué a una solución mucho más simple y robusta que está usando geodesicdesde el geopypaquete, ya que es muy probable que la uses en tu proyecto de todos modos, por lo que no es necesario instalar ningún paquete adicional.

Aquí está mi solución:

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
fuente
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
fuente
0

Hay varias formas de calcular la distancia en función de las coordenadas, es decir, latitud y longitud.

Instalar e importar

from geopy import distance
from math import sin, cos, sqrt, atan2, radians
from sklearn.neighbors import DistanceMetric
import osrm
import numpy as np

Definir coordenadas

lat1, lon1, lat2, lon2, R = 20.9467,72.9520, 21.1702, 72.8311, 6373.0
coordinates_from = [lat1, lon1]
coordinates_to = [lat2, lon2]

Usando haversine

dlon = radians(lon2) - radians(lon1)
dlat = radians(lat2) - radians(lat1)
    
a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
c = 2 * atan2(sqrt(a), sqrt(1 - a))
    
distance_haversine_formula = R * c
print('distance using haversine formula: ', distance_haversine_formula)

Usando haversine con sklearn

dist = DistanceMetric.get_metric('haversine')
    
X = [[radians(lat1), radians(lon1)], [radians(lat2), radians(lon2)]]
distance_sklearn = R * dist.pairwise(X)
print('distance using sklearn: ', np.array(distance_sklearn).item(1))

Usando OSRM

osrm_client = osrm.Client(host='http://router.project-osrm.org')
coordinates_osrm = [[lon1, lat1], [lon2, lat2]] # note that order is lon, lat
    
osrm_response = osrm_client.route(coordinates=coordinates_osrm, overview=osrm.overview.full)
dist_osrm = osrm_response.get('routes')[0].get('distance')/1000 # in km
print('distance using OSRM: ', dist_osrm)

Usando geopy

distance_geopy = distance.distance(coordinates_from, coordinates_to).km
print('distance using geopy: ', distance_geopy)
    
distance_geopy_great_circle = distance.great_circle(coordinates_from, coordinates_to).km 
print('distance using geopy great circle: ', distance_geopy_great_circle)

Salida

distance using haversine formula:  26.07547017310917
distance using sklearn:  27.847882224769783
distance using OSRM:  33.091699999999996
distance using geopy:  27.7528030550408
distance using geopy great circle:  27.839182219511834
Patel Romil
fuente