Python. Как быстро посчитать расстояние между всеми точками в массивах?

#python #массивы #scipy #spatial #haversine

Есть 733  географических точек F (массив 2*733)
Есть 22 000 географических точек S (массив 2*22 000)

Как быстро посчитать расстояние от каждой точки из F до каждой точки S.

Этот код работает очень долго, около семи минут. Почему?

import scipy.spatial
D = scipy.spatial.distance.cdist(F,S,lambda u, v: haversine(u, v))


А если вычислять евклидову метрику, то быстро.

import scipy.spatial
D2 = scipy.spatial.distance.cdist(F,S,'euclidean')


Вычислять расстояние функцией  тоже долго:

import scipy.spatial
D = scipy.spatial.distance.cdist(F,S,lambda u,v:geopy.distance.vincenty(u,v).km) 


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

import scipy.spatial
D = scipy.spatial.distance.cdist(F,S,lambda  u,v:(u-v).sum())

    

Ответы

Ответ 1

Вот рабочий пример.

Данные для примера: http://samplecsvs.s3.amazonaws.com/SacramentocrimeJanuary2006.csv

import numpy as np
import pandas as pd
from sklearn.neighbors import DistanceMetric

# http://scikit-learn.org/stable/modules/generated/sklearn.neighbors.DistanceMetric.html
dist = DistanceMetric.get_metric('haversine')

url = 'http://samplecsvs.s3.amazonaws.com/SacramentocrimeJanuary2006.csv'
df = pd.read_csv(url, usecols=['latitude','longitude'])


первые 5 записей:

In [213]: df.head()
Out[213]:
    latitude   longitude
0  38.550420 -121.391416
1  38.473501 -121.490186
2  38.657846 -121.462101
3  38.506774 -121.426951
4  38.637448 -121.384613


увеличим DataFrame в три раза, чтобы получить примерно такой же объем данных как
у вас в вопросе:

df = pd.concat([df] * 3, ignore_index=True).sample(frac=1)
# split our data set into two
F, S = np.split(df, [int(len(df)*0.033)])


получилось два DataFrame'а:

print('shape of [F]: {}'.format(F.shape))
print('shape of [S]: {}'.format(S.shape))


выдает:

shape of [F]: (750, 2)
shape of [S]: (22002, 2)


считаем расстояние:

earth_radius = 6371

D = dist.pairwise(np.radians(F), np.radians(S)) * earth_radius


PS haversine считается по формуле:

2 * arcsin(sqrt(sin^2(0.5*dx) * cos(x1) * cos(x2) * sin^2(0.5*dy)))


и возвращает результат в радианах, который надо умножить на радиус Земли (в км.),
чтобы получить расстояние в километрах

Результат:

In [208]: D
Out[208]:
array([[  3.86258192,  12.14112631,  16.1702106 , ...,  11.2399247 ,  13.71076021,
  2.59575523],
       [  3.19765341,  15.72580688,  20.35963524, ...,  15.38247654,  17.85975202,
  4.60552864],
       [  6.15251463,  19.94251384,  19.77488568, ...,  17.0637875 ,  19.15736473,
  6.55354701],
       ...,
       [  7.85715424,   9.52674314,  11.13189495, ...,   6.81544901,   9.16104377,
  5.58452969],
       [ 14.02227805,   5.08370143,   6.47846773, ...,   0.57188025,   3.04569744,
 11.81577747],
       [ 17.84804533,   3.93791815,   6.91865107, ...,   3.8743883 ,   2.94376778,
 15.72668196]])

In [209]: D.shape
Out[209]: (750, 22002)


замер скорости - около 6.5 секунд (для 16.5 миллионов пар точек) на моем стареньком
ноуте:

In [210]: %timeit dist.pairwise(np.radians(F), np.radians(S)) * earth_radius
6.45 s ± 156 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)


Проверка формулы:

In [204]: lyon = np.array([[45.7597, 4.8422]])

In [205]: paris = np.array([[48.8567, 2.3508]])

In [206]: dist.pairwise(np.radians(lyon), np.radians(paris))  * earth_radius
Out[206]: array([[ 392.21671781]])