R/projection_variability.R

In kuenm2: Detailed Development of Ecological Niche Models

#' Explores variance coming from distinct sources in model predictions
#'
#' @description
#' Calculates variance in model predictions, distinguishing
#' between the different sources of variation. Potential sources include
#' replicates, model parameterizations, and general circulation models (GCMs).
#'
#' @usage
#' projection_variability(model_projections, from_replicates = TRUE,
#'                        from_parameters = TRUE, from_gcms = TRUE,
#'                        consensus = "median", write_files = FALSE,
#'                        output_dir = NULL, return_rasters = TRUE,
#'                        progress_bar = FALSE, verbose = TRUE,
#'                        overwrite = FALSE)
#'
#' @param model_projections a `model_projections` object generated by the
#' \code{\link{project_selected}}() function. This object contains the file
#' paths to the raster projection results and the thresholds used for binarizing
#' the predictions.
#' @param from_replicates (logical) whether to compute the variance originating
#' from replicates.
#' @param from_parameters (logical) whether to compute the variance originating from
#' model parameterizations.
#' @param from_gcms (logical) whether to compute the variance originating from
#' general circulation models (GCMs)
#' @param consensus (character) (character) the consensus measure to use for
#' calculating changes. Available options are 'mean', 'median', 'range', and
#' 'stdev' (standard deviation). Default is 'median'.
#' @param write_files (logical) whether to write the raster files containing
#' the computed variance to the disk. Default is FALSE.
#' @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 return_rasters (logical) whether to return a list containing all the
#' SpatRasters with the computed changes. Default is TRUE. Setting this
#' argument to FALSE returns a NULL object.
#' @param progress_bar (logical) whether to display a progress bar during
#' processing. Default is TRUE.
#' @param verbose (logical) whether to display messages during processing.
#' Default is TRUE.
#' @param overwrite whether to overwrite SpatRaster if they already exists.
#' Only applicable if `write_files` is set to TRUE. Default is FALSE.
#'
#' @return
#' An object of class `variability_projections`. If `return_rasters = TRUE`,
#' the function returns a list containing the SpatRasters with the computed
#' variances, categorized by replicate, model, and GCMs. If `write_files = TRUE`,
#' it also returns the directory path where the computed rasters were saved to
#' disk, and the object can then be used to import these files later with the
#' `import_results()` function. If both `return_rasters = FALSE` and
#' `write_files = FALSE`, the function returns `NULL`
#'
#' @export
#'
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @importFrom terra rast nlyr mean app writeRaster
#'
#' @seealso
#' [organize_future_worldclim()], [prepare_projection()], [project_selected()],
#' [import_results()]
#'
#' @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_raw5")
#' 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_raw5")
#'
#' ## 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 = "2041-2060",
#'                          future_pscen = "ssp126",
#'                          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/maxnet3")
#' 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: Compute variance from distinct sources
#' v <- projection_variability(model_projections = p, from_replicates = FALSE)
#'
#' #terra::plot(v$Present$from_replicates)  # Variance from replicates, present projection
#' terra::plot(v$Present$from_parameters)  # From models with distinct parameters
#' #terra::plot(v$`Future_2041-2060_ssp126`$from_replicates)  # From replicates in future projection
#' terra::plot(v$`Future_2041-2060_ssp126`$from_parameters)  # From models
#' terra::plot(v$`Future_2041-2060_ssp126`$from_gcms)  # From GCMs


projection_variability <- function(model_projections,
                                   from_replicates = TRUE,
                                   from_parameters = TRUE,
                                   from_gcms = TRUE,
                                   consensus = "median",  # MAKE IT WORK WHEN CONSENSUS IS FULL MODEL
                                   write_files = FALSE,
                                   output_dir = NULL,
                                   return_rasters = TRUE,
                                   progress_bar = FALSE,
                                   verbose = TRUE,
                                   overwrite = 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(consensus, "character")) {
    stop("Argument 'consensus' must be a 'character'.")
  }

  if (length(consensus) > 1) {
    stop("Argument 'consensus' must be a single value.",
         "\nOptions are: 'median' or 'mean'.")
  }
  consensus_out <- setdiff(consensus, c("median", "mean"))
  if (length(consensus_out) > 0) {
    stop("Invalid 'consensus' provided.",
         "\nOptions are: 'median' or 'mean'.")
  }

  if (write_files & is.null(output_dir)) {
    stop("If 'write_files' = TRUE, 'output_dir' must be specified.")
  }
  if (write_files & !inherits(output_dir, "character")) {
    stop("Argument 'output_dir' must be a 'character'.")
  }

  if (!return_rasters & !write_files) {
    stop("'return_rasters' and 'write_files' cannot be both = FALSE.",
         "\nSetting 'return_rasters' to TRUE.")
    return_rasters <- TRUE
  }


  #Create directory if necessary
  if (write_files) {
    out_dir <- file.path(output_dir, "variance")
    dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
  } else {
    out_dir <- NULL
  }

  #### Get data ####
  d <- model_projections[["paths"]]


  # Get unique combinations of Time, Period, and Scenario
  uc <- unique(d[, c("Time", "Period", "Scenario")])

  #Show progress bar?
  if (progress_bar) {
    pb <- utils::txtProgressBar(0, nrow(uc), style = 3)
    progress <- function(n) utils::setTxtProgressBar(pb, n)
  }

  ####Iteration over combinations####
  res <- lapply(1:nrow(uc), function(z) {
    time <- uc$Time[z]
    period <- uc$Period[z]
    scenario <- uc$Scenario[z]

    # Filter dataframe to current combination
    d_p <- d[(is.na(d$Period) | d$Period == period) &
               (is.na(d$Scenario) | d$Scenario == scenario) &
               (is.na(d$Time) | d$Time == time), ]

    # Get all paths
    paths <- d_p$output_path

    #### By replicate ####
    if (from_replicates) {
      if (verbose) {
        message("\nCalculating variability from distinct replicates: scenario ",
                z, " of ", nrow(uc))
      }
      #Variance of means

      #### By replicates ####
      # Get variance of replicates in each gcm, than get the average across gcms
      var_rep <- terra::rast(lapply(paths, var_models_rep_by_gcm))
      var_rep <- terra::mean(var_rep)
      names(var_rep) <- "from_replicates"
    } else {#End of from_replicates
      var_rep <- NULL
    }

    #### By Model ####
    if (from_parameters) {
      if (verbose) {
        message("Calculating variability from distinct models: scenario ",
                z, " of ", nrow(uc))
      }

      # Get variance of models in each gcm, than get the average
      if (names(model_projections$thresholds[[1]])[1] == "Full_model") {
        var_model <- terra::rast(lapply(paths, var_models_model_by_gcm,
                                        "Full_model"))
      } else {
        var_model <- terra::rast(lapply(paths, var_models_model_by_gcm,
                                        consensus))
      }

      var_model <- terra::mean(var_model)
      names(var_model) <- "from_parameters"
    } else { #End of by model
      var_model <-  NULL
    }


    ####By GCM####
    if (from_gcms & period != "Present") {
      if (verbose) {
        message("Calculating variability from distinct GCMs: scenario ",
                z, " of ", nrow(uc))
      }

      var_gcm <- var_models_across_gcm(paths = paths, consensus = consensus)
      names(var_gcm) <- "from_gcms"
    } else {
      var_gcm <- NULL
    }#End of from_gcms

    all_var <- terra::rast(c("from_replicates" = var_rep, "from_parameters" = var_model,
                             "from_gcms" = var_gcm))


    #Write results
    if (write_files) {
      #Name of raster
      nr <- gsub("_NA", "", paste(time, period, scenario, sep = "_"))
      if (nr == "Present_Present_Present") {
        nr <- "Present"
      }
      terra::writeRaster(all_var,
                         filename = file.path(out_dir, paste0(nr, ".tif")),
                         overwrite = overwrite)
    }

    #Progress bar
    if (progress_bar) {
      utils::setTxtProgressBar(pb, z)
    }

    if (return_rasters) {
      return(all_var)
    } else {
      return(invisible(NULL))
    }
  }) #End of res

  if (return_rasters) {
    names(res) <- gsub("Present_Present_Present", "Present",
                       gsub("_NA", "",
                       paste(uc$Time, uc$Period, uc$Scenario, sep = "_")))
    #Append directory
    res[["root_directory"]] <- out_dir

  } else {
    res <- list("root_directory" = out_dir)
  }
  class(res) <- "variability_projections"

  #Save variability_projections object (only root directory)
  if(write_files){
    res_to_write <- res["root_directory"]
    class(res_to_write) <- "variability_projections"
    saveRDS(res_to_write,
            file.path(out_dir, "variability_projections.rds"))
  }

  return(res)
} #End of function