Merging scenario fixes #113, threshold constraints and dispersal corr…

…ections (#114) * 🐛 fix - transformed parameters added to predictor dataset #113 * Updated scripts for normalizing #113 * Further 🔥 fixes and tests #113 * Small but substantial speedup to stars to raster conversion * Improvements to dispersal constraint function * Addition of threshold constraint * Small improvements to handling of projected thresholds * Small fix remaining #113
iiasa · Mar 31, 2024 · bfbeb51 · bfbeb51
1 parent f0336bb
commit bfbeb51
Show file tree

Hide file tree

Showing 31 changed files with 589 additions and 101 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -76,6 +76,7 @@ Suggests:
     progress,
     pdp,
     rmarkdown,
+    lubridate (>= 1.9.0),
     rstan (>= 2.21.0),
     rstantools (>= 2.1.1),
     testthat (>= 3.0.0),

diff --git a/NAMESPACE b/NAMESPACE
@@ -63,6 +63,7 @@ export(add_constraint_boundary)
 export(add_constraint_connectivity)
 export(add_constraint_dispersal)
 export(add_constraint_minsize)
+export(add_constraint_threshold)
 export(add_control_bias)
 export(add_control_extrapolation)
 export(add_latent_spatial)

diff --git a/NEWS.md b/NEWS.md
@@ -11,7 +11,8 @@
 * Small bug fixes dealing with `scenario()` projections and limits, plus unit tests #104
 * Bug fixes with adding `predictor_derivate()` to scenario predictors and added unit tests #106
 * Several fixes related to different engines and priors
-# Changed default output for netcdf files to multidimensional arrays #109
+* Changed default output for netcdf files to multidimensional arrays #109
+* :fire: hot fixes for scenario scaling and normalization issue #113
 
 # ibis.iSDM 0.1.2
 

diff --git a/R/add_constraint.R b/R/add_constraint.R
@@ -57,6 +57,8 @@ NULL
 #' * \code{minsize} - Allows to specify a certain size that must be satisfied in
 #' order for a thresholded patch to be occupied. Can be thought of as a minimum
 #' size requirement. See `add_constraint_minsize()` for the required parameters.
+#' * \code{threshold} - Applies the set threshold as a constrain directly on the
+#' suitability projections. Requires a threshold to be set.
 #'
 #' @returns Adds constraints data to a [`BiodiversityScenario`] object.
 #'
@@ -105,7 +107,7 @@ methods::setMethod(
     method <- match.arg(arg = method,
                         choices = c("sdd_fixed", "sdd_nexpkernel", "kissmig", "migclim",
                                     "hardbarrier","resistance",
-                                    "boundary", "minsize",
+                                    "boundary", "minsize", "threshold",
                                     "nichelimit"), several.ok = FALSE)
 
     # Now call the respective functions individually
@@ -127,6 +129,8 @@ methods::setMethod(
                   # --- #
                   "boundary" = add_constraint_boundary(mod, ...),
                   # --- #
+                  "threshold" = add_constraint_threshold(mod, ...),
+                  # --- #
                   "minsize" = add_constraint_minsize(mod, ...)
     )
     return(o)
@@ -153,6 +157,9 @@ methods::setMethod(
 #' Parameters for \code{'method'}:
 #' * \code{sdd_fixed} - Applies a fixed uniform dispersal distance per modelling timestep.
 #' * \code{sdd_nexpkernel} - Applies a dispersal distance using a negative exponential kernel from its origin.
+#' #' The negative exponential kernel is defined as:
+#' \deqn{f(x) = \frac{1}{2 \pi a^2} e^{-\frac{x}{a}}}{fx = 1 / (2 * pi * a^2) * exp(-x / a)}
+#' where \eqn{a} is the mean dispersal distance (in m) divided by 2.
 #' * \code{kissmig} - Applies the kissmig stochastic dispersal model. Requires \code{`kissmig`} package. Applied at each modelling time step.
 #' * \code{migclim} - Applies the dispersal algorithm MigClim to the modelled objects. Requires \code{"MigClim"} package.
 #'
@@ -166,6 +173,9 @@ methods::setMethod(
 #' * \code{pcor}: [`numeric`] probability that corner cells are considered in the
 #' 3x3 neighbourhood (Default: \code{0.2}).
 #'
+#' @note
+#' Unless otherwise stated, the default unit of supplied distance values (e.g. average dispersal
+#' distance) should be in \code{"m"}.
 #' @references
 #' * Bateman, B. L., Murphy, H. T., Reside, A. E., Mokany, K., & VanDerWal, J. (2013).
 #' Appropriateness of full‐, partial‐and no‐dispersal scenarios in climate change impact
@@ -283,14 +293,15 @@ methods::setMethod(
     is.Raster(baseline_threshold), is.Raster(new_suit),
     is_comparable_raster(baseline_threshold, new_suit),
     is.numeric(value),
-    is.null(resistance) || is.Raster(resistance),
-    # Check that baseline threshold raster is binomial
-    length(unique(baseline_threshold))==2
+    is.null(resistance) || is.Raster(resistance)
   )
 
-  # Get original baseline threshold
-  ori.tr <- baseline_threshold
-  ori.tr[ori.tr>0] <- 0
+  # Check for small lon-lat values
+  if(terra::is.lonlat(baseline_threshold)){
+    if(value < 1){
+      message('Very small average dispersal vlaue provided. Check that they are in unit m!')
+    }
+  }
 
   # Set resistance layer to 0 if set to zero.
   if(is.Raster(resistance)){
@@ -309,8 +320,6 @@ methods::setMethod(
   # Now multiply the net suitability projection with this mask
   # Thus removing any grid cells outside
   out <- new_suit * ras_dis
-  # Mask with original so as to retain non-zero values
-  out <- terra::mask(out, ori.tr)
   return(out)
 }
 
@@ -328,19 +337,23 @@ methods::setMethod(
 #' @noRd
 #'
 #' @keywords internal
-.sdd_nexpkernel <- function(baseline_threshold, new_suit, value, normalize = FALSE, resistance = NULL){
+.sdd_nexpkernel <- function(baseline_threshold, new_suit, value, normalize = TRUE, resistance = NULL){
   assertthat::assert_that(
     is.Raster(baseline_threshold), is.Raster(new_suit),
     is_comparable_raster(baseline_threshold, new_suit),
     is.numeric(value),
+    is.logical(normalize),
     is.null(resistance) || is.Raster(resistance),
     # Check that baseline threshold raster is binomial
     length(unique(baseline_threshold)[,1])==2
   )
 
-  # Get original baseline threshold
-  ori.tr <- baseline_threshold
-  ori.tr[ori.tr>0] <- 0
+  # Check for small lon-lat values
+  if(terra::is.lonlat(baseline_threshold)){
+    if(value < 1){
+      message('Very small average dispersal vlaue provided. Check that they are in unit m!')
+    }
+  }
 
   # Set resistance layer to 0 if set to zero.
   if(is.Raster(resistance)){
@@ -350,23 +363,22 @@ methods::setMethod(
     baseline_threshold <- terra::mask(baseline_threshold, resistance)
   }
 
-  # Inverse of mean dispersal distance
-  alpha <- 1/value
+  # Divide alpha values by 2
+  alpha <- value/2
 
   # Grow baseline raster by using an exponentially weighted kernel
   ras_dis <- terra::gridDist(baseline_threshold, target = 1)
+  # Normalized (with a constant) negative exponential kernel
+  ras_dis <- terra::app(ras_dis, fun = function(x) (1 / (2 * pi * value ^ 2)) * exp(-x / value) )
+  # Equivalent to alpha = 1/value and
+  # ras_dis <- terra::app(ras_dis, fun = function(x) exp(-alpha * x))
   if(normalize){
-    # Normalized (with a constant) negative exponential kernel
-    ras_dis <- terra::app(ras_dis, fun = function(x) (1 / (2 * pi * value ^ 2)) * exp(-x / value) )
-  } else {
-    ras_dis <- terra::app(ras_dis, fun = function(x) exp(-alpha * x))
+    ras_dis <- predictor_transform(ras_dis, option = 'norm')
   }
 
   # Now multiply the net suitability projection with this mask Thus removing any
   # non-suitable grid cells (0) and changing the value of those within reach
   out <- new_suit * ras_dis
-  # Mask with original so as to retain non-zero values
-  out <- terra::mask(out, ori.tr)
   return(out)
 }
 
@@ -604,7 +616,7 @@ methods::setMethod(
         value > 0, msg = "Specify a value threshold (SD) and names of predictors, for which
         we do not expect the species to persist."
       )
-      if(is.null(names)) names <- character()
+      if(is.null(names)) names <- NA
       co[['adaptability']] <- list(method = method,
                                    params = c("names" = names, "value" = value,
                                               "increment" = increment))
@@ -654,9 +666,9 @@ methods::setMethod(
   for(id in names(model$biodiversity)){
     sub <- model$biodiversity[[id]]
     # Which are presence data
-    is_presence <- which(sub$observations[['observed']] > 0)
+    is_presence <- which(factor_to_numeric(sub$observations[['observed']]) > 0)
     df <- rbind(df,
-                sub$predictors[is_presence, names])
+                sub$predictors[is_presence, sub$predictors_names])
   }
   rr <- sapply(df, function(x) range(x, na.rm = TRUE))   # Calculate ranges
   rsd <- sapply(df, function(x) stats::sd(x, na.rm = TRUE))   # Calculate standard deviation
@@ -668,6 +680,7 @@ methods::setMethod(
   # Now 'clamp' all predictor values beyond these names to 0, e.g. partial out
   nd <- newdata
   for(n in names){
+    if(!(n %in% names(rr))) next() # If variable not present in model frame, skip
     # Calc min
     min_ex <- which(nd[,n] < rr[1,n])
     max_ex <- which(nd[,n] > rr[2,n])
@@ -876,3 +889,66 @@ methods::setMethod(
     return( new )
   }
 )
+
+#' Adds a threshold constraint to a scenario object
+#'
+#' @description
+#' This option adds a [`threshold()`] constraint to a scenario projection,
+#' thus effectively applying the threshold as mask to each projection step made
+#' during the scenario projection.
+#'
+#' Applying this constraint thus means that the \code{"suitability"} projection is
+#' clipped to the threshold. This method requires
+#' the `threshold()` set for a scenario object.
+#'
+#' It could be in theory possible to re calculate the threshold for each time step
+#' based on supplied parameters or even observation records. So far this option has
+#' not been necessary to implement.
+#'
+#' @note
+#' Threshold values are taken from the original fitted model.
+#'
+#' @inheritParams add_constraint
+#' @param updatevalue A [`numeric`] indicating to what the masked out values (those outside)
+#' the threshold should become (Default: \code{NA}).
+#'
+#' @family constraint
+#' @keywords scenario
+#'
+#' @examples
+#' \dontrun{
+#' # Add scenario constraint
+#' scenario(fit) |> threshold() |>
+#' add_constraint_threshold()
+#' }
+#'
+#' @name add_constraint_threshold
+NULL
+
+#' @rdname add_constraint_threshold
+#' @export
+methods::setGeneric("add_constraint_threshold",
+                    signature = methods::signature("mod"),
+                    function(mod, updatevalue = NA, ...) standardGeneric("add_constraint_threshold"))
+
+#' @rdname add_constraint_threshold
+methods::setMethod(
+  "add_constraint_threshold",
+  methods::signature(mod = "BiodiversityScenario"),
+  function(mod, updatevalue = NA, ...){
+    assertthat::assert_that(
+      inherits(mod, "BiodiversityScenario"),
+      is.na(updatevalue) || is.numeric(updatevalue)
+    )
+
+    # Add processing method #
+    # --- #
+    co <- list()
+    co[['threshold']] <- list(method = "threshold",
+                             params = c("updatevalue" = updatevalue))
+    # --- #
+    new <- mod$clone(deep = TRUE)
+    new <- new$set_constraints(co)
+    return( new )
+  }
+)