在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) 

- conradlee

---

你需要多次执行此操作吗（因此进行一些（艰苦）工作一次，并且每次需要某个点时进行简单的计算可能很有用）？您是否需要为多个点获取邻居，或者“中心”每次都是相同的点？ - Martin Thurau

我想执行这个查询数百万次。事实上，我想找到每个点的邻居。 - conradlee

地理点在应用程序的使用期间是否是静态的，还是每次执行查询时都会有所不同？ - Patrick

我不明白为什么在Postgres中使用GIST索引来处理点数据是不可取的。毕竟，您肯定不想一直重新计算每个最后一百万个点的邻居吧？ - Denis de Bernardy

这些点将保持静态。正如您所提到的，Postgres中的GIS索引可能会解决问题，但我不熟悉Postgres，并且更喜欢一个更简单的解决方案，不需要我学习和构建其他技术。 - conradlee

Answer 1

scipy

首先要明确的是，已经有现成的算法可以做这种事情，比如k-d tree。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三角剖分

另一个想法：德劳内三角剖分可能可行。对于德劳内三角剖分，任何点的最近邻都是相邻的节点。直觉上，在搜索过程中，您可以基于距离查询点的绝对值维护一个堆（优先队列）。选择最近的点，检查它是否在范围内，如果是，则添加所有邻居。我怀疑这样做不可能错过任何点，但您需要仔细考虑一下以确保......

- Eamon Nerbonne

---

我认为你的k-d树建议是正确解决方案的正确轨迹。我现在正在查看scipy.spatial.cKDTree，它似乎是一种实现，可以以直接、Pythonic的方式简单使用。 - conradlee

是的，这看起来是一个很好的实现，它甚至提供了一个范围限制的查询参数！- 链接已添加到答案中。 - Eamon Nerbonne

如果您能够包含一些使用该实现来解决问题的Python代码，我将接受您的答案... - conradlee

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)

- conradlee

---

嗨，@conradlee。你是如何计算这个距离度量的呢？我的意思是，如果我想使用2公里，那么我该如何计算max_distance的值呢？谢谢。 - pceccon

经度/纬度的米制转换实际上取决于纬度，对于纬度和经度是不同的，因此这只是一个粗略的转换。然而，您可以为您的数据想出一个有用的技巧，例如：在北纬或南纬40°时，一度经度之间的距离为85公里。 - Aaron Bramson

然而，由于纬度和经度的度数到米的转换存在差异，并且在不同的纬度处存在变化，因此该解决方案仅是近似的。然而，在scipy或sklearn KDTree实现中似乎没有使用自定义距离函数（如Haversine）的方法。 - Aaron Bramson