¿Hay una manera para llegar a C ' s OGRLayer.Intersection() o OGRLayer.Clip() en Python?

Question

Estoy trabajando con OGR en Python y notó que las bibliotecas de C tienen funciones Intersection() y Clip() útiles como parte de la clase de la capa. ¿Hay alguna manera de obtener en estas funciones en Python? Sé que estas funciones existen en el nivel de geometría pero estoy buscando específicamente en la capa. ¡Gracias!

Answer 1

Sí. Se exponen en los enlaces:

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

Supongo que necesita apoyo GEOS en GDAL y confirmado:

http://GDAL.org/OGR/classOGRLayer.html#ac189f54996c2d6fd769889ec99e0f48a

y

http://GDAL.org/OGR/classOGRLayer.html#a56d7ee3b2020e53c730d67ee4f1e2fb6

Answer 2

Una capa se compone de una o varias geometrías. Para la intersección de las capas, debe recorrer cada geometrías de la capa. Con bien es fácil, ejemplo con archivos de dos formas:

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

Es posible utilizar el mismo tipo de tratamiento para Clip(). Otra solución se proporciona mediante la creación de un poquito clipbox para sus proyectos GIS en Python