OGR Layer Intersection
There are some errors in your script but it is not the most important problem:
You cannot create a valid shapefile without specifying the geometry of the layer:
driver = ogr.GetDriverByName('ESRI Shapefile')
dstshp = driver.CreateDataSource('SomeFilename.shp')
dstlayer = dstshp.CreateLayer('mylayer',geom_type=ogr.wkbPolygon)
And you don't know a priori (upfront) the geometry of the resulting intersection layer. The intersection of two polygon layers is different from the intersection of a polygon layer and a polyline layer for example.
For that, you can get the geometry of the intersection by:
For example (with two polygons shapefiles):
layer1.GetGeomType()
3 # -> polygon
# create an empty geometry of the same type
union1=ogr.Geometry(3)
# union all the geometrical features of layer 1
for feat in layer1:
geom =feat.GetGeometryRef()
union1 = union1.Union(geom)
# same for layer2
union2=ogr.Geometry(layer2.GetGeomType())
for feat in layer2:
geom =feat.GetGeometryRef()
union2 = union2.Union(geom)
# intersection
intersection = union1.Intersection(union2)
print intersection.GetGeometryName()
'MultiPolygon'
At this stage, you can save the resulting geometry to a shapefile (without the fields of the original layers):
dstshp = driver.CreateDataSource('SomeotherFilename.shp')
dstlayer = dstshp.CreateLayer('mylayer',geom_type=ogr.wkbMultiPolygon)
But if you want to use your script (a MultiPolygon is a collection of Polygons):
driver = ogr.GetDriverByName('ESRI Shapefile')
dstshp = driver.CreateDataSource('SomeFilename.shp')
dstlayer = dstshp.CreateLayer('mylayer',geom_type=ogr.wkbPolygon)
for feature1 in layer1:
geom1 = feature1.GetGeometryRef()
attribute1 = feature1.GetField('FieldName1')
for feature2 in layer2:
geom2 = feature2.GetGeometryRef()
attribute2 = feature2.GetField('FieldName2')
# select only the intersections
if geom2.Intersects(geom1):
intersection = geom2.Intersection(geom1)
dstfeature = ogr.Feature(dstlayer.GetLayerDefn())
dstfeature.SetGeometry(intersection)
dstfeature.setField(attribute1)
dstfeature.setField(attribute2)
dstfeature.Destroy()
Don't forget to define the fields before (look at Python GDAL/OGR Cookbook:Vector Layers). And it is much easier with the module Fiona
I thought I'd add my experience here. I've been working with situations where I might need to merge hundreds or thousands of polygons and the Union
solution didn't quite cut it - the time to get the intersection of 1 polygon in 1 layer against a couple of hundred polygons in 1 other layer took ~15 secs. I managed to cut it down to 0.08ms by adding each polygon to a MultiPolygon
and avoiding doing a Union
for every polygon and Simplify
ing the entire geometry at the end of the loop:
union = ogr.Geometry(ogr.wkbMultiPolygon)
for feat in
union.AddGeometry(feat.GetGeometryRef())
union1.Simplify(0)
Actually, I found for my purposes I didn't even need the Simplify
- the intersection geometry seemed to behave as I needed it to without the extra function call.