2

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)
3
  • After a successful manual run in QGIS check Processing -> History. You will see the exact Python command that worked there. Then compare it to yours and look for potential differences / mistakes. Jan 29, 2021 at 15:56
  • I can't find it, there are other processes but not this
    – Rapu93
    Jan 29, 2021 at 16:23
  • I fixed U90 ref and a new version of your code works as you expected.
    – xunilk
    Jan 29, 2021 at 16:45

1 Answer 1

3

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.

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.