R/projection_changes.R

In kuenm2: Detailed Development of Ecological Niche Models

#' Compute changes of suitable areas between scenarios
#'
#' @description
#' This function performs map algebra operations to represent how and where
#' suitable areas change compared to the scenario in which the model was trained.
#' Changes are identified as loss (contraction), gain (expansion) and stability.
#' If multiple climate models (GCM) are used, it calculates the level of
#' agreement among them for each emission scenario.
#'
#' @usage
#' projection_changes(model_projections, reference_id = 1, consensus = "median",
#'                    include_id = NULL, user_threshold = NULL, by_gcm = TRUE,
#'                    by_change = TRUE, general_summary = TRUE,
#'                    force_resample = TRUE, write_results = TRUE,
#'                    output_dir = NULL, overwrite = FALSE,
#'                    write_bin_models = FALSE, return_raster = FALSE)
#'
#' @param model_projections a `model_projections` object generated by the
#' [project_selected()] function. This object contains the file
#' paths to the raster projection results and the thresholds used for binarization
#' of predictions.
#' @param reference_id (numeric) the reference ID for the projections
#' corresponding to the current time in `model_projections`. Default is 1.
#' See the details section for further information.
#' @param consensus (character) the consensus measure to use for calculating
#' changes. Available options are 'mean', 'median', 'range', and 'stdev'
#' (standard deviation). Default is 'median'.
#' @param include_id (numeric) a vector containing the reference IDs to include
#' when computing changes. Default is `NULL`, meaning all projections will be
#' included. See the details section for further information.
#' @param user_threshold (numeric) an optional threshold for binarizing the
#' predictions. Default is `NULL`, meaning the function will apply the
#' thresholds stored in model_projections, which were calculated earlier
#' using the omission rate from `calibration()`.
#' @param by_gcm (logical) whether to compute changes across GCMs.
#' Default is TRUE.
#' @param by_change (logical) whether to compute results separately for each
#' change, identifying areas of gain, loss, and stability for each GCM.
#' Default is TRUE.
#' @param general_summary (logical) whether to generate a general summary,
#' mapping how many GCMs project gain, loss, and stability for each scenario.
#' Default is TRUE.
#' @param force_resample (logical) whether to force the projection rasters to
#' have the same extent and resolution as the raster corresponding to the
#' `reference_id`, which represents the current projections. Default is TRUE.
#' @param write_results (logical) whether to write the raster files containing
#' the computed changes to the disk. Default is TRUE.
#' @param output_dir (character) the directory path where the resulting raster
#' files containing the computed changes will be saved. Only relevant if
#' `write_results = TRUE`.
#' @param overwrite (logical) whether to overwrite SpatRaster if they already
#' exist. Only applicable if `write_results` is set to TRUE. Default is FALSE.
#' @param write_bin_models (logical) whether to write the binarized models for
#' each GCM to the disk. Default is FALSE.
#' @param return_raster (logical) whether to return a list containing all the
#' SpatRasters with the computed changes. Default is FALSE, meaning the function
#' will return a NULL object. Setting this argument to TRUE while using multiple
#' GCMs at a large extent and fine resolution may overload the RAM.
#'
#' @details
#' When projecting a niche model to different temporal scenarios (past or
#' future), species’ areas can be classified into three categories relative to
#' the current baseline: **gain**, **loss** and **stability**. The
#' interpretation of these categories depends on the temporal direction of the projection.
#' **When projecting to future scenarios**:
#' - *Gain*: Areas that are currently unsuitable become suitable in the future.
#' - *Loss*: Areas that are currently suitable become unsuitable in the future.
#' - *Stability*: Areas that retain their current classification in the future,
#' whether suitable or unsuitable.
#'
#' **When projecting to past scenarios**:
#' - *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.
#' - *Stability*: Areas that retain their past classification in the present,
#' whether suitable or unsuitable.
#'
#' The reference scenario (current conditions) can be accessed in the paths
#' element of the model_projections object (model_projections$path). The ID will
#' differ from 1 only if there is more than one projection for the current
#' conditions.
#'
#' Specific projections can be included or excluded from the analysis using the
#' `include_id` argument. For example, setting 'include_id = c(3, 5, 7)' will
#' compute changes only for scenarios 3, 5, and 7. Conversely, setting
#' 'include_id = -c(3, 5, 7)' will exclude scenarios 3, 5, and 7 from the
#' analysis.
#'
#' @return
#' A `changes_projections` object.
#' If return_raster = TRUE, the function returns a list containing the
#' SpatRasters with the computed changes. The list includes the following
#' elements:
#'  - Binarized: binarized models for each GCM.
#'  - Results_by_gcm: computed changes for each GCM.
#'  - Results_by_change: a list where each SpatRaster represents a specific
#'  change.
#'  - Summary_changes: A general summary that indicates how many GCMs project
#'  gain, loss, and stability for each scenario
#'  - root_directory: the path to the directory where the results were saved if
#'  write_results was set to TRUE
#'
#'  If return_raster = FALSE, the function returns a NULL object.
#'
#' @export
#'
#' @importFrom
#' terra rast res ext resample writeRaster unique nlyr mean app
#'
#' @seealso
#' [organize_future_worldclim()], [prepare_projection()], [project_selected()]
#'
#' @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_raw3")
#' 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_raw3")
#'
#' ## 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/maxnet1")
#' 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)
#'
#' terra::plot(changes$Binarized)  # SpatRaster with the binarized predictions
#' terra::plot(changes$Results_by_gcm)  # SpatRaster with changes by GCM
#' changes$Results_by_change  # List of SpatRaster(s) by changes with GCM agreement
#' terra::plot(changes$Results_by_change$`Future_2081-2100_ssp585`)  # an example of the previous
#' terra::plot(changes$Summary_changes)  # SpatRaster with a general summary

 projection_changes <- function(model_projections,
                               reference_id = 1,
                               consensus = "median",
                               include_id = NULL,
                               user_threshold = NULL,
                               by_gcm = TRUE,
                               by_change = TRUE,
                               general_summary = TRUE,
                               force_resample = TRUE,
                               write_results = TRUE,
                               output_dir = NULL,
                               overwrite = FALSE,
                               write_bin_models = FALSE,
                               return_raster = FALSE) {

  #Check data
  if (missing(model_projections)) {
    stop("Argument 'model_projections' must be defined.")
  }
  if (!inherits(model_projections, "model_projections")) {
    stop("Argument 'model_projections' must be a 'model_projections' object.")
  }
  if (!inherits(reference_id, "numeric")) {
    stop("Argument 'reference_id' must be 'numeric'.")
  }
  if (!reference_id %in% model_projections$paths$id) {
    stop("'reference_id' is not present in 'model_projections$path'.")
  }
  #Get time of reference id
  time_reference <- model_projections$paths$Time[model_projections$paths$id ==
                                                   reference_id]
  if (time_reference != "Present") {
    warning("Reference scenario is ", time_reference, ", not the present time.\n",
            "To set the present time as reference scenario, check the IDS in model_projections$path")
  }
  if (!inherits(consensus, "character")) {
    stop("Argument 'consensus' must be a 'character'.")
  }
  if (length(consensus) > 1) {
    stop("Argument 'consensus' must be a unique value.",
         "\nAvailable options are: 'median', 'range', 'mean' or 'stdev'.")
  }
  consensus_out <- setdiff(consensus, c("median", "range", "mean", "stdev"))
  if (length(consensus_out) > 0) {
    stop("Invalid 'consensus' provided.",
         "\nAvailable options are: 'median', 'range', 'mean' or 'stdev'.")
  }
  if (!is.null(include_id) & !inherits(include_id, "numeric")) {
    stop("Argument 'include_id' must be NULL or 'numeric'.")
  }
  if (!is.null(user_threshold) & !inherits(user_threshold, "numeric")) {
    stop("Argument 'user_threshold' must be NULL or 'numeric'.")
  }

  if (!return_raster & !write_results) {
    stop("'return_raster' and 'write_results' cannot both be set to FALSE.
            Changing 'return_raster' or/and 'write_results' to TRUE.")
    return_raster <- TRUE
  }
  if (write_results & is.null(output_dir)) {
    stop("If 'write_results' = TRUE, 'output_dir' must be specified.")
  }
  if (write_results & !inherits(output_dir, "character")) {
    stop("Argument 'output_dir' must be a 'character'.")
  }

  if (!inherits(force_resample, "logical")) {
    stop("Argument 'force_resample' must be 'logical' (TRUE or FALSE)")
  }

  #Create directory to save results, if necessary
  if (write_results) {
    if (is.null(output_dir)) {
      stop("If write_results = TRUE, you must define output_dir")
    }
    output_dir <- file.path(output_dir, "Projection_changes/")
    if (!file.exists(output_dir)) {
    dir.create(output_dir, showWarnings = FALSE, recursive = TRUE)
    }
  } else {
    output_dir <- NULL
  }

  #Extract threshold
  if (is.null(user_threshold)) {
    thr <- model_projections[["thresholds"]][["consensus"]][[consensus]]
  } else {
    thr <- user_threshold
  }

  #Remove threshold from projection paths
  model_projections <- model_projections[["paths"]]


  #Get only specified scenarios
  if (!is.null(include_id)) {
    to_include <- include_id[include_id > 0]
    if (length(to_include) > 0) {
      model_projections <- model_projections[model_projections$id %in%
                                               to_include, ]
    }
    to_exclude <- abs(include_id[include_id < 0])
    if (length(to_exclude) > 0) {
      model_projections <- model_projections[!(model_projections$id %in%
                                               to_exclude), ]
    }
  }


  #Get raster of reference
  dir_reference <- model_projections$output_path[which(model_projections$id ==
                                                        reference_id)]
  file_reference <- file.path(dir_reference, "General_consensus.tif")
  r <- terra::rast(file_reference)[[consensus]]
  names(r) <- model_projections$Time[which(model_projections$id ==
                                                   reference_id)]


  ####Binarize models ####
  r_bin <- binarize_values(x = r, v = thr, new_value = 2)
  #Set levels
  levels(r_bin) <- data.frame(id = c(0, 2),
                              Result = c("Unsuitable", "Suitable"))
  names(r_bin) <- "Present"

  #plot(r_bin)
  #Looping throughout projections to binarize
  pp <- model_projections[model_projections$Time %in% c("Past", "Future"), ]

  #Binarize
  proj_bin <- terra::rast(lapply(1:nrow(pp), function(i) {
    #Get scenario projection
    d_i <- pp[i, ]
    dir_change <- d_i$output_path

    #Get time
    proj_time <- d_i$Time
    file_change <- file.path(dir_change, "General_consensus.tif")
    r_change <- terra::rast(file_change)[[consensus]]

    #Force resample if necessary
    if (force_resample) {
      if (any(terra::res(r) != (terra::res(r_change)) | terra::ext(r) !=
              terra::ext(r_change))) {
        r_change <- terra::resample(r_change, r_bin, method = "average")
      }
    }

    r_change <- binarize_values(x = r_change, v = thr, new_value = 1)
    levels(r_change) <- data.frame(id = c(0, 1),
                                   Result = c("Unsuitable", "Suitable"))
    return(r_change)
  }))
  #Rename binarizes scenarions
  names(proj_bin) <- paste(pp$Time, pp$Period, pp$Scenario, pp$GCM, sep = "_")
  names(proj_bin) <- gsub("_NA_", "_", names(proj_bin))

  if (write_bin_models) {
    terra::writeRaster(x = c(r_bin, proj_bin),
                      filename = file.path(output_dir, "Binarized.tif"),
                      overwrite = overwrite)
  }


  #Get single scenarios by Time and period
  sc <- unique(pp[, c("Time", "Period", "Scenario")])

  #Table to set levels in raster
  cls_future <- data.frame(id = c(1, 2, 3, 0),
                           Result = c("Gain", "Loss", "Stable suitable",
                                      "Stable unsuitable"))
  cls_past <- data.frame(id = c(1, 2, 3, 0),
                         Result = c("Loss", "Gain", "Stable suitable",
                                    "Stable unsuitable"))


  ####Identify changes by scenario####
  if (by_gcm | by_change) {
    # #Change values of r to compute gain and loss
    # r2 <- binarize_values(r, v = 1, new_value = 2)
    # r2[r2 == 0] <- 1
    # #plot(r2)

  res_by_gcm <- lapply(proj_bin, function(i) {
      #Get levels (past or future)
      if(grepl("Future", names(i))){
        cls <- cls_future
      } else if(grepl("Past", names(i))){
        cls <- cls_past
      }

      #Compute changes
      r_result <- r_bin + i
      #Get legend
      levels(r_result) <- cls
      #plot(r_result)
      return(r_result)
    })

    #Rename res
    names(res_by_gcm) <- names(proj_bin)

    #Rasterize res
    res_by_gcm <- terra::rast(res_by_gcm)

    if (write_results & by_gcm) {
      terra::writeRaster(x = res_by_gcm,
                         filename = file.path(output_dir, "Changes_by_GCM.tif"),
                         overwrite = overwrite)
    }

  } #End of by_gcm


  ####Results by change####
  if (by_change) {
    res_by_change <- lapply(1:nrow(sc), function(i) {
      sc_i <- sc[i, ]
      scenario_i <- paste(sc_i$Time, sc_i$Period, sc_i$Scenario, sep = "_")
      scenario_i <- gsub("_NA", "_", scenario_i)



      #Get levels (past or future)
      if(grepl("Future", sc_i$Time)){
        cls <- cls_future
      } else if(grepl("Past", sc_i$Time)){
        cls <- cls_past
      }

      #Subset results
      res_i <- res_by_gcm[[grep(scenario_i, names(res_by_gcm))]]
      #Looping throught changes
      #Get available changes
      out <- unique(unlist(sapply(res_i, terra::unique, simplify = T)))
      sc_out <- lapply(out, function(x) {
        #Get change value id
        out_id <- cls$id[cls$Result == x]
        res_i_bin <- (res_i == out_id) * 1

        #Sum results
        res_i_sum <- sum(res_i_bin) ####Import sum from terra
        #Levels
        l_sum <- data.frame(
          id = terra::unique(res_i_sum)$sum,
          Result = paste0(x, " in ", terra::unique(res_i_sum)$sum, " GCMs")
        )
        levels(res_i_sum) <- l_sum
        return(res_i_sum)
      })
      names(sc_out) <- out
      return(terra::rast(sc_out))
    })
    #Set names by scenarion
    names(res_by_change) <- gsub("_NA", "",
                                 paste(sc$Time, sc$Period, sc$Scenario, sep = "_"))

    #Save results
    if (write_results) {
      dir.create(file.path(output_dir, "Results_by_change"))
      sapply(names(res_by_change), function(z) {
        terra::writeRaster(x = res_by_change[[z]],
                           filename = file.path(output_dir, "Results_by_change",
                                                paste0(z, ".tif")),
                           overwrite = overwrite)
      })
    }
  } else {
    res_by_change <- NULL
  } #End of by change


  ####General summary####
  if (general_summary) {
    res_summary <- lapply(1:nrow(sc), function(i) {
      sc_i <- sc[i, ]
      scenario_i <- paste(sc_i$Time, sc_i$Period, sc_i$Scenario, sep = "_")
      scenario_i <- gsub("_NA", "_", scenario_i)

      #Get levels (past or future)
      if(grepl("Future", sc_i$Time)){
        cls <- cls_future
      } else if(grepl("Past", sc_i$Time)){
        cls <- cls_past
      }

      #Subset results
      res_i <- proj_bin[[grep(scenario_i, names(proj_bin))]]
      #Presence value in present (number of gcms + 1)
      n_gcms <- terra::nlyr(res_i) + 1
      #Change values of present
      r_present <- (r_bin == 2) * n_gcms

      #Sum rasters to get results
      res_sum <- sum(c(r_present, res_i))

      # # preparing description table
      # vals <- sort(na.omit(unique(res_sum[])))
      # #Fix vals
      # if (length(vals) != max(vals)) {
      #   vals <- seq(0, max(vals), 1)
      #   vals <- vals[vals != 2]
      # }

      vals <- seq(0, (n_gcms*2 - 1), 1)

      #Add category
      if(grepl("Past", sc_i$Time)){
        gain <- ceiling(max(vals)/2)
        g_val <- c(gain, vals[vals > gain & vals != max(vals)])
        l_val <- vals[vals < gain & vals != 0]
        gains <- paste0("Gain in ", max(vals) - g_val, " GCMs")
        gains[gains == paste0("Gain in ", n_gcms - 1, " GCMs")] <- "Gain in all GCMs"
        losses <- paste0("Loss in ", l_val, " GCMs")
        losses[losses == paste0("Loss in ", n_gcms - 1, " GCMs")] <- "Loss in all GCMs"
        descriptions <- c("Stable unsuitable",
                          losses,
                          gains,
                          "Stable suitable")

      }

      if(grepl("Future", sc_i$Time)){
        loss <- ceiling(max(vals)/2)
        l_val <- c(loss, vals[vals > loss & vals != max(vals)])
        g_val <- vals[vals < loss & vals != 0]
        gains <- paste0("Gain in ", g_val, " GCMs")
        gains[gains == paste0("Gain in ", n_gcms - 1, " GCMs")] <- "Gain in all GCMs"
        losses <- paste0("Loss in ", max(vals) - l_val, " GCMs")
        losses[losses == paste0("Loss in ", n_gcms - 1, " GCMs")] <- "Loss in all GCMs"
        descriptions <- c("Stable unsuitable", gains,
                          losses, "Stable suitable")
      }

      res_table <- data.frame(Raster_value = vals, Description = descriptions)

      #Set levels
      levels(res_sum) <- res_table
      #plot(res_sum)
      return(res_sum)
    })
    #Set names by scenario
    res_summary <- terra::rast(res_summary) #Rasterize
    names(res_summary) <- gsub("_NA", "",
                               paste(sc$Time, sc$Period, sc$Scenario, sep = "_"))

    if (write_results) {
      terra::writeRaster(x = res_summary,
                         filename = file.path(output_dir, "Changes_summary.tif"),
                         overwrite = overwrite)
    }

  } else {
    res_summary <- NULL
  } #End of general summary


  if (return_raster) {
    #Create final object
    res_final <- list(Binarized = c(r_bin, proj_bin),
                      Results_by_gcm = res_by_gcm,
                      Results_by_change = res_by_change,
                      Summary_changes = res_summary,
                      root_directory = output_dir)

  } else {
    res_final <- list(root_directory = output_dir)
  }
  class(res_final) <- "changes_projections"
  #Save changes_projections object (only root directory)
  if(write_results){
    res_to_write <- res_final["root_directory"]
    class(res_to_write) <- "changes_projections"
    saveRDS(res_to_write,
            file.path(output_dir, "changes_projections.rds"))
  }
  return(res_final)
} #End of function