0

Быстрый способ пересечения полигонов с помощью Shapely

22

У меня есть большое количество полигонов (~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 ответ(ов)

0

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

Чтобы ответить на вопрос, пожалуйста, войдите или зарегистрируйтесь