高效计算Haversine距离的最小值

Question

高效计算Haversine距离的最小值

5

我有一个包含>2.7MM坐标的数据框(dataframe)，以及一个独立的包含约2,000个坐标的列表(list)。我想要返回每一行中坐标与列表中所有坐标之间的最小距离。以下代码可以在小规模情况下（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 * 2,000 = ~11.4BN），以及运行这么多循环非常低效。

根据我的研究，似乎使用向量化的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

感谢您提前的支持！

- Walt Reed

1

告诉我们关于这个“harversine”的内容。它接受什么输入？真正的矢量化通常需要将计算降低到基本数学问题，这些问题可以由编译代码中的numpy处理。我们无法对黑匣子进行矢量化处理。 - hpaulj

"haversine" 接受两个输入参数：起始坐标和结束坐标，并计算两者之间的距离（以公里为单位）。 - Walt Reed

2

我们需要有关该软件包源代码的信息。再次发布以确认这是否是链接 - https://github.com/mapado/haversine/blob/master/haversine/__init__.py？不要假设我们已经在我们的端口安装了所有软件包。 - Divakar

抱歉，你是正确的 - 那就是我正在使用的包。 - Walt Reed

快速浏览该软件包显示它使用math.sin等标量运算。但是，看起来很容易使用numpy.sin等函数重写计算。这将是一个很好的numpy编程练习，你可以自己尝试一下。 - hpaulj

显示剩余5条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Divakar · Accepted Answer

haversine函数的本质是：

# 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广播和NumPy通用函数的向量化方法，以替换那些数学模块函数，使我们能够一次性对整个数组进行操作。

# 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模块来替换超越函数。

运行时测试和验证

方法 -

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