2

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)

4
  • Probably, it is because of paths. Take a look to this answer:link
    – Gakovî
    Jan 17, 2021 at 8:52
  • Thank you for your suggestions. I tried to make both there modifies, avoiding the space characters and adding a 'r' in front of the path, but the situation doesn't change. any other suggestions?
    – Rapu93
    Jan 17, 2021 at 9:18
  • Welcome to GIS SE! We're a little different from other sites; this isn't a discussion forum but a Q&A site. Please check out our short tour to learn about our focussed Q&A format.
    – PolyGeo
    Jan 18, 2021 at 23:15
  • See the edited answer with added new evapotranspiration string. Jan 19, 2021 at 9:47

2 Answers 2

4

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 )'
0
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.

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.