Mendapatkan persimpangan beberapa poligon secara efisien dalam Python

12

Saya ingin mendapatkan persimpangan beberapa poligon. Menggunakan shapelypaket Python , saya dapat menemukan persimpangan dua poligon menggunakan intersectionfungsi. Apakah ada fungsi efisien yang serupa untuk memperoleh persimpangan beberapa poligon?

Berikut ini cuplikan kode untuk memahami maksud saya:

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) 

Persimpangan dua lingkaran dapat ditemukan oleh circle1.intersection(circle2). Saya dapat menemukan persimpangan ketiga lingkaran dengan circle1.intersection(circle2).intersection(circle3). Namun, pendekatan ini tidak dapat dijual ke sejumlah besar poligon karena membutuhkan kode yang semakin banyak. Saya ingin fungsi yang mengambil jumlah poligon yang sewenang-wenang dan mengembalikan persimpangan mereka.

terpecah
sumber
Saya pikir mungkin menyimpan coords di kamus dan loop melalui sementara menggunakan kombinasi impor itertools. Saya akan segera memposting
ziggy
Apa yang Anda maksud dengan "persimpangan mereka"? Maksud Anda semua area yang bersinggungan dengan setidaknya satu poligon lain, atau area yang semua input berpotongan?
jpmc26
Maksud saya persimpangan semua poligon, paling tidak satu.
sempalan
Anda harus mengklarifikasi ini di atas (mungkin dengan contoh output). Saya cukup yakin sebagian besar jawaban tidak berperilaku seperti yang Anda inginkan. (Dan fakta bahwa beberapa penjawab telah salah paham adalah bukti yang cukup bahwa pertanyaan itu perlu diklarifikasi.)
jpmc26
1
@ jpmc26 Saya baru saja menambahkan pembaruan ke jawaban saya di mana rtree digunakan. Pendekatannya lebih efisien dan terukur sekarang. Semoga ini membantu!
Antonio Falciano

Jawaban:

7

Salah satu pendekatan yang mungkin dapat mempertimbangkan kombinasi pasangan poligon, persimpangan mereka dan akhirnya penyatuan semua persimpangan melalui penyatuan bertingkat (seperti yang disarankan di sini ):

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

Pendekatan yang lebih efisien harus menggunakan indeks spasial, seperti Rtree , untuk berurusan dengan banyak geometri (bukan kasus tiga lingkaran):

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
sumber
Saya tidak percaya ini melakukan apa yang diinginkan OP. Ini memberikan kembali area yang setidaknya 2 poligon menutupi, sedangkan OP hanya mencari area yang dicakup oleh semua poligon dalam set. Lihat klarifikasi dalam komentar.
jpmc26
3

Mengapa tidak menggunakan iterasi atau rekursif? sesuatu seperti :

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

Berikan kode ini kesempatan. itu cukup sederhana dalam konsep dan saya percaya membuat Anda mendapatkan apa yang Anda cari.

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

dan jika Anda ingin output disimpan sebagai penggunaan, gunakan 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)})

output ini -

masukkan deskripsi gambar di sini

ziggy
sumber
3
Saya tidak percaya ini melakukan apa yang diinginkan OP. Ini memberikan kembali area yang setidaknya 2 poligon menutupi, sedangkan OP hanya mencari area yang dicakup oleh semua poligon dalam set. Lihat klarifikasi dalam komentar. Selain itu, kdan vmerupakan pilihan yang buruk untuk nama variabel dalam dictpemahaman Anda . Variabel-variabel tersebut masing-masing merujuk pada elemen berbeda dic.items(), bukan pasangan nilai kunci. Sesuatu seperti a, bakan kurang menyesatkan.
jpmc26
1
ohh oke ya saya tidak mengerti apa yang dia maksud
ziggy
dan titik diambil dengan baik tentang k saya, v pilihan - saya hanya secara otomatis menggunakan k, v ketika perulangan melalui
kamus..tidak