QGIS raster calculator using Python console

Question

I want to calculate sum values of many rasters (from input folder) in QGIS Python using raster calculator. When I run script, QGIS crashes:

Windows fatal exception: access violation

import os
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer
work_folder = 'C:/QGIS/input_rasters'
rasters = [file for file in os.listdir(work_folder) if file.endswith('.nc')]
entries = []
i = 0
for raster in rasters:
 i = i + 1
 raster_path = os.path.join(work_folder, raster)
 input_raster = QgsRasterLayer(raster_path, raster)
 entry = QgsRasterCalculatorEntry()
 entry.ref = f"'ras@{i}'"
 entry.raster = input_raster
 entry.bandNumber = 1
 entries.append(entry)
 
expression = " + ".join(entry.ref for entry in entries)
output = 'C:/QGIS/output.tif'
calc = QgsRasterCalculator(expression, output, 'GTiff', input_raster.extent(), input_raster.width(), input_raster.height(), entries)
calc.processCalculation()

When I rewrite to two rasters, code works, but I need to do it for many rasters.

import os
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer
work_folder = 'C:/QGIS/input_rasters'
rasters = [file for file in os.listdir(work_folder) if file.endswith('.nc')]
entries = []
raster_path = os.path.join(work_folder, 'test1.nc')
input_raster = QgsRasterLayer(raster_path, raster)
entry = QgsRasterCalculatorEntry()
entry.ref = 'ras@1'
entry.raster = input_raster
entry.bandNumber = 1
entries.append(entry)
raster2_path = os.path.join(work_folder, 'test2.nc')
input2_raster = QgsRasterLayer(raster2_path, raster)
entry = QgsRasterCalculatorEntry()
entry.ref = 'ras@2'
entry.raster = input2_raster
entry.bandNumber = 1
entries.append(entry)
output = 'C:/QGIS/output/output.tif'
calc = QgsRasterCalculator('ras@1 + ras@2', output, 'GTiff', input_raster.extent(), input_raster.width(), input_raster.height(), entries)
calc.processCalculation()

• xunilk: Your code has two problems mainly. The first one is layers list. It is really a list of names; not a list of raster layers. There are problems referencing them properly with your approach in the loop. The second one is related to layers references. Please, see my answer. I hope it's not too late.

Answer 1

Your code has two problems mainly. The first one is layers list. It is really a list of names; not a list of raster layers. There are problems referencing them properly with your approach in the loop. The second one is related to layers references. Your way is for multi bands raster; not for single band.

So, following code fixes all issues. I used three *.nc raster created by me with random values between 0 and 10.

work_folder = '/home/zeito/Desktop/test4'
layers = [ QgsRasterLayer(os.path.join(work_folder, file), file[:-3]) 
 for file in os.listdir(work_folder) if file.endswith('.nc')]
entries = []
for i, raster in enumerate(layers):
 i = i + 1
 print(i, raster)
 entry = QgsRasterCalculatorEntry()
 entry.ref = 'ras' + str(i) + '@1'
 entry.raster = raster
 entry.bandNumber = 1
 entries.append(entry)
expression = " + ".join(entry.ref for entry in entries)
print(expression)
output = '/home/zeito/Desktop/test4/output.tif'
calc = QgsRasterCalculator(expression, 
 output, 
 'GTiff', 
 layers[0].extent(), 
 layers[0].width(), 
 layers[0].height(), 
 entries)
calc.processCalculation()

After running above code in Python Console of QGIS, I loaded all involved layers and the sum result can be observed in following picture. It was as expected by using Value tool plugin.

enter image description here