R/colors_for_changes.R

In kuenm2: Detailed Development of Ecological Niche Models

#' Set Colors for Change Maps
#'
#' @description
#' This functions sets the color tables associated with the `SpatRaster` object
#' resulting from `projection_changes()`. Color tables are used to associate specific colors with raster values when using `plot()`. This function defines custom colors for areas of gain, loss, and stability across scenarios.
#'
#' @param changes_projections an object of class `changes_projections`,
#' generated by `projection_changes()` or imported using `import_projections()`,
#' containing the `$Summary_changes` element.
#' @param gain_color (character) color used to define the palette for
#' representing gains. Default is "#009E73" (teal green).
#' @param loss_color (character) color used to define the palette for
#' representing losses. Default is "#D55E00" (orange-red).
#' @param stable_suitable (character) color used for representing areas that
#' remain suitable across all scenarios. Default is "#0072B2" (oxford blue).
#' @param stable_unsuitable (character) color used for representing areas that
#' remain unsuitable across all scenarios. Default is "grey".
#' @param max_alpha (numeric) opacity value (from 0 to 1) for areas where all
#' GCMs agree on the change (gain, loss, or stability). Default is 1.
#' @param min_alpha (numeric) opacity value (from 0 to 1) for areas where only
#' one GCM predicts a given change.  Default is 0.25
#'
#' @return
#' An object of class `changes_projections` with the same structure and
#' `SpatRaster`s as the input `changes_projections`, but with color tables
#' embedded in the `SpatRaster`s. These colors are used automatically when
#' visualizing the data with `plot()`.
#' @export
#' @importFrom terra coltab rast levels
#' @importFrom grDevices col2rgb rgb
#' @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_color_example")
#' 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_color_example")
#'
#' ## 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("ssp126", "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_color_example")
#' 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, by_gcm = TRUE,
#'                               by_change = TRUE, write_results = FALSE,
#'                               return_raster = TRUE)
#'
#' #Step 6: Set Colors for Change Maps
#' changes_with_colors <- colors_for_changes(changes_projections = changes)
#' terra::plot(changes_with_colors$Summary_changes)

colors_for_changes <- function(changes_projections, gain_color = "#009E73",
                               loss_color = "#D55E00",
                               stable_suitable = "#0072B2",
                               stable_unsuitable = "grey",
                               max_alpha = 1, min_alpha = 0.25){

  #### Check data ####
  if(!inherits(changes_projections, "changes_projections")){
    stop("Argument 'changes_projections' must be an object of class 'changes_projections'")
  }

  if (!inherits(gain_color, "character") || length(gain_color) > 1) {
    stop("Argument 'gain_color' must be a single 'character' value.")
  }
  if (!inherits(loss_color, "character") || length(loss_color) > 1) {
    stop("Argument 'loss_color' must be a single 'character' value.")
  }
  if (!inherits(stable_suitable, "character") || length(stable_suitable) > 1) {
    stop("Argument 'stable_suitable' must be a single 'character' value.")
  }
  if (!inherits(stable_unsuitable, "character") || length(stable_unsuitable) > 1) {
    stop("Argument 'stable_unsuitable' must be a single 'character' value.")
  }
  if(!inherits(max_alpha, "numeric") || length(max_alpha) > 1 ||
     min(max_alpha) < 0 || max(max_alpha) > 1){
    stop("Argument 'max_alpha' must be a single numeric value between 0 and 1")
  }
  if(!inherits(min_alpha, "numeric") || length(min_alpha) > 1 ||
     min(min_alpha) < 0 || max(min_alpha) > 1){
    stop("Argument 'min_alpha' must be a single numeric value between 0 and 1")
  }
  if(min_alpha >= max_alpha){
    stop("Argument 'min_alpha' can't be equal or higher than 'max_alpha'")
  }

  #Create list to save results
  r <- list()

  if("Binarized" %in% names(changes_projections)){
    r_bin <- changes_projections$Binarized
    r_bin <- lapply(r_bin, function(i){
    #Set colors
    terra::coltab(i) <- data.frame(value = 0:1,
                                       col = c(stable_unsuitable,
                                               stable_suitable))
    return(i)
    })
    r[["Binarized"]] <- terra::rast(r_bin)
  }


  if("Results_by_gcm" %in% names(changes_projections)){
    r_by_gcm <- changes_projections$Results_by_gcm
    r_by_gcm <- lapply(r_by_gcm, function(i){
      #Set colors
      terra::coltab(i) <- data.frame(value = 0:3,
                                     col = c(stable_unsuitable, gain_color,
                                             loss_color, stable_suitable))
      return(i)
    })
    r[["Results_by_gcm"]] <- terra::rast(r_by_gcm)
  }

    if("Results_by_change" %in% names(changes_projections)){
      r_by_change <- changes_projections$Results_by_change
      r_by_change <- lapply(r_by_change, function(x){
        r_changes <- lapply(names(x), function(i){
          #Setcolor
          color_i <- if(i == "Stable unsuitable"){
            stable_unsuitable} else if (i == "Stable suitable"){
              stable_suitable} else if (i == "Gain"){
                gain_color} else if (i == "Loss"){
                  loss_color}

      set_colors_by_change(x_i = x[[i]], change = i, color = color_i,
                           min_alpha = min_alpha, max_alpha = max_alpha)
          })
      names(r_changes) <- names(x)
      return(terra::rast(r_changes))
      })
      r[["Results_by_change"]] <- r_by_change
    }


  if("Summary_changes" %in% names(changes_projections)){
    #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

      #Set colors
      colnames(l)[2] <- "Description"
      l$color <- NA
      l$color[grepl("Gain", l$Description)] <- gain_color
      l$color[grepl("Loss", l$Description)] <- loss_color
      l$color[grepl(" suitable", l$Description)] <- stable_suitable
      l$color[grepl(" unsuitable", l$Description)] <- stable_unsuitable

      #Set alphas
      a <- data.frame(n_gcms = 1:max(l$n_gcms),
                      alpha = seq(from = min_alpha, to = max_alpha,
                                  length.out = max(as.numeric(l$n_gcms))))
      l <- merge(l, a)

      #Create vector of colors in rgb based on alpha value
      l$rgb <- set_rgb_adjusted(l)

      #Set col to raster
      terra::coltab(i) <- data.frame(value = l$value, col = l$rgb)
      return(i)
    })

    res <- terra::rast(res_changes)
    names(res) <- names(changes_projections$Summary_changes)

    r[["Summary_changes"]] <- res
  }

  # Add root directory
  r[["root_directory"]] <- changes_projections$root_directory

  class(r) <- "changes_projections"
  return(r)
}