This markdown script reproduces the supplementary figures associated with the paper “Atmospheric River Sequences as Indicators of Hydrologic Hazard in Historical Reanalysis and GFDL SPEAR Future Climate Projections” (https://doi.org/10.22541/essoar.167590838.86645781/v1).
Load packages & functions
## set random seed to produce consistent results
set.seed(2023)
source('_data/setup.R')
source('_scripts/create_df_functions.R')
## additional packages not in setup.R
library(patchwork)
library(pammtools)
library(lhs)
## progress bar
<- FALSE
progress
## geospatial helpers
<- colors.huc[c(6:10,1:5)]
colors.huc
## download custom boxplot function
library(assertthat)
library(devtools)
source_url('https://github.com/BoulderCodeHub/CRSSIO/blob/master/R/stat_boxplot_custom.R?raw=TRUE')
Load MERRA-2 catalogs
# ## MERRA 3-hour dataset
# load('_scripts/_checkpoints/df_3hr_1209.Rdata')
## MERRA 24-hour dataset
load('_scripts/_checkpoints/df_24hr_0209.Rdata')
## MERRA metadata
<- seq(ymd_hms('1980-1-1 00:00:00'), ymd_hms('2021-12-31 21:00:00'), by = '3 hours')
ts_merra <- unique(as.Date(ts_merra)) dates_merra
Load GFDL catalogs
## SPEAR historic
load('_scripts/_checkpoints/df_hist_1212.Rdata')
## SPEAR SSP 2-4.5
load('_scripts/_checkpoints/df_ssp245_1214.Rdata')
## SPEAR SSP 5-8.5
load('_scripts/_checkpoints/df_ssp585_1214.Rdata')
## SPEAR metadata
load('_data/GFDL/gfdl_metadata.Rdata') #ts_hist, ts_future
<- data.frame(ts = ts_future) %>%
ts_decades filter(year(ts) > 2020 & year(ts) <= 2090) %>%
pull(ts)
<-
decades data.frame(yr = 2021:2090) %>%
mutate(decade = ceiling(yr/10)-202)
<-
compare_biascorr map_dfr(.x = 1801:1810, .f = function(huc) {
<- which(grid_hucrep[] == huc)
i %>%
df_hist[[i]] filter(wateryear(date) %in% 1981:2010) %>%
select(date,ens,ivt) %>%
pivot_wider(names_from = ens, names_prefix = 'ens0', values_from = ivt) %>%
select(-date) %>%
mutate(
merra = df_24hr[[i]] %>%
filter(wateryear(date) %in% 1981:2010) %>%
pull(ivt.snap) %>% sort) %>%
apply(2, sort) %>% as.data.frame %>%
pivot_longer(starts_with('ens'), names_to = 'ens', values_to = 'gfdl') %>%
mutate(
huc4 = paste0('HUC', huc),
tag = c(paste0('(',letters[huc-1800],')'), rep(NA, nrow(.)-1)))
%>% filter(merra < 2000 & gfdl < 2000)
})
ggplot(compare_biascorr) +
geom_vline(xintercept = 0, size = 0.75) + geom_hline(yintercept = 0) +
geom_step(aes(x = merra, y = gfdl, group = ens, color = factor(ens)), direction = 'mid') +
geom_text(
aes(x = 250, y = 1200, label = tag),
family = 'Segoe UI', size = 10/.pt, fontface = 'bold') +
facet_wrap(~huc4, nrow = 2) +
scale_color_manual('Ensemble Member', labels = 1:5, values = paste0('grey', seq(20,80,15))) +
guides(color = guide_legend(override.aes = list(size = 1))) +
scale_x_origin('MERRA-2 IVT (kg/m/s)', breaks = 250*(0:10)) +
scale_y_origin('GFDL SPEAR Bias-Corrected IVT (kg/m/s)', breaks = 250*(0:10)) +
geom_parity() + coord_fixed() +
theme(
text = element_text(size = 10),
panel.grid.major = element_line(size = 0.25),
axis.line = element_blank(),
legend.position = 'bottom',
strip.background = element_rect(fill = 'grey95', size = NA),
strip.text = element_text(margin = margin(2,2,2,2)),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
ggsave('_figures/s1_biascorr.png', width = 6, height = 3.75, dpi = 600)
<-
compare_freq foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (i %in% index_ca) {
<- df_hist[[i]] %>%
temp filter(year(date) %in% 1981:2010 & month(date) %in% c(10:12,1:3)) %>%
count(ens,seq) %>% filter(seq)
<- df_24hr[[i]] %>%
merra filter(year(date) %in% 1981:2010 & month(date) %in% c(10:12,1:3)) %>%
count(seq) %>% filter(seq) %>% pull(n)
%>% mutate(merra = if (length(merra)==0) NA else merra) %>% mutate(id = i)
temp
}
}
<- map(
r2 .x = 1:5,
.f = ~summary(lm(merra ~ n, data = compare_freq %>% filter(ens==.x)))$r.squared) %>%
reduce(c)
%>%
compare_freq group_by(ens) %>%
mutate(
id = 1:length(id),
tag = case_when(id==1 ~ paste0('(', letters[ens], ')')),
rsq = case_when(id==1 ~ paste0('R^2','==', round(r2,3)[ens]))) %>%
ggplot() +
geom_vline(xintercept = 0, size = 0.75) + geom_hline(yintercept = 0) +
geom_point(aes(x = merra/30, y = n/30), size = 1) +
geom_text(
aes(x = 10, y = 81, label = tag),
family = 'Segoe UI', size = 10/.pt, fontface = 'bold') +
geom_text(
aes(x = 80, y = 20, label = rsq), hjust = 1, parse = TRUE,
family = 'Segoe UI', size = 9/.pt) +
facet_wrap(~paste('Ensemble Member',ens)) +
scale_x_origin(
'MERRA-2 Annual Average Sequence Days', breaks = 25*(0:5)) +
scale_y_origin(
'GFDL SPEAR Annual Average Sequence Days', breaks = 25*(0:5)) +
geom_abline(linetype = 'dashed', color = 'gray30') + coord_fixed() +
theme(
text = element_text(size = 10),
panel.grid.major = element_line(size = 0.25),
legend.position = 'bottom',
axis.line = element_blank(),
strip.background = element_rect(fill = 'grey95', size = NA),
strip.text = element_text(margin = margin(2,2,2,2)),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
ggsave('_figures/s2_biascorr_freq.png', width = 6, height = 4.75, dpi = 600)
<-
compare foreach (huc = 1801:1810) %do% {
<- which(grid_hucrep[] == huc)
i rbind(
%>%
df_hist[[i]] filter(seq) %>% group_by(seq.count) %>%
summarize(date = date[1], maxivt = seq.maxivt[1], duration = seq.duration[1], ens = ens[1]) %>%
filter(year(date) %in% 1981:2010 & month(date) %in% c(10:12,1:3)),
%>%
df_24hr[[i]] filter(seq) %>% group_by(seq.count) %>%
summarize(date = date[1], maxivt = seq.maxivt[1], duration = seq.duration[1], ens = NA) %>%
filter(year(date) %in% 1981:2010 & month(date) %in% c(10:12,1:3))) %>%
mutate(huc = huc)
%>%
} reduce(rbind) %>%
group_by(huc) %>%
mutate(
id = 1:length(huc),
tag = case_when(id==1 ~ paste0('(',letters[huc-1800],')'))) %>%
ungroup
%>%
compare group_by(ens,huc) %>%
arrange(maxivt) %>%
mutate(p = add_index(length(maxivt))) %>%
%>%
ungroup ggplot() +
geom_vline(xintercept = 250, size = 0.75) + geom_hline(yintercept = 0) +
geom_step(
aes(x = maxivt, y = p, group = ens, color = is.na(ens), size = is.na(ens))) +
geom_text(
aes(x = 380, y = 0.9, label = tag),
family = 'Segoe UI', size = 10/.pt, fontface = 'bold') +
facet_wrap(~paste0('HUC',huc), nrow = 2) +
scale_color_manual(
'Data Source',
values = c('grey60', 'black'),
labels = c('GFDL SPEAR', 'MERRA-2')) +
scale_size_manual(values = c(0.25, 0.5)) +
guides(
color = guide_legend(override.aes = list(size = 1)),
size = guide_none()) +
scale_x_continuous(
'Sequence Max IVT (kg/m/s)', breaks = 250*(0:10),
expand = expansion(mult = c(0,0.05))) +
scale_y_origin('Cumulative Probability') +
theme(
text = element_text(size = 10),
panel.grid.major = element_line(size = 0.25),
legend.position = 'bottom',
strip.background = element_rect(fill = 'grey95', size = NA),
strip.text = element_text(margin = margin(2,2,2,2)),
axis.line = element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
ggsave('_figures/s3_biascorr_maxivt.png', width = 6, height = 3.75, dpi = 600)
%>%
compare group_by(ens,huc) %>%
arrange(duration) %>%
mutate(p = add_index(length(duration))) %>%
%>%
ungroup ggplot() +
geom_vline(xintercept = 0, size = 0.75) + geom_hline(yintercept = 0) +
geom_step(
aes(x = duration, y = p, group = ens, color = is.na(ens), size = is.na(ens))) +
geom_text(
aes(x = 13, y = 0.9, label = tag),
family = 'Segoe UI', size = 10/.pt, fontface = 'bold') +
facet_wrap(~paste0('HUC',huc), nrow = 2) +
scale_color_manual(
'Data Source',
values = c('grey60', 'black'),
labels = c('GFDL SPEAR', 'MERRA-2')) +
scale_size_manual(values = c(0.25, 0.5)) +
guides(
color = guide_legend(override.aes = list(size = 1)),
size = guide_none()) +
scale_x_origin('Sequence Duration (days)', breaks = 30*(0:10)) +
scale_y_origin('Cumulative Probability') +
theme(
text = element_text(size = 10),
panel.grid.major = element_line(size = 0.25),
legend.position = 'bottom',
strip.background = element_rect(fill = 'grey95', size = NA),
strip.text = element_text(margin = margin(2,2,2,2)),
axis.line = element_blank(),
axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
ggsave('_figures/s4_biascorr_duration.png', width = 6, height = 3.75, dpi = 600)
## define big table of samples
<- improvedLHS(n = 100, k = 3)
samples.norm <-
samples.dist data.frame(
window = samples.norm[,1] %>% qunif(min = 3, max = 10) %>% round,
startend = samples.norm[,2] %>% qunif(min = 0.25, max = 0.75),
magcutoff = samples.norm[,3] %>% qunif(min = 200, max = 300)) %>%
mutate(magcutoff = cbind(magcutoff, startend) %>% apply(1, max))
<- samples.dist %>% rbind(c(window = 5, startend = 0.5, magcutoff = 250)) samples.dist
(this takes a while \(\to\) load from file)
# ## define progress bar
# if (progress) pb <- txtProgressBar(min = 0, max = 10*nrow(samples.dist), style = 3)
# ## historical
# timer <- Sys.time()
# clust <- parallel::makeCluster(cores)
# registerDoSNOW(clust)
# sens_hist <-
# foreach(
# huc = 1801:1810,
# .options.snow = if (progress) opts,
# .packages = c('raster', 'lubridate', 'foreach', 'tidyverse'),
# .export = c('Mean', 'lag', 'add_counter', 'create_catalog')) %:%
# foreach(samp = 1:nrow(samples.dist), .combine = 'rbind') %dopar% {
# ## create sequence catalog
# data <- df_hist[[which(grid_hucrep[]==huc)]] %>%
# filter(wateryear(date) %in% 1981:2010) %>%
# transmute(date, ivt, ens)
# sequences <-
# map(.x = 1:5, .f = function(x) {
# if (nrow(data %>% filter(ens==x)) > 0) {
# calculate_sequences(
# data %>% filter(ens==x), data %>% filter(ens==x),
# window = samples.dist$window[samp],
# startend = samples.dist$startend[samp],
# magcutoff = samples.dist$magcutoff[samp]) %>%
# mutate(ens=x)
# }
# }) %>% reduce(rbind)
# ## calculate sequence statistics
# sequences %>%
# group_by(ens) %>%
# mutate(
# nseq = length(ens),
# freq = nseq/30) %>%
# ungroup %>%
# summarize(
# maxivt.mean = mean(maxivt),
# maxivt.q95 = quantile(maxivt,0.95),
# duration.mean = mean(duration),
# duration.q95 = quantile(duration,0.95),
# nseq = nseq[1], freq = freq[1])
# }
# stopCluster(clust)
# if (progress) cat('\n')
# ## SSP 2-4.5
# clust <- parallel::makeCluster(cores)
# registerDoSNOW(clust)
# sens_ssp245 <-
# foreach(
# huc = 1801:1810,
# .options.snow = if (progress) opts,
# .packages = c('raster', 'lubridate', 'foreach', 'tidyverse'),
# .export = c('Mean', 'lag', 'add_counter', 'create_catalog')) %:%
# foreach(samp = 1:nrow(samples.dist), .combine = 'rbind') %dopar% {
# ## create sequence catalog
# data <- df_ssp245[[which(grid_hucrep[]==huc)]] %>%
# filter(wateryear(date) %in% 2061:2090) %>%
# transmute(date, ivt, ens)
# data_hist <- df_hist[[which(grid_hucrep[]==huc)]] %>%
# filter(wateryear(date) %in% 1981:2010) %>%
# transmute(date, ivt, ens)
# sequences <-
# map(.x = 1:5, .f = function(x) {
# if (nrow(data %>% filter(ens==x)) > 0) {
# calculate_sequences(
# data %>% filter(ens==x), data_hist %>% filter(ens==x),
# window = samples.dist$window[samp],
# startend = samples.dist$startend[samp],
# magcutoff = samples.dist$magcutoff[samp]) %>%
# mutate(ens=x)
# }
# }) %>% reduce(rbind)
# ## calculate sequence statistics
# sequences %>%
# group_by(ens) %>%
# mutate(
# nseq = length(ens),
# freq = nseq/30) %>%
# ungroup %>%
# summarize(
# maxivt.mean = mean(maxivt),
# maxivt.q95 = quantile(maxivt,0.95),
# duration.mean = mean(duration),
# duration.q95 = quantile(duration,0.95),
# nseq = nseq[1], freq = freq[1])
# }
# stopCluster(clust)
# if (progress) cat('\n')
# ## SSP 5-8.5
# clust <- parallel::makeCluster(cores)
# registerDoSNOW(clust)
# sens_ssp585 <-
# foreach(
# huc = 1801:1810,
# .options.snow = if (progress) opts,
# .packages = c('raster', 'lubridate', 'foreach', 'tidyverse'),
# .export = c('Mean', 'lag', 'add_counter', 'create_catalog')) %:%
# foreach(samp = 1:nrow(samples.dist), .combine = 'rbind') %dopar% {
# ## create sequence catalog
# data <- df_ssp585[[which(grid_hucrep[]==huc)]] %>%
# filter(wateryear(date) %in% 2061:2090) %>%
# transmute(date, ivt, ens)
# data_hist <- df_hist[[which(grid_hucrep[]==huc)]] %>%
# filter(wateryear(date) %in% 1981:2010) %>%
# transmute(date, ivt, ens)
# sequences <-
# map(.x = 1:5, .f = function(x) {
# if (nrow(data %>% filter(ens==x)) > 0) {
# calculate_sequences(
# data %>% filter(ens==x), data_hist %>% filter(ens==x),
# window = samples.dist$window[samp],
# startend = samples.dist$startend[samp],
# magcutoff = samples.dist$magcutoff[samp]) %>%
# mutate(ens=x)
# }
# }) %>% reduce(rbind)
# ## calculate sequence statistics
# sequences %>%
# group_by(ens) %>%
# mutate(
# nseq = length(ens),
# freq = nseq/30) %>%
# ungroup %>%
# summarize(
# maxivt.mean = mean(maxivt),
# maxivt.q95 = quantile(maxivt,0.95),
# duration.mean = mean(duration),
# duration.q95 = quantile(duration,0.95),
# nseq = nseq[1], freq = freq[1])
# }
# stopCluster(clust)
# if (progress) cat('\n')
# Sys.time() - timer
# ## save checkpoint
# save(samples.dist, sens_hist, sens_ssp245, sens_ssp585, file = '_results/sensitivity_checkpoint.Rdata')
## load checkpoint from file
load('_results/_checkpoints/sensitivity_checkpoint.Rdata')
<-
sens_hist lapply(1:10, function(i) sens_hist[[i]] %>% mutate(index = 1:nrow(.), huc = 1800+i)) %>%
reduce(rbind) %>%
mutate(freq = setNA(freq,0))
<-
sens_ssp245 lapply(1:10, function(i) sens_ssp245[[i]] %>% mutate(index = 1:nrow(.), huc = 1800+i)) %>%
reduce(rbind) %>%
mutate(freq = setNA(freq,0))
<-
sens_ssp585 lapply(1:10, function(i) sens_ssp585[[i]] %>% mutate(index = 1:nrow(.), huc = 1800+i)) %>%
reduce(rbind) %>%
mutate(freq = setNA(freq,0))
## outcome: absolute value
<-
sens rbind(
%>% mutate(source = 'hist'),
sens_hist %>% mutate(source = 'ssp245'),
sens_ssp245 %>% mutate(source = 'ssp585')) %>%
sens_ssp585 left_join(samples.dist %>% mutate(index = 1:nrow(.)), by = 'index')
## outcome: relative change
<-
sens_change_relative rbind(
/ sens_hist) %>% select(-huc, -index) %>%
(sens_ssp245 cbind(sens_hist %>% select(huc, index)) %>%
mutate(source = 'ssp245'),
/ sens_hist) %>% select(-huc, -index) %>%
(sens_ssp585 cbind(sens_hist %>% select(huc, index)) %>%
mutate(source = 'ssp585')) %>%
left_join(samples.dist %>% mutate(index = 1:nrow(.)), by = 'index')
## decide on an appropriate number of clusters
<- 3
k
## define custom kmeans function
<- function(x, ...) kmeans(x, centers = k, nstart = 50, ...)
Kmeans
## outcome: absolute value
<-
sens foreach (h = 1801:1810, .combine = 'rbind') %do% {
<- sens %>% filter(huc == h)
temp <-
good.freq !is.na(temp$freq) & !is.nan(temp$freq) & !is.infinite(temp$freq)
<-
good.maxivt !is.na(temp$maxivt.mean) & !is.nan(temp$maxivt.mean) & !is.infinite(temp$maxivt.mean)
<-
good.duration !is.na(temp$duration.mean) & !is.nan(temp$duration.mean) & !is.infinite(temp$duration.mean)
%>%
temp mutate(
cl.freq = ifelse(
Kmeans(freq[good.freq])$cluster, NA),
good.freq, cl.maxivt = ifelse(
Kmeans(maxivt.mean[good.maxivt])$cluster, NA),
good.maxivt, cl.duration = ifelse(
Kmeans(duration.mean[good.duration])$cluster, NA))
good.duration,
}
## outcome: relative change
<-
sens_change_relative foreach (h = 1801:1810, .combine = 'rbind') %do% {
<- sens_change_relative %>% filter(huc == h)
temp <-
good.freq !is.na(temp$freq) & !is.nan(temp$freq) & !is.infinite(temp$freq)
<-
good.maxivt !is.na(temp$maxivt.mean) & !is.nan(temp$maxivt.mean) & !is.infinite(temp$maxivt.mean)
<-
good.duration !is.na(temp$duration.mean) & !is.nan(temp$duration.mean) & !is.infinite(temp$duration.mean)
%>%
temp mutate(
cl.freq = ifelse(
Kmeans(freq[good.freq])$cluster, NA),
good.freq, cl.maxivt = ifelse(
Kmeans(maxivt.mean[good.maxivt])$cluster, NA),
good.maxivt, cl.duration = ifelse(
Kmeans(duration.mean[good.duration])$cluster, NA))
good.duration, }
## bootstrap L1-norm distribution
<- function(data, cl, var, boot = 1e3) {
L1_boot ## find existing parameters within clusters
<- sort(data[data$cl == cl,] %>% pull(get(var)))
x <- sort(data %>% pull(get(var)))
x.all <- (1:length(x))/length(x)
p <- (1:length(x.all))/length(x.all)
p.all <- interp1(x = p.all, y = x.all, xi = p)
x.interp
if (is.na(boot)) {
return(data.frame(x, x.interp) %>% as.matrix %>% apply(1, diff) %>% abs %>% sum)
else {
} <- c()
L1 for (b in 1:boot) {
<-
L1[b] data.frame(
x.boot = sort(sample(x.all, size = length(x), replace = TRUE)),
%>%
x.interp) %>% apply(1, diff) %>% abs %>% sum
as.matrix
}return(L1)
} }
Frequency:
if (progress) pb <- txtProgressBar(min = 0, max = 30*k, style = 3)
<- parallel::makeCluster(cores)
clust registerDoSNOW(clust)
<-
rsa_freq foreach (h = 1801:1810, .combine = 'rbind', .inorder = FALSE) %:%
foreach (cl = 1:k, .combine = 'rbind', .inorder = FALSE) %:%
foreach (
var = c('startend', 'window', 'magcutoff'),
.options.snow = if (progress) opts,
.export = 'L1_boot',
.packages = c('tidyverse', 'pracma'),
.combine = 'rbind', .inorder = FALSE) %dopar% {
## generate normalized L1 distances for every parameter & cluster
<- sens %>%
temp filter(huc==h & !is.na(cl.freq)) %>% rename(cl = cl.freq)
data.frame(
huc = h,
cl, var, d = L1_boot(temp, cl, var, NA)) %>%
cbind(
L1_boot(temp, cl, var) %>%
quantile(c(0.95, 0.99, 0.995)) %>%
%>% as.data.frame %>% setNames(c('d95', 'd99', 'd995')))
t %>%
} select(cl, var, d, d95, d99, d995, huc)
cat('\n')
<-
rsa_freq_rel foreach (h = 1801:1810, .combine = 'rbind', .inorder = FALSE) %:%
foreach (cl = 1:k, .combine = 'rbind', .inorder = FALSE) %:%
foreach (
var = c('startend', 'window', 'magcutoff'),
.options.snow = if (progress) opts,
.export = 'L1_boot',
.packages = c('tidyverse', 'pracma'),
.combine = 'rbind', .inorder = FALSE) %dopar% {
## generate normalized L1 distances for every parameter & cluster
<- sens_change_relative %>%
temp filter(huc==h & !is.na(cl.freq)) %>% rename(cl = cl.freq)
data.frame(
huc = h,
cl, var, d = L1_boot(temp, cl, var, NA)) %>%
cbind(
L1_boot(temp, cl, var) %>%
quantile(c(0.95, 0.99, 0.995)) %>%
%>% as.data.frame %>% setNames(c('d95', 'd99', 'd995')))
t %>%
} select(cl, var, d, d95, d99, d995, huc)
stopCluster(clust)
Intensity:
if (progress) pb <- txtProgressBar(min = 0, max = 30*k, style = 3)
<- parallel::makeCluster(cores)
clust registerDoSNOW(clust)
<-
rsa_maxivt foreach (h = 1801:1810, .combine = 'rbind', .inorder = FALSE) %:%
foreach (cl = 1:k, .combine = 'rbind', .inorder = FALSE) %:%
foreach (
var = c('startend', 'window', 'magcutoff'),
.options.snow = if (progress) opts,
.export = 'L1_boot',
.packages = c('tidyverse', 'pracma'),
.combine = 'rbind', .inorder = FALSE) %dopar% {
## generate normalized L1 distances for every parameter & cluster
<- sens %>%
temp filter(huc==h & !is.na(cl.maxivt)) %>% rename(cl = cl.maxivt)
data.frame(
huc = h,
cl, var, d = L1_boot(temp, cl, var, NA)) %>%
cbind(
L1_boot(temp, cl, var) %>%
quantile(c(0.95, 0.99, 0.995)) %>%
%>% as.data.frame %>% setNames(c('d95', 'd99', 'd995')))
t %>%
} select(cl, var, d, d95, d99, d995, huc)
cat('\n')
<-
rsa_maxivt_rel foreach (h = 1801:1810, .combine = 'rbind', .inorder = FALSE) %:%
foreach (cl = 1:k, .combine = 'rbind', .inorder = FALSE) %:%
foreach (
var = c('startend', 'window', 'magcutoff'),
.options.snow = if (progress) opts,
.export = 'L1_boot',
.packages = c('tidyverse', 'pracma'),
.combine = 'rbind', .inorder = FALSE) %dopar% {
## generate normalized L1 distances for every parameter & cluster
<- sens_change_relative %>%
temp filter(huc==h & !is.na(cl.maxivt)) %>% rename(cl = cl.maxivt)
data.frame(
huc = h,
cl, var, d = L1_boot(temp, cl, var, NA)) %>%
cbind(
L1_boot(temp, cl, var) %>%
quantile(c(0.95, 0.99, 0.995)) %>%
%>% as.data.frame %>% setNames(c('d95', 'd99', 'd995')))
t %>%
} select(cl, var, d, d95, d99, d995, huc)
stopCluster(clust)
Duration:
if (progress) pb <- txtProgressBar(min = 0, max = 30*k, style = 3)
<- parallel::makeCluster(cores)
clust registerDoSNOW(clust)
<-
rsa_duration foreach (h = 1801:1810, .combine = 'rbind', .inorder = FALSE) %:%
foreach (cl = 1:k, .combine = 'rbind', .inorder = FALSE) %:%
foreach (
var = c('startend', 'window', 'magcutoff'),
.options.snow = if (progress) opts,
.export = 'L1_boot',
.packages = c('tidyverse', 'pracma'),
.combine = 'rbind', .inorder = FALSE) %dopar% {
## generate normalized L1 distances for every parameter & cluster
<- sens %>%
temp filter(huc==h & !is.na(cl.duration)) %>% rename(cl = cl.duration)
data.frame(
huc = h,
cl, var, d = L1_boot(temp, cl, var, NA)) %>%
cbind(
L1_boot(temp, cl, var) %>%
quantile(c(0.95, 0.99, 0.995)) %>%
%>% as.data.frame %>% setNames(c('d95', 'd99', 'd995')))
t %>%
} select(cl, var, d, d95, d99, d995, huc)
cat('\n')
<-
rsa_duration_rel foreach (h = 1801:1810, .combine = 'rbind', .inorder = FALSE) %:%
foreach (cl = 1:k, .combine = 'rbind', .inorder = FALSE) %:%
foreach (
var = c('startend', 'window', 'magcutoff'),
.options.snow = if (progress) opts,
.export = 'L1_boot',
.packages = c('tidyverse', 'pracma'),
.combine = 'rbind', .inorder = FALSE) %dopar% {
print(c(h, cl, var))
## generate normalized L1 distances for every parameter & cluster
<- sens_change_relative %>%
temp filter(huc==h & !is.na(cl.duration)) %>% rename(cl = cl.duration)
data.frame(
huc = h,
cl, var, d = L1_boot(temp, cl, var, NA)) %>%
cbind(
L1_boot(temp, cl, var) %>%
quantile(c(0.95, 0.99, 0.995)) %>%
%>% as.data.frame %>% setNames(c('d95', 'd99', 'd995')))
t %>%
} select(cl, var, d, d95, d99, d995, huc)
stopCluster(clust)
Frequency:
<- ggplot(sens) +
g1.freq geom_smooth(
aes(x = jitter(window), y = freq, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
scale_x_continuous(breaks = 3:10, expand = expansion(mult = 0.025)) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text.x = element_blank())
<- ggplot(sens) +
g4.freq geom_smooth(
aes(x = startend, y = freq, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text = element_blank())
<- ggplot(sens) +
g7.freq geom_smooth(
aes(x = magcutoff, y = freq, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text = element_blank())
<-
lims rbind(
layer_scales(g1.freq)$y$get_limits(),
layer_scales(g4.freq)$y$get_limits(),
layer_scales(g7.freq)$y$get_limits())
<- c(apply(lims, 2, min)[1], apply(lims, 2, max)[2])
lims <- g1.freq + coord_cartesian(ylim = lims)
g1.freq <- g4.freq + coord_cartesian(ylim = lims)
g4.freq <- g7.freq + coord_cartesian(ylim = lims)
g7.freq
<- rsa_freq %>%
g1.rsa mutate(
var_name = case_when(
== 'window' ~ 'Moving\nAverage',
var == 'startend' ~ 'Start-End\nThreshold',
var == 'magcutoff' ~ 'Mag.\nFilter') %>%
var factor(levels = c('Moving\nAverage', 'Start-End\nThreshold', 'Mag.\nFilter'))) %>%
group_by(var_name, huc) %>%
summarize(d.norm = mean(d/d95), .groups = 'drop') %>%
ggplot() +
geom_col(
aes(
y = d.norm, x = var_name,
color = d.norm>1, alpha = d.norm>1,
group = factor(huc), fill = factor(huc)),
position = 'dodge', size = 0.25) +
geom_hline(yintercept = 1, linetype = 'dashed') +
scale_color_manual(values = c(NA, 'black'), na.value = NA, guide = guide_none()) +
scale_fill_manual('', values = colors.huc, guide = guide_none()) +
scale_alpha_manual(values = c(0.25, 0.75), guide = guide_none()) +
scale_y_origin('Sensitivity') +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title.x = element_blank(), legend.position = 'bottom')
<- ggplot(sens_change_relative) +
g3.freq geom_hline(yintercept = 0, color = 'grey50', size = 0.25) +
geom_smooth(
aes(x = jitter(window), y = freq-1, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
scale_x_continuous(breaks = 3:10, expand = expansion(mult = 0.025)) +
scale_y_continuous(labels = percent) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text.x = element_blank())
<- ggplot(sens_change_relative) +
g6.freq geom_hline(yintercept = 0, color = 'grey50', size = 0.25) +
geom_smooth(
aes(x = startend, y = freq-1, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
scale_y_continuous(labels = percent) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text = element_blank())
<- ggplot(sens_change_relative) +
g9.freq geom_hline(yintercept = 0, color = 'grey50', size = 0.25) +
geom_smooth(
aes(x = magcutoff, y = freq-1, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
scale_y_continuous(labels = percent) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text = element_blank())
<-
lims rbind(
layer_scales(g3.freq)$y$get_limits(),
layer_scales(g6.freq)$y$get_limits(),
layer_scales(g9.freq)$y$get_limits())
<- c(apply(lims, 2, min)[1], apply(lims, 2, max)[2])
lims <- g3.freq + coord_cartesian(ylim = lims)
g3.freq <- g6.freq + coord_cartesian(ylim = lims)
g6.freq <- g9.freq + coord_cartesian(ylim = lims)
g9.freq
<- rsa_freq_rel %>%
g2.rsa mutate(
var_name = case_when(
== 'window' ~ 'Moving\nAverage',
var == 'startend' ~ 'Start-End\nThreshold',
var == 'magcutoff' ~ 'Mag.\nFilter') %>%
var factor(levels = c('Moving\nAverage', 'Start-End\nThreshold', 'Mag.\nFilter'))) %>%
group_by(var_name, huc) %>%
summarize(d.norm = mean(d/d95), .groups = 'drop') %>%
ggplot() +
geom_col(
aes(
y = d.norm, x = var_name,
color = d.norm>1, alpha = d.norm>1,
group = factor(huc), fill = factor(huc)),
position = 'dodge', size = 0.25) +
geom_hline(yintercept = 1, linetype = 'dashed') +
scale_color_manual(values = c(NA, 'black'), na.value = NA, guide = guide_none()) +
scale_fill_manual('', values = colors.huc, guide = guide_none()) +
scale_alpha_manual(values = c(0.25, 0.75), guide = guide_none()) +
scale_y_origin('Sensitivity') +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title.x = element_blank(), legend.position = 'bottom')
Intensity:
<- ggplot(sens) +
g1.maxivt geom_smooth(
aes(x = jitter(window), y = maxivt.mean, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
scale_x_continuous(breaks = 3:10, expand = expansion(mult = 0.025)) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text.x = element_blank())
<- ggplot(sens) +
g4.maxivt geom_smooth(
aes(x = startend, y = maxivt.mean, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text = element_blank())
<- ggplot(sens) +
g7.maxivt geom_smooth(
aes(x = magcutoff, y = maxivt.mean, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text = element_blank())
<-
lims rbind(
layer_scales(g1.maxivt)$y$get_limits(),
layer_scales(g4.maxivt)$y$get_limits(),
layer_scales(g7.maxivt)$y$get_limits())
<- c(apply(lims, 2, min)[1], apply(lims, 2, max)[2])
lims <- g1.maxivt + coord_cartesian(ylim = lims)
g1.maxivt <- g4.maxivt + coord_cartesian(ylim = lims)
g4.maxivt <- g7.maxivt + coord_cartesian(ylim = lims)
g7.maxivt
<- rsa_maxivt %>%
g3.rsa mutate(
var_name = case_when(
== 'window' ~ 'Moving\nAverage',
var == 'startend' ~ 'Start-End\nThreshold',
var == 'magcutoff' ~ 'Mag.\nFilter') %>%
var factor(levels = c('Moving\nAverage', 'Start-End\nThreshold', 'Mag.\nFilter'))) %>%
group_by(var_name, huc) %>%
summarize(d.norm = mean(d/d95), .groups = 'drop') %>%
ggplot() +
geom_col(
aes(
y = d.norm, x = var_name,
color = d.norm>1, alpha = d.norm>1,
group = factor(huc), fill = factor(huc)),
position = 'dodge', size = 0.25) +
geom_hline(yintercept = 1, linetype = 'dashed') +
scale_color_manual(values = c(NA, 'black'), na.value = NA, guide = guide_none()) +
scale_fill_manual('', values = colors.huc, guide = guide_none()) +
scale_alpha_manual(values = c(0.25, 0.75), guide = guide_none()) +
scale_y_origin('Sensitivity') +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title.x = element_blank(), legend.position = 'bottom')
<- ggplot(sens_change_relative) +
g3.maxivt geom_hline(yintercept = 0, color = 'grey50', size = 0.25) +
geom_smooth(
aes(x = jitter(window), y = maxivt.mean-1, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
scale_x_continuous(breaks = 3:10, expand = expansion(mult = 0.025)) +
scale_y_continuous(labels = percent, limits = c(0,1)) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text.x = element_blank())
<- ggplot(sens_change_relative) +
g6.maxivt geom_hline(yintercept = 0, color = 'grey50', size = 0.25) +
geom_smooth(
aes(x = startend, y = maxivt.mean-1, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
scale_y_continuous(labels = percent, limits = c(0,1)) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text = element_blank())
<- ggplot(sens_change_relative) +
g9.maxivt geom_hline(yintercept = 0, color = 'grey50', size = 0.25) +
geom_smooth(
aes(x = magcutoff, y = maxivt.mean-1, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
scale_y_continuous(labels = percent, limits = c(0,1)) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text = element_blank())
<- rsa_maxivt_rel %>%
g4.rsa mutate(
var_name = case_when(
== 'window' ~ 'Moving\nAverage',
var == 'startend' ~ 'Start-End\nThreshold',
var == 'magcutoff' ~ 'Mag.\nFilter') %>%
var factor(levels = c('Moving\nAverage', 'Start-End\nThreshold', 'Mag.\nFilter'))) %>%
group_by(var_name, huc) %>%
summarize(d.norm = mean(d/d95), .groups = 'drop') %>%
ggplot() +
geom_col(
aes(
y = d.norm, x = var_name,
color = d.norm>1, alpha = d.norm>1,
group = factor(huc), fill = factor(huc)),
position = 'dodge', size = 0.25) +
geom_hline(yintercept = 1, linetype = 'dashed') +
scale_color_manual(values = c(NA, 'black'), na.value = NA, guide = guide_none()) +
scale_fill_manual('', values = colors.huc, guide = guide_none()) +
scale_alpha_manual(values = c(0.25, 0.75), guide = guide_none()) +
scale_y_origin('Sensitivity', breaks = 0:3) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title.x = element_blank(), legend.position = 'bottom')
Duration:
<- ggplot(sens) +
g1.duration geom_smooth(
aes(x = jitter(window), y = duration.mean, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
scale_x_continuous(breaks = 3:10, expand = expansion(mult = 0.025)) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text.x = element_blank())
<- ggplot(sens) +
g4.duration geom_smooth(
aes(x = startend, y = duration.mean, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text = element_blank())
<- ggplot(sens) +
g7.duration geom_smooth(
aes(x = magcutoff, y = duration.mean, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text = element_blank())
<-
lims rbind(
layer_scales(g1.duration)$y$get_limits(),
layer_scales(g4.duration)$y$get_limits(),
layer_scales(g7.duration)$y$get_limits())
<- c(apply(lims, 2, min)[1], apply(lims, 2, max)[2])
lims <- g1.duration + coord_cartesian(ylim = lims)
g1.duration <- g4.duration + coord_cartesian(ylim = lims)
g4.duration <- g7.duration + coord_cartesian(ylim = lims)
g7.duration
<- rsa_duration %>%
g5.rsa mutate(
var_name = case_when(
== 'window' ~ 'Moving\nAverage',
var == 'startend' ~ 'Start-End\nThreshold',
var == 'magcutoff' ~ 'Mag.\nFilter') %>%
var factor(levels = c('Moving\nAverage', 'Start-End\nThreshold', 'Mag.\nFilter'))) %>%
group_by(var_name, huc) %>%
summarize(d.norm = mean(d/d95), .groups = 'drop') %>%
ggplot() +
geom_col(
aes(
y = d.norm, x = var_name,
color = d.norm>1, alpha = d.norm>1,
group = factor(huc), fill = factor(huc)),
position = 'dodge', size = 0.25) +
geom_hline(yintercept = 1, linetype = 'dashed') +
scale_color_manual(values = c(NA, 'black'), na.value = NA, guide = guide_none()) +
scale_fill_manual('', values = colors.huc, guide = guide_none()) +
scale_alpha_manual(values = c(0.25, 0.75), guide = guide_none()) +
scale_y_origin('Sensitivity') +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title.x = element_blank(), legend.position = 'bottom')
<- ggplot(sens_change_relative) +
g3.duration geom_hline(yintercept = 0, color = 'grey50', size = 0.25) +
geom_smooth(
aes(x = jitter(window), y = duration.mean-1, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
scale_x_continuous(breaks = 3:10, expand = expansion(mult = 0.025)) +
scale_y_continuous(labels = percent, limits = c(0,1)) +
theme(text = element_text(family = 'Segoe UI', size = 8), axis.title = element_blank())
<- ggplot(sens_change_relative) +
g6.duration geom_hline(yintercept = 0, color = 'grey50', size = 0.25) +
geom_smooth(
aes(x = startend, y = duration.mean-1, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
scale_y_continuous(labels = percent, limits = c(0,1)) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text.y = element_blank())
<- ggplot(sens_change_relative) +
g9.duration geom_hline(yintercept = 0, color = 'grey50', size = 0.25) +
geom_smooth(
aes(x = magcutoff, y = duration.mean-1, group = huc, color = factor(huc)),
size = 0.75, se = FALSE) +
scale_color_manual('', values = colors.huc, guide = guide_none()) +
scale_y_continuous(labels = percent, limits = c(0,1)) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title = element_blank(), axis.text.y = element_blank())
<- rsa_duration_rel %>%
g6.rsa mutate(
var_name = case_when(
== 'window' ~ 'Moving\nAverage',
var == 'startend' ~ 'Start-End\nThreshold',
var == 'magcutoff' ~ 'Mag.\nFilter') %>%
var factor(levels = c('Moving\nAverage', 'Start-End\nThreshold', 'Mag.\nFilter'))) %>%
group_by(var_name, huc) %>%
summarize(d.norm = mean(d/d95), .groups = 'drop') %>%
ggplot() +
geom_col(
aes(
y = d.norm, x = var_name,
color = d.norm>1, alpha = d.norm>1,
group = factor(huc), fill = factor(huc)),
position = 'dodge', size = 0.25) +
geom_hline(yintercept = 1, linetype = 'dashed') +
scale_color_manual(values = c(NA, 'black'), na.value = NA, guide = guide_none()) +
scale_fill_manual('', values = colors.huc, guide = guide_none()) +
scale_alpha_manual(values = c(0.25, 0.75), guide = guide_none()) +
scale_y_origin('Sensitivity') +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.title.x = element_blank(), legend.position = 'bottom')
Note: Additional labels were added to this figure in Inkscape.
ggplot(sens) +
geom_histogram(aes(x = huc, fill = factor(huc))) +
scale_fill_manual(
'Hydrologic Region',
values = colors.huc, guide = guide_legend(direction = 'horizontal', nrow = 1)) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
legend.position = 'bottom',
legend.text = element_text(margin = margin(-2,-2,-2,-2)),
legend.margin = margin(-2,-2,-2,-2),
legend.box.margin = margin(-2,-2,-2,-2))
ggsave('_figures/s5_legend.png', width = 5.5, height = 1, dpi = 600)
+ g4.freq + g7.freq + g1.rsa + theme(axis.text.x = element_blank()) +
g1.freq + g6.freq + g9.freq + g2.rsa + theme(axis.text.x = element_blank()) +
g3.freq + g4.maxivt + g7.maxivt + g3.rsa + theme(axis.text.x = element_blank()) +
g1.maxivt + g6.maxivt + g9.maxivt + g4.rsa + theme(axis.text.x = element_blank()) +
g3.maxivt + g4.duration + g7.duration + g5.rsa + theme(axis.text.x = element_blank()) +
g1.duration + g6.duration + g9.duration + g6.rsa +
g3.duration theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
plot_layout(design = 'abcd\nefgh\nijkl\nmnop\nqrst\nuvwx') +
plot_annotation(tag_levels = 'a', tag_prefix = '(', tag_suffix = ')') &
theme(plot.tag = element_text(
family = 'Segoe UI', size = 10, face = 'bold', margin = margin(-2,0,-2,-10)))
ggsave('_figures/s5.png', width = 5.5, height = 7.75, dpi = 600)
<-
ndays_ref foreach (i = 1:ncell(grid_ca), .combine = 'c') %do% {
if (progress) setTxtProgressBar(pb,i)
if (i %in% index_ca) {
%>%
df_hist[[i]] filter(wateryear(date) %in% 1981:2010) %>%
filter(month(date) %in% c(10:12,1:3)) %>%
filter(seq) %>%
count(ens) %>%
summarize(nn = sum(n)/length(unique(ens))/30) %>%
pull(nn)
else NA
}
}cat('\n')
<-
ndays_ssp245 foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (is.null(df_ssp245[[i]])) rep(NA,7) else {
%>%
df_ssp245[[i]] mutate(yr = year(date)) %>%
right_join(decades, by = 'yr') %>%
filter(month(date) %in% c(10:12,1:3)) %>%
filter(seq) %>%
count(ens, decade) %>%
group_by(decade) %>%
summarize(nn = sum(n)/length(unique(ens))/10) %>%
right_join(data.frame(decade = 1:7), by = 'decade') %>%
arrange(decade) %>%
mutate(nn = setNA(nn,0)) %>%
pull(nn)
}
}cat('\n')
<-
ndays_ssp585 foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (is.null(df_ssp585[[i]])) rep(NA,7) else {
%>%
df_ssp585[[i]] mutate(yr = year(date)) %>%
right_join(decades, by = 'yr') %>%
filter(month(date) %in% c(10:12,1:3)) %>%
filter(seq) %>%
count(ens, decade) %>%
group_by(decade) %>%
summarize(nn = sum(n)/length(unique(ens))/10) %>%
right_join(data.frame(decade = 1:7), by = 'decade') %>%
arrange(decade) %>%
mutate(nn = setNA(nn,0)) %>%
pull(nn)
}
}cat('\n')
<-
ndays left_join(
%>% as.data.frame %>%
ndays_ssp245 cbind(hist = ndays_ref) %>%
cbind(huc = grid_hucrep[]) %>%
filter(huc > 0) %>%
pivot_longer(cols = starts_with('V'), names_to = 'decade') %>%
rename(ssp245 = value),
%>% as.data.frame %>%
ndays_ssp585 cbind(huc = grid_hucrep[]) %>%
filter(huc > 0) %>%
pivot_longer(cols = starts_with('V'), names_to = 'decade') %>%
rename(ssp585 = value),
by = c('huc', 'decade')) %>%
mutate(
decade = toNumber(gsub('V', '', decade)),
decade_name = case_when(
==1 ~ '2020s',
decade==2 ~ '2030s',
decade==3 ~ '2040s',
decade==4 ~ '2050s',
decade==5 ~ '2060s',
decade==6 ~ '2070s',
decade==7 ~ '2080s') %>% factor) %>%
decadepivot_longer(cols = c(ssp245, ssp585), names_to = 'sim') %>%
mutate(
sim_name = case_when(
== 'ssp245' ~ 'SSP 2-4.5 ("Middle of the Road")',
sim == 'ssp585' ~ 'SSP 5-8.5 ("Upper End")') %>%
sim %>% factor(levels = levels(.))) %>%
factor mutate(
id = 1:nrow(.),
tag = ifelse(id %in% 1:2, paste0('(', letters[id], ')'), NA))
<- ggplot(ndays) +
gdays geom_hline(yintercept = 0, color = 'grey50', size = 0.35) +
geom_line(
aes(x = decade_name, y = value-hist, group = huc, color = factor(huc)),
size = 0.75) +
geom_point(
aes(x = decade_name, y = value-hist, group = huc, color = factor(huc)),
size = 1) +
geom_text(
aes(x = 1.5, y = 40, label = tag),
size = 10/.pt, family = 'Segoe UI', fontface = 'bold') +
facet_wrap(~sim_name, nrow = 2) +
scale_color_manual(values = colors.huc, guide = guide_none()) +
scale_x_discrete(expand = expansion(mult = 0.035)) +
scale_y_continuous(
'Change in Number of Sequence Days\nRelative to WY 1981-2010',
expand = expansion(mult = 0.1)) +
theme(
text = element_text(family = 'Segoe UI', size = 10),
panel.grid.major.y = element_line(size = 0.25),
strip.background = element_rect(color = NA, fill = 'grey95'),
strip.text = element_text(size = 10, margin = margin(2,0,2,0)),
legend.position = 'bottom',
axis.line.x = element_blank(),
axis.title.x = element_blank())
<- ggplot(huc4) +
gmap geom_sf(aes(fill = factor(huc4)), alpha = 0.5, size = 0.25) +
scale_fill_manual('Hydrologic Region', values = colors.huc) +
geom_tile(
data = raster.df(grid_ca) %>% filter(!is.na(value)),
aes(x=x, y=y), color = 'grey70', fill = NA) +
geom_tile(
data = raster.df(grid_hucrep) %>% filter(value>0),
aes(x=x, y=y, fill = factor(value)), color = 'black') +
geom_point(
data = raster.df(grid_hucrep) %>% filter(value>0),
aes(x=x, y=y), color = 'black', size = 0.5) +
guides(fill = guide_legend(nrow = 1)) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
panel.border = element_rect(fill = NA, size = 0.25),
axis.line = element_blank(),
axis.text.x = element_text(angle = 90, vjust = 0.5),
axis.title = element_blank(),
legend.position = 'bottom',
legend.background = element_rect(fill = NA),
legend.margin = margin(c(0,0,0,0)))
+ get_legend(gmap) + plot_layout(heights = c(10,1))
gdays ggsave('_figures/s6_sequencedays.png', width = 6, height = 4.5, dpi = 600)
<-
stats_all map(.x = 1801:1810, .f = function(huc) {
rbind(
which(grid_hucrep[]==huc)]] %>%
df_hist[[filter(seq) %>%
filter(year(date) %in% 1981:2010) %>%
group_by(ens,seq.count) %>%
summarize(maxivt = seq.maxivt[1], duration = seq.duration[1], .groups = 'drop') %>%
mutate(source = 'hist'),
which(grid_hucrep[]==huc)]] %>%
df_ssp245[[filter(seq) %>%
group_by(ens,seq.count) %>%
summarize(maxivt = seq.maxivt[1], duration = seq.duration[1], .groups = 'drop') %>%
mutate(source = 'ssp245'),
which(grid_hucrep[]==huc)]] %>%
df_ssp585[[filter(seq) %>%
group_by(ens,seq.count) %>%
summarize(maxivt = seq.maxivt[1], duration = seq.duration[1], .groups = 'drop') %>%
mutate(source = 'ssp585')) %>%
mutate(
source_name = case_when(
== 'hist' ~ 'Historic',
source == 'ssp245' ~ 'SSP 2-4.5',
source == 'ssp585' ~ 'SSP 5-8.5') %>%
source %>% factor(levels = levels(.)[3:1]),
factor huc = huc,
huc4 = paste0('HUC',huc))
%>% reduce(rbind)
})
## add panel labels
<- stats_all %>%
stats_all group_by(huc) %>%
mutate(id = 1:length(huc)) %>%
%>%
ungroup mutate(tag = ifelse(id==1, paste0('(',letters[huc-1800],')'), NA))
ggplot() +
# stat_boxplot_custom(
# aes(x = maxivt, y = source_name, group = paste(source,ens), fill = source),
# qs = c(.05, .25, .5, .75, .95),
# outlier.size = 0.5, size = 0.25, outlier.alpha = 1, alpha = 0.5,
# show.legend = FALSE) +
geom_boxplot(
data = stats_all %>%
mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>%
group_by(source_name,source,ens,huc4) %>%
summarize(
y05 = quantile(maxivt,0.05),
y25 = quantile(maxivt,0.25),
y50 = quantile(maxivt,0.50),
y75 = quantile(maxivt,0.75),
y95 = quantile(maxivt,0.95)),
aes(
xmin = y05, xlower = y25, xmiddle = y50, xupper = y75, xmax = y95,
y = source_name, group = paste(source_name,ens), fill = source_name),
stat = 'identity', position = 'dodge',
width = 0.8, size = 0.25, alpha = 0.5, show.legend = FALSE) +
geom_point(
data = stats_all %>%
mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>%
group_by(source_name,source,ens,huc4) %>%
mutate(
y05 = quantile(maxivt,0.05),
y95 = quantile(maxivt,0.95)) %>%
%>% group_by(source) %>%
ungroup mutate(
dodge = seq(-0.32,0.32,length.out = length(unique(ens)))[ens]) %>%
%>%
ungroup filter(maxivt < y05 | maxivt > y95),
aes(x = maxivt, y = as.numeric(factor(source_name))+dodge),
size = 0.5) +
geom_text(
data = stats_all,
aes(x = 1650, y = 3.25, label = tag),
family = 'Segoe UI', size = 10/.pt, fontface = 'bold') +
facet_wrap(~huc4, ncol = 2, dir = 'v', strip.position = 'left') +
scale_fill_manual(values = scico(7, palette = 'roma')) +
scale_x_continuous('Maximum 24-Hour IVT (kg/m/s)', breaks = 250*(0:10)) +
geom_hline(yintercept = 0.4) +
ggtitle('Sequence Intensity') +
theme(
text = element_text(family = 'Segoe UI', size = 10),
panel.grid.major.x = element_line(size = 0.25),
strip.background = element_rect(color = NA, fill = 'grey95'),
strip.text = element_text(margin = margin(2,2,2,2)),
strip.placement = 'inside',
axis.title.y = element_blank(),
axis.line = element_line(color = NA))
ggsave('_figures/s7_intensity.png', width = 6, height = 7, dpi = 600)
ggplot(stats_all) +
# stat_boxplot_custom(
# aes(x = duration, y = source_name, group = paste(source,ens), fill = source),
# qs = c(.05, .25, .5, .75, .95),
# outlier.size = 0.5, size = 0.25, outlier.alpha = 1, alpha = 0.5,
# show.legend = FALSE) +
geom_boxplot(
data = stats_all %>%
mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>%
group_by(source_name,source,ens,huc4) %>%
summarize(
y05 = quantile(duration,0.05),
y25 = quantile(duration,0.25),
y50 = quantile(duration,0.50),
y75 = quantile(duration,0.75),
y95 = quantile(duration,0.95)),
aes(
xmin = y05, xlower = y25, xmiddle = y50, xupper = y75, xmax = y95,
y = source_name, group = paste(source_name,ens), fill = source_name),
stat = 'identity', position = 'dodge',
width = 0.8, size = 0.25, alpha = 0.5, show.legend = FALSE) +
geom_point(
data = stats_all %>%
mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>%
group_by(source_name,source,ens,huc4) %>%
mutate(
y05 = quantile(duration,0.05),
y95 = quantile(duration,0.95)) %>%
%>% group_by(source) %>%
ungroup mutate(
dodge = seq(-0.32,0.32,length.out = length(unique(ens)))[ens]) %>%
%>%
ungroup filter(duration < y05 | duration > y95),
aes(x = duration, y = as.numeric(factor(source_name))+dodge),
size = 0.5) +
geom_text(
aes(x = 140, y = 3.25, label = tag),
family = 'Segoe UI', size = 10/.pt, fontface = 'bold') +
facet_wrap(~huc4, ncol = 2, dir = 'v', strip.position = 'left') +
scale_fill_manual(values = scico(7, palette = 'roma')[7:5]) +
scale_x_continuous('Duration (days)', breaks = 30*(0:10)) +
geom_hline(yintercept = 0.4) +
ggtitle('Sequence Duration') +
theme(
text = element_text(family = 'Segoe UI', size = 10),
panel.grid.major.x = element_line(size = 0.25),
strip.background = element_rect(color = NA, fill = 'grey95'),
strip.text = element_text(margin = margin(2,2,2,2)),
strip.placement = 'inside',
axis.title.y = element_blank(),
axis.line = element_line(color = NA))
ggsave('_figures/s8_duration.png', width = 6, height = 7, dpi = 600)