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