QGIS Python script with Raster Calculator

Question

I have a large set of radar .tif files with several bands, and I would like to have some Python scripts to perform some basic operations that I would usually perform using the Raster Calculator in QGIS 3.22.6, namely:

• Create a new set of .tif files with only one of the bands of the input files, with the name equal to the original raster with a suffix (like here)
• Create a new set of .tif files in which the output is a single band resulting from applying (a+b+c+d)/4 to the first 4 bands of the input files.
• Create a new set of .tif files in which the output is a single band resulting from applying d/b to the 2nd and 4th input file bands.

I am currently trying to do this based on this similar question, but I don't know how to use other layers than the first one. My attempt did not produce the expected results:

from qgis.analysis import QgsRasterCalculatorEntry, QgsRasterCalculator
from qgis.core import QgsProject
lddLrs = [tree_layer.layer() for tree_layer in QgsProject.instance().layerTreeRoot().findLayers()]
path = "path/to/output"
for lyr in lddLrs:
 entries = []
 ras = QgsRasterCalculatorEntry()
 ras.ref = 'ras@1'
 ras.raster = lyr
 ras.bandNumber = 1
 entries.append(ras)
 calc = QgsRasterCalculator('("infra-red@1"+"infra-red@2"+"infra-red@3"+"infra-red@4")/4', path + lyr.name() + "_suffix.tif", 'GTiff', lyr.extent(), lyr.width(), lyr.height(), entries)
calc.processCalculation()

This other question seemed to have a similar problem, but there doesn't seem to be an answer to it yet.

Answer 1

I think it's easier to use GDAL raster calculator. I execute the tool manually once in QGIS, press Ctrl+Alt+H and copy the command. Then write the loop.

To create the output names you can use layer.source() and os.path.join with os.path.basename.

For first band, and mean:

import os
output_folder = r'/home/bera/Desktop/GIStest/rasters/'
for layer in QgsProject.instance().mapLayers().values(): #For all layers in the map
 print(layer)
 #########Extract first band
 #Create the output name
 #layer.source() 
 #'/home/bera/Desktop/GIStest/4bandraster.tif'
 #os.path.basename(layer.source()) #
 #'4bandraster.tif'
 output1 = os.path.join(output_folder, os.path.basename(layer.source()).replace('.tif', '_firstbandonly.tif'))
 #'/home/bera/Desktop/GIStest/rasters/4bandraster_suffix.tif'
 processing.run("gdal:rastercalculator", {'INPUT_A':layer,'BAND_A':1,'INPUT_B':None,
 'BAND_B':None,'INPUT_C':None,'BAND_C':None,'INPUT_D':None,'BAND_D':None,
 'INPUT_E':None,'BAND_E':None,'INPUT_F':None,'BAND_F':None,'FORMULA':'A',
 'NO_DATA':None,'RTYPE':5,'OPTIONS':'','EXTRA':'','OUTPUT':output1})
 
 ########Mean of bands 1,2,3,4:
 output2 = os.path.join(output_folder, os.path.basename(layer.source()).replace('.tif', '_mean1234.tif'))
 processing.run("gdal:rastercalculator", {'INPUT_A':layer,'BAND_A':1,'INPUT_B':layer,'BAND_B':2,'INPUT_C':layer,'BAND_C':3,'INPUT_D':layer,'BAND_D':4,'INPUT_E':None,'BAND_E':None,'INPUT_F':None,'BAND_F':None,'FORMULA':'(A+B+C+D)/4','NO_DATA':None,'RTYPE':5,'OPTIONS':'','EXTRA':'','OUTPUT':output2})
 
print("Done")

enter image description here

• Nick Carraway: Is it possible to apply raster calculator to files that are not opened in QGIS? Loading hundreds of file to QGIS takes too much memory.