python - 矢量化以计算许多距离

Question

我是 numpy/pandas 和矢量化计算的新手。我正在执行一项数据任务，我有两个数据集。数据集 1 包含具有经度和纬度的地点列表以及变量 A。数据集 2 还包含具有经度和纬度的地点列表。对于数据集 1 中的每个位置，我想计算它到数据集 2 中所有位置的距离，但我只想计算数据集 2 中小于变量 A 值的位置数。还要注意这两个数据集非常大，所以我需要使用矢量化操作来加快计算速度。

例如，我的数据集 1 可能如下所示:

id lon    lat   varA
1  20.11 19.88  100
2  20.87 18.65  90
3  18.99 20.75  120

我的数据集 2 可能如下所示:

placeid lon lat 
a       18.75 20.77
b       19.77 22.56
c       20.86 23.76
d       17.55 20.74 

然后对于数据集 1 中的 id == 1，我想计算它到数据集 2 中所有四个点(a、c、c、d)的距离，我想计算有多少距离更小比 varA 的相应值。例如计算出的四个距离分别为90、70、120、110，varA为100，则该值应为2。

我已经有一个向量化函数来计算两对坐标之间的距离。假设函数 (haversine(x,y)) 已正确实现，我有以下代码。

dataset2['count'] = dataset1.apply(lambda x: 
haversine(x['lon'],x['lat'],dataset2['lon'], dataset2['lat']).shape[0], axis 
= 1)

但是，这给出了总行数，而不是满足我要求的那些。

谁能告诉我如何使代码工作？

Answer 1

如果您可以将坐标投影到局部投影(例如 UTM )，这对于 pyproj 来说非常简单并且通常比 lon/lat 更适合测量，那么有很多使用 scipy.spatial 的方法MUCH 快得多。 df['something'] = df.apply(...) 和 np.vectorize() 都不是真正的矢量化，在幕后，它们使用循环。

ds1
    id  lon lat varA
0   1   20.11   19.88   100
1   2   20.87   18.65   90
2   3   18.99   20.75   120

ds2
    placeid lon lat
0   a   18.75   20.77
1   b   19.77   22.56
2   c   20.86   23.76
3   d   17.55   20.74


from scipy.spatial import distance

# gey coordinates of each set of points as numpy array
coords_a = ds1.values[:,(1,2)]
coords_b = ds2.values[:, (1,2)]
coords_a
#out: array([[ 20.11,  19.88],
#       [ 20.87,  18.65],
#       [ 18.99,  20.75]])

distances = distance.cdist(coords_a, coords_b)
#out: array([[ 1.62533074,  2.70148108,  3.95182236,  2.70059253],
#       [ 2.99813275,  4.06178532,  5.11000978,  3.92307278],
#       [ 0.24083189,  1.97091349,  3.54358575,  1.44003472]])

distances 实际上是每对点之间的距离。 coords_a.shape 是 (3, 2) 而 coords_b.shape 是 (4, 2)，所以结果是 (3,4)。 np.distance 的默认指标是 eculidean，但也有其他指标。为了这个例子，我们假设 vara 是:

vara = np.array([2,4.5,2])

(而不是 100 90 120)。我们需要确定第一行中 distances 中哪个值小于 2，第二行中哪个值小于 4.5，...，一种方法解决这个问题的方法是从相应的行中减去 vara 中的每个值(注意我们必须调整 vara 的大小):

vara.resize(3,1)
res = res - vara
#out: array([[-0.37466926,  0.70148108,  1.95182236,  0.70059253],
#       [-1.50186725, -0.43821468,  0.61000978, -0.57692722],
#       [-1.75916811, -0.02908651,  1.54358575, -0.55996528]])

然后将正值设置为零并将负值设置为正值将为我们提供最终数组:

res[res>0] = 0
res = np.absolute(res)
#out: array([[ 0.37466926,  0.        ,  0.        ,  0.        ],
#            [ 1.50186725,  0.43821468,  0.        ,  0.57692722],
#            [ 1.75916811,  0.02908651,  0.        ,  0.55996528]])

现在，对每一行求和:

sum_ = res.sum(axis=1)
#out:  array([ 0.37466926,  2.51700915,  2.34821989])

并计算每一行中的项目:

count = np.count_nonzero(res, axis=1)
#out: array([1, 3, 3])

这是一个完全矢量化(自定义)的解决方案，您可以根据自己的喜好进行调整，并且应该适应任何级别的复杂性。另一种解决方案是 cKDTree。代码来自文档。将它应用于您的问题应该相当容易，但如果您需要帮助，请随时提出。

x, y = np.mgrid[0:4, 0:4]
points = zip(x.ravel(), y.ravel())
tree = spatial.cKDTree(points)
tree.query_ball_point([2, 0], 1)
[4, 8, 9, 12]

query_ball_point() 找到点 x 距离 r 以内的所有点，速度快得惊人。

最后一点:不要将这些算法与经/纬度输入一起使用，尤其是当您感兴趣的区域远离赤道时，因为误差会变得很大。

更新:

要投影您的坐标，您需要将 WGS84 (lon/lat) 转换为适当的 UTM。要找出您应该计划使用哪个 utm 区域 epsg.io .

lon = -122.67598
lat = 45.52168
WGS84 = "+init=EPSG:4326"
EPSG3740 = "+init=EPSG:3740"
Proj_to_EPSG3740 = pyproj.Proj(EPSG3740)

Proj_to_EPSG3740(lon,lat)
# out: (525304.9265963673, 5040956.147893889)

您可以执行 df.apply() 并使用 Proj_to_... 来投影 df。