QGIS Python raster calculator

Question

I would like to run this code that makes calculation at every loop between different input rasters and Itot that is a raster that in the loop remains always the same.
The code doesn't give me any error but it doesn't work because at the end the writing:'invalid layer gdal provider cannot open the gdal database' appear.
I guess the problem could be in how I treat Itot in the calculation, can someone help me?
Here is the code:
import qgis
import gdal 
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer

p=[0.84, 0.83, 1.03, 1.11, 1.24, 1.25, 1.27, 1.18, 1.04, 0.96, 0.83, 0.81]
Itot = QgsRasterLayer('/Users/macbook/Desktop/TESI/CAPITANATA PROGETTO/Itot.tif')

suffix_output=['th_jan', 'th_feb', 'th_mar', 'th_apr', 'th_may', 'th_jun', 'th_jul', 'th_aug', 'th_sep', 'th_oct', 'th_nov' ,'th_dec']
suffix_input=['Tjan', 'Tfeb', 'Tmar', 'Tapr', 'Tmay', 'Tjun', 'Tjul', 'Taug', 'Tsep', 'Toct', 'Tnov', 'Tdec']

outputpath = r'/Users/macbook/Desktop/TESI/PROGETTO/PET_thornthwaite/'
inputpath = r'/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'

th = str(p[i]) +'* 16* (10* ras@1 / Itot)^(0.5 + 0.016 * Itot)'

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

Pet.processCalculation()
  iface.addRasterLayer(new_path)

xunilk · Answer

You need an Itot reference in the calculations. Following code works as expected. I marked your paths for using my own paths and corroborate adequate running of code. Observe that I only used two test raster layers. In your case, you need to change 2 for 12 in the loop, mark or delete my paths and unmark your paths.
import qgis
import gdal 
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer

p=[0.84, 0.83, 1.03, 1.11, 1.24, 1.25, 1.27, 1.18, 1.04, 0.96, 0.83, 0.81]
#Itot = QgsRasterLayer('/Users/macbook/Desktop/TESI/CAPITANATA PROGETTO/Itot.tif')
Itot = QgsRasterLayer('/home/zeito/Desktop/I/Itot.tif', 'Itot')

entries = []

itot = QgsRasterCalculatorEntry()
itot.ref = 'Itot@1'
itot.raster = Itot
itot.bandNumber = 1

entries.append( itot )

suffix_output=['th_jan', 'th_feb', 'th_mar', 'th_apr', 'th_may', 'th_jun', 'th_jul', 'th_aug', 'th_sep', 'th_oct', 'th_nov' ,'th_dec']
suffix_input=['Tjan', 'Tfeb', 'Tmar', 'Tapr', 'Tmay', 'Tjun', 'Tjul', 'Taug', 'Tsep', 'Toct', 'Tnov', 'Tdec']

#outputpath = r'/Users/macbook/Desktop/TESI/PROGETTO/PET_thornthwaite/'
outputpath = r'/home/zeito/Desktop/PROGETTO/'

#inputpath = r'/Users/macbook/Desktop/TESI/PROGETTO/raster_T_average/'
inputpath = r'/home/zeito/Desktop/raster_T_average/'

layers_path = []

layers = [ QgsRasterLayer(inputpath + suffix_input[i] + '.tif', suffix_input[i] ) for i in range(len(suffix_input)) ]

raster_ref = [ '{}@1'.format(suffix) for suffix in suffix_input ]

for i in range(2): #in your case 12 instead 2
    
    layer = QgsRasterCalculatorEntry()
    layer.ref = raster_ref[i]
    layer.raster = layers[i]
    layer.bandNumber = 1
    entries.append( layer )

new_path = outputpath + suffix_output[i] + '.tif'
    
    th = str(p[i]) +' * 16* (10* ' + layer.ref  + ' / ' + itot.ref + ' )^(0.5 + 0.016 * ' + itot.ref + ' ) '
    print(th)
    print(new_path)

cal = QgsRasterCalculator(th, 
                              new_path, 
                              'GTiff', 
                              layers[i].extent(), 
                              layers[i].width(), 
                              layers[i].height(), 
                              entries )

cal.processCalculation()
    iface.addRasterLayer(new_path)

After running above code in Python Console of QGIS, loaded layers were valid; as it can be observed in following image.

QGIS Python raster calculator

One Answer

Add your own answers!

Ask a Question