Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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!

share|improve this question
    
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 '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 7 '13 at 22:09
add comment

2 Answers

up vote 5 down vote accepted

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

share|improve this answer
    
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. –  kyle 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. –  user17574 Apr 25 '13 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 '13 at 12:24
show 2 more comments

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

share|improve this answer
    
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.