Obtención de la intersección de múltiples polígonos de manera eficiente en Python

12

Me gustaría obtener la intersección de múltiples polígonos. Usando el shapelypaquete de Python , puedo encontrar la intersección de dos polígonos usando la intersectionfunción. ¿Existe una función eficiente similar para obtener la intersección de múltiples polígonos?

Aquí hay un fragmento de código para entender lo que quiero decir:

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) 

Se puede encontrar una intersección de dos círculos por circle1.intersection(circle2). Puedo encontrar la intersección de los tres círculos por circle1.intersection(circle2).intersection(circle3). Sin embargo, este enfoque no se puede vender a una gran cantidad de polígonos, ya que requiere cada vez más código. Me gustaría una función que tome un número arbitrario de polígonos y devuelva su intersección.

astilla
fuente
Estoy pensando que tal vez almacene las coordenadas en un diccionario y repítalo mientras usa combinaciones de importación de itertools. Voy a publicar pronto
ziggy
¿Qué quieres decir con "sus intersecciones"? ¿Te refieres a todas las áreas que se cruzan con al menos otro polígono, o las áreas que se cruzan todas las entradas?
jpmc26
Me refiero a la intersección de todos los polígonos, no al menos uno.
astilla
Deberías aclarar esto arriba (quizás con un resultado de ejemplo). Estoy bastante seguro de que la mayoría de las respuestas no se comportan como lo deseas. (Y el hecho de que varios respondedores hayan entendido mal es evidencia suficiente de que la pregunta necesita aclaración)
Jpmc26
1
@ jpmc26 Acabo de agregar una actualización a mi respuesta donde se usa rtree. El enfoque es más eficiente y escalable ahora. ¡Espero que esto ayude!
Antonio Falciano

Respuestas:

7

Un posible enfoque podría ser considerar la combinación de pares de polígonos, sus intersecciones y finalmente la unión de todas las intersecciones a través de una unión en cascada (como se sugiere aquí ):

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

Un enfoque más eficiente debería usar un índice espacial, como Rtree , para tratar con muchas geometrías (no el caso de los tres círculos):

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
fuente
No creo que esto haga lo que quiere el OP. Devuelve las áreas que cubren al menos 2 polígonos, mientras que el OP solo busca áreas cubiertas por todos los polígonos del conjunto. Ver aclaración en los comentarios.
jpmc26
3

¿Por qué no usar una iteración o recursividad? algo como :

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
fuente
2

Dale una oportunidad a este código. es bastante simple en concepto y creo que te da lo que estás buscando.

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

y si desea que la salida se almacene como un archivo shape, use 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)})

esto produce

ingrese la descripción de la imagen aquí

ziggy
fuente
3
No creo que esto haga lo que quiere el OP. Devuelve las áreas que cubren al menos 2 polígonos, mientras que el OP solo busca áreas cubiertas por todos los polígonos del conjunto. Ver aclaración en los comentarios. Además, ky vson malas elecciones para nombres de variables en sus dictcomprensiones. Esas variables se refieren a elementos diferentes dic.items(), no a un par clave-valor. Algo así a, bsería menos engañoso.
jpmc26
1
ohh bien, sí, no entendí lo que quería decir
ziggy
y apunte bien sobre mis opciones k, v: solo uso automáticamente k, v al recorrer un diccionario ... no pensé mucho en ello
ziggy