Obtenir efficacement l'intersection de plusieurs polygones en Python

12

Je voudrais obtenir l'intersection de plusieurs polygones. En utilisant le shapelypackage de Python , je peux trouver l'intersection de deux polygones en utilisant la intersectionfonction. Existe-t-il une fonction efficace similaire pour obtenir l'intersection de plusieurs polygones?

Voici un extrait de code pour comprendre ce que je veux dire:

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) 

Une intersection de deux cercles peut être trouvée par circle1.intersection(circle2). Je peux trouver l'intersection des trois cercles par circle1.intersection(circle2).intersection(circle3). Cependant, cette approche n'est pas vendable à un grand nombre de polygones car elle nécessite de plus en plus de code. Je voudrais une fonction qui prend un nombre arbitraire de polygones et renvoie leur intersection.

éclat
la source
je pense peut-être stocker les coordonnées dans un dictionnaire et les parcourir tout en utilisant des combinaisons d'importation à partir d'itertools. Je
posterai
Qu'entendez-vous par «leurs intersections»? Voulez-vous dire toutes les zones qui se croisent avec au moins un autre polygone, ou les zones que toutes les entrées se croisent?
jpmc26
Je veux dire l'intersection de tous les polygones, pas au moins un.
Splinter
Vous devriez clarifier cela ci-dessus (peut-être avec un exemple de sortie). Je suis presque certain que la plupart des réponses ne se comportent pas comme vous le souhaitez. (Et le fait que plusieurs répondeurs aient mal compris est une preuve suffisante que la question doit être clarifiée.)
jpmc26
1
@ jpmc26 Je viens d'ajouter une mise à jour à ma réponse où rtree est utilisé. L'approche est désormais plus efficace et évolutive. J'espère que cela t'aides!
Antonio Falciano

Réponses:

7

Une approche possible pourrait être de considérer la combinaison de paires de polygones, leurs intersections et enfin l'union de toutes les intersections via une union en cascade (comme suggéré ici ):

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

Une approche plus efficace devrait utiliser un index spatial, comme Rtree , afin de traiter un grand nombre de géométries (pas le cas des trois cercles):

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
la source
Je ne pense pas que cela fasse ce que le PO souhaite. Il rend les zones couvertes par au moins 2 polygones, tandis que l'OP ne recherche que les zones couvertes par tous les polygones de l'ensemble. Voir la clarification dans les commentaires.
jpmc26
3

Pourquoi ne pas utiliser une itération ou une récursivité? quelque chose comme :

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
la source
2

Essayez ce code. son concept est assez simple et je crois que vous obtenez ce que vous recherchez.

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

et si vous voulez que la sortie soit stockée sous forme de fichier de formes, utilisez 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)})

cela produit -

entrez la description de l'image ici

ziggy
la source
3
Je ne pense pas que cela fasse ce que le PO souhaite. Il rend les zones couvertes par au moins 2 polygones, tandis que l'OP ne recherche que les zones couvertes par tous les polygones de l'ensemble. Voir la clarification dans les commentaires. En outre, ket vsont de mauvais choix pour les noms de variables dans vos dictcompréhensions. Ces variables font chacune référence à différents éléments de dic.items(), et non à une paire clé-valeur. Quelque chose comme a, bserait moins trompeur.
jpmc26
1
ohh ok ouais je ne comprenais pas ce qu'il voulait dire
ziggy
et point bien pris au sujet de mes choix k, v - j'utilise automatiquement k, v lors de la boucle dans un dictionnaire ... ne lui donne pas beaucoup de réflexion
ziggy