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!
6
-
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
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. – user10353 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
-
1BTW, 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
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
-
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
