半正弦波的最小的高效計算距離

Question

我有一個數據幀與> 2.7mm的座標，和一個單獨的列表的〜2000座標。我試圖返回座標每個單獨的行與列表中的每個座標之間的最小距離。以下代碼適用於小規模（200行數據框），但在計算2.7MM以上的行時，它似乎永遠運行。

from haversine import haversine 

df 
Latitude Longitude 
39.989 -89.980 
39.923 -89.901 
39.990 -89.987 
39.884 -89.943 
39.030 -89.931 

end_coords_list = [(41.342,-90.423),(40.349,-91.394),(38.928,-89.323)] 

for row in df.itertuples(): 
    def min_distance(row): 
     beg_coord = (row.Latitude, row.Longitude) 
     return min(haversine(beg_coord, end_coord) for end_coord in end_coords_list) 
    df['Min_Distance'] = df.apply(min_distance, axis=1) 

我知道問題在於計算的絕對數量正在發生（5.7MM * 2000 = 114億〜），而事實上，運行這個循環很多是非常低效的。

基於我的研究，它似乎是一個矢量化的NumPy函數可能是一個更好的方法，但我是Python和NumPy的新手，所以我不太清楚如何在這種特殊情況下實現這一點。

理想輸出：

df 
Latitude Longitude Min_Distance 
39.989 -89.980  3.7 
39.923 -89.901  4.1 
39.990 -89.987  4.2 
39.884 -89.943  5.9 
39.030 -89.931  3.1 

提前感謝！

Answer 1

在本質上haversine func是：

# convert all latitudes/longitudes from decimal degrees to radians 
lat1, lng1, lat2, lng2 = map(radians, (lat1, lng1, lat2, lng2)) 

# calculate haversine 
lat = lat2 - lat1 
lng = lng2 - lng1 

d = sin(lat * 0.5) ** 2 + cos(lat1) * cos(lat2) * sin(lng * 0.5) ** 2 
h = 2 * AVG_EARTH_RADIUS * asin(sqrt(d)) 

這裏有一個量化的方法，利用強大的NumPy broadcasting和NumPy ufuncs，以取代那些數學模塊funcs中，這樣我們就一氣呵成的整個數組操作 -

# Get array data; convert to radians to simulate 'map(radians,...)' part  
coords_arr = np.deg2rad(coords_list) 
a = np.deg2rad(df.values) 

# Get the differentiations 
lat = coords_arr[:,0] - a[:,0,None] 
lng = coords_arr[:,1] - a[:,1,None] 

# Compute the "cos(lat1) * cos(lat2) * sin(lng * 0.5) ** 2" part. 
# Add into "sin(lat * 0.5) ** 2" part. 
add0 = np.cos(a[:,0,None])*np.cos(coords_arr[:,0])* np.sin(lng * 0.5) ** 2 
d = np.sin(lat * 0.5) ** 2 + add0 

# Get h and assign into dataframe 
h = 2 * AVG_EARTH_RADIUS * np.arcsin(np.sqrt(d)) 
df['Min_Distance'] = h.min(1) 

爲了進一步提升性能，我們可以利用numexpr module來代替超凡的funcs。

---

運行時的測試和驗證

途徑 -

def loopy_app(df, coords_list): 
    for row in df.itertuples(): 
     df['Min_Distance1'] = df.apply(min_distance, axis=1) 

def vectorized_app(df, coords_list): 
    coords_arr = np.deg2rad(coords_list) 
    a = np.deg2rad(df.values) 

    lat = coords_arr[:,0] - a[:,0,None] 
    lng = coords_arr[:,1] - a[:,1,None] 

    add0 = np.cos(a[:,0,None])*np.cos(coords_arr[:,0])* np.sin(lng * 0.5) ** 2 
    d = np.sin(lat * 0.5) ** 2 + add0 

    h = 2 * AVG_EARTH_RADIUS * np.arcsin(np.sqrt(d)) 
    df['Min_Distance2'] = h.min(1) 

驗證 -

In [158]: df 
Out[158]: 
    Latitude Longitude 
0 39.989 -89.980 
1 39.923 -89.901 
2 39.990 -89.987 
3 39.884 -89.943 
4 39.030 -89.931 

In [159]: loopy_app(df, coords_list) 

In [160]: vectorized_app(df, coords_list) 

In [161]: df 
Out[161]: 
    Latitude Longitude Min_Distance1 Min_Distance2 
0 39.989 -89.980  126.637607  126.637607 
1 39.923 -89.901  121.266241  121.266241 
2 39.990 -89.987  126.037388  126.037388 
3 39.884 -89.943  118.901195  118.901195 
4 39.030 -89.931  53.765506  53.765506 

計時 -

In [163]: df 
Out[163]: 
    Latitude Longitude 
0 39.989 -89.980 
1 39.923 -89.901 
2 39.990 -89.987 
3 39.884 -89.943 
4 39.030 -89.931 

In [164]: %timeit loopy_app(df, coords_list) 
100 loops, best of 3: 2.41 ms per loop 

In [165]: %timeit vectorized_app(df, coords_list) 
10000 loops, best of 3: 96.8 µs per loop