This markdown script reproduces the figures and numerical results 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), excepting Figure 2. Please note that the figures in the markdown file may not exactly match the published version because of formatting and scale differences.
Load packages & functions
source('_data/setup.R')
source('_scripts/create_df_functions.R')
## additional packages not in setup.R
library(patchwork)
library(pammtools)
library(assertthat)
library(devtools)
library(lubridate)
## progress bar
<- FALSE
progress if (progress) pb <- txtProgressBar(min = 0, max = ncell(grid_ca), style = 3)
## geospatial helpers
<- colors.huc[c(6:10,1:5)]
colors.huc <- st_union(california)
cal
## download custom boxplot function for Figures XX & XX
source_url('https://github.com/BoulderCodeHub/CRSSIO/blob/master/R/stat_boxplot_custom.R?raw=TRUE')
Load MERRA-2 catalogs
## 3-hour dataset
load('_scripts/_checkpoints/df_3hr_1209.Rdata')
## 24-hour dataset
load('_scripts/_checkpoints/df_24hr_0209.Rdata')
## date 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)
Start-end parameter for example sequence timeseries:
<-
startend 167]] %>%
df_3hr[[filter(wateryear(ts) %in% 1981:2010) %>%
pull(ivt) %>% median(na.rm = TRUE)
cat(paste('Example start-end parameter value:', round(startend,1), 'kg/m/s'))
## Example start-end parameter value: 91.9 kg/m/s
## pull example data from historic record
<- df_3hr[[167]] %>%
df_example filter(wateryear(ts) == 2016) %>%
filter(month(ts) %in% 10:12) %>%
filter(ts >= ymd_hms('2015-10-7 00:00:00')) %>%
mutate(threshold1 = cbind(rolling5, startend) %>% apply(1, min)) %>%
mutate(threshold2 = ifelse(seq.maxrolling > 250, rolling5, threshold1))
## create timeseries plot
<- ggplot(df_example) +
g geom_hline(yintercept = 250, color = 'grey70') +
geom_text(
data = data.frame(x = max(df_example$ts), y = 250),
aes(x=x, y=y, label = '250 kg/m/s'),
# nudge_x = 3.1e5, nudge_y = 30,
hjust = 1, vjust = -0.5, nudge_x = 8.5e5,
family = 'Segoe UI', size = 8/.pt,fontface = 'italic', color = 'grey30') +
geom_hline(yintercept = startend, color = 'grey70', linetype = 'dashed') +
geom_text(
data = data.frame(x = max(df_example$ts), y = 92),
aes(x=x, y=y, label = '92 kg/m/s'),
# nudge_x = 3.6e5, nudge_y = 30,
hjust = 1, vjust = -0.5, nudge_x = 8.5e5,
family = 'Segoe UI', size = 8/.pt, fontface = 'italic', color = 'grey30') +
geom_ribbon(
aes(x = ts, ymin = threshold1, ymax = rolling5),
fill = roma.colors[4], alpha = 0.25) +
geom_ribbon(
aes(x = ts, ymin = threshold1, ymax = threshold2),
fill = roma.colors[4], alpha = 0.75) +
geom_line(aes(x = ts, y = ivt, color = '3-Hour IVT')) +
geom_line(aes(x = ts, y = rolling5, color = '5-Day Moving Average IVT'), size = 0.75) +
scale_color_manual('', values = c('grey50', 'black')) +
scale_x_datetime(
'Time (weeks)', date_breaks = '1 week',
limits = c(min(df_example$ts), max(df_example$ts)+days(10)),
expand = expansion(mult = c(0.025, 0.01))) +
scale_y_origin('IVT (kg/m/s)', breaks = 250*(0:10)) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
legend.text = element_text(size = 8),
axis.text.x = element_blank(),
legend.position = c(0.1,0.8), legend.background = element_blank(),
panel.grid.major.y = element_line(size = 0.25))
## create bar chart of AR days
<- ggplot(df_example) +
g.ar ::geom_stepribbon(
pammtoolsaes(x = ts, ymin = 0, ymax = as.numeric(ar)),
fill = roma.colors[5]) +
geom_hline(yintercept = 0, color = roma.colors[5]) +
geom_text(
data = data.frame(x = min(df_example$ts), y = 0),
aes(x=x, y=y, label = 'AR'),
hjust = 0, vjust = -0.5,
family = 'Segoe UI', size = 8/.pt, fontface = 'plain') +
scale_x_datetime(
limits = c(min(df_example$ts), max(df_example$ts)+days(10)),
expand = expansion(mult = c(0.025, 0.01))) +
theme_void()
## create bar chart of sequence days
<- ggplot(df_example) +
g.seq ::geom_stepribbon(
pammtoolsaes(x = ts, ymin = 0, ymax = as.numeric(seq)),
fill = roma.colors[4]) +
geom_hline(yintercept = 0, color = roma.colors[4]) +
geom_hline(yintercept = 1.5, color = 'white') +
geom_text(
data = data.frame(x = min(df_example$ts), y = 0),
aes(x=x, y=y, label = 'Sequence'),
hjust = 0, vjust = -0.5,
family = 'Segoe UI', size = 8/.pt, fontface = 'plain') +
scale_x_datetime(
limits = c(min(df_example$ts), max(df_example$ts)+days(10)),
expand = expansion(mult = c(0.025, 0.01))) +
theme_void()
## combine into one plot and save
ggarrange(g, g.ar, g.seq, nrow = 3, heights = c(9,1,1.5), align = 'v')
ggsave('_figures/fig1_sequence.png', width = 6, height = 2, units = 'in', dpi = 600)
Statewide average annual sequence frequency:
<-
nseq_ref foreach (i = 1:ncell(grid_ca), .combine = 'c') %do% {
if (progress) setTxtProgressBar(pb,i)
if (is.null(df_hist[[i]])) NA else {
<- df_hist[[i]] %>%
temp filter(wateryear(date) %in% 1981:2010) %>%
filter(seq) %>%
group_by(ens, seq.count) %>%
summarize(duration = seq.duration[1], .groups = 'drop')
nrow(temp)/length(unique(temp$ens))
}
}<- nseq_ref / 30
nseq_ref if (progress) cat('\n')
cat(paste(
'Mean annual frequency:',
Mean(nseq_ref) %>% round(2), 'sequences/year \n'))
cat(paste(
'Standard deviation:',
sd(nseq_ref, na.rm = TRUE) %>% round(2), 'sequences/year \n'))
## Mean annual frequency: 1.47 sequences/year
## Standard deviation: 1.06 sequences/year
Statewide percentage of years with no sequences:
<-
nzero foreach (i = 1:ncell(grid_ca), .combine = 'c') %do% {
if (progress) setTxtProgressBar(pb,i)
if (is.null(df_hist[[i]])) NA else {
<- df_hist[[i]] %>%
temp filter(wateryear(date) %in% 1981:2010) %>%
filter(seq) %>%
group_by(ens, seq.count) %>%
summarize(wy = wateryear(date[1]), duration = seq.duration[1], .groups = 'drop') %>%
count(wy,ens) %>%
right_join(
expand.grid(wy = 1981:2010, ens = unique(temp$ens)),
by = c('wy','ens')) %>%
mutate(n = setNA(n,0))
%>%
temp group_by(ens) %>%
summarize(zero = sum(n==0)) %>%
pull(zero) %>% mean
}
}if (progress) cat('\n')
cat(paste(
'Mean number of years with no sequences:',
Mean(nzero) %>% round(1), '\n'))
cat(paste(
'Mean percentage of years with no sequences:',
Mean(nzero/30) %>% percent(accuracy = 0.1)))
## Mean number of years with no sequences: 9.7
## Mean percentage of years with no sequences: 32.5%
Statewide average ARs per sequence in the historic record:
<-
nar_ref foreach (i = 1:ncell(grid_ca), .combine = 'c') %do% {
if (progress) setTxtProgressBar(pb,i)
if (!is.null(df_hist[[i]])) {
%>%
df_hist[[i]] filter(seq) %>%
filter(wateryear(date) %in% 1981:2010) %>%
group_by(ens,seq.count) %>%
summarize(nar = length(unique(ar.count[!is.na(ar.count)])), .groups = 'drop') %>%
summarize(nar.mean = Mean(nar)) %>% unlist %>% unname
else NA
}
}if (progress) cat('\n')
cat(paste(
'Mean number of ARs per sequence:',
Mean(nar_ref) %>% round(1), '\n'))
cat(paste(
'Standard deviation of ARs per sequence:',
sd(nar_ref, na.rm = TRUE) %>% round(2), '\n'))
## Mean number of ARs per sequence: 2.4
## Standard deviation of ARs per sequence: 0.24
<-
historical foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (i %in% index_ca) {
%>%
df_hist[[i]] filter(wateryear(date) %in% 1981:2010) %>%
filter(seq) %>%
group_by(seq.count) %>%
summarize(
maxivt = seq.maxivt[1], duration = seq.duration[1], ens = ens[1]) %>%
summarize(
maxivt = mean(maxivt), duration = mean(duration),
n.seq = length(seq.count)/length(unique(ens))/30) %>%
cbind(
%>%
df_hist[[i]] filter(wateryear(date) %in% 1981:2010) %>%
filter(ar) %>%
group_by(ar.count) %>%
summarize(pct.ar = mean(seq)) %>%
summarize(pct.ar = mean(pct.ar)))
else rep(NA,4)
}
}<- historical %>%
historical cbind(raster.df(grid_ca) %>% select(-value)) %>%
slice(index_ca)
<- ggplot(historical) +
g1 geom_point(aes(x=x, y=y, color = maxivt, size = n.seq)) +
geom_sf(data = california, fill = NA, size = 0.25) +
geom_point(aes(x=x, y=y, color = maxivt, size = n.seq), alpha = 0.75) +
scale_color_gradient(
'Average \nPeak 24-hr IVT \n(kg/m/s)',
low = 'white', high = 'darkred') +
scale_size_area(max_size = 4, guide = guide_none()) +
theme(
legend.position = c(0.92, 0.75),
axis.title = element_blank(),
axis.text.x = element_blank(),
# axis.text.x = element_text(angle = 90, vjust = 0.5),
text = element_text(family = 'Segoe UI', size = 8))
<- ggplot(historical) +
g2 geom_point(aes(x=x, y=y, color = duration, size = n.seq)) +
geom_sf(data = california, fill = NA, size = 0.25) +
geom_point(aes(x=x, y=y, color = duration, size = n.seq), alpha = 0.75) +
scale_color_distiller(
'Average \nDuration \n(days)',
palette = 'Blues', direction = 1, breaks = 4*(0:10)) +
scale_size_area(max_size = 4, guide = guide_none()) +
theme(
legend.position = c(0.85, 0.75),
axis.title = element_blank(),
# axis.text.y = element_blank(),
axis.text.x = element_text(angle = 90, vjust = 0.5),
text = element_text(family = 'Segoe UI', size = 8))
<- ggplot(historical) +
g3 geom_point(aes(x=x, y=y, size = n.seq)) +
scale_size_area('Annual Sequence Frequency', max_size = 4, breaks = 1:3) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
legend.position = 'bottom')
+ g2 + plot_spacer() +
g1 plot_layout(nrow = 3, heights = c(10,10,1)) +
inset_element(
get_legend(g3),
left = 0.1, bottom = 0, top = 1, right = 1,
ignore_tag = TRUE, align_to = 'full', on_top = TRUE) +
plot_annotation(tag_levels = 'a', tag_prefix = '(', tag_suffix = ')') &
theme(
plot.tag = element_text(family = 'Segoe UI', face = 'bold', size = 10))
ggsave('_figures/fig3_historic.png', width = 3, height = 6, units = 'in', dpi = 600)
Create dataframes
<-
nseq_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(seq) %>%
group_by(ens, seq.count, decade) %>%
summarize(duration = seq.duration[1], .groups = 'drop') %>%
count(ens, decade) %>%
group_by(decade) %>%
summarize(nn = sum(n)/length(unique(ens))) %>%
right_join(data.frame(decade = 1:7), by = 'decade') %>%
arrange(decade) %>%
mutate(nn = setNA(nn,0)) %>%
pull(nn)
}
}<- nseq_ssp245 / 10
nseq_ssp245 if (progress) cat('\n')
<-
nseq_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(seq) %>%
group_by(ens, seq.count, decade) %>%
summarize(duration = seq.duration[1], .groups = 'drop') %>%
count(ens, decade) %>%
group_by(decade) %>%
summarize(nn = sum(n)/length(unique(ens))) %>%
right_join(data.frame(decade = 1:7), by = 'decade') %>%
arrange(decade) %>%
mutate(nn = setNA(nn,0)) %>%
pull(nn)
}
}<- nseq_ssp585 / 10
nseq_ssp585 if (progress) cat('\n')
<-
nseq_ssp245_eoc foreach (i = 1:ncell(grid_ca), .combine = 'c') %do% {
if (progress) setTxtProgressBar(pb,i)
if (is.null(df_ssp245[[i]])) NA else {
<- df_ssp245[[i]] %>%
temp filter(wateryear(date) %in% 2061:2090) %>%
filter(seq) %>%
group_by(ens, seq.count) %>%
summarize(duration = seq.duration[1], .groups = 'drop')
nrow(temp)/length(unique(temp$ens))
}
}<- nseq_ssp245_eoc / 30
nseq_ssp245_eoc if (progress) cat('\n')
<-
nseq_ssp585_eoc foreach (i = 1:ncell(grid_ca), .combine = 'c') %do% {
if (progress) setTxtProgressBar(pb,i)
if (is.null(df_ssp585[[i]])) NA else {
<- df_ssp585[[i]] %>%
temp filter(wateryear(date) %in% 2061:2090) %>%
filter(seq) %>%
group_by(ens, seq.count) %>%
summarize(duration = seq.duration[1], .groups = 'drop')
nrow(temp)/length(unique(temp$ens))
}
}<- nseq_ssp585_eoc / 30
nseq_ssp585_eoc if (progress) cat('\n')
Change in sequence frequency by decade:
<-
frequency left_join(
%>% as.data.frame %>%
nseq_ssp245 cbind(hist = nseq_ref) %>%
cbind(huc = grid_hucrep[]) %>%
filter(huc > 0) %>%
pivot_longer(cols = starts_with('V'), names_to = 'decade') %>%
rename(ssp245 = value),
%>% as.data.frame %>%
nseq_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 (Medium)',
sim == 'ssp585' ~ 'SSP 5-8.5 (Very High)') %>%
sim %>% factor(levels = levels(.))) %>%
factor mutate(
id = 1:nrow(.),
tag = ifelse(id %in% 1:2, paste0('(', letters[id], ')'), NA))
%>%
frequency group_by(decade_name, sim) %>%
summarize(diff.mean = mean(value-hist), diff.sd = sd(value-hist), .groups = 'drop') %>%
pivot_wider(names_from = sim, values_from = c(diff.mean, diff.sd), names_sep = '.') %>%
gt() %>%
tab_header('Change in Sequence Annual Frequency') %>%
tab_spanner('SSP2-4.5', columns = ends_with('245')) %>%
tab_spanner('SSP5-8.5', columns = ends_with('585')) %>%
fmt_number(starts_with('diff'), decimals = 2) %>%
cols_label(
decade_name = 'Decade',
diff.mean.ssp245 = 'Mean',
diff.sd.ssp245 = 'Std.Dev.',
diff.mean.ssp585 = 'Mean',
diff.sd.ssp585 = 'Std.Dev') %>%
tab_options(
heading.background.color = '#d9d9d9',
column_labels.background.color = '#f2f2f2')
Change in Sequence Annual Frequency | ||||
---|---|---|---|---|
Decade | SSP2-4.5 | SSP5-8.5 | ||
Mean | Std.Dev. | Mean | Std.Dev | |
2020s | 0.13 | 0.11 | 0.25 | 0.20 |
2030s | 0.37 | 0.10 | 0.25 | 0.19 |
2040s | 0.22 | 0.19 | 0.54 | 0.15 |
2050s | 0.51 | 0.15 | 0.82 | 0.27 |
2060s | 0.59 | 0.19 | 0.78 | 0.14 |
2070s | 0.59 | 0.08 | 0.86 | 0.31 |
2080s | 0.72 | 0.16 | 1.13 | 0.17 |
Note: Additional labels were added to this figure in Inkscape.
<-
g.additive ggplot(frequency) +
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.5) +
geom_point(
aes(x = decade_name, y = value-hist, group = huc, color = factor(huc)),
size = 0.75) +
geom_text(
aes(x = 1.5, y = 1.6, 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(
# 'Frequency Change Relative to WY 1981-2010',
'Change in Sequences/Year',
# labels = percent,
expand = expansion(mult = 0.1)) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
panel.grid.major.y = element_line(size = 0.25),
strip.background = element_rect(color = NA, fill = 'grey95'),
strip.text = element_text(size = 8, margin = margin(2,0,2,0)),
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.title = element_text(margin = margin(0,0,0,0)),
legend.text = element_text(margin = margin(0,0,0,0)),
legend.margin = margin(c(0,0,0,0)))
<- gmap + guides(fill = guide_none())
gmap_blank
= 'AAAB\nCCCC'
layout + gmap_blank + get_legend(gmap) + plot_spacer() +
g.additive plot_layout(design = layout, heights = c(10,1))
ggsave('_figures/fig4.png', width = 6, height = 3.5, units = 'in', dpi = 600)
<-
nar_ssp245 foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (!is.null(df_ssp245[[i]])) {
%>%
df_ssp245[[i]] mutate(yr = year(date)) %>%
right_join(decades, by = 'yr') %>%
filter(seq) %>%
group_by(ens, decade, seq.count) %>%
summarize(nar = length(unique(ar.count[!is.na(ar.count)])), .groups = 'drop') %>%
group_by(decade) %>%
summarize(nar = mean(nar)) %>%
right_join(data.frame(decade = 1:7), by = 'decade') %>%
arrange(decade) %>%
pull(nar)
else rep(NA,7)
}
}if (progress) cat('\n')
<-
nar_ssp585 foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (!is.null(df_ssp585[[i]])) {
%>%
df_ssp585[[i]] mutate(yr = year(date)) %>%
right_join(decades, by = 'yr') %>%
filter(seq) %>%
group_by(ens, decade, seq.count) %>%
summarize(nar = length(unique(ar.count[!is.na(ar.count)])), .groups = 'drop') %>%
group_by(decade) %>%
summarize(nar = mean(nar)) %>%
right_join(data.frame(decade = 1:7), by = 'decade') %>%
arrange(decade) %>%
pull(nar)
else rep(NA,7)
}
}if (progress) cat('\n')
<-
nar left_join(
%>% as.data.frame %>%
nar_ssp245 cbind(hist = nar_ref) %>%
cbind(huc = grid_hucrep[]) %>%
filter(huc > 0) %>%
pivot_longer(cols = starts_with('V'), names_to = 'decade') %>%
rename(ssp245 = value),
%>% as.data.frame %>%
nar_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 (Medium)',
sim == 'ssp585' ~ 'SSP 5-8.5 (Very High)') %>%
sim %>% factor(levels = levels(.)))
factor
%>%
nar group_by(decade_name, sim) %>%
summarize(diff.mean = mean(value-hist), diff.sd = sd(value-hist), .groups = 'drop') %>%
pivot_wider(names_from = sim, values_from = c(diff.mean, diff.sd), names_sep = '.') %>%
gt() %>%
tab_header('Change in ARs per Sequence') %>%
tab_spanner('SSP2-4.5', columns = ends_with('245')) %>%
tab_spanner('SSP5-8.5', columns = ends_with('585')) %>%
fmt_number(starts_with('diff'), decimals = 2) %>%
cols_label(
decade_name = 'Decade',
diff.mean.ssp245 = 'Mean',
diff.sd.ssp245 = 'Std.Dev.',
diff.mean.ssp585 = 'Mean',
diff.sd.ssp585 = 'Std.Dev') %>%
tab_options(
heading.background.color = '#d9d9d9',
column_labels.background.color = '#f2f2f2')
Change in ARs per Sequence | ||||
---|---|---|---|---|
Decade | SSP2-4.5 | SSP5-8.5 | ||
Mean | Std.Dev. | Mean | Std.Dev | |
2020s | −0.02 | 0.24 | 0.07 | 0.20 |
2030s | 0.06 | 0.24 | 0.21 | 0.18 |
2040s | 0.00 | 0.11 | 0.25 | 0.24 |
2050s | 0.08 | 0.19 | 0.21 | 0.21 |
2060s | 0.29 | 0.21 | 0.46 | 0.28 |
2070s | 0.24 | 0.29 | 0.44 | 0.22 |
2080s | 0.32 | 0.12 | 0.59 | 0.27 |
<-
pctar_ref foreach (i = 1:ncell(grid_ca), .combine = 'c') %do% {
if (progress) setTxtProgressBar(pb,i)
if (!is.null(df_hist[[i]])) {
%>%
df_hist[[i]] filter(ar) %>%
filter(wateryear(date) %in% 1981:2010) %>%
group_by(ens,ar.count) %>%
summarize(pctar = sum(seq)/length(seq)>0, .groups = 'drop') %>%
pull(pctar) %>% mean
else NA
}
}if (progress) cat('\n')
<-
pctar_ssp245 foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (!is.null(df_ssp245[[i]])) {
%>%
df_ssp245[[i]] mutate(yr = year(date)) %>%
right_join(decades, by = 'yr') %>%
filter(ar) %>%
group_by(ens, decade, ar.count) %>%
summarize(pctar = sum(seq)/length(seq)>0, .groups = 'drop') %>%
group_by(decade) %>%
summarize(pctar = mean(pctar)) %>%
right_join(data.frame(decade = 1:7), by = 'decade') %>%
arrange(decade) %>%
pull(pctar)
else rep(NA,7)
}
}if (progress) cat('\n')
<-
pctar_ssp585 foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (!is.null(df_ssp585[[i]])) {
%>%
df_ssp585[[i]] mutate(yr = year(date)) %>%
right_join(decades, by = 'yr') %>%
filter(ar) %>%
group_by(ens, decade, ar.count) %>%
summarize(pctar = sum(seq)/length(seq)>0, .groups = 'drop') %>%
group_by(decade) %>%
summarize(pctar = mean(pctar)) %>%
right_join(data.frame(decade = 1:7), by = 'decade') %>%
arrange(decade) %>%
pull(pctar)
else rep(NA,7)
}
}if (progress) cat('\n')
<-
pctar left_join(
%>% as.data.frame %>%
pctar_ssp245 cbind(hist = pctar_ref) %>%
cbind(huc = grid_hucrep[]) %>%
filter(huc > 0) %>%
pivot_longer(cols = starts_with('V'), names_to = 'decade') %>%
rename(ssp245 = value),
%>% as.data.frame %>%
pctar_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 (Medium)',
sim == 'ssp585' ~ 'SSP 5-8.5 (Very High)') %>%
sim %>% factor(levels = levels(.)))
factor
%>%
pctar group_by(decade_name, sim) %>%
summarize(diff.mean = mean(value-hist), diff.sd = sd(value-hist), .groups = 'drop') %>%
pivot_wider(names_from = sim, values_from = c(diff.mean, diff.sd), names_sep = '.') %>%
gt() %>%
tab_header('Change in Percentage of ARs Within Sequences') %>%
tab_spanner('SSP2-4.5', columns = ends_with('245')) %>%
tab_spanner('SSP5-8.5', columns = ends_with('585')) %>%
fmt_percent(starts_with('diff'), decimals = 1) %>%
cols_label(
decade_name = 'Decade',
diff.mean.ssp245 = 'Mean',
diff.sd.ssp245 = 'Std.Dev.',
diff.mean.ssp585 = 'Mean',
diff.sd.ssp585 = 'Std.Dev') %>%
tab_options(
heading.background.color = '#d9d9d9',
column_labels.background.color = '#f2f2f2')
Change in Percentage of ARs Within Sequences | ||||
---|---|---|---|---|
Decade | SSP2-4.5 | SSP5-8.5 | ||
Mean | Std.Dev. | Mean | Std.Dev | |
2020s | −4.2% | 7.7% | −1.3% | 3.5% |
2030s | 1.4% | 5.8% | −1.2% | 5.8% |
2040s | −1.9% | 10.3% | 2.5% | 8.8% |
2050s | 1.1% | 7.5% | 5.7% | 5.4% |
2060s | −0.9% | 8.8% | 1.5% | 5.6% |
2070s | 1.1% | 8.3% | 2.2% | 8.4% |
2080s | 3.3% | 8.1% | 8.2% | 4.6% |
Create dataframes
<-
stats_ref foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (!is.null(df_hist[[i]])) {
%>%
df_hist[[i]] filter(wateryear(date) %in% 1981:2010) %>%
filter(seq) %>%
group_by(ens,seq.count) %>%
summarize(maxivt = seq.maxivt[1], duration = seq.duration[1], .groups = 'drop') %>%
summarize(
maxivt.med = median(maxivt),
maxivt.q95 = quantile(maxivt,0.95),
duration.med = median(duration),
duration.q95 = quantile(duration,0.95))
else rep(NA,4)
}
}if (progress) cat('\n')
<-
stats_ssp245_eoc foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (!is.null(df_ssp245[[i]])) {
%>%
df_ssp245[[i]] filter(wateryear(date) %in% 2061:2090) %>%
filter(seq) %>%
group_by(ens,seq.count) %>%
summarize(maxivt = seq.maxivt[1], duration = seq.duration[1], .groups = 'drop') %>%
summarize(
maxivt.med = median(maxivt, na.rm = TRUE),
maxivt.q95 = quantile(maxivt,0.95, na.rm = TRUE),
duration.med = median(duration, na.rm = TRUE),
duration.q95 = quantile(duration,0.95, na.rm = TRUE))
else rep(NA,4)
}
}if (progress) cat('\n')
<-
stats_ssp585_eoc foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (!is.null(df_ssp585[[i]])) {
%>%
df_ssp585[[i]] filter(wateryear(date) %in% 2061:2090) %>%
filter(seq) %>%
group_by(ens,seq.count) %>%
summarize(maxivt = seq.maxivt[1], duration = seq.duration[1], .groups = 'drop') %>%
summarize(
maxivt.med = median(maxivt, na.rm = TRUE),
maxivt.q95 = quantile(maxivt,0.95, na.rm = TRUE),
duration.med = median(duration, na.rm = TRUE),
duration.q95 = quantile(duration,0.95, na.rm = TRUE))
else rep(NA,4)
}
}if (progress) cat('\n')
<-
stats rbind(
which(grid_hucrep[]==1804)]] %>%
df_hist[[filter(wateryear(date) %in% 1981:2010) %>%
filter(seq) %>%
group_by(ens,seq.count) %>%
summarize(maxivt = seq.maxivt[1], duration = seq.duration[1], .groups = 'drop') %>%
mutate(source = 'hist'),
which(grid_hucrep[]==1804)]] %>%
df_ssp245[[filter(wateryear(date) %in% 2061:2090) %>%
filter(seq) %>%
group_by(ens,seq.count) %>%
summarize(maxivt = seq.maxivt[1], duration = seq.duration[1], .groups = 'drop') %>%
mutate(source = 'ssp245'),
which(grid_hucrep[]==1804)]] %>%
df_ssp585[[filter(wateryear(date) %in% 2061:2090) %>%
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
Filter distribution statistics to focus on sequence intensity
<-
intensity cbind(
/ stats_ref) %>%
(stats_ssp245_eoc select(starts_with('maxivt')) %>%
setNames(gsub('maxivt', 'ssp245', names(.))),
/ stats_ref) %>%
(stats_ssp585_eoc select(starts_with('maxivt')) %>%
setNames(gsub('maxivt', 'ssp585', names(.)))) %>%
cbind(raster.df(grid_ca) %>% select(-value)) %>%
slice(index_ca) %>%
pivot_longer(cols = c(-x,-y)) %>%
separate(name, into = c('sim', 'stat')) %>%
mutate(
sim_name = case_when(
== 'ssp245' ~ 'SSP 2-4.5',
sim == 'ssp585' ~ 'SSP 5-8.5') %>%
sim %>% factor(levels = levels(.)[2:1]),
factor stat_name = case_when(
== 'med' ~ 'Median',
stat == 'q95' ~ '95th Percentile') %>%
stat %>% factor(levels = levels(.)[3:1])) factor
Relative intensity change for HUC1804:
left_join(
%>%
stats_ssp245_eoc select(starts_with('maxivt')) %>%
cbind(huc = grid_hucrep[]) %>%
filter(huc == 1804) %>%
select(-huc) %>%
pivot_longer(everything(), names_to = 'stat', values_to = 'ssp245'),
%>%
stats_ref select(starts_with('maxivt')) %>%
cbind(huc = grid_hucrep[]) %>%
filter(huc == 1804) %>%
select(-huc) %>%
pivot_longer(everything(), names_to = 'stat', values_to = 'hist'),
by = 'stat') %>%
mutate(
stat = gsub('maxivt.', '', stat),
stat_name = case_when(
== 'med' ~ 'Median',
stat == 'q95' ~ '95th Percentile') %>%
stat %>% factor(levels = levels(.)[3:1])) %>%
factor mutate(fact = ssp245/hist, pct = fact-1) %>%
select(stat_name, hist, ssp245, fact, pct) %>%
%>%
gt tab_header(
'Relative Change Calculation',
subtitle = 'Sequence Intensity (Peak 24-Hour IVT, kg/m/s)') %>%
fmt_number(c(hist,ssp245), decimals = 0) %>%
fmt_number(fact, decimals = 2) %>%
fmt_percent(pct, decimals = 0) %>%
cols_label(
stat_name = 'Statistic',
hist = 'Historic',
ssp245 = 'SSP2-4.5',
fact = 'Change Factor',
pct = 'Percent Change') %>%
tab_options(
heading.background.color = '#d9d9d9',
column_labels.background.color = '#f2f2f2')
Relative Change Calculation | ||||
---|---|---|---|---|
Sequence Intensity (Peak 24-Hour IVT, kg/m/s) | ||||
Statistic | Historic | SSP2-4.5 | Change Factor | Percent Change |
Median | 539 | 576 | 1.07 | 7% |
95th Percentile | 715 | 832 | 1.16 | 16% |
Statewide average relative intensity change:
%>%
intensity group_by(sim, stat) %>%
summarize(mean = Mean(value)-1, sd = sd(value, na.rm = TRUE), .groups = 'drop') %>%
# mutate(pct.mean = round(100*(fact.mean-1),1), pct.sd = round(100*fact.sd,1)) %>%
mutate(sim_name = ifelse(sim == 'ssp245', 'SSP2-4.5', 'SSP5-8.5')) %>%
select(-sim) %>%
pivot_wider(names_from = stat, values_from = c(mean,sd), names_sep = '.') %>%
select(sim_name, ends_with('med'), ends_with('iqr'), ends_with('q95')) %>%
%>%
gt tab_header(
'Relative Change in Statewide Average Statistics',
subtitle = 'Sequence Intensity (Peak 24-Hour IVT, kg/m/s)') %>%
tab_spanner('Median', columns = ends_with('med')) %>%
tab_spanner('IQR', columns = ends_with('iqr')) %>%
tab_spanner('95th Percentile', columns = ends_with('q95')) %>%
fmt_percent(c(starts_with('mean'), starts_with('sd')), decimals = 1) %>%
cols_label(
sim_name = 'Emission Scenario',
mean.med = 'Mean', sd.med = 'Std.Dev.',
mean.q95 = 'Mean', sd.q95 = 'Std.Dev.') %>%
tab_options(
heading.background.color = '#d9d9d9',
column_labels.background.color = '#f2f2f2')
Relative Change in Statewide Average Statistics | ||||
---|---|---|---|---|
Sequence Intensity (Peak 24-Hour IVT, kg/m/s) | ||||
Emission Scenario | Median | 95th Percentile | ||
Mean | Std.Dev. | Mean | Std.Dev. | |
SSP2-4.5 | 6.5% | 4.3% | 12.7% | 6.2% |
SSP5-8.5 | 11.6% | 5.6% | 19.1% | 5.7% |
Top row (a,b):
<- ggplot(stats) +
a geom_ribbon(
aes(x = maxivt, group = source, ymin = source_name, ymax = source_name, fill = source),
color = NA, alpha = 0) +
# 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.25,
# show.legend = FALSE) +
geom_boxplot(
data = stats %>%
mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>%
group_by(source_name,source,ens) %>%
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),
stat = 'identity', position = 'dodge',
width = 0.8, size = 0.25, alpha = 0.25, show.legend = FALSE) +
geom_point(
data = stats %>%
mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>%
group_by(source_name,source,ens) %>%
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) +
scale_fill_manual(
'Climate Scenario',
labels = c('Historic', 'SSP 2-4.5', 'SSP 5-8.5'),
values = scico(7, palette = 'roma')[3:1],
guide = guide_legend(
override.aes = list(
color = scico(7, palette = 'roma')[3:1], size = 0.25, alpha = 0.5))) +
scale_x_continuous(
'Maximum 24-Hour IVT (kg/m/s)', breaks = 250*(0:10),
limits = c(250,1250), expand = expansion(mult = c(0.005,0.05))) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
panel.grid.major.x = element_line(size = 0.25),
axis.title.y = element_blank(),
legend.position = 'top',
legend.margin = margin(1,1,1,1)) + #c(0.9, 0.75)) +
ggtitle('Sequence Intensity') +
geom_segment(
x = stats %>% filter(source == 'hist') %>% pull(maxivt) %>% median,
xend = stats %>% filter(source == 'hist') %>% pull(maxivt) %>% median,
y = 2.6, yend = 3.4,
color = scico(11, palette = 'roma')[4]) +
geom_segment(
x = stats %>% filter(source == 'ssp245') %>% pull(maxivt) %>% median,
xend = stats %>% filter(source == 'ssp245') %>% pull(maxivt) %>% median,
y = 1.6, yend = 2.4,
color = scico(11, palette = 'roma')[3]) +
geom_segment(
x = stats %>% filter(source == 'ssp585') %>% pull(maxivt) %>% median,
xend = stats %>% filter(source == 'ssp585') %>% pull(maxivt) %>% median,
y = 0.6, yend = 1.4,
color = scico(11, palette = 'roma')[1])
a
<- ggplot() +
b geom_tile(
data = raster.df(grid_hucrep),
aes(x=x, y=y, fill = value==1804)) +
scale_fill_manual(
values = c('white', 'grey50'), guide = guide_none()) +
geom_tile(
data = raster.df(grid_hucrep*grid_ca/grid_ca) %>% arrange(value==1804),
aes(x=x, y=y, color = value==1804), fill = NA) +
scale_color_manual(
values = c('grey70', 'black'), na.value = NA, guide = guide_none()) +
geom_sf(data = cal, fill = NA, color = 'black', size = 0.25) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
axis.text.x = element_text(angle = 90, vjust = 0.5),
axis.title = element_blank(),
legend.background = element_rect(fill = NA))
b
Bottom row (c,d,e):
<- ggplot(stats %>% filter(source != 'ssp585')) +
c ## medians
geom_vline(
xintercept = stats %>% filter(source == 'hist') %>% pull(maxivt) %>% median,
size = 0.75, color = baker[3]) +
geom_vline(
xintercept = stats %>% filter(source == 'ssp245') %>% pull(maxivt) %>% median,
size = 0.75, color = baker[3]) +
## boxplot
# stat_boxplot_custom(
# aes(x = maxivt, y = source_name, fill = source),
# qs = c(.05, .25, .5, .75, .95),
# outlier.size = 0.5, size = 0.25, outlier.alpha = 0, alpha = 0.25,
# width = 0.5,
# show.legend = FALSE) +
geom_boxplot(
data = stats %>% filter(source != 'ssp585') %>%
mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>%
group_by(source_name,source) %>%
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), fill = source),
stat = 'identity', width = 0.5, size = 0.2, alpha = 0.25, show.legend = FALSE) +
scale_fill_manual(values = scico(7, palette = 'roma')[c(3:1)]) +
## 95th quantiles
geom_vline(
xintercept = stats %>% filter(source == 'hist') %>% pull(maxivt) %>% quantile(0.95),
size = 0.75, color = roma.colors[5], linetype = 'dashed') +
geom_vline(
xintercept = stats %>% filter(source == 'ssp245') %>% pull(maxivt) %>% quantile(0.95),
size = 0.75, color = roma.colors[5], linetype = 'dashed') +
## plot specs
scale_x_continuous(
'Maximum 24-Hour IVT (kg/m/s)', breaks = 250*(0:10),
limits = c(250,1250), expand = expansion(mult = c(0.005,0.05))) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
panel.grid.major.x = element_line(size = 0.25),
axis.title.x = element_text(hjust = 0.25),
axis.title.y = element_blank())
c
<- intensity %>%
de mutate(
id = 1:nrow(.),
tag = ifelse(id %in% 1:3, paste0('(', letters[4:6], ')'), NA),
sim_name = factor(sim_name, levels = rev(levels(sim_name)))) %>%
select(-id) %>%
filter(stat != 'iqr') %>%
ggplot() +
geom_vline(xintercept = 0, color = 'grey50', size = 0.35) +
geom_density(
aes(x = value-1, group = sim_name, fill = sim_name, color = sim_name),
alpha = 0.35, size = 0.35) +
scale_color_manual(
'Climate Scenario', values = scico(7, palette = 'roma')[2:1]) +
scale_fill_manual(
'Climate Scenario', values = scico(7, palette = 'roma')[2:1]) +
guides(color = guide_none(), fill = guide_none()) +
facet_wrap(~stat_name, nrow = 1) +
scale_x_continuous(
'Percent Change Relative to 1981-2010', labels = percent) +
scale_y_origin('Probability Density') +
theme(
text = element_text(family = 'Segoe UI', size = 8),
legend.margin = margin(0,0,0,0),
legend.position = 'top', #c(0.85,0.5),
panel.grid.major.x = element_line(size = 0.25),
strip.placement = 'outside',
strip.background = element_rect(fill = 'grey95', color = NA),
strip.text = element_text(size = 8, margin = margin(2,2,2,2)))
de
ggsave(a, filename = '_figures/fig5a.png', width = 4.25, height = 2.25, dpi = 600)
ggsave(b, filename = '_figures/fig5b.png', width = 1.75, height = 1.75, dpi = 600)
ggsave(c, filename = '_figures/fig5c.png', width = 4.25, height = 1.5, dpi = 600)
ggsave(de, filename = '_figures/fig5de.png', width = 3, height = 1.75, dpi = 600)
The individual figure panels were annotated and combined in Inkscape.
Filter distribution statistics to focus on sequence duration
<-
duration cbind(
/ stats_ref) %>%
(stats_ssp245_eoc select(starts_with('duration')) %>%
setNames(gsub('duration', 'ssp245', names(.))),
/ stats_ref) %>%
(stats_ssp585_eoc select(starts_with('duration')) %>%
setNames(gsub('duration', 'ssp585', names(.)))) %>%
cbind(raster.df(grid_ca) %>% select(-value)) %>%
slice(index_ca) %>%
pivot_longer(cols = c(-x,-y)) %>%
separate(name, into = c('sim', 'stat')) %>%
mutate(
sim_name = case_when(
== 'ssp245' ~ 'SSP 2-4.5',
sim == 'ssp585' ~ 'SSP 5-8.5') %>%
sim %>% factor(levels = levels(.)[2:1]),
factor stat_name = case_when(
== 'med' ~ 'Median',
stat == 'iqr' ~ 'Interquartile Range (IQR)',
stat == 'q95' ~ '95th Percentile') %>%
stat %>% factor(levels = levels(.)[3:1])) factor
Relative duration change for HUC1804:
left_join(
%>%
stats_ssp245_eoc select(starts_with('duration')) %>%
cbind(huc = grid_hucrep[]) %>%
filter(huc == 1804) %>%
select(-huc) %>%
pivot_longer(everything(), names_to = 'stat', values_to = 'ssp245'),
%>%
stats_ref select(starts_with('duration')) %>%
cbind(huc = grid_hucrep[]) %>%
filter(huc == 1804) %>%
select(-huc) %>%
pivot_longer(everything(), names_to = 'stat', values_to = 'hist'),
by = 'stat') %>%
mutate(
stat = gsub('duration.', '', stat),
stat_name = case_when(
== 'med' ~ 'Median',
stat == 'q95' ~ '95th Percentile') %>%
stat %>% factor(levels = levels(.)[3:1])) %>%
factor mutate(fact = ssp245/hist, pct = fact-1) %>%
select(stat_name, hist, ssp245, fact, pct) %>%
%>%
gt tab_header(
'Relative Change Calculation',
subtitle = 'Sequence Duration (days)') %>%
fmt_number(c(hist,ssp245), decimals = 0) %>%
fmt_number(fact, decimals = 2) %>%
fmt_percent(pct, decimals = 0) %>%
cols_label(
stat_name = 'Statistic',
hist = 'Historic',
ssp245 = 'SSP2-4.5',
fact = 'Change Factor',
pct = 'Percent Change') %>%
tab_options(
heading.background.color = '#d9d9d9',
column_labels.background.color = '#f2f2f2')
Relative Change Calculation | ||||
---|---|---|---|---|
Sequence Duration (days) | ||||
Statistic | Historic | SSP2-4.5 | Change Factor | Percent Change |
Median | 20 | 23 | 1.15 | 15% |
95th Percentile | 52 | 53 | 1.02 | 2% |
Statewide average relative duration change:
%>%
duration group_by(sim, stat) %>%
summarize(mean = Mean(value)-1, sd = sd(value, na.rm = TRUE), .groups = 'drop') %>%
mutate(sim_name = ifelse(sim == 'ssp245', 'SSP2-4.5', 'SSP5-8.5')) %>%
pivot_wider(names_from = stat, values_from = c(mean,sd), names_sep = '.') %>%
select(sim_name, ends_with('med'), ends_with('iqr'), ends_with('q95')) %>%
%>%
gt tab_header(
'Relative Change in Statewide Average Statistics',
subtitle = 'Sequence Duration (days)') %>%
tab_spanner('Median', columns = ends_with('med')) %>%
tab_spanner('IQR', columns = ends_with('iqr')) %>%
tab_spanner('95th Percentile', columns = ends_with('q95')) %>%
fmt_percent(c(starts_with('mean'), starts_with('sd')), decimals = 1) %>%
cols_label(
sim_name = 'Emission Scenario',
mean.med = 'Mean', sd.med = 'Std.Dev.',
mean.q95 = 'Mean', sd.q95 = 'Std.Dev.') %>%
tab_options(
heading.background.color = '#d9d9d9',
column_labels.background.color = '#f2f2f2')
Relative Change in Statewide Average Statistics | ||||
---|---|---|---|---|
Sequence Duration (days) | ||||
Emission Scenario | Median | 95th Percentile | ||
Mean | Std.Dev. | Mean | Std.Dev. | |
SSP2-4.5 | 10.9% | 10.5% | 7.6% | 12.4% |
SSP5-8.5 | 17.9% | 10.0% | 28.4% | 16.8% |
<- ggplot(stats) +
a geom_ribbon(
aes(x = duration, group = source, ymin = source_name, ymax = source_name, fill = source),
color = NA, alpha = 0) +
# 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.35,
# show.legend = FALSE) +
geom_boxplot(
data = stats %>%
mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>%
group_by(source_name,source,ens) %>%
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),
stat = 'identity', position = 'dodge',
width = 0.8, size = 0.25, alpha = 0.35, show.legend = FALSE) +
geom_point(
data = stats %>%
mutate(source_name = factor(source_name, levels = c('SSP 5-8.5','SSP 2-4.5','Historic'))) %>%
group_by(source_name,source,ens) %>%
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) +
scale_fill_manual(
'Climate Scenario',
labels = c('Historic', 'SSP 2-4.5', 'SSP 5-8.5'),
values = scico(7, palette = 'roma')[5:7],
guide = guide_legend(
override.aes = list(
color = scico(7, palette = 'roma')[5:7], size = 0.25, alpha = 0.5))) +
scale_x_continuous('Duration (days)', breaks = 30*(0:10)) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
panel.grid.major.x = element_line(size = 0.25),
axis.title.y = element_blank(),
legend.position = 'top',
legend.margin = margin(1,1,1,1)) +
ggtitle('Sequence Duration') +
geom_segment(
x = stats %>% filter(source == 'hist') %>% pull(duration) %>% median,
xend = stats %>% filter(source == 'hist') %>% pull(duration) %>% median,
y = 2.6, yend = 3.4,
color = scico(11, palette = 'roma')[8]) +
geom_segment(
x = stats %>% filter(source == 'ssp245') %>% pull(duration) %>% median,
xend = stats %>% filter(source == 'ssp245') %>% pull(duration) %>% median,
y = 1.6, yend = 2.4,
color = scico(11, palette = 'roma')[9]) +
geom_segment(
x = stats %>% filter(source == 'ssp585') %>% pull(duration) %>% median,
xend = stats %>% filter(source == 'ssp585') %>% pull(duration) %>% median,
y = 0.6, yend = 1.4,
color = scico(11, palette = 'roma')[11])
a
<- duration %>%
bc mutate(
id = 1:nrow(.),
tag = ifelse(id %in% 1:3, paste0('(', letters[2:4], ')'), NA),
sim_name = factor(sim_name, levels = rev(levels(sim_name)))) %>%
select(-id) %>%
filter(stat != 'iqr') %>%
ggplot() +
geom_vline(xintercept = 0, color = 'grey50', size = 0.35) +
geom_density(
aes(x = value-1, group = sim_name, fill = sim_name, color = sim_name),
alpha = 0.35, size = 0.35) +
scale_color_manual('Climate\nScenario', values = scico(7, palette = 'roma')[6:7]) +
scale_fill_manual('Climate\nScenario', values = scico(7, palette = 'roma')[6:7]) +
guides(color = guide_none(), fill = guide_none()) +
facet_wrap(~stat_name, nrow = 1) +
scale_x_continuous(
'Percent Change Relative to 1981-2010', labels = percent) +
scale_y_origin('Probability Density\n') +
coord_cartesian(xlim = c(NA,1.25)) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
legend.margin = margin(0,0,0,0),
legend.position = c(0.9,0.5),
panel.grid.major.x = element_line(size = 0.25),
strip.placement = 'outside',
strip.background = element_rect(fill = 'grey95', color = NA),
strip.text = element_text(margin = margin(2,2,2,2)),
plot.margin = margin(5,5,5,5))
bc
ggsave(a, filename = '_figures/fig6a.png', width = 4.25, height = 2.25, dpi = 600)
ggsave(bc, filename = '_figures/fig6bc.png', width = 4.18, height = 1.75, dpi = 600)
The individual figure panels were annotated and combined in Inkscape.
Create dataframes
<-
extremes_ref foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (is.null(df_hist[[i]])) rep(NA,5) else {
<- df_hist[[i]] %>%
temp filter(seq) %>%
filter(wateryear(date) %in% 1981:2010) %>%
group_by(ens,seq.count) %>%
summarize(duration = seq.duration[1], .groups = 'drop')
%>%
temp group_by(ens) %>%
summarize(n = sum(duration>60)/30) %>%
right_join(data.frame(ens = unique(df_hist[[i]]$ens))) %>%
mutate(n = setNA(n,0)) %>%
right_join(data.frame(ens = 1:5), by = 'ens') %>%
arrange(ens) %>%
pull(n)
}%>% apply(1, Max)
} if (progress) cat('\n')
<-
extremes_ssp245_eoc foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (is.null(df_ssp245[[i]])) rep(NA,5) else {
<- df_ssp245[[i]] %>%
temp filter(seq) %>%
filter(wateryear(date) %in% 2061:2090) %>%
group_by(ens,seq.count) %>%
summarize(duration = seq.duration[1], .groups = 'drop')
%>%
temp group_by(ens) %>%
summarize(n = sum(duration>60)/30) %>%
right_join(data.frame(ens = unique(df_ssp245[[i]]$ens))) %>%
mutate(n = setNA(n,0)) %>%
right_join(data.frame(ens = 1:5), by = 'ens') %>%
arrange(ens) %>%
pull(n)
}%>% apply(1, Max)
} if (progress) cat('\n')
<-
extremes_ssp585_eoc foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (is.null(df_ssp585[[i]])) rep(NA,5) else {
<- df_ssp585[[i]] %>%
temp filter(seq) %>%
filter(wateryear(date) %in% 2061:2090) %>%
group_by(ens,seq.count) %>%
summarize(duration = seq.duration[1], .groups = 'drop')
%>%
temp group_by(ens) %>%
summarize(n = sum(duration>60)/30) %>%
right_join(data.frame(ens = unique(df_ssp585[[i]]$ens))) %>%
mutate(n = setNA(n,0)) %>%
right_join(data.frame(ens = 1:5), by = 'ens') %>%
arrange(ens) %>%
pull(n)
}%>% apply(1, Max)
} if (progress) cat('\n')
Historical vs. projected future spatial coverage of
super-sequences:
rbind(
c('Historic', 'single', Sum(extremes_ref>0) / sum(!is.na(extremes_ref))),
c('SSP2-4.5', 'single', Sum(extremes_ssp245_eoc>0) / sum(!is.na(extremes_ssp245_eoc))),
c('SSP5-8.5', 'single', Sum(extremes_ssp585_eoc>0) / sum(!is.na(extremes_ssp585_eoc))),
c('Historic', 'multiple', Sum(extremes_ref>(1/30)) / sum(!is.na(extremes_ref))),
c('SSP2-4.5', 'multiple', Sum(extremes_ssp245_eoc>(1/30)) / sum(!is.na(extremes_ssp245_eoc))),
c('SSP5-8.5', 'multiple', Sum(extremes_ssp585_eoc>(1/30)) / sum(!is.na(extremes_ssp585_eoc)))) %>%
%>% setNames(c('source', 'super', 'pct')) %>%
as.data.frame mutate(pct = toNumber(pct)) %>%
pivot_wider(names_from = super, values_from = pct) %>%
%>%
gt tab_header(
'Super-Sequence Coverage',
subtitle = 'Percent of Grid Cells Affected over 30 Years (1981-2010 or 2021-2050)') %>%
fmt_percent(c(single, multiple), decimals = 1) %>%
cols_label(
source = 'Emission Scenario',
single = 'Cells with >0 Super-Sequences',
multiple = 'Cells with Multiple Super-Sequences') %>%
tab_options(
heading.background.color = '#d9d9d9',
column_labels.background.color = '#f2f2f2')
Super-Sequence Coverage | ||
---|---|---|
Percent of Grid Cells Affected over 30 Years (1981-2010 or 2021-2050) | ||
Emission Scenario | Cells with >0 Super-Sequences | Cells with Multiple Super-Sequences |
Historic | 73.3% | 39.4% |
SSP2-4.5 | 86.7% | 67.3% |
SSP5-8.5 | 98.8% | 93.9% |
<-
bucket_ref foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (!is.null(df_hist[[i]])) {
<- df_hist[[i]] %>%
temp filter(seq) %>%
filter(wateryear(date) %in% 1981:2010) %>%
group_by(ens,seq.count) %>%
summarize(duration = seq.duration[1], .groups = 'drop') %>%
mutate(duration_cut = cut(duration, breaks = c(0,14,30,60,200)))
summary(temp$duration_cut) / length(unique(temp$ens)) / 30) %>%
(%>% t %>% as.data.frame
as.matrix else rep(NA,6)
}
}if (progress) cat('\n')
<-
bucket_ssp245_eoc foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (!is.null(df_ssp245[[i]])) {
<- df_ssp245[[i]] %>%
temp filter(seq) %>%
filter(wateryear(date) %in% 2061:2090) %>%
group_by(ens,seq.count) %>%
summarize(duration = seq.duration[1], .groups = 'drop') %>%
mutate(duration_cut = cut(duration, breaks = c(0,14,30,60,200)))
summary(temp$duration_cut) / length(unique(temp$ens)) / 30) %>%
(%>% t %>% as.data.frame
as.matrix else rep(NA,6)
}
}if (progress) cat('\n')
<-
bucket_ssp585_eoc foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (!is.null(df_ssp585[[i]])) {
<- df_ssp585[[i]] %>%
temp filter(seq) %>%
filter(wateryear(date) %in% 2061:2090) %>%
group_by(ens,seq.count) %>%
summarize(duration = seq.duration[1], .groups = 'drop') %>%
mutate(duration_cut = cut(duration, breaks = c(0,14,30,60,200)))
summary(temp$duration_cut) / length(unique(temp$ens)) / 30) %>%
(%>% t %>% as.data.frame
as.matrix else rep(NA,6)
}
}if (progress) cat('\n')
<-
bucket rbind(
%>%
bucket_ref cbind(raster.df(grid_huc4) %>% transmute(huc4 = value)) %>%
slice(index_ca) %>%
mutate(id = index_ca, source = 'hist'),
%>%
bucket_ssp245_eoc cbind(raster.df(grid_huc4) %>% transmute(huc4 = value)) %>%
slice(index_ca) %>%
mutate(id = index_ca, source = 'ssp245'),
%>%
bucket_ssp585_eoc cbind(raster.df(grid_huc4) %>% transmute(huc4 = value)) %>%
slice(index_ca) %>%
mutate(id = index_ca, source = 'ssp585')) %>%
pivot_longer(cols = c(-id,-source,-huc4), names_to = 'bucket', values_to = 'freq') %>%
mutate(
bucket_label = case_when(
== '(0,14]' ~ 'Duration \u2264 14 Days',
bucket == '(14,30]' ~ 'Duration = 15-30 Days',
bucket == '(30,60]' ~ 'Duration = 31-60 Days',
bucket == '(60,200]' ~ 'Duration > 60 Days') %>%
bucket factor(levels = c(
'Duration \u2264 14 Days', 'Duration = 15-30 Days',
'Duration = 31-60 Days', 'Duration > 60 Days'))) %>%
pivot_wider(names_from = source, values_from = freq)
<- bucket %>%
temp filter(bucket_label == levels(bucket_label)[4]) %>%
select(hist, ssp245, ssp585) %>%
apply(2, Max)
<- (temp/temp[1])[-1]
ratio
<- bucket %>%
uppers group_by(bucket) %>%
summarize(across(c('hist', 'ssp245', 'ssp585'), Max)) %>%
select(-bucket)
<- list()
g.high for (i in 1:4) {
<-
temp ggplot(bucket %>% filter(bucket_label == levels(bucket_label)[i])) +
geom_point(aes(x = hist, y = ssp245), size = 0.75)
<- ggplot_build(temp)$layout$panel_params[[1]]$x$breaks %>% diff %>% .[1]
dx <-
g.high[[i]] ggplot(bucket %>% filter(bucket_label == levels(bucket_label)[i])) +
geom_point(aes(x = hist, y = ssp245), size = 0.75) +
facet_wrap(~bucket_label) +
scale_x_origin('Historic (1921-2010) Annual Frequency') +
scale_y_origin(
'SSP2-4.5 Annual\nFrequency (WY 2061-2090)',
breaks = seq(0, uppers$ssp245[i]+2*dx, 2*dx)) +
geom_parity() +
coord_cartesian(
clip = 'off', xlim = c(0, uppers$hist[i]), ylim = c(0, ratio[1]*1.1*uppers$hist[i])) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
panel.grid.major = element_line(size = 0.25),
axis.title = element_blank(),
strip.background = element_rect(color = NA, fill = 'grey95'),
strip.text = element_text(margin = margin(2,2,2,2)))
}1]] <- g.high[[1]] + theme(axis.title.y = element_text())
g.high[[
<- list()
g.low for (i in 1:4) {
<-
temp ggplot(bucket %>% filter(bucket_label == levels(bucket_label)[i])) +
geom_point(aes(x = hist, y = ssp585), size = 0.75)
<- ggplot_build(temp)$layout$panel_params[[1]]$x$breaks %>% diff %>% .[1]
dx <-
g.low[[i]] ggplot(bucket %>% filter(bucket_label == levels(bucket_label)[i])) +
geom_point(aes(x = hist, y = ssp585), size = 0.75) +
facet_wrap(~bucket_label) +
scale_x_origin('Historic (1921-2010)\nAnnual Frequency') +
scale_y_origin(
'SSP5-8.5 Annual\nFrequency (WY 2061-2090)',
breaks = seq(0, uppers$ssp585[i]+2*dx, 2*dx)) +
geom_parity() +
coord_cartesian(
clip = 'off', xlim = c(0, uppers$hist[i]), ylim = c(0, ratio[2]*uppers$hist[i])) +
theme(
text = element_text(family = 'Segoe UI', size = 8),
panel.grid.major = element_line(size = 0.25),
axis.title = element_blank(),
strip.background = element_rect(color = NA, fill = 'grey95'),
strip.text = element_text(margin = margin(2,2,2,2)))
}1]] <- g.low[[1]] + theme(axis.title.y = element_text())
g.low[[
<- "ABCD\nEFGH"
layout 1]] + g.high[[2]] + g.high[[3]] + g.high[[4]] + #g.high[[5]] +
g.high[[1]] + g.low[[2]] + g.low[[3]] + g.low[[4]] + #g.low[[5]] +
g.low[[plot_layout(
design = layout,
widths = rep(1,5), heights = c(3,5)) +
plot_annotation(
caption = 'Historic Annual Frequency (WY 1981-2010)',
tag_levels = list(letters[c(1,3,5,7,2,4,6,8)]), tag_prefix = '(', tag_suffix = ')', ) &
theme(
plot.tag = element_text(
family = 'Segoe UI', face = 'bold', size = 10, margin = margin(-5,-5,-5,-5)),
plot.caption = element_text(size = 8, hjust = 0.5, margin = margin(0,0,0,0)))
ggsave('_figures/fig7_returnperiod.png', width = 6, height = 3.5, units = 'in', dpi = 600)
Create dataframes
<-
exceed_ref foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (!is.null(df_hist[[i]])) {
%>%
df_hist[[i]] mutate(wy = wateryear(date)) %>%
filter(wy %in% 1981:2010) %>%
filter(seq) %>%
group_by(ens, wy, seq.count) %>%
summarize(duration = seq.duration[1], .groups = 'drop_last') %>%
summarize(exceed60 = (sum(duration>60)/length(duration)) > 0, .) %>%
right_join(expand.grid(ens = unique(df_hist[[i]]$ens), wy = 1981:2010), by = c('ens', 'wy')) %>%
mutate(exceed60 = setNA(exceed60, FALSE)) %>%
summarize(across(starts_with('exceed'), mean)) %>%
right_join(data.frame(ens = 1:5), by = 'ens') %>%
arrange(ens) %>% pull(exceed60)
else rep(NA,5)
}
}if (progress) cat('\n')
<-
exceed_ssp245_eoc foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (!is.null(df_ssp245[[i]])) {
%>%
df_ssp245[[i]] mutate(wy = wateryear(date)) %>%
filter(wy %in% 2061:2090) %>%
filter(seq) %>%
group_by(ens, wy, seq.count) %>%
summarize(duration = seq.duration[1], .groups = 'drop_last') %>%
summarize(exceed60 = (sum(duration>60)/length(duration)) > 0) %>%
right_join(expand.grid(ens = unique(df_ssp245[[i]]$ens), wy = 2061:2090), by = c('ens', 'wy')) %>%
mutate(exceed60 = setNA(exceed60, FALSE)) %>%
summarize(across(starts_with('exceed'), mean)) %>%
right_join(data.frame(ens = 1:5), by = 'ens') %>%
arrange(ens) %>% pull(exceed60)
else rep(NA,4)
}
}if (progress) cat('\n')
<-
exceed_ssp585_eoc foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
if (progress) setTxtProgressBar(pb,i)
if (!is.null(df_ssp585[[i]])) {
%>%
df_ssp585[[i]] mutate(wy = wateryear(date)) %>%
filter(wy %in% 2061:2090) %>%
filter(seq) %>%
group_by(ens, wy, seq.count) %>%
summarize(duration = seq.duration[1], .groups = 'drop_last') %>%
summarize(exceed60 = (sum(duration>60)/length(duration)) > 0) %>%
right_join(expand.grid(ens = unique(df_ssp585[[i]]$ens), wy = 2061:2090), by = c('ens', 'wy')) %>%
mutate(exceed60 = setNA(exceed60, FALSE)) %>%
summarize(across(starts_with('exceed'), mean)) %>%
right_join(data.frame(ens = 1:5), by = 'ens') %>%
arrange(ens) %>% pull(exceed60)
else rep(NA,4)
}
}if (progress) cat('\n')
Percentage of years with super-sequences:
data.frame(
`Historic` = exceed_ref %>% apply(1, Mean),
`SSP2-4.5` = exceed_ssp245_eoc %>% apply(1, Mean),
`SSP5-8.5` = exceed_ssp585_eoc %>% apply(1, Mean)) %>%
slice(index_ca) %>%
pivot_longer(everything()) %>%
group_by(name) %>%
# as_tibble() %>%
summarize(across(value, .fns = list(mean = mean, sd = sd, rp = function(x) 1/mean(x)))) %>%
# pivot_longer(-name, names_to = 'stat')
%>%
gt tab_header('Percent of Years with Super-Sequences') %>%
fmt_markdown(name) %>%
fmt_percent(c(value_mean, value_sd), decimals = 1) %>%
fmt_number(value_rp, decimals = 1) %>%
cols_label(
name = 'Emission Scenario',
value_mean = 'Mean (%)',
value_sd = 'Std. Dev. (%)',
value_rp = 'Return Period (yrs)') %>%
tab_options(
heading.background.color = '#d9d9d9',
column_labels.background.color = '#f2f2f2')
Percent of Years with Super-Sequences | |||
---|---|---|---|
Emission Scenario | Mean (%) | Std. Dev. (%) | Return Period (yrs) |
Historic |
1.9% | 1.9% | 53.0 |
SSP2.4.5 |
4.2% | 3.6% | 23.9 |
SSP5.8.5 |
10.5% | 6.3% | 9.5 |