R for Water Resources









Dave Blodgett and Laura DeCicco
2019 AWRA National Meeting, 11/4/2019
Cutthroat trout caught in the Weber River. (credit: Corey Kruitbosch
Use arrow keys or click the arrows to advance.

Introduction

These slides provide resources and examples for two common water resources applications in the R programming language.

A companion workflow to get all the data and perform processing operations is available here.


  • Introduction to examples
  • Background and resources
  • Data access
  • Data visualization
  • Spatial intersection
  • Water budget analysis
  • Trend analysis
  • Synthesis

Application 1:
Water Budget

Access and compare precip, actual ET, and streamflow.
plot of chunk wb_example_map

plot of chunk wb_example_plot Data sourced from web services using R packages shown later.

Application 2: Trends and Plotting

plot of chunk trend2

Exploration and Graphics for RivEr Trends (EGRET): A package for analysis of long-term changes in water quality and streamflow. plot of chunk trend1

Background:

Resources:

Spatial Data Access

See a nhdplusTools tutorial here.

library(nhdplusTools)
site <- list(featureSource = "nwissite", 
             featureID = "USGS-10128500")

line <- navigate_nldi(site, "UT", "")
site <- navigate_nldi(site, "UT", "nwissite")

nhdp <- subset_nhdplus(ut$nhdplus_comid, 
                       "nhdp_subset.gpkg", 
                       "download")

plot of chunk plot simple

Observations Data Access

Get USGS and EPA water data. See a dataRetrieval tutorial here.

library(dataRetrieval)
flow <- readNWISdv("10128500", "00060")

plot of chunk observationsRun

Spatial data visualization can be accomplished in several ways. This example uses base R plotting and the prettymapr and rosm packages. Output is shown on the next slide.

prettymap(title = paste("NHDPlus and NLDI data for the", outlet_name), 
          scale.label.cex = 2, scale.padin = c(0.25, 0.25),
          drawarrow = TRUE, arrow.scale = 2,
          mai = c(0.5, 0, 0, 0), { # margin for legend 
  osm.plot(nhd_bbox, type = "cartolight", quiet = TRUE, progress = "none")
  plot(gt(nhd_cat), lwd = 0.5, col = NA, border = "grey", add = TRUE)
  plot(gt(nhd_basin), lwd = 1, col = NA, border = "black", add = TRUE)
  plot(gt(nhd_fline), lwd = streamorder, col = "blue", add = TRUE)
  plot(gt(nhd_area), col = "lightblue", border = "lightblue", add = TRUE)
  plot(gt(nhd_wbody), col = "lightblue", border = "lightblue",add = TRUE)
  plot(gt(wqpsite), col = "red", pch = 46, cex = 1.25, add = TRUE)
  plot(gt(nwissite), col = "black", bg = "lightgrey", pch = 24, add = TRUE)})

See the project source code for legend code.
While extremely configurable, base-R plotting can be tedious. ggplot2 offers a different approach to plotting that some find more convenient.

plot of chunk spatial_map

Index Sites to Flowlines

nhdplusTools can attach sites to NHD flowlines using a nearest neighbor search with sites and geometry nodes.

nhdplusTools::get_flowline_index(
    points = sites,
    flines = flowline, 
    search_radius = 100, # units of flowlines
    precision = 10), # maximum node spacing
COMID REACHCODE REACH_meas offset
10093110 16020101000048 44.3 48.2
10273232 16020102000008 87.8 7.0
10274616 16020102000017 48.7 33.5
10274270 16020102000135 45.7 22.9
10092262 16020101000948 17.9 33.5

Spatial Intersection

The code below shows the anatomy of an intersectr workflow. It uses an OPeNDAP service to access remote data the same way it could a local NetCDF file. 

library(intersectr)
cells <- create_cell_geometry(x_coords, y_coords, nc_projection, catchment)

poly <- sf::st_as_sf(dplyr::select(catchment, ID = featureid))

weights <- calculate_area_intersection_weights(cells, poly)

dap_uri <- "https://cida.usgs.gov/thredds/dodsC/UofIMETDATA"
execute_intersection(nc_file = dap_uri,
                     variable_name = "precipitation_amount",
                     intersection_weights = weights,
                     cell_geometry = data_source_cells,
                     x_var = x_coord, y_var = y_coord, t_var = t_coord,
                     start_datetime = "2009-10-01 00:00:00", 
                     end_datetime = "2010-10-01 00:00:00")

intersectr is a highly scalable and adaptible workflow component. It is still an “in development” package.

plot of chunk intersectr_plotplot of chunk intersectr_plot

Precipitation sourced from gridMet
Evapotranspiration sourced from SSEBop ETa

Results of intersection shown individually on the left and differenced below.
Runoff and Infiltration?

plot of chunk intersectr_plot_real

Water Budget Analysis

plot of chunk wb

Trend Analysis with EGRET

EGRET includes a function to calculate a weighted regrestion on time, discharge, and season (WRTDS).

library(EGRET)
Daily <- readNWISDaily(siteNumber = "10128500")

Sample = readNWISSample(siteNumber = "10128500",
                        parameterCd = "00095",
                        endDate = "2018-01-01")

INFO = readNWISInfo(siteNumber = "10128500",
                    parameterCd = "00095",
                    interactive = FALSE)

eList = mergeReport(INFO = INFO,
                    Daily = Daily,
                    Sample = Sample) %>% 
    modelEstimation()

Gets daily flow, sampled concentration, and site information. mergeReport prepares data for modelEstimation which performs a Weighted Regrerssions on Time, Discharge, and Season model fit, returning timeseries and flow-concentration surfaces.

This plot is the WRTDS model result of specific conductance concentration over time and discharge. plot of chunk summary

This plot uses the model to show how June concentration changes over time at specific discharge conditions. plot of chunk summary2

Summary

plot of chunk summary3plot of chunk summary3

Questions?

Photo credit: Ricraider