Cómo producir polígonos con anillos centrales a una cierta distancia de los bordes

Con Base en una pregunta de gis.stackexchange.com, se desarrolló el código siguiente para producir polígonos con anillos centrales a una cierta distancia de los bordes (10 km en este caso). El por qué de este tipo de polígonos surge de la necesidad de una estadística zonal que excluyera los píxeles más allá de la distancia establecida.

En el código primero se obtiene una memory layer que corresponde a un buffer “negativo” (-10000 m) para cada rasgo de la capa activa; la cual no se añade al registro. Este buffer se utiliza posteriormente como capa de diferencia de la capa activa la cual se reserva como una segunda memory layer que esta vez si se añade al registro. La geometría resultante será un polígono con anillo central sólo para aquellas situaciones en que tenga sentido la operación de geoproceso (mayores de 10 km2; aunque condicionados también por su forma).

layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures() ]
 
#determine buffer -10000 m
buffer = [ feat.geometry().buffer(-10000,-1) for feat in feats ] 

#Extract CRS from route
CRS = layer.crs().postgisSrid()
 
URI = "Polygon?crs=epsg:"+str(CRS)+"&field=id:integer""&index=yes"
 
#Create polygon layer for buffer
mem_layer = QgsVectorLayer(URI,
                           "buffer",
                           "memory")
 

#Prepare mem_layer for editing
mem_layer.startEditing()

n = len(feats)
 
#Set feature for buffer
feats2 = [ QgsFeature() for i in range(n) ]

for i in range(n): 
    #Set geometry for buffer
    feats2[i].setGeometry(buffer[i])
 
    #set attributes values for buffer
    feats2[i].setAttributes([i])
  
    mem_layer.addFeature(feats2[i], True)
 
#stop editing and save changes
mem_layer.commitChanges()

feats3 = [ feat3 for feat3 in mem_layer.getFeatures() ]

difference = [ feats[i].geometry().difference(feats3[i].geometry()) for i in range(n) ] 

#Create polygon layer for difference
mem_layer2 = QgsVectorLayer(URI,
                           "difference",
                           "memory")

#Prepare mem_layer for editing
mem_layer2.startEditing()

#Set feature for difference
feats2 = [ QgsFeature() for i in range(n) ] 

for i in range(n):
    #Set geometry for difference
    feats2[i].setGeometry(difference[i])
 
    #set attributes values for difference
    feats2[i].setAttributes([i])
  
    mem_layer2.addFeature(feats2[i], True)
 
#stop editing and save changes
mem_layer2.commitChanges()

#add Map Layer to Registry
QgsMapLayerRegistry.instance().addMapLayer(mem_layer2)

La ejecución del código anterior en la Python Console de QGIS, para el shapefile de la imagen siguiente, permite evaluar su ejecución.

holes2

El resultado es el siguiente:

holes

Los polígonos sin anillo no cumplen con los requerimientos de distancia.

No obstante, si en el código cambiamos 10000 por 3000, el resultado es muy diferente. Ahora todos los polígonos si presentan anillo central.

holes3

El código funciona como se espera.

Esta entrada fue publicada en PyQGIS, QGIS, SIG, Software Libre. Guarda el enlace permanente.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s