menggunakan Shapely: menerjemahkan antara Polygons dan MultiPolygons

12

[EDIT: solusi untuk ini adalah hanya menggunakan OGR untuk membaca shapefile. Lihat contoh geografi.]

Dalam bentuk ESRI, tidak ada perbedaan antara Poligon dan MultiPoligon. Lebih jauh lagi, tidak ada perbedaan eksplisit antara lubang interior dan cincin eksterior (selain "kidal" dari poligon tertentu).

Jadi setelah membaca sebuah shapefile, saya memiliki daftar urutan koordinat yang menggambarkan cincin, tetapi tanpa proses yang lebih intensif, saya tidak dapat membedakan cincin mana yang merupakan cincin eksterior, lubang interior, atau poligon tambahan.

Tampaknya untuk rupawan 's Polygon dan MultiPolygon konstruktor, harus ada perbedaan yang jelas antara eksterior dan cincin interior, jadi bagaimana saya harus pindah dari daftar yang tidak jelas dari cincin untuk set memerintahkan poligon terpisah, dengan interior jelas ditunjuk dan cincin eksterior ?

Untuk meringkas: jika saya memiliki daftar cincin poligon, tetapi saya tidak tahu cincin mana yang lubang di bagian dalam atau merupakan poligon yang terpisah, bagaimana sebaiknya saya mengurutkannya menjadi poligon yang terpisah dengan lubang bagian dalam yang ditunjuk?

Saya mencari solusi algoritmik sederhana yang dapat saya terapkan dalam python, dapat digunakan untuk memproses ratusan poligon dalam ~ satu menit atau kurang, dan saya melakukan ini untuk melakukan sejumlah besar persimpangan.

BenjaminGolder
sumber
Pertanyaan ini tidak memiliki informasi penting tentang apa yang Anda gunakan untuk membaca Shapefile.
termasuk42
@ inc42 Saya menggunakan python untuk membaca file secara langsung.
BenjaminGolder
Ah, maka bit Shapely itu menyesatkan. Masalah Anda yang sebenarnya adalah mencari tahu cara menentukan "jenis" cincin dalam format Shapefile. :)
inc42

Jawaban:

10

Selanjutnya untuk mengembalikan jawaban tentang cara mendapatkan poligon individu, Anda kemudian dapat menjalankan persimpangan pada semua poligon untuk membuat lubang. Jika dataset Anda berisi poligon yang tumpang tindih meskipun Anda kurang beruntung.

Jelaskan lagi apa yang salah dengan pembaca shapefile yang ada?

Tidakkah lebih mudah untuk mengekspor ID fitur dan nilai M dari shapefile dan kemudian bergabung kembali ke poligon setelah menggunakan pembaca shapefile yang ada?

Untuk multipatch, Anda dapat menggunakan teknik yang sama dalam menetapkan ID poligon ke "patch ID" dan kemudian menambahkan atribut ini kembali ke fitur.

Sunting: Sementara Anda mengatakan Anda tidak ingin menggunakan OGR, kalau-kalau Anda berubah pikiran ..

import ogr
# Get the driver
driver = ogr.GetDriverByName('ESRI Shapefile')
# Open a shapefile
shapefileName = "D:/temp/myshapefile.shp"
dataset = driver.Open(shapefileName, 0)

layer = dataset.GetLayer()
for index in xrange(layer.GetFeatureCount()):
    feature = layer.GetFeature(index)
    geometry = feature.GetGeometryRef()
    #geometry for polygon as WKT, inner rings, outer rings etc. 
    print geometry

Geometri harus berupa output sebagai berikut:

POLYGON ((79285 57742,78741 54273...),(76087 55694,78511 55088,..))

Braket pertama berisi koordinat cincin eksterior, kemudian tanda kurung cincin interior. Jika Anda memiliki nilai Z poin harus dalam format 79285 57742 10 (di mana coord terakhir adalah ketinggian).

Kalau tidak, Anda bisa menggunakan fungsi Shapely Contains dan Within untuk menilai setiap poligon satu sama lain dan menerapkan indeks spasial sebelumnya - http://pypi.python.org/pypi/Rtree/ untuk mempercepat pemrosesan.

geografi
sumber
terima kasih, ini bukan yang saya cari, tetapi saya melihat beberapa klarifikasi yang bisa saya buat untuk pertanyaan saya. Untuk menjawab pertanyaan Anda: 1, karena saya menyortir mereka untuk melakukan banyak persimpangan di tempat pertama. 2, tidak ada, saya hanya perlu satu yang ringan dan yang dapat saya gunakan dengan mudah dari python, dan saya belum belajar ogr. 3, tidak, tidak akan lebih mudah dalam hal ini.
BenjaminGolder
1
Wow! Terima kasih, geografi! Sekarang saya 90% yakin bahwa ogr adalah solusinya di sini. Saya mengalami beberapa masalah instalasi dengan ogr, tetapi sepertinya mereka layak untuk dikerjakan. Jadi ogr mengatur cincin dalam output wkt? dan apakah itu baik-baik saja dengan shapefile 3d (saya menggunakan banyak 3d)?
BenjaminGolder
Sepertinya jawaban untuk pertanyaan saya adalah ya dan ya. Dan terima kasih telah memperkenalkan saya pada rtree. OGR tampaknya menyelesaikan masalah saya sepenuhnya - selama saya dapat mengonfigurasi instalasi dengan benar kali ini.
BenjaminGolder
Instalasi ulang OGR - ingat Anda akan membutuhkan binari untuk Windows dan binding Python yang cocok.
geografi
9

Pertama, gunakan ogr untuk membuka shapefile:

from osgeo import ogr
source = ogr.Open("mpolys.shp")
layers =  source.GetLayerByName("mpoly")
len(layers)
1

mengubah geometri shapefile menjadi geometri bentuk

from shapely.wkb import loads
element=layers[0] #(because lenght of layer =1, else you need "for element in layers: ...")
geom = loads(element.GetGeometryRef().ExportToWkb())
geom.geom_type
'MultiPolygon'
print geom
MULTIPOLYGON ((..... # the geometry in shapely wkt format

Untuk poligon dalam multipolygon:

poly=[]
for pol in geom:
    poly.append(pol)
poly[0]
<shapely.geometry.polygon.Polygon object at 0x00B82CB0>
poly[0].geom_type
'Polygon'
print(poly[1])
POLYGON ((.... # the geometry in shapely wkt format

Dan sekarang, Anda dapat menggunakan semua fungsi shapely ( rupawan )

Constantinius
sumber
1

Saya tidak terlalu terbiasa dengan bagaimana sebenarnya poligon disimpan dalam file bentuk, tetapi - bukankah seharusnya cincin poligon menjadi loop tertutup jika dan hanya jika koordinat awal diulang? Jadi, jika Anda membandingkan setiap koordinat berikutnya dengan koordinat awal, Anda akan menemukan titik pertama di mana poligon ditutup. Jika itu adalah koordinat terakhir dari poligon, itu adalah poligon sederhana, jika tidak, itu adalah multipoligon dan membutuhkan pemrosesan loop lainnya.

Itu mungkin 'pemrosesan lebih intensif' yang ingin Anda hindari, tetapi itu sebenarnya hanya iterasi melalui koordinat yang datang secara gratis ketika Anda harus membacanya.

relet
sumber
Saya minta maaf - pertanyaan saya agak menyesatkan dan saya sudah mengeditnya. Saya memiliki cincin poligon, dan saya tahu ketika saya memiliki daftar lebih dari satu cincin poligon, tetapi selain dari urutan poin searah jarum jam atau berlawanan arah, saya tidak tahu apakah itu cincin interior atau eksterior. Jika itu adalah cincin interior, saya tidak bisa memastikan cincin eksterior mana yang mereka miliki, selain mengukur lokasi mereka.
BenjaminGolder
Saya melihat. Saya juga senang melihat Anda menemukan jalan ke ogr. ;)
relet
-2

Seperti ditunjukkan oleh @ inc42 , pada halaman 8 ESRI Shapefile Deskripsi Teknis: Buku Putih ESRI - Juli 1998 (halaman 12 dari 34 dalam PDF) ada diskusi tentang konten rekaman poligon yang mungkin seperti apa yang Anda cari:

Poligon dapat berisi beberapa cincin luar. Urutan simpul atau orientasi untuk cincin menunjukkan sisi cincin mana yang merupakan bagian dalam poligon.

PolyGeo
sumber
p.10 memilikinya "Poligon dapat berisi beberapa cincin luar. Urutan simpul atau orientasi untuk cincin menunjukkan sisi cincin mana yang merupakan bagian dalam poligon."
termasuk42
@ inc42 diperbarui meskipun saya pikir nomor halaman harus 8 (atau 12 dari 34).
PolyGeo