¿Borrar rasgos con GDAL/OGR de Python?

Question

Tengo una característica de polígono redondo grande (5km de diámetro) y una característica de polígono redondo pequeño (3km de diámetro) que me gustaría recortar con Python GDAL/OGR para obtener el "donut". Una vez hecho esto, aplicaré el mismo proceso para 25.000 características, por eso necesito hacerlo con GDAL/OGR y no editarlo manualmente en QGIS/ArcGIS. Parece que el Funcionamiento de la pinza funciona a nivel de capas. Por lo tanto, lo que hice fue crear una capa para mi círculo de 5 km, una capa para mi círculo de 3 km y una capa de salida requerida por el prototipo de la función.

Aparentemente la operación está terminando correctamente, porque obtengo un archivo .shp con 6KB dentro, pero cuando abro este shapefile ya sea con QGIS o ArcGIS está completamente en blanco. Pensé que era la proyección, pero el archivo .prj dice que está configurado en EPSG:28992 (Amersfoort_RD_New), que es lo que necesito.

El código que escribí es el siguiente:

outDriver = ogr.GetDriverByName("ESRI Shapefile")    
outDataSource = outDriver.CreateDataSource(pathout)
outLayer = outDataSource.CreateLayer("Big", srs = outSpatialRef, geom_type=ogr.wkbPolygon)

outDriver2 = ogr.GetDriverByName("ESRI Shapefile")    
outDataSource2 = outDriver2.CreateDataSource(pathout)
outLayer2 = outDataSource2.CreateLayer("Small", srs = outSpatialRef, geom_type=ogr.wkbPolygon)

outDriver3 = ogr.GetDriverByName("ESRI Shapefile")    
outDataSource3 = outDriver3.CreateDataSource(pathout)
outLayer3 = outDataSource3.CreateLayer("Clip", srs = outSpatialRef, geom_type=ogr.wkbPolygon)

outLayer.CreateFeature(mybigfeature)
outLayer2.CreateFeature(mysmallfeature)
something = outLayer.Clip(outLayer2, outLayer3)
print something # returns zero

Entonces, no entiendo como la operación parece terminar correctamente, escribe 6KB en mi disco y luego no veo nada en QGIS. He adjuntado una imagen, mostrando las características, pero no el polígono recortado.

Supongo que falta algo: ¿el lavado de la memoria?

insertar una geometría para la salida (¿cómo si no hay retorno de Clip)?

¿cerrar los archivos?

Espero que alguien pueda detectar el problema.

---

Siguiendo el consejo de Benjamin, instalé las bibliotecas Shapely y Fiona para resolver este problema. Shapely se utiliza "para la manipulación y el análisis de objetos geométricos planares" y Fiona hace la lectura y escritura de shapefiles de una manera pitónica y muy ordenada.

Así que después de jugar con estas bibliotecas, se me ocurrió el siguiente código:

from shapely.geometry import shape, mapping
import fiona

pathb = r'PATH_TO_BIG_ROUND_FEATURES\Buffer_5km.shp'
paths = r'PATH_TO_SMALL_ROUND_FEATURES\Buffer_2km.shp'
pathout = r'PATH_TO_OUTPUT_FOLDER\clip.shp'

big = fiona.open(pathb)
small = fiona.open(paths)

polbig = big.next()
polsmall = small.next()

a = shape(polbig['geometry'])
b = shape(polsmall['geometry'])
c = a.difference(b)

source_schema = { 'geometry': 'Polygon', 'properties': {'id': 'int'}}
source_crs = {'no_defs': True, 'ellps': 'WGS84', 'datum': 'WGS84', 'proj': 'longlat'}
source_driver = "ESRI Shapefile"

with fiona.open(pathout, 'w', 'ESRI Shapefile', schema=source_schema, crs=source_crs) as out:
    out.write({
        'geometry': mapping(c),
        'properties': {'id': 123},
    })

Y esto produce la siguiente salida:

Que es exactamente lo que necesito.

Answer 1

Ya que estás usando python también podrías echar un vistazo a Shapely . Es una biblioteca para la manipulación y el análisis de objetos geométricos en el plano cartesiano.

Para su problema se procedería como

• Cargar ambas geometrías en objetos Shapely
• en el objeto del que se quiere sustraer el otro llame a object.difference(other)
• guardar el resultado en una nueva variable

Para ver un ejemplo, consulte http://toblerity.org/shapely/manual.html#object.difference

Personalmente te recomendaría usar una librería que tenga una interfaz de python, sobre todo porque cuando usas un IDE para Desarrollo tienes cosas como autocompletar y comprobación de tipo.

Answer 2

Yo utilizaría el dialecto GDAL SQLite. Aquí viene un ejemplo:

Los dos polígonos de la imagen son como WKT:

POLÍGONO (( -113 39, -113 42, -110 42, -110 39, -113 39 ))

POLÍGONO (( -112 40, -112 41, -111 41, -111 40, -112 40 ))

Puede hacer un agujero en el polígono 1 con el polígono 2 utilizando la función de Spatialite ST_Difference https://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html

Demostración con ogrinfo:

ogrinfo -ro -dialect sqlite -sql "select ST_Difference(ST_GeomFromText('POLYGON (( -113 39, -113 42, -110 42, -110 39, -113 39 ))'),ST_GeomFromText('POLYGON (( -112 40, -112 41, -111 41, -111 40, -112 40 ))'))" multipolygon.json

INFO: Open of `multipolygon.json'
      using driver `GeoJSON' successful.

Layer name: SELECT
Geometry: Unknown (any)
Feature Count: 1
Extent: (-113.000000, 39.000000) - (-110.000000, 42.000000)
Layer SRS WKT:
(unknown)
Geometry Column = ST_Difference(ST_GeomFromText('POLYGON (( -113 39, -113 42, -1
10 42, -110 39, -113 39 ))'),ST_GeomFromText('POLYGON (( -112 40, -112 41, -111
41, -111 40, -112 40 ))'))
OGRFeature(SELECT):0
  POLYGON ((-113 39,-113 42,-110 42,-110 39,-113 39),(-112 40,-111 40,-111 41,-1
12 41,-112 40))

El resultado es

POLÍGONO ((-113 39,-113 42,-110 42,-110 39,-113 39),(-112 40,-111 40,-111 41,-112 41,-112 40))

Y se ve así

Lo que hay que hacer es construir el SQL con Python y ejecutarlo.

Answer 3

Como mencioné en mi comentario, podrías dividir los polígonos interiores y los más grandes en dos capas separadas y utilizar Herramienta de diferencias (QGIS) o Herramienta de borrado (ArcGIS) para que se ejecute en un solo proceso. Si quieres, estas dos herramientas podrían configurarse para ejecutarse también en un script de python.