Effizientes Überschneiden mehrerer Polygone in Python

Ich möchte den Schnittpunkt mehrerer Polygone erhalten. Mit dem Python- shapelyPaket kann ich mithilfe der intersectionFunktion den Schnittpunkt zweier Polygone ermitteln . Gibt es eine ähnlich effiziente Funktion, um den Schnittpunkt mehrerer Polygone zu erhalten?

Hier ist ein Code-Snippet, um zu verstehen, was ich meine:

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) 

Ein Schnittpunkt zweier Kreise kann durch gefunden werden circle1.intersection(circle2). Ich kann den Schnittpunkt aller drei Kreise durch finden circle1.intersection(circle2).intersection(circle3). Dieser Ansatz ist jedoch für eine große Anzahl von Polygonen nicht verkaufbar, da er zunehmend mehr Code erfordert. Ich möchte eine Funktion, die eine beliebige Anzahl von Polygonen verwendet und deren Schnittpunkt zurückgibt.

Ein möglicher Ansatz könnte darin bestehen, die Kombination von Polygonpaaren, ihre Schnittpunkte und schließlich die Vereinigung aller Schnittpunkte über eine kaskadierte Vereinigung (wie hier vorgeschlagen ) zu berücksichtigen :

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

Ein effizienterer Ansatz sollte einen räumlichen Index wie Rtree verwenden , um mit vielen Geometrien umgehen zu können (nicht der Fall der drei Kreise):

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

Warum nicht eine Iteration oder Rekursivität verwenden? etwas wie :

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

Probieren Sie diesen Code aus. Das Konzept ist ziemlich einfach und ich glaube, Sie bekommen das, wonach Sie suchen.

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

und wenn Sie möchten, dass die Ausgabe als Shapefile gespeichert wird, verwenden Sie 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)})

diese Ausgaben -