QGIS Python script with raster calculations

Question

I have a problem with this code, that works, because it gives me the output rasters that I want, but the problem is that these outputs show wrong values because they are all 1.79769e+308.

I think the formula is right, because when I print it is what I would expect, and also the output path is good because all the output rasters are saved in the right folder.

Maybe I make some mistakes in importing the right input rasters at every loop? Because if I make calculations (by commenting the others) only with ras@1 or only with ras@2 or only with ras@3 the code works perfectly, so I guess the problem could be when I import the three rasters together at every loop.

Here is the code:

import qgis
import gdal 
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer
Rg=[6.4, 8.6, 11.4, 14.3, 16.4, 17.3, 16.7, 15.2, 12.5, 9.6, 7.0, 5.7]
d=[31,28,31,30,31,30,31,31,30,31,30,31]
suffix_output=['Harg_jan', 'Harg_feb', 'Harg_mar', 'Harg_apr', 'Harg_may', 'Harg_jun', 'Harg_jul', 'Harg_aug', 'Harg_sep', 'Harg_oct', 'Harg_nov' ,'Harg_dec']
suffix_input=['Tjan', 'Tfeb', 'Tmar', 'Tapr', 'Tmay', 'Tjun', 'Tjul', 'Taug', 'Tsep', 'Toct', 'Tnov', 'Tdec']
suffix_input2=['Tmaxjan', 'Tmaxfeb', 'Tmaxmar', 'Tmaxapr', 'Tmaxmay', 'Tmaxjun', 'Tmaxjul', 'Tmaxaug', 'Tmaxsep', 'Tmaxoct', 'Tmaxnov', 'Tmaxdec']
suffix_input3=['Tminjan', 'Tminfeb', 'Tminmar', 'Tminapr', 'Tminmay', 'Tminjun', 'Tminjul', 'Tminaug', 'Tminsep', 'Tminoct', 'Tminnov', 'Tmindec']
outputpath = r'/Users/macbook/Desktop/TESI/PROGETTO/PET_Hargreaves/'
inputpath = r'/Users/macbook/Desktop/TESI/PROGETTO/raster_T_average/'
inputpath2 = r'/Users/macbook/Desktop/TESI/PROGETTO/raster_T_max/'
inputpath3 = r'/Users/macbook/Desktop/TESI/PROGETTO/raster_T_min/'
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 )
 inputrasterfile2 = QgsRasterLayer(inputpath2 + suffix_input2[i] + ".tif")
 entries = []
 ras2 = QgsRasterCalculatorEntry()
 ras2.ref = 'ras@2'
 ras2.raster = inputrasterfile2
 ras2.bandNumber = 1
 entries.append( ras2 )
 inputrasterfile3 = QgsRasterLayer(inputpath3 + suffix_input3[i] + ".tif")
 entries = []
 ras3 = QgsRasterCalculatorEntry()
 ras3.ref = 'ras@3'
 ras3.raster = inputrasterfile3
 ras3.bandNumber = 1
 entries.append( ras3 )
 new_path = outputpath + suffix_output[i]+ '.tif'
 hargreaves = str(d[i]) + '* (0.0023 *' + str(Rg[i]) +'* (17.8 + ras@1) * (ras@2 - ras@3)^0.5)'
 print(hargreaves)
 Pet = QgsRasterCalculator(hargreaves, new_path, 'GTiff', inputrasterfile.extent(), inputrasterfile.width(), inputrasterfile.height(), entries )
 Pet.processCalculation()
 iface.addRasterLayer(new_path)

• NW_Photo_Laureate: Question Cesare. So the code works now, but the values of the output rasters are basically infinity? Are there cells in the raster that look good? Can you load the image in QGIS or plot it somehow while masking those values out. Set them to NaN for viewing purposes and see if a coherent image emerges. I've had inexplicable issues like this with gdal and numpy raster calculations. That large value may be the limit of a 64 bit or 32 bit signed floating point.
• Rapu93: Yes exactly. If I set them to NaN I see all the image black, so every pixel takes the value NaN. While if I make the same calculation only for one iteration in the raster calculator (not in the Python console) it gives me right values. But if I do the same in the console it doesn't happen.
• NW_Photo_Laureate: Just to clarify, when you say "one iteration" do you mean the for loop or from inputrasterfile to entries.append( ras )? My only observation is that you initiate a new list - entry = [] for each of the three inputraster<#>. So when you run QgsRasterCalculator(hargreaves, new_path, 'GTiff', inputrasterfile.extent(), inputrasterfile.width(), inputrasterfile.height(), entries ) then entries is just ras3 reference layer. Did you intend to reference all 3 - ras ras2 ras3? If so, remove entry = [] from ras2 and ras3 code block.