Check if a point falls within a multipolygon with Python

Question

I have tried several examples of code using libraries such as shapefile, fiona, and ogr to attempt to check whether a point (x, y) falls within the boundaries of a multipolygon created with ArcMap (and thus in shapefile format). However none of the examples work well with multipolygons, although they do fine with regular, single polygon shapefiles. Some snippets I tried are below:

# 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(polygon)
point = Point(x, y)
# point in polygon test
if poly.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'

The first example works fine with a single polygon at a time, but when I input a point within one of the shapes in my multipolygon shapefile it returns "out" even though it does fall inside one of the many parts.

For for the second example I get an error "object of type 'Geometry' has no len()" which I assume is because the geometry field can't be read as a normal, indexed list/array.

I additionally tried to replace the actual point in polygon code as suggested here to make sure it wasn't that part of the code that didn't work. And while that link's examples work fine with simple polygon shapefiles I can't get my complex multipolygon to test properly.

So I can't think of any other way to test whether a point falls within a multipolygon shapefile via python... Maybe there are other libraries out there I'm missing?

• obrl_soil: your second example looks like it might be coercing multipolygon to polygon? It might only be checking the point against the first part of the multipolygon. Try moving the point around to different parts and see if the check ever succeeds.
• spartmar: @obrl_soil Thank you for your suggestion. However, the second example never works because of the error message I described above (object of type 'Geometry' has no len())" whether I try MultiPolygon (geometry) or simply Polygon(geometry). I have tried many points in the first example and only those within the main polygon work. Hope this clarification helps.
• obrl_soil: Yeah, I think you need to replace polygon = Polygon(geometry) with some kind of try loop where it switches to polygon = MultiPolygon(geometry) if that error occurs.
• xunilk: The problem in your first example is in the first loop.

Answer 1

Shapefiles have no type MultiPolygon (type = Polygon), but they support them anyway (all rings are stored in one feature = list of polygons, look at Converting huge multipolygon to polygons)

The problem

enter image description here

If I open a MultiPolygon shapefile, the geometry is 'Polygon'

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

Solution 1 with 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 interprets the feature as a MultiPolygon and you can apply the solution presented in 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)

Solution 2 with pyshp (shapefile) and the geo_interface (GeoJSON like) protocol

This is a supplement to the answer of 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)

Solution 3 with ogr and the geo_interface protocol (Python Geo_interface applications)

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)

Solution 4 with GeoPandas as in 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]} 

The points 1,3,5,6 falls within the boundaries of the MultiPolygon

• pstatix: Slightly old thread here, but how do you call multi = shape(polys.shapeRecords()[0].shape.__geo_interface__) in Solution 2? I cannot get a call to a shape() method from shapefile.py. I have even tried shapefile.Shape(); there is a class for it but it doesnt work.
• pstatix: Further, where do you get the within() method from?
• gene: from Shapely (from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon)
• Aaron Bramson: I get this error using Solution 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'
• kevininhe: Reading this helped me to understand better gene's answer regarding geopandas's join usage: geopandas.org/en/stable/docs/user_guide/…

Answer 2

The problem in your first example is in this loop:

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

It only appends the last feature points. I tried out my approach with this shapefile:

enter image description here

I modified your code to:

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

Above code was run at the Python Console of QGIS and the result was:

enter image description here

It works perfectly and now, you can check whether a point (x, y) falls within the boundaries of each feature.

Answer 3

If you are trying to check a latitude, longitude point within a polygon, make sure you you have point object is created by the following:

from shapely.geometry.point import Point
Point(LONGITUDE, LATITUDE)
.. 
point.within(poly) # Returns true if the point within the polygon
poly.contains(point) # Returns true if the polygon contains the point. Works as well.

Point takes longitude, then latitude in the argument. Not latitude first. You can call polygon_object.contains or polygon_object.within function to check if the point is within the shape.

• mah65: Any suggestion how I can do it fast and efficient when I have many points and many polygons? Looping over all is very costly.
• Ahmed: You can use postgis to store your polygons and use one of the built-in functions of the server.
• Skippy le Grand Gourou: @mah65 If the distribution is very imbalanced (e.g. most points are in a few countries), ordering the polygons by order of most frequently seen can greatly increase performances. You can sample a number of random points from your data to find the most efficient ordering.