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")