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 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",
time_span_reference,
time_span_scenario,
dimensions_only_local = FALSE,
overtime = FALSE,
window = 30,
debug = FALSE,
external_variability = FALSE,
c2vr = NULL) {
Fabian Stenzel
committed
# check timespan consistency
nyears <- length(time_span_reference)
nyears_scen <- length(time_span_scenario)
Fabian Stenzel
committed
if ( (! nyears == window) || nyears_scen < window) {
stop(paste0("Timespan in reference is not equal to window size (", window,
"), or scenario timespan is smaller than window size."))
Fabian Stenzel
committed
# translate output names (from metric_files.yml) and folders to files_scenarios/reference lists
metric_files <- system.file(
"extdata",
"metric_files.yml",
package = "biospheremetrics"
) %>%
yaml::read_yaml()
file_extension <- get_major_file_ext(paste0(path_scen))
outputs <- metric_files$metric$ecorisk_nitrogen$output
Fabian Stenzel
committed
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
grid = paste0(path_scen, outputs$grid$name, ".", file_extension),
terr_area = paste0(path_scen, outputs$terr_area$name, ".", file_extension),
fpc = paste0(path_scen, outputs$fpc$name, ".", file_extension),
fpc_bft = paste0(path_scen, outputs$fpc_bft$name, ".", file_extension),
cftfrac = paste0(path_scen, outputs$cftfrac$name, ".", file_extension),
firec = paste0(path_scen, outputs$firec$name, ".", file_extension),
npp = paste0(path_scen, outputs$npp$name, ".", file_extension),
runoff = paste0(path_scen, outputs$runoff$name, ".", file_extension),
transp = paste0(path_scen, outputs$transp$name, ".", file_extension),
vegc = paste0(path_scen, outputs$vegc$name, ".", file_extension),
firef = paste0(path_scen, outputs$firef$name, ".", file_extension),
harvestc = paste0(path_scen, outputs$harvestc$name, ".", file_extension),
evap = paste0(path_scen, outputs$evap$name, ".", file_extension),
interc = paste0(path_scen, outputs$interc$name, ".", file_extension),
soilc = paste0(path_scen, outputs$soilc$name, ".", file_extension),
litc = paste0(path_scen, outputs$litc$name, ".", file_extension),
swc = paste0(path_scen, outputs$swc$name, ".", file_extension),
swc_vol = paste0(path_scen, outputs$swc_vol$name, ".", file_extension),
swe = paste0(path_scen, outputs$swe$name, ".", file_extension),
vegn = paste0(path_scen, outputs$vegn$name, ".", file_extension),
soilnh4 = paste0(path_scen, outputs$soilnh4$name, ".", file_extension),
soilno3 = paste0(path_scen, outputs$soilno3$name, ".", file_extension),
leaching = paste0(path_scen, outputs$leaching$name, ".", file_extension),
n2o_denit = paste0(path_scen, outputs$n2o_denit$name, ".", file_extension),
n2o_nit = paste0(path_scen, outputs$n2o_nit$name, ".", file_extension),
n2_emis = paste0(path_scen, outputs$n2_emis$name, ".", file_extension),
bnf = paste0(path_scen, outputs$bnf$name, ".", file_extension),
n_volatilization = paste0(path_scen, outputs$n_volatilization$name, ".", file_extension),
gpp = paste0(path_scen, outputs$gpp$name, ".", file_extension),
res_storage = paste0(path_scen, outputs$res_storage$name, ".", file_extension),
lakevol = paste0(path_scen, outputs$lakevol$name, ".", file_extension),
prec = paste0(path_scen, outputs$prec$name, ".", file_extension),
irrig = paste0(path_scen, outputs$irrig$name, ".", file_extension),
nfert_agr = paste0(path_scen, outputs$nfert_agr$name, ".", file_extension),
nmanure_agr = paste0(path_scen, outputs$nmanure_agr$name, ".", file_extension),
ndepos = paste0(path_scen, outputs$ndepos$name, ".", file_extension),
firen = paste0(path_scen, outputs$firen$name, ".", file_extension),
harvestn = paste0(path_scen, outputs$harvestn$name, ".", file_extension),
irrig_stor = paste0(path_scen, outputs$irrig_stor$name, ".", file_extension),
rivervol = paste0(path_scen, outputs$rivervol$name, ".", file_extension)
Fabian Stenzel
committed
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
grid = paste0(path_ref, outputs$grid$name, ".", file_extension),
terr_area = paste0(path_ref, outputs$terr_area$name, ".", file_extension),
fpc = paste0(path_ref, outputs$fpc$name, ".", file_extension),
fpc_bft = paste0(path_ref, outputs$fpc_bft$name, ".", file_extension),
cftfrac = paste0(path_ref, outputs$cftfrac$name, ".", file_extension),
firec = paste0(path_ref, outputs$firec$name, ".", file_extension),
npp = paste0(path_ref, outputs$npp$name, ".", file_extension),
runoff = paste0(path_ref, outputs$runoff$name, ".", file_extension),
transp = paste0(path_ref, outputs$transp$name, ".", file_extension),
vegc = paste0(path_ref, outputs$vegc$name, ".", file_extension),
firef = paste0(path_ref, outputs$firef$name, ".", file_extension),
harvestc = paste0(path_ref, outputs$harvestc$name, ".", file_extension),
evap = paste0(path_ref, outputs$evap$name, ".", file_extension),
interc = paste0(path_ref, outputs$interc$name, ".", file_extension),
soilc = paste0(path_ref, outputs$soilc$name, ".", file_extension),
litc = paste0(path_ref, outputs$litc$name, ".", file_extension),
swc = paste0(path_ref, outputs$swc$name, ".", file_extension),
swc_vol = paste0(path_ref, outputs$swc_vol$name, ".", file_extension),
swe = paste0(path_ref, outputs$swe$name, ".", file_extension),
vegn = paste0(path_ref, outputs$vegn$name, ".", file_extension),
soilnh4 = paste0(path_ref, outputs$soilnh4$name, ".", file_extension),
soilno3 = paste0(path_ref, outputs$soilno3$name, ".", file_extension),
leaching = paste0(path_ref, outputs$leaching$name, ".", file_extension),
n2o_denit = paste0(path_ref, outputs$n2o_denit$name, ".", file_extension),
n2o_nit = paste0(path_ref, outputs$n2o_nit$name, ".", file_extension),
n2_emis = paste0(path_ref, outputs$n2_emis$name, ".", file_extension),
bnf = paste0(path_ref, outputs$bnf$name, ".", file_extension),
n_volatilization = paste0(path_ref, outputs$n_volatilization$name, ".", file_extension),
gpp = paste0(path_ref, outputs$gpp$name, ".", file_extension),
res_storage = paste0(path_ref, outputs$res_storage$name, ".", file_extension),
lakevol = paste0(path_ref, outputs$lakevol$name, ".", file_extension),
prec = paste0(path_ref, outputs$prec$name, ".", file_extension),
irrig = paste0(path_ref, outputs$irrig$name, ".", file_extension),
nfert_agr = paste0(path_ref, outputs$nfert_agr$name, ".", file_extension),
nmanure_agr = paste0(path_ref, outputs$nmanure_agr$name, ".", file_extension),
ndepos = paste0(path_ref, outputs$ndepos$name, ".", file_extension),
firen = paste0(path_ref, outputs$firen$name, ".", file_extension),
harvestn = paste0(path_ref, outputs$harvestn$name, ".", file_extension),
irrig_stor = paste0(path_ref, outputs$irrig_stor$name, ".", file_extension),
rivervol = paste0(path_ref, outputs$rivervol$name, ".", file_extension)
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)) {
message("Loading saved data from:", save_data)
stop(
"save_data is not specified as parameter, ",
"nothing to load ... exiting"
)
# first read in all lpjml output files required for computing EcoRisks
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
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)),
Fabian Stenzel
committed
nitrogen_total = array(0, dim = c(ncells, slices)),
lat = lat,
lon = lon
message("Calculating time slice ", y, " of ", slices)
returned <- calc_ecorisk(
fpc_ref = fpc_ref,
fpc_scen = fpc_scen[, , y:(y + window - 1)],
bft_ref = bft_ref,
bft_scen = bft_scen[, , y:(y + window - 1)],
cft_ref = cft_ref,
cft_scen = cft_scen[, , y:(y + window - 1)],
state_ref = state_ref,
state_scen = state_scen[, y:(y + window - 1), ],
weighting = weighting,
lat = lat,
lon = lon,
cell_area = cell_area,
dimensions_only_local = dimensions_only_local,
nitrogen = nitrogen,
external_variability = external_variability,
c2vr = c2vr
)
ecorisk$ecorisk_total[, y] <- returned$ecorisk_total
ecorisk$vegetation_structure_change[, y] <- (
returned$vegetation_structure_change
)
ecorisk$local_change[, y] <- returned$local_change
ecorisk$global_importance[, y] <- returned$global_importance
ecorisk$ecosystem_balance[, y] <- returned$ecosystem_balance
ecorisk$c2vr[, , y] <- returned$c2vr
ecorisk$carbon_stocks[, y] <- returned$carbon_stocks
ecorisk$carbon_fluxes[, y] <- returned$carbon_fluxes
ecorisk$carbon_total[, y] <- returned$carbon_total
ecorisk$water_total[, y] <- returned$water_total
ecorisk$water_fluxes[, y] <- returned$water_fluxes
if (nitrogen) {
ecorisk$nitrogen_stocks[, y] <- returned$nitrogen_stocks
ecorisk$nitrogen_fluxes[, y] <- returned$nitrogen_fluxes
ecorisk$nitrogen_total[, y] <- returned$nitrogen_total
}
############## export and save data if requested #############
if (!(is.null(save_ecorisk))) {
message("Saving EcoRisk data to: ", save_ecorisk)
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
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
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
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) {
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
# 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)
### 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()
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)
Fabian Stenzel
committed
"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() %>%
suppressWarnings()
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)
Fabian Stenzel
committed
"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() %>%
suppressWarnings()
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)
} 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)
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
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
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))