From e17e4a05eb04259fc364446bc228923052550622 Mon Sep 17 00:00:00 2001 From: Christian Beilschmidt Date: Mon, 29 Aug 2022 10:15:54 +0200 Subject: [PATCH 1/3] get_statistics for rasterbands and md rasters --- src/raster/mdarray.rs | 114 ++++++++++++++++++++++++++++++++++++--- src/raster/mod.rs | 5 +- src/raster/rasterband.rs | 86 ++++++++++++++++++++++++++++- src/raster/tests.rs | 34 +++++++++++- 4 files changed, 226 insertions(+), 13 deletions(-) diff --git a/src/raster/mdarray.rs b/src/raster/mdarray.rs index 3bf5470aa..403bd5dc4 100644 --- a/src/raster/mdarray.rs +++ b/src/raster/mdarray.rs @@ -7,14 +7,15 @@ use gdal_sys::{ CPLErr, CSLDestroy, GDALAttributeGetDataType, GDALAttributeGetDimensionsSize, GDALAttributeH, GDALAttributeReadAsDouble, GDALAttributeReadAsDoubleArray, GDALAttributeReadAsInt, GDALAttributeReadAsIntArray, GDALAttributeReadAsString, GDALAttributeReadAsStringArray, - GDALAttributeRelease, GDALDataType, GDALDimensionGetIndexingVariable, GDALDimensionGetName, - GDALDimensionGetSize, GDALDimensionHS, GDALDimensionRelease, GDALExtendedDataTypeClass, - GDALExtendedDataTypeGetClass, GDALExtendedDataTypeGetNumericDataType, GDALExtendedDataTypeH, - GDALExtendedDataTypeRelease, GDALGroupGetAttribute, GDALGroupGetGroupNames, - GDALGroupGetMDArrayNames, GDALGroupGetName, GDALGroupH, GDALGroupOpenGroup, - GDALGroupOpenMDArray, GDALGroupRelease, GDALMDArrayGetAttribute, GDALMDArrayGetDataType, - GDALMDArrayGetDimensionCount, GDALMDArrayGetDimensions, GDALMDArrayGetNoDataValueAsDouble, - GDALMDArrayGetSpatialRef, GDALMDArrayGetTotalElementsCount, GDALMDArrayGetUnit, GDALMDArrayH, + GDALAttributeRelease, GDALDataType, GDALDatasetH, GDALDimensionGetIndexingVariable, + GDALDimensionGetName, GDALDimensionGetSize, GDALDimensionHS, GDALDimensionRelease, + GDALExtendedDataTypeClass, GDALExtendedDataTypeGetClass, + GDALExtendedDataTypeGetNumericDataType, GDALExtendedDataTypeH, GDALExtendedDataTypeRelease, + GDALGroupGetAttribute, GDALGroupGetGroupNames, GDALGroupGetMDArrayNames, GDALGroupGetName, + GDALGroupH, GDALGroupOpenGroup, GDALGroupOpenMDArray, GDALGroupRelease, + GDALMDArrayGetAttribute, GDALMDArrayGetDataType, GDALMDArrayGetDimensionCount, + GDALMDArrayGetDimensions, GDALMDArrayGetNoDataValueAsDouble, GDALMDArrayGetSpatialRef, + GDALMDArrayGetStatistics, GDALMDArrayGetTotalElementsCount, GDALMDArrayGetUnit, GDALMDArrayH, GDALMDArrayRelease, OSRDestroySpatialReference, VSIFree, }; use libc::c_void; @@ -33,6 +34,7 @@ use std::fmt::Debug; #[derive(Debug)] pub struct MDArray<'a> { c_mdarray: GDALMDArrayH, + c_dataset: GDALDatasetH, _parent: GroupOrDimension<'a>, } @@ -58,6 +60,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 }, } } @@ -72,6 +75,7 @@ impl<'a> MDArray<'a> { ) -> Self { Self { c_mdarray, + c_dataset: _dimension._md_array.c_dataset, _parent: GroupOrDimension::Dimension { _dimension }, } } @@ -341,6 +345,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 @@ -924,4 +988,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 nullptr for 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 db7acf6af..04bc436dc 100644 --- a/src/raster/mod.rs +++ b/src/raster/mod.rs @@ -8,10 +8,11 @@ mod types; mod warp; #[cfg(all(major_ge_3, minor_ge_1))] -pub use mdarray::{ExtendedDataType, Group, MDArray}; +pub use mdarray::{ExtendedDataType, 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 0c8002202..e232e15f0 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; @@ -564,6 +565,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 9d7be8fde..7811922e8 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; @@ -781,3 +783,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(); +} From 89134d2d5fbe90ed726f2f0679f8a768a95532ee Mon Sep 17 00:00:00 2001 From: Christian Beilschmidt Date: Mon, 29 Aug 2022 10:18:05 +0200 Subject: [PATCH 2/3] added changes to CHANGES.md --- CHANGES.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGES.md b/CHANGES.md index 3ede7919d..06ceca9ca 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -3,6 +3,7 @@ ## Unreleased - Add prebuild bindings for GDAL 3.5 + - - **Breaking**: Add `gdal::vector::OwnedLayer`, `gdal::vector::LayerAccess` and `gdal::vector::layer::OwnedFeatureIterator`. This requires importing `gdal::vector::LayerAccess` for using most vector layer methods. @@ -75,6 +76,10 @@ - +- Added wrapper methods for `GDALGetRasterStatistics`, `GDALComputeRasterMinMax` and `GDALMDArrayGetStatistics`. + + - + ## 0.12 - Bump Rust edition to 2021 From b30ae684af73c4cd0c86ccc542c935206904782a Mon Sep 17 00:00:00 2001 From: Christian Beilschmidt Date: Fri, 2 Sep 2022 11:52:34 +0200 Subject: [PATCH 3/3] Update src/raster/mdarray.rs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Laurențiu Nicola --- src/raster/mdarray.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/raster/mdarray.rs b/src/raster/mdarray.rs index b9264d4ac..c54578eda 100644 --- a/src/raster/mdarray.rs +++ b/src/raster/mdarray.rs @@ -1110,7 +1110,7 @@ mod tests { 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 nullptr for default behavior. + 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());