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

Question

我有一个城市的简化 map ，其中有街道作为线串，地址作为点。我需要找到从每个点到任何街道线的最近路径。我有一个执行此操作的工作脚本，但它在多项式时间内运行，因为它嵌套了 for 循环。对于 150 000 条线(匀称的 LineString)和 10 000 个点(匀称的 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。你得到几个
您需要检查最小值的点击次数，但点击次数将为
(希望) magnitud 的顺序低于检查所有段。

在 solutions dict 你会得到，对于每个点，线 id，最近的线段，
最近的点(线段的一个点)，以及到该点的距离。

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

太小了MIN_SIZE将意味着您可能在解决方案中有错误，因为如果点周围的框不与线段相交，它可能与不是最近线段的线段的框相交(很容易想到一个案例)。

A 太大 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)

并得到: