In joshualerickson/whitewater: Parallel Processing Options for Package 'dataRetrieval'

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
   fig.align='center',
  fig.height=6, 
  fig.width=8,
  dev.args = list(png = list(type = "cairo"))
)

ggplot2::theme_set(ggplot2::theme_bw())

library(whitewater)
library(ggplot2)
library(dplyr)
library(tidyr)
library(sf)
library(future)
library(ggfx)

Intro

Sometimes you just want to compare the current flow (water year or past week, whatever you want) to historical flows (figure below). The ww_statsUSGS() function does this for you! It takes the historical values for your parameter (flow in this example) using the USGS Statistics Web Service beta and returns percentiles for each site but also combines current, monthly or yearly values and conditions depending on your choice of temporalFilter. Using the web service is really helpful because it quality controls the incoming data, i.e. data that is incomplete or not qc'd will not be returned. In this vignette we'll go over the different options/arguments in the ww_statsUSGS() as well as some ways to visualize the results.

yaak_dv <- ww_dvUSGS(sites = '12304500', parameter_cd = '00060')

yaak_daily_report <- yaak_dv %>% ww_statsUSGS(temporalFilter = 'daily',
                                  days = 365)
yaak_daily_report %>% 
  pivot_longer(c('Flow', 'p25_va', 'p75_va', 'max_va', 'min_va')) %>% 
  ggplot() + 
  with_outer_glow(geom_line(aes(Date, value, color = name, alpha = name %in% c('p25_va', 'p75_va', 'max_va', 'min_va'))), sigma = .1) + 
  scale_color_manual(values = c('black', 'blue', 'red', 'yellow', 'forestgreen')) +
  scale_alpha_manual(values = c(1,.5), guide = 'none') +
  labs(y = 'Discharge (cfs)', color = '', title = 'Comparing current flow to Max, Min, 25th and 75th percentiles') + 
  theme_bw()

ww_statsUSGS()

The ww_statsUSGS() is essentially a wrapper around dataRetrieval::readNWISstat() but in addition will take the parameter values ('Flow', 'Wtemp', etc) calculated from the ww_floorIVUSGS() function and add to the data as well. Thus, the instantaneous values will look different than the daily mean values, as it should. Just like dataRetrieval, there are three different flavors of temporal windows that we can choose from: daily, monthly and yearly.

daily

The temporalFilter argument is used to generate the window of percentiles. Like in the example below we just want to get the last 30 days of daily stats so we'll put temporalFilter = 'daily' and days = 30. The days arguments is only used for daily but will go back the amount of days you input.

yaak_daily_report_30 <- yaak_dv %>% ww_statsUSGS(temporalFilter = 'daily',
                                  days = 30)
yaak_daily_report_30

yaak_daily_report_30 %>% 
  pivot_longer(c('Flow', 'p25_va', 'p75_va', 'max_va', 'min_va')) %>% 
  ggplot() + 
  with_outer_glow(geom_line(aes(Date, value, color = name, alpha = name %in% c('p25_va', 'p75_va', 'max_va', 'min_va'))), sigma = .1) + 
  scale_color_manual(values = c('black', 'blue', 'red', 'yellow', 'forestgreen')) +
  scale_alpha_manual(values = c(1,.5), guide = 'none') +
  labs(y = 'Discharge (cfs)', color = '', title = 'Comparing current flow to Max, Min, 25th and 75th percentiles') + 
  theme_bw()

If we want to compare multiple sites that fine too! In addition, you don't have to pipe a ww_dvUSGS() object into the function either. You can just input the sites you want; however, I always recommend using ww_dvUSGS() because you'll likely come back to the daily values anyway but in the example below we'll just call with inputting the sites into the sites argument. If you're going to use the sites argument you'll also need to specify the parameter_cd!

multiple_sites <- ww_statsUSGS(sites = c('12304500', '14159200'),
                               parameter_cd = '00060',
                               temporalFilter = 'daily', 
                               days = 30)

multiple_sites %>% 
  rename(`Current Flow` = 'Flow', `25%` = 'p25_va', `75%` = 'p75_va', Max = 'max_va', Min = 'min_va') %>% 
  pivot_longer(c(`Current Flow`,`25%`, `75%`, 'Max', 'Min')) %>% 
  mutate(name = factor(name, levels = c('Current Flow', 'Max', 'Min', '25%', '75%'))) %>% 
  ggplot() +
  labs(y = 'Discharge (cfs)', color = '', title = 'Comparing current flow to Max, Min, 25th and 75th percentiles') + 
  with_outer_glow(geom_line(aes(Date, value, color = name, label = name, alpha = name %in% c('Current Flow','25%', '75%', 'Max', 'Min')),text_smoothing = 30,hjust = 1, show.legend = F), sigma = .5)+ 
  scale_color_manual(values = c('black', 'blue', 'red', 'yellow', 'forestgreen')) +
  scale_alpha_manual(values = c(1,.5), guide = 'none')  + 
  facet_wrap(~Station, scales = 'free')

What's nice about this is you can run in parallel using the argument parallel = TRUE (see Parallel vignette for more details on parallel processing). We'll use the ww_current_conditions() to get the current streamflow status for the lower 48 that has a description of All-time low for this day and then we'll use these sites to view the last 30 days!

current_condition <- ww_current_conditions() %>% filter(StatisticsStatusDescription %in% c('All-time low for this day'),!TimeZoneCode %in% c('AKDT', 'AST','HST', 'GST'))  

plan(multisession(workers = 10))

atl_sites <- ww_statsUSGS(sites = unique(current_condition$SiteNumber),
                             parameter_cd = '00060', 
                             temporalFilter = 'daily',
                             days = 30,
                             parallel = TRUE,
                             verbose = FALSE)

ranking_flow <- atl_sites %>% left_join(current_condition  %>% filter(!TimeZoneCode %in% c('AKDT', 'AST','HST', 'GST')) %>% 
                                             select(site_no = "SiteNumber",
                                                    StatisticsStatusDescription,
                                                    StatisticStatusCode,
                                                    StatisticsStatusColorFill,
                                                    StatisticsStatusColorStroke,
                                                    Latitude,
                                                    Longitude))

mind <- min(atl_sites$Date, na.rm = T)
maxd <- max(atl_sites$Date, na.rm = T)



p1 <- ranking_flow %>% 
  mutate(Flow = if_else(Flow == "NaN", NA_real_, Flow)) %>% 
  filter(!is.na(Flow), !is.na(StatisticsStatusDescription)) %>% 
  ggplot(aes(Date)) + 
  with_outer_glow(geom_line(aes(y = Flow), color = '#FF0000', show.legend = F)) +
  with_inner_glow(geom_ribbon(aes(ymin = p25_va, ymax = p75_va), fill = 'grey75', alpha = 0.2)) + 
  labs(y = 'Discharge (cfs)',
       color = '',
       title = 'Comparing current flow to 25th and 75th percentiles',
       subtitle = paste0('Dates: ', mind,' to ', maxd)) +
  facet_wrap(~site_no, scales = 'free') +
  theme_void() +
  theme(text=element_text(size=10,  family="Comic Sans MS"),
        strip.text = element_blank(),
        axis.title.x = element_text(),
        axis.title.y = element_text(angle = 90,margin = unit(c(0, 3, 0, 0), 'mm')),
        plot.subtitle = element_text(margin = unit(c(0, 0, 3, 0), 'mm')),
        plot.title = element_text(margin = unit(c(0, 0, 3, 0), 'mm')))
p1

p2 <- ranking_flow %>% 
  group_by(site_no) %>% 
  slice(n = 1) %>% 
  filter(!is.na(Longitude)) %>% 
  ungroup() %>% 
  st_as_sf(coords = c('Longitude', 'Latitude'), crs = 4326) %>% filter(!is.na(StatisticsStatusDescription)) %>% 
  ggplot() + 
  with_inner_glow(borders('state',fill = 'white'),sigma = 4) + 
geom_sf(color = '#FF0000', aes(color = StatisticsStatusDescription,
                               fill = StatisticsStatusDescription),
        stroke = .25, size = 1.25, shape = 21, show.legend = F) +
  scale_fill_manual(values = "#FF0000") + 
  scale_color_manual(values = "#990000") +
  theme_void() +
  labs(
       title = 'USGS Sites with Current Streamflow Status at "All-time low"',
       subtitle = paste0('Date: ', maxd))+
  theme(text=element_text(size=10,  family="Comic Sans MS"))
p2

monthly and yearly

Now we can looked into the other temporal windows like monthly and yearly. These work the same as daily but now we just need to change the argument.

yaak_monthly_stats <- ww_statsUSGS(yaak_dv, temporalFilter = 'monthly')

yaak_monthly_stats %>% 
  filter(month == 8) %>% 
  mutate(diff = mean_value-p50_va,
         Status = if_else(diff > 0, "Above normal", "Below normal")) %>% 
  ggplot(aes(year_nu, diff, fill = Status)) + 
  with_outer_glow(geom_col(), sigma = 20) +
  scale_fill_manual(values = c("#009900","#996600")) +
  labs(y = 'Discharge (cfs)', x = 'Year', color = '', title = 'Comparing monthly mean flow in August to 50% percentile') + 
  theme_dark()

This looks interesting! What about the sites that have all-time low for today? Well, just like above we change the temporalFilter argument to monthly instead of daily and then use those site numbers to get back a tibble with mean monthly flows and statistics for each month per year per site.

atl_sites_monthly <- ww_statsUSGS(sites = unique(current_condition$SiteNumber),
                             parameter_cd = '00060', 
                             temporalFilter = 'monthly',
                             parallel = TRUE,
                             verbose = FALSE)

We can then look at a few of the sites in a heatmap style graph.

monthly_stats <- atl_sites_monthly %>% add_count(site_no,year_nu) %>% filter(n == 12)

atl_sites_monthly %>% 
  group_by(site_no) %>%   
  add_count() %>% 
  filter(n > 950) %>% 
  nest() %>% 
  .[1:4,] %>% 
  unnest() %>% 
  ungroup() %>% 
  ggplot() + 
  with_inner_glow(geom_tile(aes(month, year_nu, fill = StatisticsStatusDescription))) +
  scale_x_continuous(breaks = scales::pretty_breaks()) + 
  scale_fill_manual('Streamflow: status',
                    values = c(
                      "All-time low for this month" = "#FF0000",
                      "Much below normal" = "#BB2222",
                      "Below normal" = "#FFAA00",
                      "Normal" = "#00ff00",
                      "Above normal" = "#44dddd",
                      "Much above normal" = "#00ffff",
                      "All-time high for this month" = "#000055")) + 
  labs(title = "Heatmap of Monthly Mean Streamflow", x = 'Month', y = 'Year')+
  facet_wrap(~Station, scales = 'free')+
  theme(text=element_text(size=10),
        strip.text = element_text(size = 5),
        axis.title.x = element_text(),
        axis.title.y = element_text(angle = 90,margin = unit(c(0, 3, 0, 0), 'mm')),
        plot.subtitle = element_text(margin = unit(c(0, 0, 3, 0), 'mm')),
        plot.title = element_text(margin = unit(c(0, 0, 3, 0), 'mm')))

Or, back to the Yaak, MT station.

yaak_monthly_stats %>% 
ggplot() + 
  with_inner_glow(geom_tile(aes(month, year_nu,
                                fill = StatisticsStatusDescription,
                               color = StatisticsStatusDescription),
                            show.legend = T)) +
  scale_color_manual('Streamflow: status',
                     values = c(
                                "All-time low for this year" = "#990000",
                                "Much below normal" = "#661111",
                                "Below normal" = "#996600",
                                "Normal" = "#009900",
                                "Above normal" = "#11aaaa",
                                "Much above normal" = "#000099",
                                "All-time high for this year" = "#000000"))+ 
  scale_fill_manual('Streamflow: status',
                    values = c(
                      "All-time low for this year" = "#FF0000",
                      "Much below normal" = "#BB2222",
                      "Below normal" = "#FFAA00",
                      "Normal" = "#00ff00",
                      "Above normal" = "#44dddd",
                      "Much above normal" = "#0000FF",
                      "All-time high for this year" = "#000055"))+ 
  labs(title = "Heatmap of Monthly Mean Streamflow", x = 'Month', y = 'Year')+
  facet_wrap(~Station, scales = 'free')

Same as monthly we can do this for yearly.

yaak_yearly_stats <- ww_statsUSGS(yaak_dv, temporalFilter = 'yearly')

Now we can look at the mean value for the year against the 25th and 75th percentiles.

yaak_yearly_stats %>% 
ggplot(aes(year)) +
  with_outer_glow(geom_ribbon(aes(ymin = p25_va, ymax = p75_va), fill = '#00ff00', alpha = 0.5), sigma = 10)+
  with_outer_glow(geom_line(aes(y = mean_value), color = 'black', show.legend = F),sigma = 2, expand = .1)+
  labs(y = 'Discharge (cfs)',
       color = '',
       title = 'Comparing yearly mean flow to 25th and 75th percentiles') +
  facet_wrap(~Station, scales = 'free') +
  theme(text=element_text(size=10),
        strip.text = element_blank(),
        axis.title.x = element_text(),
        axis.title.y = element_text(angle = 90,margin = unit(c(0, 3, 0, 0), 'mm')),
        plot.subtitle = element_text(margin = unit(c(0, 0, 3, 0), 'mm')),
        plot.title = element_text(margin = unit(c(0, 0, 3, 0), 'mm')))

Or look at it in a bar graph with fill as the status.

p <- yaak_yearly_stats %>% 
ggplot(aes(year)) +
  geom_col(aes(year, mean_value, fill = StatisticsStatusDescription,
                               color = StatisticsStatusDescription), show.legend = T) +
  scale_color_manual('Streamflow: status',
                     values = c(
                                "All-time low for this year" = "#990000",
                                "Much below normal" = "#661111",
                                "Below normal" = "#996600",
                                "Normal" = "#009900",
                                "Above normal" = "#11aaaa",
                                "Much above normal" = "#000099",
                                "All-time high for this year" = "#000000"))+ 
  scale_fill_manual('Streamflow: status',
                    values = c(
                      "All-time low for this year" = "#FF0000",
                      "Much below normal" = "#BB2222",
                      "Below normal" = "#FFAA00",
                      "Normal" = "#00ff00",
                      "Above normal" = "#44dddd",
                      "Much above normal" = "#0000FF",
                      "All-time high for this year" = "#000055"))+
  labs(y = 'Mean Discharge (cfs)',
       color = '',
       x = 'Year',
       title = paste0('Yearly mean flows for ', yaak_yearly_stats %>% slice(n = 1) %>% pull(Station)))

p

Wrap-up

Hopefully this gives you an idea of how to use whitewater with the ww_statsUSGS() function? If you have any questions or issues please submit a issue on github. Thanks!

joshualerickson/whitewater documentation built on July 1, 2024, 10:14 p.m.