Fixed #25734 -- Made GDALBand min and max properties use GDALComputeRasterStatistics.

Thanks Sergey Fedoseev and Tim Graham for the review.

django/contrib/gis/gdal/prototypes/raster.py

@@ -62,8 +62,17 @@ get_band_ds = voidptr_output(std_call('GDALGetBandDataset'), [c_void_p])
 get_band_datatype = int_output(std_call('GDALGetRasterDataType'), [c_void_p])
 get_band_nodata_value = double_output(std_call('GDALGetRasterNoDataValue'), [c_void_p, POINTER(c_int)])
 set_band_nodata_value = void_output(std_call('GDALSetRasterNoDataValue'), [c_void_p, c_double])
-get_band_minimum = double_output(std_call('GDALGetRasterMinimum'), [c_void_p, POINTER(c_int)])
+get_band_statistics = void_output(std_call('GDALGetRasterStatistics'),
-get_band_maximum = double_output(std_call('GDALGetRasterMaximum'), [c_void_p, POINTER(c_int)])
+    [
+        c_void_p, c_int, c_int, POINTER(c_double), POINTER(c_double),
+        POINTER(c_double), POINTER(c_double), c_void_p, c_void_p,
+    ],
+    errcheck=False
+)
+compute_band_statistics = void_output(std_call('GDALComputeRasterStatistics'),
+    [c_void_p, c_int, POINTER(c_double), POINTER(c_double), POINTER(c_double), POINTER(c_double), c_void_p, c_void_p],
+    errcheck=False
+)
 
 # Reprojection routine
 reproject_image = void_output(std_call('GDALReprojectImage'),

django/contrib/gis/gdal/raster/band.py

@@ -1,4 +1,5 @@
-from ctypes import byref, c_int
+import math
+from ctypes import byref, c_double, c_int, c_void_p
 
 from django.contrib.gis.gdal.base import GDALBase
 from django.contrib.gis.gdal.error import GDALException
@@ -19,6 +20,14 @@ class GDALBand(GDALBase):
         self.source = source
         self._ptr = capi.get_ds_raster_band(source._ptr, index)
 
+    def _flush(self):
+        """
+        Call the flush method on the Band's parent raster and force a refresh
+        of the statistics attribute when requested the next time.
+        """
+        self.source._flush()
+        self._stats_refresh = True
+
     @property
     def description(self):
         """
@@ -47,19 +56,80 @@ class GDALBand(GDALBase):
         """
         return self.width * self.height
 
+    _stats_refresh = False
+
+    def statistics(self, refresh=False, approximate=False):
+        """
+        Compute statistics on the pixel values of this band.
+
+        The return value is a tuple with the following structure:
+        (minimum, maximum, mean, standard deviation).
+
+        If approximate=True, the statistics may be computed based on overviews
+        or a subset of image tiles.
+
+        If refresh=True, the statistics will be computed from the data directly,
+        and the cache will be updated where applicable.
+
+        For empty bands (where all pixel values are nodata), all statistics
+        values are returned as None.
+
+        For raster formats using Persistent Auxiliary Metadata (PAM) services,
+        the statistics might be cached in an auxiliary file.
+        """
+        # Prepare array with arguments for capi function
+        smin, smax, smean, sstd = c_double(), c_double(), c_double(), c_double()
+        stats_args = [
+            self._ptr, c_int(approximate), byref(smin), byref(smax),
+            byref(smean), byref(sstd), c_void_p(), c_void_p(),
+        ]
+
+        if refresh or self._stats_refresh:
+            capi.compute_band_statistics(*stats_args)
+        else:
+            # Add additional argument to force computation if there is no
+            # existing PAM file to take the values from.
+            force = True
+            stats_args.insert(2, c_int(force))
+            capi.get_band_statistics(*stats_args)
+
+        result = smin.value, smax.value, smean.value, sstd.value
+
+        # Check if band is empty (in that case, set all statistics to None)
+        if any((math.isnan(val) for val in result)):
+            result = (None, None, None, None)
+
+        self._stats_refresh = False
+
+        return result
+
     @property
     def min(self):
         """
-        Returns the minimum pixel value for this band.
+        Return the minimum pixel value for this band.
         """
-        return capi.get_band_minimum(self._ptr, byref(c_int()))
+        return self.statistics()[0]
 
     @property
     def max(self):
         """
-        Returns the maximum pixel value for this band.
+        Return the maximum pixel value for this band.
         """
-        return capi.get_band_maximum(self._ptr, byref(c_int()))
+        return self.statistics()[1]
+
+    @property
+    def mean(self):
+        """
+        Return the mean of all pixel values of this band.
+        """
+        return self.statistics()[2]
+
+    @property
+    def std(self):
+        """
+        Return the standard deviation of all pixel values of this band.
+        """
+        return self.statistics()[3]
 
     @property
     def nodata_value(self):
@@ -84,7 +154,7 @@ class GDALBand(GDALBase):
         if not isinstance(value, (int, float)):
             raise ValueError('Nodata value must be numeric.')
         capi.set_band_nodata_value(self._ptr, value)
-        self.source._flush()
+        self._flush()
 
     def datatype(self, as_string=False):
         """
@@ -149,7 +219,7 @@ class GDALBand(GDALBase):
             else:
                 return list(data_array)
         else:
-            self.source._flush()
+            self._flush()
 
 
 class BandList(list):

docs/ref/contrib/gis/gdal.txt

@@ -1413,6 +1413,35 @@ blue.
 
         The total number of pixels in this band. Is equal to ``width * height``.
 
+    .. method:: statistics(refresh=False, approximate=False)
+
+        .. versionadded:: 1.10
+
+        Compute statistics on the pixel values of this band. The return value
+        is a tuple with the following structure:
+        ``(minimum, maximum, mean, standard deviation)``.
+
+        If the ``approximate`` argument is set to ``True``, the statistics may
+        be computed based on overviews or a subset of image tiles.
+
+        If the ``refresh`` argument is set to ``True``, the statistics will be
+        computed from the data directly, and the cache will be updated with the
+        result.
+
+        If a persistent cache value is found, that value is returned. For
+        raster formats using Persistent Auxiliary Metadata (PAM) services, the
+        statistics might be cached in an auxiliary file. In some cases this
+        metadata might be out of sync with the pixel values or cause values
+        from a previous call to be returned which don't reflect the value of
+        the ``approximate`` argument. In such cases, use the ``refresh``
+        argument to get updated values and store them in the cache.
+
+        For empty bands (where all pixel values are "no data"), all statistics
+        are returned as ``None``.
+
+        The statistics can also be retrieved directly by accessing the
+        :attr:`min`, :attr:`max`, :attr:`mean`, and :attr:`std` properties.
+
     .. attribute:: min
 
         The minimum pixel value of the band (excluding the "no data" value).
@@ -1421,6 +1450,20 @@ blue.
 
         The maximum pixel value of the band (excluding the "no data" value).
 
+    .. attribute:: mean
+
+        .. versionadded:: 1.10
+
+        The mean of all pixel values of the band (excluding the "no data"
+        value).
+
+    .. attribute:: std
+
+        .. versionadded:: 1.10
+
+        The standard deviation of all pixel values of the band (excluding the
+        "no data" value).
+
     .. attribute:: nodata_value
 
         The "no data" value for a band is generally a special marker value used

docs/releases/1.10.txt

@@ -74,6 +74,11 @@ Minor features
 * Added the :meth:`GEOSGeometry.covers()
   <django.contrib.gis.geos.GEOSGeometry.covers>` binary predicate.
 
+* Added the :meth:`GDALBand.statistics()
+  <django.contrib.gis.gdal.GDALBand.statistics>` method and
+  :attr:`~django.contrib.gis.gdal.GDALBand.mean`
+  and :attr:`~django.contrib.gis.gdal.GDALBand.std` attributes.
+
 :mod:`django.contrib.messages`
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

tests/gis_tests/gdal_tests/test_raster.py

@@ -310,14 +310,34 @@ class GDALBandTests(unittest.TestCase):
         self.band = rs.bands[0]
 
     def test_band_data(self):
+        pam_file = self.rs_path + '.aux.xml'
         self.assertEqual(self.band.width, 163)
         self.assertEqual(self.band.height, 174)
         self.assertEqual(self.band.description, '')
         self.assertEqual(self.band.datatype(), 1)
         self.assertEqual(self.band.datatype(as_string=True), 'GDT_Byte')
-        self.assertEqual(self.band.min, 0)
-        self.assertEqual(self.band.max, 255)
         self.assertEqual(self.band.nodata_value, 15)
+        try:
+            self.assertEqual(
+                self.band.statistics(approximate=True),
+                (0.0, 9.0, 2.842331288343558, 2.3965567248965356)
+            )
+            self.assertEqual(
+                self.band.statistics(approximate=False, refresh=True),
+                (0.0, 9.0, 2.828326634228898, 2.4260526986669095)
+            )
+            self.assertEqual(self.band.min, 0)
+            self.assertEqual(self.band.max, 9)
+            self.assertEqual(self.band.mean, 2.8283266342289)
+            self.assertEqual(self.band.std, 2.4260526986669)
+            # Check that statistics are persisted into PAM file on band close
+            self.band = None
+            self.assertTrue(os.path.isfile(pam_file))
+        finally:
+            # Close band and remove file if created
+            self.band = None
+            if os.path.isfile(pam_file):
+                os.remove(pam_file)
 
     def test_read_mode_error(self):
         # Open raster in read mode
@@ -413,3 +433,31 @@ class GDALBandTests(unittest.TestCase):
             )
         else:
             self.assertEqual(bandmemjson.data(), list(range(25)))
+
+    def test_band_statistics_automatic_refresh(self):
+        rsmem = GDALRaster({
+            'srid': 4326,
+            'width': 2,
+            'height': 2,
+            'bands': [{'data': [0] * 4, 'nodata_value': 99}],
+        })
+        band = rsmem.bands[0]
+        # Populate statistics cache
+        self.assertEqual(band.statistics(), (0, 0, 0, 0))
+        # Change data
+        band.data([1, 1, 0, 0])
+        # Statistics are properly updated
+        self.assertEqual(band.statistics(), (0.0, 1.0, 0.5, 0.5))
+        # Change nodata_value
+        band.nodata_value = 0
+        # Statistics are properly updated
+        self.assertEqual(band.statistics(), (1.0, 1.0, 1.0, 0.0))
+
+    def test_band_statistics_empty_band(self):
+        rsmem = GDALRaster({
+            'srid': 4326,
+            'width': 1,
+            'height': 1,
+            'bands': [{'data': [0], 'nodata_value': 0}],
+        })
+        self.assertEqual(rsmem.bands[0].statistics(), (None, None, None, None))