Implement excesstopography

NAMESPACE

@@ -22,6 +22,7 @@ export(GRIDobj2mat)
 export(GRIDobj2pm)
 export(GRIDobj2xyz)
 export(and)
+export(excesstopography)
 export(ezgetnal)
 export(fillsinks)
 export(flow_accumulation)

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)
+}

src/wrap_excesstopography.c

@@ -0,0 +1,15 @@
+#include <R.h>
+#include <stddef.h>
+#include <stdint.h>
+
+#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);
+}

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))
+})