diff --git a/.cargo/config.toml b/.cargo/config.toml index e2d4b0a3..78181039 100644 --- a/.cargo/config.toml +++ b/.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" diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 6e492e6b..2370c9b6 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -68,7 +68,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" @@ -111,4 +111,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 diff --git a/.gitignore b/.gitignore index 446b55dd..34d67d08 100644 --- a/.gitignore +++ b/.gitignore @@ -4,7 +4,7 @@ /gdal-sys/target /.idea /3rd -/fixtures/tinymarble.tif.aux.xml +/fixtures/*.aux.xml # gtags GPATH diff --git a/CHANGES.md b/CHANGES.md index 1e056730..9224a66c 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -2,6 +2,11 @@ ## Unreleased + +- Added support for digital elevation model raster processing: `aspect`, `color_relief`, `hillshade`, `roughness`, `slope`, `terrain_ruggedness_index`, `topographic_position_index`. + + - + - Added `{Display|FromStr} for ResampleAlg` and `ResampleAlg::iter`. - @@ -21,7 +26,7 @@ - **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`, `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`. - diff --git a/fixtures/color-relief.clr b/fixtures/color-relief.clr new file mode 100644 index 00000000..c4115699 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/fixtures/dem-hills.tiff b/fixtures/dem-hills.tiff new file mode 100644 index 00000000..96497673 Binary files /dev/null and b/fixtures/dem-hills.tiff differ diff --git a/src/raster/mod.rs b/src/raster/mod.rs index 1e266b90..712902f1 100644 --- a/src/raster/mod.rs +++ b/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; diff --git a/src/raster/processing/dem/aspect.rs b/src/raster/processing/dem/aspect.rs new file mode 100644 index 00000000..a787cc57 --- /dev/null +++ b/src/raster/processing/dem/aspect.rs @@ -0,0 +1,125 @@ +use super::options::common_dem_options; +use crate::cpl::CslStringList; +use crate::raster::processing::dem::DemSlopeAlg; +use std::num::NonZeroUsize; + +/// Configuration options for [`aspect()`][super::aspect()]. +#[derive(Debug, Clone, Default)] +pub struct AspectOptions { + input_band: Option, + compute_edges: bool, + output_format: Option, + additional_options: CslStringList, + algorithm: Option, + zero_for_flat: Option, + trigonometric: Option, +} + +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) -> CslStringList { + let mut opts = CslStringList::default(); + + self.store_common_options_to(&mut opts); + + 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 + } +} + +#[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(()) + } +} diff --git a/src/raster/processing/dem/color_relief.rs b/src/raster/processing/dem/color_relief.rs new file mode 100644 index 00000000..230a8c9c --- /dev/null +++ b/src/raster/processing/dem/color_relief.rs @@ -0,0 +1,202 @@ +use std::num::NonZeroUsize; +use std::path::{Path, PathBuf}; + +use crate::cpl::CslStringList; +use crate::raster::processing::dem::options::common_dem_options; + +/// Configuration options for [`color_relief()`][super::color_relief()]. +#[derive(Debug, Clone)] +pub struct ColorReliefOptions { + input_band: Option, + compute_edges: bool, + output_format: Option, + additional_options: CslStringList, + color_config: PathBuf, + alpha: Option, + color_matching_mode: ColorMatchingMode, +} + +impl ColorReliefOptions { + /// Create a DEM-color-relief options set. + /// + /// `color_config` is a required text file mapping elevations to color specifiers. + /// Generally this means 4 columns per line: the elevation value, and the corresponding + /// _Red_, _Green_, and _Blue_ component (between 0 and 255). + /// The elevation value can be any floating point value, or the `nv` keyword for the no-data value. + /// + /// The elevation can also be expressed as a percentage: + /// 0% being the minimum value found in the raster, 100% the maximum value. + /// + /// An extra column can be optionally added for the alpha component. + /// If it is not specified, full opacity (255) is assumed. + /// + /// Various field separators are accepted: comma, tabulation, spaces, ':'. + /// + /// Common colors used by GRASS can also be specified by using their name, + /// instead of the RGB triplet. The supported list is: + /// _white_, _black_, _red_, _green_, _blue_, _yellow_, _magenta_, _cyan_, _aqua_, _grey/gray_, + /// _orange_, _brown_, _purple/violet_ and _indigo_. + /// + /// Note: The syntax of the color configuration file is derived from the one supported by + /// GRASS `r.colors` utility. ESRI HDR color table files (`.clr`) also match that syntax. + /// The alpha component and the support of tab and comma as separators are GDAL specific extensions. + /// + /// # Example + /// Here's an example `.clr` file showing a number of the features described above. + /// + /// ```text + /// 2600 white + /// 2000 235 220 175 + /// 50% 190 185 135 + /// 1000 240 250 150 + /// 100 50 180 50 128 + /// nv 0 0 0 0 + /// ``` + /// See: [gdaldem color-relief](https://gdal.org/programs/gdaldem.html#color-relief) for + /// details. + pub fn new>(color_config: P) -> Self { + Self { + input_band: None, + compute_edges: false, + output_format: None, + additional_options: Default::default(), + color_config: color_config.as_ref().to_path_buf(), + alpha: None, + color_matching_mode: Default::default(), + } + } + + common_dem_options!(); + + /// Add an alpha channel to the output raster + pub fn with_alpha(&mut self, state: bool) -> &mut Self { + self.alpha = Some(state); + self + } + + /// Get path to the color relief configuration file. + pub fn color_config(&self) -> PathBuf { + self.color_config.to_owned() + } + + /// Specify the color matching mode. + /// + /// See [`ColorMatchingMode`] for details. + pub fn with_color_matching_mode(&mut self, mode: ColorMatchingMode) -> &mut Self { + self.color_matching_mode = mode; + self + } + + /// Get the color matching mode to be used. + pub fn color_matching_mode(&self) -> ColorMatchingMode { + self.color_matching_mode + } + + /// Render relevant common options into [`CslStringList`] values, as compatible with + /// [`gdal_sys::GDALDEMProcessing`]. + pub fn to_options_list(&self) -> CslStringList { + let mut opts = CslStringList::default(); + + self.store_common_options_to(&mut opts); + + if self.alpha == Some(true) { + opts.add_string("-alpha").unwrap(); + } + + match self.color_matching_mode { + ColorMatchingMode::ExactColorEntry => opts.add_string("-exact_color_entry").unwrap(), + ColorMatchingMode::NearestColorEntry => { + opts.add_string("-nearest_color_entry").unwrap() + } + _ => {} + } + + opts + } +} + +/// Color relief color matching mode +#[derive(Debug, Clone, Copy, Default)] +pub enum ColorMatchingMode { + /// Colors between the given elevation values are blended smoothly. + /// + /// This is the default. + #[default] + Blended, + /// Use strict matching when searching in the color configuration file. + /// + /// If no matching color entry is found, the "0,0,0,0" RGBA quadruplet will be used. + ExactColorEntry, + /// Use the RGBA quadruplet corresponding to the closest entry in the color configuration file. + NearestColorEntry, +} + +#[cfg(test)] +mod tests { + use crate::assert_near; + use crate::cpl::CslStringList; + use crate::errors::Result; + use crate::raster::processing::dem::color_relief; + use crate::raster::StatisticsAll; + use crate::test_utils::{fixture, target}; + use crate::Dataset; + + use super::*; + + #[test] + fn test_options() -> Result<()> { + let mut proc = ColorReliefOptions::new("/dev/null"); + proc.with_input_band(2.try_into().unwrap()) + .with_compute_edges(true) + .with_alpha(true) + .with_color_matching_mode(ColorMatchingMode::NearestColorEntry) + .with_output_format("GTiff") + .with_additional_options("CPL_DEBUG=ON".parse()?); + + let expected: CslStringList = + "-compute_edges -b 2 -of GTiff CPL_DEBUG=ON -alpha -nearest_color_entry".parse()?; + assert_eq!(expected.to_string(), proc.to_options_list().to_string()); + + Ok(()) + } + + #[test] + fn test_color_relief() -> Result<()> { + let ds = Dataset::open(fixture("dem-hills.tiff"))?; + + let mut opts = ColorReliefOptions::new(fixture("color-relief.clr")); + opts.with_compute_edges(true); + + let cr = color_relief(&ds, target("dem-hills-relief.tiff"), &opts)?; + + // These numbers were generated by extracting the output from: + // gdaldem color-relief -compute_edges -alpha fixtures/dem-hills.tiff fixtures/color-relief.clr target/dest.tiff + // gdalinfo -stats target/dest.tiff + let expected = [ + StatisticsAll { + min: 0.0, + max: 255.0, + mean: 204.15542606827012, + std_dev: 57.41999604472401, + }, + StatisticsAll { + min: 0.0, + max: 255.0, + mean: 221.0581177507783, + std_dev: 33.229287115978394, + }, + StatisticsAll { + min: 0.0, + max: 255.0, + mean: 164.13047910295617, + std_dev: 60.78580825073262, + }, + ]; + for (i, e) in expected.iter().enumerate() { + let stats = cr.rasterband(i + 1)?.get_statistics(true, false)?.unwrap(); + assert_near!(StatisticsAll, stats, e, epsilon = 1e-8); + } + + Ok(()) + } +} diff --git a/src/raster/processing/dem/hillshade.rs b/src/raster/processing/dem/hillshade.rs new file mode 100644 index 00000000..44e771d9 --- /dev/null +++ b/src/raster/processing/dem/hillshade.rs @@ -0,0 +1,244 @@ +use crate::cpl::CslStringList; +use crate::raster::processing::dem::options::common_dem_options; +use crate::raster::processing::dem::DemSlopeAlg; +use std::fmt::{Display, Formatter}; +use std::num::NonZeroUsize; + +/// Configuration options for [`hillshade()`][super::hillshade()]. +#[derive(Debug, Clone, Default)] +pub struct HillshadeOptions { + input_band: Option, + compute_edges: bool, + output_format: Option, + additional_options: CslStringList, + algorithm: Option, + altitude: Option, + azimuth: Option, + scale: Option, + shading: Option, + z_factor: Option, +} + +impl HillshadeOptions { + /// Create a DEM-hillshade 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 + } + + /// Fetch the specified slope computation algorithm + pub fn algorithm(&self) -> Option { + self.algorithm + } + + /// Specify the altitude of the light, in degrees. + /// + /// `90` if the light comes from above the DEM, `0` if it is raking light. + pub fn with_altitude(&mut self, altitude: f64) -> &mut Self { + self.altitude = Some(altitude); + self + } + + /// Fetch the specified light altitude, in degrees. + pub fn altitude(&self) -> Option { + self.altitude + } + + /// Specify the azimuth of the light, in degrees: + /// + /// * `0` if it comes from the top of the raster, + /// * `90` from the east, + /// * etc. + /// + /// The default value, `315`, and should rarely be changed as it is the value generally + /// used to generate shaded maps. + pub fn with_azimuth(&mut self, azimuth: f64) -> &mut Self { + self.azimuth = Some(azimuth); + self + } + + /// Apply a elevation scaling factor. + /// + /// Routine assumes x, y and z units are identical. + /// If x (east-west) and y (north-south) units are identical, but z (elevation) units are different, + /// this scale option can be used to set the ratio of vertical units to horizontal. + /// + /// For LatLong projections near the equator, where units of latitude and units of longitude are + /// similar, elevation (z) units can be converted with the following values: + /// + /// * Elevation in feet: `370400` + /// * Elevation in meters: `111120` + /// + /// For locations not near the equator, it would be best to reproject your raster first. + pub fn with_scale(&mut self, scale: f64) -> &mut Self { + self.scale = Some(scale); + self + } + + /// Fetch the specified scaling factor. + /// + /// Returns `None` if one has not been previously set vai [`Self::with_scale`]. + pub fn scale(&self) -> Option { + self.scale + } + + /// Specify the shading mode to render with. + /// + /// See [`ShadingMode`] for mode descriptions. + pub fn with_shading_mode(&mut self, mode: ShadingMode) -> &mut Self { + self.shading = Some(mode); + self + } + + /// Fetch the specified shading mode. + pub fn shading_mode(&self) -> Option { + self.shading + } + + /// Vertical exaggeration used to pre-multiply the elevations + pub fn with_z_factor(&mut self, z_factor: f64) -> &mut Self { + self.z_factor = Some(z_factor); + self + } + + /// Fetch the applied z-factor value. + pub fn z_factor(&self) -> Option { + self.z_factor + } + + /// Render relevant common options into [`CslStringList`] values, as compatible with + /// [`gdal_sys::GDALDEMProcessing`]. + pub fn to_options_list(&self) -> CslStringList { + let mut opts = CslStringList::default(); + + self.store_common_options_to(&mut opts); + + if let Some(alg) = self.algorithm { + opts.add_string("-alg").unwrap(); + opts.add_string(&alg.to_string()).unwrap(); + } + + if let Some(scale) = self.scale { + opts.add_string("-s").unwrap(); + opts.add_string(&scale.to_string()).unwrap(); + } + + if let Some(mode) = self.shading { + opts.add_string(&format!("-{mode}")).unwrap(); + } + + if let Some(factor) = self.z_factor { + opts.add_string("-z").unwrap(); + opts.add_string(&factor.to_string()).unwrap(); + } + + if let Some(altitude) = self.altitude { + opts.add_string("-alt").unwrap(); + opts.add_string(&altitude.to_string()).unwrap(); + } + + if let Some(azimuth) = self.azimuth { + opts.add_string("-az").unwrap(); + opts.add_string(&azimuth.to_string()).unwrap(); + } + + opts + } +} + +/// Hillshade shading mode. +#[derive(Debug, Clone, Copy)] +pub enum ShadingMode { + /// Combination of slope and oblique shading. + Combined, + /// Multi-directional shading, + /// + /// A combination of hillshading illuminated from 225 deg, 270 deg, 315 deg, and 360 deg azimuth. + /// + /// See: . + Multidirectional, + /// Shading which tries to minimize effects on other map features beneath. + /// + /// Can't be used `altitude` specification + /// + /// See: . + Igor, +} + +impl Display for ShadingMode { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + let s = format!("{self:?}"); + f.write_str(&s.to_lowercase()) + } +} + +#[cfg(test)] +mod tests { + use crate::assert_near; + use crate::cpl::CslStringList; + use crate::errors::Result; + use crate::raster::processing::dem::hillshade; + use crate::raster::StatisticsAll; + use crate::test_utils::{fixture, target}; + use crate::Dataset; + + use super::*; + + #[test] + fn test_options() -> Result<()> { + let mut proc = HillshadeOptions::new(); + proc.with_input_band(2.try_into().unwrap()) + .with_algorithm(DemSlopeAlg::ZevenbergenThorne) + .with_scale(98473.0) + .with_shading_mode(ShadingMode::Igor) + .with_compute_edges(true) + .with_azimuth(330.0) + .with_altitude(45.0) + .with_z_factor(2.0) + .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 -s 98473 -igor -z 2 -alt 45 -az 330" + .parse()?; + assert_eq!(expected.to_string(), proc.to_options_list().to_string()); + + Ok(()) + } + + #[test] + fn test_hillshade() -> Result<()> { + let ds = Dataset::open(fixture("dem-hills.tiff"))?; + let scale_factor = 98473.2947; + + let mut opts = HillshadeOptions::new(); + opts.with_algorithm(DemSlopeAlg::ZevenbergenThorne) + .with_shading_mode(ShadingMode::Igor) + .with_z_factor(2.0) + .with_scale(scale_factor); + + let shade = hillshade(&ds, target("dem-hills-shade.tiff"), &opts)?; + + let stats = shade.rasterband(1)?.get_statistics(true, false)?.unwrap(); + + // These numbers were generated by extracting the output from: + // gdaldem hillshade -alg ZevenbergenThorne -s 98473.2947 -igor -z 2 fixtures/dem-hills.tiff target/dest.tiff + // gdalinfo -stats target/dest.tiff + let expected = StatisticsAll { + min: 128.0, + max: 255.0, + mean: 244.15731356401, + std_dev: 16.76881437538, + }; + + assert_near!(StatisticsAll, stats, expected, epsilon = 1e-8); + Ok(()) + } +} diff --git a/src/raster/processing/dem/mod.rs b/src/raster/processing/dem/mod.rs new file mode 100644 index 00000000..12bdb199 --- /dev/null +++ b/src/raster/processing/dem/mod.rs @@ -0,0 +1,474 @@ +//! Digital Elevation Model (DEM) processing routines. +//! +//! This module provides bindings to the algorithms in the +//! [`gdaldem` tool](https://gdal.org/programs/gdaldem.html#gdaldem). +//! +//! The routines assume an open dataset containing customary digital elevation model data. +//! This includes assumptions that `x` (east-west), `y` (north-south), and `z` (elevation) units are identical. +//! If `x` and `y` units are identical, but `z` (elevation) units are different, +//! `hillshade` and `slope` support a scale setting to set the ratio of vertical units to horizontal. +//! See [`SlopeOptions::with_scale`] for details. +//! +//! # Examples +//! +//! Examples may be found associated with the following functions: +//! +//! * [`aspect()`] +//! * [`color_relief()`] +//! * [`hillshade()`] +//! * [`roughness()`] +//! * [`slope()`] +//! * [`terrain_ruggedness_index()`] +//! * [`topographic_position_index()`] +//! + +use std::ffi::{CStr, CString}; +use std::path::{Path, PathBuf}; +use std::ptr; + +use libc::c_int; + +pub use aspect::*; +pub use color_relief::*; +use gdal_sys::{CPLErr, GDALDEMProcessing}; +pub use hillshade::*; +pub use options::{DemAlg, DemSlopeAlg}; +pub use roughness::*; +pub use slope::*; +pub use tpi::*; +pub use tri::*; + +use crate::cpl::CslStringList; +use crate::errors::Result; +use crate::utils::{_last_cpl_err, _path_to_c_string}; +use crate::Dataset; + +mod aspect; +mod color_relief; +mod hillshade; +mod options; +mod roughness; +mod slope; +mod tpi; +mod tri; + +/// Slope aspect-angle routine for DEM datasets. +/// +/// This method 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 [`AspectOptions::with_zero_for_flat`] for alternative. +/// +/// Note: Results are nonsensical if the underlying [`Dataset`] does not contain digital elevation data. +/// +/// # Example +/// +/// ```rust, no_run +/// use gdal::Dataset; +/// # fn main() -> gdal::errors::Result<()> { +/// use std::path::Path; +/// use gdal::raster::processing::dem::*; +/// let ds = Dataset::open("fixtures/dem-hills.tiff")?; +/// let mut opts = AspectOptions::new(); +/// opts +/// .with_algorithm(DemSlopeAlg::Horn) +/// .with_zero_for_flat(true); +/// let aspect_ds = aspect(&ds, Path::new("target/dem-hills-aspect.tiff"), &opts)?; +/// let stats = aspect_ds.rasterband(1)?.get_statistics(true, false)?.unwrap(); +/// println!("{stats:#?}"); +/// # Ok(()) +/// # } +/// ``` +/// The resulting output is: +/// +/// ```text +/// StatisticsAll { +/// min: 0.0, +/// max: 359.9951171875, +/// mean: 165.72752499997543, +/// std_dev: 98.5901999514453, +/// } +/// ``` +/// +/// See: [`gdaldem aspect`](https://gdal.org/programs/gdaldem.html#aspect) for details, +pub fn aspect>( + ds: &Dataset, + dest_file: P, + options: &AspectOptions, +) -> Result { + dem_eval( + ds, + dest_file.as_ref(), + DemAlg::Aspect, + &options.to_options_list(), + None, + ) +} + +/// Generate a color-relief rendering of DEM data. +/// +/// This routine outputs a 3-band (RGB) or 4-band (RGBA) raster with values computed from +/// the elevation and a text-based color configuration file. +/// +/// The color configuration file contains associations between various elevation values +/// and the corresponding desired color. See [`ColorReliefOptions::new`] for details. +/// +/// By default, the colors between the given elevation +/// values are blended smoothly and the result is a nice colorized DEM. +/// The [`ColorMatchingMode::ExactColorEntry`] or [`ColorMatchingMode::NearestColorEntry`] options +/// can be used to avoid that linear interpolation for values that don't match an index of +/// the color configuration file. See [`ColorMatchingMode`] for details. +/// +/// Note: Results are nonsensical if the underlying [`Dataset`] does not contain digital elevation data. +/// +/// # Example +/// +/// ```rust, no_run +/// use gdal::Dataset; +/// # fn main() -> gdal::errors::Result<()> { +/// use std::path::Path; +/// use gdal::raster::processing::dem::*; +/// let ds = Dataset::open("fixtures/dem-hills.tiff")?; +/// let mut opts = ColorReliefOptions::new("fixtures/color-relief.clr"); +/// opts.with_alpha(true); +/// let hs_ds = color_relief(&ds, Path::new("target/dem-hills-relief.tiff"), &opts)?; +/// // Note: Output will actually be a 4-band raster. +/// let stats = hs_ds.rasterband(1)?.get_statistics(true, false)?.unwrap(); +/// println!("{stats:#?}"); +/// # Ok(()) +/// # } +/// ``` +/// The resulting output is: +/// +/// ```text +/// StatisticsAll { +/// min: 50.0, +/// max: 255.0, +/// mean: 206.85964128690114, +/// std_dev: 52.73836661993173, +/// } +/// ``` +/// See: [`gdaldem color-relief`](https://gdal.org/programs/gdaldem.html#color-relief) for details, +/// +pub fn color_relief>( + ds: &Dataset, + dest_file: P, + options: &ColorReliefOptions, +) -> Result { + let colors = options.color_config(); + dem_eval( + ds, + dest_file.as_ref(), + DemAlg::ColorRelief, + &options.to_options_list(), + Some(colors), + ) +} + +/// Performs hill-shade rendering of DEM data. +/// +/// This routine outputs an 8-bit raster with a nice shaded relief effect. +/// It’s very useful for visualizing the terrain. +/// You can optionally specify the azimuth and altitude of the light source, +/// a vertical exaggeration factor and a scaling factor to account for +/// differences between vertical and horizontal units. +/// +/// The value `0` is used as the output no-data value. +/// +/// Note: Results are nonsensical if the underlying [`Dataset`] does not contain digital elevation data. +/// +/// # Example +/// +/// ```rust, no_run +/// use gdal::Dataset; +/// use gdal::raster::processing::dem::*; +/// # fn main() -> gdal::errors::Result<()> { +/// use std::path::Path; +/// let ds = Dataset::open("fixtures/dem-hills.tiff")?; +/// let mut opts = HillshadeOptions::new(); +/// opts +/// .with_algorithm(DemSlopeAlg::Horn) +/// .with_z_factor(4.0) +/// .with_scale(98473.0) +/// .with_shading_mode(ShadingMode::Combined); +/// let hs_ds = hillshade(&ds, Path::new("target/dem-hills-shade.tiff"), &opts)?; +/// let stats = hs_ds.rasterband(1)?.get_statistics(true, false)?.unwrap(); +/// println!("{stats:#?}"); +/// # Ok(()) +/// # } +/// ``` +/// The resulting output is: +/// +/// ```text +/// StatisticsAll { +/// min: 31.0, +/// max: 255.0, +/// mean: 234.71988886841396, +/// std_dev: 30.556285572761446, +/// } +/// ``` +/// See: [`gdaldem hillshade`](https://gdal.org/programs/gdaldem.html#hillshade) for details, +/// +pub fn hillshade>( + ds: &Dataset, + dest_file: P, + options: &HillshadeOptions, +) -> Result { + dem_eval( + ds, + dest_file.as_ref(), + DemAlg::Hillshade, + &options.to_options_list(), + None, + ) +} + +/// Roughness routine for DEM datasets. +/// +/// This processor outputs a single-band raster with values computed from the elevation. +/// Roughness is the largest inter-cell difference of a central pixel and its surrounding cell, +/// as defined in 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. +/// +/// # Example +/// +/// ```rust, no_run +/// # fn main() -> gdal::errors::Result<()> { +/// use gdal::Dataset; +/// use std::path::Path; +/// use gdal::raster::processing::dem::*; +/// let ds = Dataset::open("fixtures/dem-hills.tiff")?; +/// let roughness_ds = roughness(&ds, Path::new("target/dem-hills-roughness.tiff"), &RoughnessOptions::default())?; +/// let stats = roughness_ds.rasterband(1)?.get_statistics(true, false)?.unwrap(); +/// println!("{stats:#?}"); +/// # Ok(()) +/// # } +/// ``` +/// The resulting output is: +/// +/// ```text +/// StatisticsAll { +/// min: 0.0, +/// max: 14.361007690429688, +/// mean: 1.5128357817365072, +/// std_dev: 2.0120679959607686, +/// } +/// ``` +/// +/// See: [`gdaldem roughness`](https://gdal.org/programs/gdaldem.html#roughness) for details. +pub fn roughness>( + ds: &Dataset, + dest_file: P, + options: &RoughnessOptions, +) -> Result { + dem_eval( + ds, + dest_file.as_ref(), + DemAlg::Roughness, + &options.to_options_list(), + None, + ) +} + +/// Slope computation routine for DEM datasets. +/// +/// This method will take a DEM raster Dataset 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](SlopeOptions::with_percentage_results). +/// +/// In cases where the horizontal units differ from the vertical units, you can also supply +/// a [scaling factor](SlopeOptions::with_scale). +/// +/// 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. +/// +/// # Example +/// +/// ```rust, no_run +/// # fn main() -> gdal::errors::Result<()> { +/// use std::path::Path; +/// use gdal::Dataset; +/// use gdal::raster::processing::dem::*; +/// let ds = Dataset::open("fixtures/dem-hills.tiff")?; +/// let mut opts = SlopeOptions::new(); +/// opts +/// .with_algorithm(DemSlopeAlg::Horn) +/// .with_percentage_results(true) +/// .with_scale(98473.0); +/// let slope_ds = slope(&ds, Path::new("target/dem-hills-slope.tiff"), &opts)?; +/// let stats = slope_ds.rasterband(1)?.get_statistics(true, false)?.unwrap(); +/// println!("{stats:#?}"); +/// # Ok(()) +/// # } +/// ``` +/// The resulting output is: +/// +/// ```text +/// StatisticsAll { +/// min: 0.0, +/// max: 65.44061279296875, +/// mean: 6.171115265248098, +/// std_dev: 8.735612161193623, +/// } +/// ``` +/// +/// See: [`gdaldem slope`](https://gdal.org/programs/gdaldem.html#slope) for details, +pub fn slope>( + ds: &Dataset, + dest_file: P, + options: &SlopeOptions, +) -> Result { + dem_eval( + ds, + dest_file.as_ref(), + DemAlg::Slope, + &options.to_options_list(), + None, + ) +} + +/// Topographic Position Index (TPI) routine for DEM datasets +/// +/// This method 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. +/// +/// # Example +/// +/// ```rust, no_run +/// # fn main() -> gdal::errors::Result<()> { +/// use std::path::Path; +/// use gdal::Dataset; +/// use gdal::raster::processing::dem::*; +/// let ds = Dataset::open("fixtures/dem-hills.tiff")?; +/// let tpi_ds = topographic_position_index(&ds, Path::new("target/dem-hills-tpi.tiff"), &TpiOptions::default())?; +/// let stats = tpi_ds.rasterband(1)?.get_statistics(true, false)?.unwrap(); +/// println!("{stats:#?}"); +/// # Ok(()) +/// # } +/// ``` +/// The resulting output is: +/// +/// ```text +/// StatisticsAll { +/// min: -4.7376708984375, +/// max: 4.7724151611328125, +/// mean: 0.00012131847966825689, +/// std_dev: 0.48943078832474257, +/// } +/// ``` +/// +/// See: [`gdaldem tpi`](https://gdal.org/programs/gdaldem.html#tpi) for details. +pub fn topographic_position_index>( + ds: &Dataset, + dest_file: P, + options: &TpiOptions, +) -> Result { + dem_eval( + ds, + dest_file.as_ref(), + DemAlg::Tpi, + &options.to_options_list(), + None, + ) +} + +/// Terrain Ruggedness Index (TRI) routine for DEM datasets +/// +/// This method outputs a single-band raster with values computed from the elevation. +/// TRI stands for Terrain Ruggedness Index, which measures the difference between a +/// central pixel and its surrounding cells. +/// +/// The value `-9999` is used as the output no-data value. +/// +/// # Example +/// +/// ```rust, no_run +/// # fn main() -> gdal::errors::Result<()> { +/// use std::path::Path; +/// use gdal::Dataset; +/// use gdal::raster::processing::dem::*; +/// let ds = Dataset::open("fixtures/dem-hills.tiff")?; +/// let mut opts = TriOptions::new(); +/// opts.with_algorithm(DemTriAlg::Wilson); +/// let tri_ds = terrain_ruggedness_index(&ds, Path::new("target/dem-hills-tri.tiff"), &opts)?; +/// let stats = tri_ds.rasterband(1)?.get_statistics(true, false)?.unwrap(); +/// println!("{stats:#?}"); +/// # Ok(()) +/// # } +/// ``` +/// The resulting output is: +/// +/// ```text +/// StatisticsAll { +/// min: 0.0, +/// max: 4.983623504638672, +/// mean: 0.49063101456532326, +/// std_dev: 0.6719356336694824, +/// } +/// ``` +/// +/// See: [`gdaldem tri`](https://gdal.org/programs/gdaldem.html#tri) for details. +pub fn terrain_ruggedness_index>( + ds: &Dataset, + dest_file: P, + options: &TriOptions, +) -> Result { + dem_eval( + ds, + dest_file.as_ref(), + DemAlg::Tri, + &options.to_options_list(), + None, + ) +} + +/// Execute the processor on the given [`Dataset`]. +fn dem_eval( + src: &Dataset, + dst_file: &Path, + alg: DemAlg, + options: &CslStringList, + color_relief_config: Option, +) -> Result { + let popts = options::GdalDEMProcessingOptions::new(options)?; + let mode = CString::new(alg.to_string())?; + let dest = _path_to_c_string(dst_file)?; + let cfile = color_relief_config.and_then(|p| _path_to_c_string(&p).ok()); + let cfile_ptr = cfile.as_deref().map(CStr::as_ptr).unwrap_or(ptr::null()); + + let mut pb_usage_error: c_int = 0; + let out_ds = unsafe { + // Docs: https://github.com/OSGeo/gdal/blob/6a3584b2fea51f92022d24ad8036749ba1b98958/apps/gdaldem_lib.cpp#L3184 + GDALDEMProcessing( + dest.as_ptr(), + src.c_dataset(), + mode.as_ptr(), + cfile_ptr, + popts.as_ptr(), + &mut pb_usage_error as *mut c_int, + ) + }; + + if pb_usage_error != 0 { + Err(_last_cpl_err(CPLErr::CE_Failure)) + } else { + let out_ds = unsafe { Dataset::from_c_dataset(out_ds) }; + Ok(out_ds) + } +} diff --git a/src/raster/processing/dem/options.rs b/src/raster/processing/dem/options.rs new file mode 100644 index 00000000..76b33113 --- /dev/null +++ b/src/raster/processing/dem/options.rs @@ -0,0 +1,172 @@ +use std::fmt::{Display, Formatter}; +use std::marker::PhantomData; +use std::ptr; +use std::ptr::NonNull; + +use gdal_sys::{ + GDALDEMProcessingOptions, GDALDEMProcessingOptionsFree, GDALDEMProcessingOptionsNew, +}; + +use crate::cpl::CslStringList; +use crate::errors; +use crate::utils::_last_null_pointer_err; + +/// Payload for [`GDALDEMProcessing`]. Intended for internal use only. +pub struct GdalDEMProcessingOptions<'opts>( + NonNull, + PhantomData<&'opts CslStringList>, +); + +impl<'opts> GdalDEMProcessingOptions<'opts> { + pub fn new(opts: &'opts CslStringList) -> errors::Result { + let popts = unsafe { GDALDEMProcessingOptionsNew(opts.as_ptr(), ptr::null_mut()) }; + if popts.is_null() { + return Err(_last_null_pointer_err("GDALDEMProcessingOptionsNew")); + } + Ok(Self(unsafe { NonNull::new_unchecked(popts) }, PhantomData)) + } + + pub fn as_ptr(&self) -> *const GDALDEMProcessingOptions { + self.0.as_ptr() + } +} + +impl Drop for GdalDEMProcessingOptions<'_> { + fn drop(&mut self) { + unsafe { GDALDEMProcessingOptionsFree(self.0.as_ptr()) }; + } +} + +/// DEM processor mode, to stringify and pass to [`gdal_sys::GDALDEMProcessing`]. +#[derive(Debug, Clone, Copy)] +pub enum DemAlg { + Aspect, + ColorRelief, + Hillshade, + Roughness, + Slope, + Tpi, + Tri, +} + +impl Display for DemAlg { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + match self { + Self::ColorRelief => f.write_str("color-relief"), + _ => { + let s = format!("{self:?}").to_lowercase(); + f.write_str(&s) + } + } + } +} + +/// Slope and slope-related (aspect, hillshade) processing algorithms. +/// +/// 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 { + Horn, + ZevenbergenThorne, +} + +impl Display for DemSlopeAlg { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + f.write_fmt(format_args!("{self:?}")) + } +} + +macro_rules! common_dem_options { + () => { + /// Specify which band in the input [`Dataset`][crate::Dataset] to read from. + /// + /// Defaults to the first band. + pub fn with_input_band(&mut self, band: NonZeroUsize) -> &mut Self { + self.input_band = Some(band); + self + } + + /// Fetch the specified input band to read from. + pub fn input_band(&self) -> Option { + self.input_band + } + + /// Explicitly specify output raster format. + /// + /// This is equivalent to the `-of ` CLI flag accepted by many GDAL tools. + /// + /// The value of `format` must be the identifier of a driver supported by the runtime + /// environment's GDAL library (e.g. `COG`, `JPEG`, `VRT`, etc.). A list of these identifiers + /// is available from `gdalinfo --formats`: + /// + /// ```text + /// ❯ gdalinfo --formats + /// Supported Formats: + /// VRT -raster,multidimensional raster- (rw+v): Virtual Raster + /// DERIVED -raster- (ro): Derived datasets using VRT pixel functions + /// GTiff -raster- (rw+vs): GeoTIFF + /// COG -raster- (wv): Cloud optimized GeoTIFF generator + /// NITF -raster- (rw+vs): National Imagery Transmission Format + /// ... + /// ``` + /// + pub fn with_output_format(&mut self, format: &str) -> &mut Self { + self.output_format = Some(format.to_owned()); + self + } + + /// Fetch the specified output format driver identifier. + pub fn output_format(&self) -> Option { + self.output_format.clone() + } + + /// Compute values at image edges. + /// + /// If true, 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(&mut self, state: bool) -> &mut Self { + self.compute_edges = state; + self + } + + /// Fetch the compute edges mode. + pub fn compute_edges(&self) -> bool { + self.compute_edges + } + + /// Additional generic options to be included. + pub fn with_additional_options(&mut self, extra_options: CslStringList) -> &mut Self { + self.additional_options.extend(&extra_options); + self + } + + /// Fetch additional options. + pub fn additional_options(&self) -> &CslStringList { + &self.additional_options + } + + /// Private utility to convert common options into [`CslStringList`] options. + fn store_common_options_to(&self, opts: &mut CslStringList) { + if self.compute_edges { + 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(of) = &self.output_format { + opts.add_string("-of").unwrap(); + opts.add_string(of).unwrap(); + } + + if !self.additional_options.is_empty() { + opts.extend(&self.additional_options); + } + } + }; +} + +pub(crate) use common_dem_options; diff --git a/src/raster/processing/dem/roughness.rs b/src/raster/processing/dem/roughness.rs new file mode 100644 index 00000000..089c8418 --- /dev/null +++ b/src/raster/processing/dem/roughness.rs @@ -0,0 +1,83 @@ +use crate::cpl::CslStringList; +use crate::raster::processing::dem::options::common_dem_options; +use std::num::NonZeroUsize; + +/// Configuration options for [`roughness()`][super::roughness()]. + +#[derive(Debug, Clone, Default)] +pub struct RoughnessOptions { + input_band: Option, + compute_edges: bool, + output_format: Option, + additional_options: CslStringList, +} + +impl RoughnessOptions { + /// Create a DEM-roughness options set. + pub fn new() -> Self { + Default::default() + } + + common_dem_options!(); + + /// Render relevant options into [`CslStringList`] values, as compatible with + /// [`gdal_sys::GDALDEMProcessing`]. + pub fn to_options_list(&self) -> CslStringList { + let mut opts = CslStringList::new(); + self.store_common_options_to(&mut opts); + opts + } +} + +#[cfg(test)] +mod tests { + use crate::assert_near; + use crate::cpl::CslStringList; + use crate::errors::Result; + use crate::raster::processing::dem::roughness; + use crate::raster::processing::dem::roughness::RoughnessOptions; + use crate::raster::StatisticsAll; + use crate::test_utils::{fixture, target}; + use crate::Dataset; + + #[test] + fn test_options() -> Result<()> { + let mut proc = RoughnessOptions::new(); + proc.with_input_band(2.try_into().unwrap()) + .with_compute_edges(true) + .with_output_format("GTiff") + .with_additional_options("CPL_DEBUG=ON".parse()?); + + let expected: CslStringList = "-compute_edges -b 2 -of GTiff CPL_DEBUG=ON".parse()?; + assert_eq!(expected.to_string(), proc.to_options_list().to_string()); + + Ok(()) + } + + #[test] + fn test_roughness() -> Result<()> { + let opts = RoughnessOptions::new(); + + let ds = Dataset::open(fixture("dem-hills.tiff"))?; + + let roughness = roughness(&ds, target("dem-hills-roughness.tiff"), &opts)?; + + let stats = roughness + .rasterband(1)? + .get_statistics(true, false)? + .unwrap(); + + // These numbers were generated by extracting the output from: + // gdaldem roughness fixtures/dem-hills.tiff target/dest.tiff + // gdalinfo -stats target/dest.tiff + let expected = StatisticsAll { + min: 0.0, + max: 14.36100769043, + mean: 1.5128357817365, + std_dev: 2.0120679959608, + }; + + assert_near!(StatisticsAll, stats, expected, epsilon = 1e-10); + Ok(()) + } +} diff --git a/src/raster/processing/dem/slope.rs b/src/raster/processing/dem/slope.rs new file mode 100644 index 00000000..f1e8b8fd --- /dev/null +++ b/src/raster/processing/dem/slope.rs @@ -0,0 +1,159 @@ +use crate::cpl::CslStringList; +use crate::raster::processing::dem::options::common_dem_options; +use crate::raster::processing::dem::DemSlopeAlg; +use std::num::NonZeroUsize; + +/// Configuration options for [`slope()`][super::slope()]. +#[derive(Debug, Clone, Default)] +pub struct SlopeOptions { + input_band: Option, + compute_edges: bool, + output_format: Option, + additional_options: CslStringList, + algorithm: Option, + scale: Option, + percentage_results: Option, +} + +impl SlopeOptions { + /// Create a DEM-slope 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 + } + + /// Fetch the specified slope computation algorithm. + pub fn algorithm(&self) -> Option { + self.algorithm + } + + /// Apply a elevation scaling factor. + /// + /// Routine assumes x, y and z units are identical. + /// If x (east-west) and y (north-south) units are identical, but z (elevation) units are different, + /// this scale option can be used to set the ratio of vertical units to horizontal. + /// + /// For LatLong projections near the equator, where units of latitude and units of longitude are + /// similar, elevation (z) units can be converted with the following values: + /// + /// * Elevation in feet: `370400` + /// * Elevation in meters: `111120` + /// + /// For locations not near the equator, it would be best to reproject your raster first. + pub fn with_scale(&mut self, scale: f64) -> &mut Self { + self.scale = Some(scale); + self + } + + /// Fetch the specified scaling factor. + /// + /// Returns `None` if one has not been previously set via [`Self::with_scale`]. + pub fn scale(&self) -> Option { + self.scale + } + + /// If `state` is `true`, the slope will be expressed as percent slope. + /// + /// Otherwise, it is expressed as degrees + pub fn with_percentage_results(&mut self, state: bool) -> &mut Self { + self.percentage_results = Some(state); + self + } + + /// Render relevant common options into [`CslStringList`] values, as compatible with + /// [`gdal_sys::GDALDEMProcessing`]. + pub fn to_options_list(&self) -> CslStringList { + let mut opts = CslStringList::default(); + + self.store_common_options_to(&mut opts); + + if let Some(alg) = self.algorithm { + opts.add_string("-alg").unwrap(); + opts.add_string(&alg.to_string()).unwrap(); + } + + if let Some(scale) = self.scale { + opts.add_string("-s").unwrap(); + opts.add_string(&scale.to_string()).unwrap(); + } + + if self.percentage_results == Some(true) { + opts.add_string("-p").unwrap(); + } + + opts + } +} + +#[cfg(test)] +mod tests { + use crate::cpl::CslStringList; + use crate::errors::Result; + use crate::raster::processing::dem::slope; + use crate::raster::StatisticsAll; + use crate::test_utils::{fixture, target}; + use crate::Dataset; + use crate::{assert_near, GeoTransformEx}; + + use super::*; + + #[test] + fn test_options() -> Result<()> { + let mut proc = SlopeOptions::new(); + proc.with_input_band(2.try_into().unwrap()) + .with_algorithm(DemSlopeAlg::ZevenbergenThorne) + .with_scale(98473.0) + .with_compute_edges(true) + .with_percentage_results(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 -s 98473 -p" + .parse()?; + assert_eq!(expected.to_string(), proc.to_options_list().to_string()); + + Ok(()) + } + + #[test] + fn test_slope() -> Result<()> { + // For X & Y in degrees and Z in meters... + fn scaling_estimate(ds: &Dataset) -> f64 { + let (_, lat) = ds.geo_transform().unwrap().apply(0., 0.); + 0.5 * (111120.0 + lat.to_radians().cos() * 111120.0) + } + + let ds = Dataset::open(fixture("dem-hills.tiff"))?; + let scale_factor = scaling_estimate(&ds); + + let mut opts = SlopeOptions::new(); + opts.with_algorithm(DemSlopeAlg::Horn) + .with_percentage_results(true) + .with_scale(scale_factor); + + let slope = slope(&ds, target("dem-hills-slope.tiff"), &opts)?; + + let stats = slope.rasterband(1)?.get_statistics(true, false)?.unwrap(); + + // These numbers were generated by extracting the output from: + // gdaldem slope -alg Horn -s 98473.2947 -p fixtures/dem-hills.tiff target/dest.tiff + // gdalinfo -stats target/dest.tiff + let expected = StatisticsAll { + min: 0.0, + max: 65.440422058105, + mean: 6.1710967990449, + std_dev: 8.73558602352, + }; + + assert_near!(StatisticsAll, stats, expected, epsilon = 1e-8); + Ok(()) + } +} diff --git a/src/raster/processing/dem/tpi.rs b/src/raster/processing/dem/tpi.rs new file mode 100644 index 00000000..574a846e --- /dev/null +++ b/src/raster/processing/dem/tpi.rs @@ -0,0 +1,78 @@ +use crate::cpl::CslStringList; +use crate::raster::processing::dem::options::common_dem_options; +use std::num::NonZeroUsize; + +/// Configuration options for [`topographic_position_index()`][super::topographic_position_index()]. +#[derive(Debug, Clone, Default)] +pub struct TpiOptions { + input_band: Option, + compute_edges: bool, + output_format: Option, + additional_options: CslStringList, +} + +impl TpiOptions { + /// Create a DEM-topographic-position-index options set. + pub fn new() -> Self { + Default::default() + } + + common_dem_options!(); + + /// Render options into [`CslStringList`] values, as compatible with + /// [`gdal_sys::GDALDEMProcessing`]. + pub fn to_options_list(&self) -> CslStringList { + let mut opts = CslStringList::new(); + self.store_common_options_to(&mut opts); + opts + } +} + +#[cfg(test)] +mod tests { + use crate::assert_near; + use crate::cpl::CslStringList; + use crate::errors::Result; + use crate::raster::processing::dem::topographic_position_index; + use crate::raster::processing::dem::tpi::TpiOptions; + use crate::raster::StatisticsAll; + use crate::test_utils::{fixture, target}; + use crate::Dataset; + + #[test] + fn test_options() -> Result<()> { + let mut proc = TpiOptions::new(); + proc.with_input_band(2.try_into().unwrap()) + .with_compute_edges(true) + .with_output_format("GTiff") + .with_additional_options("CPL_DEBUG=ON".parse()?); + + let expected: CslStringList = "-compute_edges -b 2 -of GTiff CPL_DEBUG=ON".parse()?; + assert_eq!(expected.to_string(), proc.to_options_list().to_string()); + + Ok(()) + } + #[test] + fn test_tpi() -> Result<()> { + let opts = TpiOptions::new(); + + let ds = Dataset::open(fixture("dem-hills.tiff"))?; + + let slope = topographic_position_index(&ds, target("dem-hills-slope.tiff"), &opts)?; + + 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(()) + } +} diff --git a/src/raster/processing/dem/tri.rs b/src/raster/processing/dem/tri.rs new file mode 100644 index 00000000..57242fd8 --- /dev/null +++ b/src/raster/processing/dem/tri.rs @@ -0,0 +1,128 @@ +use std::fmt::{Display, Formatter}; +use std::num::NonZeroUsize; + +use crate::cpl::CslStringList; +use crate::raster::processing::dem::options::common_dem_options; + +/// Configuration options for [`terrain_ruggedness_index()`][super::terrain_ruggedness_index()]. +#[derive(Debug, Clone, Default)] +pub struct TriOptions { + input_band: Option, + compute_edges: bool, + output_format: Option, + additional_options: CslStringList, + algorithm: Option, +} + +impl TriOptions { + /// Create a DEM-terrain-ruggedness-index options set. + pub fn new() -> Self { + Default::default() + } + + common_dem_options!(); + + /// Specify the slope computation algorithm. + pub fn with_algorithm(&mut self, algorithm: DemTriAlg) -> &mut Self { + self.algorithm = Some(algorithm); + self + } + + /// Render relevant common options into [`CslStringList`] values, as compatible with + /// [`gdal_sys::GDALDEMProcessing`]. + pub fn to_options_list(&self) -> CslStringList { + let mut opts = CslStringList::default(); + + self.store_common_options_to(&mut opts); + + // Before 3.3, Wilson is the only algorithm and therefore there's no + // selection option. Rust caller can still specify Wilson, but + // we don't pass it along. + #[cfg(all(major_is_3, minor_ge_3))] + if let Some(alg) = self.algorithm { + opts.add_string("-alg").unwrap(); + opts.add_string(&alg.to_string()).unwrap(); + } + + opts + } +} + +/// Algorithm for computing Terrain Ruggedness Index (TRI). +#[derive(Debug, Clone, Copy)] +pub enum DemTriAlg { + /// The Wilson (see Wilson et al 2007, Marine Geodesy 30:3-35) algorithm uses the mean + /// difference between a central pixel and its surrounding cells. + /// This is recommended for bathymetric use cases. + Wilson, + #[cfg(all(major_is_3, minor_ge_3))] + /// The Riley algorithm (see Riley, S.J., De Gloria, S.D., Elliot, R. (1999): + /// A Terrain Ruggedness that Quantifies Topographic Heterogeneity. + /// Intermountain Journal of Science, Vol.5, No.1-4, pp.23-27) uses the square root of the + /// sum of the square of the difference between a central pixel and its surrounding cells. + /// This is recommended for terrestrial use cases. + /// + /// Only available in GDAL >= 3.3 + Riley, +} + +impl Display for DemTriAlg { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + f.write_fmt(format_args!("{self:?}")) + } +} + +#[cfg(test)] +mod tests { + use crate::assert_near; + use crate::errors::Result; + use crate::raster::processing::dem::terrain_ruggedness_index; + use crate::raster::StatisticsAll; + use crate::test_utils::{fixture, target}; + use crate::Dataset; + + use super::*; + + #[cfg(all(major_is_3, minor_ge_3))] + #[test] + fn test_options() -> Result<()> { + use crate::cpl::CslStringList; + let mut opts = TriOptions::new(); + opts.with_input_band(2.try_into().unwrap()) + .with_compute_edges(true) + .with_algorithm(DemTriAlg::Wilson) + .with_output_format("GTiff") + .with_additional_options("CPL_DEBUG=ON".parse()?); + + let expected: CslStringList = + "-compute_edges -b 2 -of GTiff CPL_DEBUG=ON -alg Wilson".parse()?; + assert_eq!(expected.to_string(), opts.to_options_list().to_string()); + + Ok(()) + } + + #[test] + fn test_tri() -> Result<()> { + let mut opts = TriOptions::new(); + opts.with_algorithm(DemTriAlg::Wilson); + + let ds = Dataset::open(fixture("dem-hills.tiff"))?; + + let tri = terrain_ruggedness_index(&ds, target("dem-hills-tri.tiff"), &opts)?; + + let stats = tri.rasterband(1)?.get_statistics(true, false)?.unwrap(); + + // These numbers were generated by extracting the output from: + // gdaldem tri -alg Wilson fixtures/dem-hills.tiff target/dest.tiff + // gdalinfo -stats target/dest.tiff + let expected = StatisticsAll { + min: 0.0, + max: 4.9836235046387, + mean: 0.49063101456532, + std_dev: 0.67193563366948, + }; + + assert_near!(StatisticsAll, stats, expected, epsilon = 1e-10); + Ok(()) + } +} diff --git a/src/raster/processing/mod.rs b/src/raster/processing/mod.rs new file mode 100644 index 00000000..6ec54fb8 --- /dev/null +++ b/src/raster/processing/mod.rs @@ -0,0 +1,3 @@ +//! GDAL processing routines. + +pub mod dem; diff --git a/src/raster/rasterband.rs b/src/raster/rasterband.rs index dbb540df..20748605 100644 --- a/src/raster/rasterband.rs +++ b/src/raster/rasterband.rs @@ -96,7 +96,7 @@ impl Dataset { /// /// # Example /// -/// ```rust +/// ```rust, no_run /// use gdal::Dataset; /// # fn main() -> gdal::errors::Result<()> { /// use gdal::raster::ResampleAlg; diff --git a/src/test_utils.rs b/src/test_utils.rs index b28dbeb6..8852fe42 100644 --- a/src/test_utils.rs +++ b/src/test_utils.rs @@ -49,12 +49,19 @@ impl AsRef for TempFixture { } /// Returns the fully qualified path to `filename` in `${CARGO_MANIFEST_DIR}/fixtures`. -pub(crate) fn fixture(filename: &str) -> PathBuf { +pub fn fixture(filename: &str) -> PathBuf { Path::new(env!("CARGO_MANIFEST_DIR")) .join("fixtures") .join(filename) } +/// Returns the fully qualified path to `filename` in `${CARGO_MANIFEST_DIR}/target`. +pub fn target(filename: &str) -> PathBuf { + Path::new(env!("CARGO_MANIFEST_DIR")) + .join("target") + .join(filename) +} + /// Scoped value for temporarily suppressing thread-local GDAL log messages. /// /// Useful for tests that expect GDAL errors and want to keep the output log clean @@ -107,3 +114,54 @@ pub fn open_gpkg_for_update(path: &Path) -> (TempPath, Dataset) { .unwrap(); (temp_path, ds) } + +/// Assert numerical difference between two expressions is less than +/// 64-bit machine epsilon or a specified epsilon. +/// +/// # Examples: +/// ```rust, no_run +/// use gdal::assert_near; +/// use std::f64::consts::{PI, E}; +/// assert_near!(PI / E, 1.1557273497909217); +/// // with specified epsilon +/// assert_near!(PI / E, 1.15572734, epsilon = 1e-8); +/// ``` +#[macro_export] +macro_rules! assert_near { + ($left:expr, $right:expr) => { + assert_near!($left, $right, epsilon = f64::EPSILON) + }; + ($left:expr, $right:expr, epsilon = $ep:expr) => { + assert!( + ($left - $right).abs() < $ep, + "|{} - {}| = {} is greater than epsilon {:.4e}", + $left, + $right, + ($left - $right).abs(), + $ep + ) + }; + ($left:expr, $right:expr, epsilon = $ep:expr, field = $field:expr) => { + assert!( + ($left - $right).abs() < $ep, + "field {}: |{} - {}| = {} is greater than epsilon {:.4e}", + $field, + $left, + $right, + ($left - $right).abs(), + $ep + ) + }; + // Pseudo-specialization + (StatisticsAll, $left:expr, $right:expr, epsilon = $ep:expr) => { + assert_near!($left.min, $right.min, epsilon = $ep, field = "min"); + assert_near!($left.max, $right.max, epsilon = $ep, field = "max"); + assert_near!($left.mean, $right.mean, epsilon = $ep, field = "mean"); + assert_near!( + $left.std_dev, + $right.std_dev, + epsilon = $ep, + field = "std_dev" + ); + }; +} diff --git a/src/vsi.rs b/src/vsi.rs index b09e9998..2b5b9a1e 100644 --- a/src/vsi.rs +++ b/src/vsi.rs @@ -2,7 +2,6 @@ //! //! This module provides safe access to a subset of the [GDAL VSI Functions](https://gdal.org/doxygen/cpl__vsi_8h.html). //! See [GDAL Virtual File Systems document](https://gdal.org/user/virtual_file_systems.html) for details. -//! use std::marker::PhantomData; use std::mem::ManuallyDrop; @@ -15,59 +14,57 @@ use crate::utils::{_last_null_pointer_err, _path_to_c_string, _pathbuf_array}; /// Read the file names from a virtual file system with optional recursion. pub fn read_dir>(path: P, recursive: bool) -> Result> { - _read_dir(path.as_ref(), recursive) -} - -fn _read_dir(path: &Path, recursive: bool) -> Result> { - let path = _path_to_c_string(path)?; - let data = if recursive { - let data = unsafe { gdal_sys::VSIReadDirRecursive(path.as_ptr()) }; - if data.is_null() { - return Err(_last_null_pointer_err("VSIReadDirRecursive")); - } - data - } else { - let data = unsafe { gdal_sys::VSIReadDir(path.as_ptr()) }; - if data.is_null() { - return Err(_last_null_pointer_err("VSIReadDir")); - } - data - }; + fn _read_dir(path: &Path, recursive: bool) -> Result> { + let path = _path_to_c_string(path)?; + let data = if recursive { + let data = unsafe { gdal_sys::VSIReadDirRecursive(path.as_ptr()) }; + if data.is_null() { + return Err(_last_null_pointer_err("VSIReadDirRecursive")); + } + data + } else { + let data = unsafe { gdal_sys::VSIReadDir(path.as_ptr()) }; + if data.is_null() { + return Err(_last_null_pointer_err("VSIReadDir")); + } + data + }; - Ok(_pathbuf_array(data)) + Ok(_pathbuf_array(data)) + } + _read_dir(path.as_ref(), recursive) } /// Creates a new VSIMemFile from a given buffer. pub fn create_mem_file>(file_name: P, data: Vec) -> Result<()> { - _create_mem_file(file_name.as_ref(), data) -} - -fn _create_mem_file(file_name: &Path, data: Vec) -> Result<()> { - let file_name = _path_to_c_string(file_name)?; - - // ownership will be given to GDAL, so it should not be automaticly dropped - let mut data = ManuallyDrop::new(data); - - let handle = unsafe { - VSIFileFromMemBuffer( - file_name.as_ptr(), - data.as_mut_ptr(), - data.len() as u64, - true as i32, - ) - }; + fn _create_mem_file(file_name: &Path, data: Vec) -> Result<()> { + let file_name = _path_to_c_string(file_name)?; + + // ownership will be given to GDAL, so it should not be automaticly dropped + let mut data = ManuallyDrop::new(data); + + let handle = unsafe { + VSIFileFromMemBuffer( + file_name.as_ptr(), + data.as_mut_ptr(), + data.len() as u64, + true as i32, + ) + }; + + if handle.is_null() { + // on error, allow dropping the data again + ManuallyDrop::into_inner(data); + return Err(_last_null_pointer_err("VSIGetMemFileBuffer")); + } - if handle.is_null() { - // on error, allow dropping the data again - ManuallyDrop::into_inner(data); - return Err(_last_null_pointer_err("VSIGetMemFileBuffer")); - } + unsafe { + VSIFCloseL(handle); + } - unsafe { - VSIFCloseL(handle); + Ok(()) } - - Ok(()) + _create_mem_file(file_name.as_ref(), data) } /// A helper struct that unlinks a mem file that points to borrowed data @@ -100,77 +97,77 @@ pub fn create_mem_file_from_ref>( file_name: P, data: &mut [u8], ) -> Result> { - _create_mem_file_from_ref(file_name.as_ref(), data) -} - -fn _create_mem_file_from_ref<'d>(file_name: &Path, data: &'d mut [u8]) -> Result> { - let file_name_c = _path_to_c_string(file_name)?; - - let handle = unsafe { - VSIFileFromMemBuffer( - file_name_c.as_ptr(), - data.as_mut_ptr(), - data.len() as u64, - false as i32, - ) - }; + fn _create_mem_file_from_ref<'d>( + file_name: &Path, + data: &'d mut [u8], + ) -> Result> { + let file_name_c = _path_to_c_string(file_name)?; + + let handle = unsafe { + VSIFileFromMemBuffer( + file_name_c.as_ptr(), + data.as_mut_ptr(), + data.len() as u64, + false as i32, + ) + }; + + if handle.is_null() { + return Err(_last_null_pointer_err("VSIGetMemFileBuffer")); + } - if handle.is_null() { - return Err(_last_null_pointer_err("VSIGetMemFileBuffer")); - } + unsafe { + VSIFCloseL(handle); + } - unsafe { - VSIFCloseL(handle); + Ok(MemFileRef::new(file_name)) } - - Ok(MemFileRef::new(file_name)) + _create_mem_file_from_ref(file_name.as_ref(), data) } /// Unlink a VSIMemFile. pub fn unlink_mem_file>(file_name: P) -> Result<()> { - _unlink_mem_file(file_name.as_ref()) -} + fn _unlink_mem_file(file_name: &Path) -> Result<()> { + let file_name_c = _path_to_c_string(file_name)?; -fn _unlink_mem_file(file_name: &Path) -> Result<()> { - let file_name_c = _path_to_c_string(file_name)?; + let rv = unsafe { VSIUnlink(file_name_c.as_ptr()) }; - let rv = unsafe { VSIUnlink(file_name_c.as_ptr()) }; + if rv != 0 { + return Err(GdalError::UnlinkMemFile { + file_name: file_name.display().to_string(), + }); + } - if rv != 0 { - return Err(GdalError::UnlinkMemFile { - file_name: file_name.display().to_string(), - }); + Ok(()) } - - Ok(()) + _unlink_mem_file(file_name.as_ref()) } /// Copies the bytes of the VSIMemFile with given `file_name`. /// Takes the ownership and frees the memory of the VSIMemFile. pub fn get_vsi_mem_file_bytes_owned>(file_name: P) -> Result> { - _get_vsi_mem_file_bytes_owned(file_name.as_ref()) -} - -fn _get_vsi_mem_file_bytes_owned(file_name: &Path) -> Result> { - let file_name = _path_to_c_string(file_name)?; + fn _get_vsi_mem_file_bytes_owned(file_name: &Path) -> Result> { + let file_name = _path_to_c_string(file_name)?; - let owned_bytes = unsafe { - let mut length: u64 = 0; - let bytes = VSIGetMemFileBuffer(file_name.as_ptr(), &mut length, true as i32); + let owned_bytes = unsafe { + let mut length: u64 = 0; + let bytes = VSIGetMemFileBuffer(file_name.as_ptr(), &mut length, true as i32); - if bytes.is_null() { - return Err(_last_null_pointer_err("VSIGetMemFileBuffer")); - } + if bytes.is_null() { + return Err(_last_null_pointer_err("VSIGetMemFileBuffer")); + } - let slice = std::slice::from_raw_parts(bytes, length as usize); - let vec = slice.to_vec(); + let slice = std::slice::from_raw_parts(bytes, length as usize); + let vec = slice.to_vec(); - VSIFree(bytes.cast::()); + VSIFree(bytes.cast::()); - vec - }; + vec + }; - Ok(owned_bytes) + Ok(owned_bytes) + } + _get_vsi_mem_file_bytes_owned(file_name.as_ref()) } /// Computes a function on the bytes of the vsi in-memory file with given `file_name`. @@ -179,27 +176,18 @@ pub fn call_on_mem_file_bytes>(file_name: P, fun: F) -> Res where F: FnOnce(&[u8]) -> R, { - _call_on_mem_file_bytes(file_name.as_ref(), fun) -} - -fn _call_on_mem_file_bytes(file_name: &Path, fun: F) -> Result -where - F: FnOnce(&[u8]) -> R, -{ - let file_name = _path_to_c_string(file_name)?; + let file_name = _path_to_c_string(file_name.as_ref())?; - unsafe { - let mut length: u64 = 0; - let bytes = VSIGetMemFileBuffer(file_name.as_ptr(), &mut length, false as i32); + let mut length: u64 = 0; + let bytes = unsafe { VSIGetMemFileBuffer(file_name.as_ptr(), &mut length, false as i32) }; - if bytes.is_null() { - return Err(_last_null_pointer_err("VSIGetMemFileBuffer")); - } + if bytes.is_null() { + return Err(_last_null_pointer_err("VSIGetMemFileBuffer")); + } - let slice = std::slice::from_raw_parts(bytes, length as usize); + let slice = unsafe { std::slice::from_raw_parts(bytes, length as usize) }; - Ok(fun(slice)) - } + Ok(fun(slice)) } #[cfg(test)] diff --git a/tests/driver-no-auto-register.rs b/tests/driver-no-auto-register.rs index ab8e67ef..a0b51694 100644 --- a/tests/driver-no-auto-register.rs +++ b/tests/driver-no-auto-register.rs @@ -1,6 +1,12 @@ -mod utils; - use gdal::{Dataset, DriverManager}; +use std::path::{Path, PathBuf}; + +/// Returns the fully qualified path to `filename` in `${CARGO_MANIFEST_DIR}/fixtures`. +fn fixture(filename: &str) -> PathBuf { + Path::new(env!("CARGO_MANIFEST_DIR")) + .join("fixtures") + .join(filename) +} #[test] /// Sequentially run tests @@ -16,21 +22,21 @@ fn test_manually_registering_drivers() { assert_eq!(DriverManager::count(), 0); - assert!(Dataset::open(fixture!("tinymarble.tif")).is_err()); + assert!(Dataset::open(fixture("tinymarble.tif")).is_err()); DriverManager::register_all(); - assert!(Dataset::open(fixture!("tinymarble.tif")).is_ok()); + assert!(Dataset::open(fixture("tinymarble.tif")).is_ok()); let driver = DriverManager::get_driver_by_name("GTiff").unwrap(); DriverManager::deregister_driver(&driver); - assert!(Dataset::open(fixture!("tinymarble.tif")).is_err()); + assert!(Dataset::open(fixture("tinymarble.tif")).is_err()); DriverManager::register_driver(&driver); - assert!(Dataset::open(fixture!("tinymarble.tif")).is_ok()); + assert!(Dataset::open(fixture("tinymarble.tif")).is_ok()); } fn test_deregister_all_but_one() { diff --git a/tests/utils.rs b/tests/utils.rs deleted file mode 100644 index 29e6368a..00000000 --- a/tests/utils.rs +++ /dev/null @@ -1,14 +0,0 @@ -#[macro_export] -macro_rules! fixture { - ($name:expr) => { - std::path::Path::new(file!()) - .parent() - .unwrap() - .parent() - .unwrap() - .join("fixtures") - .as_path() - .join($name) - .as_path() - }; -}