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 varnames data.frame with names of output files (outname) and time res.
#'        (timestep) -- can be specified to account for variable file names
#'        (default NULL -- standard names as below)
#' @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)
#'
#' @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
ecorisk_wrapper <- function(path_ref,
                            path_scen,
                            read_saved_data = FALSE,
                            save_data = NULL,
                            save_ecorisk = NULL,
                            nitrogen = TRUE,
                            weighting = "equal",
                            varnames = NULL,
                            time_span_reference,
                            time_span_scenario,
                            dimensions_only_local = FALSE,
                            overtime = FALSE,
                            window = 30,
                            debug = FALSE,
                            external_variability = FALSE,
                            c2vr = NULL) {
  if (is.null(varnames)) {
    print("variable name list not provided, using standard list, which might
          not be applicable for this case ...")
    varnames <- data.frame(
        row.names = c("grid","fpc", "fpc_bft", "cftfrac", "firec", "npp", "runoff",
                "transp", "vegc", "firef", "rh", "harvestc", "rharvestc",
                "pft_harvestc", "pft_rharvestc", "evap", "interc", "discharge",
                "soilc", "litc", "swc", "vegn", "soilnh4", "soilno3",
                "leaching", "n2o_denit", "n2o_nit", "n2_emis", "bnf",
                "n_volatilization", "gpp", "res_storage", "lakevol", "ndepos",
                "rd", "prec", "terr_area", "irrig", "nfert_agr", "nmanure_agr",
                "firen", "harvestn", "rivervol", "irrig_stor"),
        outname = c("grid.bin.json", "fpc.bin.json", "fpc_bft.bin.json",
              "cftfrac.bin.json", "firec.bin.json", "npp.bin.json",
              "runoff.bin.json", "transp.bin.json", "vegc.bin.json",
              "firef.bin.json", "rh.bin.json", "harvestc.bin.json",
              "rharvestc.bin.json", "pft_harvest.pft.bin.json",
              "pft_rharvest.pft.bin.json", "mevap.bin.json",
              "interc.bin.json", "discharge.bin.json", "soilc.bin.json",
              "litc.bin.json", "swc.bin.json", "vegn.bin.json",
              "soilnh4.bin.json", "soilno3.bin.json", "mleaching.bin.json",
              "n2o_denit.bin.json", "n2o_nit.bin.json", "n2_emis.bin.json",
              "bnf.bin.json", "n_volatilization.bin.json", "gpp.bin.json",
              "res_storage.bin.json", "lakevol.bin.json", "ndepos.bin.json",
              "rd.bin.json", "mprec.bin.json", "terr_area.bin.json",
              "irrig.bin.json", "nfert_agr.bin.json", "nmanure_agr.bin.json",
              "firen.bin.json", "harvestn.bin.json", "rivervol.bin.json",
              "irrig_stor.bin.json"),
        timestep = rep("Y", 44)
    )
  }

  nyears <- length(time_span_reference)
  nyears_scen <- length(time_span_scenario)
  if (nyears < 30 || nyears_scen < 30) {
    stop("Warning: timespan in reference or scenario is smaller than 30 years. \
          Aborting!")
  }
  # translate varnames and folders to files_scenarios/reference lists
  files_scenario <- list(
    grid = paste0(path_scen, varnames["grid", "outname"]),
    fpc = paste0(path_scen, varnames["fpc", "outname"]),
    fpc_bft = paste0(path_scen, varnames["fpc_bft", "outname"]),
    cftfrac = paste0(path_scen, varnames["cftfrac", "outname"]),
    firec = paste0(path_scen, varnames["firec", "outname"]),
    npp = paste0(path_scen, varnames["npp", "outname"]),
    runoff = paste0(path_scen, varnames["runoff", "outname"]),
    transp = paste0(path_scen, varnames["transp", "outname"]),
    vegc = paste0(path_scen, varnames["vegc", "outname"]),
    firef = paste0(path_scen, varnames["firef", "outname"]),
    rh = paste0(path_scen, varnames["rh", "outname"]),
    harvestc = paste0(path_scen, varnames["harvestc", "outname"]),
    rharvestc = paste0(path_scen, varnames["rharvestc", "outname"]),
    pft_harvestc = paste0(path_scen, varnames["pft_harvest", "outname"]),
    pft_rharvestc = paste0(path_scen, varnames["pft_rharvest", "outname"]),
    evap = paste0(path_scen, varnames["evap", "outname"]),
    interc = paste0(path_scen, varnames["interc", "outname"]),
    discharge = paste0(path_scen, varnames["discharge", "outname"]),
    soilc = paste0(path_scen, varnames["soilc", "outname"]),
    litc = paste0(path_scen, varnames["litc", "outname"]),
    swc = paste0(path_scen, varnames["swc", "outname"]),
    swc_vol = paste0(path_scen, varnames["swc_vol", "outname"]),
    vegn = paste0(path_scen, varnames["vegn", "outname"]),
    soilnh4 = paste0(path_scen, varnames["soilnh4", "outname"]),
    soilno3 = paste0(path_scen, varnames["soilno3", "outname"]),
    leaching = paste0(path_scen, varnames["leaching", "outname"]),
    n2o_denit = paste0(path_scen, varnames["n2o_denit", "outname"]),
    n2o_nit = paste0(path_scen, varnames["n2o_nit", "outname"]),
    n2_emis = paste0(path_scen, varnames["n2_emis", "outname"]),
    bnf = paste0(path_scen, varnames["bnf", "outname"]),
    n_volatilization = paste0(path_scen, varnames["n_volatilization", "outname"]),
    gpp = paste0(path_scen, varnames["gpp", "outname"]),
    res_storage = paste0(path_scen, varnames["res_storage", "outname"]),
    lakevol = paste0(path_scen, varnames["lakevol", "outname"]),
    ndepos = paste0(path_scen, varnames["ndepos", "outname"]),
    rd = paste0(path_scen, varnames["rd", "outname"]),
    prec = paste0(path_scen, varnames["prec", "outname"]),
    terr_area = paste0(path_scen, varnames["terr_area", "outname"]),
    irrig = paste0(path_scen, varnames["irrig", "outname"]),
    nfert_agr = paste0(path_scen, varnames["nfert_agr", "outname"]),
    nmanure_agr = paste0(path_scen, varnames["nmanure_agr", "outname"]),
    ndepos = paste0(path_scen, varnames["ndepos", "outname"]),
    firen = paste0(path_scen, varnames["firen", "outname"]),
    harvestn = paste0(path_scen, varnames["harvestn", "outname"]),
    irrig_stor = paste0(path_scen, varnames["irrig_stor", "outname"]),
    rivervol = paste0(path_scen, varnames["rivervol", "outname"]),
    rootmoist = paste0(path_scen, varnames["rootmoist", "outname"])
  )
  files_reference <- list(
    grid = paste0(path_ref, varnames["grid", "outname"]),
    fpc = paste0(path_ref, varnames["fpc", "outname"]),
    fpc_bft = paste0(path_ref, varnames["fpc_bft", "outname"]),
    cftfrac = paste0(path_ref, varnames["cftfrac", "outname"]),
    firec = paste0(path_ref, varnames["firec", "outname"]),
    npp = paste0(path_ref, varnames["npp", "outname"]),
    runoff = paste0(path_ref, varnames["runoff", "outname"]),
    transp = paste0(path_ref, varnames["transp", "outname"]),
    vegc = paste0(path_ref, varnames["vegc", "outname"]),
    firef = paste0(path_ref, varnames["firef", "outname"]),
    rh = paste0(path_ref, varnames["rh", "outname"]),
    harvestc = paste0(path_ref, varnames["harvestc", "outname"]),
    rharvestc = paste0(path_ref, varnames["rharvestc", "outname"]),
    pft_harvestc = paste0(path_ref, varnames["pft_harvest", "outname"]),
    pft_rharvestc = paste0(path_ref, varnames["pft_rharvest", "outname"]),
    evap = paste0(path_ref, varnames["evap", "outname"]),
    interc = paste0(path_ref, varnames["interc", "outname"]),
    discharge = paste0(path_ref, varnames["discharge", "outname"]),
    soilc = paste0(path_ref, varnames["soilc", "outname"]),
    litc = paste0(path_ref, varnames["litc", "outname"]),
    swc = paste0(path_ref, varnames["swc", "outname"]),
    swc_vol = paste0(path_ref, varnames["swc_vol", "outname"]),
    vegn = paste0(path_ref, varnames["vegn", "outname"]),
    soilnh4 = paste0(path_ref, varnames["soilnh4", "outname"]),
    soilno3 = paste0(path_ref, varnames["soilno3", "outname"]),
    leaching = paste0(path_ref, varnames["leaching", "outname"]),
    n2o_denit = paste0(path_ref, varnames["n2o_denit", "outname"]),
    n2o_nit = paste0(path_ref, varnames["n2o_nit", "outname"]),
    n2_emis = paste0(path_ref, varnames["n2_emis", "outname"]),
    bnf = paste0(path_ref, varnames["bnf", "outname"]),
    n_volatilization = paste0(path_ref, varnames["n_volatilization", "outname"]),
    gpp = paste0(path_ref, varnames["gpp", "outname"]),
    res_storage = paste0(path_ref, varnames["res_storage", "outname"]),
    lakevol = paste0(path_ref, varnames["lakevol", "outname"]),
    ndepos = paste0(path_ref, varnames["ndepos", "outname"]),
    rd = paste0(path_ref, varnames["rd", "outname"]),
    prec = paste0(path_ref, varnames["prec", "outname"]),
    terr_area = paste0(path_ref, varnames["terr_area", "outname"]),
    irrig = paste0(path_ref, varnames["irrig", "outname"]),
    nfert_agr = paste0(path_ref, varnames["nfert_agr", "outname"]),
    nmanure_agr = paste0(path_ref, varnames["nmanure_agr", "outname"]),
    ndepos = paste0(path_ref, varnames["ndepos", "outname"]),
    firen = paste0(path_ref, varnames["firen", "outname"]),
    harvestn = paste0(path_ref, varnames["harvestn", "outname"]),
    irrig_stor = paste0(path_ref, varnames["irrig_stor", "outname"]),
    rivervol = paste0(path_ref, varnames["rivervol", "outname"]),
    rootmoist = paste0(path_ref, varnames["rootmoist", "outname"])
  )

  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)) {
      print(paste("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 EcoRisk
    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
    )
    # 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))
  )
  for (y in 1:slices) {
    print(paste0("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))) {
    print(paste0("Saving EcoRisk data to: ", save_ecorisk))
    save(ecorisk, file = save_ecorisk)
  }
  #
  ###
  return(ecorisk)
}

#' Calculate the ecosystem change metric EcoRisk between 2 sets of states
#'
#' 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 1: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
#'
#' 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)
#'
#' @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) {
  file_type <- tools::file_ext(files_reference$grid)

  if (file_type %in% c("json", "clm")) {
    # read grid
    grid <- lpjmlkit::read_io(
      files_reference$grid
    )
    # calculate cell area
    cell_area <- drop(lpjmlkit::read_io(
      filename = files_reference$terr_area
      )$data) # 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)
    print("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))

    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))

    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))

    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))
    } 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))
    } 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))

    #### 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 1: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 <- 1:nstate_dimensions
    index <- 1
    # iterate over main classes (carbon pools, water fluxes ...)
    for (c in 1:nclasses) {
      classe <- metric_files$metric$ecorisk_nitrogen$metric_class[[c]]
      nsubclasses <- length(classe)
      # iterate over subclasses (vegetation carbon, soil water ...)
      for (s in 1:nsubclasses) {
        subclass <- classe[s]
        class_names[index] <- names(subclass)
        vars <- split_sign(unlist(subclass))
        for (v in 1: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)
            print(paste("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()
          } 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)
            print(paste("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()
          } 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))) {
    print(paste0("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