Return nhd plus stream network upstream of a waterbody

extract_network(
  lon = NA,
  lat = NA,
  lines = NA,
  lines_network = TRUE,
  buffer_dist = 0.01,
  maxsteps = 3,
  approve_all_dl = FALSE,
  temporary = TRUE,
  ...
)

Arguments

lon: numeric decimal degree longitude
lat: numeric decimal degree latitude
lines: sf spatial lines object to limit extent of the network search
lines_network: boolean treat lines as the complete network object. If FALSE, simply start network extraction at the terminal reach of the lines object.
buffer_dist: numeric buffer around lat-lon point in dec. deg.
maxsteps: maximum number of stream climbing iterations
approve_all_dl: logical blanket approval to download all missing data. Defaults to TRUE if session is non-interactive.
temporary: logical set FALSE to save data to a persistent rappdirs location
...: parameters passed on to sf::st_read

Value

An sf data frame with LINESTRING geometries

Details

The lon and lat arguments are used for querying the corresponding lake polygon layer which is then used to climb its intersecting stream network.

Examples

if (FALSE) {
library(mapview)
library(sf)

# headwater lakes have no upstream network
coords <- data.frame(lat = 46.32711, lon = -89.58893)
res <- extract_network(coords$lon, coords$lat, maxsteps = 9)

# fails if no lake nhdp lake found within the buffer at the query point
coords <- data.frame(lat = 43.62453, lon = -85.47164)
res <- extract_network(coords$lon, coords$lat, maxsteps = 9)

coords <- data.frame(lat = 20.79722, lon = -156.47833)
# use a non-geographic (projected) buffer size
res <- extract_network(coords$lon, coords$lat, maxsteps = 9,
  buffer_dist = units::as_units(5, "km"))

# use a projected buffer size
res <- extract_network(coords$lon, coords$lat, maxsteps = 9)

# no upstream network for lakes intersecting the Great Lakes
coords <- data.frame(lat = 44.6265, lon = -86.23121)
res <- extract_network(coords$lon, coords$lat, maxsteps = 3)

coords <- data.frame(lat = 42.96523, lon = -89.2527)
res <- extract_network(coords$lon, coords$lat, maxsteps = 9)

mapview(res)
}