Obtener la distancia entre dos puntos en función de la latitud / longitud

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

introduzca 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 devuelve es 5447.05546147. ¿Por qué?

Author: Martin Thoma, 2013-10-16
Source

3 answers

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

Puede convertir los números manualmente a radianes, o usar el radians función del módulo de matemáticas: 

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.545589351 km.

Edit: Como nota, si te has topado con este post porque solo necesitas una forma rápida y fácil de encontrar la distancia entre dos puntos, te recomiendo que utilices el enfoque recomendado en La respuesta de Kurt a continuación see ver su post para la justificación.

 100
Author: Michael0x2a,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2018-09-29 01:06:45

Actualización: 04/2018: Tenga en cuenta que la distancia de Vincenty está obsoleta desde la versión de GeoPy 1.13 - deberías usar Geopy.distancia.distance() en su lugar!

Las respuestas anteriores se basan en la fórmula Haversine, que asume que la tierra es una esfera, lo que resulta en errores de hasta aproximadamente 0.5% (de acuerdo con help(geopy.distance)). Vincenty distanceutiliza modelos elipsoidales más precisos como WGS-84, y se implementa en geopy. Para 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

Mostrará la distancia de 279.352901604 kilómetros usando el elipsoide predeterminado WGS-84. (También puede elegir .miles o una de varias otras unidades de distancia).

 73
Author: Kurt Peek,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2018-04-27 10:14:37

Para las personas (como yo) que vienen aquí a través de un motor de búsqueda y simplemente buscan una solución que funcione fuera de la caja, recomiendo instalar mpu. Instálalo a través de pip install mpu --user y úsalo así para obtener la distancia de 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()
 47
Author: Martin Thoma,
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2018-05-16 20:57:11