Raster Calculator in processing scripts crashes using PyQGIS

Question

I am currently writing a processing script for QGIS that goes as follows:

1. Takes in a digital elevation model (DEM) tiff file, a shapefile mask of the study area, and a slope cutoff value in percent
2. Clip Raster by Mask Layer: Clips the DEM to the mask of the study area
3. GDAL Slope: Calculates the slope in percent of the clipped DEM
4. Raster Calculator: Any pixels in the slope raster above the cutoff value are assigned the cutoff value. When using the Raster Calculator interface, this was achieved using this expression (example cutoff value is 20%, this is an alternative to the Con() feature in Arc):

(("raster@1" > 20) = 0) * "raster@1" + (("raster@1" < 20) = 0) * 20

Below is my code where I tried two ways to run Raster Calculator from within the script. The script works fine beforehand, and then falls over. (Note: I am using a simpler expression first just to test things.) The first one crashes, and the second one just produces this error message: An error occurred while performing the calculation.

I did some Googling about using Raster Calculator in PyQGIS for processing scripts and found some similar issues on this forum, but the solutions offered didn't work for my script.

QGIS version: 3.16.2-Hannover

Script below:

# -*- coding: utf-8 -*-
"""
***************************************************************************
* *
* This program is free software; you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation; either version 2 of the License, or *
* (at your option) any later version. *
* *
***************************************************************************
"""
from qgis.PyQt.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
 QgsProcessingException,
 QgsProcessingAlgorithm,
 QgsProcessingParameterRasterLayer,
 QgsProcessingMultiStepFeedback,
 QgsProcessingParameterVectorLayer,
 QgsProcessingParameterRasterDestination,
 QgsProcessingParameterNumber,
 QgsRasterLayer)
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis import processing
import os
class CalcLSfactor(QgsProcessingAlgorithm):
 """
 This script is a test for calculating the LS-factor raster in RUSLE.
 
 All Processing algorithms should extend the QgsProcessingAlgorithm
 class.
 """
 def tr(self, string):
 """
 Returns a translatable string with the self.tr() function.
 """
 return QCoreApplication.translate('Processing', string)
 def createInstance(self):
 return CalcLSfactor()
 def name(self):
 """
 Returns the algorithm name, used for identifying the algorithm. This
 string should be fixed for the algorithm, and must not be localised.
 The name should be unique within each provider. Names should contain
 lowercase alphanumeric characters only and no spaces or other
 formatting characters.
 """
 return 'CalcLSfactor'
 def displayName(self):
 """
 Returns the translated algorithm name, which should be used for any
 user-visible display of the algorithm name.
 """
 return self.tr('Calculate LS-factor')
 def group(self):
 """
 Returns the name of the group this algorithm belongs to. This string
 should be localised.
 """
 return self.tr('RUSLE scripts')
 def groupId(self):
 """
 Returns the unique ID of the group this algorithm belongs to. This
 string should be fixed for the algorithm, and must not be localised.
 The group id should be unique within each provider. Group id should
 contain lowercase alphanumeric characters only and no spaces or other
 formatting characters.
 """
 return 'RUSLEscripts'
 def shortHelpString(self):
 """
 Returns a localised short helper string for the algorithm. This string
 should provide a basic description about what the algorithm does and the
 parameters and outputs associated with it..
 """
 return self.tr("Calculate LS-factor")
 def initAlgorithm(self, config=None):
 """
 Here we define the inputs and output of the algorithm, along
 with some other properties.
 """
 
 self.addParameter(
 QgsProcessingParameterRasterLayer(
 'INPUT',
 'DEM raster',
 defaultValue=None
 )
 )
 
 self.addParameter(
 QgsProcessingParameterVectorLayer(
 'Studyareamask',
 'Study area mask',
 types=[QgsProcessing.TypeVectorPolygon],
 defaultValue=None
 )
 )
 
 self.addParameter(
 QgsProcessingParameterRasterDestination(
 'OUTPUT',
 'Cutoff',
 #createByDefault=True,
 #defaultValue=None
 )
 )
 
 self.addParameter(
 QgsProcessingParameterNumber(
 'CutoffPercent',
 'Cutoff Percent',
 type=QgsProcessingParameterNumber.Double,
 defaultValue=50.0)
 )
 # Intermediate outputs
 ## RB: I think to move this out because it shows up in the interface
 self.addParameter(
 QgsProcessingParameterRasterDestination(
 'CLIP_OUTPUT',
 'Clipped DEM',
 )
 )
 self.addParameter(
 QgsProcessingParameterRasterDestination(
 'SLOPE_OUTPUT',
 'Slope in percent',
 )
 )
 def processAlgorithm(self, parameters, context, model_feedback):
 """
 Here is where the processing itself takes place.
 """
 
 feedback = QgsProcessingMultiStepFeedback(2, model_feedback)
 results = {}
 outputs = {}
 
 # Clip DEM by study area mask
 alg_params = {
 'ALPHA_BAND': False,
 'CROP_TO_CUTLINE': True,
 'DATA_TYPE': 0,
 'EXTRA': '',
 'INPUT': parameters['INPUT'],
 'KEEP_RESOLUTION': True,
 'MASK': parameters['Studyareamask'],
 'MULTITHREADING': False,
 'NODATA': None,
 'OPTIONS': '',
 'SET_RESOLUTION': False,
 'SOURCE_CRS': None,
 'TARGET_CRS': None,
 'X_RESOLUTION': None,
 'Y_RESOLUTION': None,
 'OUTPUT': parameters['CLIP_OUTPUT']
 }
 
 clip_result = processing.run('gdal:cliprasterbymasklayer',
 alg_params, context=context,
 feedback=feedback, is_child_algorithm=True)
 feedback.setCurrentStep(1)
 if feedback.isCanceled():
 return {}
 
 # Calculate slope in percent
 alg_params = {
 'AS_PERCENT': True,
 'BAND': 1,
 'COMPUTE_EDGES': False,
 'EXTRA': '',
 'INPUT': clip_result['OUTPUT'],
 'OPTIONS': '',
 'SCALE': 1,
 'ZEVENBERGEN': False,
 'OUTPUT': parameters['SLOPE_OUTPUT']
 }
 slope_results = processing.run('gdal:slope',
 alg_params, context=context, feedback=feedback,
 is_child_algorithm=True)
 
 feedback.setCurrentStep(2)
 if feedback.isCanceled():
 return {}
 
 # Raster calculator
 # Define the raster layer that's going into the calculator
 slopeLyr = QgsRasterLayer(slope_results['OUTPUT'])
 entries = []
 ras = QgsRasterCalculatorEntry()
 ras.ref = 'slope@1'
 ras.bandNumber = 1
 ras.layer = slopeLyr
 entries.append(ras)
 
 # Define the expression
 ## RB: Doing a test expression first
 #exp = '(("ras@1" > ' + str(parameters['CutoffDegrees']) + ') = 0) * "ras@1" + (("ras@1" < ' + \
 #str(parameters['CutoffDegrees']) + ') = 0) * ' + str(parameters['CutoffDegrees'])
 
 #exp = '((ras@1 > ' + str(parameters['CutoffPercent']) + ') = 0) * ras@1 + ((ras@1 < ' + \
 #str(parameters['CutoffPercent']) + ') = 0) * ' + str(parameters['CutoffPercent'])
 
 exp = 'slope@1 + 25.0'
 ## Debug statements
 model_feedback.pushInfo('This is what exp looks like: ')
 model_feedback.pushInfo(exp)
 #model_feedback.pushInfo('This is what ras.layer looks like: ')
 #model_feedback.pushInfo(ras.layer)
 
 output = self.parameterAsOutputLayer(parameters, 'OUTPUT', context)
 model_feedback.pushInfo('This is what output looks like: ')
 model_feedback.pushInfo(output)
 
 ## Method one: using the processCalculation method
 ## Problem: This one seems to crash when I try it
 '''
 calc = QgsRasterCalculator(exp,
 output,
 'GTiff',
 slopeLyr.extent(),
 slopeLyr.width(),
 slopeLyr.height(),
 entries)
 res = calc.processCalculation(feedback)
 if res == QgsRasterCalculator.ParserError:
 raise QgsProcessingException(self.tr("Error parsing formula"))
 elif res == QgsRasterCalculator.CalculationError:
 raise QgsProcessingException(self.tr("An error occurred while performing the calculation"))
 '''
 
 ## Method two: This one just seems to give "an unexpected error"
 ## Without much other feedback
 alg_params = {
 'EXPRESSION': exp,
 'LAYERS':slope_results['OUTPUT'],
 'CELLSIZE':0,
 'EXTENT':None,
 'CRS':None,
 'OUTPUT':parameters['OUTPUT']
 }
 
 rasCalc_result = processing.run('qgis:rastercalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
 
 if rasCalc_result == QgsRasterCalculator.ParserError:
 raise QgsProcessingException(self.tr("Error parsing formula"))
 elif rasCalc_result == QgsRasterCalculator.CalculationError:
 raise QgsProcessingException(self.tr("An error occurred while performing the calculation"))
 
 results['Cutoff'] = outputs['RasterCalculator']['OUTPUT']
 
 
 return {'OUTPUT': rasCalc_result['OUTPUT'],
 'CLIP_OUTPUT': clip_result['OUTPUT'],
 'SLOPE_OUTPUT': slope_results['OUTPUT']}

• benaviru: Yes, I've edited that statement to reflect that. Although I am only using a simple exp = 'slope@1 + 25.0' in the actual code as a test yet it still seems to fall over.
• Khaled: For me, any processing algorithm with starting with "gdal:" crashes, processing.run('gdal:ANYTHING). Other with native:STHG, grass: STHG, perfect. If I check gdal if it is installed in cmd, it is.

Answer 1

I get a similar crash but it works when I remove the use of QgsProcessingMultiStepFeedback in the process.