Converting huge multipolygon to polygons
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)
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)
})