在大数据集上找到离每个点最近的线，可能使用 shapely 和 rtree

Question

我有一个城市的简化地图，其中街道作为线串，地址作为点。我需要找到从每个点到任何街道线的最近路径。我有一个可以执行此操作的工作脚本，但它在多项式时间内运行，因为它嵌套了 for 循环。对于 150 000 行（shapely LineString）和 10 000 个点（shapely Point），在 8 GB Ram 计算机上需要 10 个小时才能完成。

函数看起来像这样（抱歉不能完全重现）：

import pandas as pd
import shapely
from shapely import Point, LineString

def connect_nodes_to_closest_edges(edges_df , nodes_df,
                                   edges_geom,
                                   nodes_geom):
    """Finds closest line to points and returns 2 dataframes:
        edges_df
        nodes_df
    """
    for i in range(len(nodes_df)):
        point = nodes_df.loc[i,nodes_geom]
        shortest_distance = 100000
        for j in range(len(edges_df)):
            line = edges_df.loc[j,edges_geom]
            if line.distance(point) < shortest_distance:
                shortest_distance = line.distance(point)
                closest_street_index = j
                closest_line = line
                ...

然后我将结果作为新列保存在表格中，将点到线的最短路径添加为新列。

有没有办法通过添加一些功能使其更快？

例如，如果我可以过滤掉 50m 左右以外的每个点的线，这将有助于加快每次迭代？

有没有办法使用 rtree 包来加快速度？我能够找到一个答案，使用于查找多边形交集的脚本更快，但我似乎无法使其适用于最近的点到线。

Faster way of polygon intersection with shapely

https://pypi.python.org/pypi/Rtree/

对不起，如果这个问题已经回答了，但我在这里和 gis.stackexchange 上都没有找到答案

谢谢你的建议！

Answer 1

这里有一个使用rtree 库的解决方案。这个想法是建立盒子 包含对角线段，并使用该框来构建 rtree。这将是最耗时的操作。 稍后，您使用以该点为中心的框查询 rtree。你得到几个 您需要检查的最低命中数，但命中数将为 （希望）数量级低于检查所有段。

在solutions 字典中，您将获得每个点的线ID、最近的线段、 最近的点（线段的一个点），以及到该点的距离。

代码中有一些 cmets 可以帮助您。 考虑到您可以序列化 rtree 以供以后使用。事实上，我会建议构建 rtree，保存它，然后使用它。因为常量MIN_SIZE 和INFTY 的调整可能会出现异常，并且您不想丢失构建rtree 所做的所有计算。

MIN_SIZE 太小意味着您可能会在解决方案中出现错误，因为如果该点周围的框不与线段相交，它可能会与不是最近线段的线段相交（很容易思考一个案例）。

一个太大的MIN_SIZE 意味着有太多的误报，在极端情况下会使代码尝试所有的段，你将处于与以前相同的位置，或者最糟糕的是，因为你现在正在构建一个你并不真正使用的 rtree。

如果数据是来自城市的真实数据，我想您知道任何地址都会与距离小于几个街区的路段相交。这将使搜索实际上是对数的。

还有一条评论。我假设没有太大的段。由于我们使用线段作为 rtree 中框的对角线，如果一行中有一些大段，这将意味着一个巨大的框将分配给该段，并且所有地址框都会与它相交。为避免这种情况，您始终可以通过添加更多中间点来人为地增加 LineStrins 的分辨率。

import math
from rtree import index
from shapely.geometry import Polygon, LineString

INFTY = 1000000
MIN_SIZE = .8
# MIN_SIZE should be a vaule such that if you build a box centered in each 
# point with edges of size 2*MIN_SIZE, you know a priori that at least one 
# segment is intersected with the box. Otherwise, you could get an inexact 
# solution, there is an exception checking this, though.


def distance(a, b):
    return math.sqrt( (a[0]-b[0])**2 + (a[1]-b[1])**2 ) 

def get_distance(apoint, segment):
    a = apoint
    b, c = segment
    # t = <a-b, c-b>/|c-b|**2
    # because p(a) = t*(c-b)+b is the ortogonal projection of vector a 
    # over the rectline that includes the points b and c. 
    t = (a[0]-b[0])*(c[0]-b[0]) + (a[1]-b[1])*(c[1]-b[1])
    t = t / ( (c[0]-b[0])**2 + (c[1]-b[1])**2 )
    # Only if t 0 <= t <= 1 the projection is in the interior of 
    # segment b-c, and it is the point that minimize the distance 
    # (by pitagoras theorem).
    if 0 < t < 1:
        pcoords = (t*(c[0]-b[0])+b[0], t*(c[1]-b[1])+b[1])
        dmin = distance(a, pcoords)
        return pcoords, dmin
    elif t <= 0:
        return b, distance(a, b)
    elif 1 <= t:
        return c, distance(a, c)

def get_rtree(lines):
    def generate_items():
        sindx = 0
        for lid, l in lines:
            for i in xrange(len(l)-1):
                a, b = l[i]
                c, d = l[i+1]
                segment = ((a,b), (c,d))
                box = (min(a, c), min(b,d), max(a, c), max(b,d)) 
                #box = left, bottom, right, top
                yield (sindx, box, (lid, segment))
                sindx += 1
    return index.Index(generate_items())

def get_solution(idx, points):
    result = {}
    for p in points:
        pbox = (p[0]-MIN_SIZE, p[1]-MIN_SIZE, p[0]+MIN_SIZE, p[1]+MIN_SIZE)
        hits = idx.intersection(pbox, objects='raw')    
        d = INFTY
        s = None
        for h in hits: 
            nearest_p, new_d = get_distance(p, h[1])
            if d >= new_d:
                d = new_d
                s = (h[0], h[1], nearest_p, new_d)
        result[p] = s
        print s

        #some checking you could remove after you adjust the constants
        if s == None:
            raise Exception("It seems INFTY is not big enough.")

        pboxpol = ( (pbox[0], pbox[1]), (pbox[2], pbox[1]), 
                    (pbox[2], pbox[3]), (pbox[0], pbox[3]) )
        if not Polygon(pboxpol).intersects(LineString(s[1])):  
            msg = "It seems MIN_SIZE is not big enough. "
            msg += "You could get inexact solutions if remove this exception."
            raise Exception(msg)

    return result

我用这个例子测试了函数。

xcoords = [i*10.0/float(1000) for i in xrange(1000)]
l1 = [(x, math.sin(x)) for x in xcoords]
l2 = [(x, math.cos(x)) for x in xcoords]
points = [(i*10.0/float(50), 0.8) for i in xrange(50)]

lines = [('l1', l1), ('l2', l2)]

idx = get_rtree(lines)

solutions = get_solution(idx, points)

得到：

Answer 2

我一直在寻找解决方案，我找到了this，它使用了 Geopandas。 基本上，这是一种直接的方法，它考虑了点和线的边界框的重叠。 然而，由于空间索引，计算成本显着降低。