6

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!

Doug
  • 617
  • 1
  • 6
  • 10
  • 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? –  Jan 07 '13 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 07 '13 at 22:09

2 Answers2

7

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:

http://gdal.org/ogr/classOGRLayer.html#ac189f54996c2d6fd769889ec99e0f48a

and

http://gdal.org/ogr/classOGRLayer.html#a56d7ee3b2020e53c730d67ee4f1e2fb6

  • Wich versin are you using, here (version 1.9) the binding is not available. – Pablo Jan 30 '13 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. –  Jan 30 '13 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 '13 at 19:47
  • I am using gdal 1.9.2 but can't see these bindings. –  Apr 25 '13 at 20:47
  • Are you using the builds from gisinternals/sdk? –  Apr 26 '13 at 16:31
  • Building GDAL 1.10 on Ubuntu finally gave me access to those functions. – Pablo Dec 06 '13 at 12:19
  • 1
    BTW, ogr.Layer.Clip was 4x faster than Shapely and 2x faster than Shapely with a prepared_geometry pre-check. – Pablo Dec 06 '13 at 12:24
5

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

gene
  • 54,868
  • 3
  • 110
  • 187
  • 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 '13 at 19:44