17

我有数百万个地理点。对于其中的每一个，我想找到所有“相邻点”，即某个半径内的所有其他点，比如几百米。

这个问题有一个简单的 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)
4

2 回答 2

7

scipy

首先要做的事情：有预先存在的算法来做一些事情，比如kd 树。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 三角剖分，假设任何点的最近邻居都是相邻节点。直觉是，在搜索过程中，您可以根据与查询点的绝对距离来维护一个堆（优先队列）。选择最近的点，检查它是否在范围内，如果是，则添加它的所有邻居。我怀疑不可能错过这样的任何一点，但是您需要更仔细地查看以确保...

于 2011-06-16T11:54:41.460 回答
7

在 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)
于 2011-06-16T13:15:14.833 回答