Saya telah mencoba beberapa contoh kode menggunakan pustaka seperti shapefile, fiona, dan ogr untuk mencoba memeriksa apakah suatu titik (x, y) berada dalam batas-batas suatu multipoligon yang dibuat dengan ArcMap (dan dengan demikian dalam format shapefile). Namun tidak ada satu pun contoh yang bekerja dengan baik dengan multipoligon, meskipun mereka baik-baik saja dengan shapefile poligon tunggal biasa. Beberapa cuplikan yang saya coba di bawah:
# First example using shapefile and shapely:
from shapely.geometry import Polygon, Point, MultiPolygon
import shapefile
polygon = shapefile.Reader('shapefile.shp')
polygon = polygon.shapes()
shpfilePoints = []
for shape in polygon:
shpfilePoints = shape.points
polygon = shpfilePoints
poly = Polygon(poly)
point = Point(x, y)
# point in polygon test
if polygon.contains(point):
print 'inside'
else:
print 'OUT'
# Second example using ogr and shapely:
from shapely.geometry import Polygon, Point, MultiPolygon
from osgeo import ogr, gdal
driver = ogr.GetDriverByName('ESRI Shapefile')
dataset = driver.Open("shapefile.shp", 0)
layer = dataset.GetLayer()
for index in xrange(layer.GetFeatureCount()):
feature = layer.GetFeature(index)
geometry = feature.GetGeometryRef()
polygon = Polygon(geometry)
print 'polygon points =', polygon # this prints 'multipoint' + all the points fine
point = Point(x, y)
# point in polygon test
if polygon.contains(point):
print 'inside'
else:
print 'OUT'
Contoh pertama berfungsi dengan baik dengan satu poligon pada satu waktu, tetapi ketika saya memasukkan sebuah titik di dalam salah satu bentuk pada bentuk multipoligon saya, ia mengembalikan "keluar" meskipun itu termasuk dalam salah satu dari banyak bagian.
Karena untuk contoh kedua saya mendapatkan kesalahan "objek bertipe 'Geometri' tidak memiliki len ()" yang saya asumsikan karena bidang geometri tidak dapat dibaca sebagai daftar / array yang diindeks dan normal.
Saya juga mencoba mengganti titik aktual dalam kode poligon seperti yang disarankan di sini untuk memastikan bukan bagian dari kode yang tidak berfungsi. Dan sementara contoh tautan itu bekerja dengan baik dengan bentuk poligon sederhana, saya tidak bisa mendapatkan multipoligon kompleks saya untuk diuji dengan benar.
Jadi saya tidak bisa memikirkan cara lain untuk menguji apakah suatu titik termasuk dalam multipolygon shapefile via python ... Mungkin ada perpustakaan lain di luar sana yang saya lewatkan?
sumber
polygon = Polygon(geometry)
dengan semacam loop coba di mana ia beralih kepolygon = MultiPolygon(geometry)
jika kesalahan itu terjadi.Jawaban:
Shapefile tidak memiliki tipe MultiPolygon (type = Polygon), tetapi mereka tetap mendukungnya (semua cincin disimpan dalam satu fitur = daftar poligon, lihat Mengonversi multipoligon besar menjadi poligon )
Masalah
Jika saya membuka shapefile MultiPolygon, geometri adalah 'Poligon'
Solusi 1 dengan Fiona
Fiona menginterpretasikan fitur tersebut sebagai MultiPolygon dan Anda dapat menerapkan solusi yang disajikan dalam Penggabungan Tata Ruang yang Lebih Efisien dengan Python tanpa QGIS, ArcGIS, PostGIS, dll. (1)
Solusi 2 dengan pyshp (shapefile) dan protokol geo_interface (seperti GeoJSON)
Ini adalah suplemen untuk jawaban xulnik.
Solusi 3 dengan ogr dan protokol geo_interface ( aplikasi Python Geo_interface )
Solusi 4 dengan GeoPandas seperti dalam Tata Ruang yang Lebih Efisien, bergabung dengan Python tanpa QGIS, ArcGIS, PostGIS, dll. (2)
Poin 1,3,5,6 berada dalam batas MultiPolygon
sumber
multi = shape(polys.shapeRecords()[0].shape.__geo_interface__)
Solusi 2? Saya tidak bisa mendapatkan panggilan ke metode shape () darishapefile.py
. Saya bahkan sudah mencobashapefile.Shape()
; ada kelas untuk itu tetapi tidak berhasil.within()
metode ini?from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon
)File "C:\WinPython\python-3.6.5.amd64\lib\site-packages\geopandas\tools\sjoin.py", line 43, in sjoin if left_df.crs != right_df.crs:
AttributeError: 'MultiPolygon' object has no attribute 'crs'
Masalah dalam contoh pertama Anda adalah dalam loop ini:
Ini hanya menambahkan poin fitur terakhir. Saya mencoba pendekatan saya dengan shapefile ini:
Saya mengubah kode Anda menjadi:
Kode di atas dijalankan di Konsol Python QGIS dan hasilnya adalah:
Ini berfungsi dengan baik dan sekarang, Anda dapat memeriksa apakah suatu titik (x, y) berada dalam batas-batas setiap fitur.
sumber
Jika Anda mencoba untuk memeriksa titik lintang, bujur di dalam poligon, pastikan Anda memiliki objek titik yang dibuat oleh yang berikut ini:
Titik mengambil garis bujur, lalu garis lintang dalam argumen. Bukan garis lintang terlebih dahulu. Anda dapat memanggil
polygon_object.within
fungsi untuk memeriksa apakah intinya ada dalam bentuk.sumber