在大数据集上找到最接近的每一点，可能使用的是匀称和rtree

Question

我有一个城市的简化地图，其中的街道以街道和地址为点。我需要找到从每个点到任何路线的最近路径。我有一个工作脚本来做到这一点，但它在多项式时间运行，因为它已经嵌套循环。对于15万行（匀称的LineString）和10000点（匀称的Point），需要10小时才能完成8 GB Ram计算机。

功能看起来像这样（对不起，没有把它完全可重复）：

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 
       ... 

然后，我在一个表将结果保存为一个新列从点增加的最短路径，以线作为新柱。

有没有办法让它更快一些，并增加了一些功能？

如果我可以例如筛选出50米左右的每个点的线条，这将有助于加快每次迭代？

有没有办法让这个更快使用rtree包？我能够找到一个答案，使脚本更快地找到多边形的交集，但我似乎无法使它在最接近点的位置工作。

Faster way of polygon intersection with shapely

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

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

感谢您的咨询！

Answer 1

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

在solutions字典中，您将得到每个点的线段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) 

，并得到：