[processing] improve Create constant raster algorithm

Author: alexbruy

python/plugins/processing/algs/qgis/CreateConstantRaster.py

@@ -25,20 +25,27 @@
 
 __revision__ = '$Format:%H$'
 
-from osgeo import gdal
+import math
+import struct
 
-from qgis.core import (QgsProcessingParameterRasterLayer,
+from qgis.core import (Qgis,
+                       QgsRasterBlock,
+                       QgsRasterFileWriter,
+                       QgsProcessingParameterExtent,
                        QgsProcessingParameterNumber,
+                       QgsProcessingParameterCrs,
                        QgsProcessingParameterRasterDestination)
 from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
-from processing.tools.raster import RasterWriter
+from processing.tools import raster
 
 
 class CreateConstantRaster(QgisAlgorithm):
 
-    INPUT = 'INPUT'
-    OUTPUT = 'OUTPUT'
+    EXTENT = 'EXTENT'
+    TARGET_CRS = 'TARGET_CRS'
+    PIXEL_SIZE = 'PIXEL_SIZE'
     NUMBER = 'NUMBER'
+    OUTPUT = 'OUTPUT'
 
     def group(self):
         return self.tr('Raster tools')
@@ -47,10 +54,18 @@ def __init__(self):
         super().__init__()
 
     def initAlgorithm(self, config=None):
-        self.addParameter(QgsProcessingParameterRasterLayer(self.INPUT,
-                                                            self.tr('Reference layer')))
+        self.addParameter(QgsProcessingParameterExtent(self.EXTENT,
+                                                       self.tr('Desired extent')))
+        self.addParameter(QgsProcessingParameterCrs(self.TARGET_CRS,
+                                                    self.tr('Target CRS'),
+                                                    'ProjectCrs'))
+        self.addParameter(QgsProcessingParameterNumber(self.PIXEL_SIZE,
+                                                       self.tr('Pixel size'),
+                                                       QgsProcessingParameterNumber.Double,
+                                                       defaultValue=0.1))
         self.addParameter(QgsProcessingParameterNumber(self.NUMBER,
-                                                       self.tr('Constant value'), QgsProcessingParameterNumber.Double,
+                                                       self.tr('Constant value'),
+                                                       QgsProcessingParameterNumber.Double,
                                                        defaultValue=1))
         self.addParameter(QgsProcessingParameterRasterDestination(self.OUTPUT, self.tr('Constant')))
 
@@ -61,28 +76,35 @@ def displayName(self):
         return self.tr('Create constant raster layer')
 
     def processAlgorithm(self, parameters, context, feedback):
-        layer = self.parameterAsRasterLayer(parameters, self.INPUT, context)
+        extent = self.parameterAsExtent(parameters, self.EXTENT, context)
+        crs = self.parameterAsCrs(parameters, self.TARGET_CRS, context)
         value = self.parameterAsDouble(parameters, self.NUMBER, context)
+        pixelSize = self.parameterAsDouble(parameters, self.PIXEL_SIZE, context)
 
         outputFile = self.parameterAsOutputLayer(parameters, self.OUTPUT, context)
+        outputFormat = raster.formatShortNameFromFileName(outputFile)
+
+        rows = max([math.ceil(extent.height() / pixelSize) + 1, 1.0])
+        cols = max([math.ceil(extent.width() / pixelSize) + 1, 1.0])
+
+        writer = QgsRasterFileWriter(outputFile)
+        writer.setOutputProviderKey('gdal')
+        writer.setOutputFormat(outputFormat)
+        provider = writer.createOneBandRaster(Qgis.Float32, cols, rows, extent, crs)
+        provider.setNoDataValue(1, -9999)
+
+        data = [value] * cols
+        block = QgsRasterBlock(Qgis.Float32, cols, 1)
+        block.setData(struct.pack('{}f'.format(len(data)), *data))
+
+        total = 100.0 / rows if rows else 0
+        for i in range(rows):
+            if feedback.isCanceled():
+                break
+
+            provider.writeBlock(block, 1, 0, i)
+            feedback.setProgress(int(i * rows))
 
-        raster = gdal.Open(layer.source(), gdal.GA_ReadOnly)
-        geoTransform = raster.GetGeoTransform()
-
-        cellsize = (layer.extent().xMaximum() - layer.extent().xMinimum()) \
-            / layer.width()
-
-        w = RasterWriter(outputFile,
-                         layer.extent().xMinimum(),
-                         layer.extent().yMinimum(),
-                         layer.extent().xMaximum(),
-                         layer.extent().yMaximum(),
-                         cellsize,
-                         1,
-                         layer.crs(),
-                         geoTransform
-                         )
-        w.matrix.fill(value)
-        w.close()
+        provider.setEditable(False)
 
         return {self.OUTPUT: outputFile}

python/plugins/processing/tests/testdata/expected/constant_raster.tif

@@ -1223,16 +1223,15 @@ tests:
   - algorithm: qgis:createconstantrasterlayer
     name: Create constant raster
     params:
-      INPUT:
-        name: raster.tif
-        type: raster
+      EXTENT: 270736.0673250682,270899.8544675339,4458899.000550019,4459029.574521748
+      TARGET_CRS: EPSG:23030
+      PIXEL_SIZE: 10.0
       NUMBER: 3.0
     results:
       OUTPUT:
-        hash: 4cb3e82e8512cdbb75d9c39a10adc818dd6842c5dc6361fbc43dd9aa
+        hash: e453e9e36ce314d5197963ac27872a0cc3dfe43764ed586a334c66f0
         type: rasterhash
 
-
   # Case 1: Keep all fields
   - algorithm: qgis:lineintersections
     name: Line Intersection Keep All Fields from Both