Найдите ближайшую линию к каждой точке в большом наборе данных, возможно, используя shapely и rtree

У меня есть упрощенная карта города с улицами в виде линий и адресов в виде точек. Мне нужно найти ближайший путь от каждой точки до любой линии улицы. У меня есть рабочий сценарий, который делает это, но он выполняется за полиномиальное время, поскольку он вложен в цикл for. Для 150 000 линий (shapely LineString) и 10 000 точек (shapely Point) требуется 10 часов для обработки на компьютере RAM 8 ГБ.

Функция выглядит так (извините за то, что не сделала ее полностью воспроизводимой):

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? Мне удалось найти ответ, который ускоряет скрипт для поиска пересечения многоугольников, но я не могу заставить его работать для ближайшей точки к линии.

Более быстрый способ пересечения многоугольника с фигурным

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

Извините, если на этот вопрос уже был дан ответ, но я не нашел ответа ни здесь, ни на gis.stackexchange

Спасибо за совет!


person StefanK    schedule 12.09.2017    source источник
comment
Вы можете попробовать код в этом вопросе, чтобы узнать, насколько ускоряется фильтрация удаленных ссылок. Другой подход, который вы можете попробовать, - получить центральную точку линий линий (с помощью функции «shapely interpolate»), затем использовать rtree для поиска ссылок-кандидатов (поиск от точки к точке), а затем рассчитать расстояние.   -  person Alz    schedule 12.09.2017
comment
Просмотрите Как задать вопрос и минимальный воспроизводимый пример ...   -  person boardrider    schedule 13.09.2017
comment
@boardrider спасибо за ссылки, код редактировал. К сожалению, здесь сложно воспроизвести такой большой набор данных.   -  person StefanK    schedule 13.09.2017
comment
Вот где в игру вступает Minimal в минимальном воспроизводимом примере. Я посмотрел на код, и он все еще не похож на проверяемый пример   -  person boardrider    schedule 13.09.2017
comment
Вам все еще нужна помощь? Если так, нам понадобится некоторая информация о данных. 1) Сколько примерно точек у вас примерно в каждой LineString? 2) Это вычисление, которое вам нужно выполнить один раз, или вам нужно решение, которое работает для новых входящих точек и LineStrings? Я бы построил свою собственную структуру данных, подобную rtree, но вместо того, чтобы строить ее на (возможно, длинных) LineStrings, я бы использовал сегменты всех LineStrings. Вы также можете использовать библиотеку rtree напрямую, заплатив за некоторые дополнительные вычисления.   -  person eguaio    schedule 19.09.2017
comment
там есть несколько хороших советов! Что я делаю сейчас, так это разбиваю линии улиц на наименьшие сегменты всего по 2 точки, а затем запускаю вложенный цикл for с сохранением промежуточных результатов в списки, из которых я делаю фрейм данных в конце цикла. Несмотря на то, что мне нужно запустить его для выбранной области только один раз, как только область станет большой (какой-то большой город), расчет будет выполняться в течение 3 дней при текущих настройках. Не могли бы вы дать мне несколько советов по построению моего собственного r-дерева структура данных?   -  person StefanK    schedule 20.09.2017


Ответы (2)


Здесь у вас есть решение, использующее rtree библиотеку. Идея состоит в том, чтобы построить блоки, содержащие сегменты по диагонали, и использовать эти блоки для построения rtree. Это будет самая затратная по времени операция. Позже вы запрашиваете rtree с прямоугольником, центрированным в точке. Вы получаете несколько совпадений, которые вам нужно проверить по минимуму, но количество совпадений будет (надеюсь) на несколько порядков меньше, чем проверка по всем сегментам.

В solutions dict вы получите для каждой точки идентификатор линии, ближайший сегмент, ближайшую точку (точку сегмента) и расстояние до точки.

В коде есть несколько комментариев, которые могут вам помочь. Учтите, что вы можете сериализовать 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)

И получил:  Результат теста

person eguaio    schedule 20.09.2017
comment
Спасибо, что вложили в это столько труда! У меня будет время протестировать это решение позже на следующей неделе, поэтому я отмечу его как решение позже. - person StefanK; 23.09.2017
comment
Пожалуйста. Я нашел саму проблему очень интересной и потенциально пригодной для многих приложений, поэтому я хотел попробовать. Ждем результатов выступлений. - person eguaio; 23.09.2017
comment
Давай, друг мой, неизвестность убивает меня, хахаха. Знаете, мне нравятся такие тесты производительности, мне не терпелось попробовать новый код, пока вы это делаете. - person eguaio; 27.09.2017
comment
Извините за то, что не ответил раньше ... Я попытался реализовать это, но обнаружил, что это довольно сложно, так как мне пришлось бы многое изменить в моем исходном решении - написать его заново. Поэтому вместо этого я попытался отфильтровать точки, которые не находятся в квадрате 100 метров от центральной точки. Я создал столбец для точки x, для точки y и в квадратном столбце, который имеет значение True, если точка находится в квадрате 100 м. Этот шаг резко сократил общее время вычислений. Несмотря на то, что rtree будет быстрее, мое упрощенное решение можно было использовать на практике, поэтому мы остались с ним. Но все равно хорошая работа и спасибо за усилия - person StefanK; 01.12.2017
comment
К сожалению, я не получу сравнения производительности. Спасибо, что поделились своим решением и все равно ответили :) - person eguaio; 04.12.2017
comment
Кстати, большое спасибо за апвоуты, мне понятно, что это были вы :) - person eguaio; 04.12.2017

Я искал решение и нашел this, в котором используются геопанды. По сути, это простой подход, который учитывает перекрытие ограничивающих прямоугольников точек и линий. Однако благодаря пространственному индексу вычислительные затраты значительно снижаются.

person lenhhoxung    schedule 03.10.2020