Merge #292

292: get_statistics for rasterbands and md rasters r=lnicola a=ChristianBeilschmidt

- [X] I agree to follow the project's [code of conduct](https://github.com/georust/gdal/blob/master/CODE_OF_CONDUCT.md).
- [x] I added an entry to `CHANGES.md` if knowledge of this change could be valuable to users.
---



Co-authored-by: Christian Beilschmidt <christian.beilschmidt@geoengine.de>

CHANGES.md

@@ -92,6 +92,10 @@
 
   - <https://github.com/georust/gdal/pull/291>
 
+- Added wrapper methods for `GDALGetRasterStatistics`, `GDALComputeRasterMinMax` and `GDALMDArrayGetStatistics`.
+
+  - <https://github.com/georust/gdal/pull/292>
+
 ## 0.12
 
 - Bump Rust edition to 2021

src/raster/mdarray.rs

@@ -7,16 +7,16 @@ use gdal_sys::{
     CPLErr, CSLDestroy, GDALAttributeGetDataType, GDALAttributeGetDimensionsSize, GDALAttributeH,
     GDALAttributeReadAsDouble, GDALAttributeReadAsDoubleArray, GDALAttributeReadAsInt,
     GDALAttributeReadAsIntArray, GDALAttributeReadAsString, GDALAttributeReadAsStringArray,
-    GDALAttributeRelease, GDALDataType, GDALDimensionGetIndexingVariable, GDALDimensionGetName,
-    GDALDimensionGetSize, GDALDimensionHS, GDALDimensionRelease, GDALExtendedDataTypeClass,
-    GDALExtendedDataTypeGetClass, GDALExtendedDataTypeGetName,
+    GDALAttributeRelease, GDALDataType, GDALDatasetH, GDALDimensionGetIndexingVariable,
+    GDALDimensionGetName, GDALDimensionGetSize, GDALDimensionHS, GDALDimensionRelease,
+    GDALExtendedDataTypeClass, GDALExtendedDataTypeGetClass, GDALExtendedDataTypeGetName,
     GDALExtendedDataTypeGetNumericDataType, GDALExtendedDataTypeH, GDALExtendedDataTypeRelease,
     GDALGroupGetAttribute, GDALGroupGetDimensions, GDALGroupGetGroupNames,
     GDALGroupGetMDArrayNames, GDALGroupGetName, GDALGroupH, GDALGroupOpenGroup,
     GDALGroupOpenMDArray, GDALGroupRelease, GDALMDArrayGetAttribute, GDALMDArrayGetDataType,
     GDALMDArrayGetDimensionCount, GDALMDArrayGetDimensions, GDALMDArrayGetNoDataValueAsDouble,
-    GDALMDArrayGetSpatialRef, GDALMDArrayGetTotalElementsCount, GDALMDArrayGetUnit, GDALMDArrayH,
-    GDALMDArrayRelease, OSRDestroySpatialReference, VSIFree,
+    GDALMDArrayGetSpatialRef, GDALMDArrayGetStatistics, GDALMDArrayGetTotalElementsCount,
+    GDALMDArrayGetUnit, GDALMDArrayH, GDALMDArrayRelease, OSRDestroySpatialReference, VSIFree,
 };
 use libc::c_void;
 use std::ffi::CString;
@@ -34,6 +34,7 @@ use std::fmt::{Debug, Display};
 #[derive(Debug)]
 pub struct MDArray<'a> {
     c_mdarray: GDALMDArrayH,
+    c_dataset: GDALDatasetH,
     _parent: GroupOrDimension<'a>,
 }
 
@@ -65,6 +66,7 @@ impl<'a> MDArray<'a> {
     pub unsafe fn from_c_mdarray_and_group(_group: &'a Group, c_mdarray: GDALMDArrayH) -> Self {
         Self {
             c_mdarray,
+            c_dataset: _group._dataset.c_dataset(),
             _parent: GroupOrDimension::Group { _group },
         }
     }
@@ -79,6 +81,10 @@ impl<'a> MDArray<'a> {
     ) -> Self {
         Self {
             c_mdarray,
+            c_dataset: match _dimension._parent {
+                GroupOrArray::Group { _group } => _group._dataset.c_dataset(),
+                GroupOrArray::MDArray { _md_array } => _md_array.c_dataset,
+            },
             _parent: GroupOrDimension::Dimension { _dimension },
         }
     }
@@ -352,6 +358,66 @@ impl<'a> MDArray<'a> {
             Ok(Attribute::from_c_attribute(c_attribute))
         }
     }
+
+    /// Fetch statistics.
+    ///
+    /// Returns the minimum, maximum, mean and standard deviation of all pixel values in this array.
+    ///
+    /// If `force` is `false` results will only be returned if it can be done quickly (i.e. without scanning the data).
+    /// If `force` is `false` and results cannot be returned efficiently, the method will return `None`.
+    ///
+    /// When cached statistics are not available, and `force` is `true`, ComputeStatistics() is called.
+    ///
+    /// Note that file formats using PAM (Persistent Auxiliary Metadata) services will generally cache statistics in the .aux.xml file allowing fast fetch after the first request.
+    ///
+    /// This methods is a wrapper for [`GDALMDArrayGetStatistics`](https://gdal.org/api/gdalmdarray_cpp.html#_CPPv4N11GDALMDArray13GetStatisticsEbbPdPdPdPdP7GUInt6416GDALProgressFuncPv).
+    ///
+    /// TODO: add option to pass progress callback (`pfnProgress`)
+    ///
+    pub fn get_statistics(
+        &self,
+        force: bool,
+        is_approx_ok: bool,
+    ) -> Result<Option<MdStatisticsAll>> {
+        let mut statistics = MdStatisticsAll {
+            min: 0.,
+            max: 0.,
+            mean: 0.,
+            std_dev: 0.,
+            valid_count: 0,
+        };
+
+        let rv = unsafe {
+            GDALMDArrayGetStatistics(
+                self.c_mdarray,
+                self.c_dataset,
+                libc::c_int::from(is_approx_ok),
+                libc::c_int::from(force),
+                &mut statistics.min,
+                &mut statistics.max,
+                &mut statistics.mean,
+                &mut statistics.std_dev,
+                &mut statistics.valid_count,
+                None,
+                std::ptr::null_mut(),
+            )
+        };
+
+        match CplErrType::from(rv) {
+            CplErrType::None => Ok(Some(statistics)),
+            CplErrType::Warning => Ok(None),
+            _ => Err(_last_cpl_err(rv)),
+        }
+    }
+}
+
+#[derive(Debug, PartialEq)]
+pub struct MdStatisticsAll {
+    pub min: f64,
+    pub max: f64,
+    pub mean: f64,
+    pub std_dev: f64,
+    pub valid_count: u64,
 }
 
 /// Represent a mdarray in a dataset
@@ -1032,4 +1098,38 @@ mod tests {
 
         assert_eq!(md_array.unit(), "K");
     }
+
+    #[test]
+    fn test_stats() {
+        let dataset_options = DatasetOptions {
+            open_flags: GdalOpenFlags::GDAL_OF_MULTIDIM_RASTER,
+            allowed_drivers: None,
+            open_options: None,
+            sibling_files: None,
+        };
+        let dataset = Dataset::open_ex("fixtures/byte_no_cf.nc", dataset_options).unwrap();
+        let root_group = dataset.root_group().unwrap();
+        let array_name = "Band1".to_string();
+        let options = CslStringList::new(); //Driver specific options determining how the array should be opened. Pass an empty one for the default behavior.
+        let md_array = root_group.open_md_array(&array_name, options).unwrap();
+
+        assert!(md_array.get_statistics(false, true).unwrap().is_none());
+
+        assert_eq!(
+            md_array.get_statistics(true, true).unwrap().unwrap(),
+            MdStatisticsAll {
+                min: 74.0,
+                max: 255.0,
+                mean: 126.76500000000001,
+                std_dev: 22.928470838675654,
+                valid_count: 400,
+            }
+        );
+
+        // clean up aux file
+        drop(md_array);
+        drop(root_group);
+        drop(dataset);
+        std::fs::remove_file("fixtures/byte_no_cf.nc.aux.xml").unwrap();
+    }
 }

src/raster/mod.rs

@@ -8,10 +8,13 @@ mod types;
 mod warp;
 
 #[cfg(all(major_ge_3, minor_ge_1))]
-pub use mdarray::{Attribute, Dimension, ExtendedDataType, ExtendedDataTypeClass, Group, MDArray};
+pub use mdarray::{
+    Attribute, Dimension, ExtendedDataType, ExtendedDataTypeClass, Group, MDArray, MdStatisticsAll,
+};
 pub use rasterband::{
     Buffer, ByteBuffer, CmykEntry, ColorEntry, ColorInterpretation, ColorTable, GrayEntry,
-    HlsEntry, PaletteInterpretation, RasterBand, ResampleAlg, RgbaEntry,
+    HlsEntry, PaletteInterpretation, RasterBand, ResampleAlg, RgbaEntry, StatisticsAll,
+    StatisticsMinMax,
 };
 pub use rasterize::{rasterize, BurnSource, MergeAlgorithm, OptimizeMode, RasterizeOptions};
 pub use types::{GDALDataType, GdalType};

src/raster/rasterband.rs

@@ -4,8 +4,9 @@ use crate::metadata::Metadata;
 use crate::raster::{GDALDataType, GdalType};
 use crate::utils::{_last_cpl_err, _last_null_pointer_err, _string};
 use gdal_sys::{
-    self, CPLErr, GDALColorEntry, GDALColorInterp, GDALColorTableH, GDALMajorObjectH,
-    GDALPaletteInterp, GDALRWFlag, GDALRasterBandH, GDALRasterIOExtraArg,
+    self, CPLErr, GDALColorEntry, GDALColorInterp, GDALColorTableH, GDALComputeRasterMinMax,
+    GDALGetRasterStatistics, GDALMajorObjectH, GDALPaletteInterp, GDALRWFlag, GDALRasterBandH,
+    GDALRasterIOExtraArg,
 };
 use libc::c_int;
 use std::ffi::CString;
@@ -582,6 +583,87 @@ impl<'a> RasterBand<'a> {
             Ok(mask_band)
         }
     }
+
+    /// Fetch image statistics.
+    ///
+    /// Returns the minimum, maximum, mean and standard deviation of all pixel values in this band.
+    /// If approximate statistics are sufficient, the `is_approx_ok` flag can be set to true in which case overviews, or a subset of image tiles may be used in computing the statistics.
+    ///
+    /// If `force` is `false` results will only be returned if it can be done quickly (i.e. without scanning the data).
+    /// If force` is `false` and results cannot be returned efficiently, the method will return `None`.
+    ///
+    /// Note that file formats using PAM (Persistent Auxiliary Metadata) services will generally cache statistics in the .pam file allowing fast fetch after the first request.
+    ///
+    /// This methods is a wrapper for [`GDALGetRasterStatistics`](https://gdal.org/api/gdalrasterband_cpp.html#_CPPv4N14GDALRasterBand13GetStatisticsEiiPdPdPdPd).
+    ///
+    pub fn get_statistics(&self, force: bool, is_approx_ok: bool) -> Result<Option<StatisticsAll>> {
+        let mut statistics = StatisticsAll {
+            min: 0.,
+            max: 0.,
+            mean: 0.,
+            std_dev: 0.,
+        };
+
+        let rv = unsafe {
+            GDALGetRasterStatistics(
+                self.c_rasterband,
+                libc::c_int::from(is_approx_ok),
+                libc::c_int::from(force),
+                &mut statistics.min,
+                &mut statistics.max,
+                &mut statistics.mean,
+                &mut statistics.std_dev,
+            )
+        };
+
+        match CplErrType::from(rv) {
+            CplErrType::None => Ok(Some(statistics)),
+            CplErrType::Warning => Ok(None),
+            _ => Err(_last_cpl_err(rv)),
+        }
+    }
+
+    /// Compute the min/max values for a band.
+    ///
+    /// If `is_approx_ok` is `true`, then the band’s GetMinimum()/GetMaximum() will be trusted.
+    /// If it doesn’t work, a subsample of blocks will be read to get an approximate min/max.
+    /// If the band has a nodata value it will be excluded from the minimum and maximum.
+    ///
+    /// If `is_approx_ok` is `false`, then all pixels will be read and used to compute an exact range.
+    ///
+    /// This methods is a wrapper for [`GDALComputeRasterMinMax`](https://gdal.org/api/gdalrasterband_cpp.html#_CPPv4N14GDALRasterBand19ComputeRasterMinMaxEiPd).
+    ///
+    pub fn compute_raster_min_max(&self, is_approx_ok: bool) -> Result<StatisticsMinMax> {
+        let mut min_max = [0., 0.];
+
+        // TODO: The C++ method actually returns a CPLErr, but the C interface does not expose it.
+        unsafe {
+            GDALComputeRasterMinMax(
+                self.c_rasterband,
+                libc::c_int::from(is_approx_ok),
+                &mut min_max as *mut f64,
+            )
+        };
+
+        Ok(StatisticsMinMax {
+            min: min_max[0],
+            max: min_max[1],
+        })
+    }
+}
+
+#[derive(Debug, PartialEq)]
+pub struct StatisticsMinMax {
+    pub min: f64,
+    pub max: f64,
+}
+
+#[derive(Debug, PartialEq)]
+pub struct StatisticsAll {
+    pub min: f64,
+    pub max: f64,
+    pub mean: f64,
+    pub std_dev: f64,
 }
 
 impl<'a> MajorObject for RasterBand<'a> {

src/raster/tests.rs

@@ -1,7 +1,9 @@
 use crate::dataset::Dataset;
 use crate::metadata::Metadata;
 use crate::raster::rasterband::ResampleAlg;
-use crate::raster::{ByteBuffer, ColorInterpretation, RasterCreationOption};
+use crate::raster::{
+    ByteBuffer, ColorInterpretation, RasterCreationOption, StatisticsAll, StatisticsMinMax,
+};
 use crate::vsi::unlink_mem_file;
 use crate::Driver;
 use gdal_sys::GDALDataType;
@@ -801,3 +803,33 @@ fn test_color_table() {
         }
     }
 }
+
+#[test]
+fn test_raster_stats() {
+    let dataset = Dataset::open(fixture!("tinymarble.tif")).unwrap();
+    let rb = dataset.rasterband(1).unwrap();
+
+    assert!(rb.get_statistics(false, false).unwrap().is_none());
+
+    assert_eq!(
+        rb.get_statistics(true, false).unwrap().unwrap(),
+        StatisticsAll {
+            min: 12.0,
+            max: 255.0,
+            mean: 89.2526,
+            std_dev: 90.99835379412092,
+        }
+    );
+
+    assert_eq!(
+        rb.compute_raster_min_max(true).unwrap(),
+        StatisticsMinMax {
+            min: 12.0,
+            max: 255.0,
+        }
+    );
+
+    // clean up aux file
+    drop(dataset);
+    std::fs::remove_file(fixture!("tinymarble.tif.aux.xml")).unwrap();
+}