Periksa apakah suatu titik termasuk dalam multipoligon dengan Python

13

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?

Spartmar
sumber
Contoh kedua Anda terlihat seperti itu mungkin memaksa multipolygon ke polygon? Mungkin hanya memeriksa titik terhadap bagian pertama dari multipolygon. Coba pindahkan titik ke bagian yang berbeda dan lihat apakah pemeriksaan berhasil.
obrl_soil
@obrl_soil Terima kasih atas saran Anda. Namun, contoh kedua tidak pernah berfungsi karena pesan kesalahan yang saya jelaskan di atas (objek bertipe 'Geometri' tidak memiliki len ()) "apakah saya mencoba MultiPolygon (geometri) atau hanya Polygon (geometry). Saya telah mencoba banyak titik di contoh pertama dan hanya mereka yang ada dalam pekerjaan poligon utama. Semoga klarifikasi ini membantu
spartmar
Ya, saya pikir Anda perlu mengganti polygon = Polygon(geometry)dengan semacam loop coba di mana ia beralih ke polygon = MultiPolygon(geometry)jika kesalahan itu terjadi.
obrl_soil
Masalah pada contoh pertama Anda adalah pada loop pertama.
xunilk

Jawaban:

24

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

masukkan deskripsi gambar di sini

Jika saya membuka shapefile MultiPolygon, geometri adalah 'Poligon'

multipolys = fiona.open("multipol.shp")
multipolys.schema
{'geometry': 'Polygon', 'properties': OrderedDict([(u'id', 'int:10')])}
len(multipolys)
1

Solusi 1 dengan Fiona

import fiona
from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon
multipol = fiona.open("multipol.shp")
multi= multipol.next() # only one feature in the shapefile
print multi
{'geometry': {'type': 'MultiPolygon', 'coordinates': [[[(-0.5275288092189501, 0.5569782330345711), (-0.11779769526248396, 0.29065300896286816), (-0.25608194622279135, 0.01920614596670933), (-0.709346991037132, -0.08834827144686286), (-0.8629961587708066, 0.18309859154929575), (-0.734955185659411, 0.39820742637644047), (-0.5275288092189501, 0.5569782330345711)]], [[(0.19974391805377723, 0.060179257362355965), (0.5480153649167734, 0.1293213828425096), (0.729833546734955, 0.03969270166453265), (0.8143405889884763, -0.13956466069142115), (0.701664532650448, -0.38540332906530095), (0.4763124199743918, -0.5006402048655569), (0.26888604353393086, -0.4238156209987196), (0.18950064020486557, -0.2291933418693981), (0.19974391805377723, 0.060179257362355965)]], [[(-0.3764404609475033, -0.295774647887324), (-0.11523687580025621, -0.3597951344430217), (-0.033290653008962945, -0.5800256081946222), (-0.11523687580025621, -0.7413572343149808), (-0.3072983354673495, -0.8591549295774648), (-0.58898847631242, -0.6927016645326505), (-0.6555697823303457, -0.4750320102432779), (-0.3764404609475033, -0.295774647887324)]]]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'id', 1)])}

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)

points= ([pt for pt  in fiona.open("points.shp")])
for i, pt in enumerate(points):
    point = shape(pt['geometry'])
    if point.within(shape(multi['geometry'])):
         print i, shape(points[i]['geometry'])
1 POINT (-0.58898847631242 0.17797695262484)
3 POINT (0.4993597951344431 -0.06017925736235585)
5 POINT (-0.3764404609475033 -0.4750320102432779)
6 POINT (-0.3098591549295775 -0.6312419974391805)

Solusi 2 dengan pyshp (shapefile) dan protokol geo_interface (seperti GeoJSON)

Ini adalah suplemen untuk jawaban xulnik.

import shapefile
pts = shapefile.Reader("points.shp")
polys = shapefile.Reader("multipol.shp")
points = [pt.shape.__geo_interface__ for pt in pts.shapeRecords()]
multi = shape(polys.shapeRecords()[0].shape.__geo_interface__) # 1 polygon
print multi
MULTIPOLYGON (((-0.5275288092189501 0.5569782330345711, -0.117797695262484 0.2906530089628682, -0.2560819462227913 0.01920614596670933, -0.7093469910371319 -0.08834827144686286, -0.8629961587708066 0.1830985915492958, -0.734955185659411 0.3982074263764405, -0.5275288092189501 0.5569782330345711)), ((0.1997439180537772 0.06017925736235596, 0.5480153649167734 0.1293213828425096, 0.729833546734955 0.03969270166453265, 0.8143405889884763 -0.1395646606914211, 0.701664532650448 -0.3854033290653009, 0.4763124199743918 -0.5006402048655569, 0.2688860435339309 -0.4238156209987196, 0.1895006402048656 -0.2291933418693981, 0.1997439180537772 0.06017925736235596)), ((-0.3764404609475033 -0.295774647887324, -0.1152368758002562 -0.3597951344430217, -0.03329065300896294 -0.5800256081946222, -0.1152368758002562 -0.7413572343149808, -0.3072983354673495 -0.8591549295774648, -0.58898847631242 -0.6927016645326505, -0.6555697823303457 -0.4750320102432779, -0.3764404609475033 -0.295774647887324)))
for i, pt in enumerate(points):
    point = shape(pt)
    if point.within(multi): 
        print i, shape(points[i])
1 POINT (-0.58898847631242 0.17797695262484)
3 POINT (0.4993597951344431 -0.06017925736235585)
5 POINT (-0.3764404609475033 -0.4750320102432779)
6 POINT (-0.3098591549295775 -0.6312419974391805)

Solusi 3 dengan ogr dan protokol geo_interface ( aplikasi Python Geo_interface )

from osgeo import ogr
import json
def records(file):  
    # generator 
    reader = ogr.Open(file)
    layer = reader.GetLayer(0)
    for i in range(layer.GetFeatureCount()):
        feature = layer.GetFeature(i)
        yield json.loads(feature.ExportToJson())

points  = [pt for pt in records("point_multi_contains.shp")]
multipol = records("multipol.shp")
multi = multipol.next() # 1 feature
for i, pt in enumerate(points):
     point = shape(pt['geometry'])
     if point.within(shape(multi['geometry'])):
          print i, shape(points[i]['geometry'])

1 POINT (-0.58898847631242 0.17797695262484)
3 POINT (0.499359795134443 -0.060179257362356)
5 POINT (-0.376440460947503 -0.475032010243278)
6 POINT (-0.309859154929577 -0.631241997439181)

Solusi 4 dengan GeoPandas seperti dalam Tata Ruang yang Lebih Efisien, bergabung dengan Python tanpa QGIS, ArcGIS, PostGIS, dll. (2)

import geopandas
point = geopandas.GeoDataFrame.from_file('points.shp') 
poly  = geopandas.GeoDataFrame.from_file('multipol.shp')
from geopandas.tools import sjoin
pointInPolys = sjoin(point, poly, how='left')
grouped = pointInPolys.groupby('index_right')
list(grouped)
[(0.0,      geometry                               id_left  index_right id_right  

1  POINT (-0.58898847631242 0.17797695262484)       None      0.0        1.0 
3  POINT (0.4993597951344431 -0.06017925736235585)  None      0.0        1.0
5  POINT (-0.3764404609475033 -0.4750320102432779)  None      0.0        1.0 
6  POINT (-0.3098591549295775 -0.6312419974391805)  None      0.0        1.0 ]
print grouped.groups
{0.0: [1, 3, 5, 6]} 

Poin 1,3,5,6 berada dalam batas MultiPolygon

gen
sumber
Utas yang sedikit usang di sini, tetapi bagaimana Anda memanggil multi = shape(polys.shapeRecords()[0].shape.__geo_interface__)Solusi 2? Saya tidak bisa mendapatkan panggilan ke metode shape () dari shapefile.py. Saya bahkan sudah mencoba shapefile.Shape(); ada kelas untuk itu tetapi tidak berhasil.
pstatix
Selanjutnya, dari mana Anda mendapatkan within()metode ini?
pstatix
1
dari Shapely ( from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon)
gen
Saya mendapatkan kesalahan ini menggunakan Solusi 4: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'
Aaron Bramson
6

Masalah dalam contoh pertama Anda adalah dalam loop ini:

...
shpfilePoints = []
for shape in polygon:
    shpfilePoints = shape.points
...

Ini hanya menambahkan poin fitur terakhir. Saya mencoba pendekatan saya dengan shapefile ini:

masukkan deskripsi gambar di sini

Saya mengubah kode Anda menjadi:

from shapely.geometry import Polygon, Point, MultiPolygon
import shapefile 

path = '/home/zeito/pyqgis_data/polygon8.shp'

polygon = shapefile.Reader(path) 

polygon = polygon.shapes() 

shpfilePoints = [ shape.points for shape in polygon ]

print shpfilePoints

polygons = shpfilePoints

for polygon in polygons:
    poly = Polygon(polygon)
    print poly

Kode di atas dijalankan di Konsol Python QGIS dan hasilnya adalah:

masukkan deskripsi gambar di sini

Ini berfungsi dengan baik dan sekarang, Anda dapat memeriksa apakah suatu titik (x, y) berada dalam batas-batas setiap fitur.

xunilk
sumber
0

Jika Anda mencoba untuk memeriksa titik lintang, bujur di dalam poligon, pastikan Anda memiliki objek titik yang dibuat oleh yang berikut ini:

from shapely.geometry.point import Point
Point(LONGITUDE, LATITUDE)
..
poly.within(point) # Returns true if the point within the 

Titik mengambil garis bujur, lalu garis lintang dalam argumen. Bukan garis lintang terlebih dahulu. Anda dapat memanggil polygon_object.withinfungsi untuk memeriksa apakah intinya ada dalam bentuk.

Ahmed
sumber