# Check if a point falls within a multipolygon with Python

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?

• 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. Aug 26, 2016 at 22:44
• @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. Aug 26, 2016 at 22:48
• 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. Aug 27, 2016 at 6:08
• The problem in your first example is in the first loop. Aug 27, 2016 at 6:58

## 3 Answers

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

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

• 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. Aug 31, 2017 at 23:31
• Further, where do you get the `within()` method from? Aug 31, 2017 at 23:45
• from Shapely (`from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon`)
– gene
Sep 4, 2017 at 15:04
• 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'` Apr 8, 2019 at 5:59
• Reading this helped me to understand better gene's answer regarding geopandas's join usage: geopandas.org/en/stable/docs/user_guide/… Mar 25 at 19:01

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:

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:

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

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.

• Any suggestion how I can do it fast and efficient when I have many points and many polygons? Looping over all is very costly. Nov 8, 2020 at 6:05
• You can use postgis to store your polygons and use one of the built-in functions of the server. Nov 9, 2020 at 7:32