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

У меня есть большое количество полигонов (~ 100000) и попытайтесь найти разумный способ вычисления их пересекающейся области с регулярными ячейками сетки.

В настоящее время я создаю многоугольники и ячейки сетки, используя красивые (основанные на их угловых координатах). Затем, используя простой цикл for, я просматриваю каждый многоугольник и сравниваю его с соседними ячейками сетки.

Небольшой пример, иллюстрирующий ячейки полигонов/сетки.

from shapely.geometry import box, Polygon
# Example polygon 
xy = [[130.21001, 27.200001], [129.52, 27.34], [129.45, 27.1], [130.13, 26.950001]]
polygon_shape = Polygon(xy)
# Example grid cell
gridcell_shape = box(129.5, -27.0, 129.75, 27.25)
# The intersection
polygon_shape.intersection(gridcell_shape).area

(BTW: ячейки сетки имеют размеры 0,25x0,25 и многоугольники 1x1 при макс.)

На самом деле это довольно быстро для отдельного комбинированного многоугольника/сетки со скоростью около 0,003 секунды. Однако запуск этого кода на огромном количестве полигонов (каждый из которых может пересекать десятки ячеек сетки) занимает около 15 минут (до 30 + мин в зависимости от количества пересекающихся ячеек сетки) на моей машине, что неприемлемо. К сожалению, я понятия не имею, как можно написать код для пересечения многоугольников, чтобы получить область перекрытия. Есть ли у вас какие-либо советы? Есть ли альтернатива красивой?

Ответ 1

Рассмотрите возможность использования Rtree, чтобы определить, какие ячейки сетки могут пересекаться многоугольником. Таким образом, вы можете удалить цикл for, используемый с массивом lat/lons, который, вероятно, является медленной частью.

Создайте свой код примерно так:

from shapely.ops import cascaded_union
from rtree import index
idx = index.Index()

# Populate R-tree index with bounds of grid cells
for pos, cell in enumerate(grid_cells):
    # assuming cell is a shapely object
    idx.insert(pos, cell.bounds)

# Loop through each Shapely polygon
for poly in polygons:
    # Merge cells that have overlapping bounding boxes
    merged_cells = cascaded_union([grid_cells[pos] for pos in idx.intersection(poly.bounds)])
    # Now do actual intersection
    print poly.intersection(merged_cells).area

Ответ 2

Похоже, что доступное руководство пользователя Shapely устарело, но с 2013/2014 года Shapely имеет strtree.py с классом STRtree. Я использовал его, и он работает хорошо.

Вот фрагмент из docstring:

STRtree - это R-дерево, которое создается с помощью алгоритма Sort-Tile-Recursive. STRtree принимает последовательность объектов геометрии в качестве параметра инициализации. После инициализации метод запроса может быть использован для создания пространственного запроса по этим объектам.

>>> from shapely.geometry import Polygon
>>> polys = [ Polygon(((0, 0), (1, 0), (1, 1))), Polygon(((0, 1), (0, 0), (1, 0))), Polygon(((100, 100), (101, 100), (101, 101))) ]
>>> s = STRtree(polys)
>>> query_geom = Polygon(((-1, -1), (2, 0), (2, 2), (-1, 2)))
>>> result = s.query(query_geom)
>>> polys[0] in result
True