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

Question

我有一个简化的城市地图，其中的街道为线串，地址为点。我需要找到从每个点到任何一条街线的最近路径。我有一个执行此操作的脚本，但由于嵌套了循环，因此它在多项式时间内运行。对于15万行（形状为LineString）和10000点（形状为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包使其更快？我能够找到一个答案，从而使脚本可以更快地找到多边形的交点，但是我似乎无法使它适用于最接近点到线的地方。

多边形相交的更快方法

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

抱歉，如果已经回答了，但是我在这里也没有在gis.stackexchange上找到答案

谢谢你的建议！

Answer 1

在这里，您有使用rtree库的解决方案。这个想法是构建包含对角线段的框，并使用该框构建rtree。这将是最耗时的操作。稍后，您用一个以点为中心的框查询rtree。您会得到一些需要检查的最小命中数，但是命中数（希望）要比对所有分段进行查验都要少（希望）。

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

代码中有一些注释可以帮助您。考虑到您可以序列化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)

并得到：