3

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)

1 Answer 1

2

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.

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.