0

I would like to make the sum of 12 different rasters to get one raster as output. I wrote the following code but it doesn't work (or better, it works, but at the end it doesn't give me any layer, but the writing "Invalid layer, Gdal provider cannot open Gdal dataset"), any suggestions to how modifying it?

import qgis
import gdal 
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer

suffix_input=['I_jan', 'I_feb', 'I_mar', 'I_apr', 'I_may', 'I_jun', 'I_jul', 'I_aug', 'I_sep', 'I_oct', 'I_nov' ,'I_dec']
inputpath = '/Users/macbook/Desktop/TESI/PROGETTO/I/'

list=[0,1,2,3,4,5,6,7,8,9,10,11]
sum=0
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 )

  sum = 'sum + ras@1'

  calc = QgsRasterCalculator( sum, '/Users/macbook/Desktop/TESI/CAPITANATA PROGETTO/Itot.tif', 'GTiff', inputrasterfile.extent(), inputrasterfile.width(), inputrasterfile.height(), entries )
  print(calc.processCalculation())
  iface.addRasterLayer('/Users/macbook/Desktop/TESI/PROGETTO/Itot.tif')
1
  • Also, I just noticed you declare entries = [] inside a loop. So that means the list is emptied at the start of each loop iteration. Commented Jan 21, 2021 at 12:24

1 Answer 1

2
  1. If you want to iterate a list while having access to iteration number then use

for i, item in enumerate(list):

If you don't need the counter, then just do it like for item in list:

  1. I would propose that you don't sum the rasters in a loop but instead calculate the sum with a single calculation after the loop like this:

     expr = ''
     rasters = []
     for i, suffix in enumerate(suffix_input):
         inputrasterfile = QgsRasterLayer(inputpath + suffix + ".tif", suffix)
         rasters.append(inputrasterfile)
    
         expr = expr + '"{}@1"+'.format(suffix)
    
     expr = expr.rstrip('+')
    
     alg_params = {
         'CELLSIZE': 0,
         'CRS': None,
         'EXPRESSION': expr,
         'EXTENT': rasters[0].extent(),
         'LAYERS': rasters,
         'OUTPUT': '/Users/macbook/Desktop/TESI/PROGETTO/Itot.tif')
     }
    
     result = processing.run('qgis:rastercalculator', alg_params)
     iface.addRasterLayer('/Users/macbook/Desktop/TESI/PROGETTO/Itot.tif')
    
8
  • Thank you, now the code works, but at the end the raster of output Itot has all the pixels with value 1.79769e+308, can you say me where is the mistake in doing the calculations? I'm sure that input rasters have right values
    – Rapu93
    Commented Jan 21, 2021 at 13:02
  • Try to run some checks. For example: see if inputrasterfile.isValid() returns true Commented Jan 21, 2021 at 13:09
  • Maybe you need to set CRS or something else. Since I don't have your data it is hard to tell. Commented Jan 21, 2021 at 13:10
  • Maybe you want to sum a different band? Right now @1 is summed. Just tossing suggestions Commented Jan 21, 2021 at 13:11
  • @LeonPowałka Your code works perfectly. I created some raster with same prefixes and some nodata values and raster sum is also obtained where all pixel are different to nodata values. I think issues can be related to bad alignment of original raster.
    – xunilk
    Commented Jan 21, 2021 at 14:35

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.