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

Posted 2023-03-23

【中文标题】在python中查找地理数据中圆圈内的所有坐标

【英文标题】：Find all coordinates within a circle in geographic data in python

【发布时间】：2011-09-16 07:17:42

【问题描述】：

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

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

【问题讨论】：

您是否需要多次执行此操作（这样，做一些（艰苦的）工作可能会很有用，每次您需要某个点时进行更简单的计算）？您是否需要为多个点获取邻居，或者每次都是同一点的“中心”？

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

地理点在应用程序期间是静态的，还是每次执行查询时都不同？

我不明白为什么在 Postgres 中的点上使用 GIST 索引是不可取的。您肯定不想一直重新计算每一个最后一百万点的邻居吗？

点将保持不变。正如您所提到的，Postgres 中的 GIS 索引可能会解决问题，但我不知道如何使用 Postgres，并且更喜欢不需要我学习和构建其他技术的更简单的灵魂。

【参考方案1】：

scipy

​​>

首先要做的事情：有预先存在的算法来做一些事情，例如k-d tree。 Scipy 有一个 python 实现 cKDtree 可以找到给定范围内的所有点。

二分查找

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

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

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

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

德劳尼三角剖分

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

【讨论】：

我认为您对 k-d 树的建议正在寻找正确的解决方案。我现在正在查看 scipy.spatial.cKDTree ，它似乎是一个易于以直接、pythonic 方式使用的实现。

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

如果您包含一些使用该实现来解决问题的python代码，我会接受您的回答...

【参考方案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)

【讨论】：

嗨，@conradlee。你是如何计算出这个距离的？我的意思是，如果我喜欢使用 2km，例如，我将如何计算 max_distance 的值？谢谢。

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

但是，由于纬度和经度的度数->米的差异，以及不同纬度的变化，这个解决方案只是近似的。但是，似乎没有办法在 scipy 或 sklearn KDTree 实现中使用自定义距离函数（如 Haversine）。