Calculo de Indice de Vegetación NDVI de Imágenes Sentinel 2 con Python en QGIS (PyQGIS)

Los índices de vegetación se calculan a partir de las radiancias de las plantas en ciertos rangos del espectro visible e invisible. Curiosamente la vegetación tiene una mayor radiancia en el rojo y el infrarojo que en el azul. Existen varios índices de vegetación en base de diferentes fórmulas de combinación de bandas, dentro de estos índices el más conocido es el NDVI, ya que fue uno de los primeros en formularse y porque puede ser aplicado a una serie de satélites nuevos y antiguos. Este tutorial muestra el procedimiento completo para representar las bandas Rojo e Infrarojo de imágenes de Sentinel 2 con PyQGIS y posteriormente calcular el NDVI utilizando herramientas del complemento "Processing" dentro de Python.

Tutorial

Código

Este es el código completo en Python para el tutorial.

import os
from osgeo import gdal,gdal_array

os.chdir("C:\Users\Saul\Documents\Ih_HowtoCalculateVegetationIndexNDVIfromSentinel2withPyQGIS\Sentinel2Clip")
print(os.listdir(os.getcwd()))

NIR = iface.addRasterLayer('clip_RT_S2A_OPER_MSI_L1C_TL_MTI__20160506T214824_A004555_T18LTM_B08.tif','NIR')
RED = iface.addRasterLayer('clip_RT_S2A_OPER_MSI_L1C_TL_MTI__20160506T214824_A004555_T18LTM_B04.tif','RED')

import processing

NIR = processing.getObjectFromName("NIR")
RED = processing.getObjectFromName("RED")

NDVI_syntax = '(A-B)/(A+B)'

outputs_GDALOGRRASTERCALCULATOR_1=processing.runalg('gdalogr:rastercalculator', 
                                                    NIR,         #INPUT_A <ParameterRaster>
                                                    '1',         #BAND_A <ParameterString>
                                                    RED,         #INPUT_B <ParameterRaster>
                                                    '1',         #BAND_B <ParameterString>
                                                    None,        #INPUT_C <ParameterRaster>
                                                    '1',         #BAND_C <ParameterString>
                                                    None,        #INPUT_D <ParameterRaster>
                                                    '1',         #BAND_D <ParameterString>
                                                    None,        #INPUT_E <ParameterRaster>
                                                    '1',         #BAND_E <ParameterString>
                                                    None,        #INPUT_F <ParameterRaster>
                                                    '1',         #BAND_F <ParameterString>
                                                    NDVI_syntax, #FORMULA <ParameterString>
                                                    '',          #NO_DATA <ParameterString>
                                                    5,           #RTYPE <ParameterSelection>
                                                    '0',         #EXTRA <ParameterString>
                                                    None)        #OUTPUT <OutputRaster>
                                                    
NDVI = QgsRasterLayer(outputs_GDALOGRRASTERCALCULATOR_1['OUTPUT'],'NDVI')

QgsMapLayerRegistry.instance().addMapLayer(NDVI)