NDVI con la clase QgsRasterCalculator en PyQGIS

El empleo de cocientes o índices para identificar masas vegetales tiene su base en que una masa vegetal en óptimas condiciones posee una firma espectral que se caracteriza, en especial, por un claro contraste entre la banda del rojo (0.6 a 0.7 µm) y el infrarrojo cercano (0.7 a 1.1 µm). Esto se debe a que la mayor parte de la radiación solar recibida por la planta en el visible es absorbida por los pigmentos de las hojas y éstos apenas afectan a la radiación recibida en el infrarrojo cercano. En el caso del NDVI (Indice de Vegetación de Diferencias Normalizadas; por sus siglas en inglés), éste se determina dividiendo la diferencia de la banda infrarroja (C4) menos la roja (C3) entre la suma de ambas (C3+C4); expresadas en términos de reflectividad.

Con base en lo anterior, elaboré una pequeña función que ejecutada desde la Python Console toma, en ese orden, las bandas de reflectividad b3 y b4 y mediante los métodos de las clases QgsRasterCalculator y QgsRasterCalculatorEntry determina el NDVI. La función es la siguiente:

from qgis.core import *
from qgis.gui import *
from qgis.utils import  iface
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

def ndvi():

    mc = iface.mapCanvas()

    Layers = mc.layers()

    entries = []

    for lyr in Layers:
        # Define bands
        layer = QgsRasterCalculatorEntry()
        layer.ref = lyr.name() +'@1'
        layer.raster = lyr
        layer.bandNumber = 1
        entries.append( layer )

    root = 'c:/pyqgis_data/'
    
    expression = '(  ' + entries[1].ref + '   -  '  + entries[0].ref + '  ) \
    /  (  ' + entries[1].ref+'  +  ' + entries[0].ref+'  )'
    print expression
    
    calc = QgsRasterCalculator(expression, 
                                root + 'ndvi.tif', 
                                'GTiff', 
                                Layers[0].extent(), 
                                Layers[0].width(), 
                                Layers[0].height(), 
                                entries )
                                  
    calc.processCalculation()

Este código lo guardé como raster.py (por si deseo en el futuro añadir más funciones) en mi ruta de scripts registrada como PYTHONPATH en QGIS.

Para ejecutarla hay que cargar primero en la Map View de QGIS las bandas b3 y b4 (en la Map Legend deben aparecer en este orden). Luego, desde la Python Console, hay que hacer lo siguiente:

>>>from raster import ndvi
>>>ndvi()

cuya imagen después de la ejecución luce como sigue:

ndvi

El print dentro de la función para su visualización en la Python Console permite verificar que las referencias a las bandas para el cálculo del NDVI son las esperadas. En la ruta c:\pyqgis_data\ está la imagen ndvi.tif. Al cargarla a QGIS y aplicarle un estilo previamente creado para este fin (ndvi.qml) el resultado es el siguiente:

raster_calculator4

En la ventana del plugin Value Tool, ubicado el cursor en un punto arbitrario sobre el Map Canvas, aparecen unos valores con los cuales podrá comprobar que el script funciona de manera adecuada píxel a píxel aplicando manualmente la fórmula del NDVI con una calculadora manual.

Esta entrada fue publicada en Código Python, PyQGIS, SIG, Software Libre. Guarda el enlace permanente.

2 respuestas a NDVI con la clase QgsRasterCalculator en PyQGIS

  1. Pingback: NDVI con la clase QgsRasterCalculator en PyQGIS...

  2. Pingback: [ BLOG ] – El Blog de José Guerrero | All Around GIS

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