From 70afa59972cf5697299944e892d77f07d6202e76 Mon Sep 17 00:00:00 2001 From: William Kearney Date: Thu, 5 Mar 2026 16:18:32 +0100 Subject: [PATCH] Implement excesstopography Resolves #40 This PR adds 1. A binding for libtopotoolbox's `excesstopography_fsm2d` 2. An R function `excesstopography` 3. A test of `excesstopography` As in MATLAB/Python, `excesstopography` accepts either a GRIDobj or a scalar for the threshold_slope parameter. No real validation, such as checking the alignment the DEM and the threshold slope, is performed. The test checks that excesstopography runs successfully and validates that the reconstructed topography lies below the supplied DEM. Signed-off-by: William Kearney --- NAMESPACE | 1 + R/excesstopography.R | 69 ++++++++++++++++++++++++++ man/excesstopography.Rd | 44 ++++++++++++++++ src/wrap_excesstopography.c | 15 ++++++ tests/testthat/test-excesstopography.R | 12 +++++ 5 files changed, 141 insertions(+) create mode 100644 R/excesstopography.R create mode 100644 man/excesstopography.Rd create mode 100644 src/wrap_excesstopography.c create mode 100644 tests/testthat/test-excesstopography.R diff --git a/NAMESPACE b/NAMESPACE index 6f14efb..8e4893b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -22,6 +22,7 @@ export(GRIDobj2mat) export(GRIDobj2pm) export(GRIDobj2xyz) export(and) +export(excesstopography) export(ezgetnal) export(fillsinks) export(flow_accumulation) diff --git a/R/excesstopography.R b/R/excesstopography.R new file mode 100644 index 0000000..3353ca6 --- /dev/null +++ b/R/excesstopography.R @@ -0,0 +1,69 @@ +#' Excess topography +#' +#' Reconstruct surface with threshold-slope surface +#' +#' The excess topography (Blöthe et al. 2015) is computed by solving +#' an eikonal equation (Anand et al. 2023) constrained to lie below +#' the original DEM. Where the slope of the DEM is greater than the +#' threshold slope, the eikonal solver limits the output topography to +#' that slope, but where the slope of the DEM is lower that the +#' threshold slope, the output follows the DEM. +#' +#' The eikonal equation is solved using the fast sweeping method (Zhao +#' 2004), which iterates over the DEM in alternating directions and +#' updates the topography according to an upwind discretization of the +#' gradient. To constrain the solution by the original DEM, the output +#' topography is initiated with the DEM and only updates lower than +#' the DEM are accepted. +#' +#' @param dem GRIDobj; Digital elevation model +#' @param threshold_slope float | GRIDobj; The threshold +#' slope. Spatially variable slopes can be provided as a GRIDobj. +#' +#' @return GRIDobj; The reconstructed topography. The excess +#' topography is the difference between the original DEM and the +#' reconstructed one. +#' +#' @examples +#' +#' data(srtm_bigtujunga30m_utm11) +#' DEM <- GRIDobj(srtm_bigtujunga30m_utm11) +#' DEMext <- excesstopography(DEM, tan(20*pi/180)) +#' +#' @export + +excesstopography <- function(dem, threshold_slope) { + # Input validation + dem <- processgrid(dem) + provided_go <- dem$provided_go + dem <- dem$r + + # Preallocate output + ext <- dem + + # Get DEM input + dem <- get_grid_data(dem) + + # Pass grids to libtopotoolbox + if (is.numeric(threshold_slope)) { + threshold_grid <- as.single(rep(threshold_slope, length(dem$z))) + } else { + threshold_grid <- as.single(get_grid_data(processgrid(threshold_slope)$r)$z) + } + output <- single(length(dem$z)) + result <- .C( + "wrap_excesstopography", + outputR = as.single(output), + demR = as.single(dem$z), + thresholdR = threshold_grid, + cellsizeR = as.single(dem$cellsize), + dimsR = as.integer(dem$dims), + NAOK = TRUE + )$outputR + + terra::values(ext) <- result + if (provided_go) { + ext <- GRIDobj(ext) + } + return(ext) +} diff --git a/man/excesstopography.Rd b/man/excesstopography.Rd new file mode 100644 index 0000000..6524d56 --- /dev/null +++ b/man/excesstopography.Rd @@ -0,0 +1,44 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/excesstopography.R +\name{excesstopography} +\alias{excesstopography} +\title{Excess topography} +\usage{ +excesstopography(dem, threshold_slope) +} +\arguments{ +\item{dem}{GRIDobj; Digital elevation model} + +\item{threshold_slope}{float | GRIDobj; The threshold +slope. Spatially variable slopes can be provided as a GRIDobj.} +} +\value{ +GRIDobj; The reconstructed topography. The excess + topography is the difference between the original DEM and the + reconstructed one. +} +\description{ +Reconstruct surface with threshold-slope surface +} +\details{ +The excess topography (Blöthe et al. 2015) is computed by solving +an eikonal equation (Anand et al. 2023) constrained to lie below +the original DEM. Where the slope of the DEM is greater than the +threshold slope, the eikonal solver limits the output topography to +that slope, but where the slope of the DEM is lower that the +threshold slope, the output follows the DEM. + +The eikonal equation is solved using the fast sweeping method (Zhao +2004), which iterates over the DEM in alternating directions and +updates the topography according to an upwind discretization of the +gradient. To constrain the solution by the original DEM, the output +topography is initiated with the DEM and only updates lower than +the DEM are accepted. +} +\examples{ + +data(srtm_bigtujunga30m_utm11) +DEM <- GRIDobj(srtm_bigtujunga30m_utm11) +DEMext <- excesstopography(DEM, tan(20*pi/180)) + +} diff --git a/src/wrap_excesstopography.c b/src/wrap_excesstopography.c new file mode 100644 index 0000000..d2ca63f --- /dev/null +++ b/src/wrap_excesstopography.c @@ -0,0 +1,15 @@ +#include +#include +#include + +#include "topotoolbox.h" +#include "topotoolboxr.h" + +void wrap_excesstopography(float *excess, float *dem, float *threshold_slopes, + float *cellsize, int *dimsR) { + // Transformation of integers and array allocation + ptrdiff_t dims[2] = {dimsR[0], dimsR[1]}; + + // Call libtopotoolbox + excesstopography_fsm2d(excess, dem, threshold_slopes, *cellsize, dims); +} diff --git a/tests/testthat/test-excesstopography.R b/tests/testthat/test-excesstopography.R new file mode 100644 index 0000000..984c943 --- /dev/null +++ b/tests/testthat/test-excesstopography.R @@ -0,0 +1,12 @@ +test_that("Tests excesstopography", { + dem <- terra::rast(system.file("ex/elev.tif", package = "terra")) + dem <- terra::project(dem, "EPSG:32632", res = 90.0) + expect_silent(dem_ext <- excesstopography(GRIDobj(dem), tan(20 * pi / 180))) + expect_true(inherits(dem_ext, "GRIDobj")) + expect_true(all(terra::values(dem_ext$raster) <= terra::values(dem), na.rm=TRUE)) + + ts <- dem + terra::values(ts) <- runif(terra::size(ts)) + expect_silent(dem_ext <- excesstopography(GRIDobj(dem), GRIDobj(ts))) + expect_true(all(terra::values(dem_ext$raster) <= terra::values(dem), na.rm=TRUE)) +})