Calculating NDVI using QGIS Python Console

Question

I'm going to calculate NDVI using band 4 and band 3 of Landsat 5.

The code I wrote:

import os, re
import processing
from osgeo import gdal
input_B4 = r'C:\Users\test\LT05_L2SP_114034_19850820_20200918_02_T1_SR_B4.TIF'
input_B3 = r'C:\Users\test\LT05_L2SP_114034_19850820_20200918_02_T1_SR_B3.TIF'
outfile = r'C:\Users\test_raster_cal\test1.tif'
processing.run("gdal:rastercalculator", 
{'INPUT_A': input_B4,
 'BAND_A':1,
 'INPUT_B': input_B3, 
 'BAND_B':1,
 'FORMULA':'(A-B)/(A+B)',
 'NO_DATA':None,
 'RTYPE':5,
 'OPTIONS':'',
 'EXTRA':'',
 'OUTPUT':outfile
 })

In the picture below, test_cal is made using the QGIS Raster Calculator Tool Test1 is a code that I wrote.

enter image description here

enter image description here

enter image description here

enter image description here

[QGIS Raster Calculator Tool]

enter image description here

[My code]

enter image description here

When using the QGIS Tool, the value is between -1 and 1, and the NDVI made with the code I wrote is about 0 to 65.

How do I modify the code to calculate the NDVI properly?

• ChloeG: @BERA For Landsat 5, B3 and B4 seem to be the correct bands.
• ChloeG: What I'm seeing is that input_B4 and input_B3 become input4 and input3 in your code.
• xunilk: I think that there is nothing wrong with your code and the issue is relative to a bug in the visualization of limits in QGIS Map legend. Please, see my answer.

Answer 1

I think that there is nothing wrong with your code and the issue is relative to a bug in the calculation and visualization of limits in QGIS Map legend (issues with manipulation of no data values in the whole images because there is not problems in reduced areas). I tried to reproduce your results but, I got bad limits only with raster calculator instead your code.

Different results probably could be related to my QGIS 3 version (3.24.3-Tisler) on Linux. I got SR_B3 and SR_B4 bands with Google Earth Engine by using LANDSAT_PRODUCT_ID obtained from your bands (LT05_L2SP_114034_19850820_20200918_02_T1). These bands belongs to an image situated in South Korea (PATH: 114, ROW: 34, UTM_ZONE: 52).

In following picture can be observed that ndvi_sk (raster calculator) has bad limits compared with ndvi_test obtained with my version of your code.

enter image description here

The used pyqgis code looks as follows:

import os, re
import processing
from osgeo import gdal
processing.runAndLoadResults("gdal:rastercalculator", 
{'INPUT_A': '/home/zeito/pyqgis_data/SR_B4_sk.tif',
 'BAND_A':1,
 'INPUT_B': '/home/zeito/pyqgis_data/SR_B3_sk.tif', 
 'BAND_B':1,
 'FORMULA':'(A-B)/(A+B)',
 'NO_DATA':None,
 'RTYPE':5,
 'OPTIONS':'',
 'EXTRA':'',
 'OUTPUT':'ndvi_test'
 })

However, I explored tens of pixels by using Value Tool QGIS plugin, as it can be observed in the example of following picture, and all compared values (ndvi_sk and ndvi_test) were practically identical. So, probably, there is nothing wrong with your code.

enter image description here