Python - 低效的空间距离计算（如何加速）
Python - inefficient spatial distance calculation (how can it be speed up)
我目前正在 Python 中尝试一些地理编码。过程如下：我有两个具有纬度和经度值的数据框（df1 和 df2，房屋和学校），并且想为 df1 中的每个观测值找到 df2 中的最近邻居。我使用以下代码：
from tqdm import tqdm
import numpy as np
import pandas as pd
import math
def distance(lat1, long1, lat2, long2):
R = 6371 # Earth Radius in Km
dLat = math.radians(lat2 - lat1) # Convert Degrees 2 Radians
dLong = math.radians(long2 - long1)
lat1 = math.radians(lat1)
lat2 = math.radians(lat2)
a = math.sin(dLat/2) * math.sin(dLat/2) + math.sin(dLong/2) * math.sin(dLong/2) * math.cos(lat1) * math.cos(lat2)
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
d = R * c
return d
dists = []
schools =[]
for index, row1 in tqdm(df1.iterrows()):
for index, row2 in df2.iterrows():
dists.append(distance(row1.lat, row1.lng, row2.Latitude, row2.Longitude))
schools.append(min(dists))
del dists [:]
df1["school"] = pd.Series(schools)
代码有效，但需要很长时间。使用 tqdm，我的平均速度为每秒 2 次 df1 迭代。作为比较，我在 STATA 中使用 geonear 完成了整个任务，并且 df1 (950) 中的所有观察都需要 1 秒。我在 geonear 的帮助文件中读到他们使用聚类，不计算所有距离，但只计算最近的距离。然而，在我添加一个聚类函数（它也可能需要 CPU 的能力）之前，我想知道是否有人看到了一些方法来加速这个过程（我是 python 的新手并且可能有一些低效的代码会减慢进程）。或者是否有一个包可以更快地完成这个过程？
如果比 STATA 花费的时间更长，我会没问题，但不会接近 7 分钟...
提前致谢
pandas 的效率关键是对整个 dataframe/series 执行操作，而不是逐行执行。那么，让我们这样做吧。
for index, row1 in tqdm(df1.iterrows()):
for index, row2 in df2.iterrows():
在这里计算两个数据帧的笛卡尔积。这可以像这样更快地完成：
df_product = pd.merge(df1.assign(key=0, index=df1.index),
df2.assign(key=0),
on='key').drop('key', axis=1)
（代码取自 here）。我还添加了一个索引为 df1 的列，我们稍后将需要它来计算每个实体距 df1 的 min 距离。
现在使用 numpy 以矢量化方式计算所有增量、弧度纬度、a、c 和距离：
dLat = np.radians(df_product['Latitude'] - df_product['lat'])
dLong = np.radians(df_product['Longitude'] - df_product['lng'])
lat1 = np.radians(df_product['lat'])
lat2 = np.radians(df_product['Latitude'])
a = (np.sin(dLat / 2) ** 2
+ (np.sin(dLong / 2) ** 2) * np.cos(lat1) * np.cos(lat2))
c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
df_product['d'] = R * c
现在，从 df_product 开始，我们只保留包含我们之前添加的索引的列和包含距离的列。我们按索引对距离进行分组，计算相应的最小值并将它们分配给 df1['schools']，就像您在代码中所做的那样。
df1['schools'] = df_product.loc[:, ['index', 'd']].groupby('index', axis=0).min()
就是这样。对于每个数据帧中的 1000 行，对我来说一切都不到一秒钟。
您执行此操作的方式很慢，因为您使用的是 O(n²) 算法：每一行都查看其他行。 ，虽然引入了矢量化，但并没有从根本上解决效率低下的问题。
我建议将您的数据点加载到 kd-tree 中：此数据结构提供了一种在多个维度上查找最近邻居的快速方法。这样一棵树的构造在 O(n log n) 中，查询需要 O(log n)，所以总时间在 O(n log n).
如果您的数据定位到可以用平面很好地近似的地理区域，请投影您的数据，然后在二维中执行查找。否则，如果您的数据分散在全球，请投射到 spherical cartesian coordinates 并在那里执行查找。
您可以如何执行此操作的示例如下所示：
#/usr/bin/env python3
import numpy as np
import scipy as sp
import scipy.spatial
Rearth = 6371
#Generate uniformly-distributed lon-lat points on a sphere
#See: http://mathworld.wolfram.com/SpherePointPicking.html
def GenerateUniformSpherical(num):
#Generate random variates
pts = np.random.uniform(low=0, high=1, size=(num,2))
#Convert to sphere space
pts[:,0] = 2*np.pi*pts[:,0] #0-360 degrees
pts[:,1] = np.arccos(2*pts[:,1]-1) #0-180 degrees
#Convert to degrees
pts = np.degrees(pts)
#Shift ranges to lon-lat
pts[:,0] -= 180
pts[:,1] -= 90
return pts
def ConvertToXYZ(lonlat):
theta = np.radians(lonlat[:,0])+np.pi
phi = np.radians(lonlat[:,1])+np.pi/2
x = Rearth*np.cos(theta)*np.sin(phi)
y = Rearth*np.sin(theta)*np.sin(phi)
z = Rearth*np.cos(phi)
return np.transpose(np.vstack((x,y,z)))
#For each entry in qpts, find the nearest point in the kdtree
def GetNearestNeighbours(qpts,kdtree):
pts3d = ConvertToXYZ(qpts)
#See: https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.KDTree.query.html#scipy.spatial.KDTree.query
#p=2 implies Euclidean distance, eps=0 implies no approximation (slower)
return kdtree.query(pts3d,p=2,eps=0)
#Generate uniformly-distributed test points on a sphere. Note that you'll want
#to find a way to extract your pandas columns into an array of width=2, height=N
#to match this format.
df1 = GenerateUniformSpherical(10000)
df2 = GenerateUniformSpherical(10000)
#Convert df2 into XYZ coordinates. WARNING! Do not alter df2_3d or kdtree will
#malfunction!
df2_3d = ConvertToXYZ(df2)
#Build a kd-tree from df2_3D
kdtree = sp.spatial.KDTree(df2_3d, leafsize=10) #Stick points in kd-tree for fast look-up
#Return the distance to, and index of, each of df1's nearest neighbour points
distance, indices = GetNearestNeighbours(df1,kdtree)
我目前正在 Python 中尝试一些地理编码。过程如下：我有两个具有纬度和经度值的数据框（df1 和 df2，房屋和学校），并且想为 df1 中的每个观测值找到 df2 中的最近邻居。我使用以下代码：
from tqdm import tqdm
import numpy as np
import pandas as pd
import math
def distance(lat1, long1, lat2, long2):
R = 6371 # Earth Radius in Km
dLat = math.radians(lat2 - lat1) # Convert Degrees 2 Radians
dLong = math.radians(long2 - long1)
lat1 = math.radians(lat1)
lat2 = math.radians(lat2)
a = math.sin(dLat/2) * math.sin(dLat/2) + math.sin(dLong/2) * math.sin(dLong/2) * math.cos(lat1) * math.cos(lat2)
c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
d = R * c
return d
dists = []
schools =[]
for index, row1 in tqdm(df1.iterrows()):
for index, row2 in df2.iterrows():
dists.append(distance(row1.lat, row1.lng, row2.Latitude, row2.Longitude))
schools.append(min(dists))
del dists [:]
df1["school"] = pd.Series(schools)
代码有效，但需要很长时间。使用 tqdm，我的平均速度为每秒 2 次 df1 迭代。作为比较，我在 STATA 中使用 geonear 完成了整个任务，并且 df1 (950) 中的所有观察都需要 1 秒。我在 geonear 的帮助文件中读到他们使用聚类，不计算所有距离，但只计算最近的距离。然而，在我添加一个聚类函数（它也可能需要 CPU 的能力）之前，我想知道是否有人看到了一些方法来加速这个过程（我是 python 的新手并且可能有一些低效的代码会减慢进程）。或者是否有一个包可以更快地完成这个过程？
如果比 STATA 花费的时间更长，我会没问题，但不会接近 7 分钟...
提前致谢
pandas 的效率关键是对整个 dataframe/series 执行操作，而不是逐行执行。那么，让我们这样做吧。
for index, row1 in tqdm(df1.iterrows()): for index, row2 in df2.iterrows():
在这里计算两个数据帧的笛卡尔积。这可以像这样更快地完成：
df_product = pd.merge(df1.assign(key=0, index=df1.index),
df2.assign(key=0),
on='key').drop('key', axis=1)
（代码取自 here）。我还添加了一个索引为 df1 的列，我们稍后将需要它来计算每个实体距 df1 的 min 距离。
现在使用 numpy 以矢量化方式计算所有增量、弧度纬度、a、c 和距离：
dLat = np.radians(df_product['Latitude'] - df_product['lat'])
dLong = np.radians(df_product['Longitude'] - df_product['lng'])
lat1 = np.radians(df_product['lat'])
lat2 = np.radians(df_product['Latitude'])
a = (np.sin(dLat / 2) ** 2
+ (np.sin(dLong / 2) ** 2) * np.cos(lat1) * np.cos(lat2))
c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))
df_product['d'] = R * c
现在，从 df_product 开始，我们只保留包含我们之前添加的索引的列和包含距离的列。我们按索引对距离进行分组，计算相应的最小值并将它们分配给 df1['schools']，就像您在代码中所做的那样。
df1['schools'] = df_product.loc[:, ['index', 'd']].groupby('index', axis=0).min()
就是这样。对于每个数据帧中的 1000 行，对我来说一切都不到一秒钟。
您执行此操作的方式很慢，因为您使用的是 O(n²) 算法：每一行都查看其他行。
我建议将您的数据点加载到 kd-tree 中：此数据结构提供了一种在多个维度上查找最近邻居的快速方法。这样一棵树的构造在 O(n log n) 中，查询需要 O(log n)，所以总时间在 O(n log n).
如果您的数据定位到可以用平面很好地近似的地理区域，请投影您的数据，然后在二维中执行查找。否则，如果您的数据分散在全球，请投射到 spherical cartesian coordinates 并在那里执行查找。
您可以如何执行此操作的示例如下所示：
#/usr/bin/env python3
import numpy as np
import scipy as sp
import scipy.spatial
Rearth = 6371
#Generate uniformly-distributed lon-lat points on a sphere
#See: http://mathworld.wolfram.com/SpherePointPicking.html
def GenerateUniformSpherical(num):
#Generate random variates
pts = np.random.uniform(low=0, high=1, size=(num,2))
#Convert to sphere space
pts[:,0] = 2*np.pi*pts[:,0] #0-360 degrees
pts[:,1] = np.arccos(2*pts[:,1]-1) #0-180 degrees
#Convert to degrees
pts = np.degrees(pts)
#Shift ranges to lon-lat
pts[:,0] -= 180
pts[:,1] -= 90
return pts
def ConvertToXYZ(lonlat):
theta = np.radians(lonlat[:,0])+np.pi
phi = np.radians(lonlat[:,1])+np.pi/2
x = Rearth*np.cos(theta)*np.sin(phi)
y = Rearth*np.sin(theta)*np.sin(phi)
z = Rearth*np.cos(phi)
return np.transpose(np.vstack((x,y,z)))
#For each entry in qpts, find the nearest point in the kdtree
def GetNearestNeighbours(qpts,kdtree):
pts3d = ConvertToXYZ(qpts)
#See: https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.KDTree.query.html#scipy.spatial.KDTree.query
#p=2 implies Euclidean distance, eps=0 implies no approximation (slower)
return kdtree.query(pts3d,p=2,eps=0)
#Generate uniformly-distributed test points on a sphere. Note that you'll want
#to find a way to extract your pandas columns into an array of width=2, height=N
#to match this format.
df1 = GenerateUniformSpherical(10000)
df2 = GenerateUniformSpherical(10000)
#Convert df2 into XYZ coordinates. WARNING! Do not alter df2_3d or kdtree will
#malfunction!
df2_3d = ConvertToXYZ(df2)
#Build a kd-tree from df2_3D
kdtree = sp.spatial.KDTree(df2_3d, leafsize=10) #Stick points in kd-tree for fast look-up
#Return the distance to, and index of, each of df1's nearest neighbour points
distance, indices = GetNearestNeighbours(df1,kdtree)