Raster Calculation in Qgis using Python Script

Question

I want to do raster calculation over raster layer for which i am using python script as:

from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
bohLayer = qgis.utils.iface.activeLayer()
entries = []
# Define band1
boh1 = QgsRasterCalculatorEntry()
boh1.ref = 'boh@1'
boh1.raster = bohLayer
boh1.bandNumber = 1
entries.append( boh1 )
# Process calculation with input extent and resolution
calc = QgsRasterCalculator( 'boh@1/10000', 
 'E:/data/abc.tif', 
 'GTiff',
 bohLayer.extent(), 
 bohLayer.width(), 
 bohLayer.height(), 
 entries )
calc.processCalculation()

Sometimes this script run successfully and gives output but mostly it give error as

AttributeError: 'NoneType' object has no attribute 'extent'

Is there any other way to do raster calculation over raster layer with python script.

• Joseph: You could call the GDAL or SAGA raster calculator from the Processing plugin such as described in this post: Number inputs as raster calculator variables in QGIS Modeler
• Germán Carrillo: Note that you're initializing bohLayer with the QGIS Layer Tree's current layer (aka. active layer or selected layer). I guess the script fails when you have not selected your bohLayer in the Layer Tree. Let me know if that's right to tell you how to avoid relying on a selected layer.
• User18: @German Carrillo : yeah while executing script from c# program it gives an error but while i load layer in QGIS and then execute script it run. It is not finding the current layer . and if i pass the path to the input layer for BohLayer then it do not find its extent . Is any other way to do so by script ?
• User18: @Joseph : With GDAL i am using the following code : gdal_calc.py -A abc.tif --outfile=result.tif --calc="(A)/10000" but it is not working .it shows an error ImportError : numpy.core.multiarray failed to import
• Joseph: @User18 - What QGIS version and Processing version are you using?

Answer 1

I'd like to provide an alternative answer/example using GDAL: raster calculator, updated for QGIS 3. If you look in the 'processing' toolbox then GDAL> Raster Miscellaneous> Raster calculator you'll find the GUI tool. I like to test using the GUI tool first before I write the code.

First, you can call the help docs for this tool using:

processing.algorithmHelp('gdal:rastercalculator')

I find the help docs very straight forward and helpful. Here is an example of what it looks like as a standalone script:

import processing
input_raster = QgsRasterLayer('path/to/your/input/raster', 'raster') 
output_raster = 'path/to/your/output/raster'
#I find it nice to create parameters as a dictionary
parameters = {'INPUT_A' : input_raster,
 'BAND_A' : 1,
 'FORMULA' : '(A > 100)', #your expression here. Mine finds all cells with value > 100. Experiment in the GUI if needed. You can copy and paste exactly the same expression to into your code here
 'OUTPUT' : output_raster}
processing.runAndLoadResults('gdal:rastercalculator', parameters) #feed in the parameters dictionary cleanly as one argument. You can also write the parameters here as individual arguments if you want.

As you might notice, I left out a lot of optional parameters. But you if you read the help docs, just follow along as you would using the GUI tool. Lastly, make sure to capitalize the parameter names as I did in my example.