ecorisk.R

#' @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_mode = FALSE,
    suppress_warnings = 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"
    )
  } 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] +
        rowSums(cft_scen[, c(16, 32)]) * (1 - rowSums(bft_scen[, c(1:4, 7:10)]))
    )

    barren_area_scen[barren_area_scen < 0] <- 0

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

    # natural tree area fractions scaled by total natural frac
    tree_area_scen[, 1:7] <- (
      fpc_scen[, 2:8] * fpc_scen[, 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_scen[, 8:9] <- (
      cft_scen[, 16] * bft_scen[, 1:2] / rowSums(bft_scen[, 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_scen[, 10:11] <- (
      cft_scen[, 32] * bft_scen[, 7:8] / rowSums(bft_scen[, c(7, 8, 10)])
    )
    grass_area_scen <- array(0, dim = c(ncells, 20))

    # natural grass
    grass_area_scen[, 1:2] <- fpc_scen[, 9:10] * fpc_scen[, 1]

    # crops
    grass_area_scen[, 3:15] <- cft_scen[, 1:13] + cft_scen[, 17:29]

    # managed grass rf
    grass_area_scen[, 16] <- cft_scen[, 14]

    # managed grass irr
    grass_area_scen[, 17] <- cft_scen[, 30]

    # bioenergy grass
    grass_area_scen[, 18] <- cft_scen[, 15] + cft_scen[, 31]

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

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

    # evergreenness, needleleavedness, tropicalness, borealness, naturalness
    tree_attributes <- matrix(
      c(
        c(1, 0, 1, 0, 1), # 1 TrBE
        c(0, 0, 1, 0, 1), # 2 TrBR
        c(1, 1, 0, 0, 1), # 3 TeNE
        c(1, 0, 0, 0, 1), # 4 TeBE
        c(0, 0, 0, 0, 1), # 5 TeBS
        c(1, 1, 0, 1, 1), # 6 BoNE
        c(0, 0.25, 0, 1, 1), # 7 BoS (including larchs)
        c(1, 0, 1, 0, 0), # 8 TrBi tropical bioenergy rainfed
        c(0, 0, 0, 0, 0), # 9 TeBi temperate bioenergy rainfed
        c(1, 0, 1, 0, 0), # 10 TrBi tropical bioenergy irrigated
        c(0, 0, 0, 0, 0) # 11 TeBi temperate bioenergy irrigated
      ),
      nrow = 11,
      byrow = TRUE
    )

    if (weighting == "equal") {
      tree_weights <- c(0.2, 0.2, 0.2, 0.2, 0.2)

      # changed compared to Sebastian Ostberg's method
    } else if (weighting == "new") {
      tree_weights <- c(0.2, 0.2, 0.3, 0.3, 0.3) / 1.3

      # Sebastian's method (no downscaling to weightsum 1)
    } else if (weighting == "old") {
      tree_weights <- c(0.2, 0.2, 0.3, 0.3, 0.3)
    } else {
      stop("Unknown method of weighting.")
    }

    grass_attributes <- array(0, dim = c(ncells, 20, 2))
    # 1 C3grass
    # 2 C4grass
    # 3 TemperateCereals
    # 4 Rice
    # 5 Maize
    # 6 TropicalCereals
    # 7 Pulses
    # 8 TemperateRoots
    # 9 TropicalRoots
    # 10 Sunflower
    # 11 Soybean
    # 12 Groundnut
    # 13 Rapeseed
    # 14 Sugarcane
    # 15 Others
    # 16 Managed grass rainfed
    # 17 Managed grass irrigated
    # 18 Bioenergy grass
    # 19 Grass under rainfed Bioenergy trees
    # 20 Grass under irrigated Bioenergy trees

    # tropicalness
    grass_attributes[, , 1] <- rep(
      c(0, 1, 0, 1, 1, 1, 0.5, 0, 1, 0.5, 1, 1, 0.5, 1, 0.5, NA, NA, 1, NA, NA),
      each = ncells
    )

    # naturalness
    grass_attributes[, , 2] <- rep(
      c(1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
      each = ncells
    )

    # dynamic share of tropicalness for rf/irr grasslands taken from ratio of
    #   bioenergy grasses
    dyn_grass_attributes <- cbind(
      bft_scen[, 6] / rowSums(bft_scen[, 5:6]),
      bft_scen[, 12] / rowSums(bft_scen[, 11:12])
    )

    dyn_grass_attributes[!is.finite(dyn_grass_attributes)] <- 0

    # managed grass rf/irr
    grass_attributes[, 16:17, 1] <- dyn_grass_attributes

    # grass under biotrees rf/irr (taken from managed grass)
    grass_attributes[, 19:20, 1] <- dyn_grass_attributes

    if (weighting == "equal") {
      grass_weights <- c(0.2, 0.2)

      # changed compared to Sebastian Ostberg's method
    } else if (weighting == "new") {
      grass_weights <- c(0.5, 0.5)

      # Sebastian Ostbergs's method (no downscaling to weightsum 1)
    } else if (weighting == "old") {
      grass_weights <- c(0.3, 0.3)
    } else {
      stop("Unknown method of weighting.")
    }
  } else if (npfts == 11) {
    # barren = 1 - crop area - natural vegetation area +
    #   barren under bioenergy trees
    barren_area_ref <- (
      1 - rowSums(cft_ref) -
        rowSums(fpc_ref[, 2:12]) * fpc_ref[, 1] +
        rowSums(cft_ref[, c(16, 32)]) * (1 - rowSums(bft_ref[, c(4:9, 13:18)]))
    )

    barren_area_ref[barren_area_ref < 0] <- 0

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

    # natural tree area fractions scaled by total natural frac
    tree_area_ref[, 1:8] <- fpc_ref[, 2:9] * 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[, 9:10] <- (
      cft_ref[, 16] * bft_ref[, 7:8] / rowSums(bft_ref[, 4:8])
    )

    # 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[, 11:12] <- (
      cft_ref[, 32] * bft_ref[, 16:17] / rowSums(bft_ref[, 13:17])
    )

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

    # natural grass
    grass_area_ref[, 1:3] <- fpc_ref[, 10:12] * fpc_ref[, 1]

    # crops
    grass_area_ref[, 4:16] <- cft_ref[, 1:13] + cft_ref[, 17:29]

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

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

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

    # fraction of rainfed grass under bioenergy trees
    grass_area_ref[, 20] <- (
      cft_ref[, 16] * rowSums(bft_ref[, 4:6]) / rowSums(bft_ref[, 4:8])
    )

    # fraction of irrigated grass under bioenergy trees
    grass_area_ref[, 21] <- (
      cft_ref[, 32] * rowSums(bft_ref[, 13:15]) / rowSums(bft_ref[, 13:17])
    )

    # barren = 1 - crop area - natural vegetation area +
    #   barren under bioenergy trees
    barren_area_scen <- (
      1 - rowSums(cft_scen) -
        rowSums(fpc_scen[, 2:12]) * fpc_scen[, 1] +
        rowSums(cft_scen[, c(16, 32)]) *
          (1 - rowSums(bft_scen[, c(4:9, 13:18)]))
    )

    barren_area_scen[barren_area_scen < 0] <- 0

    tree_area_scen <- array(0, dim = c(ncells, 12))

    # natural tree area fractions scaled by total natural frac
    tree_area_scen[, 1:8] <- fpc_scen[, 2:9] * fpc_scen[, 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_scen[, 9:10] <- (
      cft_scen[, 16] * bft_scen[, 7:8] / rowSums(bft_scen[, 4:8])
    )

    # 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_scen[, 11:12] <- (
      cft_scen[, 32] * bft_scen[, 16:17] / rowSums(bft_scen[, 13:17])
    )

    grass_area_scen <- array(0, dim = c(ncells, 21))

    # natural grass
    grass_area_scen[, 1:3] <- fpc_scen[, 10:12] * fpc_scen[, 1]

    # crops
    grass_area_scen[, 4:16] <- cft_scen[, 1:13] + cft_scen[, 17:29]

    # managed grass rf
    grass_area_scen[, 17] <- cft_scen[, 14]

    # managed grass irr
    grass_area_scen[, 18] <- cft_scen[, 30]

    # bioenergy grass
    grass_area_scen[, 19] <- cft_scen[, 15] + cft_scen[, 31]

    # fraction of rainfed grass under bioenergy trees
    grass_area_scen[, 20] <- (
      cft_scen[, 16] * rowSums(bft_scen[, 4:6]) / rowSums(bft_scen[, 4:8])
    )

    # fraction of irrigated grass under bioenergy trees
    grass_area_scen[, 21] <- (
      cft_scen[, 32] * rowSums(bft_scen[, 13:15]) / rowSums(bft_scen[, 13:17])
    )

    # evergreenness, needleleavedness, tropicalness, borealness, naturalness
    tree_attributes <- matrix(
      c(
        c(1, 0, 1, 0, 1), # 1 TrBE
        c(0, 0, 1, 0, 1), # 2 TrBR
        c(1, 1, 0, 0, 1), # 3 TeNE
        c(1, 0, 0, 0, 1), # 4 TeBE
        c(0, 0, 0, 0, 1), # 5 TeBS
        c(1, 1, 0, 1, 1), # 6 BoNE
        c(0, 0, 0, 1, 1), # 7 BoBS
        c(0, 1, 0, 1, 1), # 8 BoNS
        c(1, 0, 1, 0, 0), # 9 TrBi tropical bioenergy rainfed
        c(0, 0, 0, 0, 0), # 10 TeBi temperate bioenergy rainfed
        c(1, 0, 1, 0, 0), # 11 TrBi tropical bioenergy irrigated
        c(0, 0, 0, 0, 0) # 12 TeBi temperate bioenergy irrigated
      ),
      nrow = 12,
      byrow = TRUE
    )

    if (weighting == "equal") {
      tree_weights <- c(0.2, 0.2, 0.2, 0.2, 0.2)

      # changed compared to Sebastian Ostberg's method
    } else if (weighting == "new") {
      tree_weights <- c(0.2, 0.2, 0.3, 0.3, 0.3) / 1.3

      # Sebastian Ostberg's method (no downscaling to weightsum 1)
    } else if (weighting == "old") {
      tree_weights <- c(0.2, 0.2, 0.3, 0.3, 0.3)
    } else {
      stop("Unknown method of weighting.")
    }

    grass_attributes <- array(0, dim = c(ncells, 21, 3))
    # 1 C4grass tropic
    # 2 C3grass temperate
    # 3 C3grass polar
    # 4 TemperateCereals
    # 5 Rice
    # 6 Maize
    # 7 TropicalCereals
    # 8 Pulses
    # 9 TemperateRoots
    # 10 TropicalRoots
    # 11 Sunflower
    # 12 Soybean
    # 13 Groundnut
    # 14 Rapeseed
    # 15 Sugarcane
    # 16 Others
    # 17 Managed grass rainfed
    # 18 Managed grass irrigated
    # 19 Bioenergy grass
    # 20 Grass under rainfed Bioenergy trees
    # 21 Grass under irrigated Bioenergy trees

    # tropicalness
    grass_attributes[, , 1] <- rep(
      c(
        1, 0, 0, 0, 1, 1, 1, 0.5, 0, 1, 0.5,
        1, 1, 0.5, 1, 0.5, NA, NA, 1, NA, NA
      ),
      each = ncells
    )

    # borealness
    grass_attributes[, , 2] <- rep(
      c(0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, NA, NA, 0, NA, NA),
      each = ncells
    )

    # naturalness
    grass_attributes[, , 3] <- rep(
      c(1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),
      each = ncells
    )

    # dynamic share of tropicalness for grass under irr biotrees
    dyn_trop_grass_attributes <- cbind(
      # dynamic share of tropicalness for rf grasslands
      bft_scen[, 1] / rowSums(bft_scen[, 1:3]),
      # dynamic share of tropicalness for irr grasslands
      bft_scen[, 10] / rowSums(bft_scen[, 10:12]),
      # dynamic share of tropicalness for grass under rf biotrees
      bft_scen[, 4] / rowSums(bft_scen[, 4:6]),
      bft_scen[, 13] / rowSums(bft_scen[, 13:15])
    )

    dyn_trop_grass_attributes[!is.finite(dyn_trop_grass_attributes)] <- 0

    # managed grass rf/irr, grass under biotrees rf/irr
    grass_attributes[, c(17, 18, 20, 21), 1] <- dyn_trop_grass_attributes

    # dynamic share of borealness for grass under irr biotrees
    dyn_boreal_grass_attributes <- cbind(
      # dynamic share of borealness for rf grasslands
      bft_scen[, 3] / rowSums(bft_scen[, 1:3]),
      # dynamic share of borealness for irr grasslands
      bft_scen[, 12] / rowSums(bft_scen[, 10:12]),
      # dynamic share of borealness for grass under rf biotrees
      bft_scen[, 6] / rowSums(bft_scen[, 4:6]),
      bft_scen[, 15] / rowSums(bft_scen[, 13:15])
    )

    dyn_boreal_grass_attributes[!is.finite(dyn_boreal_grass_attributes)] <- 0

    # managed grass rf/irr, grass under biotrees rf/irr
    grass_attributes[, c(17, 18, 20, 21), 2] <- dyn_boreal_grass_attributes

    if (weighting == "equal") {
      grass_weights <- c(0.2, 0.2, 0.2)
    } else if (weighting == "old" || weighting == "new") {
      grass_weights <- c(0.3333333, 0.3333333, 0.3333333)
    } else {
      stop("Unknown method of weighting.")
    }
  } else {
    stop("Unknown number of pfts.")
  }

  # compute vegetation_structure_change
  barren_v <- fBasics::rowMins(cbind(barren_area_ref, barren_area_scen))

  trees_v <- fBasics::rowMins(
    cbind(
      rowSums(tree_area_ref, na.rm = TRUE),
      rowSums(tree_area_scen, na.rm = TRUE)
    )
  )

  grass_v <- fBasics::rowMins(
    cbind(
      rowSums(grass_area_ref, na.rm = TRUE),
      rowSums(grass_area_scen, na.rm = TRUE)
    )
  )

  inner_sum_trees <- (
    # evergreenness
    abs(
      rowSums(tree_area_ref[, ] *
        rep(tree_attributes[, 1], each = ncells), na.rm = TRUE) -
        rowSums(tree_area_scen[, ] *
          rep(tree_attributes[, 1], each = ncells), na.rm = TRUE)
    ) * tree_weights[1] +
      # needleleavedness
      abs(
        rowSums(tree_area_ref[, ] *
          rep(tree_attributes[, 2], each = ncells), na.rm = TRUE) -
          rowSums(tree_area_scen[, ] *
            rep(tree_attributes[, 2], each = ncells), na.rm = TRUE)
      ) * tree_weights[2] +
      # tropicalness
      abs(
        rowSums(tree_area_ref[, ] *
          rep(tree_attributes[, 3], each = ncells), na.rm = TRUE) -
          rowSums(tree_area_scen[, ] *
            rep(tree_attributes[, 3], each = ncells), na.rm = TRUE)
      ) * tree_weights[3] +
      # borealness
      abs(
        rowSums(tree_area_ref[, ] *
          rep(tree_attributes[, 4], each = ncells), na.rm = TRUE) -
          rowSums(tree_area_scen[, ] *
            rep(tree_attributes[, 4], each = ncells), na.rm = TRUE)
      ) * tree_weights[4] +
      # naturalness
      abs(
        rowSums(tree_area_ref[, ] *
          rep(tree_attributes[, 5], each = ncells), na.rm = TRUE) -
          rowSums(tree_area_scen[, ] *
            rep(tree_attributes[, 5], each = ncells), na.rm = TRUE)
      ) * tree_weights[5]
  )

  if (npfts == 9) {
    inner_sum_grasses <- (
      # tropicalness
      abs(
        rowSums(grass_area_ref[, ] * grass_attributes[, , 1], na.rm = TRUE) -
          rowSums(grass_area_scen[, ] * grass_attributes[, , 1], na.rm = TRUE)
      ) * grass_weights[1] +
        # naturalness
        abs(
          rowSums(grass_area_ref[, ] * grass_attributes[, , 2], na.rm = TRUE) -
            rowSums(grass_area_scen[, ] * grass_attributes[, , 2], na.rm = TRUE)
        ) * grass_weights[2]
    )
  } else if (npfts == 11) {
    inner_sum_grasses <- (
      # tropicalness
      abs(
        rowSums(grass_area_ref[, ] * grass_attributes[, , 1], na.rm = TRUE) -
          rowSums(grass_area_scen[, ] * grass_attributes[, , 1], na.rm = TRUE)
      ) * grass_weights[1] +
        # borealness
        abs(
          rowSums(grass_area_ref[, ] * grass_attributes[, , 2], na.rm = TRUE) -
            rowSums(grass_area_scen[, ] * grass_attributes[, , 2], na.rm = TRUE)
        ) * grass_weights[2] +
        # naturalness
        abs(
          rowSums(grass_area_ref[, ] * grass_attributes[, , 3], na.rm = TRUE) -
            rowSums(grass_area_scen[, ] * grass_attributes[, , 3], na.rm = TRUE)
        ) * grass_weights[3]
    )
  } else {
    stop("Unknown number of pfts.")
  }

  vegetation_structure_change <- (
    1 - barren_v -
      trees_v * (1 - inner_sum_trees) -
      grass_v * (1 - inner_sum_grasses)
  )

  vegetation_structure_change[vegetation_structure_change < 0] <- 0

  vegetation_structure_change[!is.finite(vegetation_structure_change)] <- 0

  return(vegetation_structure_change)
}


################# further EcoRisk utility functions ##################

t_sigmoid_trafo <- function(x) {
  return(-1 / exp(3) + (1 + 1 / exp(3)) / (1 + exp(-6 * (x - 0.5))))
}


balance <- function(v1, v2) {
  return(1 - sum(v1 * v2) / (sqrt(sum(v1 * v1)) * sqrt(sum(v2 * v2))))
}


std_cellwise <- function(a) {
  return(apply(a, 1, stats::sd))
}


global_yearly_weighted_mean <- function(a, cell_area) {
  # a is matrix with dim=c(cells,years)
  # cell_area the corresponding cell_area array with dim=c(cells)
  return(
    sum(a * cell_area, na.rm = TRUE) /
      (length(a[1, ]) * sum(cell_area, na.rm = TRUE))
  )
}


globally_weighted_mean_foreach_var <- function(x, cell_area) { # nolint
  # x is matrix with dim=c(ncells,vars)
  # if vars==1 inflate x
  le <- length(x)
  if (le == length(cell_area)) dim(x) <- c(le, 1)

  # cell_area the corresponding cell_area array to x with dim=c(ncells)
  return(colSums(x * cell_area, na.rm = TRUE) / sum(cell_area, na.rm = TRUE))
}


s_change_to_var_ratio <- function(x, s) {
  return(1 / (1 + exp(-4 * (x / s - 2))))
}

aggregate_ecorisk_data_to_lower_resolution <- function(ecorisk_data_file,
                                                       new_resolution) {
  ecorisk_data_file <- "/p/projects/open/Fabian/Metrics/data/ecorisk_202311_data.RData"
  aggregation_factor <- 2 # 2 means from 0.5x0.5 to 1x1
  load(ecorisk_data_file) # load ecorisk data objects:
  # scenario data: state_scen,cft_scen,bft_scen,fpc_scen
  # reference data: state_ref,cft_ref,bft_ref,fpc_ref,
  # grid data: cell_area,lat,lon
  str(state_ref)

  state_scen <-
    state_ref <-
    fpc_scen <- fpc_ref
  bft_scen <- bft_ref
  cft_scen <- cft_ref

  di_state <- dim(state_scen)
  di_fpc <- dim(fpc_scen)
  di_bft <- dim(bft_scen)
  di_cft <- dim(cft_scen)
}


#' based on Heyder 2011 eq. 6-9; epsilon case handling from code
#'   by Sebastian Ostberg (not documented in papers)
#' @param ref mean reference state vector of dimension c(ncells,variables)
#' @param scen mean scenario state vector of dimension c(ncells,variables)
#' @param epsilon threshold for variables to be treated as 0
#'
#' @returns the length of the difference vector for each cell
state_diff_local <- function(ref, scen, epsilon = 10^-4) {
  # Ostberg code: case change_metric_lu_comparison_jun2013.c
  di <- dim(ref)
  # generally normalize the scenario state vector by the reference state
  s_scen <- scen / ref
  s_ref <- array(1, dim = di) # initialize

  # for variables in places, where ref is small (<epsilon),
  #   but scen larger (Ostberg code, line 798)
  # Sebastian set back scenario and reference vector, to keep the unscaled
  #   values (Ostberg code, line 804)
  cells_ref0 <- abs(ref) < epsilon & abs(scen) > epsilon
  s_scen[cells_ref0] <- scen[cells_ref0]
  s_ref[cells_ref0] <- ref[cells_ref0]

  # for variables in places, where ref and scen are small (<epsilon),
  #   return 0 (both are 1, difference is 0) (Ostberg code, line 809)
  s_scen[abs(ref) < epsilon & abs(scen) < epsilon] <- 1 # no change

  # normalize both state vectors by the sqrt(amount of state variables) to
  #   ensure length(s_ref)==1 (this is part of the weighting in the Ostberg
  #   code)
  s_ref <- s_ref / sqrt(di[2])
  s_scen <- s_scen / sqrt(di[2])

  # length of the local difference vector s_scen (sl2) - s_ref (sl1)
  return(sqrt(rowSums((s_scen - s_ref) * (s_scen - s_ref))))
}



#' c based on Heyder 2011 eq. 10-13
#'
#' @param ref mean reference state vector of dimension c(ncells,variables)
#' @param scen mean scenario state vector of dimension c(ncells,variables)
#' @param cell_area area of each cell as a vector of dim=c(ncells)
#' @param epsilon threshold for variables to be treated as 0
#'
#' @returns the length of the difference vector for each cell
state_diff_global <- function(ref, scen, cell_area, epsilon = 10^-4) {
  di <- dim(ref)
  ncells <- di[1]
  global_mean_ref <- globally_weighted_mean_foreach_var(ref, cell_area)
  global_mean_scen <- globally_weighted_mean_foreach_var(scen, cell_area)

  # if global mean state in ref period is 0 (e.g. for landuse vars in pnv run?)
  # take the mean scen state instead
  cells_ref0 <- abs(global_mean_ref) < epsilon & abs(global_mean_scen) > epsilon
  global_mean_ref[cells_ref0] <- global_mean_scen[cells_ref0]
  # if both are 0 take 1, then the division is defined but 0 - 0 leads
  # to no change, which is what EcoRisk should show
  cells_both0 <- (
    abs(global_mean_ref) < epsilon & abs(global_mean_scen) < epsilon
  )
  global_mean_ref[cells_both0] <- 1

  norm <- rep(global_mean_ref, each = ncells)
  dim(norm) <- dim(ref)
  s_scen <- scen / norm
  s_ref <- ref / norm

  # normalize both state vectors by the sqrt(amount of state variables) to
  #   ensure length(s_ref)==1
  # (this is part of the weighting in the Ostberg code)
  s_ref <- s_ref / sqrt(di[2])
  s_scen <- s_scen / sqrt(di[2])

  # length of the difference vector s_scen (sl2) - s_ref (sl1) for each cell
  return(sqrt(rowSums((s_scen - s_ref) * (s_scen - s_ref))))
}

calc_component <- function(ref, scen, local, cell_area) {
  di <- dim(ref)
  ncells <- di[1]
  nyears <- di[2]

  # calc mean ref and scen state
  if (length(di) == 2) {
    dim(ref) <- c(di, 1)
    dim(scen) <- c(di, 1)
  }
  ref_mean <- apply(ref, c(1, 3), mean)
  scen_mean <- apply(scen, c(1, 3), mean)


  if (local) {
    x <- t_sigmoid_trafo(state_diff_local(ref = ref_mean, scen = scen_mean))
  } else {
    x <- t_sigmoid_trafo(
      state_diff_global(ref = ref_mean, scen = scen_mean, cell_area = cell_area)
    )
  }
  # calculation of the change-to-variability ratio in my view is mathematically
  #   not correctly described in Heyder and Ostberg
  # - the way I understand it: recalculate the c/g/b value for each year of the
  #   ref period compared to the mean of the ref period as "scenario" and then
  #   calc the variability (sd) of that
  sigma_x_ref_list <- array(0, dim = c(ncells, nyears))
  for (i in seq_len(nyears)) {
    if (local) {
      sigma_x_ref_list[, i] <- t_sigmoid_trafo(
        state_diff_local(ref = ref_mean, scen = ref[, i, ])
      )
    } else {
      sigma_x_ref_list[, i] <- t_sigmoid_trafo(
        state_diff_global(
          ref = ref_mean,
          scen = ref[, i, ],
          cell_area = cell_area
        )
      )
    }
  }