2

I am trying to convert an QGIS3 plugin to the Processing Toolbox in QGIS3 but I am getting stuck on how to save my processed numpy array back to a raster file in the end.

How can I save the variable "test" back to a raster in the code example below?

class ProcessingWallHeightAscpetAlgorithm(QgsProcessingAlgorithm):

    OUTPUT_HEIGHT = 'OUTPUT_HEIGHT'
    INPUT_LIMIT = 'INPUT_LIMIT'
    INPUT = 'INPUT'

    def initAlgorithm(self, config):
        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.INPUT,
                self.tr('Input building and ground DSM'), 
                None, False))

        self.addParameter(QgsProcessingParameterNumber(
            self.INPUT_LIMIT, 
            self.tr("Lower limit for wall height (m)"), 
            QgsProcessingParameterNumber.Double,
            QVariant(3.0)))

        self.addParameter(QgsProcessingParameterRasterDestination(
            self.OUTPUT_HEIGHT,
            self.tr("Output Wall Height Raster"),
            None, False))


    def processAlgorithm(self, parameters, context, feedback):
        output_path_raster_height = self.parameterAsOutputLayer(parameters, self.OUTPUT_HEIGHT, context)
        raster = self.parameterAsRasterLayer(parameters, self.INPUT, context)
        provider = raster.dataProvider()
        filepath_dsm = str(provider.dataSourceUri())
        gdal_dsm = gdal.Open(filepath_dsm)
        raster = gdal_dsm.ReadAsArray().astype(np.float)
        limit_val = self.parameterAsDouble(parameters, self.INPUT_LIMIT, context)

        test = raster * limit_val

        results = {}
        results[self.OUTPUT_HEIGHT] = test

        return results

    def name(self):
        return 'Urban Geometry: Wall Height and Aspect'

    def displayName(self):
        return self.tr(self.name())

    def group(self):
       return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return ProcessingWallHeightAscpetAlgorithm()

2 Answers 2

1

Ok, so now I understand the logic of how to write output in the Processing Tool box. It is a combination of the "old" way and then just specify the saved filepath to the return variable in the processingAlgorithm function as examplified below:

class ProcessingWallHeightAscpetAlgorithm(QgsProcessingAlgorithm):

    OUTPUT_HEIGHT = 'OUTPUT_HEIGHT'
    INPUT_LIMIT = 'INPUT_LIMIT'
    INPUT = 'INPUT'

    def initAlgorithm(self, config):
        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.INPUT,
                self.tr('Input building and ground DSM'), 
                None, False))

        self.addParameter(QgsProcessingParameterNumber(
            self.INPUT_LIMIT, 
            self.tr("Lower limit for wall height (m)"), 
            QgsProcessingParameterNumber.Double,
            QVariant(3.0)))

        self.addParameter(QgsProcessingParameterRasterDestination(
            self.OUTPUT_HEIGHT,
            self.tr("Output Wall Height Raster"),
            None, False))


    def processAlgorithm(self, parameters, context, feedback):
        output_path_raster_height = self.parameterAsOutputLayer(parameters, self.OUTPUT_HEIGHT, context)
        raster = self.parameterAsRasterLayer(parameters, self.INPUT, context)
        provider = raster.dataProvider()
        filepath_dsm = str(provider.dataSourceUri())
        gdal_dsm = gdal.Open(filepath_dsm)
        raster = gdal_dsm.ReadAsArray().astype(np.float)
        limit_val = self.parameterAsDouble(parameters, self.INPUT_LIMIT, context)

        test = raster * limit_val

        saverasternd(gdal_dsm, output_path_raster_height, test)

        return {self.OUTPUT_HEIGHT: output_path_raster_height}

    def name(self):
        return 'Urban Geometry: Wall Height and Aspect'

    def displayName(self):
        return self.tr(self.name())

    def group(self):
       return QCoreApplication.translate('Processing', string)

    def createInstance(self):
        return ProcessingWallHeightAscpetAlgorithm()

    def saverasternd(gdal_data, filename, raster):
        rows = gdal_data.RasterYSize
        cols = gdal_data.RasterXSize

        outDs = gdal.GetDriverByName("GTiff").Create(filename, cols, rows, int(1), GDT_Float32)
        outBand = outDs.GetRasterBand(1)

        # write the data
        outBand.WriteArray(raster, 0, 0)
        # flush data to disk, set the NoData value and calculate stats
        outBand.FlushCache()
        # outBand.SetNoDataValue(-9999)

        # georeference the image and set the projection
        outDs.SetGeoTransform(gdal_data.GetGeoTransform())
        outDs.SetProjection(gdal_data.GetProjection())

Very simple when you understood the logic.

0

You are already using gdal to open the raster so you might as well use it to write the array to a raster file.

Your code could look something like:

# create empty raster
driver = gdal.GetDriverByName('GTiff')
out_tif = driver.Create(fn, test.shape[1], test.shape[0], 1, gdal.GDT_Float64)  # I am assuming you only have one band, hence the 1

# set spatial reference and geotransform (taken from the original raster)
out_tif.SetProjection(gdal_dsm.GetProjection())
out_tif.SetGeoTransform(gdal_dsm.GetGeoTransform())

# write array
band = out_tif.GetRasterBand(1)
nd_val = gdal_dsm.GetRasterBand(1).GetNoDataValue()
if nd_val:
    band.SetNoDataValue(nd_val)
band.WriteArray(test)

# save
band.FlushCache()
del out_tif, band
1
  • Yes, that would be the "old" way of saving it but I want to make use of the capabilities in the processing class in QGIS by using e.g. QgsProcessingParameterRasterDestination().
    – biglimp
    Apr 3, 2020 at 14:35

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.