3

I want to use the raster calculator to add multiple rasters in PyQGIS [raster 1 + raster 2 + raster ...}.

The rasters are from a for-loop. I do not manage to get the temporary files from inside the for-loop.

My code so far:

from qgis.core import QgsProcessing
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterRasterLayer
from qgis.core import QgsProcessingParameterNumber
from qgis.core import QgsProcessingParameterPoint
from qgis.core import QgsProcessingParameterRasterDestination
from qgis.core import QgsProcessingParameterVectorLayer
from qgis.core import QgsRasterLayer
from qgis.analysis import QgsRasterCalculatorEntry
from qgis.core import QgsRectangle
import processing


class Asd(QgsProcessingAlgorithm):

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterRasterLayer('dem', 'DEM', defaultValue=None))
        self.addParameter(QgsProcessingParameterVectorLayer('windrder', 'Windräder Punktdaten', types=[QgsProcessing.TypeVectorPoint], defaultValue=None))
        self.addParameter(QgsProcessingParameterNumber('hhe_der_anlage_m_ber_dem_boden', 'Höhe der Anlage [m über dem Boden]', type=QgsProcessingParameterNumber.Double, minValue=0, defaultValue=150))
        self.addParameter(QgsProcessingParameterNumber('hhe_der_person_m_ber_dem_boden', 'Höhe der Person [m über dem Boden]', type=QgsProcessingParameterNumber.Double, minValue=0, defaultValue=1.5))
        self.addParameter(QgsProcessingParameterNumber('maximaler_suchradius_m', 'Maximaler Suchradius [m]', type=QgsProcessingParameterNumber.Double, minValue=0, defaultValue=10))
        self.addParameter(QgsProcessingParameterRasterDestination('Output', 'output', createByDefault=True, defaultValue=None))
        self.addParameter(QgsProcessingParameterRasterDestination('raster_calculator_out', 'raster_calculator_out', createByDefault=True, defaultValue=None))

    def processAlgorithm(self, parameters, context, model_feedback):
        feedback = QgsProcessingMultiStepFeedback(7, model_feedback)
        results = {}
        outputs = {}

        # Get the input point layer
        point_layer = self.parameterAsVectorLayer(parameters, 'windrder', context)
        if not point_layer:
            raise QgsProcessingException('Point layer not found')

        # Iterate through each feature in the point layer
        for feature in point_layer.getFeatures():
            point = feature.geometry().asPoint()

            # Perform viewshed analysis for each point
            alg_params = {
                '-b': True,
                '-c': False,
                '-e': False,
                '-r': False,
                'GRASS_RASTER_FORMAT_META': '',
                'GRASS_RASTER_FORMAT_OPT': '',
                'GRASS_REGION_CELLSIZE_PARAMETER': 0,
                'GRASS_REGION_PARAMETER': None,
                'coordinates': f"{point.x()}, {point.y()}",
                'input': parameters['dem'],
                'max_distance': parameters['maximaler_suchradius_m'],
                'memory': 1000,
                'observer_elevation': parameters['hhe_der_anlage_m_ber_dem_boden'],
                'refraction_coeff': 0.24286,
                'target_elevation': parameters['hhe_der_person_m_ber_dem_boden'],
                'output': QgsProcessing.TEMPORARY_OUTPUT
            }
            result = processing.run('grass7:r.viewshed', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

            # Accumulate the resulting raster
            if 'output' in result:
                outputs[feature.id()] = result['output']


        # -- my failed approach --

        expr = ''
        rasters = []
        for i, path in enumerate(outputs.values()):
            inputrasterfile = QgsRasterLayer(path, f'raster{i}')
            rasters.append(inputrasterfile)

            expr = expr + '"{}@1"+'.format(f'raster{i}')
        expr = expr.rstrip('+')

        extents = [layer.extent() for layer in rasters]
        # Get the combined extent of all layers
        combined_extent = QgsRectangle()
        for extent in extents:
            combined_extent.combineExtentWith(extent)
        print(expr)
        
        alg_params = {
            'CELLSIZE': 0,
            'CRS': None,
            'EXPRESSION': expr,
            'EXTENT': combined_extent,
            'LAYERS': rasters,
            'OUTPUT': parameters['raster_calculator_out']
        }
        print(alg_params)

        outputs['raster_calculator'] = processing.run('qgis:rastercalculator', alg_params)
        results['raster_calculator_out'] = outputs['raster_calculator']['OUTPUT']


        return results

    def name(self):
        return 'asd'

    def displayName(self):
        return 'asd'

    def group(self):
        return ''

    def groupId(self):
        return ''

    def createInstance(self):
        return Asd()

output:

"raster0@1"+"raster1@1"+"raster2@1"+"raster3@1"+"raster4@1"+"raster5@1"+"raster6@1". 
{'CELLSIZE': 0, 'CRS': None, 'EXPRESSION': '', 'EXTENT': <QgsRectangle: 2725000 1221000, 2726000 1222000>, 'LAYERS': [<QgsRasterLayer: 'raster0' (gdal)>, <QgsRasterLayer: 'raster1' (gdal)>, <QgsRasterLayer: 'raster2' (gdal)>, <QgsRasterLayer: 'raster3' (gdal)>, <QgsRasterLayer: 'raster4' (gdal)>, <QgsRasterLayer: 'raster5' (gdal)>, <QgsRasterLayer: 'raster6' (gdal)>], 'OUTPUT': <QgsProcessingOutputLayerDefinition {'sink':TEMPORARY_OUTPUT, 'createOptions': {'fileEncoding': 'windows-1252'}}>} 

If I run the processing script, the following Error occurs:

An error occurred while performing the calculation
Execution failed after 5.36 Sekunden

Sample data is here: https://github.com/widmerc/sample_data (Sampe_Data_viewshed.gpkg and Sample_script_viewshed.txt)

So in summary: I want to iterate through a point-layer and add all resulting viewshed models to one result layer without getting an error in the script. So we can see how many points (from the point layer input) are visible from each raster cell.

12
  • Do you get any errors?
    – Taras
    Commented Feb 11 at 20:47
  • @SonofaBeach The problem is stated above. I edited the question. There is an error in the code. It would be great, if you would try the script yourself and adjust the bottom approach of my own. If you found the issue, I can (if you want) formulate a general understanding, how to use a for-loop in a processing script. Commented Feb 12 at 8:06
  • @Taras I have edited the question. The error 'An error occurred while performing the calculation' is very vague. Commented Feb 12 at 8:07
  • 1
    If I run your script, I don't get any errors. However, the final result raster layer is None. Why are you using qgis:rastercalculator? Have you tried native:rastercalc or native:virtualrastercalc?
    – Andre Geo
    Commented Feb 13 at 13:36
  • 1
    Please edit your original post to include all additional relevant information. Commented Feb 15 at 21:32

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.