在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 tree。 Scipy 有一个 python 实现 cKDtree 可以找到给定范围内的所有点。

二分查找

但是，根据您正在做的事情，实施类似的事情可能并非易事。此外，创建一棵树是相当复杂的（可能会产生相当多的开销），并且您可以通过我以前使用过的简单 hack 来摆脱困境：

1. 计算数据集的 PCA。您想要旋转数据集，以使最重要的方向是第一个，而正交（较小）的第二个方向是第二个。你可以跳过这个，只选择 X 或 Y，但它的计算成本很低，而且通常很容易实现。如果您只选择 X 或 Y，请选择方差较大的方向。
2. 按主要方向对点进行排序（将此方向称为 X）。
3. 要查找给定点的最近邻点，请通过二分搜索查找 X 中最近点的索引（如果该点已在您的集合中，您可能已经知道该索引并且不需要搜索）。迭代地查看下一个和上一个点，保持迄今为止的最佳匹配及其与搜索点的距离。当 X 的差值大于或等于迄今为止最佳匹配的距离时，您可以停止查找（实际上，通常只有很少的点）。
4. 要查找给定范围内的所有点，请执行与第 3 步相同的操作，但在 X 的差值超出范围之前不要停止。

实际上，您正在进行 O(N log(N)) 预处理，并且对于每个点大约 o(sqrt(N)) - 或更多，如果您的点分布不佳.如果这些点大致均匀分布，则 X 中比最近邻更近的点的数量将在 N 的平方根的数量级上。如果许多点在您的范围内，则效率会降低，但绝不会比蛮力差多少。

这种方法的一个优点是它可以在很少的内存分配中执行，并且大部分可以在非常好的内存局部性下完成，这意味着尽管存在明显的限制，但它的性能相当好。

德劳尼三角剖分

另一个想法：Delauney triangulation 可以工作。对于 Delauney 三角剖分，假设任何点的最近邻居都是相邻节点。直觉是，在搜索过程中，您可以根据与查询点的绝对距离来维护一个堆（优先队列）。选择最近的点，检查它是否在范围内，如果是，则添加它的所有邻居。我怀疑不可能错过任何这样的点，但您需要更仔细地查看以确保...

Answer 2

根据 Eamon 的提示，我提出了一个使用 SciPy 中实现的 btree 的简单解决方案。

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)