Writing PyQGIS script with operations on multiple set of rasters

Question

I would like to perform a loop that takes from a given folder a set of rasters and for every raster makes the following operation, that in the specific case is the calculation of the monthly evapotranspiration, that consist in multiplying every raster of monthly average temperature for a given coefficient (that is specified in the two arrays C and Rg).

The code doesn't give any error but it doesn't work, because at the end of the loop appears what follow: Invalid layer, Gdal cannot open the dataset.

Does anyone have any suggestion?

import glob, qgis
import gdal 
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer
C=[0.4, 0.37, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4,0.4]
Rg=[109, 158, 261, 376, 467, 506, 518, 448, 321, 207, 117, 94]
suffix_output=['PET_jan', 'PET_feb', 'PET_mar', 'PET_apr', 'PET_may', 'PET_jun', 'PET_jul', 'PET_aug', 'PET_sep', 'PET_oct', 'PET_nov' ,'PET_dec']
suffix_input=['Tjan', 'Tfeb', 'Tmar', 'Tapr', 'Tmay', 'Tjun', 'Tjul', 'Taug', 'Tsep', 'Toct', 'Tnov', 'Tdec']

outputpath = '/Users/macbook/Desktop/TESI/PROGETTO/'
inputpath = '/Users/macbook/Desktop/TESI/PROGETTO/raster T average/'

list=[0,1,2,3,4,5,6,7,8,9,10,11]
for i in list:
  inputrasterfile = QgsRasterLayer(inputpath + suffix_input[i] + ".tif")

  entries = []
  ras = QgsRasterCalculatorEntry()
  ras.ref = 'ras@1'
  ras.raster = inputrasterfile
  ras.bandNumber = 1
  entries.append( ras )

  new_path = outputpath + suffix_output[i]+ '.tif'

  evapotranspiration = str(C[i]) +'* ras@1 ' + '*(str(Rg[i]) + 50)' + '/(ras@1 + 15)' 

  Pet = QgsRasterCalculator(evapotranspiration, new_path, 'GTiff', inputrasterfile.extent(), inputrasterfile.width(), inputrasterfile.height(), entries )

  Pet.processCalculation()
  iface.addRasterLayer(new_path)

---

Just to understand, why these statement doesn't work? I tried to follow the former suggestion, but the software says that I can't multiply sequence by non-int of type 'str'. I just would like to understand the rule behind this calculation...

evapotranspiration = (str(k[i])  + '* ( 0.46  * ras@1 + 8')   * str(p[i]) + '-2)' *str(d[i])

Here I want to write this formula: [k(i)p(i)(0,46*raster+8)-2]*d(i)

Answer 1

The line

evapotranspiration = str(C[i]) +'* ras@1 ' + '*(str(Rg[i]) + 50)' + '/(ras@1 + 15)'

creates and invalid formula for raster calculation.

In fact, at the first iteration, evapotranspiration is equal to 0.4* ras@1 *(str(Rg[i]) + 50)/(ras@1 + 15) that is an invalid formula.

str(Rg[i]) should be outside the single quote.

Maybe the correct line is:

evapotranspiration = str(C[i]) +'* ras@1 *(' + str(Rg[i]) + ' + 50)' + '/(ras@1 + 15)'

So, at the first iteration, evapotranspiration will be equal to 0.4* ras@1 *(109 + 50)/(ras@1 + 15) that is a valid formula.

It could be also better to avoid spaces in the file path (inputpath = '/Users/macbook/Desktop/TESI/PROGETTO/raster T average/' contains spaces in the path).

EDIT:

For [k(i)p(i)(0,46*raster+8)-2]*d(i):

evapotranspiration = str(d[i]) + ' * ( ' + str(k[i]) + ' * ' + str(p[i]) + ' * ( 0.46 * ras@1 + 8 ) - 2 )'

Answer 2

For overcoming these issues, when you are coding with QgsRasterCalculator class, you need to be pretty sure that each formula is right. So, you have to print them before whatever calculation. On the other hand, you can use os.path python module for handling file and folder names adequately.

For testing my approach, in my Linux Debian system, I created similar paths to yours for my input ('/home/zeito/Desktop/raster_T_average') and output ('/home/zeito/Desktop/PROGETTO/') raster layers. Afterward, I created only 2 random raster layers (20x20 and values between 1 and 10) for simulating Tjan y Tfeb and put them in input folder.

As your formulas are wrong, I used following expression in each iteration and previously printed them for corroborating their validity:

evapotranspiration = str(C[i]) + '*' +  layer.ref + '*' + (str(Rg[i] + 50)) + '/(' + layer.ref + ' + 15)' 

Complete code looks as follows and produces two layers without any error:

from qgis.core import *
from qgis.gui import *
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from math import *
import os

inputpath = '/home/zeito/Desktop/raster_T_average'

if os.path.isdir(inputpath):
    files = os.listdir(inputpath)
else:
    print("directory not exists")

layers_paths = []
layers_refs = []

C = [0.4, 0.37, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4, 0.4,0.4]
Rg = [109, 158, 261, 376, 467, 506, 518, 448, 321, 207, 117, 94]

for file in files:
    if file[-4:] == '.tif':
        print(file)
        path = os.path.join(inputpath, file)
        layers_paths.append(path)
        layers_refs.append(file[:-4])

layers = [ QgsRasterLayer(layers_paths[i], layers_refs[i]) for i in range(len(layers_paths)) ]
 
entries = []

for i, layer in enumerate(layers):

    # Define band1
    layer = QgsRasterCalculatorEntry()
    layer.ref = layers_refs[i] + '@1'
    layer.raster = layers[i]
    layer.bandNumber = 1
    entries.append( layer )
    
    evapotranspiration = str(C[i]) + '*' +  layer.ref + '*' + (str(Rg[i] + 50)) + '/(' + layer.ref + ' + 15)' 
    
    print(evapotranspiration)
    
    calc = QgsRasterCalculator(evapotranspiration, 
                                '/home/zeito/Desktop/PROGETTO/outputfile' + str(i+1) + '.tif', 
                                'GTiff', 
                                layers[i].extent(), 
                                layers[i].width(), 
                                layers[i].height(), 
                                entries )
                                 
    calc.processCalculation()

After running above code in Python Console of QGIS, I loaded original and produced layers and they look as follows. Printed formulas (in red rectangle) can be used for corroborating, by using Value Tool QGIS plugin, if values were produced as expected.