Как я могу быстро найти многоугольник, содержащий точку, если все многоугольники являются разбиением? - PullRequest
Новая
       7

Как я могу быстро найти многоугольник, содержащий точку, если все многоугольники являются разбиением?

0 голосов
/ 09 ноября 2018

у меня есть функция 

get_polygon(polygon_collection, point):
    for polygon in polygon_collection:
        if polygon.intersects(point):
            return polygon
    return None

Этот подход работает, но он в O (n) * O (проверка одного полигона). Это наверняка можно уменьшить до O (log (n)) * O (проверка одного полигона), если построить структуру данных дерева.

Поддерживает ли формали это напрямую?

Структура коллекции полигонов

  • Коллекция содержит n MultiPolygons - следовательно, один объект может состоять из множества полигонов.
  • Каждый MultiPolygon может иметь анклавы / эксклавы
  • Обычно (> 95%) это простые многоугольники без дырок
  • Обычно (> 50%) форма довольно проста (например, квадраты)

Пример использования

Список полигонов может быть областями почтовых индексов в Германии. Это было бы несколько тысяч. Тогда у меня есть GPS-позиции от меня и некоторых друзей, это также будет несколько тысяч. Я хочу сказать, в какой области мы получили наибольшее количество точек данных.

1 Ответ

0 голосов
/ 20 ноября 2018

RT - это именно то, что вы ищете. Shapely теперь поддерживает создание и использование RTrees. Но дело в том, что стройные деревья поддерживают только прямоугольники. Итак, рецепт его использования может быть следующим:

  1. Постройте дерево с минимальными прямоугольниками, которые содержат многоугольники.
  2. Для заданной точки используйте RTree, чтобы узнать, какие прямоугольники включают эту точку. Точку нужно надуть немного, так как дерево ожидает, что поле проверяет пересечения. Это даст вам некоторые ложные срабатывания, но это нормально.
  3. Отфильтруйте ложные срабатывания, чтобы получить окончательный результат.

Если ваши полигоны не очень странные, это даст вам логарифмическое время в среднем (логарифмическое по количеству полигонов для каждой точки). Конечно, вам нужно заплатить цену за построение структуры RTree, но это также делается в логарифмическом времени, если полигоны так или иначе распределены одинаково, и rtree можно сериализовать для дальнейшего использования. (То есть добавление нового поля в дерево является логарифмическим по количеству ящиков, уже присутствующих в дереве, поэтому общая сложность меньше, чем n * log (n)).

Реализация этого заключается в следующем. 

from rtree import index
from shapely.geometry import Polygon, Point

def get_containing_box(p):
    xcoords = [x for x, _ in p.exterior.coords]
    ycoords = [y for _, y in p.exterior.coords]
    box = (min(xcoords), min(ycoords), max(xcoords), max(ycoords))   
    return box

def build_rtree(polys):
    def generate_items():
        pindx = 0
        for pol in polys:
            box = get_containing_box(pol)
            yield (pindx, box,  pol)
            pindx += 1
    return index.Index(generate_items())


def get_intersection_func(rtree_index):
    MIN_SIZE = 0.0001
    def intersection_func(point):
        # Inflate the point, since the RTree expects boxes:
        pbox = (point[0]-MIN_SIZE, point[1]-MIN_SIZE, 
                 point[0]+MIN_SIZE, point[1]+MIN_SIZE)
        hits = rtree_index.intersection(pbox, objects='raw')
        #Filter false positives:
        result = [pol for pol in hits if pol.intersects(Point(point)) ]
        return result
    return intersection_func

Пример, который я использовал для проверки: 

triangles = [
 [(8.774239095627754, 32.342041683826224),
  (8.750148703126552, 32.899346596303054),
  (4.919576457288677, 35.41040289384488)],
 [(8.485955148136222, 32.115258769399446),
  (9.263360720698277, 30.065319757618354),
  (4.562638192761559, 38.541192819415855)],
 [(2.613649959824923, 38.14802347408093),
  (7.879211442172622, 35.97223726358858),
  (0.9726266770834235, 32.12523430143625)]
]

polys = [Polygon(t) for t in triangles]
my_rtree = build_rtree(polys)
my_func = get_intersection_func(my_rtree)


my_intersections =  my_func((5, 35))
for pol in my_intersections:
    print pol.exterior.coords[:]
...