Determinando elevaciones máximas y mínimas (estadística zonal) de un ráster utilizando rpy objetos en ambiente de PyQGIS

En el GDAL/OGR Cookbook existe código para acometer la estadística zonal de un ráster empleando las librerías GDAL/OGR de Python. Sin embargo, la extensión de este contrasta con lo conciso que puede llegar a ser el mismo procedimiento con lenguaje R.

Para usar R con Python es necesario instalar en Debian los paquetes python-rpy y python-rpy2. Después de ello, utilizando el vectorial y el ráster dem de la imagen siguiente:

rcode1

desarrollé el código a continuación, en ambiente de PyQGIS, que permite obtener los valores máximo y mínimo (estadística zonal) de un raster con base en un vectorial tipo polígono e introducir en la tabla atributiva de este los valores obtenidos en sus campos respectivos.

from PyQt4.QtCore import *

import rpy2.robjects as robjects
import os

r = robjects.r

mapcanvas = iface.mapCanvas()
layers = mapcanvas.layers()

vp = layers[0].dataProvider()
vpath_tmp = vp.dataSourceUri().split('|')
vpath = vpath_tmp[0]

rp = layers[1].dataProvider()
rpath = rp.dataSourceUri()

vroot, vfile = os.path.split(vpath)
rroot, rfile = os.path.split(rpath)

r.setwd(vroot)

r.library('raster')

arg1 = 's <- shapefile(' + "\"" + vfile + "\"" + ')'
arg2 = 'd <- raster(' + "\"" + rfile + "\"" + ')'

r(arg1)
r(arg2)

print "Wait..."

demmin = r('s$demmin <- extract(d, s, fun = min, na.rm = TRUE)')
demmax = r('s$demmax <- extract(d, s, fun = max, na.rm = TRUE)')

fields = [ QgsField('min', QVariant.Double), QgsField('max', QVariant.Double) ]
vp.addAttributes( fields )
layers[0].updateFields()

idx1 = layers[0].fieldNameIndex('min')
idx2 = layers[0].fieldNameIndex('max')

n = len(demmin)

for i in range(n):
    new_values = { idx1 : demmin[i], idx2 : demmax[i] }
    vp.changeAttributeValues( {i:new_values} )

print "Done!"

Después de ejecutado el código anterior en la Python Console de QGIS se obtuvo el resultado esperado; tal como se expresa en la tabla de atributos del vectorial tipo polígono de la imagen siguiente:

rcode2

Esta entrada fue publicada en Lenguaje R, 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