OGR – OGR Layer Intersection in Python

intersectionogrpython

I am using what I would describe as the most primitive method of doing intersection with OGR. The short script below describes how I do it. Is there a best way of doing this? 

   from osgeo import ogr, osr
   shp1 = ogr.Open(file1)
   shp2 = ogr.Open(file2)

   SpatialRef = osr.SpatialReference()
   SpatialRef.SetWellKnownGeogCS('WGS84')
   # Create dst file here 
   dstshp = ogr.CreateDataSource('SomeFilename.shp')
   dstlayer = dstshp.CreateLayer()
   # define its attribute fields for dstlayer and create them

   layer1 = shp1.GetLayer(0)
   layer2 = shp2.GetLayer(0)

   for feature1 in layer1:
      geom1 = feature1.GetGeometryRef()
      attribute1 = feature1.GetField('FieldName1')
      for feature2 in layer2:
         geom2 = feature2.GetGeometryRef()
         attribute2 = feature2.GetField('FieldName2')  
         intersection = geom2.intersection(geom1)
         dstfeature = ogr.Feature(dstlayer.GetLayerDefn())
         dstfeature.SetGeometry(intersection)
         dstfeature.setField(attribute1)
         dstfeature.setField(attribute2)
         dstfeature.Destroy() # and other features must be destroyed too

Best Answer

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

Related Question