Raster calculation with Python QGIS

Question

I have this code in which I would like to perform a for loop in which every time I get one raster as input (the one called L, different for every month) and divide it always by the same raster (the one called U).
The code seems to work, but it always gives me as result an invalid raster since all the pixels have value 3.40282e+38, BUT if I do the same manually with the Raster Calculator (so not with a python script) it gives me the right result.
Can someone tell me where is the mistake in my code?
import qgis
import gdal 
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer

U90 = QgsRasterLayer('/Users/macbook/Desktop/TESI/PROGETTO/rasterUfinali/U_1990.tif', 'U90')

entries = []

u = QgsRasterCalculatorEntry()
u.ref = 'U@1'
u.raster = U90
u.bandNumber = 1

entries.append( u )

inputpath= r'/Users/macbook/Desktop/TESI/PROGETTO/L/'

suffix_input=['Ljan', 'Lfeb', 'Lmar', 'Lapr', 'Lmay', 'Ljun', 'Ljul', 'Laug', 'Lsep', 'Loct', 'Lnov', 'Ldec']

outputpath= r'/Users/macbook/Desktop/TESI/PROGETTO/lamb/'

suffix_output=['lambdajan','lambdafeb','lambdamar','lambdaapr','lambdamay','lambdajun','lambdajul','lambdaaug','lambdasep','lambdaoct','lambdanov','lambdadec']

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(12):

layer = QgsRasterCalculatorEntry()
    layer.ref = raster_ref[i]
    layer.raster = layers[i]
    layer.bandNumber = 1
    entries.append( layer )

new_path = outputpath + suffix_output[i] + '.tif'

l = layer.ref  +'/'+ u.ref  
    print(l)

lamb = QgsRasterCalculator(l, new_path, 'GTiff', layers[i].extent(), layers[i].width(), layers[i].height(), entries )

lamb.processCalculation()
    iface.addRasterLayer(new_path)

xunilk · Accepted Answer

I fixed U90 ref and following version of your code works as you expected. I used my own paths for verification purpose. Delete them an uncomment yours. It is also necessary to change 2 for 12 in corresponding loop.
import qgis
import gdal 
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer

#U90 = QgsRasterLayer('/Users/macbook/Desktop/TESI/PROGETTO/rasterUfinali/U_1990.tif', 'U90')
U90 = QgsRasterLayer('/home/zeito/Desktop/PROGETTO/rasterUfinali/U_1990.tif', 'U90')

entries = []

u90 = QgsRasterCalculatorEntry()
u90.ref = 'U90@1'
u90.raster = U90
u90.bandNumber = 1

entries.append( u90 )

#inputpath= r'/Users/macbook/Desktop/TESI/PROGETTO/L/'
inputpath= r'/home/zeito/Desktop/L/'

suffix_input=['Ljan', 'Lfeb', 'Lmar', 'Lapr', 'Lmay', 'Ljun', 'Ljul', 'Laug', 'Lsep', 'Loct', 'Lnov', 'Ldec']

#outputpath= r'/Users/macbook/Desktop/TESI/PROGETTO/lamb/'
outputpath= r'/home/zeito/Desktop/PROGETTO/lamb/'

suffix_output=['lambdajan','lambdafeb','lambdamar','lambdaapr','lambdamay','lambdajun','lambdajul','lambdaaug','lambdasep','lambdaoct','lambdanov','lambdadec']

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 change 2 for 12

layer = QgsRasterCalculatorEntry()
    layer.ref = raster_ref[i]
    layer.raster = layers[i]
    layer.bandNumber = 1
    entries.append( layer )

new_path = outputpath + suffix_output[i] + '.tif'

l = layer.ref  + ' / ' + u90.ref  
    print(l)
    print(new_path)

lamb = QgsRasterCalculator(l, 
                               new_path, 
                               'GTiff', 
                               layers[i].extent(), 
                               layers[i].width(), 
                               layers[i].height(), 
                               entries )

lamb.processCalculation()
    iface.addRasterLayer(new_path)

After running above code, valid raster layers were loaded in Map Canvas; as it can be observed in following image. I corroborated that division was correct pixel by pixel in each raster.

Raster calculation with Python QGIS

One Answer

Add your own answers!

Ask a Question