Estadísticas de un raster “Single Band” usando PyQGIS

PyQGIS tiene más de 500 clases y, en consecuencia, miles de métodos de clase listados en la documentación. Por tanto, no es difícil “perderse” explorando entre ellos; sobre todo cuando los objetos de diferentes clases tienen métodos coincidentes. En este post se van determinar las estadísticas (y otra información relevante) de un raster mediante PyQGIS. Para ello, vamos a acceder a los objetos de diferentes clases utilizando los métodos adecuados.

En primer lugar, se va a utilizar QgsInterface para acceder a los objetos de la clase QgsRasterLayer. El método ‘renderer’ de esta última nos permite acceder a los objetos de la clase QgsMultiBandColorRenderer o QgsSingleBandGrayRenderer; dependiendo de si se tiene respectivamente un raster ‘multibanda’ o ‘single band’. Estas clases tienen, a su vez, métodos coincidentes denominados ‘bandStatistics’ que pudieran hacer pensar que son los adecuados. Sin embargo, no lo son. El método ‘bandStatistics’ que se requiere para acceder a ellas es a través del ‘dataProvider’ de QgsRasterLayer.

Una vez determinado el método a emplear, se obtiene la documentación de ‘bandStatistics’ a ver que argumentos son necesarios. En la Python Console, con help, se puede tener una pista:

>>>help(QgsRasterDataProvider.bandStatistics)
Help on built-in function bandStatistics:

bandStatistics(...)
    QgsRasterInterface.bandStatistics(int, int theStats=QgsRasterBandStats.All, QgsRectangle theExtent=QgsRectangle(), int theSampleSize=0) -> QgsRasterBandStats

El primer ‘int’ es el número de banda. Yo voy a probar con una single band y, en este caso, sé que es 1 (que es lo que voy a usar en el ejemplo). Ahora, si se está automatizando un proceso éstas se pueden determinar con los métodos de la clase QgsMultiBandColorRenderer o QgsSingleBandGrayRenderer. El segundo int es de un ‘Public Types enum’ que para no complicarme mucho asumo que es QgsRasterBandStats.All (es decir, todas las posibles de la banda). El tercer argumento es la extensión sobre la cual se va a computar las estadísticas (objeto de la clase QgsRectangle que si se omite se asume que las estadísticas abarcan todo el raster). El cuarto y último argumento corresponde al número aproximado de celdas en la muestra (si es cero se tomará todo el ráster). Por tanto, es posible prescindir de los dos últimos argumentos si se desea que las estadísticas abarquen todo el ráster. Para más detalles se puede consultar la documentación de QgsRasterInterface.

El objeto que se obtiene al aplicar el método ‘bandStatistics’ pertenece a la clase QgsRasterBandStats; cuya salida con dir de sus métodos de clase es:

['All', 'Max', 'Mean', 'Min', 'None', 'Range', 'Stats', 'StdDev', 
'Sum', 'SumOfSquares', '__class__', '__delattr__', '__dict__', 
'__doc__', '__format__', '__getattribute__', '__hash__', '__init__', 
'__module__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', 
'__setattr__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', 
'bandNumber', 'contains', 'elementCount', 'extent', 'height', 
'maximumValue', 'mean', 'minimumValue', 'range', 'statsGathered', 
'stdDev', 'sum', 'sumOfSquares', 'width']

Para probar los métodos de la clase anterior se usó el siguiente código (que al inicio también lo combino con métodos de otras clases) como un script que ejecuté en el editor de la Python Console:

layer=iface.activeLayer()

print layer

ext = layer.extent()

print ext

print ext.asPolygon()

#Raster resolution
print "rasterUnitsPerPixelX = ", layer.rasterUnitsPerPixelX()
print "rasterUnitsPerPixelY = ", layer.rasterUnitsPerPixelY()

provider=layer.dataProvider()

print provider

print "UsesBand = ", layer.renderer().usesBands()
print "type = ", layer.renderer().type()

#ver =provider.hasStatistics(1,QgsRasterBandStats.All,ext,0)
#Uso solo dos argumentos en lugar de 4 y el codigo funciona!
ver =provider.hasStatistics(1,QgsRasterBandStats.All)

if ver is not False:
    #Aqui se usan todos los argumentos
    stats=provider.bandStatistics(1,QgsRasterBandStats.All,ext,0)

    print "bandNumber = ", stats.bandNumber

    print "mean = ", stats.mean

    print "stdDev = ", stats.stdDev

    print "sum = ", stats.sum

    print "sumOfSquares = ", stats.sumOfSquares
    
    print "minimumValue = ", stats.minimumValue
    
    print "maximumValue = ", stats.maximumValue
    
    print "range = ", stats.range

La salida de correr el código anterior como script es la siguiente:

Python 2.7.4 (default, Apr  6 2013, 19:54:46) [MSC v.1500 32 bit (Intel)] on zeito
## Type help(iface) for more info and list of methods.
>>>execfile(u'C:/Users/zeito/pyqgis_scripts/provider.py'.encode('mbcs'))
<qgis.core.QgsRasterLayer object at 0x207F20C0>
<qgis.core.QgsRectangle object at 0x24F52F18>
354971.34886022 4414903.32231663, 354971.34886022 4473428.40239009, 479272.40388350 4473428.40239009, 479272.40388350 4414903.32231663, 354971.34886022 4414903.32231663
rasterUnitsPerPixelX =  73.9887232281
rasterUnitsPerPixelY =  73.9887232281
<qgis.core.QgsRasterDataProvider object at 0x24F40F18>
UsesBand =  [1]
type =  singlebandgray
bandNumber =  1
mean =  1824.71800614
stdDev =  373.273662546
sum =  2424831264.0
sumOfSquares =  1.85156999563e+11
minimumValue =  1323.0
maximumValue =  3556.0
range =  2233.0

para el siguiente raster DEM:

raster

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

3 respuestas a Estadísticas de un raster “Single Band” usando PyQGIS

  1. Pingback: Estadísticas de un raster “Single ...

  2. Pingback: Cómo leer los valores de un ráster usando el módulo GDAL con PyQGIS | El Blog de José Guerrero

  3. Pingback: Funciones con PyQGIS: archivos ráster | El Blog de José Guerrero

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