ecorisk.R

  sigma_x_ref <- apply(sigma_x_ref_list, 1, stats::sd)
  c2vr <- s_change_to_var_ratio(x, sigma_x_ref)
  return(
    list(
      full = x * c2vr,
      value = x,
      var = c2vr
    )
  )
}


balance_shift <- function(ref, scen, epsilon = 10^-4) {
  # param ref with dimension c(ncells,nvars)
  # param scen with dimension c(ncells,nvars)

  if (is.null(dim(ref))) dim(ref) <- c(length(ref), 1)


  # first normalize as for local change
  s_scen <- scen / ref
  s_ref <- array(1, dim = dim(ref))

  # for variables in places, where ref is small (<epsilon), but scen larger
  # (Ostberg code, line 798/vector length calc in line 837)
  # set back scenario vector, to keep the unscaled values (Ostberg code,
  #   line 805)
  s_scen[abs(ref) < epsilon & abs(scen) > epsilon] <- (
    scen[abs(ref) < epsilon & abs(scen) > epsilon]
  )

  # for variables in places, where ref and scen are small (<epsilon),
  # set scen to 1 (both are 1, difference is 0 -- no change) (Ostberg code,
  #   line 809)

  # results in no change
  s_scen[abs(ref) < epsilon & abs(scen) < epsilon] <- 1
  # absa(_ts) in Sebastians Ostberg's code
  abs_ref <- sqrt(rowSums(s_ref * s_ref))
  # absb(_ts) in Sebastian Ostberg's code
  abs_scen <- sqrt(rowSums(s_scen * s_scen))

  # =1-angle_ts
  b1 <- 1 - (rowSums(s_ref * s_scen) / abs_ref / abs_scen) # todo: all -> 0

  # restrain to the maximum range for the acos function
  b1[b1 < 0] <- 0
  b1[b1 > 2] <- 2
  angle <- acos(1 - b1) * 360 / 2 / pi
  angle[b1 == 1] <- 0
  b <- b1 * 2
  b[angle > 60] <- 1

  return(b)
}


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

  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)


  b <- balance_shift(ref = ref_mean, scen = scen_mean)
  # 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_b_ref_list <- array(0, dim = c(ncells, nyears))
  for (i in seq_len(nyears)) {
    sigma_b_ref_list[, i] <- balance_shift(ref = ref_mean, scen = ref[, i, ])
  }
  sigma_b_ref <- apply(sigma_b_ref_list, 1, stats::sd)
  c2vr <- s_change_to_var_ratio(b, sigma_b_ref)
  return(
    list(
      full = b * c2vr,
      value = b,
      var = c2vr
    )
  )
}


#' Create modified EcoRisk data file
#'
#' Function to create a modified EcoRisk data file where each reference cell is
#' compared to the average reference biome cell. The scenario period is
#' overwritten with the original reference period and all reference cells are
#' set to the average cell of the prescribed reference biome ref_biom
#'
#' @param data_file_in path to input data
#' @param data_file_out path to save modified data to
#' @param biome_classes_in biome classes object as returned from classify_biomes
#' @param ref_biom reference biome from biome classes that all cells should
#'        be compared to
#'
#' @export
replace_ref_data_with_average_ref_biome_cell <- function(
    # nolint
    data_file_in,
    data_file_out,
    biome_classes_in,
    ref_biom) {
  if (data_file_in == data_file_out) {
    stop(
      "Same file for input and output of data, would overwrite ",
      "original data. Aborting."
    )
  }

  load(data_file_in)

  ref_cells <- which(biome_classes_in$biome_id == ref_biom)

  # first set all scen vars to the ref vars # [1:64240, 1:30, 1:10]
  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)

  # now replace all ref cells with that of the mean ref biom cell
  # FS 2022-08-10: keeping the year-to-year variation
  if (length(ref_cells) == 1) {
    av_year_state <- state_scen[ref_cells, , ]
    fpc_ref <- rep(fpc_scen[ref_cells, , ], each = di_fpc[1])
    bft_ref <- rep(bft_scen[ref_cells, , ], each = di_bft[1])
    cft_ref <- rep(cft_scen[ref_cells, , ], each = di_cft[1])
  } else {
    av_year_state <- apply(state_scen[ref_cells, , ], c(2, 3), mean)
    fpc_ref <- rep(
      apply(fpc_scen[ref_cells, , ], c(2, 3), mean),
      each = di_fpc[1]
    )
    bft_ref <- rep(
      apply(bft_scen[ref_cells, , ], c(2, 3), mean),
      each = di_bft[1]
    )
    cft_ref <- rep(
      apply(cft_scen[ref_cells, , ], c(2, 3), mean),
      each = di_cft[1]
    )
  }
  state_ref <- rep(av_year_state, each = di_state[1])
  dim(state_ref) <- di_state
  dimnames(state_ref) <- dimnames(state_scen)


  # is the same for each year, thus for the mean just take one year
  # mean_state_ref <- rep(colMeans(av_year_state), each = di_state[1])
  # FS: mean states were removed from data file, removing also here

  dim(fpc_ref) <- di_fpc
  dimnames(fpc_ref) <- dimnames(fpc_scen)
  dim(bft_ref) <- di_bft
  dimnames(bft_ref) <- dimnames(bft_scen)
  dim(cft_ref) <- di_cft
  dimnames(cft_ref) <- dimnames(cft_scen)


  # and write out the modified data
  # save(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,file = data_file_out)
  save(state_ref,
    state_scen,
    fpc_ref,
    fpc_scen,
    bft_ref,
    bft_scen,
    cft_ref,
    cft_scen,
    lat,
    lon,
    cell_area,
    file = data_file_out
  )
}


#' Create modified EcoRisk data for crosstable
#'
#' Function to create a modified EcoRisk data file where for each biome
#' the average scenario cell is compared to the average scenario cell of all
#' other biomes. This can then be used to compute a crosstable with the average
#' difference between each of them as in the SI of Ostberg et al. 2013
#' (Critical impacts of global warming on land ecosystems)
#'
#' @param data_file_in path to input data
#' @param data_file_out path to save modified data to
#' @param biome_classes_in biome classes object as returned from classify_biomes
#' @param pick_cells pick one specific cell as representative for the biome
#'        instead of computing the average state
#' @param baseline_ref logical, use reference state as baseline?
#'        default is FALSE - use scenario state
#'
#' @return c2vr array to be used in the ecorisk call
#'
#' @export
ecorisk_cross_table <- function(data_file_in,
                                data_file_out,
                                biome_classes_in,
                                pick_cells = NULL,
                                baseline_ref = FALSE) {
  if (data_file_in == data_file_out) {
    stop(
      "Same file for input and output of data, would overwrite original data. ",
      "Aborting."
    )
  }
  load(data_file_in)

  # calculate ecorisk to get the variability
  ecorisk_c2vr <- drop(ecorisk_wrapper(
    path_ref = NULL,
    path_scen = NULL,
    read_saved_data = TRUE,
    nitrogen = TRUE,
    weighting = "equal",
    save_data = data_file_in,
    save_ecorisk = NULL,
    time_span_reference = as.character(1550:1579),
    time_span_scenario = as.character(1985:2014),
    dimensions_only_local = FALSE
  )$c2vr)

  if (baseline_ref) {
    # save reference state vectors, they contain relevant data (scen can go)
    state_sav <- state_ref
    fpc_sav <- fpc_ref
    bft_sav <- bft_ref
    cft_sav <- cft_ref
  } else {
    # save scenario state vectors, they contain relevant data (ref can go)
    state_sav <- state_scen
    fpc_sav <- fpc_scen
    bft_sav <- bft_scen
    cft_sav <- cft_scen
  }


  nbiomes <- max(biome_classes_in$biome_id) # by default 19
  state_ref <- array(0, dim = c(nbiomes, nbiomes, dim(state_sav)[2:3]))
  state_scen <- state_ref
  fpc_ref <- array(0, dim = c(nbiomes, nbiomes, dim(fpc_sav)[2:3]))
  fpc_scen <- fpc_ref
  bft_ref <- array(0, dim = c(nbiomes, nbiomes, dim(bft_sav)[2:3]))
  bft_scen <- bft_ref
  cft_ref <- array(0, dim = c(nbiomes, nbiomes, dim(cft_sav)[2:3]))
  cft_scen <- cft_ref
  c2vr <- array(0, dim = c(4, nbiomes))

  # now replace all ref cells with that of the mean ref biome cell
  for (b in sort(unique(biome_classes_in$biome_id))) {
    ref_cells <- which(biome_classes_in$biome_id == b)

    if (is.null(pick_cells)) {
      if (length(ref_cells) == 1) {
        # average over cells, keeping the average year-to-year variation
        av_state <- state_sav[ref_cells, , ]
        av_fpc <- fpc_sav[ref_cells, , ]
        av_bft <- bft_sav[ref_cells, , ]
        av_cft <- cft_sav[ref_cells, , ]
        av_c2vr <- ecorisk_c2vr[, ref_cells]
      } else {
        # average over cells, keeping the average year-to-year variation
        av_state <- apply(state_sav[ref_cells, , ], c(2, 3), mean)
        av_fpc <- apply(fpc_sav[ref_cells, , ], c(2, 3), mean)
        av_bft <- apply(bft_sav[ref_cells, , ], c(2, 3), mean)
        av_cft <- apply(cft_sav[ref_cells, , ], c(2, 3), mean)
        av_c2vr <- apply(ecorisk_c2vr[, ref_cells], c(1), mean)
      }
    } else { # use pick_cells
      av_state <- state_sav[pick_cells[b], , ]
      av_fpc <- fpc_sav[pick_cells[b], , ]
      av_bft <- bft_sav[pick_cells[b], , ]
      av_cft <- cft_sav[pick_cells[b], , ]
      av_c2vr <- ecorisk_c2vr[, pick_cells[b]]
    }

    state_ref[b, , , ] <- rep(av_state, each = nbiomes)
    state_scen[, b, , ] <- rep(av_state, each = nbiomes)

    mean_state_ref <- apply(state_ref, c(1, 3), mean)
    mean_state_scen <- apply(state_scen, c(1, 3), mean)

    fpc_ref[b, , , ] <- rep(av_fpc, each = nbiomes)
    fpc_scen[, b, , ] <- rep(av_fpc, each = nbiomes)

    bft_ref[b, , , ] <- rep(av_bft, each = nbiomes)
    bft_scen[, b, , ] <- rep(av_bft, each = nbiomes)

    cft_ref[b, , , ] <- rep(av_cft, each = nbiomes)
    cft_scen[, b, , ] <- rep(av_cft, each = nbiomes)
    c2vr[, b] <- av_c2vr
  }
  dim(state_ref) <- c(nbiomes * nbiomes, dim(state_sav)[2:3])
  dim(state_scen) <- c(nbiomes * nbiomes, dim(state_sav)[2:3])
  dim(fpc_ref) <- c(nbiomes * nbiomes, dim(fpc_sav)[2:3])
  dim(fpc_scen) <- c(nbiomes * nbiomes, dim(fpc_sav)[2:3])
  dim(bft_ref) <- c(nbiomes * nbiomes, dim(bft_sav)[2:3])
  dim(bft_scen) <- c(nbiomes * nbiomes, dim(bft_sav)[2:3])
  dim(cft_ref) <- c(nbiomes * nbiomes, dim(cft_sav)[2:3])
  dim(cft_scen) <- c(nbiomes * nbiomes, dim(cft_sav)[2:3])
  dimnames(state_ref) <- list(
    cell = as.character(seq_len(nbiomes * nbiomes)),
    year = dimnames(state_sav)$year,
    class = dimnames(state_sav)$class
  )
  dimnames(state_scen) <- dimnames(state_ref)
  dimnames(fpc_ref) <- list(
    cell = as.character(seq_len(nbiomes * nbiomes)),
    band = dimnames(fpc_sav)$band,
    year = dimnames(fpc_sav)$year
  )
  dimnames(fpc_scen) <- dimnames(fpc_ref)
  dimnames(bft_ref) <- list(
    cell = as.character(seq_len(nbiomes * nbiomes)),
    band = dimnames(bft_sav)$band,
    year = dimnames(bft_sav)$year
  )
  dimnames(bft_scen) <- dimnames(bft_ref)
  dimnames(cft_ref) <- list(
    cell = as.character(seq_len(nbiomes * nbiomes)),
    band = dimnames(cft_sav)$band,
    year = dimnames(cft_sav)$year
  )
  dimnames(cft_scen) <- dimnames(cft_ref)
  dimnames(c2vr) <- list(
    component = c("vs", "lc", "gi", "eb"),
    biome = biome_classes_in$biome_names
  )


  lat <- rep(0, nbiomes * nbiomes)
  lon <- rep(1, nbiomes * nbiomes)
  cell_area <- rep(2, nbiomes * nbiomes)

  # and write out the modified data
  save(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,
    file = data_file_out
  )
  return(c2vr)
}

#' Create modified EcoRisk data combining two time series
#'
#' Function to combine two EcoRisk data files (e.g. historic and futures) into
#' one file.
#'
#' @param hist_file path to input file with historic data
#' @param scen_file path to input file with scenario data
#' @param combined_file path to save modified data to
#'
#'
#' @export
ecorisk_combine_hist_and_scen_data <- function(hist_file, scen_file, combined_file) {
  load(hist_file)
  state_scen_hist <- state_scen
  fpc_scen_hist <- fpc_scen
  bft_scen_hist <- bft_scen
  cft_scen_hist <- cft_scen
  load(scen_file)
  state_scen_scen <- state_scen
  fpc_scen_scen <- fpc_scen
  bft_scen_scen <- bft_scen
  cft_scen_scen <- cft_scen

  dimnames_state_hist <- dimnames(state_scen_hist)
  dimnames_state_scen <- dimnames(state_scen_scen)
  dimnames_state_comb <- dimnames_state_hist
  dimnames_state_comb$year <- c(dimnames_state_hist$year, dimnames_state_scen$year)
  di_state_hist <- dim(state_scen_hist)
  di_state_scen <- dim(state_scen_scen)
  di_state_comb <- di_state_hist
  di_state_comb[2] <- di_state_hist[2] + di_state_scen[2]
  state_scen <- array(0, dim = di_state_comb)
  dimnames(state_scen) <- dimnames_state_comb
  state_scen[, dimnames_state_hist$year, ] <- state_scen_hist
  state_scen[, dimnames_state_scen$year, ] <- state_scen_scen

  dimnames_cft_hist <- dimnames(cft_scen_hist)
  dimnames_cft_scen <- dimnames(cft_scen_scen)
  dimnames_cft_comb <- dimnames_cft_hist
  dimnames_cft_comb$year <- c(dimnames_cft_hist$year, dimnames_cft_scen$year)
  di_cft_hist <- dim(cft_scen_hist)
  di_cft_scen <- dim(cft_scen_scen)
  di_cft_comb <- di_cft_hist
  di_cft_comb[3] <- di_cft_hist[3] + di_cft_scen[3]
  cft_scen <- array(0, dim = di_cft_comb)
  dimnames(cft_scen) <- dimnames_cft_comb
  cft_scen[, , dimnames_cft_hist$year] <- cft_scen_hist
  cft_scen[, , dimnames_cft_scen$year] <- cft_scen_scen

  dimnames_fpc_hist <- dimnames(fpc_scen_hist)
  dimnames_fpc_scen <- dimnames(fpc_scen_scen)
  dimnames_fpc_comb <- dimnames_fpc_hist
  dimnames_fpc_comb$year <- c(dimnames_fpc_hist$year, dimnames_fpc_scen$year)
  di_fpc_hist <- dim(fpc_scen_hist)
  di_fpc_scen <- dim(fpc_scen_scen)
  di_fpc_comb <- di_fpc_hist
  di_fpc_comb[3] <- di_fpc_hist[3] + di_fpc_scen[3]
  fpc_scen <- array(0, dim = di_fpc_comb)
  dimnames(fpc_scen) <- dimnames_fpc_comb
  fpc_scen[, , dimnames_fpc_hist$year] <- fpc_scen_hist
  fpc_scen[, , dimnames_fpc_scen$year] <- fpc_scen_scen

  dimnames_bft_hist <- dimnames(bft_scen_hist)
  dimnames_bft_scen <- dimnames(bft_scen_scen)
  dimnames_bft_comb <- dimnames_bft_hist
  dimnames_bft_comb$year <- c(dimnames_bft_hist$year, dimnames_bft_scen$year)
  di_bft_hist <- dim(bft_scen_hist)
  di_bft_scen <- dim(bft_scen_scen)
  di_bft_comb <- di_bft_hist
  di_bft_comb[3] <- di_bft_hist[3] + di_bft_scen[3]
  bft_scen <- array(0, dim = di_bft_comb)
  dimnames(bft_scen) <- dimnames_bft_comb
  bft_scen[, , dimnames_bft_hist$year] <- bft_scen_hist
  bft_scen[, , dimnames_bft_scen$year] <- bft_scen_scen

  if (file.exists(combined_file)) {
    message("Output file ", combined_file, " already exists. Aborting.")
  } else {
    message("Saving data to: ", combined_file)
    save(state_ref, state_scen, fpc_ref, fpc_scen, bft_ref, bft_scen,
      cft_ref, cft_scen, lat, lon, cell_area,
      file = combined_file
    )
  }
}

################# biome (dis-)aggregation functions ##################

#' Get biome names
#'
#' Returns biome names with variable length (abbreviated, short, or full)
#'
#' @param biome_name_length integer chose from 1,2,3 for abbreviated, short,
#'                        or full biome names
#'
#' @export
get_biome_names <- function(biome_name_length = 2) {
  biome_mapping <- utils::read.csv(
    file = system.file(
      "extdata",
      "biomes.csv",
      package = "biospheremetrics"
    ),
    sep = ";"
  )

  if (biome_name_length == 1) {
    biome_class_names <- biome_mapping$abbreviation
  } else if (biome_name_length == 2) {
    biome_class_names <- biome_mapping$short_name
  } else if (biome_name_length == 3) {
    biome_class_names <- biome_mapping$name
  } else {
    stop(
      "Value for parameter biome_name_length out of range 1,2,3 - ",
      "was given as: ",
      biome_name_length
    )
  }

  return(biome_class_names)
}


#' Averages EcoRisk values across regions
#'
#' Returns the average value across either 4 regions or all (19) biomes for
#' EcoRisk and each of the subcomponents for each
#'
#' @param data List object, of which every item should be disaggregated
#' @param biome_class biome class list object as returned by classify_biomes
#' @param type string controlling whether to return  minimum, mean, maximum
#'        ("minmeanmax") or Q10,Q50,Q90 ("quantile") - default: "quantile"
#' @param classes string for into how many regions should be disaggregated
#'        "4biomes" (tropics/temperate/boreal/arctic) or "allbiomes"
#'
#' @examples
#' \dontrun{
#' disaggregate_into_biomes(
#'   ecorisk = ecorisk,
#'   biome_class = biome_classes,
#'   type = "quantile", classes = "4biomes"
#' )
#' }
#' @export
disaggregate_into_biomes <- function(data, # nolint
                                     biome_class,
                                     type = "quantile",
                                     classes = "4biomes") {
  di <- dim(data[[1]])
  comp_names <- names(data)

  if (type == "minmeanmax") {
    type_names <- c("min", "mean", "max")
  } else if (type == "quantile") {
    type_names <- c("Q10", "Q50", "Q90")
  }

  if (length(di) > 1) {
    slices <- di[2]
  } else {
    slices <- 1
  }

  if (classes == "4biomes") {
    tropics <- c(1, 2, 9, 10, 11)
    temperate <- c(3, 4, 5, 12, 13, 14)
    boreal <- c(6, 7, 8)
    arctic <- c(15, 16)
    cell_list <- list(
      tropical_cells = which(biome_class$biome_id %in% tropics),
      temperate_cells = which(biome_class$biome_id %in% temperate),
      boreal_cells = which(biome_class$biome_id %in% boreal),
      arctic_cells = which(biome_class$biome_id %in% arctic)
    )
    nclasses <- 4
  } else if (classes == "allbiomes") {
    nclasses <- max(unique(biome_class$biome_id))
  } else {
    stop(
      "Unknown parameter classes: ",
      classes,
      ", should be either '4biomes' or 'allbiomes'"
    )
  }

  data_dims <- length(data)

  data_biomes <- array(0, dim = c(nclasses, data_dims, 3, slices))

  if (classes == "4biomes") { # aggregate to trop/temp/boreal/arctic
    for (s in seq_len(slices)) {
      for (b in seq_len(nclasses)) {
        for (cc in seq_len(data_dims)) {
          if (type == "minmeanmax") {
            data_biomes[b, cc, , s] <- c(
              min(data[[cc]][cell_list[[b]], s], na.rm = TRUE),
              mean(data[[cc]][cell_list[[b]], s], na.rm = TRUE),
              max(data[[cc]][cell_list[[b]], s], na.rm = TRUE)
            )
          } else if (type == "quantile") {
            data_biomes[b, cc, , s] <- c(
              stats::quantile(
                data[[cc]][cell_list[[b]], s],
                probs = c(0.1, 0.5, 0.9),
                na.rm = TRUE
              )
            )
          } else {
            stop(paste(
              "type", type,
              "unknown. please choose either 'quantile' or 'minmeanmax'"
            ))
          } # end if
        } # end for
      } # end for
    } # end for

    biome_names <- c("tropics", "temperate", "boreal", "arctic")
    dimnames(data_biomes) <- list(
      biome_names, comp_names,
      type_names, seq_len(slices)
    )
  } else if (classes == "allbiomes") { # calculate all biomes separately
    for (s in seq_len(slices)) {
      for (b in seq_len(nclasses)) {
        for (cc in seq_len(data_dims)) {
          if (type == "minmeanmax") {
            data_biomes[b, cc, , s] <-
              c(
                min(data[[cc]][which(biome_class$biome_id == b), s], na.rm = TRUE),
                mean(data[[cc]][which(biome_class$biome_id == b), s], na.rm = TRUE),
                max(data[[cc]][which(biome_class$biome_id == b), s], na.rm = TRUE)
              )
          } else if (type == "quantile") {
            data_biomes[b, cc, , s] <- c(
              stats::quantile(
                data[[cc]][which(biome_class$biome_id == b), s],
                probs = c(0.1, 0.5, 0.9),
                na.rm = TRUE
              )
            )
          } else {
            stop(paste(
              "type", type,
              "unknown. please choose either 'quantile' or 'minmeanmax'"
            ))
          } # end if
        } # end for
      } # end for
    } # end for

    biome_names <- biome_class$biome_names
    dimnames(data_biomes) <- list(
      biome_names,
      comp_names,
      type_names,
      seq_len(slices)
    )
  } else {
    stop(
      "Unknown parameter classes: ",
      classes,
      ", should be either '4biomes' or 'allbiomes'"
    )
  }
  return(drop(data_biomes))
}


#' Calculate ecorisk with each biomes average cell
#'
#' Function to calculate ecorisk with each biomes average cell
#' as a measure of internal variability
#'
#' @param biome_classes biome classes object as returned by classify biomes,
#'                      calculated for data_file_base
#' @param data_file_base base EcoRisk to compute differences with (only ref is
#'                      relevant)
#' @param intra_biome_distrib_file file to additionally write results to
#' @param create create new modified files, or read already existing ones?
#' @param res how finegrained the distribution should be (resolution)
#' @param nitrogen include nitrogen outputs (default: TRUE)
#' @param plotting whether plots for each biome should be created
#' @param plot_folder folder to plot into
#' @param time_span_reference suitable 30 year reference period (e.g.
#'                            c(1901,1930), c(1550,1579))
#' @param ecorisk_components integer. how many subcomponents does the
#'        ecorisk object have?

#' @return data object with distibution - dim: c(biomes,ecorisk_variables,bins)
#'
#' @export
calculate_within_biome_diffs <- function(biome_classes, # nolint
                                         data_file_base,
                                         intra_biome_distrib_file,
                                         create = FALSE,
                                         nitrogen = TRUE,
                                         res = 0.05,
                                         plotting = FALSE,
                                         plot_folder,
                                         time_span_reference,
                                         ecorisk_components = 13) {
  biomes_abbrv <- get_biome_names(1)
  # nbiomes, nEcoRiskvars, nHISTclasses
  intra_biome_distrib <- array(
    0,
    dim = c(length(biome_classes$biome_names), ecorisk_components, 1 / res)
  )

  # start
  for (b in sort(unique(biome_classes$biome_id))) {
    filebase <- strsplit(data_file_base, "_data.RData")[[1]]
    message(
      "Calculating differences with biome ", b, " (",
      biome_classes$biome_names[b], ")"
    )

    data_file <- paste0(
      filebase, "_compared_to_average_", biomes_abbrv[b], "_data.RData"
    )
    ecorisk_file <- paste0(
      filebase, "_compared_to_average_", biomes_abbrv[b], "_gamma.RData"
    )

    if (create) {
      replace_ref_data_with_average_ref_biome_cell(
        data_file_in = data_file_base,
        data_file_out = data_file,
        biome_classes_in = biome_classes,
        ref_biom = b
      )
      ecorisk <- ecorisk_wrapper(
        # does not need to be specified, as data is read from file
        path_ref = NULL,
        # does not need to be specified, as data is read from file
        path_scen = NULL,
        read_saved_data = TRUE,
        nitrogen = nitrogen,
        weighting = "equal",
        save_data = data_file,
        save_ecorisk = ecorisk_file,
        time_span_reference = time_span_reference,
        time_span_scenario = time_span_reference,
        dimensions_only_local = FALSE
      )
    } else {
      # contains ecorisk list object
      load(ecorisk_file)
    }

    # compute average values per focus biom
    ref_cells <- which(biome_classes$biome_id == b)
    for (v in seq_len(10)) {
      intra_biome_distrib[b, v, ] <- graphics::hist(
        ecorisk[[v]][ref_cells],
        breaks = seq(0, 1, res), plot = FALSE
      )$counts
      intra_biome_distrib[b, v, ] <- (
        intra_biome_distrib[b, v, ] / sum(intra_biome_distrib[b, v, ])
      )
    }

    if (plotting) {
      plot_ecorisk_map(
        file = paste0(
          plot_folder, "/compare_ecorisk_to_", biomes_abbrv[b], ".png"
        ),
        focus_biome = b, biome_classes = biome_classes$biome_id,
        ecorisk_object = ecorisk,
        plot_dimension = "ecorisk_total",
        title = biome_classes$biome_names[b],
        legendtitle = "",
        eps = FALSE,
        title_size = 2,
        leg_yes = TRUE
      )

      plot_ecorisk_map(
        file = paste0(
          plot_folder, "/compare_vegetation_structure_change_to_",
          biomes_abbrv[b], ".png" # nolint
        ),
        focus_biome = b, biome_classes = biome_classes$biome_id,
        ecorisk_object = ecorisk,
        plot_dimension = "vegetation_structure_change",
        title = biome_classes$biome_names[b],
        legendtitle = "",
        eps = FALSE,
        title_size = 2,
        leg_yes = TRUE
      )

      plot_ecorisk_map(
        file = paste0(
          plot_folder, "/compare_gi_to_", biomes_abbrv[b], ".png"
        ),
        focus_biome = b,
        biome_classes = biome_classes$biome_id,
        ecorisk_object = ecorisk,
        plot_dimension = "global_importance",
        title = biome_classes$biome_names[b],
        legendtitle = "",
        eps = FALSE,
        title_size = 2,
        leg_yes = TRUE
      )

      plot_ecorisk_map(
        file = paste0(
          plot_folder, "/compare_lc_to_", biomes_abbrv[b], ".png"
        ),
        focus_biome = b,
        biome_classes = biome_classes$biome_id,
        ecorisk_object = ecorisk,
        plot_dimension = "local_change",
        title = biome_classes$biome_names[b],
        legendtitle = "",
        eps = FALSE,
        title_size = 2,
        leg_yes = TRUE
      )

      plot_ecorisk_map(
        file = paste0(
          plot_folder, "/compare_eb_to_", biomes_abbrv[b], ".png"
        ),
        focus_biome = b,
        biome_classes = biome_classes$biome_id,
        ecorisk_object = ecorisk,
        plot_dimension = "ecosystem_balance",
        title = biome_classes$biome_names[b],
        legendtitle = "",
        eps = FALSE,
        title_size = 2,
        leg_yes = TRUE
      )
    } # end if plotting
  }
  ecorisk_dimensions <- names(ecorisk)

  dim(intra_biome_distrib) <- c(
    biome = length(biome_classes$biome_names),
    variable = ecorisk_components, bin = 1 / res
  )
  dimnames(intra_biome_distrib) <- list(
    biome = biomes_abbrv, variable = ecorisk_dimensions, bin = seq(res, 1, res)
  )
  save(intra_biome_distrib, file = intra_biome_distrib_file)

  return(intra_biome_distrib)
}