Verifique si un punto cae dentro de un multipolígono con Python

Question

He probado varios ejemplos de código utilizando bibliotecas como shapefile, fiona y ogr para intentar verificar si un punto (x, y) cae 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 shapefiles de polígonos regulares simples. Algunos fragmentos que probé son los siguientes:

# Primer ejemplo usando shapefile y 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(polygon)

point = Point(x, y)
# prueba de punto en polígono
if poly.contains(point):
    print 'dentro'
else:
    print 'FUERA'

# Segundo ejemplo usando ogr y 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 'puntos del polígono =', polygon  # esto imprime 'multipunto' + todos los puntos bien

point = Point(x, y)
# prueba de punto en polígono
if polygon.contains(point):
    print 'dentro'
else:
    print 'FUERA'

El primer ejemplo funciona bien con un solo polígono a la vez, pero cuando ingreso un punto dentro de una de las formas de mi shapefile multipolígono, devuelve "fuera" aunque sí cae dentro de una de las muchas partes.

Para el segundo ejemplo, obtengo un error "object of type 'Geometry' has no len()" que supongo es porque el campo de geometría no se puede leer como una lista/matriz normal e indexada.

Además, intenté reemplazar el código real de punto en polígono como se sugiere aquí para asegurarme de que no era esa parte del código la que no funcionaba. Y aunque los ejemplos de ese enlace funcionan bien con shapefiles de polígonos simples, no puedo hacer que mi multipolígono complejo se pruebe correctamente.

Así que no puedo pensar en ninguna otra forma de probar si un punto cae dentro de un shapefile multipolígono a través de python... ¿Quizás hay otras bibliotecas por ahí que me estoy perdiendo?

Answer 1

Los shapefiles no tienen el tipo MultiPolygon (tipo = Polígono), pero de todos modos los admiten (todos los anillos se almacenan en una característica = lista de polígonos, ver Converting huge multipolygon to polygons)

El problema

Si abro un shapefile MultiPolygon, la geometría es 'Polígono'

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

Solución 1 con Fiona

import fiona
from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon
multipol = fiona.open("multipol.shp")
multi= multipol.next() # solo una característica en el 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 interpreta la característica como un MultiPolygon y puedes aplicar la solución presentada en More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (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)

Solución 2 con pyshp (shapefile) y el protocolo geo_interface (similar a GeoJSON)

Esto es un suplemento a la respuesta de 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 polígono
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)

Solución 3 con ogr y el protocolo geo_interface (Aplicaciones Python Geo_interface)

from osgeo import ogr
import json
def records(file):  
    # generador 
    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 característica
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)

Solución 4 con GeoPandas como en More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (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]} 

Los puntos 1, 3, 5, 6 caen dentro de los límites del MultiPolygon

Answer 2

El problema en tu primer ejemplo está en este bucle:

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

Solo añade los puntos de la última característica. Probé mi enfoque con este archivo de formas:

Modifiqué tu código a:

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

El código anterior se ejecutó en la Consola de Python de QGIS y el resultado fue:

Funciona perfectamente y ahora puedes verificar si un punto (x, y) cae dentro de los límites de cada característica.

Answer 3

Si estás intentando comprobar si un punto de latitud, longitud está dentro de un polígono, asegúrate de que tengas el objeto punto creado de la siguiente manera:

from shapely.geometry.point import Point
Point(LONGITUDE, LATITUDE)
..    
point.within(poly)   # Devuelve verdadero si el punto está dentro del polígono
poly.contains(point) # Devuelve verdadero si el polígono contiene el punto. Funciona también.

El Point recibe primero la longitud y luego la latitud como argumento. No la latitud primero. Puedes llamar a la función polygon_object.contains o polygon_object.within para comprobar si el punto está dentro de la forma.