Initial approach to binding DEM routines.

.cargo/config.toml

@@ -1,5 +1,9 @@
 [alias]
-# Run doctests, displaying compiler output
-dto = "test --doc -- --show-output"
+# Run doctests, displaying compiler output.
+# Due to this issue:
+#    https://github.com/rust-lang/cargo/pull/9705#issuecomment-1226149265
+# the following is required for full output during documentation development debugging:
+#    RUSTDOCFLAGS="-Z unstable-options --nocapture" cargo +nightly test --doc
+dto = "test --doc -- --show-output --nocapture"
 # Run clippy, raising warnings to errors
 nowarn = "clippy --all-targets -- -D warnings"

CHANGES.md

@@ -2,6 +2,10 @@
 
 ## Unreleased
 
+ - Added support for digital elevation model raster processing.
+
+   - <https://github.com/georust/gdal/pull/456> 
+
  - **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<String>` instead of `Result<Option<String>>`, 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`, `CslStringList::add_name_value`.

src/cpl.rs

@@ -29,7 +29,7 @@ use crate::utils::_string;
 ///
 /// # Example
 ///
-/// ```rust
+/// ```rust, no_run
 /// use gdal::cpl::CslStringList;
 ///
 /// let mut sl1 = CslStringList::new();

src/raster/mod.rs

@@ -76,6 +76,7 @@
 
 #[cfg(all(major_ge_3, minor_ge_1))]
 mod mdarray;
+pub mod processing;
 mod rasterband;
 mod rasterize;
 mod types;

src/raster/processing/dem/aspect.rs

@@ -0,0 +1,140 @@
+use crate::cpl::CslStringList;
+use crate::raster::processing::dem::{DemAlg, DemCommonOptions, DemSlopeAlg, GdalDemProcessor};
+
+/// Slope-aspect angle processor for DEM datasets.
+///
+/// This processor outputs a 32-bit float raster with values between 0° and 360°
+/// representing the azimuth that slopes are facing. The definition of the azimuth is such that:
+///
+/// * 0° means that the slope is facing the North,
+/// * 90° it's facing the East,
+/// * 180° it's facing the South and;
+/// * 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.
+///
+/// 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, Default)]
+pub struct DemAspectProcessor {
+    common_options: DemCommonOptions,
+    algorithm: Option<DemSlopeAlg>,
+    zero_for_flat: Option<bool>,
+    trigonometric: Option<bool>,
+}
+
+impl DemAspectProcessor {
+    /// Create a aspect-from-DEM processor
+    pub fn new() -> Self {
+        Default::default()
+    }
+
+    /// Specify the slope computation algorithm.
+    pub fn with_algorithm(&mut self, algorithm: DemSlopeAlg) -> &mut Self {
+        self.algorithm = Some(algorithm);
+        self
+    }
+
+    /// Return `0` for flat areas with `slope=0`, instead of `-9999`.
+    ///
+    /// See: [`zero_for_flat`](https://gdal.org/programs/gdaldem.html#cmdoption-zero_for_flat)
+    pub fn with_zero_for_flat(&mut self, state: bool) -> &mut Self {
+        self.zero_for_flat = Some(state);
+        self
+    }
+
+    /// Return trigonometric angle instead of azimuth. Thus 0° means East, 90° North, 180° West, 270° South.
+    pub fn with_trigonometric_angles(&mut self, state: bool) -> &mut Self {
+        self.trigonometric = Some(state);
+        self
+    }
+}
+
+impl GdalDemProcessor for DemAspectProcessor {
+    fn common_options(&self) -> &DemCommonOptions {
+        &self.common_options
+    }
+
+    fn common_options_mut(&mut self) -> &mut DemCommonOptions {
+        &mut self.common_options
+    }
+
+    fn to_options(&self) -> CslStringList {
+        let mut opts = self.common_options.to_options();
+
+        if let Some(alg) = self.algorithm {
+            opts.add_string("-alg").unwrap();
+            opts.add_string(&alg.to_string()).unwrap();
+        }
+
+        if self.zero_for_flat == Some(true) {
+            opts.add_string("-zero_for_flat").unwrap();
+        }
+
+        if self.trigonometric == Some(true) {
+            opts.add_string("-trigonometric").unwrap();
+        }
+
+        opts
+    }
+
+    fn dem_algorithm(&self) -> DemAlg {
+        DemAlg::Aspect
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::assert_near;
+    use crate::cpl::CslStringList;
+    use crate::errors::Result;
+    use crate::raster::StatisticsAll;
+    use crate::test_utils::fixture;
+    use crate::Dataset;
+
+    #[test]
+    fn options() -> Result<()> {
+        let mut proc = DemAspectProcessor::new();
+        proc.with_input_band(2.try_into().unwrap())
+            .with_algorithm(DemSlopeAlg::ZevenbergenThorne)
+            .with_compute_edges(true)
+            .with_zero_for_flat(true)
+            .with_trigonometric_angles(true)
+            .with_options("-of GTiff".parse()?);
+
+        let expected: CslStringList =
+            "-compute_edges -b 2 -of GTiff -alg ZevenbergenThorne -zero_for_flat -trigonometric"
+                .parse()?;
+        assert_eq!(expected.to_string(), proc.to_options().to_string());
+
+        Ok(())
+    }
+
+    #[test]
+    fn aspect() -> Result<()> {
+        let mut proc = DemAspectProcessor::new();
+        proc.with_algorithm(DemSlopeAlg::Horn)
+            .with_zero_for_flat(true);
+
+        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 aspect -alg Horn -zero_for_flat fixtures/dem-hills.tiff target/dest.tiff
+        //    gdalinfo -stats target/dest.tiff
+        let expected = StatisticsAll {
+            min: 0.0,
+            max: 359.9951171875,
+            mean: 165.72752499998,
+            std_dev: 98.590199951445,
+        };
+        assert_near!(StatisticsAll, stats, expected, epsilon = 1e-10);
+        Ok(())
+    }
+}