This markdown script reproduces Figure 2 from 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).

Setup

source('_data/setup.R')

## load required packages
require(prism); prism_set_dl_dir('D:/Research/_data/PRISM/files')
require(patchwork)
# require(mapview)
# require(glue)
## 24-hour dataset
load('_scripts/_checkpoints/df_24hr_0209.Rdata')

## date metadata 
ts_merra <- seq(ymd_hms('1980-1-1 00:00:00'), ymd_hms('2021-12-31 21:00:00'), by = '3 hours')
dates_merra <- unique(as.Date(ts_merra))

## progress bar
progress <- FALSE

## geospatial 
cal <- st_union(california)

Precipitation

Create dataframe

Load file names

# ## note: do not need to run this again (takes a long time)
# get_precip_dailys(
#   type = 'ppt', 
#   minDate = '1981-01-01', 
#   maxDate = '2021-12-31', 
#   keepZip = FALSE)
# ## load file names
# ppt <-
#   prism_archive_subset(
#     type = 'ppt', temp_period = 'daily',
#     dates = dates_merra)
# 
# ## filter to wet season
# dates_precip <- 
#   lapply(ppt, function(x) x %>% str_split('_') %>% .[[1]] %>% .[5] %>% ymd) %>% 
#   reduce(c)
# ppt <- ppt[month(dates_precip) %in% c(10:12,1:3)]
# dates_precip <- dates_precip[month(dates_precip) %in% c(10:12,1:3)]

Subset to California cells

# ## determine which cells I actually need
# precip_raster <- ppt[1] %>% pd_to_file %>% rast %>% terra::crop(extent(grid_ca)) %>% terra::mask(california)
# precip_match <- grid_ca %>% rast %>% 
#   project(precip_raster, method = 'near') %>%
#   .[] %>% c %>% 
#   cbind(precip_raster %>% terra::xyFromCell(1:ncell(.))) %>% 
#   data.frame %>% setNames(c('grid_merra', 'x', 'y')) %>% 
#   mutate(
#     value_ppt = c(precip_raster[]),
#     grid_ppt = 1:nrow(.), 
#     grid_merra = ifelse(is.nan(grid_merra), NA, grid_merra)) %>% 
#   mutate(
#     grid_ppt = ifelse(is.na(value_ppt), NA, grid_ppt),
#     valid = !is.na(grid_merra) & !is.na(grid_ppt)) %>% 
#   select(-value_ppt)
# # ggplot(precip_match) + geom_raster(aes(x=x, y=y, fill=valid)) + geom_sf(data = cal, fill = NA)
# precip_valid <- precip_match$valid

Extract precipitation values

# # takes about five minutes to run locally (parallelizing takes longer)
# if (progress) pb <- txtProgressBar(min = 0, max = length(ppt), style = 3)
# df.precip <- 
#   foreach (d = 1:length(ppt), .combine = 'cbind') %do% {
#     if (progress) setTxtProgressBar(pb,d)
#     ppt[d] %>% 
#       pd_to_file %>% rast %>% terra::crop(extent(grid_ca)) %>% 
#       .[] %>% c %>% .[precip_valid]
#   } %>% t
# 
# ## gut check
# dim(df.precip)  
# c(length(dates_precip), sum(precip_valid))

Find AR days and sequence days

# if (progress) pb <- txtProgressBar(min = 0, max = ncell(grid_ca), style = 3)
# ar_binary <-
#   foreach (i = 1:ncell(grid_ca), .combine = 'cbind') %do% {
#     if (progress) setTxtProgressBar(pb,i)
#     if (i %in% index_ca) {
#       df_24hr[[i]] %>% 
#         filter(date %in% dates_precip) %>% 
#         select(ar)
#     } else rep(NA,length(dates_precip))
#   }  
# seq_binary <-
#   foreach (i = 1:ncell(grid_ca), .combine = 'cbind') %do% {
#     if (progress) setTxtProgressBar(pb,i)
#     if (i %in% index_ca) {
#       df_24hr[[i]] %>% 
#         filter(date %in% dates_precip) %>% 
#         select(seq)
#     } else rep(NA,length(dates_precip))
#   }  

Create precip.extreme dataframe

# if (progress) pb <- txtProgressBar(min = 0, max = length(precip_valid), style = 3)
# cl <- makeCluster(cores)
# registerDoSNOW(cl)
# precip.extreme <-
#   foreach (
#     j = 1:length(precip_valid),
#     .options.snow = if (progress) opts,
#     .packages = c('tidyverse'),
#     .combine = 'rbind') %dopar% {
#       if (precip_valid[j]) {
#         temp <-
#           data.frame(
#             precip = unname(df.precip[,match(j,precip_match$grid_ppt[precip_valid])]),
#             ar = ar_binary[,precip_match$grid_merra[j]],
#             seq = seq_binary[,precip_match$grid_merra[j]]) %>%
#           count(ar, seq, precip.extreme = precip > quantile(precip, 0.95, na.rm = TRUE)) %>%
#           group_by(ar, seq) %>%
#           mutate(p.all = sum(n)) %>%
#           ungroup %>%
#           mutate(p.all = p.all/sum(p.all)*2) %>%
#           filter(precip.extreme) %>%
#           mutate(p.precip = n/sum(n))
#         c(temp %>% filter(ar) %>% summarize(across(starts_with('p.'), sum)) %>% unlist,
#           temp %>% filter(seq) %>% summarize(across(starts_with('p.'), sum)) %>% unlist) %>%
#           setNames(c('ar.all', 'ar.precip', 'seq.all', 'seq.precip'))
#       } else rep(NA,4)
#     }
# stopCluster(cl)

# ## save out
# save(
#   dates_precip, precip_match, df.precip, precip.extreme,
#   file = '_data/ARTMIP/ppt_0905.Rdata')

Load from checkpoint

## load from checkpoint
load('_data/ARTMIP/ppt_0905.Rdata')

Plot results

g.precip <- precip_match %>% 
  cbind(precip.extreme) %>% filter(valid) %>% 
  select(x, y, starts_with('ar'), starts_with('seq')) %>%
  pivot_longer(c(-x,-y)) %>% 
  separate(name, into = c('arseq', 'extreme'), sep = '\\.') %>% 
  filter(extreme == 'precip') %>%
  ggplot() + 
  geom_raster(aes(x=x, y=y, fill = value)) + 
  geom_sf(data = california, fill = NA, size = 0.2) + 
  facet_grid(arseq ~ extreme) + 
  scale_fill_scico(
    'Percentage of Days     ', palette = 'lapaz', direction = -1, 
    limits = c(0,1), labels = percent) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    axis.title = element_blank(),
    axis.text.x = element_text(angle = 90, vjust = 0.5),
    axis.text.y = element_blank(), axis.line.y = element_blank(), axis.ticks.y = element_blank(),
    strip.background.y = element_blank(), strip.text.y = element_blank(),
    plot.background = element_rect(fill = 'white', color = NA),
    strip.background.x = element_blank(), strip.text.x = element_blank(),
    plot.margin = margin(5,2,5,2),
    legend.position = 'bottom',
    legend.key.width = unit(0.5, 'in'), 
    legend.margin = margin(0,0,0,0))
# g.precip

Runoff

Load data

load('_data/streamflow/_checkpoints/huc_runoff_0928.Rdata')
df.runoff <- huc_runoff %>% 
  select(date, huc8, runoff) %>% 
  pivot_wider(id_cols = 'date', names_from = 'huc8', values_from = 'runoff') %>% 
  filter(month(date) %in% c(10:12,1:3) & year(date) %in% 1980:2021) %>% 
  select(-date)

huc8 <- st_read('_data/WBD/WBDHU8.shp', quiet = TRUE)

Find AR days and sequence days

if (progress) pb <- txtProgressBar(min = 0, max = ncell(grid_ca), style = 3)
datelength <- length(dates_merra[month(dates_merra) %in% c(10:12,1:3) & year(dates_merra) %in% 1980:2021])

ar_binary <-
  foreach (i = 1:ncell(grid_ca), .combine = 'cbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (i %in% index_ca) {
      df_24hr[[i]] %>% 
        filter(month(date) %in% c(10:12,1:3) & year(date) %in% 1980:2021) %>% 
        select(ar)
    } else rep(NA, datelength)
  }  
seq_binary <-
  foreach (i = 1:ncell(grid_ca), .combine = 'cbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (i %in% index_ca) {
      df_24hr[[i]] %>% 
        filter(month(date) %in% c(10:12,1:3) & year(date) %in% 1980:2021) %>% 
        select(seq)
    } else rep(NA, datelength)
  }

Match HUC8 subbasins to MERRA-2 cells

huc_intersects <- huc8 %>% 
  st_intersects(grid_ca %>% rasterToPolygons %>% st_as_sf %>% st_transform(st_crs(huc8)))
# ggplot() + 
#   geom_tile(data = raster.df(grid_ca) %>% filter(!is.na(value)), aes(x=x, y=y), fill = NA, color = 'blue') + 
#   geom_sf(data = huc8, fill = NA, color = 'darkred')

Create runoff.extreme dataframe

if (progress) pb <- txtProgressBar(min = 0, max = nrow(huc8), style = 3)
runoff.extreme <- 
  foreach (h = 1:nrow(huc8), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,h)
    if (length(huc_intersects[h][[1]]) == 0) {
      return(rep(NA,4))
    } else if (length(huc_intersects[h][[1]]) == 1) {
      temp <- 
        data.frame(
          runoff = unlist(df.runoff[,names(df.runoff) == huc8$huc8[h]]),
          ar = ar_binary[,index_ca[huc_intersects[h][[1]]]],
          seq = seq_binary[,index_ca[huc_intersects[h][[1]]]])
    } else {
      temp <- 
        data.frame(
          runoff = unlist(df.runoff[,names(df.runoff) == huc8$huc8[h]]),
          ar = ar_binary[,index_ca[huc_intersects[h][[1]]]] %>% apply(1, function(x) mean(x)>0.1),
          seq = seq_binary[,index_ca[huc_intersects[h][[1]]]] %>% apply(1, function(x) mean(x)>0.1))
    }
    temp <- temp %>% 
      count(ar, seq, runoff.extreme = runoff > quantile(runoff, 0.95, na.rm = TRUE)) %>% 
      group_by(ar, seq) %>%
      mutate(p.all = sum(n)) %>% 
      ungroup %>% 
      mutate(p.all = p.all/sum(p.all)*2) %>% 
      filter(runoff.extreme) %>% 
      mutate(p.runoff = n/sum(n))
    c(temp %>% filter(ar) %>% summarize(across(starts_with('p.'), sum)) %>% unlist,
      temp %>% filter(seq) %>% summarize(across(starts_with('p.'), sum)) %>% unlist) %>% 
      setNames(c('ar.all', 'ar.runoff', 'seq.all', 'seq.runoff'))
  }

Plot results

huc8.plot <- huc8 %>% 
  select(huc8) %>% 
  cbind(runoff.extreme) %>% 
  pivot_longer(c(-huc8,-geometry)) %>% 
  separate(name, into = c('arseq', 'extreme'), sep = '\\.') %>% 
  filter(!is.na(value))
huc8.union <- st_union(huc8.plot)
g.runoff <- huc8.plot %>% 
  filter(extreme == 'runoff') %>% 
  ggplot() +
  geom_sf(aes(fill = value), size = 0.2) + 
  geom_sf(data = cal, fill = NA, color = 'grey70', size = 0.2) + 
  geom_sf(data = huc8.union, fill = NA, size = 0.2) + 
  facet_grid(arseq ~ extreme) + 
  scale_fill_scico(
    'Percentage of Days     ', palette = 'lapaz', direction = -1, 
    limits = c(0,1), labels = percent) + 
  coord_sf(xlim = layer_scales(g.precip)$x$range$range, ylim = layer_scales(g.precip)$y$range$range) +
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    axis.title = element_blank(),
    axis.text.x = element_text(angle = 90, vjust = 0.5),
    axis.text.y = element_blank(), axis.line.y = element_blank(), axis.ticks.y = element_blank(),
    strip.background.y = element_blank(), strip.text.y = element_blank(),
    plot.background = element_rect(fill = 'white', color = NA),
    strip.background.x = element_blank(), strip.text.x = element_blank(),
    plot.margin = margin(5,2,5,2),
    legend.position = 'bottom',
    legend.key.width = unit(0.5, 'in'), 
    legend.margin = margin(0,0,0,0))
# g.runoff

Soil moisture

Create dataframe

Load file names

# sm_filelist <- list.files('_data/WLDAS/files', full.names = TRUE)
# sm_datelist <- str_sub(sm_filelist, start = 28, end = -13) %>% ymd
# 
# sm_filelist <- sm_filelist[month(sm_datelist) %in% c(10:12,1:3) & year(sm_datelist) %in% 1981:2021]
# sm_datelist <- sm_datelist[month(sm_datelist) %in% c(10:12,1:3) & year(sm_datelist) %in% 1981:2021]

Subset to California

# sm_raster <- sm_filelist[1] %>% rast %>% terra::crop(extent(grid_ca)) %>% terra::mask(california)
# sm_match <- grid_ca %>% rast %>% 
#   project(sm_raster, method = 'near') %>%
#   .[] %>% c %>% 
#   cbind(sm_raster %>% terra::xyFromCell(1:ncell(.))) %>% 
#   data.frame %>% setNames(c('grid_merra', 'x', 'y')) %>% 
#   mutate(
#     value_sm = c(sm_raster[]),
#     grid_sm = 1:nrow(.), 
#     grid_merra = ifelse(is.nan(grid_merra), NA, grid_merra)) %>% 
#   mutate(
#     grid_sm = ifelse(is.na(value_sm), NA, grid_sm),
#     valid = !is.na(grid_merra) & !is.na(grid_sm)) %>% 
#   select(-value_sm)
# sm_valid <- sm_match$valid

Extract soil moisture values

# if (progress) pb <- txtProgressBar(min = 0, max = length(sm_filelist), style = 3)
# df.sm <- 
#   foreach (d = 1:length(sm_filelist), .combine = 'cbind') %do% {
#     if (progress) setTxtProgressBar(pb,d)
#     sm_filelist[d] %>% rast %>% terra::crop(extent(grid_ca)) %>% 
#       .[] %>% c %>% .[sm_valid]
#   } %>% t
# 
# ## gut check
# dim(df.sm)  
# c(length(sm_datelist), sum(sm_valid))

Find AR days and sequence days

# if (progress) pb <- txtProgressBar(min = 0, max = ncell(grid_ca), style = 3)
# ar_binary <-
#   foreach (i = 1:ncell(grid_ca), .combine = 'cbind') %do% {
#     if (progress) setTxtProgressBar(pb,i)
#     if (i %in% index_ca) {
#       df_24hr[[i]] %>% 
#         filter(date %in% sm_datelist) %>% 
#         select(ar)
#     } else rep(NA, length(sm_datelist))
#   }  
# seq_binary <-
#   foreach (i = 1:ncell(grid_ca), .combine = 'cbind') %do% {
#     if (progress) setTxtProgressBar(pb,i)
#     if (i %in% index_ca) {
#       df_24hr[[i]] %>% 
#         filter(date %in% sm_datelist) %>% 
#         select(seq)
#     } else rep(NA, length(sm_datelist))
#   }

Create sm.extreme dataframe

# if (progress) pb <- txtProgressBar(min = 0, max = length(sm_valid), style = 3)
# cl <- makeCluster(cores)
# registerDoSNOW(cl)
# sm.extreme <- 
#   foreach (
#     k = 1:length(sm_valid), 
#     .options.snow = if (progress) opts,
#     .packages = c('tidyverse'),
#     .combine = 'rbind') %dopar% {
#       if (sm_valid[k]) {
#         temp <- 
#           data.frame(
#             sm = unname(df.sm[,match(k,sm_match$grid_sm[sm_valid])]),
#             ar = ar_binary[,sm_match$grid_merra[k]],
#             seq = seq_binary[,sm_match$grid_merra[k]]) %>% 
#           count(ar, seq, sm.extreme = sm > quantile(sm, 0.95, na.rm = TRUE)) %>% 
#           group_by(ar, seq) %>%
#           mutate(p.all = sum(n)) %>% 
#           ungroup %>% 
#           mutate(p.all = p.all/sum(p.all)*2) %>% 
#           filter(sm.extreme) %>% 
#           mutate(p.sm = n/sum(n))
#         c(temp %>% filter(ar) %>% summarize(across(starts_with('p.'), sum)) %>% unlist,
#           temp %>% filter(seq) %>% summarize(across(starts_with('p.'), sum)) %>% unlist) %>% 
#           setNames(c('ar.all', 'ar.sm', 'seq.all', 'seq.sm'))
#       } else rep(NA,4)
#     }
# stopCluster(cl)

# ## save out
# save(
#   sm_filelist, sm_datelist, sm_match, df.sm, sm.extreme,
#   file = '_data/ARTMIP/sm_0905.Rdata')

Load from checkpoint

load('_data/ARTMIP/sm_0905.Rdata')

Plot results

g.sm <- sm_match %>% 
  cbind(sm.extreme) %>% filter(valid) %>% 
  select(x, y, starts_with('ar'), starts_with('seq')) %>%
  pivot_longer(c(-x,-y)) %>% 
  separate(name, into = c('arseq', 'extreme'), sep = '\\.') %>% 
  filter(extreme == 'sm') %>% 
  ggplot() + 
  geom_raster(aes(x=x, y=y, fill = value)) + 
  geom_sf(data = california, fill = NA, size = 0.2) + 
  facet_grid(arseq ~ extreme) + 
  scale_fill_scico(
    'Percentage of Days     ', palette = 'lapaz', direction = -1, 
    limits = c(0,1), labels = percent) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    axis.title = element_blank(),
    axis.text.x = element_text(angle = 90, vjust = 0.5),
    axis.text.y = element_blank(), axis.line.y = element_blank(), axis.ticks.y = element_blank(),
    strip.background.y = element_blank(), strip.text.y = element_blank(),
    plot.background = element_rect(fill = 'white', color = NA),
    strip.background.x = element_blank(), strip.text.x = element_blank(),
    plot.margin = margin(5,2,5,2),
    legend.position = 'bottom',
    legend.key.width = unit(0.5, 'in'), 
    legend.margin = margin(0,0,0,0))
# g.sm

FIG 2: Hydrologic impacts of sequences

Plot all AR & sequence days

if (progress) pb <- txtProgressBar(min = 0, max = ncell(grid_ca), style = 3)
arseq <- 
  foreach (i = 1:ncell(grid_ca), .combine = 'rbind') %do% {
    if (progress) setTxtProgressBar(pb,i)
    if (i %in% index_ca) {
      c('ar.all' = df_24hr[[i]] %>% 
          filter(month(date) %in% c(10:12,1:3) & year(date) %in% 1981:2021) %>% 
          count(ar) %>% mutate(pct = n/sum(n)) %>% filter(ar) %>% pull(pct),
        'seq.all' = df_24hr[[i]] %>% 
          filter(month(date) %in% c(10:12,1:3) & year(date) %in% 1981:2021) %>% 
          count(seq) %>% mutate(pct = n/sum(n)) %>% filter(seq) %>% pull(pct))
    } else rep(NA,2)
  }

Combine all plots into full figure

Note: Additional labels were added to this figure in Inkscape.

g.all <- raster.df(grid_ca) %>% cbind(arseq) %>% 
  filter(!is.na(value)) %>% select(-value) %>% 
  pivot_longer(ends_with('all'), names_to = 'arseq') %>% 
  mutate(empty = NA) %>% 
  ggplot() + 
  geom_raster(aes(x=x, y=y, fill = value)) + 
  facet_grid(arseq ~ empty) +
  geom_sf(data = california, fill = NA, size = 0.2) + 
  scale_fill_scico(
    'Percentage of Days     ', palette = 'lapaz', direction = -1, 
    limits = c(0,1), labels = percent) + 
  theme(
    text = element_text(family = 'Segoe UI', size = 8),
    axis.title = element_blank(),
    axis.text.x = element_text(angle = 90, vjust = 0.5),
    strip.background.y = element_blank(), strip.text.y = element_blank(),
    plot.background = element_rect(fill = 'white', color = NA),
    strip.background.x = element_blank(), strip.text.x = element_blank(),
    plot.margin = margin(5,2,5,2),
    legend.position = 'bottom',
    legend.key.width = unit(0.5, 'in'), 
    legend.margin = margin(0,0,0,0))
g.all + g.precip + g.runoff + g.sm + guide_area() + 
  plot_layout(design = 'abcd\neeee', heights = c(10,1), guides = 'collect')
ggsave('_figures/fig2_new.png', width = 5.5, height = 3.5, units = 'in', dpi = 300)

Numerical results

AR days & sequence days

Statewide wet-season AR day percentage:

pct_ar <- arseq %>% data.frame %>% pull(ar.all)

## print mean & standard deviation
cat(paste0(
  'Distribution Mean: ', 
  Mean(pct_ar*100) %>% round(1), '%\n'))
cat(paste0(
  'Distribution Standard Deviation: ', 
  sd(pct_ar*100, na.rm = TRUE) %>% round(1), '%'))
## Distribution Mean: 8.8%
## Distribution Standard Deviation: 5.6%


Statewide wet-season sequence day percentage:

pct_seq <- arseq %>% data.frame %>% pull(seq.all)

## print mean & standard deviation
cat(paste0(
  'Distribution Mean: ', 
  Mean(pct_seq*100) %>% round(1), '%\n'))
cat(paste0(
  'Distribution Standard Deviation: ', 
  sd(pct_seq*100, na.rm = TRUE) %>% round(1), '%'))
## Distribution Mean: 16.4%
## Distribution Standard Deviation: 11.7%


How many days are AR days and sequence days in the North Coast?

temp <- 
  grid_ca %>% raster.df %>% 
  cbind(arseq) %>% select(-value) %>% 
  rasterFromXYZ(crs = crs(grid_ca)) %>% 
  raster::extract(huc4 %>% st_transform(crs(grid_ca)), fun = mean, na.rm = TRUE) %>% 
  cbind(huc4 %>% st_drop_geometry %>% select(huc4, name))

cat('Maximum percentage of AR days and sequence days: \n')
temp %>% 
  filter(ar.all == max(ar.all) | seq.all == max(seq.all)) %>% 
  transmute(Name = name, AR.pct = percent(ar.all), seq.pct = percent(seq.all))
cat('\n')

cat('Minimum percentage of AR days and sequence days: \n')
temp %>% 
  filter(ar.all == min(ar.all) | seq.all == min(seq.all)) %>% 
  transmute(Name = name, AR.pct = percent(ar.all), seq.pct = percent(seq.all))
## Maximum percentage of AR days and sequence days: 
##                                  Name AR.pct seq.pct
## 1 Klamath-Northern California Coastal    18%     31%
## 
## Minimum percentage of AR days and sequence days: 
##                        Name AR.pct seq.pct
## 1 Northern Mojave-Mono Lake     3%      5%

Extreme precipitation


Summary of extreme precipitation days:

temp <- 
  precip_match %>% cbind(precip.extreme) %>% 
  select(x, y, ar.precip, seq.precip) %>% 
  rasterFromXYZ(crs = crs(grid_ca)) %>% 
  raster::extract(huc4 %>% st_transform(crs(grid_ca)), fun = mean, na.rm = TRUE) %>% 
  cbind(huc4 %>% st_drop_geometry %>% select(huc4, name))

temp %>% 
  mutate(diff = ar.precip - seq.precip) %>% 
  gt %>% 
  tab_header('Percentage of Extreme Precipitation Days') %>% 
  fmt_percent(c(ar.precip, seq.precip, diff), decimals = 1) %>% 
  cols_label(
    huc4 = 'Hydrologic Region',
    name = 'Name', 
    ar.precip = 'AR Days',
    seq.precip = 'Sequence Days',
    diff = 'Difference') %>% 
  tab_options(
    heading.background.color = '#d9d9d9', 
    column_labels.background.color = '#f2f2f2')

# lab <- glue("{temp$name}: {percent(temp$ar.precip)}")
# mapview(cbind(huc4, temp), zcol = 'ar.precip', label = lab)
# 
# lab <- glue("{temp$name}: {percent(temp$seq.precip)}")
# mapview(cbind(huc4, temp), zcol = 'seq.precip', label = lab)
Percentage of Extreme Precipitation Days
AR Days Sequence Days Hydrologic Region Name Difference
67.1% 67.4% 1801 Klamath-Northern California Coastal −0.3%
57.9% 55.1% 1802 Sacramento 2.8%
30.6% 29.3% 1803 Tulare-Buena Vista Lakes 1.3%
42.9% 43.7% 1804 San Joaquin −0.8%
59.3% 71.8% 1805 San Francisco Bay −12.5%
41.5% 56.3% 1806 Central California Coastal −14.8%
26.4% 32.0% 1807 Southern California Coastal −5.6%
34.9% 26.9% 1808 North Lahontan 8.0%
18.4% 11.4% 1809 Northern Mojave-Mono Lake 7.0%
21.8% 14.7% 1810 Southern Mojave-Salton Sea 7.1%

Extreme runoff


Summary of extreme runoff days:

temp <- 
  huc8 %>% cbind(runoff.extreme) %>% 
  st_drop_geometry %>% 
  select(huc8, ar.runoff, seq.runoff) %>% 
  group_by(huc4 = toNumber(str_sub(huc8, end = 4))) %>% 
  summarize(across(-huc8, Mean)) %>% 
  left_join(huc4 %>% st_drop_geometry %>% select(huc4, name), by = 'huc4')

temp %>% 
  mutate(diff = ar.runoff - seq.runoff) %>%
  gt %>% 
  tab_header('Percentage of Extreme Runoff Days') %>% 
  fmt_percent(c(ar.runoff, seq.runoff, diff), decimals = 1) %>% 
  cols_label(
    huc4 = 'Hydrologic Region',
    name = 'Name', 
    ar.runoff = 'AR Days',
    seq.runoff = 'Sequence Days',
    diff = 'Difference') %>% 
  tab_options(
    heading.background.color = '#d9d9d9', 
    column_labels.background.color = '#f2f2f2')

# lab <- glue("{temp$name}: {percent(temp$ar.runoff)}")
# mapview(cbind(huc4, temp), zcol = 'ar.runoff', label = lab)
# 
# lab <- glue("{temp$name}: {percent(temp$seq.runoff)}")
# mapview(cbind(huc4, temp), zcol = 'seq.runoff', label = lab)
Percentage of Extreme Runoff Days
Hydrologic Region AR Days Sequence Days Name Difference
1801 57.5% 90.1% Klamath-Northern California Coastal −32.7%
1802 43.6% 84.7% Sacramento −41.1%
1803 23.8% 53.6% Tulare-Buena Vista Lakes −29.7%
1804 26.0% 72.6% San Joaquin −46.6%
1805 56.3% 91.5% San Francisco Bay −35.2%
1806 38.8% 85.1% Central California Coastal −46.4%
1807 38.2% 59.2% Southern California Coastal −20.9%
1808 20.4% 39.1% North Lahontan −18.7%
1809 19.9% 35.1% Northern Mojave-Mono Lake −15.2%
1810 23.8% 47.4% Southern Mojave-Salton Sea −23.7%

Extreme soil moisture


Summary of extreme soil moisture days:

temp <- 
  sm_match %>% cbind(sm.extreme) %>% 
  select(x, y, ar.sm, seq.sm) %>% 
  rasterFromXYZ(crs = crs(grid_ca)) %>% 
  raster::extract(huc4 %>% st_transform(crs(grid_ca)), fun = mean, na.rm = TRUE) %>% 
  cbind(huc4 %>% st_drop_geometry %>% select(huc4, name))

temp %>% 
  mutate(diff = ar.sm - seq.sm) %>% 
  gt %>% 
  tab_header('Percentage of Extreme Soil Moisture Days') %>% 
  fmt_percent(c(ar.sm, seq.sm, diff), decimals = 1) %>% 
  cols_label(
    huc4 = 'Hydrologic Region',
    name = 'Name', 
    ar.sm = 'AR Days',
    seq.sm = 'Sequence Days',
    diff = 'Difference') %>% 
  tab_options(
    heading.background.color = '#d9d9d9', 
    column_labels.background.color = '#f2f2f2')

# lab <- glue("{temp$name}: {percent(temp$ar.sm)}")
# mapview(cbind(huc4, temp), zcol = 'ar.sm', label = lab)
# 
# lab <- glue("{temp$name}: {percent(temp$seq.sm)}")
# mapview(cbind(huc4, temp), zcol = 'seq.sm', label = lab)
Percentage of Extreme Soil Moisture Days
AR Days Sequence Days Hydrologic Region Name Difference
33.1% 72.5% 1801 Klamath-Northern California Coastal −39.4%
25.0% 59.5% 1802 Sacramento −34.5%
11.9% 36.2% 1803 Tulare-Buena Vista Lakes −24.3%
16.7% 47.7% 1804 San Joaquin −31.0%
27.8% 73.5% 1805 San Francisco Bay −45.7%
20.3% 66.8% 1806 Central California Coastal −46.5%
10.9% 37.1% 1807 Southern California Coastal −26.2%
12.3% 25.5% 1808 North Lahontan −13.2%
5.8% 11.4% 1809 Northern Mojave-Mono Lake −5.6%
7.0% 14.0% 1810 Southern Mojave-Salton Sea −7.0%