He intentado varios ejemplos de código utilizando bibliotecas como shapefile, fiona y ogr para intentar verificar si un punto (x, y) se encuentra dentro de los límites de un multipolígono creado con ArcMap (y, por lo tanto, en formato shapefile). Sin embargo, ninguno de los ejemplos funciona bien con multipolígonos, aunque funcionan bien con archivos de forma de polígonos individuales regulares. Algunos fragmentos que probé están a continuación:
# 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'
El primer ejemplo funciona bien con un solo polígono a la vez, pero cuando ingreso un punto dentro de una de las formas en mi archivo de formas multipolígono, regresa "fuera" a pesar de que cae dentro de una de las muchas partes.
Para el segundo ejemplo, recibo un error "el objeto de tipo 'Geometry' no tiene len ()", lo que supongo es porque el campo de geometría no puede leerse como una lista / matriz indexada normal.
También intenté reemplazar el punto real en el código de polígono como se sugirió aquí para asegurarme de que no era esa parte del código que no funcionaba. Y aunque los ejemplos de ese enlace funcionan bien con archivos de forma de polígonos simples, no puedo hacer que mi multipolígono complejo se pruebe correctamente.
Así que no puedo pensar en otra forma de probar si un punto cae dentro de un archivo de forma multipolígono a través de Python ... ¿Tal vez hay otras bibliotecas que me faltan?
fuente
polygon = Polygon(geometry)con algún tipo de bucle de prueba donde cambiapolygon = MultiPolygon(geometry)si se produce ese error.Respuestas:
Los archivos de forma no tienen tipo MultiPolygon (type = Polygon), pero los admiten de todos modos (todos los anillos se almacenan en una característica = lista de polígonos, mira Convertir gran polígono en polígonos )
El problema
Si abro un archivo de forma MultiPolygon, la geometría es 'Polygon'
Solución 1 con Fiona
Fiona interpreta la función como un MultiPolygon y puede aplicar la solución presentada en la unión espacial más eficiente en Python sin QGIS, ArcGIS, PostGIS, etc. (1)
Solución 2 con pyshp (shapefile) y el protocolo geo_interface (como GeoJSON)
Este es un suplemento a la respuesta de xulnik.
Solución 3 con ogr y el protocolo geo_interface ( aplicaciones Python Geo_interface )
Solución 4 con GeoPandas como en combinación espacial más eficiente en Python sin QGIS, ArcGIS, PostGIS, etc. (2)
Los puntos 1,3,5,6 caen dentro de los límites del MultiPolygon
fuente
multi = shape(polys.shapeRecords()[0].shape.__geo_interface__)en la Solución 2? No puedo recibir una llamada a un método shape () desdeshapefile.py. Incluso lo he intentadoshapefile.Shape(); hay una clase para eso pero no funciona.within()método?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'El problema en su primer ejemplo está en este bucle:
Solo agrega los últimos puntos de características. Probé mi enfoque con este shapefile:
Modifiqué tu código para:
El código anterior se ejecutó en la Consola Python de QGIS y el resultado fue:
Funciona perfectamente y ahora, puede verificar si un punto (x, y) cae dentro de los límites de cada entidad.
fuente
Si está tratando de verificar un punto de latitud y longitud dentro de un polígono, asegúrese de tener un objeto de punto creado de la siguiente manera:
El punto toma longitud, luego latitud en el argumento. No la latitud primero. Puede llamar a la
polygon_object.withinfunción para verificar si el punto está dentro de la forma.fuente