Effizientes Überschneiden mehrerer Polygone in Python

11

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.

Splitter
quelle
Ich denke, vielleicht speichern Sie die Koordinaten in einem Wörterbuch und durchlaufen sie, während Sie aus itertools Importkombinationen verwenden. Ich werde bald
posten
Was meinst du mit "ihren Schnittpunkten"? Meinen Sie alle Bereiche, die sich mit mindestens einem anderen Polygon schneiden, oder die Bereiche, die sich alle Eingaben schneiden?
jpmc26
Ich meine den Schnittpunkt aller Polygone, nicht mindestens eines.
Splitter
Sie sollten dies oben klarstellen (möglicherweise mit einer Beispielausgabe). Ich bin mir ziemlich sicher, dass sich die meisten Antworten nicht so verhalten, wie Sie es wünschen. (Und die Tatsache, dass mehrere Antwortende missverstanden haben, ist ein Beweis genug, dass die Frage geklärt werden muss.)
jpmc26
1
@ jpmc26 Ich habe gerade ein Update zu meiner Antwort hinzugefügt, in dem rtree verwendet wird. Der Ansatz ist jetzt effizienter und skalierbarer. Hoffe das hilft!
Antonio Falciano

Antworten:

6

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
Antonio Falciano
quelle
Ich glaube nicht, dass dies das tut, was das OP will. Es gibt die Bereiche zurück, die mindestens 2 Polygone abdecken, während das OP nur nach Bereichen sucht, die von allen Polygonen im Satz abgedeckt werden . Siehe Klarstellung in den Kommentaren.
jpmc26
3

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
julsbreakdown
quelle
2

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 -

Geben Sie hier die Bildbeschreibung ein

zickig
quelle
3
Ich glaube nicht, dass dies das tut, was das OP will. Es gibt die Bereiche zurück, die mindestens 2 Polygone abdecken, während das OP nur nach Bereichen sucht, die von allen Polygonen im Satz abgedeckt werden . Siehe Klarstellung in den Kommentaren. Darüber hinaus kund vsind schlechte Entscheidungen für Variablennamen in Ihrem dictVerständnis. Diese Variablen beziehen sich jeweils auf verschiedene Elemente von dic.items(), nicht auf ein Schlüssel-Wert-Paar. So etwas a, bwäre weniger irreführend.
jpmc26
1
ohh okay ja ich habe nicht verstanden was er meinte
zickig
und Punkt gut über meine k, v Entscheidungen getroffen - ich benutze nur automatisch k, v, wenn ich ein Wörterbuch
durchlaufe. Ich