TransWikia.com

Save numpy array as raster in a Processing toolbox QGIS3

Geographic Information Systems Asked by biglimp on April 5, 2021

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

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

Answered by Marcelo Villa-Piñeros on April 5, 2021

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.

Answered by biglimp on April 5, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP