""" *************************************************************************** gdalcalc.py --------------------- Date : Janaury 2015 Copyright : (C) 2015 by Giovanni Manghi Email : giovanni dot manghi at naturalgis dot pt *************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * *************************************************************************** """ __author__ = "Giovanni Manghi" __date__ = "January 2015" __copyright__ = "(C) 2015, Giovanni Manghi" from qgis.core import ( QgsProcessingException, QgsProcessingParameterDefinition, QgsProcessingParameterRasterLayer, QgsProcessingParameterBand, QgsProcessingParameterNumber, QgsProcessingParameterEnum, QgsProcessingParameterExtent, QgsProcessingParameterString, QgsProcessingParameterRasterDestination, ) from processing.algs.gdal.GdalAlgorithm import GdalAlgorithm from processing.algs.gdal.GdalUtils import GdalUtils from processing.tools.system import isWindows class gdalcalc(GdalAlgorithm): INPUT_A = "INPUT_A" INPUT_B = "INPUT_B" INPUT_C = "INPUT_C" INPUT_D = "INPUT_D" INPUT_E = "INPUT_E" INPUT_F = "INPUT_F" BAND_A = "BAND_A" BAND_B = "BAND_B" BAND_C = "BAND_C" BAND_D = "BAND_D" BAND_E = "BAND_E" BAND_F = "BAND_F" FORMULA = "FORMULA" # TODO QGIS 4.0 : Rename EXTENT_OPT to EXTENT EXTENT_OPT = "EXTENT_OPT" EXTENT_OPTIONS = ["ignore", "fail", "union", "intersect"] # TODO QGIS 4.0 : Rename EXTENT to PROJWIN or CUSTOM_EXTENT EXTENT = "PROJWIN" OUTPUT = "OUTPUT" NO_DATA = "NO_DATA" OPTIONS = "OPTIONS" CREATION_OPTIONS = "CREATION_OPTIONS" EXTRA = "EXTRA" RTYPE = "RTYPE" TYPE = [ "Byte", "Int16", "UInt16", "UInt32", "Int32", "Float32", "Float64", "CInt16", "CInt32", "CFloat32", "CFloat64", "Int8", ] def __init__(self): super().__init__() def initAlgorithm(self, config=None): self.addParameter( QgsProcessingParameterRasterLayer( self.INPUT_A, self.tr("Input layer A"), optional=False ) ) self.addParameter( QgsProcessingParameterBand( self.BAND_A, self.tr("Number of raster band for A"), parentLayerParameterName=self.INPUT_A, ) ) self.addParameter( QgsProcessingParameterRasterLayer( self.INPUT_B, self.tr("Input layer B"), optional=True ) ) self.addParameter( QgsProcessingParameterBand( self.BAND_B, self.tr("Number of raster band for B"), parentLayerParameterName=self.INPUT_B, optional=True, ) ) self.addParameter( QgsProcessingParameterRasterLayer( self.INPUT_C, self.tr("Input layer C"), optional=True ) ) self.addParameter( QgsProcessingParameterBand( self.BAND_C, self.tr("Number of raster band for C"), parentLayerParameterName=self.INPUT_C, optional=True, ) ) self.addParameter( QgsProcessingParameterRasterLayer( self.INPUT_D, self.tr("Input layer D"), optional=True ) ) self.addParameter( QgsProcessingParameterBand( self.BAND_D, self.tr("Number of raster band for D"), parentLayerParameterName=self.INPUT_D, optional=True, ) ) self.addParameter( QgsProcessingParameterRasterLayer( self.INPUT_E, self.tr("Input layer E"), optional=True ) ) self.addParameter( QgsProcessingParameterBand( self.BAND_E, self.tr("Number of raster band for E"), parentLayerParameterName=self.INPUT_E, optional=True, ) ) self.addParameter( QgsProcessingParameterRasterLayer( self.INPUT_F, self.tr("Input layer F"), optional=True ) ) self.addParameter( QgsProcessingParameterBand( self.BAND_F, self.tr("Number of raster band for F"), parentLayerParameterName=self.INPUT_F, optional=True, ) ) self.addParameter( QgsProcessingParameterString( self.FORMULA, self.tr( "Calculation in gdalnumeric syntax using +-/* or any numpy array functions (i.e. logical_and())" ), "A*2", optional=False, ) ) self.addParameter( QgsProcessingParameterNumber( self.NO_DATA, self.tr("Set output NoData value"), type=QgsProcessingParameterNumber.Type.Double, defaultValue=None, optional=True, ) ) if GdalUtils.version() >= 3030000: extent_opt_param = QgsProcessingParameterEnum( self.EXTENT_OPT, self.tr("Handling of extent differences"), options=[o.title() for o in self.EXTENT_OPTIONS], defaultValue=0, ) extent_opt_param.setHelp( self.tr( "This option determines how to handle rasters with different extents" ) ) self.addParameter(extent_opt_param) if GdalUtils.version() >= 3030000: extent_param = QgsProcessingParameterExtent( self.EXTENT, self.tr("Output extent"), optional=True ) extent_param.setHelp(self.tr("Custom extent of the output raster")) self.addParameter(extent_param) self.addParameter( QgsProcessingParameterEnum( self.RTYPE, self.tr("Output raster type"), options=self.TYPE, defaultValue=5, ) ) # backwards compatibility parameter # TODO QGIS 4: remove parameter and related logic options_param = QgsProcessingParameterString( self.OPTIONS, self.tr("Additional creation options"), defaultValue="", optional=True, ) options_param.setFlags( options_param.flags() | QgsProcessingParameterDefinition.Flag.Hidden ) options_param.setMetadata({"widget_wrapper": {"widget_type": "rasteroptions"}}) self.addParameter(options_param) creation_options_param = QgsProcessingParameterString( self.CREATION_OPTIONS, self.tr("Additional creation options"), defaultValue="", optional=True, ) creation_options_param.setFlags( creation_options_param.flags() | QgsProcessingParameterDefinition.Flag.FlagAdvanced ) creation_options_param.setMetadata( {"widget_wrapper": {"widget_type": "rasteroptions"}} ) self.addParameter(creation_options_param) extra_param = QgsProcessingParameterString( self.EXTRA, self.tr("Additional command-line parameters"), defaultValue=None, optional=True, ) extra_param.setFlags( extra_param.flags() | QgsProcessingParameterDefinition.Flag.FlagAdvanced ) self.addParameter(extra_param) self.addParameter( QgsProcessingParameterRasterDestination(self.OUTPUT, self.tr("Calculated")) ) def name(self): return "rastercalculator" def displayName(self): return self.tr("Raster calculator") def group(self): return self.tr("Raster miscellaneous") def groupId(self): return "rastermiscellaneous" def commandName(self): return "gdal_calc" def getConsoleCommands(self, parameters, context, feedback, executing=True): out = self.parameterAsOutputLayer(parameters, self.OUTPUT, context) self.setOutputValue(self.OUTPUT, out) formula = self.parameterAsString(parameters, self.FORMULA, context) if self.NO_DATA in parameters and parameters[self.NO_DATA] is not None: noData = self.parameterAsDouble(parameters, self.NO_DATA, context) else: noData = None arguments = [ "--overwrite", f'--calc "{formula}"', "--format", GdalUtils.getFormatShortNameFromFilename(out), ] rtype = self.parameterAsEnum(parameters, self.RTYPE, context) if ( self.TYPE[rtype] in ["CInt16", "CInt32", "CFloat32", "CFloat64"] and GdalUtils.version() < 3050300 ): raise QgsProcessingException( self.tr("{} data type requires GDAL version 3.5.3 or later").format( self.TYPE[rtype] ) ) if self.TYPE[rtype] == "Int8" and GdalUtils.version() < 3070000: raise QgsProcessingException( self.tr("Int8 data type requires GDAL version 3.7 or later") ) arguments.append("--type " + self.TYPE[rtype]) if noData is not None: arguments.append("--NoDataValue") arguments.append(noData) layer_a = self.parameterAsRasterLayer(parameters, self.INPUT_A, context) if layer_a is None: raise QgsProcessingException( self.invalidRasterError(parameters, self.INPUT_A) ) layer_a_details = GdalUtils.gdal_connection_details_from_layer(layer_a) def all_equal(iterator): iterator = iter(iterator) try: first = next(iterator) except StopIteration: return True return all(first == x for x in iterator) # Check GDAL version for projwin and extent options (GDAL 3.3 is required) if GdalUtils.version() < 3030000 and self.EXTENT in parameters.keys(): raise QgsProcessingException( self.tr( "The custom output extent option (--projwin) is only available on GDAL 3.3 or later" ) ) if GdalUtils.version() < 3030000 and self.EXTENT_OPT in parameters.keys(): raise QgsProcessingException( self.tr( "The output extent option (--extent) is only available on GDAL 3.3 or later" ) ) # --projwin and --extent option are mutually exclusive if ( self.EXTENT in parameters.keys() and parameters[self.EXTENT] is not None ) and ( self.EXTENT_OPT in parameters.keys() and parameters[self.EXTENT_OPT] != 0 ): raise QgsProcessingException( self.tr( "The custom output extent option (--projwin) and output extent option (--extent) are mutually exclusive" ) ) # If extent option is defined, pixel size and SRS of all input raster must be the same if self.EXTENT_OPT in parameters.keys() and parameters[self.EXTENT_OPT] != 0: pixel_size_X, pixel_size_Y, srs = [], [], [] for input_layer in [ self.INPUT_A, self.INPUT_B, self.INPUT_C, self.INPUT_D, self.INPUT_E, self.INPUT_F, ]: if input_layer in parameters and parameters[input_layer] is not None: layer = self.parameterAsRasterLayer( parameters, input_layer, context ) pixel_size_X.append(layer.rasterUnitsPerPixelX()) pixel_size_Y.append(layer.rasterUnitsPerPixelY()) srs.append(layer.crs().authid()) if not ( all_equal(pixel_size_X) and all_equal(pixel_size_Y) and all_equal(srs) ): raise QgsProcessingException( self.tr( "For all output extent options, the pixel size (resolution) and SRS (Spatial Reference System) of all the input rasters must be the same" ) ) extent = self.EXTENT_OPTIONS[ self.parameterAsEnum(parameters, self.EXTENT_OPT, context) ] if extent != "ignore": arguments.append(f"--extent={extent}") bbox = self.parameterAsExtent(parameters, self.EXTENT, context, layer_a.crs()) if not bbox.isNull(): arguments.append("--projwin") arguments.append(str(bbox.xMinimum())) arguments.append(str(bbox.yMaximum())) arguments.append(str(bbox.xMaximum())) arguments.append(str(bbox.yMinimum())) arguments.append("-A") arguments.append(layer_a_details.connection_string) if self.parameterAsString(parameters, self.BAND_A, context): arguments.append( "--A_band " + self.parameterAsString(parameters, self.BAND_A, context) ) if self.INPUT_B in parameters and parameters[self.INPUT_B] is not None: layer_b = self.parameterAsRasterLayer(parameters, self.INPUT_B, context) if layer_b is None: raise QgsProcessingException( self.invalidRasterError(parameters, self.INPUT_B) ) input_b_details = GdalUtils.gdal_connection_details_from_layer(layer_b) arguments.append("-B") arguments.append(input_b_details.connection_string) if self.parameterAsString(parameters, self.BAND_B, context): arguments.append( "--B_band " + self.parameterAsString(parameters, self.BAND_B, context) ) if self.INPUT_C in parameters and parameters[self.INPUT_C] is not None: layer_c = self.parameterAsRasterLayer(parameters, self.INPUT_C, context) if layer_c is None: raise QgsProcessingException( self.invalidRasterError(parameters, self.INPUT_C) ) input_c_details = GdalUtils.gdal_connection_details_from_layer(layer_c) arguments.append("-C") arguments.append(input_c_details.connection_string) if self.parameterAsString(parameters, self.BAND_C, context): arguments.append( "--C_band " + self.parameterAsString(parameters, self.BAND_C, context) ) if self.INPUT_D in parameters and parameters[self.INPUT_D] is not None: layer_d = self.parameterAsRasterLayer(parameters, self.INPUT_D, context) if layer_d is None: raise QgsProcessingException( self.invalidRasterError(parameters, self.INPUT_D) ) input_d_details = GdalUtils.gdal_connection_details_from_layer(layer_d) arguments.append("-D") arguments.append(input_d_details.connection_string) if self.parameterAsString(parameters, self.BAND_D, context): arguments.append( "--D_band " + self.parameterAsString(parameters, self.BAND_D, context) ) if self.INPUT_E in parameters and parameters[self.INPUT_E] is not None: layer_e = self.parameterAsRasterLayer(parameters, self.INPUT_E, context) if layer_e is None: raise QgsProcessingException( self.invalidRasterError(parameters, self.INPUT_E) ) input_e_details = GdalUtils.gdal_connection_details_from_layer(layer_e) arguments.append("-E") arguments.append(input_e_details.connection_string) if self.parameterAsString(parameters, self.BAND_E, context): arguments.append( "--E_band " + self.parameterAsString(parameters, self.BAND_E, context) ) if self.INPUT_F in parameters and parameters[self.INPUT_F] is not None: layer_f = self.parameterAsRasterLayer(parameters, self.INPUT_F, context) if layer_f is None: raise QgsProcessingException( self.invalidRasterError(parameters, self.INPUT_F) ) input_f_details = GdalUtils.gdal_connection_details_from_layer(layer_f) arguments.append("-F") arguments.append(input_f_details.connection_string) if self.parameterAsString(parameters, self.BAND_F, context): arguments.append( "--F_band " + self.parameterAsString(parameters, self.BAND_F, context) ) options = self.parameterAsString(parameters, self.CREATION_OPTIONS, context) # handle backwards compatibility parameter OPTIONS if self.OPTIONS in parameters and parameters[self.OPTIONS] not in (None, ""): options = self.parameterAsString(parameters, self.OPTIONS, context) if options: parts = options.split("|") for p in parts: arguments.append("--co " + p) if self.EXTRA in parameters and parameters[self.EXTRA] not in (None, ""): extra = self.parameterAsString(parameters, self.EXTRA, context) arguments.append(extra) arguments.append("--outfile") arguments.append(out) return [ self.commandName() + (".bat" if isWindows() else ".py"), GdalUtils.escapeAndJoin(arguments), ]