Converting huge multipolygon to polygons

Question

I have a shapefile with some huge multipolygons, with 100.000's of parts. What would be the easiest way to split them to singlepart polygons? I'm looking for something like QGIS ”Multipart to singlepart” function, but the file is to large for QGIS to handle. I'm guessing that there is probably already some Python module that can do it for me. Any tips?

Answer 1

Shapefiles have no type MultiPolygon (type = Polygon), but they support them anyway (all rings are stored in one polygon = list of polygons, look at GDAL: ESRI Shapefile)

It is easier with Fiona and Shapely:

import fiona
from shapely.geometry import shape, mapping

# open the original MultiPolygon file
with fiona.open('multipolygons.shp') as source:
    # create the new file: the driver and crs are the same
    # for the schema the geometry type is "Polygon" instead
    output_schema = dict(source.schema)  # make an independant copy
    output_schema['geometry'] = "Polygon"

    with fiona.open('output.shp', 'w', 
                    driver=source.driver,
                    crs=source.crs,
                    schema=output_schema) as output:

        # read the input file
        for multi in source:

           # extract each Polygon feature
           for poly in shape(multi['geometry']):

              # write the Polygon feature
              output.write({
                  'properties': multi['properties'],
                  'geometry': mapping(poly)
              })

Answer 2

from GDAL mailing list using python

import os
from osgeo import ogr

def multipoly2poly(in_lyr, out_lyr):
    for in_feat in in_lyr:
        geom = in_feat.GetGeometryRef()
        if geom.GetGeometryName() == 'MULTIPOLYGON':
            for geom_part in geom:
                addPolygon(geom_part.ExportToWkb(), out_lyr)
        else:
            addPolygon(geom.ExportToWkb(), out_lyr)

def addPolygon(simplePolygon, out_lyr):
    featureDefn = out_lyr.GetLayerDefn()
    polygon = ogr.CreateGeometryFromWkb(simplePolygon)
    out_feat = ogr.Feature(featureDefn)
    out_feat.SetGeometry(polygon)
    out_lyr.CreateFeature(out_feat)
    print 'Polygon added.'

from osgeo import gdal
gdal.UseExceptions()
driver = ogr.GetDriverByName('ESRI Shapefile')
in_ds = driver.Open('data/multipoly.shp', 0)
in_lyr = in_ds.GetLayer()
outputshp = 'data/poly.shp'
if os.path.exists(outputshp):
    driver.DeleteDataSource(outputshp)
out_ds = driver.CreateDataSource(outputshp)
out_lyr = out_ds.CreateLayer('poly', geom_type=ogr.wkbPolygon)
multipoly2poly(in_lyr, out_lyr)