Newer
Older
# 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) {
# TODO: compare length time_span_reference and time_span_scenario
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"]),
Fabian Stenzel
committed
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"]),
Fabian Stenzel
committed
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"
)
# first read in all lpjml output files required for computing EcoRisks
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
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))
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
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
}
############## 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(
)
# 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
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
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 ################
Fabian Stenzel
committed
# 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)
ref = state_ref,
scen = state_scen,
Fabian Stenzel
committed
ref = state_ref,
scen = state_scen,
local = FALSE,
cell_area = cell_area
) # global importance
ref = state_ref,
scen = state_scen
) # ecosystem balance
ref = state_ref[, , non_nitrogen_dimensions],
scen = state_scen[, , non_nitrogen_dimensions],
ref = state_ref[, , non_nitrogen_dimensions],
scen = state_scen[, , non_nitrogen_dimensions],
ref = state_ref[, , non_nitrogen_dimensions],
scen = state_scen[, , non_nitrogen_dimensions]
) # ecosystem balance
}
if (dimensions_only_local == TRUE) {
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
# 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")],
# 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,
# 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")],
}
} else { # local == FALSE
cf <- (
calc_component(
ref = state_ref[, , c("carbon_influx", "carbon_outflux")],
scen = state_scen[, , c("carbon_influx", "carbon_outflux")],
ref = state_ref[, , c("carbon_influx", "carbon_outflux")],
scen = state_scen[, , c("carbon_influx", "carbon_outflux")],
local = FALSE,
cell_area = cell_area
ref = state_ref[, , c("carbon_influx", "carbon_outflux")],
scen = state_scen[, , c("carbon_influx", "carbon_outflux")]
# 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 +
ref = state_ref[, , c("vegetation_carbon_pool", "soil_carbon_pool")],
scen = state_scen[, , c("vegetation_carbon_pool", "soil_carbon_pool")],
ref = state_ref[, , c("vegetation_carbon_pool", "soil_carbon_pool")],
scen = state_scen[, , c("vegetation_carbon_pool", "soil_carbon_pool")]
# 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
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")]
# 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 +
ref = state_ref[, , c("water_influx", "water_outflux")],
scen = state_scen[, , c("water_influx", "water_outflux")],
local = FALSE,
cell_area = cell_area
ref = state_ref[, , c("water_influx", "water_outflux")],
scen = state_scen[, , c("water_influx", "water_outflux")]
# 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 +
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
ref = state_ref[, , c("water_influx", "water_outflux", "soil_water_pool")],
scen = state_scen[, , c("water_influx", "water_outflux", "soil_water_pool")]
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 +
ref = state_ref[, , c("vegetation_nitrogen_pool", "soil_mineral_nitrogen_pool")],
scen = state_scen[, , c("vegetation_nitrogen_pool", "soil_mineral_nitrogen_pool")],
ref = state_ref[, , c("vegetation_nitrogen_pool", "soil_mineral_nitrogen_pool")],
scen = state_scen[, , c("vegetation_nitrogen_pool", "soil_mineral_nitrogen_pool")]
# 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 +
ref = state_ref[, , c("nitrogen_influx", "nitrogen_outflux")],
scen = state_scen[, , c("nitrogen_influx", "nitrogen_outflux")],
local = FALSE,
cell_area = cell_area
ref = state_ref[, , c("nitrogen_influx", "nitrogen_outflux")],
scen = state_scen[, , c("nitrogen_influx", "nitrogen_outflux")]
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
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")]
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", ]
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,
carbon_total = ct,
water_stocks = NA,
water_total = wt,
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,
carbon_total = ct,
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(
Fabian Stenzel
committed
files_reference$grid
Fabian Stenzel
committed
cell_area <- drop(lpjmlkit::read_io(
filename = files_reference$terr_area
Fabian Stenzel
committed
ncells <- length(lat)
nyears <- length(time_span_scenario)
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
### 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))
Fabian Stenzel
committed
#### new input reading ###
metric_files <- system.file(
"extdata",
"metric_files.yml",
package = "biospheremetrics"
) %>%
yaml::read_yaml()
Fabian Stenzel
committed
nclasses <- length(metric_files$metric$ecorisk_nitrogen$metric_class)
nstate_dimensions <- 0
for (i in seq_len(nclasses)) {
nstate_dimensions <- nstate_dimensions +
length(metric_files$metric$ecorisk_nitrogen$metric_class[[i]])
}
state_ref <- array(0, dim = c(ncells, nyears, nstate_dimensions))
state_scen <- array(0, dim = c(ncells, nyears, nstate_dimensions))
class_names <- seq_len(nstate_dimensions)
Fabian Stenzel
committed
index <- 1
# iterate over main classes (carbon pools, water fluxes ...)
Fabian Stenzel
committed
classe <- metric_files$metric$ecorisk_nitrogen$metric_class[[c]]
nsubclasses <- length(classe)
# iterate over subclasses (vegetation carbon, soil water ...)
Fabian Stenzel
committed
subclass <- classe[s]
class_names[index] <- names(subclass)
vars <- split_sign(unlist(subclass))
Fabian Stenzel
committed
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]
))
Fabian Stenzel
committed
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()
Fabian Stenzel
committed
} else {
stop(paste("Couldn't read in:", path_scen_file, " - stopping!"))
Fabian Stenzel
committed
}
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]
))
Fabian Stenzel
committed
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()
Fabian Stenzel
committed
} else {
stop(paste("Couldn't read in:", path_ref_file, " - stopping!"))
Fabian Stenzel
committed
}
# 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
Fabian Stenzel
committed
}
Fabian Stenzel
committed
}
index <- index + 1
Fabian Stenzel
committed
}
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)
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
} 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"
#' )