Creating a Python loop to add up automatically generated rasters on QGIS 3

Question

I'm creating a QGIS 3 program that extracts flow paths from a DTM, then slightly modifies the DTM to extract new flow paths that I want to add to the previous flow path raster.

I have used the graphical modeler to carry out the different steps and then extracted the associated Python code but I couldn't find the way to unite them into one Python code, that could open the initial DTM and create a loop that could go through these steps from 50 to 100 times.

#opening of the DTM

import processing

petitmnt = 'C:\\Users\\peter\\Documents\\OTEIS QGIS\\PetitMNT.asc'
rlayer = iface.addRasterLayer(petitmnt,'PetitMNT')  

#extraction of the flow paths (code from the graphical modeler)

from qgis.core import QgsProcessing
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterRasterLayer
from qgis.core import QgsProcessingParameterRasterDestination


class FlowPaths(QgsProcessingAlgorithm):

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterRasterLayer('PetitMNT', 'MNT', 
defaultValue=None))
      
  self.addParameter(QgsProcessingParameterRasterDestination('Prethalwegs', 'PreThalwegs', createByDefault=True, defaultValue=None))
       
 self.addParameter(QgsProcessingParameterRasterDestination('Mnt_traite', 'MNT_Traite', createByDefault=True, defaultValue=None))
    
self.addParameter(QgsProcessingParameterRasterDestination('Thalwegs_traites', 'Thalwegs_Traites', createByDefault=True, defaultValue=None))

    def processAlgorithm(self, parameters, context, model_feedback):
        # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
        # overall progress through the model
        feedback = QgsProcessingMultiStepFeedback(3, model_feedback)
        results = {}
        outputs = {}

        # Fill sinks
        alg_params = {
            'DEM': parameters['MNT'],
            'MINSLOPE': 0.01,
            'RESULT': parameters['Mnt_traite']
    }
        outputs['FillSinks'] = processing.run('saga:fillsinks', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['Mnt_traite'] = outputs['FillSinks']['RESULT']

        feedback.setCurrentStep(1)
        if feedback.isCanceled():
            return {}

        # Catchment area (recursive) finds the flow paths and creates a new raster
    alg_params = {
        'ACCU_MATERIAL': None,
        'ACCU_TARGET': None,
        'CONVERGENCE': 1.1,
        'ELEVATION': outputs['FillSinks']['RESULT'],
        'FLOW_UNIT': 1,
        'METHOD': 0,
        'NO_NEGATIVES': True,
        'SINKROUTE': None,
        'STEP': 1,
        'TARGETS': None,
        'VAL_INPUT': None,
        'WEIGHTS': None,
        'ACCU_LEFT': QgsProcessing.TEMPORARY_OUTPUT,
        'ACCU_RIGHT': QgsProcessing.TEMPORARY_OUTPUT,
        'ACCU_TOTAL': QgsProcessing.TEMPORARY_OUTPUT,
        'FLOW': parameters['Prethalwegs'],
        'FLOW_LENGTH': QgsProcessing.TEMPORARY_OUTPUT,
        'VAL_MEAN': QgsProcessing.TEMPORARY_OUTPUT,
        'WEIGHT_LOSS': QgsProcessing.TEMPORARY_OUTPUT
    }
    outputs['CatchmentAreaRecursive'] = processing.run('saga:catchmentarearecursive', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
    results['Prethalwegs'] = outputs['CatchmentAreaRecursive']['FLOW']

    feedback.setCurrentStep(2)
    if feedback.isCanceled():
        return {}

    # Raster calculator that eliminates values under 10
    alg_params = {
        'FORMULA': 'ifelse(gt(a, 9), a, 0)',
        'GRIDS': outputs['CatchmentAreaRecursive']['FLOW'],
        'RESAMPLING': 3,
        'TYPE': 7,
        'USE_NODATA': False,
        'XGRIDS': [],
        'RESULT': parameters['Thalwegs_traites']
    }
    outputs['RasterCalculator'] = processing.run('saga:rastercalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
    results['Thalwegs_traites'] = outputs['RasterCalculator']['RESULT']
    return results

#creation of an automatically generated raster containing zeros or ones

import numpy as np

proportion_1 = 0.1 #Adjust
outfile = r'C:\Users\peter\Documents\OTEIS QGIS\Création du 
bruitage\RasterAlea.asc' #Adjust path

rl = QgsProject.instance().mapLayersByName('MNTSudCassagnoles')[0] #Adjust to match your raster layer name

e = rl.extent()

h = rl.height()
w = rl.width()
xres = rl.rasterUnitsPerPixelX()

arr = np.random.choice(2, (w,h), p=[1-proportion_1,proportion_1])

with open(outfile,'w') as file:
    header= ['ncols {}'.format(h),
             'nrows {}'.format(w), 
             'xllcorner {}'.format(e.xMinimum()), 
             'yllcorner {}'.format(e.yMinimum()), 
             'cellsize {}'.format(xres), 
             'nodata_value -9999']
    for row in header:
        file.write(row+'\n')
    for row in arr.tolist():
        file.write(' '.join([str(val) for val in row]))

#modification of the DTM using the previously generated raster (adding the random raster to the DTM

class DTM_modif(QgsProcessingAlgorithm):

def initAlgorithm(self, config=None):
    self.addParameter(QgsProcessingParameterRasterLayer('MNT', 'MNT', defaultValue=None))
    self.addParameter(QgsProcessingParameterRasterLayer('array', 'RasterAlea', defaultValue=None))
    
self.addParameter(QgsProcessingParameterRasterDestination('Rasteralea_20cm', 'RasterAlea_20cm', createByDefault=True, defaultValue=None))
    self.addParameter(QgsProcessingParameterRasterDestination('Mnt_modifie', 'MNT_Modifie', createByDefault=True, defaultValue=None))

def processAlgorithm(self, parameters, context, model_feedback):
    # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
    # overall progress through the model
    feedback = QgsProcessingMultiStepFeedback(2, model_feedback)
    results = {}
    outputs = {}

    # Raster calculator that allows to modify the random raster
    alg_params = {
        'CELLSIZE': 0,
        'CRS': 'ProjectCrs',
        'EXPRESSION': '(\"RasterAlea@1\" > 0) * 0.2',
        'EXTENT': None,
        'LAYERS': parameters['array'],
        'OUTPUT': parameters['Rasteralea_20cm']
    }
    outputs['RasterCalculator'] = processing.run('qgis:rastercalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
    results['Rasteralea_20cm'] = outputs['RasterCalculator']['OUTPUT']

    feedback.setCurrentStep(1)
    if feedback.isCanceled():
        return {}

    # Raster calculator that adds up the DTM and the random raster
    alg_params = {
        'FORMULA': 'a + b ',
        'GRIDS': parameters['MNT'],
        'RESAMPLING': 3,
        'TYPE': 7,
        'USE_NODATA': False,
        'XGRIDS': outputs['RasterCalculator']['OUTPUT'],
        'RESULT': parameters['Mnt_modifie']
    }
    outputs['RasterCalculator'] = processing.run('saga:rastercalculator', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
    results['Mnt_modifie'] = outputs['RasterCalculator']['RESULT']
    return results

Answer 1

It's useless to try to modify this code since there is an easier solution to this problem. You can run the plugin algorithms from the Toolbox, and simply go into the 'History' file, then simply extract the code from there. Finally, just change the paths to your input and output rasters. The loop can be easily added with basic python code.