Merge pull request #50586 from nicogodet/gdalcalcextent

Add user interface for `--extent` gdalcalc option
qgis · Nov 17, 2022 · 60e2c98 · 60e2c98
2 parents b60ff3c + 94c0998
commit 60e2c98
Show file tree

Hide file tree

Showing 2 changed files with 79 additions and 5 deletions.
diff --git a/python/plugins/processing/algs/gdal/gdalcalc.py b/python/plugins/processing/algs/gdal/gdalcalc.py
@@ -50,6 +50,10 @@ class gdalcalc(GdalAlgorithm):
     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'
@@ -141,6 +145,15 @@ def initAlgorithm(self, config=None):
                 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'),
@@ -192,11 +205,6 @@ def groupId(self):
     def commandName(self):
         return 'gdal_calc'
 
-    def processAlgorithm(self, parameters, context, feedback):
-        if GdalUtils.version() < 3030000 and self.EXTENT in parameters.keys():
-            raise QgsProcessingException(self.tr('The output extent option is only available on GDAL 3.3 or later'))
-        return GdalAlgorithm.processAlgorithm(self, parameters, context, feedback)
-
     def getConsoleCommands(self, parameters, context, feedback, executing=True):
 
         out = self.parameterAsOutputLayer(parameters, self.OUTPUT, context)
@@ -223,6 +231,38 @@ def getConsoleCommands(self, parameters, context, feedback, executing=True):
         if layer is None:
             raise QgsProcessingException(self.invalidRasterError(parameters, self.INPUT_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.crs())
         if not bbox.isNull():
             arguments.append('--projwin')

diff --git a/python/plugins/processing/tests/GdalAlgorithmsRasterTest.py b/python/plugins/processing/tests/GdalAlgorithmsRasterTest.py
@@ -685,6 +685,40 @@ def testGdalCalc(self):
                     ['gdal_calc.py',
                      '--overwrite --calc "{}" --format JPEG --type Float32 --projwin 1.0 4.0 3.0 2.0 -A {} --A_band 1 --outfile {}'.format(formula, source, output)])
 
+                # Inputs A and B share same pixel size and CRS
+                self.assertEqual(
+                    alg.getConsoleCommands({'INPUT_A': source,
+                                            'BAND_A': 1,
+                                            'INPUT_B': source,
+                                            'BAND_B': 1,
+                                            'FORMULA': formula,
+                                            'EXTENT_OPT': 3,
+                                            'OUTPUT': output}, context, feedback),
+                    ['gdal_calc.py',
+                     '--overwrite --calc "{}" --format JPEG --type Float32 --extent=intersect -A {} --A_band 1 -B {} --B_band 1 --outfile {}'.format(formula, source, source, output)])
+
+                # Test mutually exclusive --extent and --projwin. Should raise an exception
+                self.assertRaises(
+                    QgsProcessingException,
+                    lambda: alg.getConsoleCommands({'INPUT_A': source,
+                                                    'BAND_A': 1,
+                                                    'FORMULA': formula,
+                                                    'PROJWIN': extent,
+                                                    'EXTENT_OPT': 3,
+                                                    'OUTPUT': output}, context, feedback))
+
+                # Inputs A and B do not share same pixel size and CRS. Should raise an exception
+                source2 = os.path.join(testDataPath, 'raster.tif')
+                self.assertRaises(
+                    QgsProcessingException,
+                    lambda: alg.getConsoleCommands({'INPUT_A': source,
+                                                    'BAND_A': 1,
+                                                    'INPUT_B': source2,
+                                                    'BAND_B': 1,
+                                                    'FORMULA': formula,
+                                                    'EXTENT_OPT': 3,
+                                                    'OUTPUT': output}, context, feedback))
+
             # check that formula is not escaped and formula is returned as it is
             formula = 'A * 2'  # <--- add spaces in the formula
             self.assertEqual(