How to Calculate the Vegetation Index NDVI from Sentinel 2 Imagery with PyQGIS

Vegetation indexes are calculated from the plant radiation in certain ranges of the visible and infrared spectrum. There are many indexes based on different band combination formulas, one of the most common indexes is the Normalized Deviation Vegetation Index (NDVI) because it was of the first vegetation index and it can be applied to images from current and old satellites. This tutorial shows the complete procedure to represent in QGIS the red and near infrared (NIR) images from a clipped Sentinel 2 image with PyQGIS and then calculate the NDVI using the processing package.

Tutorial

 

Input data

You can download the input data from this link.

Code

This is the complete code in Python:

import os
from osgeo import gdal,gdal_array

os.chdir("/Users/saulmontoya/Documents/Ih_VegetationIndexNDVISentinel2PyQGIS/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)