在Pandas中计算地理密度的高效方法？

Question

我有一个包含美国快餐店经纬度数据的大列表。对于每个快餐店，我想知道在5英里范围内有多少其他快餐店。我可以使用Geopy和Pandas计算如下（DataFrame中的每一行表示不同的快餐店）：

import pandas as pd
import geopy.distance

df = pd.DataFrame({'Fast Food Place':[1,2,3], 'Lat':[33,34,35], 'Lon':[42,43,44]})

for index1, row1 in df.iterrows():
    num_fastfood = 0

    for index2, row2 in df.iterrows():
        # calculate distance in miles between longitude and latitude
        dist = geopy.distance.VincentyDistance(row1[['Lat','Lon']],
                                               row2[['Lat','Lon']]).miles

        # if fast food is within 5 miles, increment num_fastfood
        if dist < 5: # if less than five miles
            num_fastfood = num_fastfood + 1

    df.loc[index1, 'num_fastfood_5miles'] = num_fastfood - 1 # (subtract 1 to exclude self)

但是在非常大的数据集上（例如50,000行），这种方法非常慢。我考虑使用KDTree进行搜索，但想知道其他人是否有更快捷的方法？

- user1566200

---

KDTrees在这项任务中很难被击败。有没有不使用它的特别原因？ - Paul Brodersen

@Paul 没有特别的原因，只是好奇。我需要花一点时间来回忆如何使用sklearn的KDTree设置。类似于 tree = KDTree(my_lat_long) 然后循环遍历 nnDist 吗？ - user1566200

不要这样，使用 query_ball_tree 来获取半径内的所有点：tree = KDTree(my_lat_long); within_5 = tree.query_ball_tree(tree, radius=5)。然后展开嵌套列表并计数。 - Paul Brodersen

@Paul没有意识到这个存在，谢谢。如果我这样做，它只看经纬度，但我需要结合geopy来获取它的英里“真实”距离。 - user1566200

只需预先计算出5英里相当于多少度数，然后使用纬度/经度即可。 - Paul Brodersen

@Paul 不错的观点 - 大约是0.0724度。这对我来说足够准确了。我认为现在在SKLearn中是query_radius，我没有看到query_ball_tree。感谢您的帮助。如果您回答了我的问题，我会标记它为已解决。 - user1566200

Answer 1

使用scipy.spatial.cKDTree进行实现：

from scipy.spatial import cKDTree

def find_neighbours_within_radius(xy, radius):
    tree = cKDTree(xy)
    within_radius = tree.query_ball_tree(tree, r=radius)
    return within_radius

def flatten_nested_list(nested_list):
    return [item for sublist in nested_list for item in sublist]

def total_neighbours_within_radius(xy, radius):
    neighbours = find_neighbours_within_radius(xy, radius)
    return len(flatten_nested_list(neighbours))

- Paul Brodersen

---

我在使用时遇到了错误 TypeError: query_ball_tree() takes at least 2 positional arguments (1 given)。更新：算了，我发现我需要单独提供经纬度距离。 - user1566200

抱歉，cKDTree.query_ball_tree的关键字是r而不是radius。已修复代码。 - Paul Brodersen