Obtendo interseção de vários polígonos com eficiência em Python

12

Gostaria de obter a interseção de vários polígonos. Usando o shapelypacote do Python , posso encontrar a interseção de dois polígonos usando a intersectionfunção Existe uma função eficiente semelhante para obter a interseção de vários polígonos?

Aqui está um trecho de código para entender o que quero dizer:

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) 

Uma interseção de dois círculos pode ser encontrada por circle1.intersection(circle2). Eu posso encontrar a interseção dos três círculos por circle1.intersection(circle2).intersection(circle3). No entanto, essa abordagem não é vendável para um grande número de polígonos, pois requer cada vez mais código. Eu gostaria de uma função que pega um número arbitrário de polígonos e retorna sua interseção.

lasca
fonte
estou pensando que talvez armazene as cordas em um dicionário e passe por ele enquanto usa combinações de importação do itertools. Vou postar em breve
ziggy
O que você quer dizer com "suas interseções"? Você quer dizer todas as áreas que se cruzam com pelo menos um outro polígono ou as áreas que todas as entradas se cruzam?
perfil completo de jpmc26
Quero dizer a interseção de todos os polígonos, não pelo menos um.
Splinter
Você deve esclarecer isso acima (talvez com um exemplo de saída). Estou bastante certo de que a maioria das respostas não se comporta como você deseja. (Eo fato de vários answerers ter entendido mal é prova suficiente de que a questão precisa de esclarecimento.)
jpmc26
1
@ jpmc26 Acabei de adicionar uma atualização à minha resposta em que rtree é usado. A abordagem é mais eficiente e escalável agora. Espero que isto ajude!
Antonio Falciano

Respostas:

7

Uma abordagem possível seria considerar a combinação de pares de polígonos, suas interseções e, finalmente, a união de todas as interseções por meio de uma união em cascata (como sugerido aqui ):

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

Uma abordagem mais eficiente deve usar um índice espacial, como Rtree , para lidar com muitas geometrias (não o caso dos três 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
fonte
Não acredito que isso faça o que o OP quer. Ele devolve as áreas que pelo menos 2 polígonos cobrem, enquanto o OP está procurando apenas áreas cobertas por todos os polígonos no conjunto. Veja esclarecimentos nos comentários.
jpmc26
3

Por que não usar uma iteração ou recursividade? 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
fonte
2

Teste este código. é bastante simples em conceito e acredito que você recebe o que procura.

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

e se você deseja que a saída seja armazenada como um arquivo shapefile, 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)})

isso gera -

insira a descrição da imagem aqui

ziggy
fonte
3
Não acredito que isso faça o que o OP quer. Ele devolve as áreas que pelo menos 2 polígonos cobrem, enquanto o OP está procurando apenas áreas cobertas por todos os polígonos no conjunto. Veja esclarecimentos nos comentários. Além disso, ke vsão más escolhas para nomes de variáveis em suas dictcompreensões. Essas variáveis ​​se referem a diferentes elementos de dic.items(), e não a um par de valores-chave. Algo como a, bseria menos enganador.
jpmc26
1
ohh bem sim, eu não entendia o que ele queria dizer
Ziggy
e ponto bem aceite sobre meus k, v escolhas - eu só uso automaticamente k, v quando looping através de um dictionary..didnt dar-lhe muito pensamento
Ziggy