7

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
8

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

  • Thanks for the comments and for assistance in improving the script. Much appreciated and it works. I will go through the entire OGR Cookbook. I found it to be helpful. – user2757128 Jan 16 '14 at 8:57
  • I have edited my own with your suggestions. The inclusion of 'if geom2.Intersects(geom1):' before doing the actual instersection has reduced the computational time by a big margin (more than half). – user2757128 Jan 16 '14 at 9:06
  • I don't follow the purpose of doing a Union first. What do you gain? – Phlucious Mar 16 '18 at 17:57
  • I'm not sure how the union works in the for loop, I keep getting a NoneType error. Having created union1=ogr.Geometry(3), the line union1 = union1.Union(geom) keeps resulting in an AttributeError: 'NoneType' object has no attribute 'Union' – grego Oct 16 at 19:49
2

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 Simplifying 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.

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.