Быстрый способ пересечения полигонов с помощью Shapely
У меня есть большое количество полигонов (~100000), и я пытаюсь найти умный способ расчета их пересекающейся площади с ячейками регулярной сетки.
В настоящее время я создаю полигоны и ячейки сетки, используя библиотеку Shapely (на основе их координат углов). Затем, используя простой цикл for
, я перебираю каждый полигон и сравниваю его с соседними ячейками сетки.
Вот небольшой пример для иллюстрации полигонов и ячеек сетки:
from shapely.geometry import box, Polygon
# Пример полигона
xy = [[130.21001, 27.200001], [129.52, 27.34], [129.45, 27.1], [130.13, 26.950001]]
polygon_shape = Polygon(xy)
# Пример ячейки сетки
gridcell_shape = box(129.5, -27.0, 129.75, 27.25)
# Пересечение
polygon_shape.intersection(gridcell_shape).area
(Кстати, размеры ячеек сетки составляют 0.25x0.25, а максимальный размер полигонов — 1x1.)
На данный момент это довольно быстро для одной пары "полигон-ячейка сетки" — около 0.003 секунд. Однако, запуск этого кода на большом количестве полигонов (каждый из которых может пересекаться с десятками ячеек сетки) занимает около 15+ минут (иногда до 30+ минут в зависимости от числа пересекающихся ячеек сетки) на моем компьютере, что является неприемлемым. К сожалению, у меня нет представления о том, как можно написать код для нахождения пересечения полигонов, чтобы получить площадь перекрытия. Есть ли у вас советы? Есть ли альтернатива Shapely?
1 ответ(ов)
Вот еще одна версия ответа с более точным контролем над IoU (Intersection over Union):
def merge_intersecting_polygons(list_of_polygons, image_width, image_height):
"""Сливает пересекающиеся полигоны с использованием библиотеки shapely.
Слияние происходит только если Intersection over Union больше 0.5.
Ускорение достигается с помощью STRTree.
"""
# создаем полигоны shapely
shapely_polygons = []
for polygon in list_of_polygons:
shapely_polygons.append(Polygon(polygon))
# создаем STRTree
tree = STRtree(shapely_polygons)
# сливаем полигоны
merged_polygons = []
for i, polygon in enumerate(shapely_polygons):
# находим пересекающиеся полигоны
intersecting_polygons = tree.query(polygon)
# сливаем пересекающиеся полигоны
for intersecting_polygon in intersecting_polygons:
if polygon != intersecting_polygon:
# вычисляем intersection over union
intersection = polygon.intersection(intersecting_polygon).area
union = polygon.union(intersecting_polygon).area
iou = intersection / union
if iou > 0.5:
# сливаем полигоны
polygon = polygon.union(intersecting_polygon)
# добавляем слитый полигон в список
merged_polygons.append(polygon)
# удаляем дубликаты
merged_polygons = list(set(merged_polygons))
# конвертируем полигоны shapely в список полигонов
list_of_polygons = []
for polygon in merged_polygons:
list_of_polygons.append(np.array(polygon.exterior.coords).tolist())
return list_of_polygons
В этом коде реализована функция, которая сливает пересекающиеся полигоны, используя библиотеку shapely
. Обратите внимание, что слияние происходит только в том случае, если IoU больше 0.5, что позволяет получить более точные результаты. Ускорение процесса обеспечивается использованием структуры данных STRTree
, что позволяет быстро находить пересечения между полигонами.
Наиболее эффективный способ применения функции к массиву NumPy
Как извлечь частоту, связанную с FFT значениями в Python?
Цветовой график 2D массива в matplotlib
Преобразование байтового массива обратно в массив numpy
Взвешенный процентиль с помощью numpy