Climate Inputs
var_info <- data.table(variable = c('dlwrf', 'nbdsf', 'nddsf', 'vbdsf',
'vddsf', 'prate', 'pres', 'hgt',
'ugrd', 'vgrd', 'sh', 'tmp'),
description = c('downward long wave radiation',
'near infrared beam downward solar radiation',
'near IR diffuse downward solar radiation',
'visible beam downward solar radiation',
'visible diffuse downward solar radiation',
'precipitation rate',
'atmospheric pressure',
'geopotential height',
'zonal wind',
'meridional wind',
'specific humidity',
'air temperature'),
units = c('W m-2', 'W m-2', 'W m-2', 'W m-2',
'W m-2', 'kgH2O m-2 s-1', 'Pa', 'm',
'm s-1', 'm s-', 'kgH2O kgAir-1', 'K'), stringsAsFactors=FALSE)
list.files(file.path(BASE_DIR, 'ED-outputs', 'met_results'), "constant_yr", full.names = TRUE) %>%
lapply(function(x){
d <- read.csv(x, stringsAsFactors = FALSE)
met <- gsub(pattern = '-constant_yr.csv|NARR-ED2_met-',replacement = "", basename(x))
d$met <- met
return(d)
}) %>%
do.call(what = "rbind") %>%
as.data.table() ->
annual_data
annnual_data <- melt.data.table(annual_data, id.vars = c("year", "met"), variable.name = "variable",
value.name = "value")
annnual_data <- annnual_data[var_info, on = 'variable']
annnual_data$met <- as.character(annnual_data$met)
annnual_data %>%
ggplot(aes(year, value, color = met)) +
geom_line() +
facet_wrap("description", scales = 'free', strip.position = "top",
labeller = label_wrap_gen(width = 20, multi_line = TRUE)) +
labs(title = 'Annual Average Met Data', y = NULL)
list.files(file.path(BASE_DIR, 'ED-outputs', 'met_results'), "constant_month", full.names = TRUE) %>%
lapply(function(x){
d <- read.csv(x, stringsAsFactors = FALSE)
met <- gsub(pattern = '-constant_month.csv|NARR-ED2_met-',replacement = "", basename(x))
d$met <- met
return(d)
}) %>%
do.call(what = "rbind") %>%
as.data.table() ->
monthly_data
monthly_data <- melt.data.table(monthly_data, id.vars = c("time", "met", "month"), variable.name = "variable",
value.name = "value")
monthly_data <- monthly_data[var_info, on = 'variable']
monthly_data <- na.omit(monthly_data)
monthly_data$month <- factor(x = monthly_data$month, levels = toupper(month.abb), ordered = TRUE)
monthly_data$met <- as.character(monthly_data$met)
mon <- data.table(month = toupper(month.abb), mon_time = 1:12)
monthly_data <- monthly_data[mon, on = 'month']
monthly_data %>%
.[,.(value = mean(value)), by = c("met", "month", "variable", "description", "mon_time")] %>%
ggplot(aes(mon_time, value, color = met)) +
geom_line() +
facet_wrap("description", scales = 'free', strip.position = "top",
labeller = label_wrap_gen(width = 20, multi_line = TRUE)) +
labs(title = 'Monthly Met', y = NULL, x = "Month")
Baseline Scenario
ED_outputs[scn == "harvest_0"] %>%
ggplot() +
geom_line(aes(year, value, color = met)) +
facet_wrap("variable", scales = "free") +
THEME +
labs(title = "Control Run Idealized Climate Inputs",
y = "MgC ha-1 year-1",
x = "Year")
Stability Metrics
ED_outputs %>%
ln_disturbance() %>%
add_severity_lables() ->
disturbance_df
disturbance_df %>%
ggplot() +
geom_line(aes(year, value, color = severity, line = met), alpha = 0.5, size = 1.25) +
geom_hline(yintercept = 0, linetype = 2, size = 1) +
facet_wrap("variable", scales = "free") +
labs(x = 'Years since disturbance', y = "ln(disturbance/control)") +
scale_color_manual(values = FORTE_SEVERITY_COLORS) +
theme_bw(base_size = 20) +
theme(legend.position = "none")
# Clculate the resistance of the trough and when that trough occurs
# Args:
# data: ED data!
# Return: a data frame of the trough resistance value
through_resistance <- function(data){
d <- ln_disturbance(data)
dd <- d[, .(year, scn, variable, met, value)]
info <- unique(data[, .(variable, description)])
rslt <- do.call(what = "rbind",
args = split(dd, interaction(dd$scn, dd$variable, dd$met), drop = TRUE) %>%
lapply(function(x){
index <- which.min(x$value)
out <- x[index, ]
out <- out[ , .(scn, trough_year = year, variable, met, trough_resistance = value)]
return(out)
}))
out <- rslt[info, on = "variable"]
out$unit <- "unitless"
return(out)
}
ED_outputs %>%
through_resistance() %>%
add_severity_lables() %>%
na.omit() %>%
rename(value = trough_resistance) %>%
mutate(type = "trough resistance") ->
trough_resistance
ED_outputs %>%
ln_disturbance() %>%
na.omit() %>%
add_severity_lables() %>%
filter(year == 0) %>%
mutate(type = "yr resistance") ->
yr_resistance
ggplot(data = trough_resistance) +
geom_boxplot(aes(variable, value, fill = severity)) +
scale_fill_manual(values = FORTE_SEVERITY_COLORS) +
THEME +
labs(y = "Resistance", x = "Carbon Flux", title = "Resistance")
disturbance_df %>%
left_join(select(trough_resistance, met, variable, severity, trough_year)) %>%
dplyr::filter(year >= trough_year) %>%
dplyr::mutate(year = year - trough_year) %>%
select(-trough_year) ->
resillience_data1
## Joining, by = c("variable", "met", "severity")
resillience_data1 %>%
dplyr::mutate(label = paste0(description, "(", unit, ")")) %>%
ggplot(aes(year, value, color = severity, line = met)) +
geom_line() +
facet_wrap('variable', scales = 'free',
labeller = label_wrap_gen(width = 30, multi_line = TRUE), ncol = 2) +
THEME +
labs(y = 'ln(Disturbance/Baseline)',
x = 'Year post trough resistance',
title = "Time series post the trough resistance") +
scale_color_manual(values = FORTE_SEVERITY_COLORS) +
guides(linetype = FALSE)
df_list <- split(resillience_data1,
interaction(resillience_data1$scn, resillience_data1$variable, resillience_data1$met, sep = "~"),
drop = TRUE)
df_list1 <- df_list
fit_data <- function(d){
x <- d$year
y <- d$value
lin.mod <- lm(y ~ x)
segmented.mod <- lin.mod
tryCatch({
segmented.mod <- segmented(lin.mod, seg.Z = ~x)
},
error = function(e){lin.mod})
out <- list("lin.mod" = lin.mod, "seg.mod" = segmented.mod)
return(out)
}
extract_fit_data <- function(l, n){
ln_slope <- l[["lin.mod"]]$coefficients[2]
# If it was determined that there was no segment
if(length(l[["seg.mod"]]) == 12){
seg_slope1 <- l[["seg.mod"]]$coefficients[2]
seg_slope2 <- NA
seg_brk <- NA
} else {
slopes <- slope(l[["seg.mod"]])
seg_slope1 <- slopes$x[1,1] # Extract the slope of the first segment
seg_slope2 <- slopes$x[2,1] # Extract the slope of the second segment
seg_brk <- l[["seg.mod"]]$psi[1,3] # Extract the year of the break point
}
# Extract the data about the fit.
info <- unlist(strsplit(n, split = "~"))
# Format the slope into a data table
dt <- data.table(scn = info[1],
variable = info[2],
met = info[3],
lin_slope = ln_slope,
seg_slope1 = seg_slope1,
seg_slope2 = seg_slope2,
seg_brk = seg_brk)
}
fits <- lapply(df_list, fit_data)
fits1 <- fits
mapply(FUN = extract_fit_data, l = fits, n = names(fits), SIMPLIFY = FALSE) %>%
do.call(what = "rbind") %>%
add_severity_lables() ->
compare_slope_fits1
compare_slope_fits1 %>%
add_severity_lables() %>%
select(scn, variable, met, lin_slope, seg_slope1, severity) %>%
ggplot() +
geom_boxplot(aes(variable, seg_slope1, fill = severity)) +
scale_fill_manual(values = FORTE_SEVERITY_COLORS) +
THEME +
labs(y = "Resistance", x = "Carbon Flux", title = "Resilience")
compare_slope_fits1[trough_resistance, on = c("scn", "variable", "met")] %>%
mutate(severity = paste0(severity, " %")) ->
wide_data
wide_data %>%
ggplot(aes(value,seg_slope1, color = severity)) +
geom_point(size = 3) +
facet_wrap('variable', scales = "free") +
theme_bw() +
THEME +
labs(x = "Resistance (initial response)",
y = "Resilience (rate of recovery)") +
scale_color_manual(values = FORTE_SEVERITY_COLORS) +
theme_bw(base_size = 20) +
theme(legend.position = "none")
Met & Metric Interactions
annual_data %>%
select(prate, tmp, met) %>%
distinct() %>%
dplyr::left_join(wide_data, by = "met") %>%
# dplyr::filter(variable %in% c("GPP", "Rh")) %>%
ggplot(aes(tmp, value, color = severity)) +
geom_point(size = 3) +
facet_wrap("variable", scales = "free") +
THEME +
scale_color_manual(values = FORTE_SEVERITY_COLORS) +
labs(x = "Average Annual Temperature (K)", y = "Resistance (unitless)") +
theme(legend.position = "none") +
theme_bw(base_size = 20) +
theme(legend.position = "none")
