From e844edc8435075a1098e4c11641ca659aa65b9c0 Mon Sep 17 00:00:00 2001
From: Bock <abock@usgs.gov>
Date: Thu, 30 Jun 2022 10:48:47 -0600
Subject: [PATCH] Added functions for WB events
---
workspace/R/NHD_navigate.R | 238 +++++++++++++++++++++++++++++++++++++
1 file changed, 238 insertions(+)
diff --git a/workspace/R/NHD_navigate.R b/workspace/R/NHD_navigate.R
index 9ff7e7c..2253a59 100644
--- a/workspace/R/NHD_navigate.R
+++ b/workspace/R/NHD_navigate.R
@@ -775,3 +775,241 @@ split_tt <- function(in_POI_ID, split_DF, tt_diff, max_DA){
}
+#' Retrieves waterbodies from NHDArea layer or NHD Hi-Res for ResOpsUs
+#' locations if absent from NHD waterbody
+#' @param nhd (sf data.frame) flowlines for given VPU
+#' @param WBs (sf data.frame) waterbodys for discretization within VPU
+#' @param data_paths (list) data paths to data layers
+#'
+#' @return (list) wbs - sf data frame for NHDArea and HR waterbodies
+#' wb_table - table of flowlines and outlet info for each
+#' feature in wb
+HR_Area_coms <- function(nhd, WBs, data_paths){
+ # Read in resops table
+ resops_wb_df <- read.csv(data_paths$resops_NID_CW)
+
+ # Pull out rows for VPU that are NHDArea
+ nhd_area_resops <- resops_wb_df %>%
+ filter(FlowLcomid %in% nhd$COMID, COMI_srclyr == "NHDAREA")
+
+ # Pull out rows for VPU that are NHD HR
+ nhd_hr_wb_resops <- resops_wb_df %>%
+ filter(FlowLcomid %in% nhd$COMID, COMI_srclyr == "HR ID AVAILABLE")
+
+ # Get reachcodes for waterbody outlets, so we have an idea of which
+ # NHD HR 4-digit geodatabase we may need to retrieve
+ RC <- filter(nhd, COMID %in% c(nhd_area_resops$FlowLcomid,
+ nhd_hr_wb_resops$FlowLcomid))$REACHCODE
+
+ # If no NHDArea or HR waterbodies needed return NULL
+ if (nrow(nhd_area_resops) == 0 & nrow(nhd_hr_wb_resops) == 0){
+ return(NA)
+ }
+
+ # If NHDArea feature needed retrieve from National GDB
+ if (nrow(nhd_area_resops) > 0){
+ nhd_area_vpu <- read_sf(data_paths$nhdplus_gdb, "NHDArea") %>%
+ filter(COMID %in% nhd_area_resops$NHDAREA) %>%
+ mutate(source = "NHDv2Area")
+
+ wb <- nhd_area_vpu
+ }
+
+ # If NHDHR feature needed
+ if (nrow(nhd_hr_wb_resops) > 0){
+ # HUC04 we need to download
+ huc04 <- substr(RC, 1, 4)
+
+ # Download NHD HR HUC04 if we dont' have it, other wise load and
+ # Bind NHDHR waterbodies into one layer
+ hr_wb <- do.call("rbind", lapply(unique(huc04), function(x){
+
+ if(!file.exists(file.path(data_paths$nhdplus_dir, vpu,
+ paste0("NHDPLUS_H_", x, "_HU4_GDB.gdb")))){
+ download_nhdplushr(data_paths$nhdplus_dir, unique(huc04))
+ }
+
+ # Format to similar to NHDArea/Waterbody layers
+ read_sf(file.path(data_paths$nhdplus_dir, substr(vpu, 1, 2),
+ paste0("NHDPLUS_H_", x, "_HU4_GDB.gdb")), "NHDWaterbody")})) %>%
+ filter(NHDPlusID %in% nhd_hr_wb_resops$HR_NHDPLUSID) %>%
+ rename(COMID = NHDPlusID, GNIS_NAME = GNIS_Name, RESOLUTION = Resolution, AREASQKM = AreaSqKm, ELEVATION = Elevation,
+ FTYPE = FType, FCODE = FCode, FDATE = FDate, REACHCODE = ReachCode) %>%
+ select(-c(Permanent_Identifier, VisibilityFilter, VPUID)) %>%
+ st_zm(.) %>%
+ st_as_sf() %>%
+ mutate(source = "NHDHR")
+
+ # Bind or create new object
+ if(exists("wb")){
+ wb <- data.table::rbindlist(list(wb, hr_wb), fill = TRUE) %>%
+ st_as_sf()
+ } else {
+ wb <- hr_wb
+ }
+ }
+
+ # get the outlt rows from the table
+ resops_outlet <- filter(resops_wb_df, NHDAREA %in% wb$COMID | HR_NHDPLUSID %in% wb$COMID)
+
+ # Create table of all flowlines that intersect the waterbody
+ nhd_wb <- st_intersects(nhd, wb)
+ comid <- nhd[lengths(nhd_wb) > 0,]$COMID
+ nhd_wb_all <- nhd_wb[lengths(nhd_wb) > 0] %>%
+ purrr::set_names(comid) %>%
+ stack() %>%
+ # Ind is the default name for the set_names
+ rename(comid = ind, nhd_wb = values) %>%
+ mutate(wb_comid = wb[nhd_wb,]$COMID,
+ outlet = ifelse(comid %in% resops_outlet$FlowLcomid, "outlet", NA),
+ comid = as.integer(as.character(comid))) %>%
+ left_join(select(resops_wb_df, DAM_ID, DAM_NAME, FlowLcomid), by = c("comid" = "FlowLcomid")) %>%
+ left_join(select(st_drop_geometry(wb), COMID, source), by = c("wb_comid" = "COMID"))
+
+ return(list(wb_table = nhd_wb_all, wb = wb))
+}
+
+#' Creates wb inlet and outlet events for splitting in hyRefactor
+#' for waterbodies derived from NHDArea and NHDHR waterbodies
+#' @param WBs (sf data.frame) return from HR_Area_coms
+#' @param nhd_wb (sf data.frame) flowlines that intersect waterbodies
+#' @param type (character) whether to derive inlet or outlet points
+#'
+#' @return (sf data.frame) dataframe of WB inlet and outlet points to split
+#'
+WB_event <- function(WBs, nhd_wb, type){
+ # split into features and table
+ WBs_table <- WBs_VPU_areaHR$wb_table
+ WBs_layer <- WBs_VPU_areaHR$wb
+
+ if (type == "outlet"){
+ # get the outlet comid from the ResOps Table
+ outlet_fl <- filter(nhd_wb, COMID %in% filter(WBs_table, outlet == "outlet")$comid)
+
+ # Get the downstream flowline for continuity
+ ds_fl <- filter(nhd_wb, DnHydroseq %in% outlet_fl$Hydroseq,
+ LevelPathI %in% outlet_fl$LevelPathI) %>%
+ rbind(outlet_fl) %>%
+ group_by(WBAREACOMI) %>%
+ # union together
+ summarize(do_union = T)
+
+ WBs_area <- filter(WBs_layer, source == "NHDv2Area")
+ WBs_HR <- filter(WBs_layer, source == "NHDHR")
+
+ # For NHDArea waterbodies (better congruity with th flowline st)
+ if (nrow(WBs_area) > 0){
+
+ # Intersect unioned FL with waterbody and get intersecting point
+ outlet_pnts <- sf::st_intersection(ds_fl, WBs_area) %>%
+ # Cast to point
+ st_cast("POINT") %>%
+ group_by(WBAREACOMI) %>%
+ # Get furthest downstream point of intersection
+ filter(row_number() == max(row_number(), na.rm = T)) %>%
+ ungroup()
+
+ # Derive event from point to use for splitting within hyRefactor
+ wb_events <- get_flowline_index(nhd_wb, outlet_pnts) %>%
+ inner_join(select(st_drop_geometry(nhd_wb), COMID, WBAREACOMI), by = "COMID") %>%
+ mutate(event_type = type) %>%
+ cbind(select(outlet_pnts, geom)) %>%
+ st_as_sf()
+ }
+
+ # For NHD HR waterbody outlets its a bit more complicated
+ if (nrow(WBs_HR) > 0){
+ # Get the flowlines intersecting the HR waterbody and find one with the
+ # max DA
+ outlet_wb_int <- nhd_wb[lengths(st_intersects(nhd_wb, WBs_HR)) > 0,] %>%
+ group_by(WBAREACOMI) %>%
+ filter(TotDASqKM == max(TotDASqKM)) %>%
+ ungroup()
+
+ # get the ds flo with the same levepath (JIC)
+ ds_fl <- filter(nhd_wb, DnHydroseq %in% outlet_wb_int$Hydroseq,
+ LevelPathI %in% outlet_wb_int$LevelPathI)
+
+ outlet_fl <- rbind(outlet_wb_int, ds_fl)
+
+ # Cast flowlines within NHDHR waterbody to point
+ WB_FL_pnts <- outlet_wb_int %>%
+ st_cast("POINT") %>%
+ group_by(WBAREACOMI) %>%
+ mutate(pnt_id = row_number()) %>%
+ ungroup()
+
+ # Determine which points intersect waterbody
+ WB_outlet_pnts <- WB_FL_pnts[lengths(st_intersects(WB_FL_pnts, WBs_HR)) > 0,] %>%
+ st_drop_geometry() %>%
+ group_by(WBAREACOMI) %>%
+ mutate(within_wb_id = row_number()) %>%
+ filter(within_wb_id >= max(within_wb_id)) %>%
+ ungroup() %>%
+ select(WBAREACOMI, orig_pnt_id = pnt_id, within_wb_id)
+
+ # Deriv new linestring by concating points from most upstream point
+ # within waterbody to downstream so we can split at FL/waterbody
+ # nexus
+ outlet_FL <- WB_FL_pnts %>%
+ inner_join(WB_outlet_pnts, by = "WBAREACOMI") %>%
+ select(WBAREACOMI, pnt_id, orig_pnt_id, within_wb_id) %>%
+ filter(pnt_id >= orig_pnt_id) %>%
+ group_by(WBAREACOMI) %>%
+ summarize(do_union = F) %>%
+ st_cast("LINESTRING")
+
+ # now run the intersection with the modified flowline
+ outlet_pnts <- sf::st_intersection(outlet_FL, WBs_HR) %>%
+ st_cast("POINT") %>%
+ group_by(WBAREACOMI) %>%
+ filter(row_number() == max(row_number(), na.rm = T)) %>%
+ ungroup()
+
+ # Derive the events
+ if(exists("wb_events")){
+ hr_events <- get_flowline_index(nhd_wb, outlet_pnts) %>%
+ inner_join(select(st_drop_geometry(nhd_wb), COMID, WBAREACOMI), by = "COMID") %>%
+ mutate(event_type = type) %>%
+ cbind(select(outlet_pnts, geom)) %>%
+ st_as_sf()
+
+ wb_events <- rbind(wb_events, hr_events)
+ } else {
+ wb_events <- get_flowline_index(nhd_wb, outlet_pnts) %>%
+ inner_join(select(st_drop_geometry(nhd_wb), COMID, WBAREACOMI), by = "COMID") %>%
+ mutate(event_type = type) %>%
+ cbind(select(outlet_pnts, geom)) %>%
+ st_as_sf()
+ }
+
+ }
+
+ # For inlet points its alot easier for both NHDARea and NHDHR
+ } else {
+ start_pts <- get_node(nhd_wb, position = "start") %>%
+ cbind(st_drop_geometry(nhd_wb))
+
+ inlet_FL <- nhd_wb[lengths(st_intersects(start_pts, WBs_layer)) == 0,] %>%
+ rbind(filter(nhd_wb, Hydroseq %in% .$DnHydroseq,
+ LevelPathI %in% .$LevelPathI))
+
+ inlet_ls <- inlet_FL %>%
+ group_by(LevelPathI) %>%
+ summarize(do_union = T)
+
+ inlet_pnts <- sf::st_intersection(inlet_ls, WBs_layer) %>%
+ st_cast("POINT") %>%
+ group_by(LevelPathI) %>%
+ filter(row_number() == min(row_number())) %>%
+ ungroup()
+
+ wb_events <- get_flowline_index(nhd_wb, inlet_pnts) %>%
+ inner_join(select(st_drop_geometry(nhd_wb), COMID, WBAREACOMI, LevelPathI), by = "COMID") %>%
+ mutate(event_type = type) %>%
+ inner_join(select(inlet_pnts, LevelPathI), by = "LevelPathI") %>%
+ select(-LevelPathI) %>%
+ st_as_sf()
+ }
+ return(wb_events)
+}
\ No newline at end of file
--
GitLab