Aproximación rápida de Haversine(Python/Pandas)

numpy gis (5)

Cada fila en un marco de datos de Pandas contiene coordenadas lat / lng de 2 puntos. ¡Usar el código de Python a continuación, calcular la distancia entre estos 2 puntos para muchas (millones) de filas lleva mucho tiempo!

Teniendo en cuenta que los 2 puntos están separados por menos de 50 millas y la precisión no es muy importante, ¿es posible hacer el cálculo más rápido?

from math import radians, cos, sin, asin, sqrt def haversine(lon1, lat1, lon2, lat2): """ Calculate the great circle distance between two points on the earth (specified in decimal degrees) """ # convert decimal degrees to radians lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2]) # haversine formula dlon = lon2 - lon1 dlat = lat2 - lat1 a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2 c = 2 * asin(sqrt(a)) km = 6367 * c return km for index, row in df.iterrows(): df.loc[index, ''distance''] = haversine(row[''a_longitude''], row[''a_latitude''], row[''b_longitude''], row[''b_latitude''])

Algunas de estas respuestas son "redondeando" el radio de la tierra. Si los compara con otras calculadoras de distancia (como geopy ), estas funciones estarán desactivadas.

Puede cambiar R=3959.87433 para la constante de conversión a continuación si desea obtener la respuesta en millas.

Si quieres kilómetros, usa R= 6372.8 .

lon1 = -103.548851 lat1 = 32.0004311 lon2 = -103.6041946 lat2 = 33.374939 def haversine(lat1, lon1, lat2, lon2): R = 3959.87433 # this is in miles. For Earth radius in kilometers use 6372.8 km dLat = radians(lat2 - lat1) dLon = radians(lon2 - lon1) lat1 = radians(lat1) lat2 = radians(lat2) a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2 c = 2*asin(sqrt(a)) return R * c print(haversine(lat1, lon1, lat2, lon2))

Aquí hay una versión numérica vectorizada de la misma función:

import numpy as np def haversine_np(lon1, lat1, lon2, lat2): """ Calculate the great circle distance between two points on the earth (specified in decimal degrees) All args must be of equal length. """ 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)) km = 6367 * c return km

Las entradas son todas matrices de valores, y deberían poder hacer millones de puntos al instante. El requisito es que las entradas sean ndarrays pero las columnas de tu tabla de pandas funcionarán.

Por ejemplo, con valores generados aleatoriamente:

>>> import numpy as np >>> import pandas >>> lon1, lon2, lat1, lat2 = np.random.randn(4, 1000000) >>> df = pandas.DataFrame(data={''lon1'':lon1,''lon2'':lon2,''lat1'':lat1,''lat2'':lat2}) >>> km = haversine_np(df[''lon1''],df[''lat1''],df[''lon2''],df[''lat2''])

El bucle a través de matrices de datos es muy lento en python. Numpy proporciona funciones que operan en arreglos completos de datos, lo que le permite evitar los bucles y mejorar drásticamente el rendimiento.

Este es un ejemplo de vectorization .

En caso de que se permita el uso de scikit-learn, le daría una oportunidad al siguiente:

from sklearn.neighbors import DistanceMetric dist = DistanceMetric.get_metric(''haversine'') # example data lat1, lon1 = 36.4256345, -5.1510261 lat2, lon2 = 40.4165, -3.7026 lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2]) X = [[lat1, lon1], [lat2, lon2]] kms = 6367 print(kms * dist.pairwise(X))

Simplemente por un ejemplo ilustrativo, tomé la versión numpy en la respuesta de @ballsdotballs y también hice una implementación de C complementaria para ser llamada a través de ctypes . Dado que numpy es una herramienta altamente optimizada, hay pocas posibilidades de que mi código C sea tan eficiente, pero debería ser un tanto cercano. La gran ventaja aquí es que al ejecutar un ejemplo con tipos C, puede ayudarlo a ver cómo puede conectar sus propias funciones personales de C a Python sin demasiada sobrecarga. Esto es especialmente bueno cuando solo quieres optimizar una pequeña parte de un cálculo más grande escribiendo esa pequeña parte en alguna fuente de C en lugar de Python. El simple uso de numpy solucionará el problema la mayor parte del tiempo, pero para aquellos casos en los que realmente no necesita todos los numpy y no desea agregar el acoplamiento para requerir el uso de tipos de datos numpy en algún código, es muy útil para saber cómo bajar a la biblioteca de ctypes y hacerlo usted mismo.

Primero vamos a crear nuestro archivo fuente C, llamado haversine.c :

#include <stdlib.h> #include <stdio.h> #include <math.h> int haversine(size_t n, double *lon1, double *lat1, double *lon2, double *lat2, double *kms){ if ( lon1 == NULL || lon2 == NULL || lat1 == NULL || lat2 == NULL || kms == NULL){ return -1; } double km, dlon, dlat; double iter_lon1, iter_lon2, iter_lat1, iter_lat2; double km_conversion = 2.0 * 6367.0; double degrees2radians = 3.14159/180.0; int i; for(i=0; i < n; i++){ iter_lon1 = lon1[i] * degrees2radians; iter_lat1 = lat1[i] * degrees2radians; iter_lon2 = lon2[i] * degrees2radians; iter_lat2 = lat2[i] * degrees2radians; dlon = iter_lon2 - iter_lon1; dlat = iter_lat2 - iter_lat1; km = pow(sin(dlat/2.0), 2.0) + cos(iter_lat1) * cos(iter_lat2) * pow(sin(dlon/2.0), 2.0); kms[i] = km_conversion * asin(sqrt(km)); } return 0; } // main function for testing int main(void) { double lat1[2] = {16.8, 27.4}; double lon1[2] = {8.44, 1.23}; double lat2[2] = {33.5, 20.07}; double lon2[2] = {14.88, 3.05}; double kms[2] = {0.0, 0.0}; size_t arr_size = 2; int res; res = haversine(arr_size, lon1, lat1, lon2, lat2, kms); printf("%d/n", res); int i; for (i=0; i < arr_size; i++){ printf("%3.3f, ", kms[i]); } printf("/n"); }

Tenga en cuenta que estamos tratando de mantener con las convenciones de C. Pase explícitamente argumentos de datos por referencia, usando size_t para una variable de tamaño, y esperando que nuestra función haversine funcione mutando una de las entradas pasadas de manera que contenga los datos esperados en la salida. La función en realidad devuelve un entero, que es un indicador de éxito / falla que podría ser usado por otros consumidores de nivel C de la función.

Vamos a necesitar encontrar una manera de manejar todos estos pequeños problemas específicos de C dentro de Python.

A continuación, pongamos nuestra versión numpy de la función junto con algunas importaciones y algunos datos de prueba en un archivo llamado haversine.py :

Escogí hacer lats y lons (en grados) que se eligen aleatoriamente entre 0 y 50, pero no importa demasiado para esta explicación.

Lo siguiente que debemos hacer es compilar nuestro módulo C de tal manera que Python pueda cargarlo dinámicamente. Estoy usando un sistema Linux (puedes encontrar ejemplos para otros sistemas muy fácilmente en Google), por lo que mi objetivo es compilar haversine.c en un objeto compartido, como por ejemplo:

gcc -shared -o haversine.so -fPIC haversine.c -lm

También podemos compilar en un ejecutable y ejecutarlo para ver qué funciones main del programa C muestra:

> gcc haversine.c -o haversine -lm > ./haversine 0 1964.322, 835.278,

Ahora que hemos compilado el objeto compartido haversine.so , podemos usar ctypes para cargarlo en Python y necesitamos proporcionar la ruta al archivo para hacerlo:

lib_path = "/path/to/haversine.so" # Obviously use your real path here. haversine_lib = ctypes.CDLL(lib_path)

Ahora haversine_lib.haversine actúa casi como una función de Python, excepto que es posible que tengamos que hacer algunos cálculos de tipo manual para asegurarnos de que las entradas y salidas se interpreten correctamente.

numpy realidad proporciona algunas herramientas útiles para esto y la que numpy.ctypeslib aquí es numpy.ctypeslib . Vamos a construir un tipo de puntero que nos permitirá pasar numpy.ndarrays a estas funciones ctypes ctypes, ya que a través de ellos eran punteros. Aquí está el código:

arr_1d_double = np.ctypeslib.ndpointer(dtype=np.double, ndim=1, flags=''CONTIGUOUS'') haversine_lib.haversine.restype = ctypes.c_int haversine_lib.haversine.argtypes = [ctypes.c_size_t, arr_1d_double, arr_1d_double, arr_1d_double, arr_1d_double, arr_1d_double]

Observe que le decimos a la función haversine_lib.haversine que interprete sus argumentos de acuerdo con los tipos que deseamos.

Ahora, para probarlo desde Python, lo que queda es simplemente hacer una variable de tamaño, y una matriz que será mutada (como en el código C) para contener los datos de resultados, entonces podemos llamarlo:

size = len(lat1) output = np.empty(size, dtype=np.double) print "=====" print output t2 = time.time() res = haversine_lib.haversine(size, lon1, lat1, lon2, lat2, output) t3 = time.time() print t3 - t2, res print type(output), output

Poniéndolo todo junto en el bloque __main__ de haversine.py , el archivo completo ahora se ve así:

import time import ctypes import numpy as np from math import radians, cos, sin, asin, sqrt def haversine(lon1, lat1, lon2, lat2): """ Calculate the great circle distance between two points on the earth (specified in decimal degrees) """ # convert decimal degrees to radians lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2]) # haversine formula dlon = lon2 - lon1 dlat = lat2 - lat1 a = (np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2) c = 2 * np.arcsin(np.sqrt(a)) km = 6367 * c return km if __name__ == "__main__": lat1 = 50.0 * np.random.rand(1000000) lon1 = 50.0 * np.random.rand(1000000) lat2 = 50.0 * np.random.rand(1000000) lon2 = 50.0 * np.random.rand(1000000) t0 = time.time() r1 = haversine(lon1, lat1, lon2, lat2) t1 = time.time() print t1-t0, r1 lib_path = "/home/ely/programming/python/numpy_ctypes/haversine.so" haversine_lib = ctypes.CDLL(lib_path) arr_1d_double = np.ctypeslib.ndpointer(dtype=np.double, ndim=1, flags=''CONTIGUOUS'') haversine_lib.haversine.restype = ctypes.c_int haversine_lib.haversine.argtypes = [ctypes.c_size_t, arr_1d_double, arr_1d_double, arr_1d_double, arr_1d_double, arr_1d_double] size = len(lat1) output = np.empty(size, dtype=np.double) print "=====" print output t2 = time.time() res = haversine_lib.haversine(size, lon1, lat1, lon2, lat2, output) t3 = time.time() print t3 - t2, res print type(output), output

Para ejecutarlo, que ejecutará y ctypes versiones de Python y ctypes separado e imprimirá algunos resultados, solo podemos hacer

python haversine.py

que muestra:

0.111340045929 [ 231.53695005 3042.84915093 169.5158946 ..., 1359.2656769 2686.87895954 3728.54788207] ===== [ 6.92017600e-310 2.97780954e-316 2.97780954e-316 ..., 3.20676686e-001 1.31978329e-001 5.15819721e-001] 0.148446083069 0 <type ''numpy.ndarray''> [ 231.53675618 3042.84723579 169.51575588 ..., 1359.26453029 2686.87709456 3728.54493339]

Como se esperaba, la versión numpy es ligeramente más rápida (0.11 segundos para vectores con una longitud de 1 millón), pero nuestra versión rápida y sucia de ctypes no se queda atrás: un respetable 0.148 segundos con los mismos datos.

Comparemos esto con una solución ingenua de bucle for en Python:

from math import radians, cos, sin, asin, sqrt def slow_haversine(lon1, lat1, lon2, lat2): n = len(lon1) kms = np.empty(n, dtype=np.double) for i in range(n): lon1_v, lat1_v, lon2_v, lat2_v = map( radians, [lon1[i], lat1[i], lon2[i], lat2[i]] ) dlon = lon2_v - lon1_v dlat = lat2_v - lat1_v a = (sin(dlat/2)**2 + cos(lat1_v) * cos(lat2_v) * sin(dlon/2)**2) c = 2 * asin(sqrt(a)) kms[i] = 6367 * c return kms

Cuando coloco esto en el mismo archivo de Python que los otros y cronometro en los mismos datos de millones de elementos, constantemente veo un tiempo de alrededor de 2.65 segundos en mi máquina.

Entonces, al cambiar rápidamente a ctypes , mejoramos la velocidad en un factor de aproximadamente 18. Para muchos cálculos que pueden beneficiarse del acceso a datos simples y contiguos, a menudo se ven ganancias mucho más altas incluso que esto.

Solo para ser súper claro, no estoy respaldando esto como una mejor opción que solo usar numpy . Este es precisamente el problema por el que se construyó el numpy para resolverlo, por lo que ctypes tu propio código ctypes cada vez que sea (a) tiene sentido incorporar tipos de datos numpy en tu aplicación y (b) existe una manera fácil de asignar tu código a un numpy equivalente numpy , no es muy eficiente.

Pero sigue siendo muy útil saber cómo hacer esto para aquellas ocasiones en las que prefiera escribir algo en C y llamarlo en Python, o en situaciones en numpy que no sea práctico depender de numpy (en un sistema integrado donde no se puede instalar numpy , por numpy ejemplo).

Una extensión trivial de la solución vectorizada de @ derricw , puede usar numba para mejorar el rendimiento en ~ 2x prácticamente sin cambios en su código. Para cálculos numéricos puros, esto probablemente debería usarse para pruebas comparativas / pruebas frente a soluciones posiblemente más eficientes.

from numba import njit @njit def haversine_nb(lon1, lat1, lon2, lat2): lon1, lat1, lon2, lat2 = np.radians(lon1), np.radians(lat1), np.radians(lon2), np.radians(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 return 6367 * 2 * np.arcsin(np.sqrt(a))

Benchmarking versus la función de Pandas:

%timeit haversine_pd(df[''lon1''], df[''lat1''], df[''lon2''], df[''lat2'']) # 1 loop, best of 3: 1.81 s per loop %timeit haversine_nb(df[''lon1''].values, df[''lat1''].values, df[''lon2''].values, df[''lat2''].values) # 1 loop, best of 3: 921 ms per loop

Código completo de referencia:

import pandas as pd, numpy as np from numba import njit def haversine_pd(lon1, lat1, lon2, lat2): 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 return 6367 * 2 * np.arcsin(np.sqrt(a)) @njit def haversine_nb(lon1, lat1, lon2, lat2): lon1, lat1, lon2, lat2 = np.radians(lon1), np.radians(lat1), np.radians(lon2), np.radians(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 return 6367 * 2 * np.arcsin(np.sqrt(a)) np.random.seed(0) lon1, lon2, lat1, lat2 = np.random.randn(4, 10**7) df = pd.DataFrame(data={''lon1'':lon1,''lon2'':lon2,''lat1'':lat1,''lat2'':lat2}) km = haversine_pd(df[''lon1''], df[''lat1''], df[''lon2''], df[''lat2'']) km_nb = haversine_nb(df[''lon1''].values, df[''lat1''].values, df[''lon2''].values, df[''lat2''].values) assert np.isclose(km.values, km_nb).all() %timeit haversine_pd(df[''lon1''], df[''lat1''], df[''lon2''], df[''lat2'']) # 1 loop, best of 3: 1.81 s per loop %timeit haversine_nb(df[''lon1''].values, df[''lat1''].values, df[''lon2''].values, df[''lat2''].values) # 1 loop, best of 3: 921 ms per loop