Как рассчитать расстояние в милях для пары широты долготы в рамке данных панд? [Дубликат]

Ниже приведена векторная версия с одной и той же функцией:

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

Входы - это все массивы значений, и они должны иметь возможность делать миллионы точек мгновенно. Требование состоит в том, что входы представляют собой ndarrays, но столбцы вашей таблицы pandas будут работать.

Например, со случайно генерируемыми значениями:

>>> 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'])

Циклизация через массивы данных очень медленный в python. Numpy предоставляет функции, которые работают со всеми массивами данных, что позволяет избежать циклизации и значительно повысить производительность.

Это пример векторизации .

Чисто для иллюстративного примера, я взял версию numpy в ответ от @ballsdotballs, а также сделал реализацию C-компаньона, которая будет вызвана через ctypes. Поскольку numpy является таким высокооптимизированным инструментом, мало шансов, что мой C-код будет таким же эффективным, но он должен быть несколько близким. Большим преимуществом здесь является то, что, просматривая пример с помощью типов C, он может помочь вам понять, как вы можете подключить свои собственные функции C к Python без чрезмерных затрат. Это особенно приятно, когда вы просто хотите оптимизировать небольшой кусок большего вычисления, написав эту небольшую часть в некотором источнике C, а не в Python. Просто использование numpy решит проблему большую часть времени, но для тех случаев, когда вам действительно не нужно все numpy, и вы не хотите добавлять муфту, чтобы требовать использования типов данных numpy во всем некоторый код, очень удобно знать, как спуститься во встроенную библиотеку ctypes и сделать это самостоятельно.

Сначала давайте создадим наш исходный C-файл под названием 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");
}

Обратите внимание, что мы стараемся придерживаться конвенций C. Явно передавая аргументы данных по ссылке, используя size_t для переменной размера, и ожидая, что наша функция haversine будет работать путем изменения одного из переданных входов, так что она будет содержать ожидаемые данные при выходе. Функция фактически возвращает целое число, которое является флагом успеха / сбоя, который может использоваться другими потребителями этой функции.

Нам понадобится найти способ справиться со всеми этими маленькими C-специфическими проблемами внутри Python.

Затем давайте поместим нашу numpy версию функции вдоль с некоторыми импортными данными и некоторыми тестовыми данными в файл с именем haversine.py:

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

Я решил сделать lats и lons (в градусах), которые случайным образом выбираются между 0 и 50, но это не слишком много для этого объяснения.

Следующее, что нам нужно сделать, - это скомпилировать наш C-модуль таким образом, чтобы он мог динамически загружаться Python. Я использую систему Linux (вы можете легко найти примеры для других систем в Google), поэтому моя цель состоит в компиляции haversine.c в общий объект, например:

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

Мы можем также скомпилируйте исполняемый файл и запустите его, чтобы увидеть, что отображает функция main программы C:

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

Теперь, когда мы скомпилировали общий объект haversine.so, мы можем использовать ctypes для загрузки это в Python, и нам нужно указать путь к файлу, чтобы сделать это:

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

Теперь haversine_lib.haversine действует почти так же, как функция Python, за исключением того, что нам может понадобиться сделать какой-то ручной тип marshaling, чтобы убедиться, что входы и выходы интерпретируются правильно.

numpy на самом деле предоставляет некоторые полезные инструменты для этого, и тот, который я буду использовать здесь, - numpy.ctypeslib. Мы собираемся создать тип указателя , который позволит нам передавать numpy.ndarrays в эти ctypes -груженные функции, поскольку они были указателями. Вот код:

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] 

Обратите внимание, что мы передаем прокси-сервер функции haversine_lib.haversine для интерпретации его аргументов в соответствии с типами, которые мы хотим.

Теперь, чтобы проверить его из Python остается только создать переменную размера и массив, который будет мутировать (как и в коде C), чтобы содержать данные результата, тогда мы можем назвать это:

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

Объединяя все это в блоке __main__ в haversine.py, весь файл теперь выглядит следующим образом:

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

Чтобы запустить его, который будет запускаться и время Python и ctypes и напечатать некоторые результаты, мы можем просто сделать

python haversine.py

, который отображает:

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]

Как и ожидалось, версия numpy немного быстрее (0.11 секунды для векторов длиной 1 миллион), но наша быстрая и грязная версия ctypes не сутулится: респектабельные 0,148 секунды на одни и те же данные.

Давайте сравним это с наивным for-loop решением в 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

Когда я помещаю это в тот же файл Python, что и другие, и время на том же стане ионного элемента, я постоянно вижу время около 2,65 секунды на моей машине.

Таким образом, быстро переключаясь на ctypes, мы улучшаем скорость примерно в 18 раз. Для многих расчетов, которые могут принести пользу от доступа к голым, непрерывным данным, вы часто видите выигрыш гораздо выше, чем это.

Просто, чтобы быть предельно ясным, я вовсе не одобряю это как лучший вариант, чем просто использовать numpy. Это именно проблема, с которой numpy был создан для решения, и поэтому homebrewing ваш собственный код ctypes всякий раз, когда он (a) имеет смысл включать numpy типы данных в ваше приложение и (b) существует простой способ для сопоставления вашего кода с эквивалентом numpy не очень эффективно.

Но по-прежнему очень полезно знать, как это сделать в тех случаях, когда вы предпочитаете писать что-то в C, но назовите его в Python , или ситуации, когда зависимость от numpy непрактична (во встроенной системе, где numpy не может быть установлена, например).

В случае, если использование scikit-learn разрешено, я бы дал следующий шанс:

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

Некоторые из этих ответов «округляют» радиус земли. Если вы проверите их с другими калькуляторами расстояния (например, география ), эти функции будут отключены.

Вы можете отключить R=3959.87433 для константы преобразования ниже, если вы хотите ответьте в милях.

Если вам нужны километры, используйте 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))