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

Question

我有一张简化的城市地图，其中街道作为线串，地址作为点。我需要找到从每个点到任何街道线最近的路径。我有一个可以这样做的工作脚本，但它在多项式时间内运行，因为它嵌套了for循环。对于150000行（shapely LineString）和10000点（shapely Point），在8GB 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包来加快这个过程？我能够找到一个答案，使寻找多边形交点的脚本更快，但我似乎无法使它适用于最接近直线的点。

多边形与Shapely相交的更快方式

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

如果已经回答了，很抱歉，但我在这里或gis上都找不到答案。堆栈交换

谢谢你的建议！

Answer 1

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

Answer 2

这里有一个使用rtree库的解决方案。其思想是构建包含对角线段的框，并使用该框构建rtree。这将是最耗时的操作。稍后，使用以点为中心的框查询rtree。你会得到一些你需要检查的最小点击数，但是点击数（希望）会比检查所有片段的数量少几个magnitud。

在解决方案中，您将获得每个点的行id、最近的段、最近的点（段中的一个点）以及到该点的距离。

代码中有一些注释可以帮助您。考虑到您可以序列化rtree以供以后使用。事实上，我建议构建rtree，保存它，然后使用它。因为调整常量MIN_SIZE和INFTY的例外情况可能会增加，并且您不希望丢失构建rtree所做的所有计算。

太小的MIN_SIZE将意味着您可能会在解决方案中出现错误，因为如果点周围的框不与线段相交，它可能会与不是最近线段的线段框相交（很容易认为是这种情况）。

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

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

还有一条评论。我假设没有太大的部分。由于我们使用段作为rtree中方框的对角线，因此如果在一条直线中有一些较大的段，这将意味着将为该段分配一个较大的方框，并且所有地址方框都将与之相交。要避免这种情况，始终可以通过添加更多中间点来人为提高线条的分辨率。

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)