Расчет расстояния между параболическими отводами

Question

У меня есть два массива с широким и длинным. Я хочу рассчитать расстояние между каждой парой lat и long с каждой другой парой lat и long в массиве. Вот мои два массива.

lat_array 

array([ 0.33356456, 0.33355585, 0.33355585, 0.33401788, 0.33370132, 
     0.33370132, 0.33370132, 0.33371075, 0.33371075, 0.33370132, 
     0.33370132, 0.33370132, 0.33356488, 0.33356488, 0.33370132, 
     0.33370132, 0.33370132, 0.33401788, 0.33362632, 0.33362632, 
     0.33364007, 0.33370132, 0.33401788, 0.33401788, 0.33358399, 
     0.33358399, 0.33358399, 0.33370132, 0.33370132, 0.33362632, 
     0.33370132, 0.33370132, 0.33370132, 0.33370132, 0.33370132, 
     0.33356488, 0.33356456, 0.33391071, 0.33370132, 0.33356488, 
     0.33356488, 0.33356456, 0.33356456, 0.33356456, 0.33362632, 
     0.33364804, 0.3336314 , 0.33370132, 0.33370132, 0.33370132, 
     0.33364034, 0.33359921, 0.33370132, 0.33360397, 0.33348863, 
     0.33370132]) 
long_array 

array([ 1.27253229, 1.27249141, 1.27249141, 1.27259085, 1.2724337 , 
     1.2724337 , 1.2724337 , 1.27246931, 1.27246931, 1.2724337 , 
     1.2724337 , 1.2724337 , 1.27254305, 1.27254305, 1.2724337 , 
     1.2724337 , 1.2724337 , 1.27259085, 1.27250461, 1.27250461, 
     1.27251211, 1.2724337 , 1.27259085, 1.27259085, 1.27252134, 
     1.27252134, 1.27252134, 1.2724337 , 1.2724337 , 1.27250461, 
     1.2724337 , 1.2724337 , 1.2724337 , 1.2724337 , 1.2724337 , 
     1.27254305, 1.27253229, 1.27266808, 1.2724337 , 1.27254305, 
     1.27254305, 1.27253229, 1.27253229, 1.27253229, 1.27250461, 
     1.27250534, 1.27250184, 1.2724337 , 1.2724337 , 1.2724337 , 
     1.27251339, 1.27223739, 1.2724337 , 1.2722575 , 1.27237575, 
     1.2724337 ]) 

После преобразования в радианы. Теперь я хочу расстояние между первой парой lat и длинными с оставшимися парами lat и long и так далее. И хотите напечатать пары и соответствующее расстояние.

Это то, что я делаю в python.

distance = [] 
R = 6371.0 

for i in range(len(lat_array)): 
    for j in (i+1,len(lat_array)): 
     dlon = long_array[j]-long_array[i] 
     dlat = lat_array[j]-lat_array[i] 
     a = sin(dlat/2)**2 + cos(lat_array[i]) * cos(lat_array[j]) *  
      sin(dlon/2)**2 
     c = 2 * atan2(sqrt(a), sqrt(1 - a)) 

     distance.append(R * c) 

Это дает мне ошибку IndexError: index 56 is out of bounds for axis 0 with size 56 Где я делаю неправильно? И как сделать вычисление быстрее, если массив большой? Пожалуйста помоги.

Answer 1

У вас есть опечатка в коде. Изменение

for j in (i+1,len(lat_array)): 

в

for j in range(i+1,len(lat_array)): 

В противном случае вы переборе кортеж, состоящий из двух элементов i+1 и len(lat_array). Вторая причина ошибки.

Answer 2

Предполагая lat и lng как широты & долгот массивов и те имеют данные в радианах, вот один Векторизованное решения, основанное на this other solution -

# Elementwise differentiations for lattitudes & longitudes 
dflat = lat[:,None] - lat 
dflng = lng[:,None] - lng 

# Finally Calculate haversine using its distance formula 
d = np.sin(dflat/2)**2 + np.cos(lat[:,None])*np.cos(lat) * np.sin(dflng/2)**2 
hav_dists = 2 * 6371 * np.arcsin(np.sqrt(d)) 

Теперь, вышеупомянутый подход дал бы нам выход для всех пар независимо их порядка. Таким образом, у нас были бы два выходных выхода для двух пар: (point1,point2) & (point2,point1), хотя расстояния были бы одинаковыми. Таким образом, чтобы сохранить в памяти и, надеюсь, более высокая производительность, вы можете создавать уникальные парные идентификаторы с np.triu_indices и изменить ранее перечисленный подход, как так -

# Elementwise differentiations for lattitudes & longitudes, 
# but not repeat for the same paired elements 
N = lat.size 
idx1,idx2 = np.triu_indices(N,1) 
dflat = lat[idx2] - lat[idx1] 
dflng = lng[idx2] - lng[idx1] 

# Finally Calculate haversine using its distance formula 
d = np.sin(dflat/2)**2 + np.cos(lat[idx2])*np.cos(lat[idx1]) * np.sin(dflng/2)**2 
hav_dists = 2 * 6371 * np.arcsin(np.sqrt(d)) 

---

Функция определение -

def original_app(lat,lng): 
    distance = [] 
    R = 6371.0 
    for i in range(len(lat)): 
     for j in range(i+1,len(lat)): 
      dlon = lng[j]-lng[i] 
      dlat = lat[j]-lat[i] 
      a = sin(dlat/2)**2 + cos(lat[i]) * cos(lat[j]) * sin(dlon/2)**2 
      c = 2 * atan2(sqrt(a), sqrt(1 - a)) 
      distance.append(R * c) 
    return distance 

def vectorized_app1(lat,lng):        
    dflat = lat[:,None] - lat 
    dflng = lng[:,None] - lng 
    d = np.sin(dflat/2)**2 + np.cos(lat[:,None])*np.cos(lat) * np.sin(dflng/2)**2 
    return 2 * 6371 * np.arcsin(np.sqrt(d)) 

def vectorized_app2(lat,lng):        
    N = lat.size 
    idx1,idx2 = np.triu_indices(N,1) 
    dflat = lat[idx2] - lat[idx1] 
    dflng = lng[idx2] - lng[idx1] 
    d =np.sin(dflat/2)**2+np.cos(lat[idx2])*np.cos(lat[idx1])*np.sin(dflng/2)**2 
    return 2 * 6371 * np.arcsin(np.sqrt(d)) 

Проверьте выход - тест

In [78]: lat 
Out[78]: array([ 0.33356456, 0.33355585, 0.33355585, 0.33401788, 0.33370132]) 

In [79]: lng 
Out[79]: array([ 1.27253229, 1.27249141, 1.27249141, 1.27259085, 1.2724337 ]) 

In [80]: original_app(lat,lng) 
Out[80]: 
[0.2522702110418014, 
0.2522702110418014, 
2.909533226553249, 
1.0542204712876762, 
0.0, 
3.003834632906676, 
0.9897592295963831, 
3.003834632906676, 
0.9897592295963831, 
2.2276138997714474] 

In [81]: vectorized_app1(lat,lng) 
Out[81]: 
array([[ 0.  , 0.25227021, 0.25227021, 2.90953323, 1.05422047], 
     [ 0.25227021, 0.  , 0.  , 3.00383463, 0.98975923], 
     [ 0.25227021, 0.  , 0.  , 3.00383463, 0.98975923], 
     [ 2.90953323, 3.00383463, 3.00383463, 0.  , 2.2276139 ], 
     [ 1.05422047, 0.98975923, 0.98975923, 2.2276139 , 0.  ]]) 

In [82]: vectorized_app2(lat,lng) 
Out[82]: 
array([ 0.25227021, 0.25227021, 2.90953323, 1.05422047, 0.  , 
     3.00383463, 0.98975923, 3.00383463, 0.98975923, 2.2276139 ]) 

Продолжительность -

In [83]: lat = np.random.randn(1000) 

In [84]: lng = np.random.randn(1000) 

In [85]: %timeit original_app(lat,lng) 
1 loops, best of 3: 2.11 s per loop 

In [86]: %timeit vectorized_app1(lat,lng) 
1 loops, best of 3: 263 ms per loop 

In [87]: %timeit vectorized_app2(lat,lng) 
1 loops, best of 3: 224 ms per loop 

Таким образом, для работы, похоже, vectorized_app2 может быть путем!

Answer 3

Поскольку в настоящее время это лучший результат Google для «парного дальнего расстояния», я добавлю свои два цента: эту проблему можно решить очень быстро, если у вас есть доступ к scikit-learn. Если посмотреть на sklearn.metrics.pairwise_distances, вы заметите, что метрика «haversine» не поддерживается, однако она реализована в sklearn.neighbors.DistanceMetric.

Это означает, что вы можете сделать следующее:

from sklearn.neighbors import DistanceMetric 

def sklearn_haversine(lat, lon): 
    haversine = DistanceMetric.get_metric('haversine') 
    latlon = np.hstack((lat[:, np.newaxis], lon[:, np.newaxis])) 
    dists = haversine.pairwise(latlon) 
    return 6371 * dists 

Обратите внимание, что конкатенация lat и lon необходимо только, так как они являются отдельными массивами. Если вы передадите их в виде комбинированного массива формы (n_samples, 2), вы можете позвонить в haversine.pairwise прямо на них. Кроме того, умножение на 6371 также необходимо, только если вам нужны расстояния в километрах. Например. этот шаг не был бы необходим, если бы вы хотели просто найти ближайшую пару точек.

Проверка:

In [87]: lat = np.array([ 0.33356456, 0.33355585, 0.33355585, 0.33401788, 0.33370132]) 

In [88]: lng = np.array([ 1.27253229, 1.27249141, 1.27249141, 1.27259085, 1.2724337 ]) 

In [89]: sklearn_haversine(lat, lng) 
Out[89]: 
array([[ 0.  , 0.25227021, 0.25227021, 2.90953323, 1.05422047], 
     [ 0.25227021, 0.  , 0.  , 3.00383463, 0.98975923], 
     [ 0.25227021, 0.  , 0.  , 3.00383463, 0.98975923], 
     [ 2.90953323, 3.00383463, 3.00383463, 0.  , 2.2276139 ], 
     [ 1.05422047, 0.98975923, 0.98975923, 2.2276139 , 0.  ]]) 

Производительность:

In [91]: lat = np.random.randn(1000) 

In [92]: lng = np.random.randn(1000) 

In [93]: %timeit original_app(lat,lng) 
1 loops, best of 3: 1.46 s per loop 

In [94]: %timeit vectorized_app1(lat,lng) 
10 loops, best of 3: 86.7 ms per loop 

In [95]: %timeit vectorized_app2(lat,lng) 
10 loops, best of 3: 75.7 ms per loop 

In [96]: %timeit sklearn_haversine(lat,lng) 
10 loops, best of 3: 76 ms per loop 

В заключение вы можете получить выход Divakar-х vectorized_app1 со скоростью vectorized_app2 в более короткий и простой код.