Python：为大量位置生成距离矩阵

Posted 2023-03-05

【中文标题】Python：为大量位置生成距离矩阵

【英文标题】：Python: generate distance matrix for large number of locations

【发布时间】：2019-10-09 15:10:21

【问题描述】：

我想根据 500 个位置的纬度和经度，使用 Haversine 公式生成一个距离矩阵 500X500。

这是 10 个位置的示例数据“coordinate.csv”：

Name,Latitude,Longitude
depot1,35.492807,139.6681689
depot2,33.6625572,130.4096027
depot3,35.6159881,139.7805445
customer1,35.622632,139.732631
customer2,35.857287,139.821461
customer3,35.955313,139.615387
customer4,35.16073,136.926239
customer5,36.118163,139.509548
customer6,35.937351,139.909783
customer7,35.949508,139.676462

得到距离矩阵后，我想根据距离矩阵找到离每个客户最近的仓库，然后将输出（从每个客户到壁橱仓库的距离和最近仓库的名称）保存到 Pandas DataFrame。

预期输出：

// Distance matrix
[ [..],[..],[..],[..],[..],[..],[..],[..],[..],[..] ]

// Closet depot to each customer (just an example)
Name,Latitude,Longitude,Distance_to_closest_depot,Closest_depot
depot1,35.492807,139.6681689,,
depot2,33.6625572,130.4096027,,
depot3,35.6159881,139.7805445,,
customer1,35.622632,139.732631,10,depot1
customer2,35.857287,139.821461,20,depot3
customer3,35.955313,139.615387,15,depot2
customer4,35.16073,136.926239,12,depot3
customer5,36.118163,139.509548,25,depot1
customer6,35.937351,139.909783,22,depot2
customer7,35.949508,139.676462,15,depot1

【问题讨论】：

所以你已经解释了你想要做什么。请注意，Stack Overflow 不是一个关于人们计划做什么的博客。但是，如果您有任何问题，请咨询How to ask a good question，因为目前没有问题。不要忘记发布您拥有的代码，指定您卡在哪里。

【参考方案1】：

有几个库函数可以帮助您：

来自scipy 的cdist 可用于使用您喜欢的任何距离度量生成距离矩阵。 还有一个haversine 函数可以传递给cdist。

之后，只需从距离矩阵中找到逐行最小值并将它们添加到您的 DataFrame 中。完整代码如下：

import pandas as pd
from scipy.spatial.distance import cdist
from haversine import haversine


df = pd.read_clipboard(sep=',')
df.set_index('Name', inplace=True)
customers = df[df.index.str.startswith('customer')]
depots = df[df.index.str.startswith('depot')]

dm = cdist(customers, depots, metric=haversine)
closest = dm.argmin(axis=1)
distances = dm.min(axis=1)

customers['Closest Depot'] = depots.index[closest]
customers['Distance'] = distances

结果：

            Latitude   Longitude Closest Depot    Distance
Name                                                      
customer1  35.622632  139.732631        depot3    4.393506
customer2  35.857287  139.821461        depot3   27.084212
customer3  35.955313  139.615387        depot3   40.565820
customer4  35.160730  136.926239        depot1  251.466152
customer5  36.118163  139.509548        depot3   60.945377
customer6  35.937351  139.909783        depot3   37.587862
customer7  35.949508  139.676462        depot3   38.255776

根据评论，我创建了一个替代解决方案，它使用平方距离矩阵。我认为原始解决方案更好，因为问题表明我们只想为每个客户找到最近的站点，因此无需计算客户之间和站点之间的距离。但是，如果您出于其他目的需要平方距离矩阵，请按照以下方式创建它：

import pandas as pd
import numpy as np
from scipy.spatial.distance import squareform, pdist
from haversine import haversine


df = pd.read_clipboard(sep=',')
df.set_index('Name', inplace=True)

dm = pd.DataFrame(squareform(pdist(df, metric=haversine)), index=df.index, columns=df.index)
np.fill_diagonal(dm.values, np.inf)  # Makes it easier to find minimums

customers = df[df.index.str.startswith('customer')]
depots = df[df.index.str.startswith('depot')]
customers['Closest Depot'] = dm.loc[depots.index, customers.index].idxmin()
customers['Distance'] = dm.loc[depots.index, customers.index].min()

最终的结果和之前一样，只是你现在有了一个平方距离矩阵。如果您愿意，可以在提取最小值后将 0 放回对角线上：

np.fill_diagonal(dm.values, 0)

【讨论】：

感谢您的回答。我可以使用 pd.read_csv("coordinate.csv") 代替 pd.read_clipboard 吗？

是的，pd.read_clipboard 只是我将您的数据读入 DataFrame 的方式，但 pd.read_csv("coordinate.csv") 应该适合您。

我希望 dm 是一个平方距离矩阵格式dm = [ [ ],[ ],...] 与 0 值对角线，因为从 A 到 A 的距离为 0。我该如何做到这一点并用你的其他代码来实现它同样的输出？

@belle - 我已经编辑了我的答案以提供第二种方法。

对不起，因为我的距离矩阵是来自谷歌地图距离矩阵api的数组（600x600），所以我想知道如何通过从该数组中读取数据来获取输出。你能再帮帮忙吗？非常感谢您的帮助。

【参考方案2】：

如果您需要一个非常大的矩阵并且可以使用带有 CUDA 的 NVIDIA GPU，您可以使用这个 numba 函数：

from numba import cuda
import math

@cuda.jit
def haversine_gpu_distance_matrix(p, G):
  i,j = cuda.grid(2)
  if i < p.shape[0] == G.shape[0] and j < p.shape[0] == G.shape[1]:
    if i == j:
      G[i][j] = 0
    else:
      longit_a = math.radians(p[i][0])
      latit_a = math.radians(p[i][1])
      longit_b = math.radians(p[j][0])
      latit_b =  math.radians(p[j][1])
      dist_longit_add = longit_b - longit_a
      dist_latit_sub = latit_b - latit_a
      dist_latit_add = latit_b + latit_a
      pre_comp = math.sin(dist_latit_sub/2)**2
      area = pre_comp + ((1 - pre_comp - math.sin(dist_latit_add/2)**2) * math.sin(dist_longit_add/2)**2)
      central_angle = 2 * math.asin(math.sqrt(area))
      radius = 3958
      G[i][j] = math.fabs(central_angle * radius)

您可以使用以下命令调用此函数：

# 10k [lon, lat] elements, replace this with your [lon, lat] array
# if your data is in a Pandas DataFrame, please convert it to a numpy array
geo_array = np.ones((10000, 2)) 
# allocate an empty distance matrix to fill when the function is called
dm_global_mem = cuda.device_array((geo_array.shape[0], geo_array.shape[0]))
# move the data in geo_array onto the GPU
geo_array_global_mem = cuda.to_device(geo_array)

# specify kernel dimensions, this can/should be further optimized for your hardware
threadsperblock = (16, 16)
blockspergrid_x = math.ceil(geo_array.shape[0] / threadsperblock[0])
blockspergrid_y = math.ceil(geo_array.shape[1] / threadsperblock[1])
blockspergrid = (blockspergrid_x, blockspergrid_y)

# run the function, which will transform dm_global_mem inplace
haversine_gpu_distance_matrix[blockspergrid, threadsperblock](geo_array_global_mem, dm_global_mem)

请注意，这可以针对您的硬件进一步优化。在 10k 个地理坐标对（即 100M 距离测量）上的 g4dn.xlarge 实例上的运行时间在编译后不到 0.01 秒。半径值设置为距离矩阵以英里为单位，您可以将其更改为 6371 以表示米。