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 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 suppressWarnings 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 = FALSE,
                            replace_input_file_names = NULL,
                            external_variability = FALSE,
                            c2vr = NULL,
                            suppressWarnings = TRUE) {
  # check timespan consistency
  nyears <- length(time_span_reference)
  nyears_scen <- length(time_span_scenario)
  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_names <- metric_files$file_name
  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 = debug,
      suppressWarnings = suppressWarnings
    )
    # 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)
  slices <- (nyears_scen - window + 1)
  # todo: add dimnames already here and then get rid of the idexing
  ecorisk <- list(
    ecorisk_total = array(0, dim = c(ncells, slices)),
    vegetation_structure_change = array(0, dim = c(ncells, slices)),
    local_change = array(0, dim = c(ncells, slices)),
    global_importance = array(0, dim = c(ncells, slices)),
    ecosystem_balance = array(0, dim = c(ncells, slices)),
    c2vr = array(0, dim = c(4, ncells, slices)),
    carbon_stocks = array(0, dim = c(ncells, slices)),
    carbon_fluxes = array(0, dim = c(ncells, slices)),
    carbon_total = array(0, dim = c(ncells, slices)),
    water_total = array(0, dim = c(ncells, slices)),
    water_fluxes = array(0, dim = c(ncells, slices)),
    nitrogen_stocks = array(0, dim = c(ncells, slices)),
    nitrogen_fluxes = array(0, dim = c(ncells, slices)),
    nitrogen_total = array(0, dim = c(ncells, slices)),
    lat = lat,
    lon = lon
  )
  for (y in seq_len(slices)) {
    message("Calculating time slice ", y, " of ", slices)
    returned <- calc_ecorisk(
      fpc_ref = fpc_ref,
      fpc_scen = fpc_scen[, , y:(y + window - 1)],
      bft_ref = bft_ref,
      bft_scen = bft_scen[, , y:(y + window - 1)],
      cft_ref = cft_ref,
      cft_scen = cft_scen[, , y:(y + window - 1)],
      state_ref = state_ref,
      state_scen = state_scen[, y:(y + window - 1), ],
      weighting = weighting,
      lat = lat,
      lon = lon,
      cell_area = cell_area,
      dimensions_only_local = dimensions_only_local,
      nitrogen = nitrogen,
      external_variability = external_variability,
      c2vr = c2vr
    )
    ecorisk$ecorisk_total[, y] <- returned$ecorisk_total
    ecorisk$vegetation_structure_change[, y] <- (
      returned$vegetation_structure_change
    )
    ecorisk$local_change[, y] <- returned$local_change
    ecorisk$global_importance[, y] <- returned$global_importance
    ecorisk$ecosystem_balance[, y] <- returned$ecosystem_balance
    ecorisk$c2vr[, , y] <- returned$c2vr
    ecorisk$carbon_stocks[, y] <- returned$carbon_stocks
    ecorisk$carbon_fluxes[, y] <- returned$carbon_fluxes
    ecorisk$carbon_total[, y] <- returned$carbon_total
    ecorisk$water_total[, y] <- returned$water_total
    ecorisk$water_fluxes[, y] <- returned$water_fluxes
    if (nitrogen) {
      ecorisk$nitrogen_stocks[, y] <- returned$nitrogen_stocks
      ecorisk$nitrogen_fluxes[, y] <- returned$nitrogen_fluxes
      ecorisk$nitrogen_total[, y] <- returned$nitrogen_total
    }
  }

  ############## export and save data if requested #############
  if (!(is.null(save_ecorisk))) {
    message("Saving EcoRisk data to: ", save_ecorisk)
    save(ecorisk, file = save_ecorisk)
  }
  #
  ###
  return(ecorisk)
}

#' Calculate the ecosystem change metric EcoRisk between 2 sets of states
#' This function is called by the wrapper function (ecorisk_wrapper),
#' unless you know what you are doing, don't use this function directly.
#'
#' Function to calculate the ecosystem change metric EcoRisk, based on
#' gamma/vegetation_structure_change
#' work from Sykes (1999), Heyder (2011), and Ostberg (2015,2018).
#' This is a reformulated version in R, not producing 100% similar values
#' than the C/bash version from Ostberg et al. 2018, but similar the methodology
#'
#' @param fpc_ref reference run data for fpc
#' @param fpc_scen scenario run data for fpc
#' @param bft_ref reference run data for fpc_bft
#' @param bft_scen scenario run data for fpc_bft
#' @param cft_ref reference run data for cftfrac
#' @param cft_scen scenario run data for cftfrac
#' @param state_ref reference run data for state variables
#' @param state_scen scenario run data for state variables
#' @param weighting apply "old" (Ostberg-like), "new", or "equal" weighting of
#'        vegetation_structure_change weights (default "equal")
#' @param lat latitude array
#' @param lon longitude array
#' @param cell_area cellarea array
#' @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 nitrogen include nitrogen outputs (default: TRUE)
#' @param external_variability include external change_to_variability_ratio?
#'        (default: FALSE)
#' @param c2vr list with external change_to_variability_ratios for each
#'        component (default: NULL)
#'
#' @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)
#'
#' @export
calc_ecorisk <- function(fpc_ref,
                         fpc_scen,
                         bft_ref,
                         bft_scen,
                         cft_ref,
                         cft_scen,
                         state_ref,
                         state_scen,
                         weighting = "equal",
                         lat,
                         lon,
                         cell_area,
                         dimensions_only_local = FALSE,
                         nitrogen = TRUE,
                         external_variability = FALSE,
                         c2vr = NULL) {
  if (external_variability && is.null(c2vr)) {
    stop("external_variability enabled, but not supplied (c2vr). Aborting.")
  }
  di_ref <- dim(fpc_ref)
  di_scen <- dim(fpc_scen)
  ncells <- di_ref[1]
  nyears <- di_ref[3]
  if (di_ref[3] != di_scen[3]) {
    stop("Dimension year does not match between fpc_scen and fpc_ref.")
  }
  # calc vegetation_structure_change and variability of
  #   vegetation_structure_change within
  # reference period S(vegetation_structure_change,
  #   sigma_vegetation_structure_change)
  fpc_ref_mean <- apply(fpc_ref, c(1, 2), mean)
  bft_ref_mean <- apply(bft_ref, c(1, 2), mean)
  cft_ref_mean <- apply(cft_ref, c(1, 2), mean)


  sigma_vegetation_structure_change_ref_list <- array(
    0,
    dim = c(ncells, nyears)
  )
  # calculate for every year of the reference period,
  #   vegetation_structure_change between that year and the average reference
  #   period year
  # this gives the variability of vegetation_structure_change within the
  #   reference period
  for (y in seq_len(nyears)) {
    sigma_vegetation_structure_change_ref_list[, y] <- calc_delta_v( # nolint
      fpc_ref = fpc_ref_mean,
      fpc_scen = fpc_ref[, , y],
      bft_ref = bft_ref_mean,
      bft_scen = bft_ref[, , y],
      cft_ref = cft_ref_mean,
      cft_scen = cft_ref[, , y],
      weighting = weighting
    )
  }

  # calculate the std deviation over the reference period for each gridcell
  vegetation_structure_changesd <- apply(
    sigma_vegetation_structure_change_ref_list,
    c(1),
    stats::sd
  )

  # calculate vegetation_structure_change between average reference and average
  #   scenario period
  vegetation_structure_change <- calc_delta_v(
    fpc_ref = fpc_ref_mean,
    fpc_scen = apply(fpc_scen, c(1, 2), mean),
    bft_ref = bft_ref_mean,
    bft_scen = apply(bft_scen, c(1, 2), mean),
    cft_ref = cft_ref_mean,
    cft_scen = apply(cft_scen, c(1, 2), mean),
    weighting = weighting
  )
  #
  ####
  ############## calc EcoRisk components ################
  # dimensions in the state vector
  # 1 "vegetation_carbon_pool"
  # 2 "soil_carbon_pool"
  # 3 "carbon_influx"
  # 4 "carbon_outflux"
  # 5 "soil_water_pool"
  # 6 "water_influx"
  # 7 "water_outflux"
  # 8 "other"
  # 9 "vegetation_nitrogen_pool"
  # 10 "soil_mineral_nitrogen_pool"
  # 11 "nitrogen_influx"
  # 12 "nitrogen_outflux"

  delta_var <- s_change_to_var_ratio(
    vegetation_structure_change,
    vegetation_structure_changesd
  )
  nitrogen_dimensions <- c("vegetation_nitrogen_pool",
                           "soil_mineral_nitrogen_pool",
                           "nitrogen_influx",
                           "nitrogen_outflux")
  all_dimensions <- dimnames(state_scen)$class
  non_nitrogen_dimensions <- setdiff(all_dimensions, nitrogen_dimensions)
  if (nitrogen) {
    lc_raw <- calc_component(
      ref = state_ref,
      scen = state_scen,
      local = TRUE,
      cell_area = cell_area
    ) # local change

    gi_raw <- calc_component(
      ref = state_ref,
      scen = state_scen,
      local = FALSE,
      cell_area = cell_area
    ) # global importance

    eb_raw <- calc_ecosystem_balance(
      ref = state_ref,
      scen = state_scen
    ) # ecosystem balance
  } else {
    lc_raw <- calc_component(
      ref = state_ref[, , non_nitrogen_dimensions],
      scen = state_scen[, , non_nitrogen_dimensions],
      local = TRUE,
      cell_area = cell_area
    ) # local change

    gi_raw <- calc_component(
      ref = state_ref[, , non_nitrogen_dimensions],
      scen = state_scen[, , non_nitrogen_dimensions],
      local = FALSE,
      cell_area = cell_area
    ) # global importance

    eb_raw <- calc_ecosystem_balance(
      ref = state_ref[, , non_nitrogen_dimensions],
      scen = state_scen[, , non_nitrogen_dimensions]
    ) # ecosystem balance
  }
  if (dimensions_only_local == TRUE) {
    # carbon stocks (local change)
    cs <- calc_component(
      ref = state_ref[, , c("vegetation_carbon_pool", "soil_carbon_pool")],
      scen = state_scen[, , c("vegetation_carbon_pool", "soil_carbon_pool")],
      local = TRUE,
      cell_area = cell_area
    )$full
    # carbon fluxes (local change)
    cf <- calc_component(
      ref = state_ref[, , c("carbon_influx", "carbon_outflux")],
      scen = state_scen[, , c("carbon_influx", "carbon_outflux")],
      local = TRUE,
      cell_area = cell_area
    )$full
    # total carbon (local change)
    ct <- calc_component(
      ref = state_ref[, , c("vegetation_carbon_pool",
                            "soil_carbon_pool",
                            "carbon_influx",
                            "carbon_outflux")],
      scen = state_scen[, , c("vegetation_carbon_pool",
                              "soil_carbon_pool",
                              "carbon_influx",
                              "carbon_outflux")],
      local = TRUE,
      cell_area = cell_area
    )$full
    # water fluxes (local change)
    wf <- calc_component(
      ref = state_ref[, , c("water_influx",
                            "water_outflux")],
      scen = state_scen[, , c("water_influx",
                              "water_outflux")],
      local = TRUE,
      cell_area = cell_area
    )$full
    # total water (local change)
    wt <- calc_component(
      ref = state_ref[, , c("water_influx",
                            "water_outflux",
                            "soil_water_pool")],
      scen = state_scen[, , c("water_influx",
                              "water_outflux",
                              "soil_water_pool")],
      local = TRUE,
      cell_area = cell_area
    )$full

    # nitrogen stocks (local change)
    if (nitrogen) {
      ns <- calc_component(
        ref = state_ref[, , c("vegetation_nitrogen_pool",
                              "soil_mineral_nitrogen_pool")],
        scen = state_scen[, , c("vegetation_nitrogen_pool",
                                "soil_mineral_nitrogen_pool")],
        local = TRUE,
        cell_area = cell_area
      )$full
      # nitrogen fluxes (local change)
      nf <- calc_component(
        ref = state_ref[, , c("nitrogen_influx", "nitrogen_outflux")],
        scen = state_scen[, , c("nitrogen_influx", "nitrogen_outflux")],
        local = TRUE,
        cell_area = cell_area
      )$full
      # total nitrogen (local change)
      nt <- calc_component(
        ref = state_ref[, , c("vegetation_nitrogen_pool",
                              "soil_mineral_nitrogen_pool",
                              "nitrogen_influx",
                              "nitrogen_outflux")],
        scen = state_scen[, , c("vegetation_nitrogen_pool",
                                "soil_mineral_nitrogen_pool",
                                "nitrogen_influx",
                                "nitrogen_outflux")],
        local = TRUE,
        cell_area = cell_area
      )$full
    }
  } else { # local == FALSE
    cf <- (
      calc_component(
        ref = state_ref[, , c("carbon_influx",
                              "carbon_outflux")],
        scen = state_scen[, , c("carbon_influx",
                                "carbon_outflux")],
        local = TRUE,
        cell_area = cell_area
      )$full + # carbon fluxes
        calc_component(
          ref = state_ref[, , c("carbon_influx",
                                "carbon_outflux")],
          scen = state_scen[, , c("carbon_influx",
                                  "carbon_outflux")],
          local = FALSE,
          cell_area = cell_area
        )$full +
        calc_ecosystem_balance(
          ref = state_ref[, , c("carbon_influx",
                                "carbon_outflux")],
          scen = state_scen[, , c("carbon_influx",
                                  "carbon_outflux")]
        )$full
    ) / 3

    # carbon stocks
    cs <- (
      calc_component(
        ref = state_ref[, , c("vegetation_carbon_pool",
                              "soil_carbon_pool")],
        scen = state_scen[, , c("vegetation_carbon_pool",
                                "soil_carbon_pool")],
        local = TRUE,
        cell_area = cell_area
      )$full +
        calc_component(
          ref = state_ref[, , c("vegetation_carbon_pool",
                                "soil_carbon_pool")],
          scen = state_scen[, , c("vegetation_carbon_pool",
                                  "soil_carbon_pool")],
          local = FALSE,
          cell_area = cell_area
        )$full +
        calc_ecosystem_balance(
          ref = state_ref[, , c("vegetation_carbon_pool",
                                "soil_carbon_pool")],
          scen = state_scen[, , c("vegetation_carbon_pool",
                                  "soil_carbon_pool")]
        )$full
    ) / 3

    # carbon total
    ct <- (
      calc_component(
        ref = state_ref[, , c("vegetation_carbon_pool",
                              "soil_carbon_pool",
                              "carbon_influx",
                              "carbon_outflux")],
        scen = state_scen[, , c("vegetation_carbon_pool",
                                "soil_carbon_pool",
                                "carbon_influx",
                                "carbon_outflux")],
        local = TRUE,
        cell_area = cell_area
      )$full +
        calc_component(
          ref = state_ref[, , c("vegetation_carbon_pool",
                                "soil_carbon_pool",
                                "carbon_influx",
                                "carbon_outflux")],
          scen = state_scen[, , c("vegetation_carbon_pool",
                                  "soil_carbon_pool",
                                  "carbon_influx",
                                  "carbon_outflux")],
          local = FALSE,
          cell_area = cell_area
        )$full +
        calc_ecosystem_balance(
          ref = state_ref[, , c("vegetation_carbon_pool",
                                "soil_carbon_pool",
                                "carbon_influx",
                                "carbon_outflux")],
          scen = state_scen[, , c("vegetation_carbon_pool",
                                  "soil_carbon_pool",
                                  "carbon_influx",
                                  "carbon_outflux")]
        )$full
    ) / 3

    # water fluxes
    wf <- (
      calc_component(
        ref = state_ref[, , c("water_influx",
                              "water_outflux")],
        scen = state_scen[, , c("water_influx",
                                "water_outflux")],
        local = TRUE,
        cell_area = cell_area
      )$full +
        calc_component(
          ref = state_ref[, , c("water_influx",
                                "water_outflux")],
          scen = state_scen[, , c("water_influx",
                                  "water_outflux")],
          local = FALSE,
          cell_area = cell_area
        )$full + calc_ecosystem_balance(
          ref = state_ref[, , c("water_influx",
                                "water_outflux")],
          scen = state_scen[, , c("water_influx",
                                  "water_outflux")]
        )$full
    ) / 3

    # water total
    wt <- (
      calc_component(
        ref = state_ref[, , c("water_influx",
                              "water_outflux",
                              "soil_water_pool")],
        scen = state_scen[, , c("water_influx",
                                "water_outflux",
                                "soil_water_pool")],
        local = TRUE,
        cell_area = cell_area
      )$full +
        calc_component(
          ref = state_ref[, , c("water_influx",
                                "water_outflux",
                                "soil_water_pool")],
          scen = state_scen[, , c("water_influx",
                                  "water_outflux",
                                  "soil_water_pool")],
          local = FALSE,
          cell_area = cell_area
        )$full +
        calc_ecosystem_balance(
          ref = state_ref[, , c("water_influx",
                                "water_outflux",
                                "soil_water_pool")],
          scen = state_scen[, , c("water_influx",
                                  "water_outflux",
                                  "soil_water_pool")]
        )$full
    ) / 3

    if (nitrogen) {
      # nitrogen stocks (local change)
      ns <- (
        calc_component(
          ref = state_ref[, , c("vegetation_nitrogen_pool",
                                "soil_mineral_nitrogen_pool")],
          scen = state_scen[, , c("vegetation_nitrogen_pool",
                                  "soil_mineral_nitrogen_pool")],
          local = TRUE,
          cell_area = cell_area
        )$full +
          calc_component(
            ref = state_ref[, , c("vegetation_nitrogen_pool",
                                  "soil_mineral_nitrogen_pool")],
            scen = state_scen[, , c("vegetation_nitrogen_pool",
                                    "soil_mineral_nitrogen_pool")],
            local = FALSE, cell_area = cell_area
          )$full +
          calc_ecosystem_balance(
            ref = state_ref[, , c("vegetation_nitrogen_pool",
                                  "soil_mineral_nitrogen_pool")],
            scen = state_scen[, , c("vegetation_nitrogen_pool",
                                    "soil_mineral_nitrogen_pool")]
          )$full
      ) / 3

      # nitrogen fluxes (local change)
      nf <- (
        calc_component(
          ref = state_ref[, , c("nitrogen_influx",
                                "nitrogen_outflux")],
          scen = state_scen[, , c("nitrogen_influx",
                                  "nitrogen_outflux")],
          local = TRUE,
          cell_area = cell_area
        )$full +
          calc_component(
            ref = state_ref[, , c("nitrogen_influx",
                                  "nitrogen_outflux")],
            scen = state_scen[, , c("nitrogen_influx",
                                    "nitrogen_outflux")],
            local = FALSE,
            cell_area = cell_area
          )$full +
          calc_ecosystem_balance(
            ref = state_ref[, , c("nitrogen_influx",
                                  "nitrogen_outflux")],
            scen = state_scen[, , c("nitrogen_influx",
                                    "nitrogen_outflux")]
          )$full
      ) / 3

      nt <- (
        calc_component(
          ref = state_ref[, , c("vegetation_nitrogen_pool",
                                "soil_mineral_nitrogen_pool",
                                "nitrogen_influx",
                                "nitrogen_outflux")],
          scen = state_scen[, , c("vegetation_nitrogen_pool",
                                  "soil_mineral_nitrogen_pool",
                                  "nitrogen_influx",
                                  "nitrogen_outflux")],
          local = TRUE,
          cell_area = cell_area
        )$full +
          calc_component(
            ref = state_ref[, , c("vegetation_nitrogen_pool",
                                  "soil_mineral_nitrogen_pool",
                                  "nitrogen_influx",
                                  "nitrogen_outflux")],
            scen = state_scen[, , c("vegetation_nitrogen_pool",
                                    "soil_mineral_nitrogen_pool",
                                    "nitrogen_influx",
                                    "nitrogen_outflux")],
            local = FALSE,
            cell_area = cell_area
          )$full +
          calc_ecosystem_balance(
            ref = state_ref[, , c("vegetation_nitrogen_pool",
                                  "soil_mineral_nitrogen_pool",
                                  "nitrogen_influx",
                                  "nitrogen_outflux")],
            scen = state_scen[, , c("vegetation_nitrogen_pool",
                                    "soil_mineral_nitrogen_pool",
                                    "nitrogen_influx",
                                    "nitrogen_outflux")]
          )$full
      ) / 3
    }
  }
  if (external_variability) {
    delta <- vegetation_structure_change * c2vr["vs", ] # vegetation_structure
    lc <- lc_raw$value * c2vr["lc", ]
    gi <- gi_raw$value * c2vr["gi", ]
    eb <- eb_raw$value * c2vr["eb", ]
  } else {
    delta <- vegetation_structure_change * delta_var # vegetation_structure
    lc <- lc_raw$value * lc_raw$var
    gi <- gi_raw$value * gi_raw$var
    eb <- eb_raw$value * eb_raw$var
    c2vr <- rbind(delta_var, lc_raw$var, gi_raw$var, eb_raw$var) # dim=(4,ncell)
    dimnames(c2vr) <- list(component = c("vs", "lc", "gi", "eb"),
                           cell = 0:(ncells - 1))
  }

  # calc total EcoRisk as the average of the 4 components
  ecorisk_full <- (delta + lc + gi + eb) / 4 # check for NAs

  if (nitrogen) {
    ecorisk <- list(
      ecorisk_total = ecorisk_full,
      vegetation_structure_change = delta,
      local_change = lc,
      global_importance = gi,
      ecosystem_balance = eb,
      c2vr = c2vr,
      carbon_stocks = cs,
      carbon_fluxes = cf,
      carbon_total = ct,
      water_fluxes = wf,
      water_stocks = NA,
      water_total = wt,
      nitrogen_stocks = ns,
      nitrogen_fluxes = nf,
      nitrogen_total = nt
    )
  } else {
    ecorisk <- list(
      ecorisk_total = ecorisk_full,
      vegetation_structure_change = delta,
      local_change = lc,
      global_importance = gi,
      ecosystem_balance = eb,
      c2vr = c2vr,
      carbon_stocks = cs,
      carbon_fluxes = cf,
      carbon_total = ct,
      water_fluxes = wf,
      water_stocks = NA,
      water_total = wt,
      nitrogen_stocks = NA,
      nitrogen_fluxes = NA,
      nitrogen_total = NA
    )
  }
  ###
  return(ecorisk)
}

#' Read in output data from LPJmL to calculate the ecosystem change metric
#' EcoRisk. This function is called by the wrapper function (ecorisk_wrapper),
#' unless you know what you are doing, don't use this function directly.
#'
#' Utility function to read in output data from LPJmL for calculation of EcoRisk
#'
#' @param files_reference folder of reference run
#' @param files_scenario folder of scenario run
#' @param save_file file to save read in data to (default NULL)
#' @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 nitrogen include nitrogen outputs for pools and fluxes into EcoRisk
#'                 calculation (default FALSE)
#' @param debug write out all nitrogen state variables (default FALSE)
#' @param suppressWarnings suppress writing of Warnings, default: TRUE
#'
#' @return list data object containing arrays of state_ref, mean_state_ref,
#'         state_scen, mean_state_scen, fpc_ref, fpc_scen, bft_ref, bft_scen,
#'         cft_ref, cft_scen, lat, lon, cell_area
#'
#' @export
read_ecorisk_data <- function(
    files_reference, # nolint
    files_scenario,
    save_file = NULL,
    time_span_reference,
    time_span_scenario,
    nitrogen,
    debug = FALSE,
    suppressWarnings = TRUE) {
  file_type <- tools::file_ext(files_reference$grid)

  if (file_type %in% c("json", "clm")) {
    # read grid
    grid <- lpjmlkit::read_io(
      files_reference$grid
    ) %>% suppressWarnings()
    # calculate cell area
    cell_area <- drop(lpjmlkit::read_io(
      filename = files_reference$terr_area
    )$data) %>% suppressWarnings() # in m2
    lat <- grid$data[, , 2]
    lon <- grid$data[, , 1]
    ncells <- length(lat)
    nyears <- length(time_span_scenario)
    nyears_ref <- length(time_span_reference)

    ### read in lpjml output
    # for vegetation_structure_change (fpc,fpc_bft,cftfrac)
    message("Reading in fpc, fpc_bft, cftfrac")

    cft_scen <- aperm(lpjmlkit::read_io(
      files_scenario$cftfrac,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)), c(1, 3, 2)) %>%
      suppressWarnings()

    bft_scen <- aperm(lpjmlkit::read_io(
      files_scenario$fpc_bft,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)), c(1, 3, 2)) %>%
      suppressWarnings()

    fpc_scen <- aperm(lpjmlkit::read_io(
      files_scenario$fpc,
      subset = list(year = as.character(time_span_scenario))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)), c(1, 3, 2)) %>%
      suppressWarnings()

    if (file.exists(files_reference$cftfrac)) {
      cft_ref <- aperm(lpjmlkit::read_io(
        files_reference$cftfrac,
        subset = list(year = as.character(time_span_reference))
      ) %>%
        lpjmlkit::transform(to = c("year_month_day")) %>%
        lpjmlkit::as_array(aggregate = list(month = sum)), c(1, 3, 2)) %>%
        suppressWarnings()
    } else {
      cft_ref <- cft_scen * 0
    }

    if (file.exists(files_reference$fpc_bft)) {
      bft_ref <- aperm(lpjmlkit::read_io(
        files_reference$fpc_bft,
        subset = list(year = as.character(time_span_reference))
      ) %>%
        lpjmlkit::transform(to = c("year_month_day")) %>%
        lpjmlkit::as_array(aggregate = list(month = sum)), c(1, 3, 2)) %>%
        suppressWarnings()
    } else {
      bft_ref <- bft_scen * 0
    }

    fpc_ref <- aperm(lpjmlkit::read_io(
      files_reference$fpc,
      subset = list(year = as.character(time_span_reference))
    ) %>%
      lpjmlkit::transform(to = c("year_month_day")) %>%
      lpjmlkit::as_array(aggregate = list(month = sum)), c(1, 3, 2)) %>%
      suppressWarnings()

    #### new input reading ###
    metric_files <- system.file(
      "extdata",
      "metric_files.yml",
      package = "biospheremetrics"
    ) %>%
      yaml::read_yaml()
    nclasses <- length(metric_files$metric$ecorisk_nitrogen$metric_class)
    nstate_dimensions <- 0
    for (i in seq_len(nclasses)) {
      nstate_dimensions <- nstate_dimensions +
        length(metric_files$metric$ecorisk_nitrogen$metric_class[[i]])
    }
    state_ref <- array(0, dim = c(ncells, nyears_ref, nstate_dimensions))
    state_scen <- array(0, dim = c(ncells, nyears, nstate_dimensions))
    class_names <- seq_len(nstate_dimensions)
    index <- 1
    # iterate over main classes (carbon pools, water fluxes ...)
    for (c in seq_len(nclasses)) {
      classe <- metric_files$metric$ecorisk_nitrogen$metric_class[[c]]
      nsubclasses <- length(classe)
      # iterate over subclasses (vegetation carbon, soil water ...)
      for (s in seq_len(nsubclasses)) {
        subclass <- classe[s]
        class_names[index] <- names(subclass)
        vars <- split_sign(unlist(subclass))
        for (v in seq_len(length(vars[, 1]))) {
          path_scen_file <- files_scenario[[vars[v, "variable"]]]
          if (file.exists(path_scen_file)) {
            header_scen <- lpjmlkit::read_meta(filename = path_scen_file)
            message(
              "Reading in ", path_scen_file, " with unit ", header_scen$unit,
              " -> as part of ", class_names[index]
            )
            var_scen <- lpjmlkit::read_io(
              path_scen_file,
              subset = list(year = as.character(time_span_scenario))
            ) %>%
              lpjmlkit::transform(to = c("year_month_day")) %>%
              lpjmlkit::as_array(aggregate = list(month = sum, band = sum), )%>%
              drop() %>%
              suppressWarnings()
          } else {
            stop(paste("Couldn't read in:", path_scen_file, " - stopping!"))
          }
          path_ref_file <- files_reference[[vars[v, "variable"]]]
          if (file.exists(path_ref_file)) {
            header_ref <- lpjmlkit::read_meta(path_ref_file)
            message(
              "Reading in ", path_ref_file, " with unit ", header_ref$unit,
              " -> as part of ", class_names[index]
            )
            var_ref <- lpjmlkit::read_io(
              path_ref_file,
              subset = list(year = as.character(time_span_reference))
            ) %>%
              lpjmlkit::transform(to = c("year_month_day")) %>%
              lpjmlkit::as_array(aggregate = list(month = sum, band = sum)) %>%
              drop() %>%
              suppressWarnings()
          } else {
            stop(paste("Couldn't read in:", path_ref_file, " - stopping!"))
          }
          # if (window > 30){
          #  if (vars[v,"sign"] == "+"){
          #    state_scen[,,index,] <- state_scen[,,index,] + var_scen
          #    state_ref[,,index,] <- state_ref[,,index,] + var_ref
          #  } else { # vars[v,"sign"] == "-"
          #    state_scen[,,index,] <- state_scen[,,index,] - var_scen
          #    state_ref[,,index,] <- state_ref[,,index,] - var_ref
          #  }
          # }else{
          if (vars[v, "sign"] == "+") {
            state_scen[, , index] <- state_scen[, , index] + var_scen
            state_ref[, , index] <- state_ref[, , index] + var_ref
          } else { # vars[v,"sign"] == "-"
            state_scen[, , index] <- state_scen[, , index] - var_scen
            state_ref[, , index] <- state_ref[, , index] - var_ref
          }
          # }
        }
        index <- index + 1
      }
    }

    dimnames(state_scen) <- list(cell = 0:(ncells - 1),
                  year = as.character(time_span_scenario), class = class_names)
    dimnames(state_ref) <- list(cell = 0:(ncells - 1),
                  year = as.character(time_span_reference), class = class_names)
  } else if (file_type == "nc") { # to be added
    stop(
      "nc reading has not been updated to latest functionality. ",
      "Please contact Fabian Stenzel"