1

I would like to perform this raster calculation where when ras@1 >0 gives me 1 as pixel value of the output raster, while when ras@1 <0 the values are given by the other formula ras@1/(1-2.71828^(ras@2)), but it doesn't work because it gives me in any case the result of the other formula, even when ras@1 is higher than zero. Can someone give me some suggestion?

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

inputpath = '/Users/macbook/Desktop/TESI/PROGETTO/sigm/sigmath90.tif'
inputpath2 = '/Users/macbook/Desktop/TESI/PROGETTO/lambdamin/lambminth90.tif'
outputpath = '/Users/macbook/Desktop/TESI/PROGETTO/alfamax/alphamaxth90.tif'

entries = []
inputrasterfile = QgsRasterLayer(inputpath)  
ras = QgsRasterCalculatorEntry()
ras.ref = 'ras@1'
ras.raster = inputrasterfile
ras.bandNumber = 1
entries.append( ras )

inputrasterfile2 = QgsRasterLayer(inputpath2)  
ras2 = QgsRasterCalculatorEntry()
ras2.ref = 'ras@2'
ras2.raster = inputrasterfile2
ras2.bandNumber = 1
entries.append( ras2 )


alpha_max='(ras@1 < 1)* ras@1 / (1-2.71828^(ras@2)) + (ras@1 >>1)*1'

new_path = outputpath

alfamax= QgsRasterCalculator(alpha_max, new_path, 'GTiff', inputrasterfile.extent(), inputrasterfile.width(), inputrasterfile.height(), entries )
alfamax.processCalculation()
iface.addRasterLayer(new_path)
0

1 Answer 1

1

This statement:

I would like to perform this raster calculation where when ras@1 >0 gives me 1 as pixel value of the output raster, while when ras@1 <0 the values are given by the other formula ras@1/(1-2.71828^(ras@2))

is not compatible with your alpha_max code formula:

alpha_max='(ras@1 < 1)* ras@1 / (1-2.71828^(ras@2)) + (ras@1 >>1)*1'

With adequate layer references, your formula should be changed in this way:

alpha_max = '(' + ras1.ref + ' > 0)*1  + ( ' + ras1.ref + ' <= 0)*( ' + ras1.ref + ' / (1 - 2.71828^' + ras2.ref + '))'

So, following version of your code works as expected (I used my own paths and commented yours).

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

#inputpath = '/Users/macbook/Desktop/TESI/PROGETTO/sigm/sigmath90.tif'
inputpath = '/home/zeito/Desktop/PROGETTO/sigm/sigmahar90.tif'

#inputpath2 = '/Users/macbook/Desktop/TESI/PROGETTO/lambdamin/lambminth90.tif'
inputpath2 = '/home/zeito/Desktop/PROGETTO/lambdamin/lambminhar90.tif'

#outputpath = '/Users/macbook/Desktop/TESI/PROGETTO/alfamax/alphamaxth90.tif'
outputpath = '/home/zeito/Desktop/PROGETTO/alfamax/alphamaxhar90.tif'

entries = []
inputrasterfile = QgsRasterLayer(inputpath)  
ras1 = QgsRasterCalculatorEntry()
ras1.ref = 'ras1@1'
ras1.raster = inputrasterfile
ras1.bandNumber = 1
entries.append( ras1 )

inputrasterfile2 = QgsRasterLayer(inputpath2)  
ras2 = QgsRasterCalculatorEntry()
ras2.ref = 'ras2@1'
ras2.raster = inputrasterfile2
ras2.bandNumber = 1
entries.append( ras2 )

alpha_max = '(' + ras1.ref + ' > 0)*1  + ( ' + ras1.ref + ' <= 0)*( ' + ras1.ref + ' / (1 - 2.71828^' + ras2.ref + '))'
print(alpha_max)

new_path = outputpath

alfamax= QgsRasterCalculator(alpha_max, 
                             new_path, 
                             'GTiff', 
                             inputrasterfile.extent(), 
                             inputrasterfile.width(), 
                             inputrasterfile.height(), 
                             entries )

alfamax.processCalculation()
iface.addRasterLayer(new_path)

I ran above code in QGIS Python Console and pixel values of resulting loaded layer were corroborated with original raster values. Its pixel values are as expected (according to conditional used formulas).

enter image description here

Not the answer you're looking for? Browse other questions tagged or ask your own question.