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!
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:
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