Эффективное пересечение нескольких полигонов в Python

12

Я хотел бы получить пересечение нескольких полигонов. Используя shapelyпакет Python , я могу найти пересечение двух полигонов, используя intersectionфункцию. Существует ли подобная эффективная функция для получения пересечения нескольких многоугольников?

Вот фрагмент кода, чтобы понять, что я имею в виду:

from shapely.geometry import Point

coord1 = ( 0,0 )
point1 = Point(coord1)
circle1 = point1.buffer(1)

coord2 = ( 1,1 )
point2 = Point(coord2)
circle2 = point2.buffer(1)

coord3 = ( 1,0 )
point3 = Point(coord3)
circle3 = point3.buffer(1) 

Пересечение двух кругов можно найти с помощью circle1.intersection(circle2). Я могу найти пересечение всех трех кругов circle1.intersection(circle2).intersection(circle3). Однако этот подход не подходит для большого числа полигонов, так как требует все больше кода. Я хотел бы функцию, которая принимает произвольное количество многоугольников и возвращает их пересечение.

осколок
источник
Я думаю, может быть, сохранить координаты в словаре и циклически проходить по нему при использовании из itertools импорта комбинаций. Я скоро
отправлю
Что вы подразумеваете под "их пересечениями"? Вы имеете в виду все области, которые пересекаются хотя бы с одним другим многоугольником, или области, которые пересекают все входные данные?
jpmc26
Я имею в виду пересечение всех многоугольников, а не хотя бы одного.
осколок
Вы должны пояснить это выше (возможно, с примером вывода). Я почти уверен, что большинство ответов ведут себя не так, как вы хотите. (И тот факт, что несколько опрошенных неправильно поняли, является достаточным доказательством того, что вопрос нуждается в уточнении.)
jpmc26
1
@ jpmc26 Я только что добавил обновление в свой ответ, где используется rtree. Подход теперь более эффективен и масштабируем. Надеюсь это поможет!
Антонио Фальчано

Ответы:

7

Одним из возможных подходов может быть рассмотрение комбинации пар многоугольников, их пересечений и, наконец, объединение всех пересечений посредством каскадного объединения (как предложено здесь ):

from shapely.geometry import Point
from shapely.ops import cascaded_union
from itertools import combinations

circles = [
    Point(0,0).buffer(1),
    Point(1,0).buffer(1),
    Point(1,1).buffer(1),
]

intersection = cascaded_union(
    [a.intersection(b) for a, b in combinations(circles, 2)]
)
print intersection

Более эффективный подход должен использовать пространственный индекс, такой как Rtree , для работы с множеством геометрий (не в случае трех окружностей):

from shapely.geometry import Point
from shapely.ops import cascaded_union
from rtree import index

circles = [
    Point(0,0).buffer(1),
    Point(1,0).buffer(1),
    Point(1,1).buffer(1),
]
intersections = []
idx = index.Index()

for pos, circle in enumerate(circles):
    idx.insert(pos, circle.bounds)

for circle in circles:
    merged_circles = cascaded_union([circles[pos] for pos in idx.intersection(circle.bounds) if circles[pos] != circle])
    intersections.append(circle.intersection(merged_circles))

intersection = cascaded_union(intersections)
print intersection
Антонио Фальчано
источник
Я не верю, что это делает то, что хочет ОП. Он возвращает области, которые покрывают по крайней мере 2 полигона, тогда как OP ищет только области, покрытые всеми полигонами в наборе. Смотри разъяснения в комментариях.
jpmc26
3

Почему бы не использовать итерацию или рекурсивность? что-то типа :

from shapely.geometry import Point

def intersection(circle1, circle2):
    return circle1.intersection(circle2)

coord1 = ( 0,0 )
point1 = Point(coord1)
circle1 = point1.buffer(1)

coord2 = ( 1,1 )
point2 = Point(coord2)    
circle2 = point2.buffer(1)


coord3 = ( 1,0 )
point3 = Point(coord3)
circle3 = point3.buffer(1)
circles = [circle1, circle2, circle3]
intersectionResult = None

for j, circle  in enumerate(circles[:-1]):

    #first loop is 0 & 1
    if j == 0:
        circleA = circle
        circleB = circles[j+1]
     #use the result if the intersection
    else:
        circleA = intersectionResult
        circleB = circles[j+1]
    intersectionResult = intersection(circleA, circleB)

result= intersectionResult
julsbreakdown
источник
2

Дайте этот код выстрел. это довольно просто в концепции, и я верю, что вы получите то, что ищете.

from shapely.geometry import Point
from itertools import combinations
dic ={}
dic['coord1']=Point(0,0).buffer(1)
dic['coord2']=Point(1,1).buffer(1)
dic['coord3']=Point(1,0).buffer(1)
inter = {k[0]+v[0]:k[1].intersection(v[1]) for k,v in combinations(dic.items(),2)}
print inter

и если вы хотите, чтобы вывод был сохранен как шейп-файл, используйте fiona:

from shapely.geometry import Point,mapping
import fiona
from itertools import combinations
schema = {'geometry': 'Polygon', 'properties': {'Place': 'str'}}
dic ={}
dic['coord1']=Point(0,0).buffer(1)
dic['coord2']=Point(1,1).buffer(1)
dic['coord3']=Point(1,0).buffer(1)
inter = {k[0]+v[0]:k[1].intersection(v[1]) for k,v in combinations(dic.items(),2)}
print inter
with fiona.open(r'C:\path\abid', "w", "ESRI Shapefile", schema) as output:
    for x,y in inter.items():
        output.write({'properties':{'Place':x},'geometry':mapping(y)})

это выводит -

введите описание изображения здесь

Зигги
источник
3
Я не верю, что это делает то, что хочет ОП. Он возвращает области, которые покрывают по крайней мере 2 полигона, тогда как OP ищет только области, покрытые всеми полигонами в наборе. Смотри разъяснения в комментариях. Кроме того, kи vплохой выбор для имен переменных в вашем dictпонимании. Каждая из этих переменных относится к разным элементам dic.items(), а не к паре ключ-значение. Нечто подобное a, bбудет менее обманчивым.
jpmc26
1
ооо да, да, я не понял, что он имел в виду
зигги
и укажи хорошо взятый мой выбор k, v - я просто автоматически использую k, v при циклическом просмотре словаря .. не задумывался об этом
зигги