diff --git a/CHANGES.md b/CHANGES.md index a4320a74c..6c37e00bf 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -92,6 +92,10 @@ - +- Added wrapper methods for `GDALGetRasterStatistics`, `GDALComputeRasterMinMax` and `GDALMDArrayGetStatistics`. + + - + ## 0.12 - Bump Rust edition to 2021 diff --git a/src/raster/mdarray.rs b/src/raster/mdarray.rs index 0237c67e5..c54578eda 100644 --- a/src/raster/mdarray.rs +++ b/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> { + 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(); + } } diff --git a/src/raster/mod.rs b/src/raster/mod.rs index 4784a098f..103cabfa2 100644 --- a/src/raster/mod.rs +++ b/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}; diff --git a/src/raster/rasterband.rs b/src/raster/rasterband.rs index dbe8d3521..46d41d9d4 100644 --- a/src/raster/rasterband.rs +++ b/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> { + 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 { + 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> { diff --git a/src/raster/tests.rs b/src/raster/tests.rs index e46913924..3f6a1c722 100644 --- a/src/raster/tests.rs +++ b/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(); +}