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

con*_*lee 16 python gis distance geospatial geo

我有数百万个地理点.对于其中的每一个,我想找到所有"相邻点",即在某个半径内的所有其他点,比如说几百米.

对这个问题有一个天真的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) 
Run Code Online (Sandbox Code Playgroud)

con*_*lee 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)
Run Code Online (Sandbox Code Playgroud)

  • 嗨,@ corradlee.你怎么弄清楚这个距离的衡量标准?我的意思是,如果我喜欢使用2km,我怎么能算出max_distance的值呢?谢谢. (2认同)
  • 然而,由于纬度和经度的度数->米存在差异,并且不同纬度的变化,该解决方案只是近似的。然而,似乎没有办法在 scipy 或 sklearn KDTree 实现中使用自定义距离函数(如 Haversine)。 (2认同)

Eam*_*nne 6

SciPy的

首先要做的事情是:有预先存在的算法可以做某事,比如kd树.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三角剖分

另一个想法:Delauney三角测量可以工作.对于Delauney三角剖分,给出任何点的最近邻居都是相邻节点.直觉是在搜索过程中,您可以根据与查询点的绝对距离来维护堆(优先级队列).选择最近的点,检查它是否在范围内,如果是,则添加其所有邻居.我怀疑不可能错过这样的任何一点,但你需要更仔细地看一下才能确定......