From 8f5904560a7ab51d480410ca5228df1a28ec2295 Mon Sep 17 00:00:00 2001 From: Daniel Wiesmann Date: Tue, 24 Nov 2015 14:48:00 +0000 Subject: [PATCH] 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 | 13 ++- django/contrib/gis/gdal/raster/band.py | 84 ++++++++++++++++++-- docs/ref/contrib/gis/gdal.txt | 43 ++++++++++ docs/releases/1.10.txt | 5 ++ tests/gis_tests/gdal_tests/test_raster.py | 52 +++++++++++- 5 files changed, 186 insertions(+), 11 deletions(-) diff --git a/django/contrib/gis/gdal/prototypes/raster.py b/django/contrib/gis/gdal/prototypes/raster.py index 12167b4774..4058e1912a 100644 --- a/django/contrib/gis/gdal/prototypes/raster.py +++ b/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_maximum = double_output(std_call('GDALGetRasterMaximum'), [c_void_p, POINTER(c_int)]) +get_band_statistics = void_output(std_call('GDALGetRasterStatistics'), + [ + 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'), diff --git a/django/contrib/gis/gdal/raster/band.py b/django/contrib/gis/gdal/raster/band.py index bd273b2985..a573ef53c8 100644 --- a/django/contrib/gis/gdal/raster/band.py +++ b/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): diff --git a/docs/ref/contrib/gis/gdal.txt b/docs/ref/contrib/gis/gdal.txt index 2f3cb9f3d9..23728a354a 100644 --- a/docs/ref/contrib/gis/gdal.txt +++ b/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 diff --git a/docs/releases/1.10.txt b/docs/releases/1.10.txt index 26e8849f3f..1c8b0b072d 100644 --- a/docs/releases/1.10.txt +++ b/docs/releases/1.10.txt @@ -74,6 +74,11 @@ Minor features * Added the :meth:`GEOSGeometry.covers() ` binary predicate. +* Added the :meth:`GDALBand.statistics() + ` method and + :attr:`~django.contrib.gis.gdal.GDALBand.mean` + and :attr:`~django.contrib.gis.gdal.GDALBand.std` attributes. + :mod:`django.contrib.messages` ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/tests/gis_tests/gdal_tests/test_raster.py b/tests/gis_tests/gdal_tests/test_raster.py index 8b9021bfc7..3f281d277e 100644 --- a/tests/gis_tests/gdal_tests/test_raster.py +++ b/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))