[processing] port raster layer statistics

Author: alexbruy

python/plugins/processing/algs/help/qgis.yaml

@@ -467,9 +467,7 @@ qgis:rasterlayerhistogram: >
   The raster layer must have a single band.
 
 qgis:rasterlayerstatistics: >
-  This algorithm computes basic statistics from the values in a raster layer.
-
-  The raster layer must have a single band.
+  This algorithm computes basic statistics from the values in a given band of the raster layer.
 
 qgis:refactorfields: >
   This algorithm allows editing the structure of the attributes table of a vector layer. Fields can be modified in their type and name, using a fields mapping.

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

@@ -99,6 +99,7 @@
 from .RandomPointsExtent import RandomPointsExtent
 from .RandomPointsLayer import RandomPointsLayer
 from .RandomPointsPolygons import RandomPointsPolygons
+from .RasterLayerStatistics import RasterLayerStatistics
 from .RegularPoints import RegularPoints
 from .ReverseLineDirection import ReverseLineDirection
 from .Ruggedness import Ruggedness
@@ -144,7 +145,6 @@
 # from .HubDistanceLines import HubDistanceLines
 # from .HubLines import HubLines
 # from .GeometryConvert import GeometryConvert
-# from .RasterLayerStatistics import RasterLayerStatistics
 # from .StatisticsByCategories import StatisticsByCategories
 # from .FieldsCalculator import FieldsCalculator
 # from .FieldPyculator import FieldsPyculator
@@ -270,6 +270,7 @@ def getAlgs(self):
                 RandomPointsExtent(),
                 RandomPointsLayer(),
                 RandomPointsPolygons(),
+                RasterLayerStatistics(),
                 RegularPoints(),
                 ReverseLineDirection(),
                 Ruggedness(),

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

@@ -16,7 +16,6 @@
 *                                                                         *
 ***************************************************************************
 """
-from builtins import str
 
 __author__ = 'Victor Olaya'
 __date__ = 'January 2013'
@@ -26,30 +25,31 @@
 
 __revision__ = '$Format:%H$'
 
-import math
 import codecs
 
-from qgis.core import (QgsApplication,
-                       QgsProcessingUtils)
+from qgis.core import (QgsRectangle,
+                       QgsRasterBandStats,
+                       QgsProcessingParameterRasterLayer,
+                       QgsProcessingParameterNumber,
+                       QgsProcessingParameterFileDestination,
+                       QgsProcessingOutputHtml,
+                       QgsProcessingOutputNumber)
 from processing.algs.qgis.QgisAlgorithm import QgisAlgorithm
-from processing.core.parameters import ParameterRaster
-from processing.core.outputs import OutputNumber
-from processing.core.outputs import OutputHTML
-from processing.tools import raster
 
 
 class RasterLayerStatistics(QgisAlgorithm):
 
     INPUT = 'INPUT'
+    BAND = 'BAND'
+    OUTPUT_HTML_FILE = 'OUTPUT_HTML_FILE'
 
     MIN = 'MIN'
     MAX = 'MAX'
+    RANGE = 'RANGE'
     SUM = 'SUM'
     MEAN = 'MEAN'
-    COUNT = 'COUNT'
-    NO_DATA_COUNT = 'NO_DATA_COUNT'
     STD_DEV = 'STD_DEV'
-    OUTPUT_HTML_FILE = 'OUTPUT_HTML_FILE'
+    SUM_OF_SQUARES = 'SUM_OF_SQUARES'
 
     def group(self):
         return self.tr('Raster tools')
@@ -58,16 +58,22 @@ def __init__(self):
         super().__init__()
 
     def initAlgorithm(self, config=None):
-        self.addParameter(ParameterRaster(self.INPUT, self.tr('Input layer')))
-
-        self.addOutput(OutputHTML(self.OUTPUT_HTML_FILE, self.tr('Statistics')))
-        self.addOutput(OutputNumber(self.MIN, self.tr('Minimum value')))
-        self.addOutput(OutputNumber(self.MAX, self.tr('Maximum value')))
-        self.addOutput(OutputNumber(self.SUM, self.tr('Sum')))
-        self.addOutput(OutputNumber(self.MEAN, self.tr('Mean value')))
-        self.addOutput(OutputNumber(self.COUNT, self.tr('valid cells count')))
-        self.addOutput(OutputNumber(self.COUNT, self.tr('No-data cells count')))
-        self.addOutput(OutputNumber(self.STD_DEV, self.tr('Standard deviation')))
+        self.addParameter(QgsProcessingParameterRasterLayer(self.INPUT,
+                                                            self.tr('Input layer')))
+        self.addParameter(QgsProcessingParameterNumber(self.BAND,
+                                                       self.tr('Band number'),
+                                                       QgsProcessingParameterNumber.Integer,
+                                                       1, False, 1, 999))
+        self.addParameter(QgsProcessingParameterFileDestination(self.OUTPUT_HTML_FILE, self.tr('Statistics'), self.tr('HTML files (*.html)'), None, True))
+        self.addOutput(QgsProcessingOutputHtml(self.OUTPUT_HTML_FILE, self.tr('Statistics')))
+
+        self.addOutput(QgsProcessingOutputNumber(self.MIN, self.tr('Minimum value')))
+        self.addOutput(QgsProcessingOutputNumber(self.MAX, self.tr('Maximum value')))
+        self.addOutput(QgsProcessingOutputNumber(self.RANGE, self.tr('Range')))
+        self.addOutput(QgsProcessingOutputNumber(self.SUM, self.tr('Sum')))
+        self.addOutput(QgsProcessingOutputNumber(self.MEAN, self.tr('Mean value')))
+        self.addOutput(QgsProcessingOutputNumber(self.STD_DEV, self.tr('Standard deviation')))
+        self.addOutput(QgsProcessingOutputNumber(self.SUM_OF_SQUARES, self.tr('Sum of the squares')))
 
     def name(self):
         return 'rasterlayerstatistics'
@@ -76,62 +82,41 @@ def displayName(self):
         return self.tr('Raster layer statistics')
 
     def processAlgorithm(self, parameters, context, feedback):
-        outputFile = self.getOutputValue(self.OUTPUT_HTML_FILE)
-        uri = self.getParameterValue(self.INPUT)
-        layer = QgsProcessingUtils.mapLayerFromString(uri, context)
-        values = raster.scanraster(layer, feedback)
-
-        n = 0
-        nodata = 0
-        mean = 0
-        M2 = 0
-        sum = 0
-        minvalue = None
-        maxvalue = None
-
-        for v in values:
-            if v is not None:
-                sum += v
-                n = n + 1
-                delta = v - mean
-                mean = mean + delta / n
-                M2 = M2 + delta * (v - mean)
-                if minvalue is None:
-                    minvalue = v
-                    maxvalue = v
-                else:
-                    minvalue = min(v, minvalue)
-                    maxvalue = max(v, maxvalue)
-            else:
-                nodata += 1
-
-        variance = M2 / (n - 1)
-        stddev = math.sqrt(variance)
+        layer = self.parameterAsRasterLayer(parameters, self.INPUT, context)
+        band = self.parameterAsInt(parameters, self.BAND, context)
+        outputFile = self.parameterAsFileOutput(parameters, self.OUTPUT_HTML_FILE, context)
+
+        stat = layer.dataProvider().bandStatistics(band, QgsRasterBandStats.All, QgsRectangle(), 0)
 
         data = []
-        data.append('Valid cells: ' + str(n))
-        data.append('No-data cells: ' + str(nodata))
-        data.append('Minimum value: ' + str(minvalue))
-        data.append('Maximum value: ' + str(maxvalue))
-        data.append('Sum: ' + str(sum))
-        data.append('Mean value: ' + str(mean))
-        data.append('Standard deviation: ' + str(stddev))
-
-        self.createHTML(outputFile, data)
-
-        self.setOutputValue(self.COUNT, n)
-        self.setOutputValue(self.NO_DATA_COUNT, nodata)
-        self.setOutputValue(self.MIN, minvalue)
-        self.setOutputValue(self.MAX, maxvalue)
-        self.setOutputValue(self.SUM, sum)
-        self.setOutputValue(self.MEAN, mean)
-        self.setOutputValue(self.STD_DEV, stddev)
+        data.append(self.tr('Analyzed file: {} (band {})').format(layer.source(), band))
+        data.append(self.tr('Minimum value: {}').format(stat.minimumValue))
+        data.append(self.tr('Maximum value: {}').format(stat.maximumValue))
+        data.append(self.tr('Range: {}').format(stat.range))
+        data.append(self.tr('Sum: {}').format(stat.sum))
+        data.append(self.tr('Mean value: {}').format(stat.mean))
+        data.append(self.tr('Standard deviation: {}').format(stat.stdDev))
+        data.append(self.tr('Sum of the squares: {}').format(stat.sumOfSquares))
+
+        results = {self.MIN: stat.minimumValue,
+                   self.MAX: stat.maximumValue,
+                   self.RANGE: stat.range,
+                   self.SUM: stat.sum,
+                   self.MEAN: stat.mean,
+                   self.STD_DEV: stat.stdDev,
+                   self.SUM_OF_SQUARES: stat.sumOfSquares}
+
+        if outputFile:
+            self.createHTML(outputFile, data)
+            results[self.OUTPUT_HTML_FILE] = outputFile
+
+        return results
 
     def createHTML(self, outputFile, algData):
         with codecs.open(outputFile, 'w', encoding='utf-8') as f:
-            f.write('<html><head>')
+            f.write('<html><head>\n')
             f.write('<meta http-equiv="Content-Type" content="text/html; \
-                    charset=utf-8" /></head><body>')
+                    charset=utf-8" /></head><body>\n')
             for s in algData:
-                f.write('<p>' + str(s) + '</p>')
-            f.write('</body></html>')
+                f.write('<p>' + str(s) + '</p>\n')
+            f.write('</body></html>\n')

python/plugins/processing/tests/testdata/expected/raster_statistics.html

@@ -0,0 +1,11 @@
+<html><head>
+<meta http-equiv="Content-Type" content="text/html;                     charset=utf-8" /></head><body>
+<p>Analyzed file: /home/alex/devel/qgis/python/plugins/processing/tests/testdata/dem.tif (band 1)</p>
+<p>Minimum value: 85.0</p>
+<p>Maximum value: 243.0</p>
+<p>Range: 158.0</p>
+<p>Sum: 19213301.982429504</p>
+<p>Mean value: 147.17197994967066</p>
+<p>Standard deviation: 43.9618116337985</p>
+<p>Sum of the squares: 252304334.52061242</p>
+</body></html>

python/plugins/processing/tests/testdata/qgis_algorithm_tests.yaml

@@ -2758,3 +2758,23 @@ tests:
       OUTPUT:
         hash: f09384c64f56286ec4146a7b9a679cea7c6711ec4c7d77eec054e364
         type: rasterhash
+
+  - algorithm: qgis:rasterlayerstatistics
+    name: Raster layer statistics
+    params:
+      INPUT:
+        name: dem.tif
+        type: raster
+      BAND: 1
+    results:
+      OUTPUT_HTML_FILE:
+        name: raster_statistics.html
+        type: regex
+        rules:
+          - 'Minimum value: 85.0'
+          - 'Maximum value: 243.0'
+          - 'Range: 158.0'
+          - 'Sum: 19213301.982429504'
+          - 'Mean value: 147.17197994967066'
+          - 'Standard deviation: 43.9618116337985'
+          - 'Sum of the squares: 252304334.52061242'