R/binarize_changes.R

In kuenm2: Detailed Development of Ecological Niche Models

#' Binarize changes based on the agreement among GCMs
#'
#' @description
#' Highlights areas where a specified number of Global Climate Models (GCMs)
#' agree on a given outcome in the projected scenario. The function can identify
#' areas classified as suitable, unsuitable, stable-suitable, stable-unsuitable,
#' gained, or lost.
#'
#' @usage binarize_changes(changes_projections, outcome = "suitable", n_gcms)
#'
#' @param changes_projections an object of class `changes_projections`,
#' generated by `projection_changes()` or imported using `import_projections()`,
#' containing the `$Summary_changes` element.
#' @param outcome (character) The outcome to binarize. Available options are
#' "suitable", "unsuitable", "stable-suitable", "stable-unsuitable", "gain" or
#' "loss". Default is "suitable". See details.
#' @param n_gcms (numeric) The minimum number of GCMs that must agree on the
#' specified outcome for a cell to be included in that category.
#'
#' @details
#' The interpretation of the outcomes depends on the temporal direction of the
#' projection.
#' **When projecting to future scenarios**:
#' - *suitable*: Areas that remain suitable (stable-suitable) or become suitable
#' (gain) in the future.
#' - *unsuitable*: Areas that remain unsuitable (stable-unsuitable) or become
#' unsuitable (loss) in the future.
#' - *gain*: Areas that are currently unsuitable become suitable in the future.
#' - *loss*: Areas that are currently suitable become unsuitable in the future.
#' - *stable-suitable or stable-unsuitable*: Areas that retain their current
#' classification in the future, whether suitable or unsuitable.
#'
#' **When projecting to past scenarios**:
#' - *suitable*: Areas that remain suitable (stable-suitable) or become
#' unsuitable (loss) in the present.
#' - *unsuitable*: Areas that remain unsuitable (stable-unsuitable) or become
#' suitable (gain) in the present
#' - *gain*: Areas that were unsuitable in the past are now suitable in the
#' present.
#' - *loss*: Areas that were suitable in the past are now unsuitable in the
#' present.
#' - *stable-suitable or stable-unsuitable*: Areas that retain their current
#' classification in the future, whether suitable or unsuitable.
#'
#'
#' @returns
#' A `SpatRaster` or a list of `SpatRaster` objects (one per scenario) with the
#' binarized outcomes. For example, if `outcome = "suitable"` and `n_gcms = 3`,
#' cells with a value of 1 indicate areas where three or more GCMs agree that
#' the area is suitable for the species in that scenario.
#' @export
#'
#' @importFrom terra levels classify
#' @examples
#' # Step 1: Organize variables for current projection
#' ## Import current variables (used to fit models)
#' var <- terra::rast(system.file("extdata", "Current_variables.tif",
#'                                package = "kuenm2"))
#'
#' ## Create a folder in a temporary directory to copy the variables
#' out_dir_current <- file.path(tempdir(), "Current_raw_bin")
#' dir.create(out_dir_current, recursive = TRUE)
#'
#' ## Save current variables in temporary directory
#' terra::writeRaster(var, file.path(out_dir_current, "Variables.tif"))
#'
#'
#' # Step 2: Organize future climate variables (example with WorldClim)
#' ## Directory containing the downloaded future climate variables (example)
#' in_dir <- system.file("extdata", package = "kuenm2")
#'
#' ## Create a folder in a temporary directory to copy the future variables
#' out_dir_future <- file.path(tempdir(), "Future_raw_bin")
#'
#' ## Organize and rename the future climate data (structured by year and GCM)
#' ### 'SoilType' will be appended as a static variable in each scenario
#' organize_future_worldclim(input_dir = in_dir, output_dir = out_dir_future,
#'                           name_format = "bio_", static_variables = var$SoilType)
#'
#' # Step 3: Prepare data to run multiple projections
#' ## An example with maxnet models
#' ## Import example of fitted_models (output of fit_selected())
#' data(fitted_model_maxnet, package = "kuenm2")
#'
#' ## Prepare projection data using fitted models to check variables
#' pr <- prepare_projection(models = fitted_model_maxnet,
#'                          present_dir = out_dir_current,
#'                          future_dir = out_dir_future,
#'                          future_period = c("2081-2100"),
#'                          future_pscen = c("ssp585"),
#'                          future_gcm = c("ACCESS-CM2", "MIROC6"),
#'                          raster_pattern = ".tif*")
#'
#' # Step 4: Run multiple model projections
#' ## A folder to save projection results
#' out_dir <- file.path(tempdir(), "Projection_results/maxnet_bin")
#' dir.create(out_dir, recursive = TRUE)
#'
#' ## Project selected models to multiple scenarios
#' p <- project_selected(models = fitted_model_maxnet, projection_data = pr,
#'                       out_dir = out_dir)
#'
#' # Step 5: Identify areas of change in projections
#' ## Contraction, expansion and stability
#' changes <- projection_changes(model_projections = p, write_results = FALSE,
#'                               return_raster = TRUE)
#'
#' # Step 6: Binarize changes
#' future_suitable <- binarize_changes(changes_projections = changes,
#'                                     outcome = "suitable",
#'                                     n_gcms = 1)
#' terra::plot(future_suitable)

binarize_changes <- function(changes_projections,
                             outcome = "suitable",
                             n_gcms) {
  #Check data
  if (missing(changes_projections)) {
    stop("Argument 'changes_projections' must be defined.")
  }
  if (!inherits(changes_projections, "changes_projections")) {
    stop("Argument 'changes_projections' must be a 'changes_projections' object.")
  }

  if (!inherits(outcome, "character")) {
    stop("Argument 'outcome' must be a 'character'.")
  }

  if(length(outcome) > 1){
    stop("Argument 'outcome' must be a single character")
  }

  if(!(outcome %in% c("suitable", "unsuitable", "gain", "loss",
                        "stable-suitable", "stable-unsuitable"))){
    stop("Argument 'outcome' must be 'suitable', 'unsuitable', 'gain', 'loss',
'stable-suitable', or 'stable-unsuitable'")
  }

  if (!inherits(n_gcms, "numeric")) {
    stop("Argument 'n_gcms' must be a 'numeric' object.")
  }


  #Get Summary changes
  r_changes <- changes_projections$Summary_changes

  res_changes <- lapply(r_changes, function(i){
    #Extract levels
    l <- terra::levels(i)[[1]]
    colnames(l)[1] <- "value"
    l$event <- sub(" in.*", "", l[,2])
    l$n_gcms <- as.numeric(gsub("\\D+", "", l[,2]))
    l$n_gcms[l$n_gcms == "" |
               is.na(l$n_gcms)] <- max(as.numeric(l$n_gcms), na.rm = TRUE) + 1

    # Check the direction of projection
    if(grepl("Future", names(i))){
      d <- "Future"
    } else if(grepl("Past", names(i))){
      d <- "Past"
    }

    # Extract outcomes
    if(outcome == "suitable" && d == "Future"){
      v <- l[l$event %in% c("gain", "stable, suitable") &
               l$n_gcms >= n_gcms,] }
    if(outcome == "suitable" && d == "Past"){
      v <- l[l$event %in% c("loss", "stable, suitable") &
               l$n_gcms >= n_gcms,] }

    if(outcome == "unsuitable" && d == "Future"){
      v <- l[l$event %in% c("loss", "stable, unsuitable") &
               l$n_gcms >= n_gcms,] }
    if(outcome == "unsuitable" && d == "Past"){
      v <- l[l$event %in% c("gain", "stable, unsuitable") &
               l$n_gcms >= n_gcms,] }

    if(outcome == "gain"){
      v <- l[l$event == "gain" &
               l$n_gcms >= n_gcms,]
    }
    if(outcome == "loss"){
      v <- l[l$event == "loss" &
               l$n_gcms >= n_gcms,]
    }
    if(outcome == "stable-suitable"){
      v <- l[l$event == "stable, suitable" &
               l$n_gcms >= n_gcms,]
    }
    if(outcome == "stable-unsuitable"){
      v <- l[l$event == "stable, unsuitable" &
               l$n_gcms >= n_gcms,]
    }

    # Reclassify
    if(nrow(v) > 0){
      v$new_value <- 1
      rv <- terra::classify(i, v[, c("value", "new_value")], others = 0)
    } else {
      rv <- i * 0
    }

    names(rv) <- outcome
   return(rv)
  }) # End of lapply

  res <- rast(res_changes)
  names(res) <- names(r_changes)
  return(res)
} #End of function