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.

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.