In USGS-R/dataRetrieval: Retrieval Functions for USGS and EPA Hydrology and Water Quality Data

library(knitr)

opts_chunk$set(
  warning = FALSE,
  message = FALSE,
  out.width = "75%",
  fig.align = "center",
  comment = "#>"
)

This post will demonstrate dataRetrieval functions to query the Network Linked Data Index (NLDI).

The NLDI provides a information backbone to navigate the NHDPlusV2 network and discover features indexed to the network. This process of feature discovery mirrors web-based navigation tools like Google Maps.

Navigation Analogy

Imagine you want to know about UCSB using Google Maps. In asking about the entity 'UCSB' we get information related to the feature (links, operating hours, photos) and things nearby.

knitr::include_graphics("ucsb.png")

Imagine then that you want to travel from UCSB to UCLA and know you need gas along the way. In Google Maps you might use UCSB as your origin and ask for directions to UCLA that include gas stations along your route.

include_graphics("nldi-gas-station-example.png")

The NLDI works to provide the same functionality for the hydrologic network allowing users to discover known features, navigate from those along a route, and find other features on that route. The primary difference being navigation occurs over the hydrologic network rather then a transportation network. This service is now available as part of the dataRetrieval package through the findNLDI function.

Like any routing service, there are three types of information you can provide findNLDI. These include:

(1) a feature to discover (2) direction(s) to travel from that feature (3) types of features to search for along the way

Each of these is discussed below using the following packages:

library(dplyr) # Data frame manipulation
library(ggplot2) # Plotting
library(patchwork) # Arranging plots
library(dataRetrieval) # The star of the show!

What's available?

First, we need to know what features are indexed in the NLDI. The most current offerings can be found using get_nldi_sources, and new features are regularly added. At the time of writing (r Sys.Date()), r nrow(get_nldi_sources()) data sets have been indexed to the NHDPlus and cataloged in the NLDI.

get_nldi_sources()

kable(get_nldi_sources())

Feature/Origin discovery

Features can be requested in two primary ways: Using the the native data set identifier, and using a location. The core feature set of the NLDI include the NHD flowlines, USGS NWIS locations, WQP locations, and HUC12 pour points. Each of these are available as arguments in findNLDI and can be used to request a feature object.

By Identifier

As an illustrative example, NHDPlus features can be requested by their COMID from the NHDPlusV2 data set.

findNLDI(comid = 101)

The returned simple features object contains the native data set identifier ("identifier"), sourceName of the native data set, and the indexed NHD COMID (in this case a duplicate since an NHD feature was requested). In the example above, we see the geometry column is of type LINESTRING. To keep dataRetrieval lightweight, the sf package is not a dependency. Instead, if sf is not installed - or no_sf = TRUE - only the sourceName, comid, and identifier will be returned.

findNLDI(comid = 101, no_sf = TRUE)

To provide another example, we can request the NLDI representation of USGS NWIS gauge 11120000 in both a sf and "non-sf" way. Features indexed to the NHDPlus are returned as POINT objects. If sf is enabled, the sourceName, identifier, X, Y and geometry (sfc) are returned. If sf is not available, the geometry is dropped but the X and Y values are retained.

# local sf installation
findNLDI(nwis = "11120000")

# No sf use/installation
findNLDI(nwis = "11120000", no_sf = TRUE)

Any NLDI feature found with get_nldi_source can be requested by passing a type/ID pair as a list to the origin argument. This will allow the networking capabilities offered in dataRetrieval to grow naturally with the NLDI itself. For example, we can use the origin argument to request features that don't offer a specific parameter.

# Water Data Exchange 2.0 Site CA_45206
findNLDI(origin = list("wade" = "CA_45206"))

Location (longitude/latitude)

If you don't know a feature ID, a longitude/latitude (X,Y) pair or a sf/sfc POINT object can be passed to the location argument. Doing so will identify the NHDPlus catchment the location fall within and will return the associated NHDPlusV2 flowline.

# Request by coordinates
findNLDI(location = c(-115, 40))

# Request by sf/sfc POINT object
ucsb <- sf::st_sfc(sf::st_point(c(-119.8458, Y = 34.4146)), crs = 4326)
findNLDI(location = ucsb)

Navigation

From any feature (comid, huc12, nwis, wqp, origin) or location, four modes of navigation are available and include:

(1) UT: Upper Tributary (2) UM: Upper Mainstream (3) DM: Downstream Tributary (4) DD: Downstream Diversions

A example view of these navigation types can be seen below for NWIS site 11109000.

x <- findNLDI(
  nwis = 11109000,
  nav = c("UT", "DM", "UM", "DD")
)

ggplot() +
  geom_sf(data = x$UT, aes(color = "UT")) +
  geom_sf(data = x$UM, aes(color = "UM")) +
  geom_sf(data = x$DD, aes(color = "DD")) +
  geom_sf(data = x$DM, aes(color = "DM")) +
  geom_sf(data = x$origin, col = "green", size = 2) +
  geom_sf_label(
    data = x$origin,
    aes(label = x$origin$identifier),
    nudge_x = -.08,
    nudge_y = .03
  ) +
  theme_void() +
  scale_color_manual(
    name = "Navigation Direction",
    breaks = c("UT", "UM", "DD", "DM"),
    values = c(
      "UT" = "blue",
      "UM" = "red",
      "DM" = "orange",
      "DD" = "purple"
    )
  ) +
  theme(
    legend.position = c(.05, .7),
    legend.justification = c(0, 0)
  )

One or more modes of navigation can be supplied to the nav argument. For example we can ask to navigate along the upper mainstem (UM) from COMID 101.

summarize.nldi <- function(input) {
  data.frame(
    name = names(input),
    class = sapply(input, class)[1],
    row.names = NULL
  ) %>%
    mutate(feature_count = ifelse(class == "sf", sapply(input, nrow),
      sapply(input, length)
    ))
}


findNLDI(comid = 101, nav = "UM") %>%
  summarize.nldi()

Or along the upper mainstem (UM) and upper tributary (UT) of COMID 101.

findNLDI(comid = 101, nav = c("UM", "UT")) %>%
  summarize.nldi()

In both cases the returned named list includes the origin and the flowlines along the requested navigation. If sf is not enabled, the returned object for a flowpath navigation is a vector of COMIDs.

findNLDI(comid = 101, nav = c("UM", "DM"), no_sf = TRUE) %>%
  summarize.nldi()

Searching along the Navigation

Like the gas station example, any of the features listed in get_nldi_sources can be searched for along the network, for example, we can find all NWIS gauges, on the upper tributary, of COMID 101.

findNLDI(comid = 101, nav = "UT", find = "nwis") %>%
  summarize.nldi()

Of course, more than one resource can be requested, for example, lets replicate the previous search, this time adding Water Quality Points to the returned list:

findNLDI(comid = 101, nav = "UT", find = c("nwis", "wqp")) %>%
  summarize.nldi()

Note that flowlines are no longer the default return for navigation once a new feature is requested. To retain flowlines, the must be explicitly requested.

findNLDI(comid = 101, nav = "UT", find = c("nwis", "flowlines")) %>%
  summarize.nldi()

Upstream Basin Boundary

The Upstream Basin Boundary is a unique object that can be found for any feature by adding "basin" to find. Basins are only geometries with no specific attribute data. Therefore basins can only be returned if sf is installed. Otherwise, the result will be a 0 column data.frame

# with sf
findNLDI(comid = 101, find = "basin") %>%
  summarize.nldi()

# No sf
findNLDI(comid = 101, find = "basin", no_sf = TRUE) %>%
  summarize.nldi()

Distance Constraints

In some cases, particularly for DM and DD navigation, the network can extend for hundreds of kilometers. You can limit (or extend) the distance of your search using the distance_km argument. As the name implies the value provided should be the maximum kilometers you want to search for features. The default for distance_km is 100.

# Default 100 km
findNLDI(comid = 101, nav = "DM", find = c("nwis", "wqp")) %>%
  summarize.nldi()

# Extended 200 km search
findNLDI(comid = 101, nav = "DM", find = c("nwis", "wqp"), distance_km = 200) %>%
  summarize.nldi()

Basic `dataRetrieval` integration

Last, as this functionality is being added to the dataRetrieval package, lets see a basic example of how the NLDI tools provide a discovery mechanism for working with the dataRetrieval tools. Here we will take a location that is near Fountain Creek in Colorado Springs, Colorado.

In this example we will use that location as the origin, navigate upstream along the mainstem, search for NWIS gauges, and use the identified siteIDs to query streamflow records from January 1^st^, 2020 to the current day.

# Upstream nwis, flowlines, and basin
fountainCreek <- findNLDI(
  location = c(-104.780837, 38.786796),
  nav = "UM",
  find = c("nwis", "basin", "flowlines")
)

summarize.nldi(fountainCreek)

# Identify NLDI sites with daily values "dv"
#     and record streamflow ("00060")
#     and recorded flows in 2020
find <- whatNWISdata(sites = gsub(
  "USGS-", "",
  fountainCreek$UM_nwissite$identifier
)) %>%
  filter(
    data_type_cd == "dv",
    parm_cd == "00060",
    end_date > as.Date("2020-01-01")
  ) %>%
  mutate(identifier = paste0("USGS-", site_no)) %>%
  inner_join(fountainCreek$UM_nwissite, by = "identifier") %>%
  sf::st_as_sf()

# Extract Streamflow for identified sites
Q <- readNWISdv(find$site_no,
  parameterCd = "00060",
  startDate = "2020-01-01"
) %>%
  renameNWISColumns()

# Plot!
ggplot() +
  geom_line(
    data = Q,
    aes(x = Date, y = Flow, col = site_no),
    size = .5
  ) +
  facet_wrap(~site_no, nrow = 4) +
  theme_minimal() +
  scale_color_brewer(palette = "Set1") +
  theme(legend.position = "none") +
  ggplot() +
  geom_sf(data = fountainCreek$basin, col = NA) +
  geom_sf(data = fountainCreek$UM_flowlines, col = "blue", alpha = .5) +
  geom_sf(data = find, aes(col = site_no)) +
  scale_color_brewer(palette = "Set1") +
  theme_void() +
  labs(
    title = "2020 Streamflow",
    caption = "Fountain Creek, Colorado"
  ) +
  theme(
    legend.position = "none",
    plot.title = element_text(face = "bold", hjust = .5)
  )