在Python中计算给定纬度和经度数据的距离矩阵的有效方法

我有纬度和经度的数据，我需要计算两个包含位置的数组之间的距离矩阵。我用过这个This https://stackoverflow.com/questions/120283/working-with-latitude-longitude-values-in-java/123305#123305获取给定纬度和经度的两个位置之间的距离。

这是我的代码的示例：

import numpy as np
import math

def get_distances(locs_1, locs_2):
    n_rows_1 = locs_1.shape[0]
    n_rows_2 = locs_2.shape[0]
    dists = np.empty((n_rows_1, n_rows_2))
    # The loops here are inefficient
    for i in xrange(n_rows_1):
        for j in xrange(n_rows_2):
            dists[i, j] = get_distance_from_lat_long(locs_1[i], locs_2[j])
    return dists


def get_distance_from_lat_long(loc_1, loc_2):

    earth_radius = 3958.75

    lat_dif = math.radians(loc_1[0] - loc_2[0])
    long_dif = math.radians(loc_1[1] - loc_2[1])
    sin_d_lat = math.sin(lat_dif / 2)
    sin_d_long = math.sin(long_dif / 2)
    step_1 = (sin_d_lat ** 2) + (sin_d_long ** 2) * math.cos(math.radians(loc_1[0])) * math.cos(math.radians(loc_2[0])) 
    step_2 = 2 * math.atan2(math.sqrt(step_1), math.sqrt(1-step_1))
    dist = step_2 * earth_radius

    return dist

我的预期输出是这样的：

>>> locations_1 = np.array([[34, -81], [32, -87], [35, -83]])
>>> locations_2 = np.array([[33, -84], [39, -81], [40, -88], [30, -80]])
>>> get_distances(locations_1, locations_2)
array([[ 186.13522573,  345.46610882,  566.23466349,  282.51056676],
       [ 187.96657622,  589.43369894,  555.55312473,  436.88855214],
       [ 149.5853537 ,  297.56950329,  440.81203371,  387.12153747]])

性能对我来说很重要，我能做的一件事就是使用Cython加快循环速度，但如果我不必去那里就好了。

有没有一个模块可以做这样的事情？或者还有其他解决方案吗？

您使用的半正矢方程中有很多次优的东西。您可以修剪其中的一部分，并最大限度地减少需要计算的正弦、余弦和平方根的数量。以下是我能想到的最好的，在我的系统上，在 1000 和 2000 个元素的两个随机数组上，它的运行速度比 Ophion 的代码（就矢量化而言基本相同）快 5 倍：

def spherical_dist(pos1, pos2, r=3958.75):
    pos1 = pos1 * np.pi / 180
    pos2 = pos2 * np.pi / 180
    cos_lat1 = np.cos(pos1[..., 0])
    cos_lat2 = np.cos(pos2[..., 0])
    cos_lat_d = np.cos(pos1[..., 0] - pos2[..., 0])
    cos_lon_d = np.cos(pos1[..., 1] - pos2[..., 1])
    return r * np.arccos(cos_lat_d - cos_lat1 * cos_lat2 * (1 - cos_lon_d))

如果你“按原样”向它提供两个数组，它会抱怨，但这不是一个错误，而是一个功能。基本上，该函数计算球体在最后一个维度上的距离，并在其余维度上广播。所以你可以获得你想要的东西：

>>> spherical_dist(locations_1[:, None], locations_2)
array([[ 186.13522573,  345.46610882,  566.23466349,  282.51056676],
       [ 187.96657622,  589.43369894,  555.55312473,  436.88855214],
       [ 149.5853537 ,  297.56950329,  440.81203371,  387.12153747]])

但它也可以用于计算两个点列表之间的距离，即：

>>> spherical_dist(locations_1, locations_2[:-1])
array([ 186.13522573,  589.43369894,  440.81203371])

或者在两个单点之间：

>>> spherical_dist(locations_1[0], locations_2[0])
186.1352257300577

这受到了 gufunc 工作方式的启发，一旦你习惯了它，我发现它是一种美妙的“瑞士军刀”编码风格，它可以让你在许多不同的设置中重用单个函数。