# written by Fabian Stenzel
# 2022-2023 - stenzel@pik-potsdam.de
################# BioCol calc functions ###################
#' Calculate BioCol based on file lists from a PNV run and LU run of LPJmL.
#' Do not use this function directly, unless you are instructed to do so, there
#' is a wrapper called calc_biocol() which is for use of endusers.
#'
#' Function to calculate BioCol based on a PNV run and LU run of LPJmL
#' @param files_scenario list with variable names and corresponding file paths
#' (character string) of the scenario LPJmL run. All needed files are
#' provided in XXX. E.g.: list(npp = "/temp/npp.bin.json")
#' @param files_baseline list with variable names and corresponding file paths
#' (character string) of the baseline LPJmL run. All needed files are
#' provided in XXX. E.g.: list(npp = "/temp/npp.bin.json"). If not
#' needed for the applied method, set to NULL.
#' @param files_reference list with npp file path (character string) of the
#' reference LPJmL run (usually Holocene/preindustrial).
#' E.g.: list(npp = "/temp/npp.bin.json"). If NULL uses baseline npp.
#' @param time_span_scenario time span to be used for the scenario run, defined
#' as a character vector, e.g. `as.character(1982:2011)` (required)
#' @param time_span_baseline time span to be used for the baseline run, defined
#' as a character vector, e.g. `as.character(1901:1930)`. Can differ in offset
#' and length from `time_span_scenario`! If `NULL` value of `time_span_scenario`
#' is used
#' @param time_span_reference time span to read reference npp from, using
#' index years 10:39 from potential npp input if set to NULL (default: NULL)
#' @param gridbased logical are pft outputs gridbased or pft-based?
#' @param read_saved_data flag whether to read previously saved data
#' instead of reading it in from output files (default FALSE)
#' @param save_data whether to save input data to file (default FALSE)
#' @param data_file file to save/read input data to/from (default NULL)
#' @param include_fire boolean include firec in calculation of BioCol?
#' (default TRUE)
#' @param external_fire instead of reading in firec for fire emissions, read in
#' this external firec file from a separate spitfire run with disabled
#' lighning. this will then include only human induced fires
#' (default FALSE)
#' @param external_wood_harvest include external wood harvest from LUH2_v2h
#' (default FALSE)
#' @param grass_scaling whether to scale pasture harvest according to
#' data given via grass_harvest_file (default FALSE)
#' @param npp_threshold lower threshold for npp (to mask out non-lu areas
#' according to Haberl et al. 2007). Below BioCol will be set to 0.
#' (default: 20 gC/m2)
#' @param epsilon minimum value for npp, below which it will be set to 0
#' @param grass_harvest_file file containing grazing data to rescale the
#' grassland harvests according to Herrero et al. 2013. File contains:
#' grazing_data list object with $name and $id of 29 world regions, and
#' $Herrero_2000_kgDM_by_region containing for each of these regions and
#' mapping_lpj67420_to_grazing_regions array with a mapping between 67420
#' LPJmL cells and the 29 regions
#' @param external_fire_file path to external file with human induced fire
#' fraction c(cell,month,year) since 1500
#' @param external_wood_harvest_file path to R-file containing processed
#' timeline of maps for LUH2_v2h woodharvest
#' @param suppress_warnings suppress warnings when reading files (default: TRUE)
#'
#' @return list data object containing BioCol and components as arrays:
#' biocol_overtime, biocol_overtime_abs, biocol_overtime_abs_frac_piref,
#' biocol_overtime_abs_frac, biocol_overtime_pos,
#' biocol_overtime_pos_frac_piref,biocol_overtime_pos_frac,
#' biocol_overtime_frac_piref, biocol_overtime_frac, npp_harv_overtime,
#' npp_luc_overtime,npp_act_overtime, npp_pot_overtime,npp_eco_overtime,
#' harvest_grasslands_overtime, harvest_bioenergy_overtime,
#' harvest_cft_overtime, rharvest_cft_overtime, fire_overtime,
#' timber_harvest_overtime, wood_harvest_overtime, biocol, biocol_frac,
#' npp, biocol_frac_piref, npp_potential, npp_ref, harvest_cft,
#' rharvest_cft, biocol_harvest, biocol_luc, lat, lon, cellarea
#'
#' @export
read_calc_biocol <- function(
# nolint
files_scenario,
files_baseline,
files_reference = NULL,
time_span_scenario,
time_span_baseline = NULL,
time_span_reference = NULL,
gridbased = TRUE,
read_saved_data = FALSE,
save_data = FALSE,
data_file = NULL,
include_fire = FALSE,
external_fire = FALSE,
external_wood_harvest = FALSE,
grass_scaling = FALSE,
npp_threshold = 20,
epsilon = 0.001, # gC/m2
grass_harvest_file = NULL,
external_fire_file = NULL,
external_wood_harvest_file = NULL,
suppress_warnings = TRUE) {
if (is.null(files_reference)) {
files_reference <- list(npp = baseline_npp_file)
}
if (is.null(time_span_baseline)) {
time_span_baseline <- time_span_scenario
}
if (is.null(time_span_reference)) {
time_span_reference <- time_span_scenario[3:12]
}
if (grass_scaling && !file.exists(grass_harvest_file)) {
stop(
paste0("Grass harvest scaling enabled, but grass_harvest_file \
does not exist in: ", grass_harvest_file)
)
}
if (external_wood_harvest && !file.exists(external_wood_harvest_file)) {
stop(
paste0("External wood harvest enabled, but external_wood_harvest_file \
does not exist in: ", external_wood_harvest_file)
)
}
if (external_fire && !file.exists(external_fire_file)) {
stop(
paste0("External fire fraction file enabled, but external_fire_file \
does not exist in: ", external_fire_file)
)
}
# reading required data
if (read_saved_data) {
if (file.exists(data_file)) {
message("Reading in data from previously saved data file")
load(data_file)
wood_harvest[is.na(wood_harvest)] <- 0
} else {
stop(
paste0(
"data_file: '",
data_file,
"' does not exist but is required since reading is set to FALSE."
)
)
}
if (save_data) {
save_data <- FALSE
message(
"Both read_saved_data and save_data have been set to TRUE. ",
"Overwriting with the same data does not make sense, saving ",
"disabled. "
)
}
} else {
message("Reading in data from outputs")
file_type <- tools::file_ext(files_baseline$grid)
if (file_type %in% c("json", "clm")) {
# read grid
grid <- lpjmlkit::read_io(
files_baseline$grid
) %>% suppressWarnings()
# calculate cell area
cellarea <- drop(lpjmlkit::read_io(
filename = files_baseline$terr_area
)$data) %>% suppressWarnings() # in m2
lat <- grid$data[, , 2]
lon <- grid$data[, , 1]
npp <- abind::adrop(lpjmlkit::read_io(
files_scenario$npp,
subset = list(year = as.character(time_span_scenario))
) %>%
lpjmlkit::transform(to = c("year_month_day")) %>%
lpjmlkit::as_array(aggregate = list(month = sum)) %>%
# drop() %>% suppressWarnings()
suppressWarnings(), drop = "band") # gC/m2
npp[npp < epsilon] <- 0
if (!is.null(files_reference)) {
npp_ref <- abind::adrop(lpjmlkit::read_io(
files_reference$npp,
subset = list(year = as.character(time_span_reference))
) %>%
lpjmlkit::transform(to = c("year_month_day")) %>%
lpjmlkit::as_array(aggregate = list(month = sum)) %>%
# drop() %>% suppressWarnings()
suppressWarnings(), drop = "band")
npp_ref[npp_ref < epsilon] <- 0
}
pftnpp <- lpjmlkit::read_io(
files_scenario$pft_npp,
subset = list(year = as.character(time_span_scenario))
) %>%
lpjmlkit::transform(to = c("year_month_day")) %>%
lpjmlkit::as_array(aggregate = list(month = sum)) %>%
suppressWarnings()
pftnpp[pftnpp < epsilon] <- 0
harvest <- lpjmlkit::read_io(
files_scenario$pft_harvestc,
subset = list(year = as.character(time_span_scenario))
) %>%
lpjmlkit::transform(to = c("year_month_day")) %>%
lpjmlkit::as_array(aggregate = list(month = sum)) %>%
suppressWarnings()
rharvest <- lpjmlkit::read_io(
files_scenario$pft_rharvestc,
subset = list(year = as.character(time_span_scenario))
) %>%
lpjmlkit::transform(to = c("year_month_day")) %>%
lpjmlkit::as_array(aggregate = list(month = sum)) %>%
suppressWarnings()
timber <- abind::adrop(lpjmlkit::read_io(
files_scenario$timber_harvestc,
subset = list(year = as.character(time_span_scenario))
) %>%
lpjmlkit::transform(to = c("year_month_day")) %>%
lpjmlkit::as_array(aggregate = list(month = sum)) %>%
# drop() %>%
suppressWarnings(), drop = "band")
if (include_fire) {
# read fire in monthly res. if possible, then multiply with monthly
# human/total ignition frac and aggregate to yearly. Otherwise aggregate
# human/total ignition frac to yearly and multiply with yearly firec
fire_raw <- lpjmlkit::read_io(
files_scenario$firec,
subset = list(year = as.character(time_span_scenario))
) %>%
lpjmlkit::transform(to = c("year_month_day")) %>%
lpjmlkit::as_array(aggregate = list(band = sum)) %>%
# drop() %>%
suppressWarnings()
if (external_fire) {
load(external_fire_file) # frac = c(cell,month,year)
}
if ("month" %in% names(dim(fire_raw))) {
if (external_fire) {
fire <- apply(
fire_raw *
lpjmlkit::asub(frac,
year = time_span_scenario,
drop = FALSE
),
c("cell", "year"),
sum,
na.rm = TRUE
) # gC/m2
rm(frac)
} else {
fire <- apply(
fire_raw,
c("cell", "year"),
sum,
na.rm = TRUE
) # gC/m2
}
rm(fire_raw)
} else {
if (external_fire) {
frac_yearly <- apply(
lpjmlkit::asub(frac,
year = time_span_scenario,
drop = FALSE
),
c("cell", "year"),
mean,
na.rm = TRUE
) # gC/m2
fire <- fire_raw * frac_yearly
rm(frac_yearly, frac)
}
}
gc()
} else {
fire <- timber * 0
}
if (external_wood_harvest) {
load(external_wood_harvest_file) # wh_lpj in kgC
# from kgC to gC/m2
wood_harvest <- (
lpjmlkit::asub(wh_lpj, year = time_span_scenario, drop = FALSE) *
10^3 / cellarea
)
# the division can lead to NAs
wood_harvest[is.na(wood_harvest)] <- 0
rm(wh_lpj)
gc()
} else {
wood_harvest <- fire * 0
}
cftfrac <- 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)) %>%
suppressWarnings()
npp_potential <- abind::adrop(lpjmlkit::read_io(
files_baseline$npp,
subset = list(year = as.character(time_span_baseline))
) %>%
lpjmlkit::transform(to = c("year_month_day")) %>%
lpjmlkit::as_array(aggregate = list(month = sum)) %>%
# drop() %>% suppressWarnings()
suppressWarnings(), drop = "band") # gC/m2
npp_potential[npp_potential < epsilon] <- 0
fpc <- lpjmlkit::read_io(
files_scenario$fpc,
subset = list(year = as.character(time_span_scenario))
) %>%
lpjmlkit::transform(to = c("year_month_day")) %>%
lpjmlkit::as_array(subset = list(band = "natural stand fraction")) %>%
suppressWarnings()
pftbands <- lpjmlkit::read_meta(files_scenario$fpc)$nbands - 1
} 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,
")"
)
}
bp_bands <- c(15, 16, 31, 32)
grass_bands <- c(14, 30)
nat_bands <- seq_len(pftbands)
if (!gridbased) { # needs to be scaled with standfrac
pftnpp[, , nat_bands] <- lpjmlkit::asub(pftnpp,
band = nat_bands,
drop = FALSE
) *
lpjmlkit::asub(fpc, band = rep(
"natural stand fraction",
pftbands
), drop = FALSE)
pftnpp[, , -c(nat_bands)] <- lpjmlkit::asub(pftnpp,
band = -c(nat_bands),
drop = FALSE
) * cftfrac
harvest <- harvest * cftfrac
rharvest <- rharvest * cftfrac
}
pftnpp_grasslands <- apply(
lpjmlkit::asub(pftnpp, band = (pftbands + grass_bands), drop = FALSE),
c("cell", "year"),
sum
) # gC/m2 only from grassland bands
pftnpp_cft <- apply(
lpjmlkit::asub(pftnpp,
band = -c(nat_bands, pftbands + grass_bands, pftbands + bp_bands),
drop = FALSE
),
c("cell", "year"),
sum
) # gC/m2 not from grassland and bioenergy bands
pftnpp_bioenergy <- apply(
lpjmlkit::asub(pftnpp, band = pftbands + bp_bands, drop = FALSE),
c("cell", "year"),
sum
) # gC/m2 only from bioenergy bands
pftnpp_nat <- apply(
lpjmlkit::asub(pftnpp, band = nat_bands, drop = FALSE),
c("cell", "year"),
sum
) # gC/m2
if (is.null(files_reference)) {
pi_window <- 3:32
npp_ref <- lpjmlkit::asub(npp_potential, year = pi_window, drop = FALSE)
} # npp_ref
# lpjmlkit::asub(INARRAY, band = BANDS, drop = FALSE)
harvest_grasslands <- apply(
lpjmlkit::asub(harvest, band = grass_bands, drop = FALSE),
c("cell", "year"),
sum
) # gC/m2 only from grassland bands
harvest_bioenergy <- apply(
lpjmlkit::asub(harvest, band = bp_bands, drop = FALSE),
c("cell", "year"),
sum
) # gC/m2 only from bioenergy bands
harvest_cft <- apply(
lpjmlkit::asub(harvest, band = -c(grass_bands, bp_bands), drop = FALSE),
c("cell", "year"),
sum
) # gC/m2 not from grassland and bioenergy bands
rharvest_cft <- apply(
lpjmlkit::asub(rharvest, band = -c(grass_bands, bp_bands), drop = FALSE),
c("cell", "year"),
sum
) # gC/m2 not from grassland and bioenergy bands
if (save_data) {
if (!file.exists(data_file)) {
message("Writing data file: ", data_file)
} else {
message(
"Data file (",
data_file,
") already exists, old file renamed to: ",
data_file,
"_sav"
)
file.rename(data_file, paste0(data_file, "_sav"))
}
save(npp_potential,
npp,
npp_ref,
pftnpp_cft,
pftnpp_nat,
pftnpp_grasslands,
pftnpp_bioenergy,
harvest_cft,
rharvest_cft,
fire,
timber,
fpc,
cftfrac,
harvest_grasslands,
harvest_bioenergy,
wood_harvest,
lat,
lon,
cellarea,
file = data_file
)
}
}
message("Calculating data")
if (grass_scaling) {
load(grass_harvest_file)
nregs <- length(grazing_data$name)
lpj_grass_harvest_region <- array(0, dim = nregs)
lpj_grass_harvest_2000 <- rowMeans(
harvest_grasslands[, (1995 - start_year + 1):(2005 - start_year + 1)]
) * cellarea / 1000 * 2 # from gC/m2 to kgDM
grassland_scaling_factor_cellwise <- array(1, dim = grid$ncells)
for (r in seq_len(nregs)) {
lpj_grass_harvest_region[r] <- sum(
lpj_grass_harvest_2000[which(mapping_lpj67420_to_grazing_regions == r)]
)
}
scaling_factor <- (
grazing_data$Herrero_2000_kgDM_by_region / lpj_grass_harvest_region
)
for (r in seq_len(nregs)) {
grassland_scaling_factor_cellwise[
which(mapping_lpj67420_to_grazing_regions == r)
] <- scaling_factor[r]
}
harvest_grasslands <- harvest_grasslands * rep(
grassland_scaling_factor_cellwise,
times = length(harvest_grasslands[1, ])
)
}
npp_act_overtime <- colSums(npp * cellarea) / 10^15 # gC/m2 to GtC
npp_pot_overtime <- colSums(npp_potential * cellarea) / 10^15 # gC/m2 to GtC
npp_eco_overtime <- colSums(pftnpp_nat * cellarea) / 10^15 # gC/m2 to GtC
npp_luc_overtime <- npp_pot_overtime - npp_act_overtime
harvest_cft_overtime <- colSums(
harvest_cft * cellarea
) / 10^15 # gC/m2 to GtC
rharvest_cft_overtime <- colSums(
rharvest_cft * cellarea
) / 10^15 # gC/m2 to GtC
harvest_grasslands_overtime <- colSums(
harvest_grasslands * cellarea
) / 10^15 # gC/m2 to GtC
harvest_bioenergy_overtime <- colSums(
harvest_bioenergy * cellarea
) / 10^15 # gC/m2 to GtC
timber_harvest_overtime <- colSums(
timber * cellarea
) / 10^15 # gC/m2 to GtC
fire_overtime <- colSums(
fire * cellarea
) / 10^15 # gC/m2 to GtC
wood_harvest_overtime <- colSums(
wood_harvest * cellarea
) / 10^15 # gC/m2 to GtC
if (include_fire) {
npp_harv_overtime <- harvest_cft_overtime + rharvest_cft_overtime +
harvest_grasslands_overtime + harvest_bioenergy_overtime +
timber_harvest_overtime + fire_overtime + wood_harvest_overtime
} else {
npp_harv_overtime <- harvest_cft_overtime + rharvest_cft_overtime +
harvest_grasslands_overtime + harvest_bioenergy_overtime +
timber_harvest_overtime + wood_harvest_overtime
}
biocol_overtime <- npp_harv_overtime + npp_luc_overtime
biocol_overtime_frac_piref <- (
biocol_overtime / mean(colSums(npp_ref * cellarea) / 10^15)
)
biocol_overtime_frac <- (
biocol_overtime / npp_pot_overtime
)
biocol_luc <- npp_potential - npp
# pick a PI window that excludes onset effects, but is reasonable early
if (include_fire) {
biocol_harvest <- (
harvest_cft + rharvest_cft + harvest_grasslands + harvest_bioenergy +
timber + fire + wood_harvest
)
} else {
biocol_harvest <- (
harvest_cft + rharvest_cft + harvest_grasslands + harvest_bioenergy +
timber + wood_harvest
)
}
biocol <- biocol_harvest + biocol_luc
# set to 0 below lower threshold of NPP
biocol[abs(npp_potential) < npp_threshold] <- 0
# actual NPPpot as ref
biocol_frac <- biocol / npp_potential
# NPPpi as ref
biocol_frac_piref <- biocol / rowMeans(npp_ref)
# take the abs of biocol and sum that up for overtime
biocol_overtime_abs <- colSums(abs(biocol * cellarea)) / 10^15
biocol_overtime_abs_frac_piref <- biocol_overtime_abs * 10^15 /
mean(colSums(npp_ref * cellarea))
biocol_overtime_abs_frac <- biocol_overtime_abs / npp_pot_overtime
# take the abs of biocol and sum that up for overtime
biocol_pos <- biocol
biocol_pos[biocol_pos < 0] <- 0
biocol_overtime_pos <- colSums(biocol_pos * cellarea) / 10^15
biocol_overtime_pos_frac_piref <- biocol_overtime_pos * 10^15 /
mean(colSums(npp_ref * cellarea))
biocol_overtime_pos_frac <- biocol_overtime_pos / npp_pot_overtime
return(list(
biocol_overtime = biocol_overtime,
biocol_overtime_abs = biocol_overtime_abs,
biocol_overtime_abs_frac_piref = biocol_overtime_abs_frac_piref,
biocol_overtime_abs_frac = biocol_overtime_abs_frac,
biocol_overtime_pos = biocol_overtime_pos,
biocol_overtime_pos_frac_piref = biocol_overtime_pos_frac_piref,
biocol_overtime_pos_frac = biocol_overtime_pos_frac,
biocol_overtime_frac_piref = biocol_overtime_frac_piref,
biocol_overtime_frac = biocol_overtime_frac,
npp_harv_overtime = npp_harv_overtime,
npp_luc_overtime = npp_luc_overtime,
npp_act_overtime = npp_act_overtime,
npp_pot_overtime = npp_pot_overtime,
npp_eco_overtime = npp_eco_overtime,
harvest_grasslands_overtime = harvest_grasslands_overtime,
harvest_bioenergy_overtime = harvest_bioenergy_overtime,
harvest_cft_overtime = harvest_cft_overtime,
rharvest_cft_overtime = rharvest_cft_overtime,
fire_overtime = fire_overtime,
timber_harvest_overtime = timber_harvest_overtime,
wood_harvest_overtime = wood_harvest_overtime,
biocol = biocol,
biocol_frac = biocol_frac,
npp = npp,
biocol_frac_piref = biocol_frac_piref,
npp_potential = npp_potential,
npp_ref = npp_ref,
harvest_cft = harvest_cft,
rharvest_cft = rharvest_cft,
biocol_harvest = biocol_harvest,
biocol_luc = biocol_luc,
lat = lat,
lon = lon,
cellarea = cellarea
)) # , biocol_luc_piref = biocol_luc_piref))
}
#' Calculate BioCol
#'
#' Wrapper function to calculate BioCol
#'
#' @param path_lu folder of landuse scenario run
#' @param path_pnv folder of pnv reference run
#' @param start_year first year of simulations
#' @param stop_year last year of simulations
#' @param reference_npp_time_span time span to read reference npp from, using
#' index years 10:39 from potential npp input if set to NULL (default: NULL)
#' @param reference_npp_file file to read reference npp from, using
#' potential npp input if set to NULL (default: NULL)
#' @param gridbased logical are pft outputs gridbased or pft-based?
#' @param read_saved_data flag whether to read previously saved data
#' instead of reading it in from output files (default FALSE)
#' @param save_data whether to save input data to file (default FALSE)
#' @param data_file file to save/read input data to/from (default NULL)
#' @param include_fire boolean include firec in calculation of BioCol?
#' (default TRUE)
#' @param external_fire instead of reading in firec for fire emissions, read in
#' this external firec file from a separate spitfire run with disabled
#' lighning. this will then include only human induced fires
#' (default FALSE)
#' @param external_wood_harvest include external wood harvest from LUH2_v2h
#' (default FALSE)
#' @param grass_scaling whether to scale pasture harvest according to
#' data given via grass_harvest_file (default FALSE)
#' @param npp_threshold lower threshold for npp (to mask out non-lu areas
#' according to Haberl et al. 2007). Below BioCol will be set to 0.
#' (default: 20 gC/m2)
#' @param grass_harvest_file file containing grazing data to rescale the
#' grassland harvests according to Herrero et al. 2013. File contains:
#' grazing_data list object with $name and $id of 29 world regions, and
#' $Herrero_2000_kgDM_by_region containing for each of these regions and
#' mapping_lpj67420_to_grazing_regions array with a mapping between 67420
#' LPJmL cells and the 29 regions
#' @param external_fire_file path to external file with human induced fire
#' fraction c(cell,month,year) since 1500
#' @param external_wood_harvest_file path to R-file containing processed
#' timeline of maps for LUH2_v2h woodharvest
#' @param replace_input_file_names list with alternative names for output
#' identifiers to replace the ones in inst/ext_files/metric_files.yml.
#' e.g. list(npp="mnpp") would replace the expected output for npp with
#' mnpp followed by the automatically detected file extension (.bin.json)
#' @param suppress_warnings suppress warnings when reading files (default: TRUE)
#'
#' @return list data object containing BioCol and components as arrays: biocol,
#' biocol_overtime, biocol_overtime_piref, biocol_frac, npp_potential,
#' biocol_overtime_abs_frac_piref, biocol_frac_piref, npp_act_overtime,
#' npp_pot_overtime, npp_eco_overtime, npp_ref, harvest_cft_overtime,
#' npp_luc_overtime, rharvest_cft_overtime, fire_overtime,
#' timber_harvest_overtime, harvest_cft, rharvest_cft,
#' wood_harvest_overtime, biocol_harvest, biocol_luc, lat, lon, cellarea
#'
#' @examples
#' \dontrun{
#' calc_biocol(
#' path_lu = run_folder,
#' path_pnv = pnv_folder,
#' gridbased = TRUE,
#' start_year = 1980,
#' stop_year = 2014,
#' reference_npp_time_span = 1510:1539,
#' read_saved_data = FALSE,
#' save_data = FALSE,
#' npp_threshold = 20,
#' )
#' }
#'
#' @md
#' @export
calc_biocol <- function(
path_lu,
path_pnv,
start_year,
stop_year,
reference_npp_time_span = NULL,
reference_npp_file = NULL,
gridbased = TRUE,
read_saved_data = FALSE,
save_data = FALSE,
data_file = NULL,
include_fire = FALSE,
external_fire = FALSE,
external_wood_harvest = FALSE,
grass_scaling = FALSE,
npp_threshold = 20,
grass_harvest_file = NULL,
external_fire_file = NULL,
external_wood_harvest_file = NULL,
replace_input_file_names = NULL,
suppress_warnings = TRUE) {
metric_files <- system.file(
"extdata",
"metric_files.yml",
package = "biospheremetrics"
) %>%
yaml::read_yaml()
# translate output names (from metric_files.yml)
# and folders to files_scenarios/reference lists
file_extension <- get_major_file_ext(paste0(path_lu))
files_scenario <- list()
files_baseline <- list()
for (output in names(metric_files$metric$biocol$output)) {
# Iterate over all outputs
if (is.null(replace_input_file_names[[output]])) {
for (file in metric_files$file_name[[output]]) {
full_file_path_lu <- paste0(path_lu, file, ".", file_extension)
if (file.exists(full_file_path_lu)) {
files_scenario[[output]] <- full_file_path_lu
}
full_file_path_pnv <- paste0(path_pnv, file, ".", file_extension)
if (file.exists(full_file_path_pnv)) {
files_baseline[[output]] <- full_file_path_pnv
}
}
if (is.null(files_scenario[[output]])) {
stop(
"None of the default file names for ", output,
" were found in ", path_lu, "please check or define manually",
" using argument 'replace_input_file_names'. Stopping."
)
}
if (is.null(files_baseline[[output]])) {
stop(
"None of the default file names for ", output,
" were found in ", path_pnv, "please check or define manually",
" using argument 'replace_input_file_names'. Stopping."
)
}
} else {
files_scenario[[output]] <- paste0(
path_lu,
replace_input_file_names[[output]], ".", file_extension
)
files_baseline[[output]] <- paste0(
path_pnv,
replace_input_file_names[[output]], ".", file_extension
)
}
}
if (is.null(reference_npp_file)) reference_npp_file <- files_baseline$npp
files_reference <- list(
npp = reference_npp_file
)
return(
read_calc_biocol(
files_scenario = files_scenario,
files_baseline = files_baseline,
files_reference = files_reference,
time_span_scenario = as.character(start_year:stop_year),
time_span_baseline = as.character(start_year:stop_year),
time_span_reference = reference_npp_time_span,
gridbased = gridbased,
read_saved_data = read_saved_data,
save_data = save_data,
data_file = data_file,
include_fire = include_fire,
external_fire = external_fire,
external_wood_harvest = external_wood_harvest,
grass_scaling = grass_scaling,
npp_threshold = npp_threshold,
grass_harvest_file = grass_harvest_file,
external_fire_file = external_fire_file,
external_wood_harvest_file = external_wood_harvest_file,
suppress_warnings = suppress_warnings
)
)
}