- add_areasqkm
- add_lengthmap
- aggregate_catchments
- aggregate_network
- clean_geometry
- collapse_flowlines
- download_elev
- download_fdr_fac
- flowpaths_to_linestrings
- get_minimal_network
- get_row_col
- map_outlet_ids
- prep_cat_fdr_fac
- reconcile_catchment_divides
- reconcile_collapsed_flowlines
- refactor_nhdplus
- split_catchment_divide
- split_flowlines
- trace_upstream
- union_linestrings
- union_linestrings_geos
- union_polygons
September 22, 2022
Package: hyRefactor
Type: Package
Title: Hydrologic Network Refactoring Tools Based on HY_Features
Version: 0.4.8
Authors@R: c(person("David", "Blodgett", role = c("aut", "cre"), email = "dblodgett@usgs.gov"),
person(given = "Mike", family = "Johnson", role = "aut", comment = c(ORCID = "0000-0002-5288-8350")))
Description: A collection of tools for normalizing and changing the size of catchments.
URL: https://github.com/dblodgett-usgs/hyrefactor
BugReports: https://github.com/dblodgett-usgs/hyrefactor/issues
Depends:
R (>= 3.5.0)
Imports:
dplyr,
sf,
lwgeom,
units,
data.table,
terra,
nhdplusTools,
tidyr,
pbapply,
rvest,
httr,
methods,
rlang,
rmapshaper,
utils
Suggests: testthat, knitr, rmarkdown
Remotes: dblodgett-usgs/nhdplusTools
License: CC0
Encoding: UTF-8
LazyData: true
RoxygenNote: 7.2.1
VignetteBuilder: knitr
Compute km2 area Short hand for safely computing area in sqkm and returning as numeric vector.
Compute km2 area Short hand for safely computing area in sqkm and returning as numeric vector.
add_areasqkm(x)| Argument | Description |
|---|---|
x |
sf object |
numeric vector
library(sf)
nc = st_read(system.file("shape/nc.shp", package="sf"))
add_areasqkm(nc[1,])Add Length Map to Refactored Network
This function replicates the member_COMID column of a refactored network but adds a new notation Following each COMID is '.' which is proceeded by the fraction of that COMID making up the new flowpath. For example 101.1 would indicate 100 % of COMID 101 is in the new ID. Equally 101.05 would indicate 50 % of COMID 101 is present in the new ID'd flowpath
add_lengthmap(flowpaths, length_table)| Argument | Description |
|---|---|
flowpaths |
a refactored flowpath network containing an member_COMID column |
length_table |
a table of NHDPlus COMIDs and LENGTH to use as weights. Can be found with nhdplusTools::get_vaa("lengthkm") |
sf object
path <- system.file("extdata/walker_reconcile.gpkg", package = "hyRefactor")
fps <- add_lengthmap(flowpaths = sf::read_sf(path),
length_table = nhdplusTools::get_vaa("lengthkm"))Aggregate Catchments
Aggregates catchments according to a set of outlet catchments.
Network aggregation is completed using: See aggregate_network .
aggregate_catchments(
flowpath,
divide,
outlets,
zero_order = NULL,
coastal_cats = NULL,
da_thresh = NA,
only_larger = FALSE,
post_mortem_file = NA,
keep = NULL
)| Argument | Description |
|---|---|
flowpath |
sf data.frame Flowpaths as generated by refactor_nhdplus |
divide |
sf data.frame Reconciled catchment divides as generated by reconcile_catchment_divides |
outlets |
data.frame with "ID" and "type" columns. "ID" must be identifiers from flowpath and divide data.frames. "type" should be "outlet", or "terminal". "outlet" will include the specified ID. "terminal" will be treated as a terminal node with nothing downstream. |
zero_order |
list of vectors containing IDs to be aggregated into 0-order catchments. |
coastal_cats |
sf data.frame with coastal catchments to be used with zero order. |
da_thresh |
numeric See aggregate_network |
only_larger |
boolean See aggregate_network |
post_mortem_file |
rda file to dump environment to in case of error |
keep |
logical passed along to clean_geometry |
source(system.file("extdata", "walker_data.R", package = "hyRefactor"))
outlets <- data.frame(ID = c(31, 3, 5, 1, 45, 92),
type = c("outlet", "outlet", "outlet", "terminal", "outlet", "outlet"),
stringsAsFactors = FALSE)
aggregated <- aggregate_catchments(flowpath = walker_fline_rec,
divide = walker_catchment_rec,
outlets = outlets)
plot(aggregated$cat_sets$geom, lwd = 3, border = "red")
plot(walker_catchment_rec$geom, lwd = 1.5, border = "green", col = NA, add = TRUE)
plot(walker_catchment$geom, lwd = 1, add = TRUE)
plot(walker_flowline$geom, lwd = .7, col = "blue", add = TRUE)
plot(aggregated$cat_sets$geom, lwd = 3, border = "black")
plot(aggregated$fline_sets$geom, lwd = 3, col = "red", add = TRUE)
plot(walker_flowline$geom, lwd = .7, col = "blue", add = TRUE)Aggregate Network
Aggregates a catchment network according to a set of outlet.
aggregate_network(
flowpath,
outlets,
da_thresh = NA,
only_larger = FALSE,
post_mortem_file = NA
)| Argument | Description |
|---|---|
flowpath |
sf data.frame Flowpaths with ID, toID, LevelPathID, and Hydroseq attributes. |
outlets |
data.frame with "ID" and "type" columns. "ID" must be identifiers from fowpath and divide data.frames. "type" should be "outlet", or "terminal". "outlet" will include the specified ID. "terminal" will be treated as a terminal node with nothing downstream. |
da_thresh |
numeric Defaults to NA. A threshold total drainage area in the units of the TotDASqKM field of the flowpath data.frame. When automatically adding confluences to make the network valid, tributary catchments under this threshold will be lumped with the larger tributaries rather than being added to the set of output catchments. |
only_larger |
boolean Defaults to TRUE. If TRUE when adding confluences to make the network valid, only tributaries larger than the one with an upstream outlet will be added. e.g. if a tributary is required in the model this will add main stems that the tributary contributes to. Note that the NHDPlus treats divergences as part of the main stem, so the da_thresh may still be needed to eliminate small tributary catchments introduced by divergences near confluences. |
post_mortem_file |
rda file to dump environment to in case of error |
This function operates on the catchment network as a node-edge graph. The outlet types are required to ensure that graph searches start from the appropriate nodes and includes the appropriate catchments. Outlets such as gages should be treated as "outlet" outlets. While it may be possible for the algorithm to determine terminal outlets, at this time, it is required that they be specified explicitely as "terminal" outlet types.
The function checks supplied outlets to make sure they connect downstream. Checks verify that the outlet of the levelpath (main stem of a total catchment) of each supplied outlet is in the supplied outlet set. If the outlet of a levelpath is not in the supplied set, it is added along with other catchments that contribute to the same receiving catchment. These checks ensure that all output catchments have one and only one input and output nexus and that all catchments are well-connected.
source(system.file("extdata", "walker_data.R", package = "nhdplusTools"))
fline <- dplyr::right_join(dplyr::select(walker_flowline, COMID),
nhdplusTools::prepare_nhdplus(walker_flowline, 0, 0, 0, FALSE))
fline <- dplyr::select(fline, ID = COMID, toID = toCOMID,
LevelPathID = LevelPathI, Hydroseq)
outlets <- data.frame(ID = c(5329357, 5329317, 5329365, 5329303, 5329435, 5329817),
type = c("outlet", "outlet", "outlet", "terminal", "outlet", "outlet"),
stringsAsFactors = FALSE)
aggregated <- aggregate_network(fline, outlets)
aggregated <- aggregate_network(fline, outlets)
outlets <- dplyr::filter(fline, ID %in% outlets$ID)
outlets <- nhdplusTools::get_node(outlets)
plot(aggregated$fline_sets$geom, lwd = 3, col = "red")
plot(walker_flowline$geom, lwd = .7, col = "blue", add = TRUE)
plot(outlets$geometry, add = TRUE)Clean Catchment Geometry
Fixes geometry issues present in catchments that originate in the
CatchmentSP layers, or from the reconcile_catchments hyRefactor preocess.
These include, but are not limited to disjoint polygon fragments, artifacts
from the 30m DEM used to generate the catchments, and non-valid geometry topolgies.
A goal of this functions is also to provide means to reduce the data column
of the catchments by offering a topology preserving simplification
through ms_simplify .
Generally a "keep" parameter of .9 seems appropriate for the resolution of
the data but can be modified in function
clean_geometry(catchments, ID = "ID", keep = 0.9, crs = 5070, sys = NULL)| Argument | Description |
|---|---|
catchments |
catchments geometries to fix |
ID |
name of uniquely identifying column |
keep |
proportion of points to retain in geometry simplification (0-1; default 0.05). See ms_simplify . If NULL, then no simplification will be executed. |
crs |
integer or object compatible with sf::st_crs coordinate reference. Should be a projection that supports area-calculations. |
sys |
logical should the mapshaper system library be used. If NULL the system library will be used if available. |
sf object
Collapse NHDPlus Network
Refactors the NHDPlus flowline network, eliminating short and non-confluence flowlines. The aim of this function is to create flowpaths that describe a network of catchments that combines complex hydrology near confluences into upstream catchments and removes very short flowlines along mainstem flow-paths.
collapse_flowlines(
flines,
thresh,
add_category = FALSE,
mainstem_thresh = NULL,
exclude_cats = NULL,
warn = TRUE
)| Argument | Description |
|---|---|
flines |
data.frame with COMID, toCOMID, LENGTHKM, Hydroseq, and LevelPathI columns |
thresh |
numeric a length threshold (km). Flowlines shorter than this will be collapsed with up or downstream flowlines. |
add_category |
boolean if combination category is desired in output, set to TRUE |
mainstem_thresh |
numeric threshold for combining inter-confluence mainstems |
exclude_cats |
integer vector of COMIDs to be excluded from collapse modifications. |
warn |
boolean controls whether warning an status messages are printed |
A refactored network with merged up and down flowlines.
The refactor_nhdplus function implements a complete
workflow using collapse_flowlines().
Download Elevation and Derivatives
Download Elevation and Derivatives
download_elev(product, out_dir, regions = NULL)| Argument | Description |
|---|---|
product |
character DEM, hydroDEM, or FDRFAC. |
out_dir |
path to directory to store output. |
regions |
character vector of two digit hydrologic |
Download FDR FAC
Download FDR FAC
download_fdr_fac(out_dir, regions = NULL)| Argument | Description |
|---|---|
out_dir |
path to directory to store output. |
regions |
character vector of two digit hydrologic |
Convert MULITLINESTINGS to LINESTRINGS
Convert MULITLINESTINGS to LINESTRINGS
flowpaths_to_linestrings(flowpaths)| Argument | Description |
|---|---|
flowpaths |
a flowpath sf object |
a sf object
Get Minimal Network
Given a set of outlets, will generate a minimal network by
calling aggregate_network and adding nhdplus attributes to the result.
If geometry is included with the network, it will be merged and returned.
get_minimal_network(flowpath, outlets)| Argument | Description |
|---|---|
flowpath |
sf data.frame Flowpaths with ID, toID, LevelPathID, Hydroseq and LENGTHKM and AreaSqKM attributes. |
outlets |
data.frame with "ID" and "type" columns. "ID" must be identifiers from fowpath and divide data.frames. "type" should be "outlet", or "terminal". "outlet" will include the specified ID. "terminal" will be treated as a terminal node with nothing downstream. |
a data.frame (potentially including an sfc list column) with
attributes generated by add_plus_network_attributes
and a list column "set" containing members of each output flowpath.
source(system.file("extdata", "walker_data.R", package = "nhdplusTools"))
fline <- walker_flowline
outlets <- data.frame(ID = c(5329357, 5329317, 5329365, 5329435, 5329817),
type = c("outlet", "outlet", "outlet", "outlet", "outlet"))
#' Add toCOMID
fline <- nhdplusTools::get_tocomid(fline, add = TRUE)
# get attributes set
fline <- dplyr::select(fline, ID = comid, toID = tocomid,
LevelPathID = levelpathi, hydroseq = hydroseq,
AreaSqKM = areasqkm, LENGTHKM = lengthkm)
min_net <- get_minimal_network(fline, outlets)
plot(sf::st_geometry(fline), col = "blue")
plot(sf::st_geometry(min_net), lwd = 2, add = TRUE)
plot(sf::st_geometry(nhdplusTools::get_node(min_net)), add = TRUE)Get Row and Column
Get Row and Column
get_row_col(fdr, start, fac_matrix)| Argument | Description |
|---|---|
fdr |
flow direction grid |
start |
matrix (row, col) |
fac_matrix |
flow accumulation matrix |
Map outlets from COMID to ID for aggregate catchments
given reconciled flowlines and a set of source outlets, returns a set of outlets with reconciled IDs suitable for use with aggregate_catchments.
map_outlet_ids(source_outlets, reconciled)| Argument | Description |
|---|---|
source_outlets |
data.frame with COMID and type columns |
reconciled |
data.frame as returned by refactor workflow |
Prep catchment with FDR/FAC
Prep catchment with FDR/FAC
prep_cat_fdr_fac(cat, fdr, fac)| Argument | Description |
|---|---|
cat |
catchment (sf object) |
fdr |
flow direction grid |
fac |
flow accumulation grid |
Reconcile Catchment Divides
Reconciles catchment divides according to the output of
reconcile_collapsed_flowlines and refactor_nhdplus
reconcile_catchment_divides(
catchment,
fline_ref,
fline_rec,
fdr = NULL,
fac = NULL,
para = 2,
cache = NULL,
min_area_m = 800,
snap_distance_m = 100,
simplify_tolerance_m = 40,
vector_crs = 5070,
fix_catchments = TRUE,
keep = 0.9
)| Argument | Description |
|---|---|
catchment |
sf data.frame NHDPlus Catchment or CatchmentSP layers for included COMIDs |
fline_ref |
sf data.frame flowlines as returned by refactor_nhdplus and reconcile_collapsed_flowlines |
fline_rec |
sf data.frame flowpaths as returned by reconcile_collapsed_flowlines |
fdr |
character path to D8 flow direction |
fac |
character path to flow accumulation |
para |
integer numer of cores to use for parallel execution |
cache |
path .rda to cache incremental outputs |
min_area_m |
minimum area in m^2 to filter out slivers (caution, use with care!!) |
snap_distance_m |
distance in meters to snap SpatRaster generated geometry to polygon geometry |
simplify_tolerance_m |
dTolerance in meters for simplification of grid-cell based polygons |
vector_crs |
integer or object compatible with sf::st_crs coordinate reference. Should be a projection that supports area-calculations. |
fix_catchments |
logical. should catchment geometries be rectified? |
keep |
Only applicable if fix_catchments = TRUE. Defines the proportion of points to retain in geometry simplification (0-1; default 0.05). See ms_simplify . Set to NULL to skip simplification. |
Note that all inputs must be passed in the same projection.
Catchment divides that have been split and collapsed according to input flowpaths
The refactor_nhdplus function implements a complete
workflow using reconcile_collapsed_flowlines() and can be used in prep
for this function.
Reconcile Collapsed Flowlines
Reconciles output of collapse_flowlines giving a unique ID to each new flowpath and providing a mapping to NHDPlus COMIDs.
reconcile_collapsed_flowlines(flines, geom = NULL, id = "COMID")| Argument | Description |
|---|---|
flines |
data.frame with COMID, toCOMID, LENGTHKM, LevelPathI, Hydroseq, and TotDASqKM columns |
geom |
sf data.frame for flines |
id |
character id collumn name. |
reconciled flowpaths with new ID, toID, LevelPathID, and Hydroseq identifiers. Note that all the identifiers are new integer IDs. LevelPathID and Hydroseq are consistent with the LevelPathID and Hydroseq from the input NHDPlus flowlines.
The refactor_nhdplus function implements a complete
workflow using reconcile_collapsed_flowlines().
Refactor NHDPlus
A complete network refactor workflow has been packaged into this function. Builds a set of normalized catchment-flowpaths from input flowline features. See details and vignettes for more information.
refactor_nhdplus(
nhdplus_flines,
split_flines_meters,
split_flines_cores,
collapse_flines_meters,
collapse_flines_main_meters,
out_refactored,
out_reconciled,
three_pass = FALSE,
purge_non_dendritic = TRUE,
exclude_cats = NULL,
events = NULL,
warn = TRUE
)| Argument | Description |
|---|---|
nhdplus_flines |
data.frame raw nhdplus flowline features as derived from the national seamless geodatabase. |
split_flines_meters |
numeric the maximum length flowpath desired in the output. |
split_flines_cores |
numeric the number of processing cores to use while splitting flowlines. |
collapse_flines_meters |
numeric the minimum length of inter-confluence flowpath desired in the output. |
collapse_flines_main_meters |
numeric the minimum length of between-confluence flowpaths. |
out_refactored |
character where to write a geopackage containing the split and collapsed flowlines. |
out_reconciled |
character where to write a geopackage containing the reconciled flowpaths. |
three_pass |
boolean whether to perform a three pass collapse or single pass. |
purge_non_dendritic |
boolean passed on to prepare_nhdplus |
exclude_cats |
integer vector of COMIDs to be excluded from collapse modifications. |
events |
data.frame containing events as generated by nhdplusTools::get_flowline_index() |
warn |
boolean controls whether warning an status messages are printed |
This is a convenient wrapper function that implements three phases
of the network refactor workflow: split, collapse, reconcile. See the
NHDPlus Refactor vignette for details of these three steps by running:
vignette("refactor_nhdplus", package = "hyRefactor")
In addition to prepare_nhdplus from the nhdplusTools package,
The following three functions are used in the refactor_nhdplus workflow.
source(system.file("extdata",
"sample_flines.R",
package = "nhdplusTools"))
nhdplus_flowlines <- sf::st_zm(sample_flines)
refactor_nhdplus(nhdplus_flines = nhdplus_flowlines,
split_flines_meters = 2000,
split_flines_cores = 2,
collapse_flines_meters = 500,
collapse_flines_main_meters = 500,
out_refactored = "temp.gpkg",
out_reconciled = "temp_rec.gpkg",
three_pass = TRUE,
purge_non_dendritic = FALSE,
warn = FALSE)
unlink("temp.gpkg")
unlink("temp_rec.gpkg")Split Catchment Divides
A catchment-divide splitting algorithm that works with a D8 flow direction grid and the output of nhdplus_refactor. See Vignette for examples.
split_catchment_divide(
catchment,
fline,
fdr,
fac,
lr = FALSE,
min_area_m = 800,
snap_distance_m = 100,
simplify_tolerance_m = 40,
vector_crs = NULL
)| Argument | Description |
|---|---|
catchment |
sf data.frame with one catchment divide |
fline |
sf data.frame with one or more flowline segments in upstream downstream order. |
fdr |
character path to flow direction that fully covers the catchment |
fac |
character path to flow accumulation that fuller covers the catchment |
lr |
boolean should catchments be split along the left/right bank? |
min_area_m |
minimum area in m^2 to filter out slivers (caution, use with care!!) |
snap_distance_m |
distance in meters to snap SpatRaster generated geometry to polygon geometry |
simplify_tolerance_m |
dTolerance in meters for simplification of grid-cell based polygons |
vector_crs |
any object compatible with sf::st_crs. Used for vector-based calculations in case that fdr projection is not suitable (e.g. lon/lat) -- must result in units of meters. |
Split catchment divides as an sfc geometry.
Split Flowlines
A wrapper for split_lines that works on nhdplus attributes
split_flowlines(flines, max_length = NULL, events = NULL, para = 0, avoid = NA)| Argument | Description |
|---|---|
flines |
data.frame with COMID, toCOMID, LENGTHKM and LINESTRING sf column in "meters" projection |
max_length |
maximum segment length to return |
events |
data.frame containing events as generated by nhdplusTools::get_flowline_index() if an identifier attribute is included, it will be passed through in the output table. |
para |
numeric how many threads to use in parallel computation |
avoid |
vector of ids to avoid |
All the flowlines with some split apart. COMIDs are returned as strings with a semantic part number appended. That is .1, .2, ... .10, .11, etc. are appended and must be treated as one would treat a semantic version. .1 is the most upstream and the sequence increases in the downstream direction.
The refactor_nhdplus function implements a complete
workflow using split_flowlines().
source(system.file("extdata", "new_hope_data.R", package = "hyRefactor"))
new_hope_flowline <-
dplyr::right_join(dplyr::select(new_hope_flowline, COMID, REACHCODE, FromMeas, ToMeas),
suppressWarnings(nhdplusTools::prepare_nhdplus(
new_hope_flowline, 0, 0, 0, FALSE, warn = FALSE)),
by = "COMID")
split <- split_flowlines(suppressWarnings(sf::st_cast(sf::st_transform(
new_hope_flowline, 5070), "LINESTRING")),
max_length = 2000, events = new_hope_events)Trace Upstream
Trace Upstream
trace_upstream(start_point, cat, fdr, fac_matrix, fdr_matrix)| Argument | Description |
|---|---|
start_point |
row col index |
cat |
catchment |
fdr |
flow direction grid |
fac_matrix |
flow accumulation matrix |
fdr_matrix |
flow direction matrix |
sfc
Fast LINESTRING union
Wayyyy faster then either data.table, or sf based line merging
union_linestrings(lines, ID)| Argument | Description |
|---|---|
lines |
lines to merge |
ID |
ID to merge over |
an sf object
DEPRECATED: Fast LINESTRING union
Wayyyy faster then either data.table, or sf based line merging
union_linestrings_geos(lines, ID)| Argument | Description |
|---|---|
lines |
lines to merge |
ID |
ID to merge over |
an sf object
Fast POLYGON Union
This is significantly faster then sf::st_union or summarize
union_polygons(poly, ID)| Argument | Description |
|---|---|
poly |
sf POLYGON object |
ID |
the column name over which to union geometries |
sf object