0

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.

These are the three main steps for which I have extracted the Python code using the graphical modeler

#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
3
  • 1
    Welcome to Geographic Information Systems! Welcome to GIS SE! We're a little different from other sites; this isn't a discussion forum but a Q&A site. Your questions should as much as possible describe not just what you want to do, but precisely what you have tried and where you are stuck trying that. Please check out our short tour for more about how the site works.
    – Ian Turton
    Commented Jun 30, 2020 at 9:50
  • 2
    Generally, we're looking for some code to show us how far you got and where you are stuck
    – Ian Turton
    Commented Jun 30, 2020 at 9:51
  • I've added the code I have for now but it's very long for not much, I'd just like to know if it is possible to create the type of loop I'm looking for Commented Jun 30, 2020 at 13:25

1 Answer 1

0

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.

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.