向量化计算多个距离

Question

向量化计算多个距离

6

我对numpy/pandas和向量化计算都很陌生。我正在处理一个数据任务，其中有两个数据集。数据集1包含一系列带有经度、纬度和变量A的地点列表。数据集2也包含一系列带有经度和纬度的地点列表。对于数据集1中的每个地点，我想计算它与数据集2中所有地点的距离，但我只想得到小于变量A值的数据集2中地点的数量。请注意，这两个数据集都非常大，因此我需要使用向量化操作来加快计算速度。

例如，我的数据集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)

然而，这只给出了总行数，但并没有满足我的要求的行数。

有人能指导我如何使代码运行吗？

- macintosh81

3

是否指在向量化 Haversine 计算方面进行翻译？请查看您发布的右侧相关链接，例如 https://dev59.com/vFsW5IYBdhLWcg3w4aX7?rq=1。 - NaN

你能发布一下你正在使用的 Haversine 函数吗？ - DJK

@macintosh81 如果我的回答有用，请考虑接受/点赞它。 - Alz

3个回答

1

理解我所说的：

源数据框：

In [160]: d1
Out[160]:
   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

In [161]: d2
Out[161]:
  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

矢量化的 haversine 函数：

def haversine(lat1, lon1, lat2, lon2, to_radians=True, earth_radius=6371):
    if to_radians:
        lat1, lon1, lat2, lon2 = pd.np.radians([lat1, lon1, lat2, lon2])

    a = pd.np.sin((lat2-lat1)/2.0)**2 + \
        pd.np.cos(lat1) * pd.np.cos(lat2) * pd.np.sin((lon2-lon1)/2.0)**2

    return earth_radius * 2 * pd.np.arcsin(np.sqrt(a))

解决方案：

x = d2.assign(x=1) \
      .merge(d1.loc[d1['id']==1, ['lat','lon']].assign(x=1),
             on='x', suffixes=['','2']) \
      .drop(['x'], 1)

x['dist']  = haversine(x.lat, x.lon, x.lat2, x.lon2)

产生：

In [163]: x
Out[163]:
  placeid    lon    lat   lat2   lon2        dist
0       a  18.75  20.77  19.88  20.11  172.924852
1       b  19.77  22.56  19.88  20.11  300.078600
2       c  20.86  23.76  19.88  20.11  438.324033
3       d  17.55  20.74  19.88  20.11  283.565975

筛选：

In [164]: x.loc[x.dist < d1.loc[d1['id']==1, 'varA'].iat[0]]
Out[164]:
Empty DataFrame
Columns: [placeid, lon, lat, lat2, lon2, dist]
Index: []

让我们更改d1，这样几行就能满足条件:

In [171]: d1.loc[0, 'varA'] = 350

In [172]: d1
Out[172]:
   id    lon    lat  varA
0   1  20.11  19.88   350   # changed: 100 --> 350 
1   2  20.87  18.65    90
2   3  18.99  20.75   120

In [173]: x.loc[x.dist < d1.loc[d1['id']==1, 'varA'].iat[0]]
Out[173]:
  placeid    lon    lat   lat2   lon2        dist
0       a  18.75  20.77  19.88  20.11  172.924852
1       b  19.77  22.56  19.88  20.11  300.078600
3       d  17.55  20.74  19.88  20.11  283.565975

- MaxU - stand with Ukraine

嘿，感谢您的输入。这确实很有帮助！然而，这段代码似乎没有为数据集1中的每个id创建计数？此外，这对于大型数据集是否能够快速工作呢？ - macintosh81

1

使用 scipy.spatial.distance.cdist 函数，并将您自定义的距离算法作为 metric 参数传入。

h = lambda u, v: haversine(u['lon'], u['lat'], v['lon'], v['lat'])
dist_mtx = scipy.spatial.distance.cdist(dataset1, dataset2, metric = h)

然后要检查该地区的号码，只需广播它。

dataset2['count'] = np.sum(dataset1['A'][:, None] > dist_mtx, axis = 0)

- Daniel F

嗨，感谢您的输入！在使用空间距离方法时，我遇到了IndexError：只有整数、切片（:）、省略号（...）、numpy.newaxis（None）和整数或布尔数组是有效索引。请问我做错了什么吗？ - macintosh81

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

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

- Alz · Accepted Answer

如果您可以将坐标投影到本地投影（例如UTM），这在使用pyproj时非常简单，通常比经度/纬度更有利于测量，那么使用scipy.spatial会有一个更加快速的方法。无论是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)。默认的度量标准是 eculidean，但也有其他的度量标准。为了这个例子，假设 vara 是：

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

我们需要确定第一行中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（经度/纬度）转换为适当的 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。