Calcular las estadísticas zonales de una parte de los polígonos

Question

Estoy trabajando con un shapefile que contiene los polígonos de distrito del mapa político de la India y un archivo raster que contiene los datos de elevación del mismo país. Mi objetivo es calcular la altitud media y la varianza en un radio de 10 km desde las fronteras de cada distrito. Utilizando Estadísticas zonales He podido calcular las estadísticas que busco a nivel de distrito, el paso siguiente que necesito es descartar los píxeles situados en la parte central.

He intentado deshacerme de las partes internas de los distritos con añadir anillos pero el usuario debe dibujar los polígonos internos por sí mismo. Esto no es lo que busco: Quiero tener un polígono interno con exactamente los mismos ángulos que el original pero con lados más pequeños, calculados automáticamente para toda la India.

Cómo puedo crear, alternativamente, una máscara en el shapefile para la zona que no quiero (o que quiero, la que más se ajuste al caso) tener en cuenta y luego utilizarla para calcular la media y la varianza con el método Estadísticas zonales ¿un plugin?

La respuesta ideal me proporcionaría un código PyQGIS pero está bien también usando iconos en QGIS.

Answer 1

He probado el siguiente código con un shapefile poligonal (como capa activa) donde, primero, se produce un buffer "negativo" (-10000 m) como capa de memoria. Después, esta capa de memoria se utiliza como capa de diferencia para la capa activa para producir un polígono sin la parte interna (se obtiene un anillo central).

layer = iface.activeLayer()

feat = layer.getFeatures().next()

#determining buffer -10000 m
buffer = feat.geometry().buffer(-10000,-1)

#Extract CRS from layer
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()

#Set feature for buffer
feat2 = QgsFeature()

#Set geometry for buffer
feat2.setGeometry(buffer)

#set attributes values for buffer
feat2.setAttributes([1])

mem_layer.addFeature(feat2, True)

#stop editing and save changes
mem_layer.commitChanges()

#preparing layer for difference
feat3 = mem_layer.getFeatures().next()

difference = feat.geometry().difference(feat3.geometry())

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

#Prepare mem_layer2 for editing
mem_layer2.startEditing()

#Set feature for difference
feat2 = QgsFeature()

#Set geometry for difference
feat2.setGeometry(difference)

#set attributes values for difference
feat2.setAttributes([1])

mem_layer2.addFeature(feat2, True)

#stop editing and save changes
mem_layer2.commitChanges()

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

La siguiente imagen se obtuvo tras ejecutar el código en la consola Python de QGIS:

Funciona como se esperaba.

Nota de edición:

La siguiente versión del script funciona para todas las características de una capa vectorial poligonal.

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 layer
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 = []

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

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

Se ejecutó con una capa vectorial en la que, en algunos rasgos, la distancia a las fronteras es inferior a 10 km; como puede observarse en la siguiente imagen: