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()
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
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP