1

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']}
1
  • 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.
    – benaviru
    Jan 20, 2022 at 22:00

1 Answer 1

0

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

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.