diff --git a/workspace/00_get_data.Rmd b/workspace/00_get_data.Rmd
index 241ef2b9ac46cc3bb8585746f688aed60fbac8e1..cdeef4f85971abec57863c3b5022835c7196a3dd 100644
--- a/workspace/00_get_data.Rmd
+++ b/workspace/00_get_data.Rmd
@@ -213,9 +213,7 @@ if(!file.exists(nhdplus_dir)) {
, " -o", data_dir), ignore.stdout = T)}
)
- nhdplus_path(nhdplus_gdb)
-
- suppressWarnings(staged_nhd <- stage_national_data())
+ suppressWarnings(staged_nhd <- stage_national_data(nhdplus_data = nhdplus_gdb))
x <- tryCatch(
download_nhdplusv2(islands_dir, url = paste0("https://dmap-data-commons-ow.s3.amazonaws.com/NHDPlusV21/",
@@ -227,9 +225,7 @@ if(!file.exists(nhdplus_dir)) {
, " -o", islands_dir), ignore.stdout = T)}
)
- nhdplus_path(islands_gdb)
-
- suppressWarnings(staged_nhdplus_islands <- stage_national_data())
+ suppressWarnings(staged_nhdplus_islands <- stage_national_data(nhdplus_data = islands_gdb))
rm(gLz, xWalk)
@@ -304,8 +300,6 @@ waterbodies_path <- file.path(nhdplus_dir, "nhdplus_waterbodies.rds")
if(!file.exists(waterbodies_path)) {
message("formatting NHDPlus watebodies...")
- nhdplus_path(nhdplus_gdb)
-
# Read the feature class
wbSF <- read_sf(nhdplus_gdb, "NHDWaterbody") %>%
st_transform(5070)
diff --git a/workspace/01_Gage_Selection.Rmd b/workspace/01_Gage_Selection.Rmd
index 8a9a6ecc6c84c10d8129b9b466c295ab4fb0e7f9..30546dfa191eb278ad7b1ed6e6727da40e4e99b2 100644
--- a/workspace/01_Gage_Selection.Rmd
+++ b/workspace/01_Gage_Selection.Rmd
@@ -36,7 +36,7 @@ source("R/gage_mod.R")
# Load user-defined data paths
data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
# Load NHDPlus Data if precprocessed to RDS format
-nhdplus_path(data_paths$nhdplus_gdb)
+
staged_nhd <- stage_national_data(nhdplus_data = data_paths$nhdplus_gdb,
output_path = data_paths$nhdplus_dir)
diff --git a/workspace/R/config.R b/workspace/R/config.R
index 4d06f0ce56da7096418a4a79e1b858c7f42a900a..014d0127349ca49d6710ad3404b3cdf041420252 100644
--- a/workspace/R/config.R
+++ b/workspace/R/config.R
@@ -54,8 +54,7 @@ options(scipen = 9999)
crs <- 5070
data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
-nhdplus_path(data_paths$nhdplus_gdb)
-suppressWarnings(staged_nhd <- stage_national_data())
+suppressWarnings(staged_nhd <- stage_national_data(nhdplus_data = data_paths$nhdplus_gdb))
# Defined and used broadly
nhd_flowline <- "reference_flowline"
diff --git a/workspace/R/utils.R b/workspace/R/utils.R
index e7727d973be4676a3fa4ade12c5d182f739e5be5..0c8535912efc3143b8d2b42b84f695c7a056f41f 100644
--- a/workspace/R/utils.R
+++ b/workspace/R/utils.R
@@ -10,117 +10,6 @@ if(!(require(hyfabric, quietly = TRUE) && packageVersion("hyfabric") == hyfabric
repos = NULL, type = "source")
}
-fix_headwaters <- function(nhd_rds, out_gpkg, new_atts = NULL,
- nhdpath = nhdplusTools::nhdplus_path(),
- par = 6) {
-
- if(!file.exists(out_gpkg)) {
- nhd <- readRDS(nhd_rds)
-
- # Need to replace headwater flowlines with the trimmed burn line event.
- # Some headwater flowlines are not fully within the catchment they go with.
- ble <- sf::read_sf(nhdpath, "BurnLineEvent")
-
- ble <- dplyr::left_join(dplyr::select(sf::st_drop_geometry(nhd), COMID),
- ble, by = "COMID")
-
- ble <- sf::st_zm(sf::st_as_sf(ble))
- nhd <- sf::st_zm(nhd)
-
- sf::st_geometry(nhd)[!sf::st_is_empty(sf::st_geometry(ble)) & (nhd$StartFlag == 1)] <-
- sf::st_geometry(ble)[!sf::st_is_empty(sf::st_geometry(ble)) & (nhd$StartFlag == 1)]
-
- if(!is.null(new_atts)) {
- new_atts <- data.table::fread(new_atts,
- integer64 = "character",
- showProgress = FALSE)
-
- if("weight" %in% names(new_atts)) new_atts <- rename(new_atts,
- ArbolateSu = weight)
-
- nhd <- select(nhd, COMID, FromNode, ToNode,
- StartFlag, StreamCalc, Divergence, DnMinorHyd) %>%
- left_join(new_atts, by = c("COMID" = "comid")) %>%
- nhdplusTools::align_nhdplus_names()
- }
-
- nhd <- sf::st_sf(nhd)
-
- m_to_km <- (1/1000)
-
- nhd$LENGTHKM <- as.numeric(sf::st_length(sf::st_transform(nhd, 5070))) * m_to_km
-
- custom_net <- select(sf::st_drop_geometry(nhd),
- COMID, FromNode, ToNode, Divergence)
-
- custom_net <- nhdplusTools::get_tocomid(custom_net, remove_coastal = FALSE)
-
- custom_net <- select(custom_net, comid, override_tocomid = tocomid)
-
- nhd <- left_join(nhd, custom_net, by = c("COMID" = "comid"))
-
- nhd <- mutate(nhd, override_tocomid = ifelse(toCOMID == 0, override_tocomid, toCOMID))
-
- # is a headwater and for sure flows to something.
- check <- !nhd$COMID %in% nhd$override_tocomid &
- !(nhd$override_tocomid == 0 | is.na(nhd$override_tocomid) |
- !nhd$override_tocomid %in% nhd$COMID)
- check_direction <- dplyr::filter(nhd, check)
-
- if(!all(check_direction$override_tocomid[check_direction$override_tocomid != 0] %in% nhd$COMID))
- stop("this won't work")
-
- matcher <- match(check_direction$override_tocomid, nhd$COMID)
-
- matcher <- nhd[matcher, ]
-
- fn_list <- Map(list,
- split(check_direction, seq_len(nrow(check_direction))),
- split(matcher, seq_len(nrow(matcher))))
-
- cl <- parallel::makeCluster(16)
-
- # check and fix these.
- new_geom <- pbapply::pblapply(cl = cl,
- X = fn_list,
- FUN = function(x) {
- nhdplusTools::fix_flowdir(x[[1]]$COMID,
- fn_list = list(flowline = x[[1]],
- network = x[[2]],
- check_end = "end"))
- })
-
- parallel::stopCluster(cl)
-
- # One is empty. Deal with it and remove from the replacement.
- error_index <- sapply(new_geom, inherits, what = "try-error")
- errors <- filter(nhd, COMID %in% check_direction$COMID[error_index])
- check[which(nhd$COMID %in% errors$COMID)] <- FALSE
-
- if(!all(sapply(sf::st_geometry(errors), sf::st_is_empty))) {
- stop("Errors other than empty geometry from fix flowdir")
- }
-
- new_geom <- new_geom[!error_index]
- new_geom <- do.call(c, new_geom)
- st_geometry(nhd)[check] <- new_geom
-
- nhd <- dplyr::filter(nhd, COMID %in% new_atts$comid)
-
- nhd <- select(nhd, -override_tocomid)
-
- # remove remnant attributes
- nhd <- select(nhd, -FromNode, -ToNode, -StartFlag -StreamCalc, -Divergence, -DnMinorHyd)
-
- warning("naive duplication!")
- nhd <- nhd[!duplicated(nhd$COMID), ]
-
- sf::write_sf(nhd, out_gpkg, "reference_flowlines")
-
- }
- return(out_gpkg)
-}
-
## Is available in sbtools now, but will leave till on cran.
#' Retrieves files from SB in facet file_structure
#' @param sbitem (character) ScienceBase item ID
@@ -595,3 +484,78 @@ generate_lookup_table = function(nl,
}
+# deprecated function from nhdplusTools -- should remove from package but here for compatibility
+stage_national_data <- function(include = c("attribute",
+ "flowline",
+ "catchment"),
+ output_path = NULL,
+ nhdplus_data = NULL,
+ simplified = TRUE) {
+
+ if (is.null(output_path)) {
+ output_path <- dirname(nhdplus_data)
+ warning(paste("No output path provided, using:", output_path))
+ }
+
+ if (is.null(nhdplus_data)) {
+ stop("must provide nhdplus_data")
+ }
+
+ allow_include <- c("attribute", "flowline", "catchment")
+
+ if (!all(include %in% allow_include)) {
+ stop(paste0("Got invalid include entries. Expect one or more of: ",
+ paste(allow_include, collapse = ", "), "."))
+ }
+
+ outlist <- list()
+
+ if (any(c("attribute", "flowline") %in% include)) {
+
+ out_path_attributes <- file.path(output_path,
+ "nhdplus_flowline_attributes.rds")
+ out_path_flines <- file.path(output_path, "nhdplus_flowline.rds")
+
+ if (!(file.exists(out_path_flines) | file.exists(out_path_attributes))) {
+ fline <- sf::st_zm(sf::read_sf(nhdplus_data,
+ get_flowline_layer_name(nhdplus_data)))
+ }
+
+ if ("attribute" %in% include) {
+ if (file.exists(out_path_attributes)) {
+ warning("attributes file exists")
+ } else {
+ saveRDS(sf::st_set_geometry(fline, NULL), out_path_attributes)
+ }
+ outlist["attributes"] <- out_path_attributes
+ }
+
+ if ("flowline" %in% include) {
+ if (file.exists(out_path_flines)) {
+ warning("flowline file exists")
+ } else {
+ saveRDS(fline, out_path_flines)
+ }
+ outlist["flowline"] <- out_path_flines
+ }
+ }
+
+ if (exists("fline")) rm(fline)
+
+ if ("catchment" %in% include) {
+ out_path_catchments <- file.path(output_path, "nhdplus_catchment.rds")
+ if (file.exists(out_path_catchments)) {
+ warning("catchment already exists.")
+ } else {
+
+ layer_name <- get_catchment_layer_name(simplified, nhdplus_data)
+
+ saveRDS(sf::st_zm(sf::read_sf(nhdplus_data, layer_name)),
+ out_path_catchments)
+ }
+ outlist["catchment"] <- out_path_catchments
+ }
+ assign("national_data", outlist, envir = nhdplusTools_env)
+
+ return(outlist)
+}
\ No newline at end of file