8

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?

  • Do you have the option of loading the layer into PostGIS? – Get Spatial Aug 28 '14 at 20:44
  • What is your layer stored in now? It could be that changing your storage format would allow it to work in QGIS, just accounting for differences in efficiencies. – Get Spatial Aug 28 '14 at 20:46
  • 1
    It's in a shapefile, so @Barrett's solution below is perfect! – leo Aug 29 '14 at 5:38
11

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

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.