Package {hydrogeofetch}

Type:	Package
Title:	Hydrologic Geospatial Fabric Extraction Tool Chain
Version:	2.0.0
Description:	Traverses and works with National Hydrography Dataset Plus (NHDPlus) data. All methods implemented in 'hydrogeofetch' are available in the NHDPlus documentation available from the US Environmental Protection Agency https://www.epa.gov/waterdata/basic-information. Previously published as 'nhdplusTools'.
URL:	https://doi-usgs.github.io/nhdplusTools/ https://github.com/doi-usgs/nhdplusTools/
BugReports:	https://github.com/doi-usgs/nhdplusTools/issues/
Depends:	R (≥ 4.1)
Imports:	hydroloom, dataRetrieval, dplyr, sf, units, jsonlite, httr2, xml2, utils, data.table, methods, arrow, tools, zip, memoise, digest
Suggests:	testthat, knitr, rmarkdown, ggmap, ggplot2, lwgeom, gifski, leaflet, httptest2, StreamCatTools, terra, maptiles, mapsf
License:	CC0
Encoding:	UTF-8
VignetteBuilder:	knitr
Config/testthat/parallel:	true
Config/testthat/edition:	3
LazyData:	true
LazyDataCompression:	xz
Language:	en-US
Config/roxygen2/version:	8.0.0
NeedsCompilation:	no
Packaged:	2026-08-17 13:07:00 UTC; root
Author:	David Blodgett [aut, cre], Mike Johnson [ctb], Marc Weber [ctb], Josh Erickson [ctb], Lauren Koenig [ctb]
Maintainer:	David Blodgett <dblodgett@usgs.gov>
Repository:	CRAN
Date/Publication:	2026-08-20 14:40:08 UTC

Add mainstem identifiers

Description

Joins geoconnex mainstem identifiers onto a table that contains NHDPlusV2 (comid/featureid) or NHDPlusHR (nhdplusid) identifiers using a lookup table.

Usage

add_mainstems(x, join_col = NULL, join_col_type = NULL)

Arguments

Details

Source data comes from the ref_rivers GitHub release here. The csv source data is downloaded once, converted to parquet, and cached in the user data dir indicated by hydrogeofetch_data_dir.

The NHDPlusV2 lookup csv is roughly 120 MB and the NHDPlusHR lookup roughly 50 MB, so the mainstem functions have no examples. The first call for a given type downloads and converts the table; later calls in the same session read a cached parquet file and are fast.

old_dir <- hydrogeofetch_data_dir()
hydrogeofetch_data_dir(file.path(tempdir(), "mainstems"))

add_mainstems(data.frame(comid = c(2804607, 2804621)))

hydrogeofetch_data_dir(old_dir)

Value

x with mainstem_uri and mainstemid columns added.

Align NHD Dataset Names

Description

this function takes any NHDPlus dataset and aligns the attribute names with those used in hydrogeofetch.

Usage

align_nhdplus_names(x)

Arguments

Value

data.frame renamed sf object

Examples

source(system.file("extdata/new_hope_data.R", package = "hydrogeofetch"))

names(new_hope_flowline)

names(new_hope_flowline) <- tolower(names(new_hope_flowline))

new_hope_flowline <- align_nhdplus_names(new_hope_flowline)

names(new_hope_flowline)

Check mainstem identifiers for supersession

Description

Checks whether geoconnex mainstem identifiers have been superseded by a newer reference release. See update_mainstems to resolve a superseded mainstem id to its replacement.

Usage

check_mainstems(x)

Arguments

Value

logical vector the same length as x. TRUE indicates the mainstem id has been superseded.

Examples

check_mainstems(c(2086165, 2086637))

discover geoconnex reference feature layers

Description

Queries the geoconnex.us reference feature server for available layers and attributes.

Usage

discover_geoconnex_reference()

Value

data.frame containing layers available and fields that are available to query.

Examples

  discover_geoconnex_reference()

Discover NHDPlus ID

Description

Multipurpose function to find a COMID of interest.

Note that NHDPlusV2 uses "featureid" for catchment polygons and "comid" for flowline linestrings. These two identifiers are the same where a flowline/catchment pair exists. In some cases, a catchment will not have a flowline and in others, a flowline will not have a catchment.

If a point is provided, and raindrop is false, the comid/featureid integer of the catchment the point is in is returned. This options uses a web service here: https://api.water.usgs.gov/fabric/pygeoapi/collections/catchmentsp

If a point is provided, and raindrop is true, the response is the result of a call to get_raindrop_trace.

If no point is provided, the raindrop argument is ignored and the result is a comid integer derived from a call to get_nldi_feature.

Usage

discover_nhdplus_id(point = NULL, nldi_feature = NULL, raindrop = FALSE)

Arguments

Value

integer COMID or list containing COMID and raindrop trace.

Examples

point <- sf::st_sfc(sf::st_point(c(-76.874, 39.482)), crs = 4326)
discover_nhdplus_id(point)

discover_nhdplus_id(point, raindrop = TRUE)

nldi_nwis <- list(featureSource = "nwissite", featureID = "USGS-08279500")
discover_nhdplus_id(nldi_feature = nldi_nwis)

Download mainstem lookup table from ref_rivers

Description

Downloads and caches a geoconnex mainstem identifier lookup table on your computer.

Usage

download_mainstem_lookup(
  type,
  path = get_mainstem_lookup_path(type),
  force = FALSE
)

Arguments

Details

Source data comes from the ref_rivers GitHub release here. The csv source data is downloaded once, converted to parquet, and cached in the user data dir indicated by hydrogeofetch_data_dir.

The NHDPlusV2 lookup csv is roughly 120 MB and the NHDPlusHR lookup roughly 50 MB, so the mainstem functions have no examples. The first call for a given type downloads and converts the table; later calls in the same session read a cached parquet file and are fast.

old_dir <- hydrogeofetch_data_dir()
hydrogeofetch_data_dir(file.path(tempdir(), "mainstems"))

add_mainstems(data.frame(comid = c(2804607, 2804621)))

hydrogeofetch_data_dir(old_dir)

Value

character path to cached data

Download NHD

Description

Download NHD

Usage

download_nhd(nhd_dir, hu_list, download_files = TRUE)

Arguments

Value

character Paths to geodatabases created.

Examples

hu <- get_huc(sf::st_sfc(sf::st_point(c(-73, 42)), crs = 4326),
                            type = "huc08")

(hu <- substr(hu$huc8, 1, 2))

download_nhd(tempdir(), c(hu, "0203"), download_files = FALSE)

Download NHDPlus HiRes

Description

Download NHDPlus HiRes

Usage

download_nhdplushr(nhd_dir, hu_list, download_files = TRUE, archive = FALSE)

Arguments

Value

character Paths to geodatabases created.

Examples

hu <- get_huc(sf::st_sfc(sf::st_point(c(-73, 42)), crs = 4326),
                            type = "huc08")
if(inherits(hu, "sf")) {
(hu <- substr(hu$huc8, 1, 2))

download_nhdplushr(tempdir(), c(hu, "0203"), download_files = FALSE)

download_nhdplushr(tempdir(), c(hu, "0203"), download_files = FALSE, archive = TRUE)
}

Download seamless National Hydrography Dataset Version 2 (NHDPlusV2)

Description

This function downloads and decompresses staged seamless NHDPlusV2 data. The following requirements are needed: p7zip (MacOS), 7zip (windows) Please see: https://www.epa.gov/waterdata/get-nhdplus-national-hydrography-dataset-plus-data for more information and metadata about this data.

Default downloads lower-48 only. Pass the island archive URL to 'url' to get Hawaii, Puerto Rico, the Virgin Islands, and the Pacific Islands instead. No Alaska data are available.

The lower-48 archive is roughly 8 GB and extraction needs 7zip installed, so this function has no example. Pass any writable directory as outdir; the archive is downloaded there, extracted in place, and the path to the geodatabase returned:

download_nhdplusv2(file.path(tempdir(), "nhdplusv2"))

Usage

download_nhdplusv2(
  outdir,
  url = paste0("https://dmap-data-commons-ow.s3.amazonaws.com/NHDPlusV21/",
    "Data/NationalData/NHDPlusV21_NationalData_Seamless", "_Geodatabase_Lower48_07.7z"),
  progress = TRUE
)

Arguments

Value

character path to the local geodatabase

Download the seamless Reach File (RF1) Database

Description

This function downloads and decompresses staged RF1 data. See: https://water.usgs.gov/GIS/metadata/usgswrd/XML/erf1_2.xml for metadata.

The archive is roughly 46 MB and the server is slow, so this function has no example. Pass any writable directory as outdir; the gzipped e00 is downloaded there, decompressed in place, and the path to the e00 returned:

download_rf1(file.path(tempdir(), "rf1"))

Usage

download_rf1(
  outdir,
  url = "https://water.usgs.gov/GIS/dsdl/erf1_2.e00.gz",
  progress = TRUE
)

Arguments

Value

character path to the local e00 file

Download NHDPlusVAA data from HydroShare

Description

downloads and caches NHDPlusVAA data on your computer

Usage

download_vaa(
  path = get_vaa_path(updated_network),
  force = FALSE,
  updated_network = FALSE
)

Arguments

Details

The VAA data is a aggregate table of information from the NHDPlusV2 elevslope.dbf(s), PlusFlowlineVAA.dbf(s); and NHDFlowlines. All data originates from the EPA NHDPlus Homepage here. To see the location of cached data on your machine use get_vaa_path. To view aggregate data and documentation, see here

Value

character path to cached data

Download the seamless Watershed Boundary Dataset (WBD)

Description

This function downloads and decompresses staged seamless WBD data. Please see: https://prd-tnm.s3.amazonaws.com/StagedProducts/Hydrography/WBD/National/GDB/WBD_National_GDB.xml for metadata.

The national archive is roughly 3 GB, so this function has no example. The "hydrogeofetch Data Access Overview" article works through a download in its Watershed Boundary Dataset section.

Usage

download_wbd(
  outdir,
  url = paste0("https://prd-tnm.s3.amazonaws.com/StagedProducts/",
    "Hydrography/WBD/National/GDB/WBD_National_GDB.zip"),
  progress = TRUE
)

Arguments

Value

character path to the local geodatabase

Get 3DHP Data

Description

Calls the 3DHP_all web service and returns sf data.frames for the selected layers. See https://3dhp.nationalmap.gov/arcgis/rest/services/usgs_3dhp_all/FeatureServer for source data documentation.

Usage

get_3dhp(
  AOI = NULL,
  ids = NULL,
  type = NULL,
  universalreferenceid = NULL,
  t_srs = NULL,
  buffer = 0.5,
  page_size = 2000
)

Arguments

Details

The returned object(s) will have the same Spatial Reference System (SRS) as the input AOI. If a individual or set of IDs are used to query, then the default CRS of EPSG:4269 is preserved. In all cases, a user-defined SRS can be passed to t_srs which will override all previous SRS (either input or default). All buffer and distance operations are handled internally using in EPSG:5070 Albers Equal Area projection

Value

a simple features (sf) object or valid types if no type supplied

Examples

AOI <- sf::st_as_sfc(sf::st_bbox(c(xmin = -89.56684, ymin = 42.99816,
                                   xmax = -89.24681, ymax = 43.17192),
                                 crs = "+proj=longlat +datum=WGS84 +no_defs"))

# get flowlines and hydrolocations
flowlines <- get_3dhp(AOI = AOI, type = "flowline")
hydrolocation <- get_3dhp(AOI = AOI, type = "hydrolocation")
waterbody <- get_3dhp(AOI = AOI, type = "waterbody")

if(!is.null(waterbody) & !is.null(flowlines) & !is.null(hydrolocation)) {
plot(sf::st_geometry(waterbody), col = "lightblue", border = "lightgrey")
plot(sf::st_geometry(flowlines), col = "blue", add = TRUE)
plot(sf::st_geometry(hydrolocation), col = "grey", pch = "+", add = TRUE) }

# given mainstem ids from any source, can query for them in ids.

SU <- get_3dhp(ids = "https://geoconnex.us/ref/mainstems/194408",
               type = "flowline")

if(!is.null(SU))
  plot(sf::st_geometry(SU), col = "blue")

# get all the waterbodies along the Susquehanna river
SU_wb <- get_3dhp(ids = unique(SU$waterbodyid3dhp), type = "waterbody")

if(!is.null(SU_wb)) {
plot(sf::st_geometry(SU_wb[grepl("Otsego", SU_wb$gnisidlabel),]),
     col = "blue", border = "NA") }

# given a workunitid, can query for all features in that work unit with
# ids = "workunitid:300585" -- not run here as a work unit is tens of
# thousands of features.

# given universalreferenceid (reachcodes), can query for them but only
# for hydrolocations. This is useful for looking up mainstem ids.

if(!is.null(hydrolocation)) {
get_3dhp(universalreferenceid =
           head(unique(hydrolocation$universalreferenceid), 5),
         type = "hydrolocation")
}

Index Points to the 3DHP Network

Description

Addresses (snaps) point features to the 3D Hydrography Program (3DHP) network using the HydroAdd3D web service. For each input point the service determines the nearest 3DHP flowline, the snapped location and its elevation, the mainstem of that flowline, the measure along the mainstem, and associated GNIS attributes.

Usage

get_3dhp_index(points, t_srs = NULL, convert_missing = TRUE, batch_size = 100)

Arguments

Details

Points are sent to the service in batches. The 'mainstemid' is a geoconnex mainstem uri (e.g. 'https://geoconnex.us/ref/mainstems/312091') that can be passed to get_3dhp or used with the functions in add_mainstems.

See https://apps.usgs.gov/hydroadd3d for the web service.

Value

sf data.frame of POINT (XYZ) features in 't_srs', one row per input point. The geometry is the snapped location on the 3DHP flowline: the X and Y coordinates are the snapped position and the Z coordinate is the snapped elevation in meters. Points that did not snap carry their original input location as the geometry. Attribute columns:

• 'source_id' - input row, echoed for joins

• 'zsnap' - snapped elevation in meters (the geometry Z coordinate)

• 'mainstemid' - geoconnex mainstem uri

• 'm' - measure along the mainstem, 0 = downstream, 100 = upstream

• 'gnisidlabel' / 'gnisid' - GNIS name and id

• 'featuredate' - date the 3DHP feature was loaded

• 'snapdistance' - distance in meters from the original to the snapped point

• 'snapdate' - timestamp of the addressing run

• 'snapped' - logical, added when 'convert_missing = TRUE', flagging points that addressed successfully

Examples

points <- sf::st_sf(
  id = c(1, 2, 3),
  geometry = sf::st_sfc(sf::st_point(c(-76.86934, 39.49328)),
                        sf::st_point(c(-76.91711, 39.40884)),
                        sf::st_point(c(-76.88081, 39.36354)),
                        crs = 4326))

indexed <- get_3dhp_index(points)

if(inherits(indexed, "sf")) {

# points that did not snap to a flowline
indexed[!indexed$snapped, ]

# pull the flowline a point was addressed to
get_3dhp(ids = indexed$mainstemid[indexed$snapped][1], type = "flowline")

}

Return RPU or VPU boundaries

Description

Return RPU or VPU boundaries

Usage

get_boundaries(type = "vpu")

Arguments

Value

An object of class "sf"

Get Catchment Characteristics

Description

Downloads (subsets of) catchment characteristics from a cloud data store. See get_characteristics_metadata for available characteristics.

When source = "usgs" (the default), data is retrieved from: Wieczorek, M.E., Jackson, S.E., and Schwarz, G.E., 2018, Select Attributes for NHDPlus Version 2.1 Reach Catchments and Modified Network Routed Upstream Watersheds for the Conterminous United States (ver. 3.0, January 2021): U.S. Geological Survey data release, doi:10.5066/F7765D7V.

When source = "streamcat", data is retrieved from the EPA StreamCat dataset via the StreamCatTools package (must be installed separately). The aoi parameter controls the area of interest for StreamCat queries.

Usage

get_catchment_characteristics(
  varname,
  ids,
  reference_fabric = "nhdplusv2",
  source = "usgs",
  aoi = "cat"
)

Arguments

Value

data.frame with columns characteristic_id, comid, characteristic_value, and percent_nodata. NULL if no requested variables were found or the data store is unavailable.

Examples

  get_catchment_characteristics("CAT_BFI", c(5329343, 5329427))

Get catchment characteristics metadata table

Description

Download and cache table of catchment characteristics.

When source = "usgs" (the default), returns metadata from: Wieczorek, M.E., Jackson, S.E., and Schwarz, G.E., 2018, Select Attributes for NHDPlus Version 2.1 Reach Catchments and Modified Network Routed Upstream Watersheds for the Conterminous United States (ver. 3.0, January 2021): U.S. Geological Survey data release, doi:10.5066/F7765D7V.

When source = "streamcat", returns metric metadata from the EPA StreamCat dataset accessed via the StreamCatTools package (must be installed separately): Weber, Marc H, Hill, Ryan A., Brookes, Allen F. 2024, StreamCatTools: Tools to work with the StreamCat API within R and access the full suite of StreamCat and LakeCat metrics, https://usepa.github.io/StreamCatTools

Usage

get_characteristics_metadata(search, source = "usgs", cache = TRUE)

Arguments

Value

data.frame of characteristic metadata with columns ID, description, units, datasetLabel, datasetURL, themeLabel, themeURL, watershedType, sbid, end, s3_url, and http_url. NULL if the metadata service is unavailable.

Examples

# the metadata table is cached in hydrogeofetch_data_dir(); point it at a
# temporary directory so this example does not write to user space.
old_dir <- hydrogeofetch_data_dir()
hydrogeofetch_data_dir(file.path(tempdir(check = TRUE), "meta_demo"))

get_characteristics_metadata()

hydrogeofetch_data_dir(old_dir)

Estimate drainage area from HUC and catchment data

Description

Combines HUC12 areas upstream of HUC outlets with NHDPlusV2 catchment areas for the portion of the basin between the outlet and HUC12 outlets to produce a drainage area estimate. Non-contributing areas captured in HUC12 boundaries are included.

Usage

get_drainage_area_estimates(
  start,
  catchments = FALSE,
  nhdplushr = TRUE,
  local_navigation = FALSE,
  huc12_data = NULL,
  huc12_outlets = NULL,
  waterbody_data = NULL,
  catchment_data = NULL,
  HU_inclusion_override = NULL,
  outlet_split_threshold_m = 100
)

Arguments

Details

By default, network navigation is performed via the NLDI web service and flowline attributes are retrieved from the NHDPlusV2 OGC API. When local_navigation = TRUE, the NHDPlusV2 Value Added Attributes (get_vaa) are used for network navigation instead, and only HUC12 pour points are fetched from the NLDI.

HUC drainage area is used upstream of the nearest HUC12 outlet. Between the outlet (e.g. gage) and the HUC outlet(s), catchment areas are used. For large upstream areas the largest HUC level is used to define connectedness, but drainage estimates are derived from HUC12 because that is where the non-contributing area attribute lives.

Three pairs of drainage area estimates are returned: one using only NLDI-identified HUC12s (HUC12-level), one using HUC10-level queries, and one using HUC08-level queries for basins spanning multiple HUC08s. Each pair includes a total and a contributing-only estimate derived from the ncontrb_a (non-contributing acres) attribute on HUC12 features.

Value

list with elements:

da_huc12_sqkm

numeric. Total DA using NLDI-identified HUC12s only.

da_huc10_sqkm

numeric or NA. Total DA using HUC10-level queries. NA when basin is within a single HUC10.

da_huc08_sqkm

numeric or NA. Total DA using HUC08-level queries. NA when basin is within a single HUC08.

contrib_da_huc12_sqkm

numeric. Contributing DA (HUC12-only).

contrib_da_huc10_sqkm

numeric or NA. Contributing DA (HUC10-level).

contrib_da_huc08_sqkm

numeric or NA. Contributing DA (HUC08-level).

network_da_sqkm

numeric. Network-derived total DA for comparison.

basin_da_sqkm

numeric or NA. Area of drainage_basin. Not the same quantity as network_da_sqkm: that one is the NHDPlusV2 totdasqkm attribute for the outlet flowline and covers the whole outlet catchment, while this is a geometric area with the portion below the submitted point removed. The two differ by that sliver plus the usual drift between attribute and geometry. NA when drainage_basin is NULL.

nhdplushr_network_dasqkm

numeric or NA. Drainage area from NHDPlusHR catchments upstream of the matched HR flowline. NA with a warning when the HR web service is unavailable or fails.

nhdplushr_boundary

sfc_GEOMETRY or NULL. Dissolved boundary of upstream NHDPlusHR catchments. NULL when the HR estimate is unavailable.

start_feature

sf data.frame. The resolved NLDI start feature.

hu12_by_huc12

sf data.frame. NLDI-identified upstream HUC12 polygons (EPSG:5070).

hu12_by_huc10

sf data.frame or NULL. Upstream HUC12 polygons (HUC10 query). NULL when basin is within a single HUC10.

hu12_by_huc08

sf data.frame or NULL. Upstream HUC12 polygons (HUC08 query). NULL when basin is within a single HUC08.

extra_catchments

sf data.frame. Catchments between outlet and HUC12 outlets.

split_catchment

sf data.frame. Split catchment at HUC12 outlet(s).

drainage_basin

sf data.frame or NULL. NHDPlusV2 catchments for the upstream network dissolved into a single boundary (EPSG:5070, with a dasqkm column), with the outlet catchment replaced by the portion above the submitted point. No HUC boundaries are involved. Available when catchments = TRUE or catchment_data was supplied, and always for headwater starts, where the upstream network is a single catchment; NULL otherwise. Also NULL when a split of the outlet catchment was due but the service returned none, since the basin would otherwise carry area below the submitted point.

all_network

data.frame. Full upstream flowline attributes.

all_catchments

sf data.frame or NULL. NHDPlusV2 catchment polygons for the full upstream network. NULL when catchments = FALSE.

outlet_flowline_measure

numeric or NULL. Flowline measure (0–100) for the gage on its outlet flowline. NULL for non-gage starts.

outlet_split_catchment

sf data.frame or NULL. Split catchment at the gage point. Contains "catchment" and "splitCatchment" rows with areas. NULL when no split is needed (gage at outlet or below threshold).

hu12_outlet

sf data.frame. HUC12 pour points upstream of outlet.

Examples

# Black Earth Creek
start <- list(featureSource = "nwissite", featureID = "USGS-05406500")
result <- get_drainage_area_estimates(start)
result$da_huc10_sqkm
result$network_da_sqkm

Get Elevation Along Path (experimental)

Description

Uses a cross section retrieval web services to retrieve elevation along a path.

Usage

get_elev_along_path(points, num_pts, res = 1, status = TRUE)

Arguments

Value

sf data.frame containing points retrieved. Names include "id", "distance_m", "elevation_m", "spatial_ref", "geometry", and ".group". .group tracks which input point each set of output points belongs to.

Examples

point1 <- sf::st_sfc(sf::st_point(x = c(-105.9667, 36.17602)), crs = 4326)
point2 <- sf::st_sfc(sf::st_point(x = c(-105.97768, 36.17526)), crs = 4326)
point3 <- sf::st_sfc(sf::st_point(x = c(-105.98869, 36.17450)), crs = 4326)

points <- sf::st_as_sf(c(point1, point2, point3))

(xs <- get_elev_along_path(points, 100))

if(inherits(xs, "sf")) {

bbox <- sf::st_bbox(xs) + c(-0.005, -0.005, 0.005, 0.005)

hydrogeofetch::plot_nhdplus(bbox = bbox, cache_data = FALSE)

plot(sf::st_transform(sf::st_geometry(xs), 3857), pch = ".", add = TRUE, col = "red")
plot(sf::st_transform(sf::st_sfc(point1, crs = 4326), 3857), add = TRUE)
plot(sf::st_transform(sf::st_sfc(point2, crs = 4326), 3857), add = TRUE)
plot(sf::st_transform(sf::st_sfc(point3, crs = 4326), 3857), add = TRUE)

plot(xs$distance_m, xs$elevation_m)
}

Find gagesII Features

Description

Subsets the gagesII dataset by location (POINT), area (POLYGON), or set of IDs. See <doi:10.5066/P96CPHOT> for documentation of source data.

Usage

get_gagesII(AOI = NULL, id = NULL, t_srs = NULL, buffer = 0.5, basin = FALSE)

Arguments

Details

The returned object(s) will have the same Spatial Reference System (SRS) as the input AOI. If a individual or set of IDs are used to query, then the default server CRS of EPSG:4326 is preserved. In all cases, a user-defined SRS can be passed to t_srs which will override all previous SRS (either input or default). All buffer and distance operations are handled internally using an EPSG:5070 Albers Equal Area projection

Value

a simple features (sf) object, or a data.frame when skip_geometry = TRUE

get geoconnex reference feature layers

Description

Queries the geoconnex reference feature server for features of interest.

Usage

get_geoconnex_reference(
  AOI,
  type = NULL,
  t_srs = NULL,
  buffer = 0.5,
  status = TRUE
)

Arguments

Value

sf data.frame containing requested reference features

Examples

dplyr::distinct(discover_geoconnex_reference()[c("id", "title")])

AOI <- sf::st_as_sfc(sf::st_bbox(c(xmin = -89.56684, ymin = 42.99816,
                                   xmax = -89.24681, ymax = 43.17192),
                                 crs = "+proj=longlat +datum=WGS84 +no_defs"))

get_geoconnex_reference(AOI, type = "hu04")

get_geoconnex_reference("https://geoconnex.us/ref/mainstems/315626", type = "hu04", )

AOI <- sf::st_sfc(sf::st_point(c(-89.56684, 42.99816)),
                  crs = "+proj=longlat +datum=WGS84 +no_defs")

get_geoconnex_reference(AOI, type = "hu04", buffer = 100000)

Get NHDPlus HiRes Data

Description

Use to remove unwanted detail NHDPlusHR data See get_nhdplushr for examples.

Usage

get_hr_data(
  gdb,
  layer = NULL,
  min_size_sqkm = NULL,
  simp = NULL,
  proj = NULL,
  rename = TRUE
)

Arguments

Value

sf data.frame containing requested data

Find WBD HUC unit subsets

Description

Subsets WBD features by location (POINT), area (POLYGON), or set of HUC IDs.

Usage

get_huc(AOI = NULL, id = NULL, t_srs = NULL, buffer = 0.5, type = NULL)

Arguments

Details

The returned object(s) will have the same Spatial Reference System (SRS) as the input AOI. If a individual or set of IDs are used to query, then the default server CRS of EPSG:4326 is preserved. In all cases, a user-defined SRS can be passed to t_srs which will override all previous SRS (either input or default). All buffer and distance operations are handled internally using an EPSG:5070 Albers Equal Area projection

Value

a simple features (sf) object, or a data.frame when skip_geometry = TRUE

File path to mainstem lookup table cache

Description

hydrogeofetch will download and cache a parquet file translating geoconnex mainstem identifiers to NHDPlusV2 comids (type = "nhdpv2") or NHDPlusHR NHDPlusIDs (type = "nhdphr"). Will use the user data dir indicated by hydrogeofetch_data_dir.

Usage

get_mainstem_lookup_path(type)

Arguments

Details

Source data comes from the ref_rivers GitHub release here. The csv source data is downloaded once, converted to parquet, and cached in the user data dir indicated by hydrogeofetch_data_dir.

The NHDPlusV2 lookup csv is roughly 120 MB and the NHDPlusHR lookup roughly 50 MB, so the mainstem functions have no examples. The first call for a given type downloads and converts the table; later calls in the same session read a cached parquet file and are fast.

old_dir <- hydrogeofetch_data_dir()
hydrogeofetch_data_dir(file.path(tempdir(), "mainstems"))

add_mainstems(data.frame(comid = c(2804607, 2804621)))

hydrogeofetch_data_dir(old_dir)

Value

character file path

Find NHDPlusV2 Areas

Description

Subsets NHDPlusV2 Area features by location (POINT), area (POLYGON), or set of IDs. See download_nhdplusv2 for source data documentation.

Usage

get_nhdarea(AOI = NULL, id = NULL, t_srs = NULL, buffer = 0.5)

Arguments

Details

The returned object(s) will have the same Spatial Reference System (SRS) as the input AOI. If a individual or set of IDs are used to query, then the default server CRS of EPSG:4326 is preserved. In all cases, a user-defined SRS can be passed to t_srs which will override all previous SRS (either input or default). All buffer and distance operations are handled internally using an EPSG:5070 Albers Equal Area projection

Value

a simple features (sf) object, or a data.frame when skip_geometry = TRUE

Get NHDPlusHR Data

Description

Calls the NHDPlus_HR web service and returns sf data.frames for the selected layers. See https://hydro.nationalmap.gov/arcgis/rest/services/NHDPlus_HR/MapServer for source data documentation.

Usage

get_nhdphr(
  AOI = NULL,
  ids = NULL,
  type = NULL,
  reachcode = NULL,
  t_srs = NULL,
  buffer = 0.5,
  page_size = 2000
)

Arguments

Details

The returned object(s) will have the same Spatial Reference System (SRS) as the input AOI. If a individual or set of IDs are used to query, then the default CRS of EPSG:4269 is preserved. In all cases, a user-defined SRS can be passed to t_srs which will override all previous SRS (either input or default). All buffer and distance operations are handled internally using in EPSG:5070 Albers Equal Area projection

Value

a simple features (sf) object or valid types if no type supplied

Examples

AOI <- sf::st_as_sfc(sf::st_bbox(c(xmin = -89.56684, ymin = 42.99816,
                                   xmax = -89.24681, ymax = 43.17192),
                                 crs = "+proj=longlat +datum=WGS84 +no_defs"))

# get flowlines and hydrolocations
flowlines <- get_nhdphr(AOI = AOI, type = "networknhdflowline")
point <- get_nhdphr(AOI = AOI, type = "nhdpoint")
waterbody <- get_nhdphr(AOI = AOI, type = "nhdwaterbody")

if(!is.null(waterbody) & !is.null(flowlines) & !is.null(point)) {
plot(sf::st_geometry(waterbody), col = "lightblue", border = "lightgrey")
plot(sf::st_geometry(flowlines), col = "blue", add = TRUE)
plot(sf::st_geometry(point), col = "grey", pch = "+", add = TRUE) }

# given universalreferenceid (reachcodes), can query for them but only
# for hydrolocations. This is useful for looking up mainstem ids.

get_nhdphr(reachcode = "13020101021927", type = "networknhdflowline")

Get National Hydrography Dataset V2 Subsets (Multirealization)

Description

Subsets NHDPlusV2 features by location (POINT), area (POLYGON), or set of COMIDs. Multi realizations are supported allowing you to query for flowlines, catchments, or outlets.

NOTE: not all flowlines have catchments and not all catchment have flowlines. Flowlines without catchments will have a 0 areasqkm attribute. Catchments without flowlines terminate in a sink and will have a negative valued comid.

Usage

get_nhdplus(
  AOI = NULL,
  comid = NULL,
  nwis = NULL,
  realization = "flowline",
  streamorder = NULL,
  t_srs = NULL,
  properties = NULL,
  skip_geometry = FALSE
)

Arguments

Details

The returned object(s) will have the same Spatial Reference System (SRS) as the input AOI. If a individual or set of IDs are used to query, then the default server CRS of EPSG:4326 is preserved. In all cases, a user-defined SRS can be passed to t_srs which will override all previous SRS (either input or default). All buffer and distance operations are handled internally using an EPSG:5070 Albers Equal Area projection

Value

sf single, or named list of, simple feature objects. NULL if no realization succeeded (e.g. degraded service).

Examples

 point <- sf::st_sfc(sf::st_point(c(-119.845, 34.4146)), crs = 4326)
 get_nhdplus(point)
 get_nhdplus(point, realization = "catchment")
 get_nhdplus(point, realization = "all")
 get_nhdplus(comid = 101)
 get_nhdplus(nwis  = c(11120000, 11120500))
 area <- sf::st_as_sfc(sf::st_bbox(c(xmin = -119.8851, xmax =-119.8361,
 ymax = 34.42439, ymin = 34.40473), crs = 4326))
 get_nhdplus(area)
 get_nhdplus(area, realization = "flowline", streamorder = 3)
 

Get NHDPlus HiRes

Description

Get NHDPlus HiRes

Usage

get_nhdplushr(
  hr_dir,
  out_gpkg = NULL,
  layers = c("NHDFlowline", "NHDPlusCatchment"),
  pattern = ".*GDB.gdb$",
  check_terminals = TRUE,
  overwrite = FALSE,
  keep_cols = NULL,
  ...
)

Arguments

Details

NHDFlowline is joined to value added attributes prior to being returned. Names are not modified from the NHDPlusHR geodatabase. Set layers to "NULL" to get all layers.

Value

sf data.frames containing output that may also be written to a geopackage for later use.

Examples

# a small NHDPlusHR sample, downloaded to a temporary directory
work_dir <- file.path(tempdir(check = TRUE), "hr_example")
try(source(system.file("extdata/nhdplushr_data.R", package = "hydrogeofetch")))

if(file.exists(file.path(work_dir, "03_sub.gpkg"))) {

get_nhdplushr(work_dir, pattern = "03_sub.gpkg")

get_nhdplushr(work_dir, pattern = "03_sub.gpkg", layers = NULL)

get_nhdplushr(work_dir, pattern = "03_sub.gpkg",
              layers = "NHDFlowline",
              min_size_sqkm = 10, simp = 10, proj = 5070)

}

# For real data, download subregions first -- see download_nhdplushr().

Get NLDI Basin Boundary

Description

Get a basin boundary for a given NLDI feature.

Usage

get_nldi_basin(nldi_feature, simplify = TRUE, split = FALSE)

Arguments

Details

Only resolves to the nearest NHDPlus catchment divide. See: https://waterdata.usgs.gov/blog/nldi-intro/ for more info on the nldi.

Value

sf data.frame with result basin boundary

Examples

library(sf)
library(dplyr)

nldi_nwis <- list(featureSource = "nwissite", featureID = "USGS-05428500")

site <- get_nldi_feature(nldi_nwis)

basin <- get_nldi_basin(nldi_feature = nldi_nwis)

plot(st_geometry(basin))

basin

basin2 <- get_nldi_basin(nldi_feature = nldi_nwis,
                         simplify = FALSE, split = TRUE)

if(inherits(basin, "sf") & inherits(basin2, "sf")) {

length(st_coordinates(basin))
length(st_coordinates(basin2))

plot(st_geometry(st_buffer(st_transform(site, 5070),
                          units::set_units(3000, "m"))), border = NA)

plot(st_geometry(st_transform(site, 5070)), add = TRUE)
plot(st_geometry(st_transform(basin2, 5070)), add = TRUE)

plot(st_geometry(st_transform(basin, 5070)), border = "red", add = TRUE)

}

Get NLDI Feature

Description

Get a single feature from the NLDI

Usage

get_nldi_feature(nldi_feature)

Arguments

Value

sf data.frame with one feature

Examples

get_nldi_feature(list("featureSource" = "nwissite", featureID = "USGS-05428500"))

Get NLDI Index

Description

uses the Network Linked Data Index to retrieve an estimated network location for the given point. If not within a grid cell of a flowline, will use a raindrop trace service to find the nearest downslope flowline location.

Usage

get_nldi_index(location)

Arguments

Value

sf data.frame of estimated network locations for the given point. NULL if the NLDI service is unavailable.

Examples

index <- get_nldi_index(c(-89.276, 42.988))

if(inherits(index, "sf")) {

plot_nhdplus(
  bbox = sf::st_bbox(
    sf::st_buffer(
      sf::st_transform(index[1,], 5070), units::set_units(1000, "m")
      )
    )
)
plot(sf::st_geometry(sf::st_transform(index, 3857)), add = TRUE)

}

Discover USGS NWIS Stream Gages

Description

Returns a POINT feature class of active, stream network, NWIS gages for an Area of Interest. If a POINT feature is used as an AOI, then the returned sites within the requested buffer, are sorted by distance (in meters) from that POINT.

Usage

get_nwis(AOI = NULL, t_srs = NULL, buffer = 20000)

Arguments

Details

The returned object(s) will have the same Spatial Reference System (SRS) as the input AOI. If a individual or set of IDs are used to query, then the default server CRS of EPSG:4326 is preserved. In all cases, a user-defined SRS can be passed to t_srs which will override all previous SRS (either input or default). All buffer and distance operations are handled internally using an EPSG:5070 Albers Equal Area projection

Value

a simple features (sf) object, or a data.frame when skip_geometry = TRUE

Get Raindrop Trace

Description

Uses a raindrop trace web service to trace the nhdplus digital elevation model to the nearest downslope flowline.

Usage

get_raindrop_trace(point, direction = "down")

Arguments

Value

sf data.frame containing raindrop trace and requested portion of flowline.

Examples

point <- sf::st_sfc(sf::st_point(x = c(-89.2158, 42.9561)), crs = 4326)

(trace <- get_raindrop_trace(point))

if(inherits(trace, "sf")) {
bbox <- sf::st_bbox(trace) + c(-0.005, -0.005, 0.005, 0.005)

hydrogeofetch::plot_nhdplus(bbox = bbox, cache_data = FALSE)

plot(sf::st_transform(sf::st_sfc(point, crs = 4326), 3857), add = TRUE)
plot(sf::st_transform(sf::st_geometry(trace)[1], 3857), add = TRUE, col = "red")
plot(sf::st_transform(sf::st_geometry(trace)[2], 3857), add = TRUE, col = "black")
}

Get split catchment

Description

Uses a catchment splitting web service to retrieve the portion of a catchment upstream of the point provided.

Usage

get_split_catchment(point, upstream = TRUE)

Arguments

Details

This service works within the conterminous US NHDPlusV2 domain. If the point provided falls on an NHDPlusV2 flowline as retrieved from get_raindrop_trace the catchment will be split across the flow line. IF the point is not along the flowline a small sub catchment will typically result. As a result, most users of this function will want to use get_raindrop_trace prior to calls to this function.

An attempt is made to eliminate polygon shards if they exist in the output. However, there is a chance that this function will return a multipolygon data.frame.

Value

sf data.frame containing the local catchment, the split portion and optionally the total drainage basin.

Examples

point <- sf::st_sfc(sf::st_point(x = c(-89.2158, 42.9561)), crs = 4326)

trace <- get_raindrop_trace(point)

if(inherits(trace, "sf")) {

(snap_point <- sf::st_sfc(sf::st_point(trace$intersection_point[[1]]),
                          crs = 4326))

(catchment <- get_split_catchment(snap_point))

if(inherits(catchment, "sf")) {
bbox <- sf::st_bbox(catchment) + c(-0.005, -0.005, 0.005, 0.005)

hydrogeofetch::plot_nhdplus(bbox = bbox, cache_data = FALSE)

if(length(sf::st_geometry(catchment)) > 1) {
  plot(sf::st_transform(sf::st_geometry(catchment)[2], 3857), add = TRUE, col = "black")
}
plot(sf::st_transform(sf::st_geometry(catchment)[1], 3857), add = TRUE, col = "red")
plot(sf::st_transform(sf::st_sfc(point, crs = 4326), 3857), add = TRUE, col = "white")
}

(catchment <- get_split_catchment(snap_point, upstream = FALSE))


if(inherits(catchment, "sf")) {
bbox <- sf::st_bbox(catchment) + c(-0.005, -0.005, 0.005, 0.005)

hydrogeofetch::plot_nhdplus(bbox = bbox, cache_data = FALSE)

plot(sf::st_transform(sf::st_geometry(catchment)[1], 3857), add = TRUE, col = "red")
if(length(sf::st_geometry(catchment)) > 1) {
  plot(sf::st_transform(sf::st_geometry(catchment)[2], 3857), add = TRUE, col = "black")
}
plot(sf::st_transform(sf::st_sfc(point, crs = 4326), 3857), add = TRUE, col = "white")
}

pour_point <- sf::st_sfc(sf::st_point(x = c(-89.25619, 42.98646)), crs = 4326)

(catchment <- get_split_catchment(pour_point, upstream = FALSE))

if(inherits(catchment, "sf")) {
bbox <- sf::st_bbox(catchment) + c(-0.005, -0.005, 0.005, 0.005)

hydrogeofetch::plot_nhdplus(bbox = bbox, cache_data = FALSE)

plot(sf::st_transform(sf::st_geometry(catchment)[1], 3857), add = TRUE, col = "red")
if(length(sf::st_geometry(catchment)) > 1) {
  plot(sf::st_transform(sf::st_geometry(catchment)[2], 3857), add = TRUE, col = "black")
}
plot(sf::st_transform(sf::st_sfc(pour_point, crs = 4326), 3857), add = TRUE, col = "white")
}

}

Get tocomid

Description

Given flowlines with fromnode and tonode attributes, will return a toid attribute that is the result of joining tonode and fromnode attributes. In the case that a terminalpa attribute is included, the join is executed by terminalpa group. This is done grouped by terminalpathID because duplicate node ids have been encountered across basins in some datasets. If 'remove_coastal' is 'TRUE' (the default) either ftype or fcode are required. Uses the add_toids function.

Usage

get_tocomid(
  x,
  return_dendritic = TRUE,
  missing = 0,
  remove_coastal = TRUE,
  add = TRUE
)

Arguments

Value

data.frame containing comid and tocomid attributes or all attributes provided with comid and tocomid in the first and second columns..

Examples

source(system.file("extdata", "sample_flines.R", package = "hydrogeofetch"))

tocomid <- get_tocomid(sample_flines)

tocomid <- get_tocomid(sample_flines, return_dendritic = FALSE)

NHDPlusV2 Attribute Subset

Description

Return requested NHDPlusV2 Attributes.

Calling this function downloads the attribute table to get_vaa_path if it is not already cached there. The base table is roughly 260 MB and the updated network table another 170 MB, so this function has no example. To try it without writing to your user data directory, point hydrogeofetch_data_dir at a temporary path first:

old_dir <- hydrogeofetch_data_dir()
hydrogeofetch_data_dir(file.path(tempdir(), "vaa"))

get_vaa("slope")
get_vaa(c("slope", "lengthkm"))
get_vaa("reachcode", updated_network = TRUE)

hydrogeofetch_data_dir(old_dir)

Usage

get_vaa(
  atts = NULL,
  path = get_vaa_path(),
  download = TRUE,
  updated_network = FALSE
)

Arguments

Details

The VAA data is a aggregate table of information from the NHDPlusV2 elevslope.dbf(s), PlusFlowlineVAA.dbf(s); and NHDFlowlines. All data originates from the EPA NHDPlus Homepage here. To see the location of cached data on your machine use get_vaa_path. To view aggregate data and documentation, see here

Value

data.frame containing requested VAA data

Available NHDPlusV2 Attributes

Description

Find variables available from the NHDPlusV2 attribute data.frame

Calling this function downloads the attribute table to get_vaa_path if it is not already cached there. The table is roughly 260 MB, so this function has no example.

old_dir <- hydrogeofetch_data_dir()
hydrogeofetch_data_dir(file.path(tempdir(), "vaa"))

get_vaa_names()

hydrogeofetch_data_dir(old_dir)

Usage

get_vaa_names(updated_network = FALSE)

Arguments

Details

The VAA data is a aggregate table of information from the NHDPlusV2 elevslope.dbf(s), PlusFlowlineVAA.dbf(s); and NHDFlowlines. All data originates from the EPA NHDPlus Homepage here. To see the location of cached data on your machine use get_vaa_path. To view aggregate data and documentation, see here

Value

character vector

File path to value added attribute (vaa) Cache

Description

hydrogeofetch will download and cache a parquet file with NHDPlusV2 attribute data sans geometry. This function returns the file path to the cached file. Will use the user data dir indicated by hydrogeofetch_data_dir.

Usage

get_vaa_path(updated_network = FALSE)

Arguments

Details

The VAA data is a aggregate table of information from the NHDPlusV2 elevslope.dbf(s), PlusFlowlineVAA.dbf(s); and NHDFlowlines. All data originates from the EPA NHDPlus Homepage here. To see the location of cached data on your machine use get_vaa_path. To view aggregate data and documentation, see here

Value

character file path

Examples

get_vaa_path()

get_vaa_path(updated_network = TRUE)

Find NHDPlusV2 Water Bodies

Description

Subsets NHDPlusV2 waterbody features by location (POINT), area (POLYGON), or set of IDs. See download_nhdplusv2 for source data documentation.

Usage

get_waterbodies(AOI = NULL, id = NULL, t_srs = NULL, buffer = 0.5)

Arguments

Details

The returned object(s) will have the same Spatial Reference System (SRS) as the input AOI. If a individual or set of IDs are used to query, then the default server CRS of EPSG:4326 is preserved. In all cases, a user-defined SRS can be passed to t_srs which will override all previous SRS (either input or default). All buffer and distance operations are handled internally using an EPSG:5070 Albers Equal Area projection

Value

a simple features (sf) object, or a data.frame when skip_geometry = TRUE

Get Waterbody Outlet

Description

Get Waterbody Outlet

Usage

get_wb_outlet(lake_id, network)

Arguments

Value

sf data.frame with single record of network COMID associated with most-downstream reach in the NHD Waterbody

Examples

source(system.file("extdata/sample_data.R", package = "hydrogeofetch"))

fline <- sf::read_sf(sample_data, "NHDFlowline_Network")
wtbdy <- sf::read_sf(sample_data, "NHDWaterbody")

lake_COMID <- wtbdy$COMID[wtbdy$GNIS_NAME=='Lake Mendota 254']

get_wb_outlet(13293262, fline)

Get Cross Section From Point (experimental)

Description

Uses a cross section retrieval web services to retrieve a cross section given a point and specified width. Orientation is determined based on direction of the flowline found near point. This function uses a 10m National Elevation Dataset request on the back end.

Usage

get_xs_point(point, width, num_pts)

Arguments

Value

sf data.frame containing points retrieved.

Examples

point <- sf::st_sfc(sf::st_point(x = c(-105.97218, 36.17592)), crs = 4326)

(xs <- get_xs_point(point, 300, 100))

if(inherits(xs, "sf")) {

bbox <- sf::st_bbox(xs) + c(-0.005, -0.005, 0.005, 0.005)

hydrogeofetch::plot_nhdplus(bbox = bbox, cache_data = FALSE)

plot(sf::st_transform(sf::st_geometry(xs), 3857), pch = ".", add = TRUE, col = "red")
plot(sf::st_transform(sf::st_sfc(point, crs = 4326), 3857), add = TRUE)

plot(xs$distance_m, xs$elevation_m)
}

Get Cross Section Endpoints (experimental)

Description

Uses a cross section retrieval web services to retrieve a cross section between two endpoints.

Usage

get_xs_points(point1, point2, num_pts, res = 1)

Arguments

Value

sf data.frame containing points retrieved.

Examples

point1 <- sf::st_sfc(sf::st_point(x = c(-105.9667, 36.17602)), crs = 4326)
point2 <- sf::st_sfc(sf::st_point(x = c(-105.97768, 36.17526)), crs = 4326)

(xs <- get_xs_points(point1, point2, 100))

if(inherits(xs, "sf")) {

bbox <- sf::st_bbox(xs) + c(-0.005, -0.005, 0.005, 0.005)

hydrogeofetch::plot_nhdplus(bbox = bbox, cache_data = FALSE)

plot(sf::st_transform(sf::st_geometry(xs), 3857), pch = ".", add = TRUE, col = "red")
plot(sf::st_transform(sf::st_sfc(point1, crs = 4326), 3857), add = TRUE)
plot(sf::st_transform(sf::st_sfc(point2, crs = 4326), 3857), add = TRUE)

plot(xs$distance_m, xs$elevation_m)
}

Clear cached hydrogeofetch data

Description

Removes files that hydrogeofetch has written to hydrogeofetch_data_dir. Use this to reclaim space or to force downloaded data to be refreshed.

Downloaded datasets are kept until removed with this function – nothing expires them automatically – so it is worth calling periodically if you use get_vaa or add_mainstems on many different areas.

Usage

hydrogeofetch_cache_clear(what = "all", ask = interactive())

Arguments

Value

character vector of paths removed, invisibly.

Examples

# pointed at an empty temporary directory here so the example does not
# remove anything you have actually cached
old_dir <- hydrogeofetch_data_dir()
hydrogeofetch_data_dir(file.path(tempdir(), "hgf_cache_demo"))

hydrogeofetch_cache_clear(ask = FALSE)

hydrogeofetch_data_dir(old_dir)

hydrogeofetch cache settings

Description

Provides an interface to adjust hydrogeofetch 'memoise' cache.

Mode and timeout can also be set using environment variables. 'HYDROGEOFETCH_MEMOISE_CACHE' and 'HYDROGEOFETCH_MEMOISE_TIMEOUT' are used unless overridden with this function. The old 'NHDPLUSTOOLS_*' names are still recognized as a fallback.

The default mode is 'memory', so cached responses live only as long as the R session. Set 'filesystem' to persist them in hydrogeofetch_data_dir, and see hydrogeofetch_cache_clear to remove them again.

Usage

hydrogeofetch_cache_settings(mode = NULL, timeout = NULL)

Arguments

Value

list containing settings at time of calling. If inputs are NULL, current settings. If settings are altered, previous setting values.

Get or set hydrogeofetch data directory

Description

if left unset, will return the user data dir as returned by 'tools::R_user_dir' for this package.

Usage

hydrogeofetch_data_dir(dir = NULL)

Arguments

Value

character path of data directory (silent when setting)

Examples

hydrogeofetch_data_dir()

# set it somewhere else, then put it back
old_dir <- hydrogeofetch_data_dir()

hydrogeofetch_data_dir(file.path(tempdir(check = TRUE), "demo"))

hydrogeofetch_data_dir()

hydrogeofetch_data_dir(old_dir)

Make isolated NHDPlusHR region a standalone dataset

Description

Cleans up and prepares NHDPlusHR regional data for use as complete NHDPlus data. The primary modification applied is to ensure that any flowpath that exits the domain is labeled as a terminal path and attributes are propagated upstream such that the domain is independently complete.

Usage

make_standalone(flowlines)

Arguments

Value

sf data.frame containing standalone network

Examples

library(dplyr)
library(sf)
source(system.file("extdata/nhdplushr_data.R", package = "hydrogeofetch"))

(outlet <- filter(hr_data$NHDFlowline, Hydroseq == min(Hydroseq)))
nrow(filter(hr_data$NHDFlowline, TerminalPa == outlet$Hydroseq))

hr_data$NHDFlowline <- make_standalone(hr_data$NHDFlowline)

(outlet <- filter(hr_data$NHDFlowline, Hydroseq == min(Hydroseq)))
nrow(filter(hr_data$NHDFlowline, TerminalPa == outlet$Hydroseq))

source(system.file("extdata/nhdplushr_data.R", package = "hydrogeofetch"))

# Remove mainstem and non-dendritic stuff.
subset <- filter(hr_data$NHDFlowline,
                        StreamLeve > min(hr_data$NHDFlowline$StreamLeve) &
                          StreamOrde == StreamCalc)

subset <- subset_nhdplus(subset$COMID, nhdplus_data = hr_gpkg)$NHDFlowline

plot(sf::st_geometry(hr_data$NHDFlowline))

flowline_mod <- make_standalone(subset)

terminals <- unique(flowline_mod$TerminalPa)

colors <- sample(hcl.colors(length(terminals), palette = "Zissou 1"))

for(i in 1:length(terminals)) {
  fl <- flowline_mod[flowline_mod$TerminalPa == terminals[i], ]
  plot(st_geometry(fl), col = colors[i], lwd = 2, add = TRUE)
}

ol <- filter(flowline_mod, TerminalFl == 1 & TerminalPa %in% terminals)

plot(st_geometry(ol), lwd = 2, add = TRUE)

Make Interactive Map of NHDPlus

Description

Given a list of outlets, get their basin boundaries and network and return a leaflet map in EPSG:4326.

Usage

map_nhdplus(
  outlets = NULL,
  bbox = NULL,
  streamorder = NULL,
  nhdplus_data = NULL,
  gpkg = NULL,
  flowline_only = NULL,
  plot_config = NULL,
  overwrite = TRUE,
  cache_data = NULL,
  return_map = FALSE
)

Arguments

Details

map_nhdplus supports several input specifications. An unexported function "as_outlet" is used to convert the outlet formats as described below.

1. if outlets is omitted, the bbox input is required and all nhdplus data in the bounding box is plotted.

2. If outlets is a list of integers, it is assumed to be NHDPlus IDs (comids) and all upstream tributaries are plotted.

3. if outlets is an integer vector, it is assumed to be all NHDPlus IDs (comids) that should be plotted. Allows custom filtering.

4. If outlets is a character vector, it is assumed to be NWIS site ids.

5. if outlets is a list containing only characters, it is assumed to be a list of nldi features and all upstream tributaries are plotted.

6. if outlets is a data.frame with point geometry, a point in polygon match is performed and upstream with tributaries from the identified catchments is plotted.

See plot_nhdplus for details on plot configuration.

Value

data.frame or leaflet map (see return_map)

Examples

map_nhdplus("05428500")

map_nhdplus("05428500", streamorder = 2)

map_nhdplus(list(13293970, 13293750))

source(system.file("extdata/sample_data.R", package = "hydrogeofetch"))

map_nhdplus(list(13293970, 13293750), streamorder = 3, nhdplus_data = sample_data)

#return leaflet object
map_nhdplus("05428500", return_map = TRUE)

Navigate NLDI

Description

Navigate the Network Linked Data Index network.

Usage

navigate_nldi(
  nldi_feature,
  mode = "upstreamMain",
  data_source = "flowlines",
  distance_km = 10
)

Arguments

Value

sf data.frame with result

Examples

library(sf)
library(dplyr)

nldi_nwis <- list(featureSource = "nwissite", featureID = "USGS-05428500")

navigate_nldi(nldi_feature = nldi_nwis,
              mode = "upstreamTributaries")$UT |>
  st_geometry() |>
  plot()

navigate_nldi(nldi_feature = nldi_nwis,
              mode = "UM")$UM |>
  st_geometry() |>
  plot(col = "blue", add = TRUE)



nwissite <- navigate_nldi(nldi_feature = nldi_nwis,
                          mode = "UT",
                          data_source = "nwissite")$UT_nwissite

st_geometry(nwissite) |>
  plot(col = "green", add = TRUE)

nwissite

NHDPlus Data Path

Description

Allows specification of a custom path to a source dataset. Typically this will be the national seamless dataset in geodatabase or geopackage format.

Usage

nhdplus_path(path = NULL, warn = FALSE)

Arguments

Value

0 (invisibly) if set successfully, character path if no input.

Examples

nhdplus_path("/data/NHDPlusV21_National_Seamless.gdb")

nhdplus_path("/data/NHDPlusV21_National_Seamless.gdb", warn=FALSE)

nhdplus_path()

Plot NHDPlus

Description

Given a list of outlets, get their basin boundaries and network and return a plot in EPSG:3857 Web Mercator Projection.

Usage

plot_nhdplus(
  outlets = NULL,
  bbox = NULL,
  streamorder = NULL,
  nhdplus_data = NULL,
  gpkg = NULL,
  plot_config = NULL,
  basemap = "Esri.NatGeoWorldMap",
  zoom = NULL,
  add = FALSE,
  actually_plot = TRUE,
  overwrite = TRUE,
  flowline_only = NULL,
  cache_data = NULL
)

Arguments

Details

plot_nhdplus supports several input specifications. An unexported function "as_outlet" is used to convert the outlet formats as described below.

1. if outlets is omitted, the bbox input is required and all nhdplus data in the bounding box is plotted.

2. If outlets is a list of integers, it is assumed to be NHDPlus IDs (comids) and all upstream tributaries are plotted.

3. if outlets is an integer vector, it is assumed to be all NHDPlus IDs (comids) that should be plotted. Allows custom filtering.

4. If outlets is a character vector, it is assumed to be NWIS site ids.

5. if outlets is a list containing only characters, it is assumed to be a list of nldi features and all upstream tributaries are plotted.

6. if outlets is a data.frame with point geometry, a point in polygon match is performed and upstream with tributaries from the identified catchments is plotted.

The plot_config parameter is a list with names "basin", "flowline", "outlets", "network_wtbd", and "off_network_wtbd". The following shows the defaults that can be altered.

1. basin

       list(lwd = 1, col = NA, border = "black")

2. flowline

       list(lwd = 1, col = "blue")

3. outlets

      list(default = list(col = "black", border = NA, pch = 19, cex = 1),
           nwissite = list(col = "grey40", border = NA, pch = 17, cex = 1),
           huc12pp = list(col = "white", border = "black", pch = 22, cex = 1),
           wqp = list(col = "red", border = NA, pch = 20, cex = 1))

4. network_wtbd

       list(lwd = 1, col = "lightblue", border = "black")

5. off_network_wtbd

       list(lwd = 1, col = "darkblue", border = "black")

If adding additional layers to the plot, data must be projected to EPSG:3857 with 'sf::st_transform(x, 3857)' prior to adding to the plot.

Value

data.frame plot data is returned invisibly in NAD83 Lat/Lon.

Examples

# Beware plot_nhdplus caches data to the default location.
# If you do not want data in "user space" change the default.
old_dir <- hydrogeofetch::hydrogeofetch_data_dir()
hydrogeofetch_data_dir(tempdir())

plot_nhdplus("05428500")

plot_nhdplus("05428500", streamorder = 2)

plot_nhdplus(list(13293970, 13293750))

source(system.file("extdata/sample_data.R", package = "hydrogeofetch"))

plot_nhdplus(list(13293970, 13293750), streamorder = 3, nhdplus_data = sample_data)

plot_nhdplus(list(list("comid", "13293970"),
                  list("nwissite", "USGS-05428500"),
                  list("huc12pp", "070900020603"),
                  list("huc12pp", "070900020602")),
             streamorder = 2,
             nhdplus_data = sample_data)


plot_nhdplus(sf::st_as_sf(data.frame(x = -89.36083,
                                     y = 43.08944),
                          coords = c("x", "y"), crs = 4326),
             streamorder = 2,
             nhdplus_data = sample_data)

plot_nhdplus(list(list("comid", "13293970"),
                  list("nwissite", "USGS-05428500"),
                  list("huc12pp", "070900020603"),
                  list("huc12pp", "070900020602")),
             streamorder = 2,
             nhdplus_data = sample_data,
             plot_config = list(basin = list(lwd = 2),
                                outlets = list(huc12pp = list(cex = 1.5),
                                               comid = list(col = "green"))))

bbox <- sf::st_bbox(c(xmin = -89.43, ymin = 43, xmax = -89.28, ymax = 43.1),
                    crs = "+proj=longlat +datum=WGS84 +no_defs")

fline <- sf::read_sf(sample_data, "NHDFlowline_Network")
comids <- hydroloom::navigate_hydro_network(fline, 13293970, "UT")

plot_nhdplus(comids)

#' # With Local Data
plot_nhdplus(bbox = bbox, nhdplus_data = sample_data)

# With downloaded data
plot_nhdplus(bbox = bbox, streamorder = 3)

# Can also plot on top of the previous!
plot_nhdplus(bbox = bbox, nhdplus_data = sample_data,
             plot_config = list(flowline = list(lwd = 0.5)))
plot_nhdplus(comids, nhdplus_data = sample_data, streamorder = 3, add = TRUE,
             plot_config = list(flowline = list(col = "darkblue")))

hydrogeofetch::hydrogeofetch_data_dir(old_dir)

Prep NHDPlus Data

Description

Function to prep NHDPlus data for use by hydrogeofetch functions

Usage

prepare_nhdplus(
  flines,
  min_network_size = 0,
  min_path_length = 0,
  min_path_size = 0,
  purge_non_dendritic = TRUE,
  warn = TRUE,
  error = TRUE,
  skip_toCOMID = FALSE,
  align_names = TRUE
)

Arguments

Value

data.frame with prepared NHDPlus network attributes.

Query USGS Hydro ESRI Rest Server

Description

Query the USGS Hydro ESRI Rest Server for spatial data by location, area, or ID.

Usage

query_usgs_arcrest(
  AOI = NULL,
  ids = NULL,
  type = NULL,
  service = NULL,
  where = NULL,
  t_srs = NULL,
  buffer = 0.5,
  page_size = 2000
)

Arguments

Details

The returned object(s) will have the same Spatial Reference System (SRS) as the input AOI. If a individual or set of IDs are used to query, then the default CRS of EPSG:4269 is preserved. In all cases, a user-defined SRS can be passed to t_srs which will override all previous SRS (either input or default). All buffer and distance operations are handled internally using in EPSG:5070 Albers Equal Area projection

Value

a simple features (sf) object or valid types if no type supplied

Query USGS Water OGC API Features

Description

Query the USGS Water OGC API for spatial data by location, area, or ID.

Usage

query_usgs_oafeat(
  AOI = NULL,
  ids = NULL,
  type = NULL,
  filter = NULL,
  t_srs = NULL,
  buffer = 0.5,
  properties = NULL,
  skip_geometry = FALSE
)

Arguments

Details

The returned object(s) will have the same Spatial Reference System (SRS) as the input AOI. If a individual or set of IDs are used to query, then the default server CRS of EPSG:4326 is preserved. In all cases, a user-defined SRS can be passed to t_srs which will override all previous SRS (either input or default). All buffer and distance operations are handled internally using an EPSG:5070 Albers Equal Area projection

Value

a simple features (sf) object, or a data.frame when skip_geometry = TRUE

Rescale Catchment Characteristics

Description

Given catchment characteristics to retrieve or process will aggregate and / or split the characteristics according to a lookup table.

Usage

rescale_catchment_characteristics(
  vars,
  lookup_table,
  refactored_areas = NULL,
  catchment_characteristics = NULL,
  catchment_areas = NULL,
  source = "usgs",
  aoi = "cat"
)

Arguments

Details

NOTE: Since this algorithm works on catchment characteristics that are spatial averages, when splitting, the average condition is apportioned evenly to each split. In some cases, such as with land cover or elevation, this may not be appropriate and source data should be used to derive new characteristics. In addition, this function handles catchment areas for split catchments but makes no adjustments for the length of flowlines in those catchments. Therefore, requests for length-weighted mean values may not be appropriate when working with split catchments.

Value

data.frame with one row per unique 'id' in 'lookup_table'. Always contains 'id', 'areasqkm_sum' (total catchment area contributing to that 'id'), and 'lengthkm_sum' (total flowline length). Each characteristic named in 'vars' adds a column called '<characteristic_id>_<statistic>', where the statistic suffix is one of "area_wtd", "length_wtd", "sum", "min", or "max" following the 'summary_statistic' requested for it – for example 'CAT_IMPV11_area_wtd' or 'CAT_BASIN_AREA_sum'. Where the source reports a no-data fraction, a matching 'percent_nodata_<characteristic_id>_area_wtd' column is included. Returns 'NULL' if 'catchment_characteristics' was not supplied and could not be retrieved from the web service.

Examples

vars <- data.frame(characteristic_id = c("CAT_IMPV11","CAT_BASIN_AREA"),
                   summary_statistic = c("area_weighted_mean","sum"))
lookup_table <- data.frame(id = rep(10012268, 2),
                           comid = c(4146596, 4147382),
                           member_comid = c(4146596, 4147382))
rescale_catchment_characteristics(vars, lookup_table)

vars <- data.frame(characteristic_id = c("CAT_ELEV_MIN","CAT_ELEV_MAX"),
                   summary_statistic = c("min","max"))
lookup_table <- data.frame(id = rep(10012268, 2),
                           comid = c(4146596, 4147382),
                           member_comid = c(4146596, 4147382))
rescale_catchment_characteristics(vars, lookup_table)

vars <- data.frame(characteristic_id = c("CAT_EWT","CAT_TWI", "CAT_BASIN_AREA"),
                   summary_statistic = c("area_weighted_mean", "area_weighted_mean","sum"))
lookup_table <- data.frame(id = c(10012268, 10012268, 10024047, 10024048),
                           comid = c(4146596, 4147382, 4147396, 4147396),
                           member_comid = c("4146596", "4147382", "4147396.1", "4147396.2"))
comid_areas <- data.frame(featureid = c("4146596", "4147382", "4147396.1", "4147396.2"),
                               areasqkm = c(0.9558, 11.9790, 6.513294, 1.439999))
rescale_catchment_characteristics(vars, lookup_table, refactored_areas = comid_areas)

 

RPU Boundaries Raster Processing Unit boundaries

Description

RPU Boundaries Raster Processing Unit boundaries

Usage

rpu_boundaries

Format

An object of class "sf"

Subset NHDPlus

Description

Saves a subset of the National Seamless database or other hydrogeofetch compatible data based on a specified collection of COMIDs. This function uses get_nhdplus for the "download" data source but returns data consistent with local data subsets in a subset file.

Usage

subset_nhdplus(
  comids = NULL,
  output_file = NULL,
  nhdplus_data = NULL,
  bbox = NULL,
  simplified = TRUE,
  overwrite = FALSE,
  return_data = TRUE,
  status = TRUE,
  flowline_only = NULL,
  streamorder = NULL,
  out_prj = 4269
)

Arguments

Details

This function relies on the National Seamless Geodatabase or Geopackage. It can be downloaded here.

The "download" option of this function should be considered preliminary and subject to revision. It does not include as many layers and may not be available permanently.

Value

character path to the saved subset geopackage

Examples

source(system.file("extdata/sample_data.R", package = "hydrogeofetch"))

nhdplus_path(sample_data)

sample_flines <- sf::st_zm(sf::read_sf(nhdplus_path(), "NHDFlowline_Network"))

plot(sf::st_geometry(sample_flines),
     lwd = 3)

start_point <- sf::st_sfc(sf::st_point(c(-89.362239, 43.090266)),
                          crs = 4326)

plot(start_point, cex = 1.5, lwd = 2, col = "red", add = TRUE)

start_comid <- discover_nhdplus_id(start_point)

comids <- hydroloom::navigate_hydro_network(sample_flines, start_comid, "UT")

plot(sf::st_geometry(dplyr::filter(sample_flines, COMID %in% comids)),
     add=TRUE, col = "red", lwd = 2)

output_file <- tempfile(fileext = ".gpkg")

subset_nhdplus(comids = comids,
               output_file = output_file,
               nhdplus_data = sample_data,
               overwrite = TRUE,
               status = TRUE)

sf::st_layers(output_file)

catchment <- sf::read_sf(output_file, "CatchmentSP")

plot(sf::st_geometry(catchment), add = TRUE)

waterbody <- sf::read_sf(output_file, "NHDWaterbody")

plot(sf::st_geometry(waterbody),
     col = rgb(0, 0, 1, alpha = 0.5), add = TRUE)

# Cleanup temp
unlink(output_file)

# Download Option:
subset_nhdplus(comids = comids,
               output_file = output_file,
               nhdplus_data = "download",
               overwrite = TRUE,
               status = TRUE, flowline_only = FALSE)

sf::st_layers(output_file)

# NHDPlusHR
source(system.file("extdata/nhdplushr_data.R", package = "hydrogeofetch"))

up_ids <- hydroloom::navigate_hydro_network(hr_data$NHDFlowline, 15000500028335, "UT")

sub_gpkg <- file.path(work_dir, "sub.gpkg")
sub_nhdhr <- subset_nhdplus(up_ids, output_file = sub_gpkg,
                            nhdplus_data = hr_gpkg, overwrite = TRUE)

sf::st_layers(sub_gpkg)
names(sub_nhdhr)

plot(sf::st_geometry(hr_data$NHDFlowline), lwd = 0.5)
plot(sf::st_geometry(sub_nhdhr$NHDFlowline), lwd = 0.6, col = "red", add = TRUE)

unlink(output_file)
unlink(sub_gpkg)

Subset by Raster Processing Unit

Description

Given flowlines and an rpu_code, performs a network-safe subset such that the result can be used in downstream processing. Has been tested to work against the entire NHDPlusV2 domain and satisfies a number of edge cases.

Usage

subset_rpu(fline, rpu, run_make_standalone = TRUE, strict = FALSE)

Arguments

Value

data.frame containing subset network

Examples

source(system.file("extdata/sample_data.R", package = "hydrogeofetch"))

sample_flines <- sf::read_sf(sample_data, "NHDFlowline_Network")

subset_rpu(sample_flines, rpu = "07b")

Subset by Vector Processing Unit

Description

Calls subset_rpu for all raster processing units for the requested vector processing unit.

Usage

subset_vpu(fline, vpu, include_null_rpuid = TRUE, run_make_standalone = TRUE)

Arguments

Value

data.frame containing subset network

Examples

source(system.file("extdata/sample_data.R", package = "hydrogeofetch"))

sample_flines <- sf::read_sf(sample_data, "NHDFlowline_Network")

subset_vpu(sample_flines, "07")

Update superseded mainstem identifiers

Description

Given point features with a known geoconnex mainstem identifier, checks whether the mainstem has been superseded and, if so, re-indexes the point to the replacement mainstem using index_points_to_lines. Points whose mainstem has no replacement are left unchanged and flagged "unresolved". See check_mainstems for a lightweight superseded/current check.

Usage

update_mainstems(x, mainstem_col = NULL, search_radius = NULL)

Arguments

Value

x with mainstem_col updated in place and a mainstem_update_status column added with values "unchanged", "updated", or "unresolved".

Examples

pt <- sf::st_sf(mainstemid = 2086165,
                geometry = sf::st_sfc(sf::st_point(c(-75.567, 43.176)),
                                      crs = 4326))
update_mainstems(pt)

VPU Boundaries Vector Processing Unit boundaries

Description

VPU Boundaries Vector Processing Unit boundaries

Usage

vpu_boundaries

Format

An object of class "sf"