ecorisk.R

# written by Fabian Stenzel, based on work by Sebastian Ostberg
# 2022-2023 - stenzel@pik-potsdam.de

################# EcoRisk calc functions  ###################

#' Wrapper for calculating the ecosystem change metric EcoRisk
#'
#' Function to read in data for ecorisk, and call the calculation function once,
#' if overtime is FALSE, or for each timeslice of length window years, if
#' overtime is TRUE
#'
#' @param path_ref folder of reference run
#' @param path_scen folder of scenario run
#' @param read_saved_data whether to read in previously saved data
#'        (default: FALSE)
#' @param save_data file to save read in data to (default NULL)
#' @param save_ecorisk file to save EcoRisk data to (default NULL)
#' @param nitrogen include nitrogen outputs for pools and fluxes into EcoRisk
#'        calculation (default FALSE)
#' @param weighting apply "old" (Ostberg-like), "new", or "equal" weighting of
#'        vegetation_structure_change weights (default "equal")
#' @param time_span_reference vector of years to use as scenario period
#' @param time_span_scenario vector of years to use as scenario period
#' @param dimensions_only_local flag whether to use only local change component
#'        for water/carbon/nitrogen fluxes and pools, or use an average of
#'        local change, global change and ecosystem balance (default FALSE)
#' @param overtime logical: calculate ecorisk as time-series? (default: FALSE)
#' @param window integer, number of years for window length (default: 30)
#' @param debug_mode write out all nitrogen state variables (default FALSE)
#' @param replace_input_file_names list with alternative names for output
#'        identifiers to replace the ones in inst/ext_files/metric_files.yml.
#'        e.g. list(npp="mnpp") would replace the expected output for npp with
#'        mnpp followed by the automatically detected file extension (.bin.json)
#' @param suppress_warnings suppress warnings - default: TRUE
#' @param external_variability use externally supplied variability for the
#'        reference period? experimental! (default: FALSE)
#' @param c2vr external variability array
#'
#' @return list data object containing arrays of ecorisk_total,
#'         vegetation_structure_change, local_change, global_importance,
#'         ecosystem_balance, carbon_stocks, carbon_fluxes, water_fluxes
#'         (+ nitrogen_stocks and nitrogen_fluxes)
#'
#' @examples
#' \dontrun{
#' ecorisk_wrapper(
#'   path_ref = pnv_folder,
#'   path_scen = run_folder,
#'   read_saved_data = FALSE,
#'   nitrogen = TRUE,
#'   save_data = NULL,
#'   save_ecorisk = NULL,
#'   time_span_reference = c(1550:1579),
#'   time_span_scenario = c(1987:2016)
#' )
#' }
#'
#' @md
#' @export
ecorisk_wrapper <- function(path_ref,
                            path_scen,
                            read_saved_data = FALSE,
                            save_data = NULL,
                            save_ecorisk = NULL,
                            nitrogen = TRUE,
                            weighting = "equal",
                            time_span_reference,
                            time_span_scenario,
                            dimensions_only_local = FALSE,
                            overtime = FALSE,
                            window = 30,
                            debug_mode = FALSE,
                            replace_input_file_names = NULL,
                            external_variability = FALSE,
                            c2vr = NULL,
                            suppress_warnings = TRUE) {
  # check timespan consistency
  nyears <- length(time_span_reference)
  nyears_scen <- length(time_span_scenario)

  # prepare slices for overtime calculation
  slices <- (nyears_scen - window + 1)
  window_half <- round(window / 2.0)
  slice_years <- time_span_scenario[
    (window_half + 1):(nyears_scen - window_half + 1)
  ]

  if ((!nyears == window) || nyears_scen < window) {
    stop(
      "Timespan in reference is not equal to window size (",
      window,
      "), or scenario timespan is smaller than window size."
    )
  }

  if (!read_saved_data) {
    # translate output names (from metric_files.yml) and
    # folders to files_scenarios/reference lists
    metric_files <- system.file(
      "extdata",
      "metric_files.yml",
      package = "biospheremetrics"
    ) %>%
      yaml::read_yaml()

    file_extension <- get_major_file_ext(paste0(path_scen))
    files_scenario <- list()
    files_reference <- list()

    for (output in names(metric_files$metric$ecorisk_nitrogen$output)) {
      # Iterate over all outputs
      if (is.null(replace_input_file_names[[output]])) {
        for (file in metric_files$file_name[[output]]) {
          full_file_path_lu <- paste0(path_scen, file, ".", file_extension)
          if (file.exists(full_file_path_lu)) {
            files_scenario[[output]] <- full_file_path_lu
          }
          full_file_path_pnv <- paste0(path_ref, file, ".", file_extension)
          if (file.exists(full_file_path_pnv)) {
            files_reference[[output]] <- full_file_path_pnv
          }
        }
        if (is.null(files_scenario[[output]])) {
          stop(
            "None of the default file names for ", output,
            " were found in ", path_scen, " - please check or define manually",
            " using argument 'replace_input_file_names'. Stopping."
          )
        }
        if (is.null(files_reference[[output]])) {
          stop(
            "None of the default file names for ", output,
            " were found in ", path_ref, " - please check or define manually",
            " using argument 'replace_input_file_names'. Stopping."
          )
        }
      } else {
        files_scenario[[output]] <- paste0(
          path_scen,
          replace_input_file_names[[output]], ".", file_extension
        )
        files_reference[[output]] <- paste0(
          path_ref,
          replace_input_file_names[[output]], ".", file_extension
        )
      }
    }
  } # end if !read_saved_data
  if (overtime && (window != nyears)) stop("Overtime is enabled, but window \
                  length (", window, ") does not match the reference nyears.")

  if (read_saved_data) {
    if (!is.null(save_data)) {
      message("Loading saved data from:", save_data)
      load(file = save_data)
    } else {
      stop(
        "save_data is not specified as parameter, ",
        "nothing to load ... exiting"
      )
    }
  } else {
    # first read in all lpjml output files required for computing EcoRisks
    returned_vars <- read_ecorisk_data(
      files_reference = files_reference,
      files_scenario = files_scenario,
      save_file = save_data,
      nitrogen = nitrogen,
      time_span_reference = time_span_reference,
      time_span_scenario = time_span_scenario,
      debug_mode = debug_mode,
      suppress_warnings = suppress_warnings
    )
    if (overtime == FALSE && slices > 1) {
      print("Overtime == FALSE, but too many time slices.
            I assume this is just a read all data operation. Exiting.")
      return(NULL)
    }

    # extract variables from return list object and give them proper names
    state_ref <- returned_vars$state_ref
    state_scen <- returned_vars$state_scen
    fpc_ref <- returned_vars$fpc_ref
    fpc_scen <- returned_vars$fpc_scen
    bft_ref <- returned_vars$bft_ref
    bft_scen <- returned_vars$bft_scen
    cft_ref <- returned_vars$cft_ref
    cft_scen <- returned_vars$cft_scen
    lat <- returned_vars$lat
    lon <- returned_vars$lon
    cell_area <- returned_vars$cell_area
    rm(returned_vars)
  }

  ncells <- length(cell_area)
  cell_ids <- 0:(ncells - 1)

  ecorisk <- list(
    ecorisk_total = array(0,
      dim = c(ncells, slices),