Work on scenario improvements and mechanistic functions #67 #68

iiasa · Oct 1, 2023 · fb173c3 · fb173c3
1 parent dac8f5f
commit fb173c3
Show file tree

Hide file tree

Showing 21 changed files with 788 additions and 51 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -142,6 +142,7 @@ Collate:
     'pseudoabsence.R'
     'scenario.R'
     'similarity.R'
+    'simulate_population_steps.R'
     'summary.R'
     'threshold.R'
     'train.R'

diff --git a/NAMESPACE b/NAMESPACE
@@ -139,6 +139,7 @@ export(sanitize_names)
 export(scenario)
 export(sel_predictors)
 export(similarity)
+export(simulate_population_steps)
 export(spartial)
 export(thin_observations)
 export(threshold)
@@ -203,6 +204,7 @@ exportMethods(rm_priors)
 exportMethods(scenario)
 exportMethods(sel_predictors)
 exportMethods(similarity)
+exportMethods(simulate_population_steps)
 exportMethods(threshold)
 exportMethods(train)
 exportMethods(validate)

diff --git a/NEWS.md b/NEWS.md
@@ -1,12 +1,15 @@
 # ibis.iSDM 0.1.0 (current dev branch)
+
 #### New features
 * Added a small convenience wrapper to add model outputs to another model `add_predictors_model()`
-* Started and added mechanistic SDM vignette #67
+* Started adding mechanistic SDM vignette #67
+* Wrapper for *steps* implemented via `simulate_population_steps()` #68
 
 #### Minor improvements and bug fixes
-
+* Minor bug fixes in `scenario()` object, and MigClim and Kissmig wrappers.
 
 # ibis.iSDM 0.0.9
+
 #### New features
 * Added new vignette on available functions for data preparation #67
 * Addition of small `mask()` function that emulates the for `terra`.
@@ -21,6 +24,7 @@
 * Improved error messages and handling of formula's.
 
 # ibis.iSDM 0.0.8
+
 #### New features
 * Implemented min size constraint (`add_constraint_minsize()`) #56
 * Added a function for estimating partial effects of ensembles `ensemble_spartial()`.
@@ -34,6 +38,7 @@
 * Allow to specify location of biodiversity point records in `threshold()`.
 
 # ibis.iSDM 0.0.7
+
 #### New features
 * Added method proximity to `add_control_bias()` to place lower weights on points closer to another.
 * Added helper functions `get_data()` and the option to apply `threshold()` directly on BiodiversityScenarios.

diff --git a/R/add_constraint.R b/R/add_constraint.R
@@ -61,12 +61,16 @@ NULL
 #'   can be particular useful to limit the influence of increasing marginal
 #'   responses and avoid biologically unrealistic projections.
 #'
-#' **Boundary**:
+#' **Boundary and size**:
 #' * \code{boundary} - Applies a hard boundary constraint on the projection, thus disallowing an expansion of a range outside
 #'   the provide layer. Similar as specifying projection limits (see
 #'   [`distribution`]), but can be used to specifically constrain a projection
 #'   within a certain area (e.g. a species range or an island).
 #'
+#' * \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.
+#'
 #' @returns Adds constraints data to a [`BiodiversityScenario`] object.
 #' @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 modelling. Diversity and Distributions, 19(10), 1224-1234.
@@ -106,7 +110,8 @@ methods::setMethod(
     # Match method
     method <- match.arg(arg = method,
                         choices = c("sdd_fixed", "sdd_nexpkernel", "kissmig", "migclim",
-                                    "hardbarrier","resistance", "boundary",
+                                    "hardbarrier","resistance",
+                                    "boundary", "minsize",
                                     "nichelimit"), several.ok = FALSE)
 
     # Now call the respective functions individually
@@ -126,8 +131,10 @@ methods::setMethod(
                   # --- #
                   "nichelimit" = add_constraint_adaptability(mod, method = "nichelimit", ...),
                   # --- #
-                  "boundary" = add_constraint_boundary(mod, ...)
-                  )
+                  "boundary" = add_constraint_boundary(mod, ...),
+                  # --- #
+                  "minsize" = add_constraint_minsize(mod, ...)
+    )
     return(o)
   }
 )
@@ -352,7 +359,7 @@ methods::setMethod(
     is.vector(params) || is.list(params),
     is.null(resistance) || is.logical(resistance) || is.Raster(resistance),
     # Check that baseline threshold raster is binomial
-    length(unique(baseline_threshold))==2
+    length(unique(baseline_threshold)[,1])==2
   )
 
   check_package('kissmig')
@@ -366,12 +373,15 @@ methods::setMethod(
   # Simulate kissmig for a given threshold and suitability raster
   km <- kissmig::kissmig(O = terra_to_raster( baseline_threshold ),
                          # Rescale newsuit to 0-1
-                         S = predictor_transform(new_suit, 'norm'),
+                         S = predictor_transform(new_suit, 'norm') |>
+                           terra_to_raster(),
                          it = as.numeric( params['iteration'] ),
                          type = params['type'],
                          pext = as.numeric(params['pext']),
                          pcor = as.numeric(params['pcor'])
                         )
+  # Convert to terra again
+  km <- terra::rast(km)
   if(is.factor(km)) km <- terra::as.int(km)
 
   # Now multiply the net suitability projection with this mask Thus removing any

diff --git a/R/add_constraint_MigClim.R b/R/add_constraint_MigClim.R
@@ -63,7 +63,7 @@ NULL
 #'
 #' @name add_constraint_MigClim
 #' @aliases add_constraint_MigClim
-#' @family constrain
+#' @family constraint
 #' @keywords scenario
 #' @exportMethod add_constraint_MigClim
 #' @export

diff --git a/R/bdproto-biodiversityscenario.R b/R/bdproto-biodiversityscenario.R
@@ -48,6 +48,9 @@ BiodiversityScenario <- bdproto(
     # Thresholds
     tr <- self$get_threshold()
 
+    # Any other simulation outputs modules
+    simmods <- self$get_simulation()
+
     message(paste0(
       ifelse(is.Waiver(self$limits),"Spatial-temporal scenario:","Spatial-temporal scenario (limited):"),
                    '\n  Used model: ',ifelse(is.Waiver(fit) || isFALSE(fit), text_red('None'), class(fit)[1] ),
@@ -57,6 +60,7 @@ BiodiversityScenario <- bdproto(
                    ifelse(!is.Waiver(cs)||!is.Waiver(tr), "\n --------- ", ""),
                    ifelse(is.Waiver(cs),"", paste0("\n  Constraints:    ", text_green(paste(paste0(names(cs),' (',cs,')'),collapse = ', ')) ) ),
                    ifelse(is.Waiver(tr),"", paste0("\n  Threshold:      ", round(tr[1], 3),' (',names(tr[1]),')') ),
+                   ifelse(is.Waiver(simmods),"", paste0("\n  Simulations:    ", text_green(paste(paste0(names(simmods),' (',simmods[[1]][[1]],')'),collapse = ', ')) ) ),
                    "\n --------- ",
                    "\n  Scenarios fitted: ", ifelse(is.Waiver(self$scenarios),text_yellow('None'), text_green('Yes'))
       )
@@ -194,6 +198,16 @@ BiodiversityScenario <- bdproto(
       bdproto(NULL, self, constraints = x)
     }
   },
+  # Get simulation options and parameters
+  get_simulation = function(self){
+    if('simulation' %notin% names(self)) return( new_waiver() )
+    return( self$simulation )
+  },
+  # Set simulation objects
+  set_simulation = function(self, x ){
+    # We only take a single simulation so far
+    bdproto(NULL, self, simulation = x)
+  },
   # Get Predictors
   get_predictors = function(self){
     return(self$predictors)
@@ -232,6 +246,8 @@ BiodiversityScenario <- bdproto(
       if(getOption('ibis.setupmessages')) myLog('[Scenario]','red','No scenarios found')
       invisible()
     } else {
+      assertthat::assert_that(what %in% names( self$get_data() ),
+                              msg = paste(what, "not found in scenario projections?"))
       # Get unique number of data values. Surely there must be an easier val
       vals <- self$get_data()[what] |> stars:::pull.stars() |> as.vector() |> unique()
       vals <- base::length(stats::na.omit(vals))

diff --git a/R/project.R b/R/project.R
@@ -14,6 +14,7 @@ NULL
 #' @details In the background the function \code{x$project()} for the respective
 #' model object is called, where \code{x} is fitted model object. For specifics
 #' on the constraints, see the relevant \code{constrain} functions, respectively:
+#'
 #' * [`add_constraint()`] for generic wrapper to add any of the available constrains.
 #' * [`add_constraint_dispersal()`] for specifying dispersal constraint on the temporal projections at each step.
 #' * [`add_constraint_MigClim()`] Using the \pkg{MigClim} R-package to simulate dispersal in projections.
@@ -33,6 +34,11 @@ NULL
 #' computation of thresholds at every step based on the threshold used for the
 #' main model (threshold values are taken from there).
 #'
+#' It is also possible to make a complementary simulation with the \code{steps}
+#' package, which can be provided via [`simulate_population_steps()`] to the
+#' [`BiodiversityScenario-class`] object. Similar as with thresholds, estimates
+#' values will then be added to the outputs.
+#'
 #' Finally this function also allows temporal stabilization across prediction
 #' steps via enabling the parameter \code{stabilize} and checking the
 #' \code{stablize_method} argument. Stabilization can for instance be helpful in
@@ -175,7 +181,9 @@ methods::setMethod(
     if(inherits(baseline_threshold, 'SpatRaster')){
       baseline_threshold <- baseline_threshold[[grep(layer, names(baseline_threshold))]]
     }
+    # Optional constraints or simulations if specified
     scenario_constraints <- mod$get_constraints()
+    scenario_simulations <- mod$get_simulation()
 
     # Create a template for use
     template <- emptyraster( mod$get_predictors()$get_data() )
@@ -208,13 +216,15 @@ methods::setMethod(
     assertthat::assert_that(nrow(df)>0,
                             utils::hasName(df,'x'), utils::hasName(df,'y'), utils::hasName(df,'time'),
                             msg = "Error: Projection data and training data are not of equal size and format!")
+
     df <- subset(df, select = c("x", "y", "time", mod_pred_names) )
     df$time <- to_POSIXct( df$time )
     # Convert all units classes to numeric or character to avoid problems
     df <- units::drop_units(df)
 
     # ------------------ #
-    if(getOption('ibis.setupmessages')) myLog('[Scenario]','green','Starting suitability projections for ', length(unique(df$time)), ' timesteps.')
+    if(getOption('ibis.setupmessages')) myLog('[Scenario]','green','Starting suitability projections for ',
+                                              length(unique(df$time)), ' timesteps.')
 
     # Now for each unique element, loop and project in order
     proj <- terra::rast()
@@ -496,6 +506,23 @@ methods::setMethod(
       }
     }
 
+    # Optional simulation steps
+    if(!is.Waiver(scenario_simulations)){
+      if("steps" %in% scenario_simulations$simulation$method ){
+        pops <- .simulate_steps(proj, scenario_simulations)
+        if(!is.null(pops)){
+          assertthat::assert_that(is.Raster(pops))
+          pops <- stars::st_as_stars(pops,
+                                     crs = sf::st_crs(new_crs)
+          ); names(pops) <- 'population'
+        }
+      }
+    } else {
+      pops <- NULL
+    }
+
+    # --------------------------------- #
+    # # # # # # # # # # # # # # # # # # #
     # Finally convert to stars and rename
     proj <- stars::st_as_stars(proj,
                                crs = sf::st_crs(new_crs)
@@ -510,8 +537,24 @@ methods::setMethod(
       if(all(!stars::st_get_dimension_values(proj, 3) != stars::st_get_dimension_values(proj_thresh, 3 ))){
         proj_thresh <- stars::st_set_dimensions(proj_thresh, 3, values = stars::st_get_dimension_values(proj, 3))
       }
-      proj <- stars:::c.stars(proj, proj_thresh)
+      # Try
+      new <- try({ c(proj, proj_thresh) },silent = TRUE)
+      if(inherits(new, "try-error")){
+        # Replace dimensions and check again.
+        # This error likely originates from some part of the modelling routing
+        # reprojecting input layers
+        stars::st_dimensions(proj_thresh) <- stars::st_dimensions(proj)
+        new <- try({ c(proj, proj_thresh) },silent = TRUE)
+      }
+      proj <- new; rm(new)
+    }
+    if(!is.null(pops)){
+      stars::st_dimensions(pops) <- stars::st_dimensions(proj)
+
+      proj <- try({ c(proj, pops) },silent = TRUE)
     }
+    # # # # # # # # # # # # # # # # # # #
+    # --------------------------------- #
 
     # Return output by adding it to the scenario object
     bdproto(NULL, mod,