In mikejohnson51/nwmHistoric: National Water Model Historic Data

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4,
  #eval=nzchar(Sys.getenv("BUILD_VIGNETTES"))
  warning = FALSE,
  message = FALSE,
  error  = TRUE
)

oldoption <- options(scipen = 9999)
options(scipen = 9999)

library(ggplot2)
library(viridis)
library(gridExtra)
library(nwmHistoric)

Feature Identification

The NWM reanalysis product is indexed by feature ID. Identifying ID's is not always easy. Tools exist in dataRetirval and nhdplusTools to aid in feautre discovery.

Point Search

NHD COMIDs can be determined from point locations as the NHD catchment containing the requested point. discover_nhd can help find the NHDPlus catchment containing a defined point:

library(AOI)
library(dataRetrieval)
library(nhdplusTools)
library(nwmHistoric)

# Create POINT object by geocoding "Goleta"
(pt = AOI::geocode("Goleta", pt = T))

# Find COMID
(id = dataRetrieval::findNLDI(location = pt))

Area Search

discover_nhd can also accept polygon features as a search constraint, and will return all NHD catchments contained in or intersecting the input polygon using the National Water Census Geoserver.

# Goleta footprint
(AOI = AOI::aoi_get("Santa Barbara"))

# Find COMIDS
(ids = nhdplusTools::get_nhdplus(AOI)$comid)

Finding NWIS gages

Analogous functionality for finding USGS NWIS site IDs is provided with discover_nwis. This function only returns NWIS sites that record streamflow (USGS parameter code '00060') and are collocated with an NHD catchment represented in the NWM.

(nwis = nhdplusTools::get_nwis(AOI))

Geometric Discovery

The NHDPlusV2 data model loosely conforms to the HY_Features Conceptual Model with a mapping shared here. nwmHistoric extends the ability of discover_nhd to retrieve catchment-divide, flowline, and outlet represnetations of a COMID.

# Return catchments
catch = nhdplusTools::get_nhdplus(AOI, realization = "catchment")

# Return flowlines
fl    = nhdplusTools::get_nhdplus(AOI, realization = "flowline")

# Return outlets
out   = nhdplusTools::get_nhdplus(AOI, realization = "outlet")

library(ggplot2)

ggplot() +
  geom_sf(data = catch, aes(fill = "Catchment"), color = "black") +
  geom_sf(data = fl,  aes(color = 'Flowline'), show.legend  = "line") +
  geom_sf(data = out, aes(color = 'Outlet'), show.legend  = "point") + 
  geom_sf(data = nwis, aes(color = 'NWIS'),show.legend = "point")+
  scale_fill_manual(values = c("Catchment" = "gray99"), 
                    guide = guide_legend(override.aes = list(linetype = 1, border = "solid", shape = NA))) +
  scale_colour_manual(values = c("Flowline" = "blue","Outlet" = "red", "NWIS" = "darkgreen"), guide = guide_legend(override.aes = list(linetype = c("solid", "blank", "blank"), shape = c(NA,16,16)))) + theme_minimal() +
  labs(fill = "", color = "") + 
  theme_void() +
  theme(legend.position = 'bottom')

Putting it all together

Lets look at one integrative example. The aim is to identify a self-contained (e.g watershed) catchment and extract reanalysis records along the mainstem for 2015. The dataRetrival package provides an interface to the Network Linked Data Index which allow users to define a starting point (as a USGS siteID or NHDPlus COMID), and traverse the hydrographic network to find the upstream (or downstream) geometries and indexed elements.

Here we use the NLDI to trace the upstream network and the nwmHistoric to extract the relevant streamflow forecasts.

loc = findNLDI(nwis = '05428500', 
               nav = c("UT", "UM"), 
               find = c("flowline", "basin"))

# Find 2015 Flows along the Mainstem
nldi_flows <- readNWMdata(comid = loc$UM_flowlines$nhdplus_comid, 
                          startDate = "2015-01-01", 
                          endDate   = "2015-12-31")

# nwmHistric find_outlets functions finds outlets from a flowpath representation
um_outlets = nhdplusTools::get_node(loc$UM_flowlines)
um_outlets$order = factor(nrow(um_outlets):1)

# Make a map
g1 = ggplot() +
  geom_sf(data = loc$basin, col = 'gray90') +
  geom_sf(data = loc$UT_flowlines, col = 'lightblue') + 
  geom_sf(data = loc$UM_flowlines, col = 'blue', size = .5) + 
  geom_sf(data = um_outlets, aes(color = order), size = 1.25) + 
  scale_color_manual(values= viridis::viridis(24)) +
  theme_void() + theme(legend.position="none")


# Plots the flows
g2 = ggplot(data = nldi_flows, aes(x = time_utc, y = flow, color = comid)) +
  geom_line(size = .25) + 
  scale_color_manual(values= viridis::viridis(24)) +
  ggpubr::theme_pubr() +
  labs(y = "Qsim (cms)",
       x = paste0("DateTime (", base::format(nldi_flows$time[1], format="%Z"), ")"), 
       title = paste(nldi_flows$model[1],"Flow Records")) +
  theme(legend.position = 'NA')

gridExtra::grid.arrange(g1,g2, nrow = 1)