Union intersected polygons with python and save the result in a new shapefile

Question

I'm trying to make an union with overlapping polygons and save the result in a new shapefile.

from osgeo import ogr
import os
test = ogr.Open(r"F:\Python\Shapefile\BuiltBuf7000.shp")
layer = test.GetLayer()
print layer.GetGeomType()
3 # -> polygons
union_poly = ogr.Geometry(ogr.wkbPolygon)
for feature in layer:
geom =feature.GetGeometryRef()
union_poly = union_poly.Union(geom)

print union_poly

When I run the code I got this result :

Now, I would like to save the result in a new shapefile, Does someone know how to do it ?

Answer 1

First thing, you need to create an empty shapefile:

# Create datasource.
driver = ogr.GetDriverByName('ESRI Shapefile')
try:
    ds_out = driver.CreateDataSource(filepath)
except:
    raise Exception("File already exists.")

# Grab current layer's SRS and geometry type.
srs = layer_in.GetSpatialRef()
geom_type = layer_in.GetGeomType()

# Create empty layer.
layer_out = ds_out.CreateLayer(name = *name you want*, srs = srs, geom_type = geom_type)

# Copy fields <-- not required.
layer_defn = layer_in.GetLayerDefn()
n_fields = layer_defn.GetFieldCount()
for i in range(n_fields):
    layer_out.CreateField(layer_defn.GetFieldDefn(i))

or, you know, just load an existing empty shapefile:

driver = ogr.GetDriverByName('ESRI Shapefile')
ds_out = driver.Open(filepath)
layer_out = ds_out.GetLayer()

To save your created polygon:

# Create feature.
layer_defn = layer_out.GetLayerDefn()
feat_out = ogr.Feature(layer_defn)

# Load geometry in feature, and feature in layer.
feat_out.SetGeometry(union_poly)
layer_out.CreateFeature(feat_out)

# Set attributes <-- not required.
feat_to.SetField(field_name, value)

Don't forget to feat_out.Destroy() and ds_out.Destroy() afterwards!