Merge pull request #456 from metasim/feature/dem

GDAL DEM Processing 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"

.github/workflows/ci.yml

@@ -69,7 +69,7 @@ jobs:
       - name: Build (--all-features)
         run: cargo build --all-features
       - name: Run tests (--all-features)
-        run: cargo test --all-features
+        run: cargo test --all-features -- --nocapture
 
   ubuntu_lts:
     name: "ci ubuntu-lts"
@@ -112,4 +112,4 @@ jobs:
       - name: Build (--all-features)
         run: cargo build --all-features
       - name: Run tests (--all-features)
-        run: cargo test --all-features
+        run: cargo test --all-features -- --nocapture

.gitignore

@@ -4,7 +4,7 @@
 /gdal-sys/target
 /.idea
 /3rd
-/fixtures/tinymarble.tif.aux.xml
+/fixtures/*.aux.xml
 
 # gtags
 GPATH

CHANGES.md

@@ -2,30 +2,35 @@
 
 ## Unreleased
 
+- Added support for digital elevation model raster processing: `aspect`, `color_relief`, `hillshade`, `roughness`, `slope`, `terrain_ruggedness_index`, `topographic_position_index`.
+
+   - <https://github.com/georust/gdal/pull/456> 
+
 - Added pre-built bindings for GDAL 3.8
 
-  - <https://github.com/georust/gdal/pull/466>
+   - <https://github.com/georust/gdal/pull/466>
+
 
 - Added `{Display|FromStr} for ResampleAlg` and `ResampleAlg::iter`.
 
-    - <https://github.com/georust/gdal/pull/462>
+   - <https://github.com/georust/gdal/pull/462>
 
- - **Breaking**: Replaced `TryFrom<&[(&str, &str); N]> for CslStringList` with `impl FromIterator<CslStringListEntry> for CslStringList`, `impl FromIterator<String> for CslStringList` and `impl<'a> FromIterator<&'a str> for CslStringList`
- - Added `Extend<CslStringListEntry> for CslStringList`, and `CslStringList::merge`
+- **Breaking**: Replaced `TryFrom<&[(&str, &str); N]> for CslStringList` with `impl FromIterator<CslStringListEntry> for CslStringList`, `impl FromIterator<String> for CslStringList` and `impl<'a> FromIterator<&'a str> for CslStringList`
+- Added `Extend<CslStringListEntry> for CslStringList`, and `CslStringList::merge`
 
-    - <https://github.com/georust/gdal/pull/457>
+  - <https://github.com/georust/gdal/pull/457>
 
 - **Breaking**: `SpatialRef::set_axis_mapping_strategy` now takes `&mut self`
 
-  - <https://github.com/georust/gdal/pull/461>
+   - <https://github.com/georust/gdal/pull/461>
 
 - **Breaking**: `Dataset::raster_count` now returns an `usize` and `Dataset::rasterband` now takes `usize` instead of `isize`
 
-  - <https://github.com/georust/gdal/pull/434>
+   - <https://github.com/georust/gdal/pull/434>
 
 - **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`.
+- Added `CslStringList::get_field`, `CslStringList::find_string`, `CslStringList::partial_find_string`, `CslStringList::find_string_case_sensitive`, `CslStringList::into_ptr`, `CslStringList::add_name_value`.
 
    - <https://github.com/georust/gdal/pull/455>
 

fixtures/color-relief.clr

@@ -0,0 +1,6 @@
+100%   white
+75%   235:220:175
+50%   190 185 135
+20%    240 250 150
+0%      50  180  50
+nv     0   0   0   0

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,124 @@
+use std::num::NonZeroUsize;
+
+use super::options::{common_dem_options, CommonOptions};
+use crate::cpl::CslStringList;
+use crate::errors;
+use crate::raster::processing::dem::DemSlopeAlg;
+
+/// Configuration options for [`aspect()`][super::aspect()].
+#[derive(Debug, Clone, Default)]
+pub struct AspectOptions {
+    common_options: CommonOptions,
+    algorithm: Option<DemSlopeAlg>,
+    zero_for_flat: Option<bool>,
+    trigonometric: Option<bool>,
+}
+
+impl AspectOptions {
+    /// Create a DEM-aspect options set.
+    pub fn new() -> Self {
+        Default::default()
+    }
+
+    common_dem_options!();
+
+    /// 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
+    }
+
+    /// Render relevant common options into [`CslStringList`] values, as compatible with
+    /// [`gdal_sys::GDALDEMProcessing`].
+    pub fn to_options_list(&self) -> errors::Result<CslStringList> {
+        let mut opts = CslStringList::default();
+
+        self.store_common_options_to(&mut opts)?;
+
+        if let Some(alg) = self.algorithm {
+            opts.add_string("-alg")?;
+            opts.add_string(alg.to_gdal_option())?;
+        }
+
+        if self.zero_for_flat == Some(true) {
+            opts.add_string("-zero_for_flat")?;
+        }
+
+        if self.trigonometric == Some(true) {
+            opts.add_string("-trigonometric")?;
+        }
+
+        Ok(opts)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::assert_near;
+    use crate::cpl::CslStringList;
+    use crate::errors::Result;
+    use crate::raster::processing::dem::aspect;
+    use crate::raster::StatisticsAll;
+    use crate::test_utils::{fixture, target};
+    use crate::Dataset;
+
+    use super::*;
+
+    #[test]
+    fn test_options() -> Result<()> {
+        let mut proc = AspectOptions::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_output_format("GTiff")
+            .with_additional_options("CPL_DEBUG=ON".parse()?);
+
+        let expected: CslStringList =
+            "-compute_edges -b 2 -of GTiff CPL_DEBUG=ON -alg ZevenbergenThorne -zero_for_flat -trigonometric"
+                .parse()?;
+        assert_eq!(expected.to_string(), proc.to_options_list()?.to_string());
+
+        Ok(())
+    }
+
+    #[test]
+    fn test_aspect() -> Result<()> {
+        let mut opts = AspectOptions::new();
+        opts.with_algorithm(DemSlopeAlg::Horn)
+            .with_zero_for_flat(true);
+
+        let ds = Dataset::open(fixture("dem-hills.tiff"))?;
+
+        let aspect = aspect(&ds, target("dem-hills-aspect.tiff"), &opts)?;
+
+        let stats = aspect.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(())
+    }
+}