From a38e30e19668898317413f7a61b3fab9e21fc4bb Mon Sep 17 00:00:00 2001 From: "Simeon H.K. Fitch" Date: Wed, 25 Oct 2023 09:13:29 -0400 Subject: [PATCH] Added Topographic Position Index processor. --- CHANGES.md | 4 + src/raster/processing/dem/aspect.rs | 28 ++--- src/raster/processing/dem/mod.rs | 10 +- src/raster/processing/dem/roughness.rs | 18 +-- src/raster/processing/dem/slope.rs | 30 ++--- src/raster/processing/dem/tpi.rs | 166 +++++++++++++++++++++++++ 6 files changed, 214 insertions(+), 42 deletions(-) create mode 100644 src/raster/processing/dem/tpi.rs diff --git a/CHANGES.md b/CHANGES.md index b4f81689f..ddef1ba60 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -2,6 +2,10 @@ ## Unreleased + - Added support for digital elevation model raster processing. + + - <> + - **Breaking**: `CslStringListIterator` returns a `CslStringListEntry` instead of `(String, String)` in order to differentiate between `key=value` entries vs `flag` entries. - **Breaking**: `CslStringList::fetch_name_value` returns `Option` instead of `Result>`, better reflecting the semantics of GDAL C API. - Added `CslStringList::get_field`, `CslStringList::find_string`, `CslStringList::partial_find_string`, `CslStringList::find_string_case_sensitive`, `CslStringList::into_ptr` diff --git a/src/raster/processing/dem/aspect.rs b/src/raster/processing/dem/aspect.rs index 7dc94a300..467763980 100644 --- a/src/raster/processing/dem/aspect.rs +++ b/src/raster/processing/dem/aspect.rs @@ -3,7 +3,7 @@ use std::path::{Path, PathBuf}; use crate::cpl::CslStringList; use crate::errors::*; -use crate::raster::processing::dem::{DEMSlopeAlg, GdalDemProcessor}; +use crate::raster::processing::dem::{DemSlopeAlg, GdalDemProcessor}; use crate::raster::processing::{merge, GdalProcessor}; use crate::Dataset; @@ -18,15 +18,15 @@ use crate::Dataset; /// * 270° it's facing the West (provided that the top of your input raster is north oriented). /// /// By default, the aspect value `-9999` is used as the no-data value to indicate undefined aspect in flat -/// areas with slope=0. See [`DEMAspectProcessor::with_zero_for_flat`] for alternative. +/// areas with slope=0. See [`DemAspectProcessor::with_zero_for_flat`] for alternative. /// /// Note: Results are nonsensical if the underlying [`Dataset`] does not contain digital elevation data. /// /// See: [`gdaldem aspect`](https://gdal.org/programs/gdaldem.html#aspect) for details, #[derive(Debug, Clone)] -pub struct DEMAspectProcessor { +pub struct DemAspectProcessor { input_band: Option, - algorithm: Option, + algorithm: Option, destination: Option, compute_edges: Option, zero_for_flat: Option, @@ -34,7 +34,7 @@ pub struct DEMAspectProcessor { extra_options: Option, } -impl DEMAspectProcessor { +impl DemAspectProcessor { /// Create a aspect-from-DEM processor pub fn new() -> Self { Self { @@ -59,7 +59,7 @@ impl DEMAspectProcessor { } /// Specify the slope computation algorithm. - pub fn with_algorithm(self, algorithm: DEMSlopeAlg) -> Self { + pub fn with_algorithm(self, algorithm: DemSlopeAlg) -> Self { Self { algorithm: Some(algorithm), ..self @@ -114,13 +114,13 @@ impl DEMAspectProcessor { } } -impl Default for DEMAspectProcessor { +impl Default for DemAspectProcessor { fn default() -> Self { Self::new() } } -impl GdalProcessor for DEMAspectProcessor { +impl GdalProcessor for DemAspectProcessor { fn to_options(&self) -> CslStringList { let mut opts = CslStringList::new(); @@ -158,7 +158,7 @@ impl GdalProcessor for DEMAspectProcessor { } } -impl GdalDemProcessor for DEMAspectProcessor { +impl GdalDemProcessor for DemAspectProcessor { fn mode(&self) -> &'static str { "aspect" } @@ -179,9 +179,9 @@ mod tests { #[test] fn options() -> Result<()> { - let dsp = DEMAspectProcessor::new() + let dsp = DemAspectProcessor::new() .with_input_band(2.try_into().unwrap()) - .with_algorithm(DEMSlopeAlg::ZevenbergenThorne) + .with_algorithm(DemSlopeAlg::ZevenbergenThorne) .with_compute_edges(true) .with_zero_for_flat(true) .with_trigonometric_angles(true) @@ -197,13 +197,13 @@ mod tests { #[test] fn aspect() -> Result<()> { - let dsp = DEMAspectProcessor::new() - .with_algorithm(DEMSlopeAlg::Horn) + let proc = DemAspectProcessor::new() + .with_algorithm(DemSlopeAlg::Horn) .with_zero_for_flat(true); let ds = Dataset::open(fixture("dem-hills.tiff"))?; - let slope = dsp.eval(&ds)?; + let slope = proc.eval(&ds)?; let stats = slope.rasterband(1)?.get_statistics(true, false)?.unwrap(); diff --git a/src/raster/processing/dem/mod.rs b/src/raster/processing/dem/mod.rs index 0f8d2e038..27bc31a08 100644 --- a/src/raster/processing/dem/mod.rs +++ b/src/raster/processing/dem/mod.rs @@ -1,10 +1,12 @@ mod aspect; mod roughness; mod slope; +mod tpi; pub use aspect::*; pub use roughness::*; pub use slope::*; +pub use tpi::*; use crate::cpl::CslStringList; use crate::errors::Result; @@ -64,16 +66,16 @@ pub(crate) trait GdalDemProcessor: GdalProcessor { /// The literature suggests `ZevenbergenThorne` to be more suited to smooth landscapes, /// whereas `Horn` performs better on rougher terrain. #[derive(Debug, Clone, Copy)] -pub enum DEMSlopeAlg { +pub enum DemSlopeAlg { Horn, ZevenbergenThorne, } -impl Display for DEMSlopeAlg { +impl Display for DemSlopeAlg { fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { match self { - DEMSlopeAlg::Horn => f.write_str("Horn"), - DEMSlopeAlg::ZevenbergenThorne => f.write_str("ZevenbergenThorne"), + DemSlopeAlg::Horn => f.write_str("Horn"), + DemSlopeAlg::ZevenbergenThorne => f.write_str("ZevenbergenThorne"), } } } diff --git a/src/raster/processing/dem/roughness.rs b/src/raster/processing/dem/roughness.rs index 23a87e8b4..cf73b1293 100644 --- a/src/raster/processing/dem/roughness.rs +++ b/src/raster/processing/dem/roughness.rs @@ -19,14 +19,14 @@ use crate::Dataset; /// /// See: [`gdaldem roughness`](https://gdal.org/programs/gdaldem.html#roughness) for details. #[derive(Debug, Clone)] -pub struct DEMRoughnessProcessor { +pub struct DemRoughnessProcessor { input_band: Option, destination: Option, compute_edges: Option, extra_options: Option, } -impl DEMRoughnessProcessor { +impl DemRoughnessProcessor { /// Create a aspect-from-DEM processor pub fn new() -> Self { Self { @@ -77,13 +77,13 @@ impl DEMRoughnessProcessor { } } -impl Default for DEMRoughnessProcessor { +impl Default for DemRoughnessProcessor { fn default() -> Self { Self::new() } } -impl GdalProcessor for DEMRoughnessProcessor { +impl GdalProcessor for DemRoughnessProcessor { fn to_options(&self) -> CslStringList { let mut opts = CslStringList::new(); @@ -108,7 +108,7 @@ impl GdalProcessor for DEMRoughnessProcessor { } } -impl GdalDemProcessor for DEMRoughnessProcessor { +impl GdalDemProcessor for DemRoughnessProcessor { fn mode(&self) -> &'static str { "roughness" } @@ -129,7 +129,7 @@ mod tests { #[test] fn options() -> Result<()> { - let dsp = DEMRoughnessProcessor::new() + let dsp = DemRoughnessProcessor::new() .with_input_band(2.try_into().unwrap()) .with_compute_edges(true) .with_options("-of GTiff".parse()?); @@ -141,12 +141,12 @@ mod tests { } #[test] - fn aspect() -> Result<()> { - let dsp = DEMRoughnessProcessor::new(); + fn roughness() -> Result<()> { + let proc = DemRoughnessProcessor::new(); let ds = Dataset::open(fixture("dem-hills.tiff"))?; - let slope = dsp.eval(&ds)?; + let slope = proc.eval(&ds)?; let stats = slope.rasterband(1)?.get_statistics(true, false)?.unwrap(); diff --git a/src/raster/processing/dem/slope.rs b/src/raster/processing/dem/slope.rs index b9f1e240e..259c27c56 100644 --- a/src/raster/processing/dem/slope.rs +++ b/src/raster/processing/dem/slope.rs @@ -3,7 +3,7 @@ use std::path::{Path, PathBuf}; use crate::cpl::CslStringList; use crate::errors::*; -use crate::raster::processing::dem::{DEMSlopeAlg, GdalDemProcessor}; +use crate::raster::processing::dem::{DemSlopeAlg, GdalDemProcessor}; use crate::raster::processing::{merge, GdalProcessor}; use crate::Dataset; @@ -11,10 +11,10 @@ use crate::Dataset; /// /// This processor will take a DEM raster and output a 32-bit float raster with slope values. /// You have the option of specifying the type of slope value you want: -/// [degrees or percent slope](DEMSlopeProcessor::with_percentage_results). +/// [degrees or percent slope](DemSlopeProcessor::with_percentage_results). /// /// In cases where the horizontal units differ from the vertical units, you can also supply -/// a [scaling factor](DEMSlopeProcessor::with_scale). +/// a [scaling factor](DemSlopeProcessor::with_scale). /// /// The value `-9999` is used as the output no-data value. /// @@ -22,9 +22,9 @@ use crate::Dataset; /// /// See: [`gdaldem slope`](https://gdal.org/programs/gdaldem.html#slope) for details, #[derive(Debug, Clone)] -pub struct DEMSlopeProcessor { +pub struct DemSlopeProcessor { input_band: Option, - algorithm: Option, + algorithm: Option, scale: Option, destination: Option, compute_edges: Option, @@ -32,7 +32,7 @@ pub struct DEMSlopeProcessor { extra_options: Option, } -impl DEMSlopeProcessor { +impl DemSlopeProcessor { /// Create a slope-from-DEM processor pub fn new() -> Self { Self { @@ -57,7 +57,7 @@ impl DEMSlopeProcessor { } /// Specify the slope computation algorithm. - pub fn with_algorithm(self, algorithm: DEMSlopeAlg) -> Self { + pub fn with_algorithm(self, algorithm: DemSlopeAlg) -> Self { Self { algorithm: Some(algorithm), ..self @@ -122,13 +122,13 @@ impl DEMSlopeProcessor { } } -impl Default for DEMSlopeProcessor { +impl Default for DemSlopeProcessor { fn default() -> Self { Self::new() } } -impl GdalProcessor for DEMSlopeProcessor { +impl GdalProcessor for DemSlopeProcessor { fn to_options(&self) -> CslStringList { let mut opts = CslStringList::new(); @@ -167,7 +167,7 @@ impl GdalProcessor for DEMSlopeProcessor { } } -impl GdalDemProcessor for DEMSlopeProcessor { +impl GdalDemProcessor for DemSlopeProcessor { fn mode(&self) -> &'static str { "slope" } @@ -188,9 +188,9 @@ mod tests { #[test] fn options() -> Result<()> { - let dsp = DEMSlopeProcessor::new() + let dsp = DemSlopeProcessor::new() .with_input_band(2.try_into().unwrap()) - .with_algorithm(DEMSlopeAlg::ZevenbergenThorne) + .with_algorithm(DemSlopeAlg::ZevenbergenThorne) .with_scale(37.0) .with_compute_edges(true) .with_percentage_results(true) @@ -205,14 +205,14 @@ mod tests { #[test] fn slope() -> Result<()> { - let dsp = DEMSlopeProcessor::new() - .with_algorithm(DEMSlopeAlg::Horn) + let proc = DemSlopeProcessor::new() + .with_algorithm(DemSlopeAlg::Horn) .with_percentage_results(true) .with_scale(111120.0); let ds = Dataset::open(fixture("dem-hills.tiff"))?; - let slope = dsp.eval(&ds)?; + let slope = proc.eval(&ds)?; let stats = slope.rasterband(1)?.get_statistics(true, false)?.unwrap(); diff --git a/src/raster/processing/dem/tpi.rs b/src/raster/processing/dem/tpi.rs new file mode 100644 index 000000000..5f495d580 --- /dev/null +++ b/src/raster/processing/dem/tpi.rs @@ -0,0 +1,166 @@ +use std::num::NonZeroUsize; +use std::path::{Path, PathBuf}; + +use crate::cpl::CslStringList; +use crate::errors::*; +use crate::raster::processing::dem::GdalDemProcessor; +use crate::raster::processing::{merge, GdalProcessor}; +use crate::Dataset; + +/// Topographic Position Index (TPI) processor for DEM datasets +/// +/// This processor outputs a single-band raster with values computed from the elevation. +/// A Topographic Position Index is defined as the difference between a central pixel and the +/// mean of its surrounding cells (see Wilson et al 2007, Marine Geodesy 30:3-35). +// +/// The value `-9999` is used as the output no-data value. +/// +/// Note: Results are nonsensical if the underlying [`Dataset`] does not contain digital elevation data. +/// +/// See: [`gdaldem tpi`](https://gdal.org/programs/gdaldem.html#tpi) for details. +#[derive(Debug, Clone)] +pub struct DemTpiProcessor { + input_band: Option, + destination: Option, + compute_edges: Option, + extra_options: Option, +} + +impl DemTpiProcessor { + /// Create a aspect-from-DEM processor + pub fn new() -> Self { + Self { + input_band: None, + destination: None, + compute_edges: None, + extra_options: None, + } + } + + /// Specify which band in the input [`Dataset`] to read from. + /// + /// Defaults to the first band. + pub fn with_input_band(self, band: NonZeroUsize) -> Self { + Self { + input_band: Some(band), + ..self + } + } + + /// Specify a file destination for the generated Dataset. + /// + /// Note: If none is specified, an anonymous in-memory file is used instead. + pub fn with_destination>(self, path: P) -> Self { + Self { + destination: Some(path.as_ref().to_path_buf()), + ..self + } + } + + /// Compute slope values at image edges. + /// + /// This causes interpolation of values at image edges or if a no-data value is found + /// in the 3x3 processing window. + pub fn with_compute_edges(self, state: bool) -> Self { + Self { + compute_edges: Some(state), + ..self + } + } + + /// Additional options. Applied last. + pub fn with_options(self, extra_options: CslStringList) -> Self { + Self { + extra_options: Some(extra_options), + ..self + } + } +} + +impl Default for DemTpiProcessor { + fn default() -> Self { + Self::new() + } +} + +impl GdalProcessor for DemTpiProcessor { + fn to_options(&self) -> CslStringList { + let mut opts = CslStringList::new(); + + if self.compute_edges == Some(true) { + opts.add_string("-compute_edges").unwrap(); + } + + if let Some(band) = self.input_band { + opts.add_string("-b").unwrap(); + opts.add_string(&band.to_string()).unwrap(); + } + + if let Some(extra) = &self.extra_options { + merge(&mut opts, extra); + } + + opts + } + + fn eval(&self, ds: &Dataset) -> Result { + GdalDemProcessor::dem_eval(self, ds) + } +} + +impl GdalDemProcessor for DemTpiProcessor { + fn mode(&self) -> &'static str { + "tpi" + } + + fn destination(&self) -> Option { + self.destination.clone() + } +} + +#[cfg(test)] +mod tests { + use crate::assert_near; + use crate::cpl::CslStringList; + use crate::raster::StatisticsAll; + use crate::test_utils::fixture; + + use super::*; + + #[test] + fn options() -> Result<()> { + let dsp = DemTpiProcessor::new() + .with_input_band(2.try_into().unwrap()) + .with_compute_edges(true) + .with_options("-of GTiff".parse()?); + + let expected: CslStringList = "-compute_edges -b 2 -of GTiff".parse()?; + assert_eq!(expected, dsp.to_options()); + + Ok(()) + } + + #[test] + fn tpi() -> Result<()> { + let proc = DemTpiProcessor::new(); + + let ds = Dataset::open(fixture("dem-hills.tiff"))?; + + let slope = proc.eval(&ds)?; + + let stats = slope.rasterband(1)?.get_statistics(true, false)?.unwrap(); + + // These numbers were generated by extracting the output from: + // gdaldem tpi fixtures/dem-hills.tiff target/dest.tiff + // gdalinfo -stats target/dest.tiff + let expected = StatisticsAll { + min: -4.7376708984375, + max: 4.7724151611328, + mean: 0.00012131847966826, + std_dev: 0.48943078832474 + }; + + assert_near!(StatisticsAll, stats, expected, epsilon = 1e-10); + Ok(()) + } +}