#' @title Prep NHDPlus Data
#' @description Function to prep NHDPlus data for use by nhdplusTools functions
#' @param flines data.frame NHDPlus flowlines including:
#' COMID, LENGTHKM, FTYPE (or FCODE), TerminalFl, FromNode, ToNode, TotDASqKM,
#' StartFlag, StreamOrde, StreamCalc, TerminalPa, Pathlength,
#' and Divergence variables.
#' @param min_network_size numeric Minimum size (sqkm) of drainage network
#' to include in output.
#' @param min_path_length numeric Minimum length (km) of terminal level
#' path of a network.
#' @param min_path_size numeric Minimum size (sqkm) of outlet level
#' path of a drainage basin. Drainage basins with an outlet drainage area
#' smaller than this will be removed.
#' @param purge_non_dendritic logical Should non dendritic paths be removed
#' or not.
#' @param warn logical controls whether warning an status messages are printed
#' @param error logical controls whether to return potentially invalid data with a warning rather than an error
#' @param skip_toCOMID logical if TRUE, toCOMID will not be added to output.
#' @param align_names logical
#' @return data.frame ready to be used with the refactor_flowlines function.
#' @importFrom dplyr select filter left_join group_split group_by bind_rows
#' @export
#' @examples
#'
#' source(system.file("extdata", "sample_flines.R", package = "nhdplusTools"))
#'
#' prepare_nhdplus(sample_flines,
#' min_network_size = 10,
#' min_path_length = 1,
#' warn = FALSE)
#'
prepare_nhdplus <- function(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) {
req_names <- c("COMID", "LENGTHKM", "TerminalFl",
"Hydroseq", "TotDASqKM", "StartFlag", "TerminalPa",
"Pathlength", "Divergence",
"LevelPathI")
dend <- c("StreamOrde", "StreamCalc")
sz <- c("TotDASqKM", "Hydroseq", "LevelPathI")
tocomid <- c("ToNode", "FromNode")
dont_return <- tolower(c("ToNode", "FromNode", "TerminalFl",
"StartFlag", "StreamOrde", "StreamCalc",
"TerminalPa", "Pathlength", "Divergence"))
if(purge_non_dendritic) {
req_names <- c(req_names, dend)
}
if(min_path_size > 0 | min_network_size > 0 | min_path_length > 0) {
req_names <- c(req_names, sz)
}
if(!skip_toCOMID) {
req_names <- c(req_names, tocomid)
}
assign("prepare_nhdplus_attributes",
unique(req_names),
envir = nhdplusTools_env)
flines <- check_names(flines, "prepare_nhdplus", tolower = TRUE)
if ("sf" %in% class(flines)) {
if (warn) warning("removing geometry")
flines <- st_drop_geometry(flines)
}
orig_rows <- nrow(flines)
flines <- filter_coastal(flines)
flines <- flines[ , names(flines) %in% tolower(req_names)]
if (!any(flines$terminalfl == 1)) {
warning("Got NHDPlus data without a Terminal catchment. Attempting to find it.")
if (all(flines$terminalpa == flines$terminalpa[1])) {
out_ind <- which(flines$hydroseq == min(flines$hydroseq))
flines$terminalfl[out_ind] <- 1
} else {
stop("Multiple networks without terminal flags found. Can't proceed.")
}
}
if (purge_non_dendritic) {
flines <- filter(flines, .data$streamorde == .data$streamcalc)
} else {
if(!skip_toCOMID) {
flines[["fromnode"]][which(flines$divergence == 2)] <- NA
}
}
if(min_path_size > 0) {
remove_paths <- group_by(flines, .data$levelpathi)
remove_paths <- filter(remove_paths, .data$hydroseq == min(.data$hydroseq))
remove_paths <- filter(remove_paths, .data$totdasqkm < min_path_size &
.data$totdasqkm >= 0)$levelpathi
flines <- filter(flines, !.data$levelpathi %in% remove_paths)
}
if(min_network_size > 0 | min_path_length > 0) {
terminal_filter <- flines$terminalfl == 1 &
flines$totdasqkm < min_network_size
start_filter <- flines$startflag == 1 &
flines$pathlength < min_path_length
if (any(terminal_filter, na.rm = TRUE) | any(start_filter, na.rm = TRUE)) {
tiny_networks <- rbind(filter(flines, terminal_filter),
filter(flines, start_filter))
flines <- filter(flines, !flines$terminalpa %in%
unique(tiny_networks$terminalpa))
}
}
if (warn) {
warning(paste("Removed", orig_rows - nrow(flines),
"flowlines that don't apply.\n",
"Includes: Coastlines, non-dendritic paths, \nand networks",
"with drainage area less than",
min_network_size, "sqkm, and drainage basins smaller than",
min_path_size))
}
if(!skip_toCOMID) {
if(nrow(flines) > 0) {
flines <- get_tocomid(flines, remove_coastal = FALSE)
}
if (!all(flines[["terminalfl"]][which(flines$tocomid == 0)] == 1)) {
warn <- paste("FromNode - ToNode imply terminal flowlines that are not\n",
"flagged terminal. Can't assume NA toCOMIDs go to the ocean.")
if(error) {
stop(warn)
} else {
warning(warn)
}
}
}
out <- flines[ , !names(flines) %in% dont_return]
if(align_names) {
out <- align_nhdplus_names(out)
}
out
}
filter_coastal <- function(flines) {
if("fcode" %in% names(flines)) {
flines <- filter(flines, (.data$fcode != 56600))
} else if("ftype" %in% names(flines)) {
flines <- filter(flines, (.data$ftype != "Coastline" | .data$ftype != 566))
} else {
stop("must provide fcode and/or ftype to filter coastline.")
}
}
#' 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 \link[hydroloom]{add_toids} function.
#' @param x data.frame with comid, tonode, fromnode, and (optionally)
#' divergence and terminalpa attributes.
#' @param return_dendritic logical if TRUE, a divergence attribute is required
#' (2 indicates diverted path, 1 is main) and diverted paths will be treated
#' as headwaters. If this is FALSE, the return value is a data.frame including
#' the comid and tocomid attributes.
#' @param missing integer value to use for terminal nodes.
#' @param remove_coastal logical remove coastal features prior to generating
#' tocomid values? ftype or fcode are required if `TRUE`. fcode == 56600 or
#' fcode == "Coastline" will be removed.
#' @param add logical if TRUE, a tocomid column will be added, otherwise
#' a data.frame with two columns will be returned.
#' @return data.frame containing comid and tocomid attributes or all
#' attributes provided with comid and tocomid in the first and second columns..
#' @importFrom hydroloom add_toids
#' @export
#' @examples
#' source(system.file("extdata", "sample_flines.R", package = "nhdplusTools"))
#'
#' tocomid <- get_tocomid(sample_flines)
#'
#' tocomid <- get_tocomid(sample_flines, return_dendritic = FALSE)
#'
get_tocomid <- function(x, return_dendritic = TRUE, missing = 0,
remove_coastal = TRUE, add = TRUE) {
x <- check_names(x, "get_tocomid", tolower = TRUE)
hy_g <- get_hyg(x, add, "comid")
x <- st_drop_geometry(x)
if(remove_coastal)
x <- filter_coastal(x)
x <- rename(x, id = "comid")
if("terminalpa" %in% names(x)) {
order <- data.frame(id = x$id)
x <- group_split(group_by(x, .data$terminalpa))
x <- bind_rows(lapply(x, add_toids, return_dendritic = return_dendritic))
x <- left_join(order, x, by = c("id"))
} else {
x <- add_toids(x, return_dendritic)
}
x <- rename(x, comid = "id", tocomid = "toid")
if(!is.na(missing) & missing != 0) {
repl <- x[["tocomid"]] == 0
x[["tocomid"]][repl] <- rep(missing, sum(repl))
}
if(add) {
if(!is.null(hy_g)) {
x <- sf::st_sf(left_join(x, hy_g, by = "comid"))
}
x[ , c("comid", "tocomid",
names(x)[!names(x) %in% c("comid", "tocomid")])]
} else {
as.data.frame(select(x, "comid", "tocomid"))
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.