Быстрый способ пересечения полигонов с помощью 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