Return to Answer

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)

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)

Solution 4 with GeoPandas as in More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (2)

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)

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)

Solution 4 with GeoPandas as in More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (2)

 import fiona  from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon  multi = fiona.open("multipol.shp") multipol = multi.next() # only one feature in the shapefile  print multipol  {'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)])} 

for i, pt in enumerate(points): point = shape(pt['geometry']) if point.within(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) 

This is a Supplement to the answer of xulnik.

Solution 3 with GeoPandas as in More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (2)

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)])} 

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) 

This is a supplement to the answer of xulnik.

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 fiona # creation of a MultiPolygon shapefile polys = fiona.open("multipoly.shp") print len(polys) 3 # 3 polygons in the shapefile polys.schema {'geometry': 'Polygon', 'properties': OrderedDict([(u'id', 'int:10')])} # creation of a MultiPolygon from shapely.geometry import shape,mapping, Point, Polygon, MultiPolygon multi = MultiPolygon([shape(i['geometry']) for i in polys]) # schema of the new shapefile mult = polys.schema mult['geometry'] = 'MultiPolygon' print mult {'geometry': 'MultiPolygon', 'properties': OrderedDict([(u'id', 'int:10')])} with fiona.open("multipol.shp", 'w',driver='ESRI Shapefile', schema=mult) as output:  output.write({'geometry': mapping(multi), 'properties':{'id': 1}}) 

Now if I open the new shapefile, the geometry is 'Polygon' but with only one feature (= list of the 3 polygons)

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

Solution 1 expanding the list of polygons:

polygons = [poly for poly in multi for multi in multipolys] len(polygons) 3 

And now 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')] # iterate through points for i, pt in enumerate(points): point = shape(pt['geometry']) #iterate through polygons for j, poly in enumerate(polygons):  if point.within(poly):  print i, shape(points[i]['geometry']) #, polygons[j] 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 : recreating the MultiPolygon with Fiona

multi = MultiPolygon([poly for poly in multi for multi in multipolys]) for i, pt in enumerate(points): point = shape(pt['geometry']) if point.within(multi): 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 3 with pyshp (shapefile) and the geo_interface protocol

This is a Supplement to the answer of xulnik

pts = shapefile.Reader("points.shp") polys = shapefile.Reader("multipol.shp") points = [pt.shape.__geo_interface__ for pt in pts.shapeRecords()] multi = [shape(pol.shape.__geo_interface__) for pol in poly.shapeRecords()] for i, pt in enumerate(points): point = shape(pt) if point.within(multi[0]): #first multipolygon 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 4 with GeoPandas as in More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (2)

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

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  multi = fiona.open("multipol.shp")   multipol = multi.next() # only one feature in the shapefile print multipol  {'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)

for i, pt in enumerate(points): point = shape(pt['geometry']) if point.within(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 GeoPandas as in More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (2)

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