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 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,
                            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."
    )
  }

  # 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
    files_scenario[[output]] <- paste0(path_scen, metric_files$file_name[[output]][1], ".", file_extension)
    files_reference[[output]] <- paste0(path_ref, metric_files$file_name[[output]][1], ".", file_extension)
  }
  
  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)
  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_change
    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_change
    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,ncells)
    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)

    ### 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, 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"
    )
  } else {
    stop("Unrecognized file type (", file_type, ")")
  }

  if (!(is.null(save_file))) {
    message("Saving data to: ", save_file)
    save(state_ref, state_scen, fpc_ref, fpc_scen,
      bft_ref, bft_scen, cft_ref, cft_scen, lat, lon, cell_area,
      file = save_file
    )
  }
  return(
    list(
      state_ref = state_ref,
      state_scen = state_scen,
      fpc_ref = fpc_ref,
      fpc_scen = fpc_scen,
      bft_ref = bft_ref,
      bft_scen = bft_scen,
      cft_ref = cft_ref,
      cft_scen = cft_scen,
      lat = lat,
      lon = lon,
      cell_area = cell_area
    )
  )
}

#' Calculates changes in vegetation structure (vegetation_structure_change)
#'
#' Utility function to calculate changes in vegetation structure
#' (vegetation_structure_change) for calculation of EcoRisk
#'
#' @param fpc_ref reference fpc array (dim: [ncells,npfts+1])
#' @param fpc_scen scenario fpc array (dim: [ncells,npfts+1])
#' @param bft_ref reference bft array (dim: [ncells,nbfts])
#' @param bft_scen scenario bft array (dim: [ncells,nbfts])
#' @param cft_ref reference cft array (dim: [ncells,ncfts])
#' @param cft_scen scenario cft array (dim: [ncells,ncfts])
#' @param weighting apply "old" (Ostberg-like), "new", or "equal" weighting of
#'                  vegetation_structure_change weights (default "equal")
#'
#' @return vegetation_structure_change array of size ncells with the
#'         vegetation_structure_change value [0,1] for each cell
#'
#' @examples
#' \dontrun{
#' 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 = "equal"
#' )
#' }
#' @export
calc_delta_v <- function(fpc_ref, # nolint
                         fpc_scen,
                         bft_ref,
                         bft_scen,
                         cft_ref,
                         cft_scen,
                         weighting = "equal") {
  di <- dim(fpc_ref)
  ncells <- di[1]
  npfts <- di[2] - 1

  fpc_ref[fpc_ref < 0] <- 0
  fpc_scen[fpc_scen < 0] <- 0
  bft_ref[bft_ref < 0] <- 0
  bft_scen[bft_scen < 0] <- 0
  cft_ref[cft_ref < 0] <- 0
  cft_scen[cft_scen < 0] <- 0

  if (npfts == 9) {
    # barren = 1 - crop area - natural vegetation area +
    #   barren under bioenergy trees
    barren_area_ref <- (
      1 - rowSums(cft_ref) -
        rowSums(fpc_ref[, 2:10]) * fpc_ref[, 1] +
        rowSums(cft_ref[, c(16, 32)]) * (1 - rowSums(bft_ref[, c(1:4, 7:10)]))
    )

    barren_area_ref[barren_area_ref < 0] <- 0

    tree_area_ref <- array(0, dim = c(ncells, 11))

    # natural tree area fractions scaled by total natural frac
    tree_area_ref[, 1:7] <- (
      fpc_ref[, 2:8] * fpc_ref[, 1]
    )

    # fraction of rainfed tropical and temperate BE trees scaled by total
    #   rainfed bioenergy tree area and relative fpc of bioenergy trees and
    #   grass under bioenergy trees
    tree_area_ref[, 8:9] <- (
      cft_ref[, 16] * bft_ref[, 1:2] / rowSums(bft_ref[, c(1, 2, 4)])
    )

    # fraction of irrigated tropical and temperate BE trees scaled by total
    #   irrigated bioenergy tree area and relative fpc of bioenergy trees and
    #   grass under bioenergy trees
    tree_area_ref[, 10:11] <- (
      cft_ref[, 32] * bft_ref[, 7:8] / rowSums(bft_ref[, c(7, 8, 10)])
    )

    grass_area_ref <- array(0, dim = c(ncells, 20))

    # natural grass
    grass_area_ref[, 1:2] <- fpc_ref[, 9:10] * fpc_ref[, 1]

    # crops
    grass_area_ref[, 3:15] <- cft_ref[, 1:13] + cft_ref[, 17:29]

    # managed grass rf
    grass_area_ref[, 16] <- cft_ref[, 14]

    # managed grass irr
    grass_area_ref[, 17] <- cft_ref[, 30]

    # bioenergy grass
    grass_area_ref[, 18] <- cft_ref[, 15] + cft_ref[, 31]

    # fraction of rainfed grass under bioenergy trees
    grass_area_ref[, 19] <- (
      cft_ref[, 16] * bft_ref[, 4] / rowSums(bft_ref[, c(1, 2, 4)])
    )

    # fraction of irrigated grass under bioenergy trees
    grass_area_ref[, 20] <- (
      cft_ref[, 32] * bft_ref[, 10] / rowSums(bft_ref[, c(7, 8, 10)])
    )

    # barren
    barren_area_scen <- (
      1 - rowSums(cft_scen) -
        rowSums(fpc_scen[, 2:10]) * fpc_scen[, 1] +