在Python中查找地理数据中圆内的所有坐标

Question

我有数百万个地理点。对于每一个点，我想找到所有的“邻近点”，也就是在一定半径内的其他点，比如几百米。

这个问题有一个简单的解决办法，时间复杂度是O(N^2)——就是计算所有点对之间的距离。但是，由于我处理的是地理距离，应该有更快的方法来解决这个问题。

我想在Python中实现这个功能。一个想到的解决方案是使用某种数据库（比如带有GIS扩展的MySQL，或者PostGIS），希望这样的数据库能通过某种索引高效地完成上述操作。不过，我更希望有一种更简单的方法，不需要我去搭建和学习这些技术。

几点说明

• 我将执行“查找邻居”的操作数百万次
• 数据将保持静态
• 因为这个问题在某种意义上比较简单，我希望能看到解决它的Python代码。

用Python代码来说，我想要类似这样的东西：

points = [(lat1, long1), (lat2, long2) ... ] # this list contains millions lat/long tuples
points_index = magical_indexer(points)
neighbors = []
for point in points:
    point_neighbors = points_index.get_points_within(point, 200) # get all points within 200 meters of point
    neighbors.append(point_neighbors) 

Answer 1

scipy

首先，得说说，已经有一些现成的算法可以用来做这类事情，比如 k-d 树。Scipy 提供了一个 Python 实现的 cKDtree，可以用来找到在特定范围内的所有点。

二分查找

不过，根据你要做的事情，实现这样的算法可能并不简单。而且，创建一棵树的过程相对复杂（可能会消耗不少资源），你也可以试试我之前用过的一个简单方法：

1. 计算数据集的主成分分析（PCA）。你想把数据集旋转，使得最重要的方向在前面，第二个方向（相对不那么重要的）在后面。你可以跳过这一步，直接选择 X 或 Y，但这样做计算上比较便宜，通常也容易实现。如果你选择 X 或 Y，建议选择方差更大的方向。
2. 按照主要方向（我们称之为 X）对点进行排序。
3. 要找到某个点的最近邻，先通过二分查找找到在 X 方向上最近的点的索引（如果这个点已经在你的集合中，你可能已经知道这个索引，不需要再查找）。然后逐个查看前后相邻的点，保持目前为止最好的匹配和它与搜索点的距离。当 X 方向的差值大于或等于目前最佳匹配的距离时，你可以停止查找（实际上，通常只需要查看很少的点）。
4. 要找到在给定范围内的所有点，做和第 3 步一样的事情，只是不要停止，直到 X 方向的差值超过范围。

实际上，你是在进行 O(N log(N)) 的预处理，对于每个点大约是 o(sqrt(N)) - 或者更多，如果你的点分布得不好。如果点大致均匀分布，X 方向上比最近邻更近的点数量大约是 N 的平方根。虽然如果很多点都在你的范围内，这种方法效率会降低，但通常也不会比暴力搜索差太多。

这个方法的一个优点是，它在内存分配上非常少，且大部分操作都能很好地利用内存局部性，这意味着尽管有明显的局限性，它的性能还是相当不错的。

Delauney 三角剖分

还有一个想法：可以使用 Delauney 三角剖分。在 Delauney 三角剖分中，任何点的最近邻都是相邻的节点。直观上说，在搜索过程中，你可以根据与查询点的绝对距离维护一个堆（优先队列）。选择最近的点，检查它是否在范围内，如果在，就添加它的所有邻居。我 怀疑 这样做不可能漏掉任何点，但你需要更仔细地研究一下才能确认……

Answer 2

在Eamon的提示下，我想出了一个简单的解决办法，使用了在SciPy中实现的b树。

from scipy.spatial import cKDTree
from scipy import inf

max_distance = 0.0001 # Assuming lats and longs are in decimal degrees, this corresponds to 11.1 meters
points = [(lat1, long1), (lat2, long2) ... ]
tree = cKDTree(points)

point_neighbors_list = [] # Put the neighbors of each point here

for point in points:
    distances, indices = tree.query(point, len(points), p=2, distance_upper_bound=max_distance)
    point_neighbors = []
    for index, distance in zip(indices, distances):
        if distance == inf:
            break
        point_neighbors.append(points[index])
    point_neighbors_list.append(point_neighbors)