I am working with OGR in Python and noticed that the C libraries have useful Intersection() and Clip() functions as part of the Layer class. Is there any way to get at these functions in Python? I know that these functions exist on the Geometry level but I'm looking specifically at the Layer. Thanks!

  • Unfortunately the problem with Intersection at the geometry level, as far as I understand it, is that you lose any reference to the associated attributes. The geometry based intersection method produces a collection of geometries with no original feature level attributes. Anyone know a solution for that?
    – user14032
    Jan 7, 2013 at 21:51
  • At a higher level, handling multiple records for each attribute requires aggregation. I can think of a good PostGIS solution with SQL.
    – Mike T
    Jan 7, 2013 at 22:09

2 Answers 2


Yes. They are exposed in the bindings:

>>> from osgeo import ogr

>>> help(ogr.Layer.Intersection)
Help on method Intersection in module osgeo.ogr:

Intersection(self, *args, **kwargs) unbound osgeo.ogr.Layer method
    Intersection(self, Layer method_layer, Layer result_layer, char options = None, 
    GDALProgressFunc callback = None, void callback_data = None) -> OGRErr

>>> help(ogr.Layer.Clip)
Help on method Clip in module osgeo.ogr:

Clip(self, *args, **kwargs) unbound osgeo.ogr.Layer method
Clip(self, Layer method_layer, Layer result_layer, char options = None, 
    GDALProgressFunc callback = None, void callback_data = None) -> OGRErr

I am guessing you need GEOS support built into GDAL and confirmed:




  • Wich versin are you using, here (version 1.9) the binding is not available.
    – Pablo
    Jan 30, 2013 at 13:31
  • I'm on the trunk, 1.10, but those have been around for a while. If GEOS support isn't enabled, they may not be built in the bindings. I am not sure though.
    – user10353
    Jan 30, 2013 at 16:45
  • Hey thanks for this, I'm going to try and figure out my setup and see what I'm missing. At the moment I also can't see the bindings.
    – Doug
    Feb 14, 2013 at 19:47
  • I am using gdal 1.9.2 but can't see these bindings.
    – user17574
    Apr 25, 2013 at 20:47
  • 1
    BTW, ogr.Layer.Clip was 4x faster than Shapely and 2x faster than Shapely with a prepared_geometry pre-check.
    – Pablo
    Dec 6, 2013 at 12:24

A layer is composed of one or several geometries. For the intersection of layers, you must iterate through each layer geometries. With shapely it is easy, example with two shapefiles:

from osgeo import ogr
from shapely.wkb import loads
from shapely.geometry import *
# first layer, a polygon shapefile
first = Polygon()
# open shapefile
source1 = ogr.Open("test1.shp")
layer1 = source1.GetLayer()
# combination of all the geometries of the layer in a single shapely object
for element in layer1:
   geom = loads(element.GetGeometryRef().ExportToWkb())
   first = first.union(geom)    
# second layer, a polygon shapefile
two = Polygon()
source2 = ogr.Open("test2.shp")
layer2 = source2.GetLayer()
for element in layer2:
   geom = loads(element.GetGeometryRef().ExportToWkb())
   two = two.union(geom)

# intersection between the two layers 
print first.intersection(two).wkt

It is possible to use the same type of treatment for Clip(). Another solution is provided by Creating a little clipbox for your GIS projects in Python

  • Thanks, this is actually what I ended up doing. The only problem however is that the polygon I want to take the intersection over is a global shapefile. Basically I have an area of interest and want to Clip the global polygon. Using Shapely this way ends up taking quite a long time compared to using OGR.
    – Doug
    Feb 14, 2013 at 19:44

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.