Parallel Processing Options for Package 'dataRetrieval'.

whitewater

The goal of whitewater is to provide sequential and parallel processing for USGS stations in a tidy-style format. This package allows user to plan() their choice of parallel processing and then use the argument parallel = TRUE in whitewater function calls. The package also puts every output in a tibble with data munging of sites, parameter and stat codes, which results in a tidy style data frame.

Attention!

Due to potentially crashing (USGS Water Services)[https://waterservices.usgs.gov/] REST services parallel processing is kept to 120 requests/min. By following this rate limit, we can still benefit from parallel processing but also being mindful/respectful to the USGS Water Services via {dataRetrieval} and REST services. Thank you!

Installation

You can install the development version of whitewater from GitHub with:

# install.packages("devtools")
devtools::install_github("joshualerickson/whitewater")

Example

This is a basic example which shows you how to solve a common problem: get daily values of discharge for multiple sites (all active sites in Pacific Northwest (Region 17)) using parallel processing. Please see furrr and future for more details on parallel processing methods.

Running in parallel

library(whitewater)
library(tidyverse)
library(sf)
library(future)
library(dataRetrieval)

huc17_sites <- dataRetrieval::whatNWISdata(huc = 17,
                                           siteStatus = 'active',
                                           service = 'dv',
                                           parameterCd = '00060',
                                           drainAreaMax = 2000)
cat("# of sites: ", nrow(huc17_sites))
#> # of sites:  675

st_as_sf(huc17_sites, coords = c('dec_long_va', 'dec_lat_va')) %>% 
  ggplot() + 
  geom_sf() +
  borders('state', xlim = c(-130, -110), ylim = c(20, 50)) + 
  theme_bw()

#need to call future::plan()

##### Remember, please use 10 or less workers #####
plan(multisession(workers = 11))

#running on 11 cores

system.time({
pnw_dv <- ww_dvUSGS(huc17_sites$site_no,
                    parameter_cd = '00060',
                    wy_month = 10,
                    parallel = TRUE,
                    verbose = FALSE)
})
#>    user  system elapsed 
#>   39.96    2.47  559.56

nrow(pnw_dv)
#> [1] 11752061

pnw_dv
#> # A tibble: 11,752,061 x 28
#>    Station site_no drainage_area   lat  long altitude agency_cd Date        Flow
#>    <chr>   <chr>           <dbl> <dbl> <dbl>    <dbl> <chr>     <date>     <dbl>
#>  1 DONNER~ 103960~          124.  42.8 -119.     4254 USGS      1911-04-01   160
#>  2 DONNER~ 103960~          124.  42.8 -119.     4254 USGS      1911-04-02   166
#>  3 DONNER~ 103960~          124.  42.8 -119.     4254 USGS      1911-04-03   173
#>  4 DONNER~ 103960~          124.  42.8 -119.     4254 USGS      1911-04-04   180
#>  5 DONNER~ 103960~          124.  42.8 -119.     4254 USGS      1911-04-05   187
#>  6 DONNER~ 103960~          124.  42.8 -119.     4254 USGS      1911-04-06   232
#>  7 DONNER~ 103960~          124.  42.8 -119.     4254 USGS      1911-04-07   207
#>  8 DONNER~ 103960~          124.  42.8 -119.     4254 USGS      1911-04-08   182
#>  9 DONNER~ 103960~          124.  42.8 -119.     4254 USGS      1911-04-09   157
#> 10 DONNER~ 103960~          124.  42.8 -119.     4254 USGS      1911-04-10   132
#> # ... with 11,752,051 more rows, and 19 more variables: Flow_cd <chr>,
#> #   STILLING.WELL_Flow <dbl>, STILLING.WELL_Flow_cd <chr>,
#> #   .Estimated.By.Regression._Flow <dbl>,
#> #   .Estimated.By.Regression._Flow_cd <chr>, ..2.._Flow <dbl>,
#> #   ..2.._Flow_cd <chr>, .Discharge.1921.to.1952._Flow <dbl>,
#> #   .Discharge.1921.to.1952._Flow_cd <chr>, year <dbl>, month <dbl>, day <int>,
#> #   doy <dbl>, wy_doy <dbl>, month_day <chr>, wy <int>, month_abb <fct>, ...

Now we can use other ww_ functions to filter the data by water year, month, water year and month, as well as stat reporting (percentiles comparing current readings).

Water Year

Same as above, we can just call parallel = TRUE to run in parallel since we’ll be getting peak flows from dataRetrieval::readNWISpeak().

system.time({
pnw_wy <- suppressMessages(ww_wyUSGS(pnw_dv,
                                     parallel = TRUE,
                                     verbose = FALSE))
})
#>    user  system elapsed 
#>   16.86    0.72  418.25

pnw_wy
#> # A tibble: 32,663 x 27
#>    Station   site_no    wy peak_va peak_dt    drainage_area   lat  long altitude
#>    <chr>     <chr>   <int>   <dbl> <date>             <dbl> <dbl> <dbl>    <dbl>
#>  1 ABIQUA C~ 142007~  2014      NA NA                  45.3  45.0 -123.      193
#>  2 ABIQUA C~ 142007~  2015    3330 2014-12-21          45.3  45.0 -123.      193
#>  3 ABIQUA C~ 142007~  2016    5980 2015-12-07          45.3  45.0 -123.      193
#>  4 ABIQUA C~ 142007~  2017    3740 2017-03-09          45.3  45.0 -123.      193
#>  5 ABIQUA C~ 142007~  2018    4290 2017-10-22          45.3  45.0 -123.      193
#>  6 ABIQUA C~ 142007~  2019    3700 2019-04-07          45.3  45.0 -123.      193
#>  7 ABIQUA C~ 142007~  2020    2600 2020-01-28          45.3  45.0 -123.      193
#>  8 ABIQUA C~ 142007~  2021    5360 2020-12-20          45.3  45.0 -123.      193
#>  9 ABIQUA C~ 142007~  2022      NA NA                  45.3  45.0 -123.      193
#> 10 AHTANUM ~ 125025~  1904      NA NA                 107.   46.5 -120.      940
#> # ... with 32,653 more rows, and 18 more variables: obs_per_wy <int>,
#> #   wy_count <int>, Flow_sum <dbl>, Flow_max <dbl>, Flow_min <dbl>,
#> #   Flow_mean <dbl>, Flow_median <dbl>, Flow_stdev <dbl>, Flow_coef_var <dbl>,
#> #   Flow_max_dnorm <dbl>, Flow_min_dnorm <dbl>, Flow_mean_dnorm <dbl>,
#> #   Flow_med_dnorm <dbl>, Flow_max_sdnorm <dbl>, Flow_min_sdnorm <dbl>,
#> #   Flow_mean_sdnorm <dbl>, Flow_med_sdnorm <dbl>, Flow_sd_norm <dbl>

Without using parallel

If you just want a few sites (or one) and not use parallel processing, go for it! You’ll still get the advantages of filtering and stats. In addition, You don’t always have to pipe a ww_dvUSGS() object into the ww_*() and can just use the sites argument. In the example below I’ll do this but IMO its nice to start with a ww_dvUSGS() object because you’ll likely come back to it.

withlacoochee_temp_and_flow <- ww_wyUSGS(sites="02319394",
                          parameter_cd = c("00010", "00060"))
#> v 'water year' was successfully downloaded.
#> > now starting to gather peak flows using dataRetrieval::readNWISpeak
#> v 02319394 'peak flows' were successfully downloaded.

withlacoochee_temp_and_flow %>% 
  pivot_longer(c('Wtemp_max', 'Flow_max')) %>% 
ggplot(aes(wy, value)) +
  geom_point() + 
  geom_line() +
  theme_bw() +
  labs(x = 'Water Year') +
  facet_wrap(~name, scale = 'free')

Stats

Sometimes you just want to compare the current flow (water year or past week, whatever you want) to historical flows. The ww_statsUSGS() function does this for you! It takes the historical values for your parameter (flow in this example) and returns percentiles (dataRetrieval::readNWISstat()) but also combines the current values.

yaak_river_dv <- ww_dvUSGS('12304500')
#> v Yaak River near Troy MT 'daily' was successfully downloaded.

yaak_daily_report <- ww_statsUSGS(yaak_river_dv, 
                                  temporalFilter = 'daily',
                                  days = 365)
#> v Yaak River near Troy MT 'NWIS Stat' for Temporal Filter (daily) was successfully downloaded.
yaak_daily_report %>% 
  pivot_longer(c('Flow', 'p25_va', 'p75_va')) %>% 
  ggplot() + 
  geom_line(aes(Date, value, color = name, alpha = name %in% c('p25_va', 'p75_va'))) + 
  scale_color_manual(values = c('black', 'red', 'blue')) +
  scale_alpha_manual(values = c(1,.5), guide = 'none') +
  labs(y = 'Discharge (cfs)', color = '', title = 'Comparing current flow to 25th and 75th percentiles') + 
  theme_bw()