From bf49008a0d642c8ad2553fd16f437ef81b34a8fe Mon Sep 17 00:00:00 2001
From: David Blodgett <dblodgett@usgs.gov>
Date: Sat, 30 Nov 2024 21:31:22 -0600
Subject: [PATCH] waterbody outlet pois for #153
---
workspace/00_get_data.Rmd | 3 +-
workspace/R/02_POI_creation_functions.R | 750 ++++++++++++++++++-
workspace/R/NHD_navigate.R | 611 ---------------
workspace/_targets/02_POI_creation_targets.R | 21 +-
workspace/cache/data_paths.json | 1 +
5 files changed, 771 insertions(+), 615 deletions(-)
diff --git a/workspace/00_get_data.Rmd b/workspace/00_get_data.Rmd
index 9c1d2af..3e53781 100644
--- a/workspace/00_get_data.Rmd
+++ b/workspace/00_get_data.Rmd
@@ -26,7 +26,8 @@ if(!dir.exists("data")) {dir.create("data")}
if(!dir.exists("bin")) {dir.create("bin")}
data_dir <- "data"
-out_list <- list("data_dir" = data_dir)
+out_list <- list("data_dir" = data_dir,
+ nhdplushr_dir = file.path(data_dir, "nhdplushr"))
out_file <- file.path("cache", "data_paths.json")
sevenz <- "7z"
diff --git a/workspace/R/02_POI_creation_functions.R b/workspace/R/02_POI_creation_functions.R
index a1a371e..5de2321 100644
--- a/workspace/R/02_POI_creation_functions.R
+++ b/workspace/R/02_POI_creation_functions.R
@@ -379,4 +379,752 @@ create_hilarri_pois <- function(tmp_POIs, hilarri_data, wb_poi_lst, flowline, po
list(tmp_POIs = tmp_POIs, wb_pois = wb_poi_lst)
-}
\ No newline at end of file
+}
+
+#' create waterbody outlet pois
+#' @param tmp_POIs POIs from previous step
+#' @param NHDPlus_NHDArea NHDArea features from NHDPlusV2
+#' @param wb_poi_lst wb_pois from hilarri output list
+#' @param resops_wb_df wb_resops from resops output list
+#' @param flowline flowlines to attach pois to
+#' @param events events from gage output list
+#' @param nav_gpkg name of gpkg for current vpu
+#' @param data_paths list of data sources
+#' @param proj_crs projected crs for overall project
+#' @param gages gages from gage_info_gpkg
+#' @param poi_name name for poi set
+create_wb_outlet_pois <- function(tmp_POIs,
+ NHDPlus_NHDArea,
+ wb_poi_lst,
+ resops_wb_df,
+ flowline,
+ events,
+ nav_gpkg,
+ data_paths,
+ proj_crs,
+ gages,
+ poi_name) {
+
+ # Hillari IDs
+ hilarri_pois <- filter(tmp_POIs, !is.na(Type_hilarri))
+
+ # Get NHDARea polygons already define
+ nhd_area_vpu <- NHDPlus_NHDArea %>%
+ select(COMID, AREASQKM) %>%
+ mutate(source = "NHDv2Area") %>%
+ st_compatibalize(wb_poi_lst)
+
+ # Non-nhdv2 waterbodies from resopsus
+ other_wbs <- filter(resops_wb_df, !wb_id %in% wb_poi_lst$wb_id,
+ source != "REACH") %>%
+ select(grand_id, member_comid = wbareacomi, resops_flowlcomid = flowlcomid,
+ source, accounted) %>%
+ left_join(select(st_drop_geometry(hilarri_pois),
+ COMID, hilarriid = Type_hilarri),
+ by = c("resops_flowlcomid" = "COMID")) %>%
+ group_by(resops_flowlcomid) %>%
+ filter(row_number() == 1) %>%
+ ungroup() %>%
+ left_join(select(nhd_area_vpu, COMID) %>%
+ mutate(COMID = as.character(COMID)),
+ by = c("member_comid" = "COMID")) %>%
+ st_as_sf() %>%
+ mutate(wb_id = ifelse(source == "NHDAREA", member_comid, NA)) %>%
+ st_drop_geometry()
+
+ # copy empty geometries to wb_list for use in wb POI function
+ if(nrow(other_wbs) > 0){
+ wb_poi_list_fn <- data.table::rbindlist(list(wb_poi_lst,
+ other_wbs),
+ fill = TRUE) %>%
+ distinct() %>%
+ st_as_sf()
+ } else {
+ wb_poi_list_fn <- wb_poi_lst
+ }
+
+ # Final waterbody list with all requirements
+ final_wb_list <- filter(wb_poi_list_fn,
+ area_sqkm > wb_area_thresh |
+ !is.na(resops_flowlcomid) |
+ !is.na(hilarriid))
+
+ new_path <- "?prefix=StagedProducts/Hydrography/NHDPlusHR/VPU/Current/GDB/"
+ assign("nhdhr_file_list", new_path, envir = nhdplusTools:::nhdplusTools_env)
+
+ # fix something here
+ wb_layers <- wbout_POI_creaton(filter(flowline, poi == 1), final_wb_list,
+ data_paths, tmp_POIs, proj_crs, wb_poi_lst, poi_name)
+
+ wb_events <- poi_missing_wb <- wb_missing_poi <- double_wb <- double_wb_pois <- NULL
+
+ if(!is.data.frame(wb_layers)){
+ # write_sf(wb_layers$nhd_wb_table, temp_gpkg, "wb_flowpaths")
+ #
+ wbs <- wb_layers$WBs
+ wb_pois <- wb_layers$POIs %>% filter(!is.na(Type_WBOut))
+ missing_wb <- filter(wbs, !wb_id %in% wb_pois$Type_WBOut)
+ print(nrow(missing_wb))
+
+ double_wb <- wbs %>% group_by(wb_id) %>% filter(n() > 1)
+
+ double_wb_pois <- wb_pois %>%
+ group_by(Type_WBOut) %>%
+ filter(n() > 1)
+
+ wb_missing_poi <- filter(wbs,
+ !wb_id %in% wb_pois$Type_WBOut)
+
+ poi_missing_wb <- filter(wb_pois, !Type_WBOut %in% wb_layers$WBs$wb_id)
+
+ #*************************************
+ if(all(!is.na(wb_layers$events))){
+ wb_events <- wb_layers$events %>%
+ select(COMID, REACHCODE, REACH_meas, POI_identifier = wb_id,
+ event_type, nexus) %>%
+ st_compatibalize(., events)
+
+ events <- rbind(events, wb_events)
+ }
+
+ tmp_POIs <- wb_layers$POIs
+ # Write out waterbodies
+ WB_VPU_pois <- wb_layers$WBs
+ } else {
+ tmp_POIs <- wb_layers
+ WB_VPU_pois <- final_wb_list
+ }
+
+ #12047928
+ wb_out_col <- wb_poi_collapse(tmp_POIs, flowline, events, gages)
+ tmp_POIs <- wb_out_col$POIs
+
+ if(all(!is.na(wb_out_col$events_ret))){
+ events <- wb_out_col$events_ret
+ }
+
+ list(tmp_POIs = tmp_POIs, wb_out_col = wb_out_col,
+ all_events = events, wb_events = wb_events,
+ double_wb = double_wb, double_wb_pois = double_wb_pois,
+ wb_missing_poi = wb_missing_poi, poi_missing_wb = poi_missing_wb,
+ updated_flowline = wb_layers$nhd_wb_table)
+}
+
+#' Creates Waterbody POIs, calls a few other functions
+#' @param nhd (sf data.frame) flowline data.frame
+#' @param WBs_VPU (sf data.frame)
+#' @param data_paths (list) data paths to data layers
+#' @param crs (integer) CRS to use (epsg code)
+#'
+#' @return (sf data.frame) dataframe of waterbody outlet POIs
+#'
+wbout_POI_creaton <- function(nhd, wb_table, data_paths, tmp_POIs, crs, wb_poi_lst, poi_name){
+
+ wb_components_table <- st_drop_geometry(wb_table) %>%
+ dplyr::mutate(member_comid = strsplit(member_comid, ",")) %>%
+ tidyr::unnest(cols = member_comid) %>%
+ mutate(member_comid = as.integer(member_comid)) %>%
+ distinct()
+
+ # R17, some wierd stuff going on in the waterbodies
+ if("947070203" %in% wb_components_table$member_comid){
+ wb_components_table <- wb_components_table %>%
+ filter(!wb_id %in% c(1503847, 1513298))
+ }
+
+ # Bring wb_id to NHD data frame
+ if(("wb_id") %in% colnames(nhd)){
+ nhd <- select(nhd, -wb_id)
+ }
+
+ nhd <- nhd %>%
+ left_join(distinct(st_drop_geometry(wb_components_table), wb_id, member_comid) %>%
+ mutate(member_comid = as.numeric(member_comid)),
+ by = c("WBAREACOMI" = "member_comid"))
+
+
+ # Resops that are either NHDArea or MHD HR source and are missing in the
+ # nhd attributiong table
+ resops_wb_other_sources <- filter(wb_table, source ==
+ "HR ID AVAILABLE" | (source == "NHDAREA" &
+ !wb_id %in% nhd$wb_id) &
+ accounted == 0)
+
+ # Resops that are not accounted for and NHDv2 waterbody or NHDArea sourced
+ resops_main <- filter(wb_table, !is.na(resops_flowlcomid) &
+ !wb_id %in% resops_wb_other_sources$wb_id,
+ accounted == 0)
+
+ # Remaining WB that don't have outlet COMID defined and aren't in ResOPsUS
+ WB_outlet <- filter(wb_table, is.na(resops_flowlcomid),
+ !wb_id %in% tmp_POIs$Type_WBOut)
+
+ # Size-based waterbody outlet POIs ot in NHDArea, HR, or
+ # manually defined with ResOPS
+ wbout_COMIDs <- nhd %>%
+ filter(wb_id %in% WB_outlet$wb_id) %>%
+ group_by(wb_id) %>%
+ slice(which.min(Hydroseq)) %>%
+ ungroup() %>%
+ switchDiv(., nhd) %>%
+ select(COMID, wb_id) %>%
+ mutate(wbtype = "size") %>%
+ st_drop_geometry() %>%
+ distinct()
+
+ # ResOpsUS defined outlets
+ resops_main_outlets <- select(st_drop_geometry(resops_main), COMID = resops_flowlcomid,
+ wb_id) %>%
+ mutate(wbtype = "resops")
+
+ # combine together
+ wbout_coms <- distinct(rbind(wbout_COMIDs, resops_main_outlets))
+
+ # Get NHDArea and HI-Res waterbodies for ResOpsUS locaitons
+ WBs_VPU_areaHR <- HR_Area_coms(nhd, resops_wb_other_sources, data_paths, crs)
+
+ # if there are no HR/Area waterbodies generated
+ if(all(is.na(WBs_VPU_areaHR))){
+ tmp_POIs <- POI_creation(wbout_coms,
+ nhd, "wb_id") %>%
+ addType(., tmp_POIs, "wb_id", nexus = TRUE) %>%
+ mutate(Type_WBOut = ifelse(is.na(Type_WBOut), Type_wb_id, Type_WBOut)) %>%
+ select(-Type_wb_id) %>%
+ left_join(select(st_drop_geometry(resops_main), GRAND_ID, resops_FlowLcomid),
+ by = c("COMID" = "resops_FlowLcomid")) %>%
+ mutate(Type_resops = ifelse(!is.na(GRAND_ID), GRAND_ID, Type_resops)) %>%
+ select(-GRAND_ID)
+ } else {
+ # WBarea out comids
+ wbareaout_coms <- select(st_drop_geometry(nhd), COMID, Hydroseq) %>%
+ inner_join(filter(WBs_VPU_areaHR$nhd_wb_table, !is.na(outlet)),
+ by = c("COMID" = "comid")) %>%
+ group_by(wb_id) %>%
+ filter(Hydroseq == min(Hydroseq)) %>% #, row_number() == 1
+ ungroup() %>%
+ inner_join(select(WBs_VPU_areaHR$wb, -source), by = "wb_id") %>%
+ filter(source != "NHDHR") %>%
+ select(COMID, wb_id) %>%
+ mutate(wbtype = "resops")
+
+ if(nrow(wbareaout_coms) > 0){
+ wbout_coms <- rbind(wbout_coms, wbareaout_coms)
+ }
+
+ tmp_POIs <- POI_creation(select(wbout_coms, COMID, wb_id),
+ nhd, poi_name) %>%
+ addType(., tmp_POIs, poi_name, nexus = TRUE) %>%
+ mutate(Type_WBOut = ifelse(is.na(Type_WBOut) & (!nexus | is.na(nexus)),
+ Type_wb_id, Type_WBOut)) %>%
+ select(-Type_wb_id) %>%
+ left_join(select(st_drop_geometry(resops_main), grand_id, resops_flowlcomid),
+ by = c("COMID" = "resops_flowlcomid")) %>%
+ mutate(Type_resops = ifelse(!is.na(grand_id) & (!nexus | is.na(nexus)),
+ grand_id, Type_resops)) %>%
+ select(-grand_id)
+ }
+
+ # If the return is not NA
+ if(inherits(WBs_VPU_areaHR, "list")){
+ # Attribute flowlines that intersect NHDARea and/or NHD HR waterbodies with
+ # waterbody COMIDs
+ # For NHDHR, two waterbodies might be back to back and share an intersecting
+ # flowline, this may cause problems below, so deal with this case
+ WBs_VPU_areaHR$nhd_wb_table <- WBs_VPU_areaHR$nhd_wb_table %>%
+ left_join(select(st_drop_geometry(resops_wb_other_sources),
+ resops_flowlcomid, member_comid, grand_id),
+ by = c("comid" = "resops_flowlcomid",
+ "wb_comid" = "member_comid")) %>%
+ group_by(comid) %>%
+ filter(case_when(n() == 2 ~ !is.na(grand_id),
+ TRUE ~ n() == 1)) %>%
+ ungroup()
+
+ nhd <- nhd %>%
+ left_join(distinct(WBs_VPU_areaHR$nhd_wb_table, comid, wb_id_hr = wb_id),
+ by = c("COMID" = "comid")) %>%
+ mutate(wb_id = ifelse(!is.na(wb_id_hr), wb_id_hr, wb_id)) %>%
+ select(-wb_id_hr)
+
+ # NHD flowlines within the HR waterbodies
+ nhd_WBs <- filter(nhd, wb_id %in%
+ filter(WBs_VPU_areaHR$nhd_wb_table, source == "NHDHR")$wb_id)
+
+ if(nrow(nhd_WBs) > 0){
+ # Create outlet events for splitting HR and NHDArea waterbodies
+ wb_outlet_events <- WB_event(WBs_VPU_areaHR, nhd_WBs, "outlet") %>%
+ st_compatibalize(., tmp_POIs) %>%
+ mutate(nexus = TRUE, Type_resops = 2)
+
+ tmp_POIs <- data.table::rbindlist(list(tmp_POIs,
+ select(wb_outlet_events, COMID, Type_WBOut = wb_id, Type_resops, nexus)), fill = TRUE) %>%
+ st_as_sf()
+ } else {
+ wb_outlet_events <- NA
+ }
+
+ WBs_VPU_areaHR$wb <- WBs_VPU_areaHR$wb %>%
+ sfheaders::sf_remove_holes(.) %>%
+ st_compatibalize(wb_poi_lst)
+
+ WBs_VPU_areaHR$wb <- nhdplusTools::st_compatibalize(
+ WBs_VPU_areaHR$wb, wb_table)
+
+ wb_table_filtered <- wb_table %>%
+ filter(!st_is_empty(wb_table))
+
+ # Append NHDArea and NHDHR waterbodies to waterbody layer
+ WBs_VPU <- data.table::rbindlist(list(wb_table_filtered, WBs_VPU_areaHR$wb),
+ fill = TRUE) %>%
+ st_as_sf()
+
+ st_geometry(WBs_VPU) <- st_make_valid(st_geometry(WBs_VPU))
+
+ return(list(POIs = tmp_POIs, events = wb_outlet_events, WBs = WBs_VPU,
+ nhd_wb_table = nhd))
+
+ } else {
+ return(list(POIs = tmp_POIs, events = NA, WBs = wb_table,
+ nhd_wb_table = nhd))
+ }
+}
+
+### USEDIN wbout_POI_creaton and wbin_POIcreation
+#' Retrieves waterbodies from NHDArea layer or NHD Hi-Res for ResOpsUs
+#' locations if absent from NHD waterbody
+#' @param nhd (sf data.frame) Fowlines for given VPU
+#' @param WBs (sf data.frame) waterbodiess 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, crs){
+
+ # Pull out rows for VPU that are NHDArea
+ nhd_area_resops <- WBs %>%
+ filter(resops_flowlcomid %in% nhd$COMID, source == "NHDAREA")
+
+ # Pull out rows for VPU that are NHD HR
+ nhd_hr_wb_resops <- WBs %>%
+ filter(member_comid != 65000300139130) %>% # R09, in Canada
+ filter(resops_flowlcomid %in% nhd$COMID, source == "HR ID AVAILABLE") %>%
+ st_drop_geometry()
+
+ # 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% nhd_hr_wb_resops$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$member_comid) %>%
+ mutate(source = "NHDv2Area")
+
+ wb <- st_transform(nhd_area_vpu, crs)
+ }
+
+ # 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){
+ print(x)
+ vpu <- unique(substr(x, 1, 2))
+ nhd_file <- list.files(data_paths$nhdplushr_dir, pattern = paste0(".*",x,".*\\.gdb$"),
+ recursive = TRUE, full.names = TRUE,include.dirs = TRUE)
+
+ if(length(nhd_file) < 1) {
+ download_nhdplushr(data_paths$nhdplus_dir, unique(huc04))
+ }
+
+ nhd_file <- list.files(data_paths$nhdplushr_dir, pattern = paste0(".*",x,".*\\.gdb$"),
+ recursive = TRUE, full.names = TRUE,include.dirs = TRUE)
+
+ if(length(nhd_file) != 1) {
+ stop("something wrong with nhdplushr file in ", x)
+ }
+
+ # Format to similar to NHDArea/Waterbody layers
+ read_sf(nhd_file, "NHDWaterbody")})) #%>%
+
+ # convert to lower case and designate new shape field
+ names(hr_wb) <- tolower(names(hr_wb))
+ st_geometry(hr_wb) <- "shape"
+
+ hr_wb <- filter(hr_wb, permanent_identifier %in% WBs$member_comid) %>%
+ rename(GNIS_NAME = gnis_name, GNIS_ID = gnis_id,
+ Permanent_Identifier = permanent_identifier, AREASQKM = areasqkm,
+ FTYPE = ftype, FCODE = fcode, FDATE = fdate,
+ REACHCODE = reachcode) %>%
+ #select(-c(permanent_identifier, visibilityfilter, vpuid)) %>%
+ st_zm(.) %>%
+ st_as_sf() %>%
+ mutate(source = "NHDHR",
+ wb_id = ifelse(!is.na(GNIS_ID), GNIS_ID, Permanent_Identifier)) %>%
+ st_transform(crs)
+
+ # Bind or create new object
+ if(exists("wb")){
+ hr_wb <- st_compatibalize(hr_wb, 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 <- read.csv("data/reservoir_data/ISTARF-CONUS.csv") %>%
+ # filter(WBAREACOMI %in% wb$COMID | HR_NHDPLUSID %in% wb$COMID)
+
+ # Clear out columns not needed
+ WBs_sub <- WBs %>%
+ select(-any_of(c("COMID", "GNIS_ID")))
+
+ WBs_fin <- st_drop_geometry(WBs_sub) %>%
+ mutate(member_comid = member_comid) %>%
+ inner_join(select(wb, Permanent_Identifier, GNIS_ID, GNIS_NAME2 = GNIS_NAME),
+ by = c("member_comid" = "Permanent_Identifier")) %>%
+ mutate(wb_id = ifelse(GNIS_ID %in% c(NA, " "), member_comid, GNIS_ID),
+ GNIS_NAME = ifelse(is.na(GNIS_NAME), GNIS_NAME2, GNIS_NAME)) %>%
+ select(-c(GNIS_ID, GNIS_NAME2)) %>%
+ st_as_sf()
+
+ # Create table of all flowlines that intersect the waterbody
+ nhd_wb <- st_intersects(st_transform(nhd, st_crs(wb)), 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,]$Permanent_Identifier,
+ outlet = ifelse(comid %in% WBs_fin$resops_flowlcomid, "outlet", NA),
+ comid = as.integer(as.character(comid))) %>%
+ #left_join(select(wb_table, DAM_ID, DAM_NAME, resops_FlowLcomid),
+ # by = c("comid" = "resops_FlowLcomid")) %>%
+ left_join(select(st_drop_geometry(wb), wb_id, Permanent_Identifier, GNIS_ID,
+ source),
+ by = c("wb_comid" = "Permanent_Identifier"))
+
+ return(list(nhd_wb_table = nhd_wb_all, wb = WBs_fin))
+}
+
+### USEDIN wbout_POI_creaton and wbin_POIcreation
+#' 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$nhd_wb_table
+ WBs_layer <- WBs$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, Hydroseq %in% outlet_fl$DnHydroseq,
+ LevelPathI %in% outlet_fl$LevelPathI) %>%
+ rbind(outlet_fl) %>%
+ arrange(desc(Hydroseq)) %>%
+ group_by(LevelPathI) %>%
+ # union together
+ summarize(do_union = TRUE) %>%
+ st_cast("LINESTRING")
+
+ WBs_HR <- filter(WBs_layer, source == "HR ID AVAILABLE")
+
+ # 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(st_transform(nhd_wb, st_crs(WBs_HR)),
+ WBs_HR)) > 0,] %>%
+ group_by(wb_id) %>%
+ 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(wb_id) %>%
+ mutate(pnt_id = row_number()) %>%
+ ungroup()
+
+ outlet_pnts <- do.call("rbind", lapply(unique(WB_FL_pnts$wb_id), function(x){
+ fl <- filter(WB_FL_pnts, wb_id == x)
+ wb <- filter(WBs_HR, wb_id == x)
+
+ # Determine which points intersect waterbody
+ WB_outlet_pnts <- fl[lengths(st_intersects(fl, wb)) > 0,] %>%
+ st_drop_geometry() %>%
+ group_by(wb_id) %>%
+ mutate(within_wb_id = row_number()) %>%
+ filter(within_wb_id >= max(within_wb_id)) %>%
+ ungroup() %>%
+ select(wb_id, 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 <- fl %>%
+ inner_join(WB_outlet_pnts, by = "wb_id") %>%
+ select(wb_id, pnt_id, orig_pnt_id, within_wb_id) %>%
+ filter(pnt_id >= orig_pnt_id) %>%
+ group_by(wb_id) %>%
+ summarize(do_union = F) %>%
+ st_cast("LINESTRING") %>%
+ filter(wb_id %in% wb$wb_id)
+
+ outlet_pnt <- sf::st_intersection(outlet_FL, WBs_HR) %>%
+ st_cast("MULTIPOINT") %>%
+ st_cast("POINT") %>% # was point
+ group_by(wb_id) %>%
+ filter(row_number() == max(row_number(), na.rm = T)) %>%
+ ungroup() %>%
+ mutate(id = row_number()) %>%
+ filter(wb_id == wb_id.1)
+ }))
+
+ # Derive the events
+ if(exists("wb_events")){
+ hr_events <- get_flowline_index(st_transform(nhd_wb, st_crs(outlet_pnts)),
+ outlet_pnts) %>%
+ inner_join(select(st_drop_geometry(nhd_wb), COMID, wb_id), by = "COMID") %>%
+ filter(wb_id %in% WBs_HR$wb_id) %>%
+ mutate(event_type = type) %>%
+ cbind(select(outlet_pnts, geom)) %>%
+ st_as_sf() %>%
+ st_transform(st_crs(wb_events)) %>%
+ st_compatibalize(., wb_events)
+
+ wb_events <- rbind(wb_events, hr_events) %>%
+ distinct()
+ } else {
+ wb_events <- do.call("rbind", lapply(unique(outlet_pnts$wb_id), function(x){
+ outlet_pnt <- filter(outlet_pnts, wb_id == x)
+ get_flowline_index(nhd_wb, outlet_pnt)})) %>%
+ left_join(distinct(st_drop_geometry(outlet_pnts), resops_flowlcomid, wb_id),
+ by = c("COMID" = "resops_flowlcomid")) %>%
+ 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(
+ st_transform(start_pts, st_crs(WBs_layer)), WBs_layer)) == 0,] %>%
+ rbind(filter(nhd_wb, Hydroseq %in% .$DnHydroseq,
+ LevelPathI %in% .$LevelPathI)) %>%
+ arrange(desc(Hydroseq)) %>%
+ unique()
+
+ inlet_ls <- inlet_FL %>%
+ group_by(LevelPathI) %>%
+ st_cast("POINT") %>%
+ summarize(do_union = F) %>%
+ st_cast("LINESTRING") %>%
+ ungroup()
+
+ inlet_pnts_int <- sf::st_intersection(
+ st_transform(inlet_ls, st_crs(WBs_layer)), WBs_layer) #%>%
+ #st_cast("LINESTRING")
+
+ int_point <- inlet_pnts_int[st_geometry_type(inlet_pnts_int) %in% c("POINT"),]
+ if(nrow(int_point) > 0){
+ inlet_fl <- inlet_pnts_int[st_geometry_type(inlet_pnts_int) %in% c("LINESTRING", "MULTILINESTRING"),] %>%
+ st_cast("LINESTRING")
+ inlet_pnts_int <- rbind(int_point, inlet_fl)
+ }
+
+ inlet_pnts <- inlet_pnts_int %>%
+ group_by(LevelPathI) %>%
+ st_cast("POINT") %>%
+ mutate(id = row_number()) %>%
+ filter(id == min(id), !is.na(wb_id)) %>%
+ ungroup()
+
+ wb_events <- get_flowline_index(st_transform(nhd_wb, st_crs(inlet_pnts)),
+ inlet_pnts) %>%
+ inner_join(select(st_drop_geometry(nhd_wb), COMID, wb_id, LevelPathI), by = "COMID") %>%
+ mutate(event_type = type) %>%
+ inner_join(select(inlet_pnts, LevelPathI), by = "LevelPathI") %>%
+ #group_by(COMID, LEVELPATHI) %>%
+ mutate(nexus = T) %>%
+ rename(POI_identifier = wb_id) %>%
+ select(-c(id, offset)) %>%
+ st_as_sf()
+ }
+ return(wb_events)
+}
+
+### used as part of waterbody poi vreation in the workflow
+#' Collapses POIs us/ds of waterbody POIs
+#' @param (sf data.frame) rolling POI data frame
+#' @param nhd (sf data.frame) nhd flowlines
+#' @param events (sf data.frame) waterbody events
+#' @param gages gages from gage selection
+#'
+#' @return (sf data.frame) dataframe of wb inlet POIs collapsed
+wb_poi_collapse <- function(tmp_POIs, nhd, events, gages){
+ gage_dist_node <- function(x, wb_ds_ds, gage_add, events){
+ print (x)
+ wb_out_fl <- distinct(filter(wb_ds_ds, COMID == x))
+ gage_ds <- filter(wb_ds_ds, Hydroseq %in% wb_out_fl$Hydroseq |
+ Hydroseq %in% wb_out_fl$DnHydroseq)
+
+ gage_reach <- gage_ds %>%
+ group_by(REACHCODE) %>%
+ summarize(do_union = TRUE,
+ total_length = sum(LENGTHKM))
+
+ #print(nrow(gage_reach))
+ gage_event <- filter(gage_add, COMID %in% gage_ds$COMID) %>%
+ filter(reach_meas == min(reach_meas))
+ wb_event <- filter(events, COMID == x, event_type == "outlet") %>%
+ unique()
+
+ if(nrow(gage_reach) == 0){
+ print("no gages")
+ }
+
+ if(nrow(gage_event) == 0){
+ return("no events")
+ } else if(gage_event$COMID != wb_out_fl$COMID) {
+ gage_reach <- gage_reach %>%
+ filter(REACHCODE == unique(gage_event$reachcode)) %>%
+ mutate(gage_dist = ifelse(gage_event$nexus == TRUE,
+ total_length * (1 - (gage_event$reach_meas/100)),
+ total_length)) %>%
+ mutate(gage_comid = gage_event$COMID,
+ wbout_comid = x)
+ } else if(gage_event$COMID == wb_out_fl$COMID){
+ if(nrow(wb_event) >0){
+ wb_out_meas <- min(wb_event$REACH_meas)
+ wb_RC <- wb_event$REACHCODE
+ } else {
+ wb_out_meas <- min(wb_out_fl$FromMeas)
+ wb_RC <- wb_out_fl$REACHCODE
+ }
+
+ # gage reach
+ gage_reach <- gage_reach %>%
+ filter(REACHCODE == gage_event$reachcode) %>%
+ mutate(gage_dist = total_length * (1 - (gage_event$reach_meas/100)))
+
+ # wb info
+ wb_reach <- gage_reach %>%
+ filter(REACHCODE == unique(wb_RC)) %>%
+ mutate(wb_dist = total_length * (1 - (wb_out_meas /100)))
+
+ gage_reach <- gage_reach %>%
+ mutate(gage_dist = abs(wb_reach$wb_dist - gage_dist),
+ gage_comid = gage_event$COMID,
+ wbout_comid = x)
+ }
+ }
+
+ # Previously identified streamgages within Gage_Selection.Rmd
+ streamgages_VPU <- gages %>%
+ rename(COMID = comid) %>%
+ filter(COMID %in% nhd$COMID) %>%
+ switchDiv(., nhd)
+
+ # get waterbody outlets
+ wb_out <- filter(tmp_POIs, !is.na(Type_WBOut), is.na(Type_Gages))
+
+ wb_out_nexus <- filter(wb_out, nexus)
+ wb_out_nhd <- filter(nhd, COMID %in% wb_out$COMID)
+ wb_ds_ds <- wb_out_nhd %>%
+ rbind(filter(nhd, LevelPathI %in% .$LevelPathI & Hydroseq %in% .$DnHydroseq))
+
+ # get gages
+ gage_wb <- filter(tmp_POIs, !is.na(Type_Gages) & is.na(Type_WBOut)) #%>%
+ # filter(COMID %in% wb_ds_ds$COMID)
+ #gage_nhd <- filter(nhd, COMID %in% gage_wb$COMID)
+ #gage_node <- filter(gage_wb, nexus == FALSE)
+ #gage_nexus <- filter(gage_wb, nexus == TRUE)
+
+ gage_add <- filter(streamgages_VPU, site_no %in% gage_wb$Type_Gages) %>%
+ select(COMID, reachcode, reach_meas, site_no) %>%
+ inner_join(select(st_drop_geometry(gage_wb), site_no = Type_Gages, nexus),
+ by = "site_no") %>%
+ distinct()
+
+ # 8693865
+ output <- lapply(wb_out$COMID,
+ gage_dist_node, wb_ds_ds, gage_add, events)
+ output_length <- output[lengths(output) > 1]
+
+ if(length(output_length) == 0){
+ return(list(POIs = tmp_POIs, events_ret = NA))
+ }
+
+ output_full <- do.call("rbind", output[lengths(output) > 1]) %>%
+ filter(gage_dist < 1) %>%
+ left_join(select(st_drop_geometry(nhd), COMID, TotDASqKM_WB = TotDASqKM),
+ by = c("wbout_comid" = "COMID")) %>%
+ left_join(select(st_drop_geometry(nhd), COMID, TotDASqKM_gage = TotDASqKM),
+ by = c("gage_comid" = "COMID")) %>%
+ mutate(DAR = TotDASqKM_WB / TotDASqKM_gage) %>%
+ filter(gage_dist < 1, DAR > .975) %>%
+ st_drop_geometry()
+
+ gage_POI <- filter(tmp_POIs, COMID %in% output_full$gage_comid) %>%
+ select(COMID, Type_HUC12_ds = Type_HUC12, Type_Gages_ds = Type_Gages,
+ Type_TE_ds = Type_TE, nexus) %>%
+ st_drop_geometry() %>%
+ group_by(COMID) %>%
+ summarise(Type_HUC12_ds = last(na.omit(Type_HUC12_ds)),
+ Type_Gages_ds = last(na.omit(Type_Gages_ds)),
+ Type_TE_ds = last(na.omit(Type_TE_ds)),
+ nexus = last(na.omit(nexus)))
+
+ WB_POI <- filter(tmp_POIs, COMID %in% output_full$wbout_comid, !is.na(Type_WBOut)) %>%
+ inner_join(select(output_full, gage_comid, wbout_comid), by = c("COMID" = "wbout_comid")) %>%
+ inner_join(select(gage_POI, -nexus), by = c("gage_comid" = "COMID")) %>%
+ mutate(Type_HUC12 = ifelse(!is.na(Type_HUC12_ds), Type_HUC12_ds, Type_HUC12),
+ Type_Gages = ifelse(!is.na(Type_Gages_ds), Type_Gages_ds, Type_Gages),
+ Type_TE = ifelse(!is.na(Type_TE_ds), Type_TE_ds, Type_TE)) %>%
+ select(-c(Type_HUC12_ds, Type_Gages_ds, Type_TE_ds))
+
+ tmp_POIs_fin <- filter(tmp_POIs, !COMID %in% c(WB_POI$COMID, WB_POI$gage_comid)) %>%
+ rbind(select(WB_POI, -gage_comid))
+
+ if(any(gage_POI$nexus == TRUE)){
+ gage_events <- filter(gage_POI, nexus == TRUE)
+ events <- filter(events, !POI_identifier %in% gage_events$Type_Gages_ds)
+ }
+
+
+ return(list(POIs = tmp_POIs_fin, events_ret = events))
+}
diff --git a/workspace/R/NHD_navigate.R b/workspace/R/NHD_navigate.R
index 9d6764b..b00ffa9 100644
--- a/workspace/R/NHD_navigate.R
+++ b/workspace/R/NHD_navigate.R
@@ -508,477 +508,7 @@ split_elev <- function(in_POI_ID, split_DF, elev_diff, max_DA){
return(split_elev_DF)
}
-### USEDIN wbout_POI_creaton and wbin_POIcreation
-#' Retrieves waterbodies from NHDArea layer or NHD Hi-Res for ResOpsUs
-#' locations if absent from NHD waterbody
-#' @param nhd (sf data.frame) Fowlines for given VPU
-#' @param WBs (sf data.frame) waterbodiess 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, crs){
-
- # Pull out rows for VPU that are NHDArea
- nhd_area_resops <- WBs %>%
- filter(resops_flowlcomid %in% nhd$COMID, source == "NHDAREA")
-
- # Pull out rows for VPU that are NHD HR
- nhd_hr_wb_resops <- WBs %>%
- filter(member_comid != 65000300139130) %>% # R09, in Canada
- filter(resops_flowlcomid %in% nhd$COMID, source == "HR ID AVAILABLE") %>%
- st_drop_geometry()
-
- # 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% nhd_hr_wb_resops$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$member_comid) %>%
- mutate(source = "NHDv2Area")
-
- wb <- st_transform(nhd_area_vpu, crs)
- }
-
- # 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){
- print(x)
- vpu <- unique(substr(x, 1, 2))
-
- 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))
- }
-
- gdb <- list.files(file.path(data_paths$nhdplus_dir, substr(vpu, 1, 2)),
- pattern = paste0("_", x, "_.+gdb"), full.names = TRUE)
-
- # Format to similar to NHDArea/Waterbody layers
- read_sf(gdb, "NHDWaterbody")})) #%>%
-
- # convert to lower case and designate new shape field
- names(hr_wb) <- tolower(names(hr_wb))
- st_geometry(hr_wb) <- "shape"
-
- hr_wb <- filter(hr_wb, permanent_identifier %in% WBs$member_comid) %>%
- rename(GNIS_NAME = gnis_name, GNIS_ID = gnis_id,
- Permanent_Identifier = permanent_identifier, AREASQKM = areasqkm,
- FTYPE = ftype, FCODE = fcode, FDATE = fdate,
- REACHCODE = reachcode) %>%
- #select(-c(permanent_identifier, visibilityfilter, vpuid)) %>%
- st_zm(.) %>%
- st_as_sf() %>%
- mutate(source = "NHDHR",
- wb_id = ifelse(!is.na(GNIS_ID), GNIS_ID, Permanent_Identifier)) %>%
- st_transform(crs)
-
- # Bind or create new object
- if(exists("wb")){
- hr_wb <- st_compatibalize(hr_wb, 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 <- read.csv("data/reservoir_data/ISTARF-CONUS.csv") %>%
- # filter(WBAREACOMI %in% wb$COMID | HR_NHDPLUSID %in% wb$COMID)
-
- # Clear out columns not needed
- WBs_sub <- WBs %>%
- select(-any_of(c("COMID", "GNIS_ID")))
-
- WBs_fin <- st_drop_geometry(WBs_sub) %>%
- mutate(member_comid = member_comid) %>%
- inner_join(select(wb, Permanent_Identifier, GNIS_ID, GNIS_NAME2 = GNIS_NAME),
- by = c("member_comid" = "Permanent_Identifier")) %>%
- mutate(wb_id = ifelse(GNIS_ID %in% c(NA, " "), member_comid, GNIS_ID),
- GNIS_NAME = ifelse(is.na(GNIS_NAME), GNIS_NAME2, GNIS_NAME)) %>%
- select(-c(GNIS_ID, GNIS_NAME2)) %>%
- st_as_sf()
-
- # Create table of all flowlines that intersect the waterbody
- nhd_wb <- st_intersects(st_transform(nhd, st_crs(wb)), 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,]$Permanent_Identifier,
- outlet = ifelse(comid %in% WBs_fin$resops_flowlcomid, "outlet", NA),
- comid = as.integer(as.character(comid))) %>%
- #left_join(select(wb_table, DAM_ID, DAM_NAME, resops_FlowLcomid),
- # by = c("comid" = "resops_FlowLcomid")) %>%
- left_join(select(st_drop_geometry(wb), wb_id, Permanent_Identifier, GNIS_ID,
- source),
- by = c("wb_comid" = "Permanent_Identifier"))
-
- return(list(nhd_wb_table = nhd_wb_all, wb = WBs_fin))
-}
-
-### USEDIN wbout_POI_creaton and wbin_POIcreation
-#' 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$nhd_wb_table
- WBs_layer <- WBs$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, Hydroseq %in% outlet_fl$DnHydroseq,
- LevelPathI %in% outlet_fl$LevelPathI) %>%
- rbind(outlet_fl) %>%
- arrange(desc(Hydroseq)) %>%
- group_by(LevelPathI) %>%
- # union together
- summarize(do_union = TRUE) %>%
- st_cast("LINESTRING")
-
- WBs_HR <- filter(WBs_layer, source == "HR ID AVAILABLE")
-
- # 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(st_transform(nhd_wb, st_crs(WBs_HR)),
- WBs_HR)) > 0,] %>%
- group_by(wb_id) %>%
- 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(wb_id) %>%
- mutate(pnt_id = row_number()) %>%
- ungroup()
-
- outlet_pnts <- do.call("rbind", lapply(unique(WB_FL_pnts$wb_id), function(x){
- fl <- filter(WB_FL_pnts, wb_id == x)
- wb <- filter(WBs_HR, wb_id == x)
-
- # Determine which points intersect waterbody
- WB_outlet_pnts <- fl[lengths(st_intersects(fl, wb)) > 0,] %>%
- st_drop_geometry() %>%
- group_by(wb_id) %>%
- mutate(within_wb_id = row_number()) %>%
- filter(within_wb_id >= max(within_wb_id)) %>%
- ungroup() %>%
- select(wb_id, 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 <- fl %>%
- inner_join(WB_outlet_pnts, by = "wb_id") %>%
- select(wb_id, pnt_id, orig_pnt_id, within_wb_id) %>%
- filter(pnt_id >= orig_pnt_id) %>%
- group_by(wb_id) %>%
- summarize(do_union = F) %>%
- st_cast("LINESTRING") %>%
- filter(wb_id %in% wb$wb_id)
-
- outlet_pnt <- sf::st_intersection(outlet_FL, WBs_HR) %>%
- st_cast("MULTIPOINT") %>%
- st_cast("POINT") %>% # was point
- group_by(wb_id) %>%
- filter(row_number() == max(row_number(), na.rm = T)) %>%
- ungroup() %>%
- mutate(id = row_number()) %>%
- filter(wb_id == wb_id.1)
- }))
-
- # Derive the events
- if(exists("wb_events")){
- hr_events <- get_flowline_index(st_transform(nhd_wb, st_crs(outlet_pnts)),
- outlet_pnts) %>%
- inner_join(select(st_drop_geometry(nhd_wb), COMID, wb_id), by = "COMID") %>%
- filter(wb_id %in% WBs_HR$wb_id) %>%
- mutate(event_type = type) %>%
- cbind(select(outlet_pnts, geom)) %>%
- st_as_sf() %>%
- st_transform(st_crs(wb_events)) %>%
- st_compatibalize(., wb_events)
-
- wb_events <- rbind(wb_events, hr_events) %>%
- distinct()
- } else {
- wb_events <- do.call("rbind", lapply(unique(outlet_pnts$wb_id), function(x){
- outlet_pnt <- filter(outlet_pnts, wb_id == x)
- get_flowline_index(nhd_wb, outlet_pnt)})) %>%
- left_join(distinct(st_drop_geometry(outlet_pnts), resops_flowlcomid, wb_id),
- by = c("COMID" = "resops_flowlcomid")) %>%
- 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(
- st_transform(start_pts, st_crs(WBs_layer)), WBs_layer)) == 0,] %>%
- rbind(filter(nhd_wb, Hydroseq %in% .$DnHydroseq,
- LevelPathI %in% .$LevelPathI)) %>%
- arrange(desc(Hydroseq)) %>%
- unique()
-
- inlet_ls <- inlet_FL %>%
- group_by(LevelPathI) %>%
- st_cast("POINT") %>%
- summarize(do_union = F) %>%
- st_cast("LINESTRING") %>%
- ungroup()
-
- inlet_pnts_int <- sf::st_intersection(
- st_transform(inlet_ls, st_crs(WBs_layer)), WBs_layer) #%>%
- #st_cast("LINESTRING")
-
- int_point <- inlet_pnts_int[st_geometry_type(inlet_pnts_int) %in% c("POINT"),]
- if(nrow(int_point) > 0){
- inlet_fl <- inlet_pnts_int[st_geometry_type(inlet_pnts_int) %in% c("LINESTRING", "MULTILINESTRING"),] %>%
- st_cast("LINESTRING")
- inlet_pnts_int <- rbind(int_point, inlet_fl)
- }
-
- inlet_pnts <- inlet_pnts_int %>%
- group_by(LevelPathI) %>%
- st_cast("POINT") %>%
- mutate(id = row_number()) %>%
- filter(id == min(id), !is.na(wb_id)) %>%
- ungroup()
- wb_events <- get_flowline_index(st_transform(nhd_wb, st_crs(inlet_pnts)),
- inlet_pnts) %>%
- inner_join(select(st_drop_geometry(nhd_wb), COMID, wb_id, LevelPathI), by = "COMID") %>%
- mutate(event_type = type) %>%
- inner_join(select(inlet_pnts, LevelPathI), by = "LevelPathI") %>%
- #group_by(COMID, LEVELPATHI) %>%
- mutate(nexus = T) %>%
- rename(POI_identifier = wb_id) %>%
- select(-c(id, offset)) %>%
- st_as_sf()
- }
- return(wb_events)
-}
-
-#' Creates Waterbody POIs, calls a few other functions
-#' @param nhd (sf data.frame) flowline data.frame
-#' @param WBs_VPU (sf data.frame)
-#' @param data_paths (list) data paths to data layers
-#' @param crs (integer) CRS to use (epsg code)
-#'
-#' @return (sf data.frame) dataframe of waterbody outlet POIs
-#'
-wbout_POI_creaton <- function(nhd, wb_table, data_paths, tmp_POIs, crs){
-
- wb_components_table <- st_drop_geometry(wb_table) %>%
- dplyr::mutate(member_comid = strsplit(member_comid, ",")) %>%
- tidyr::unnest(cols = member_comid) %>%
- mutate(member_comid = as.integer(member_comid)) %>%
- distinct()
-
- # R17, some wierd stuff going on in the waterbodies
- if("947070203" %in% wb_components_table$member_comid){
- wb_components_table <- wb_components_table %>%
- filter(!wb_id %in% c(1503847, 1513298))
- }
-
- # Bring wb_id to NHD data frame
- if(("wb_id") %in% colnames(nhd)){
- nhd <- select(nhd, -wb_id)
- }
-
- nhd <- nhd %>%
- left_join(distinct(st_drop_geometry(wb_components_table), wb_id, member_comid) %>%
- mutate(member_comid = as.numeric(member_comid)),
- by = c("WBAREACOMI" = "member_comid"))
-
-
- # Resops that are either NHDArea or MHD HR source and are missing in the
- # nhd attributiong table
- resops_wb_other_sources <- filter(wb_table, source ==
- "HR ID AVAILABLE" | (source == "NHDAREA" &
- !wb_id %in% nhd$wb_id) &
- accounted == 0)
-
- # Resops that are not accounted for and NHDv2 waterbody or NHDArea sourced
- resops_main <- filter(wb_table, !is.na(resops_flowlcomid) &
- !wb_id %in% resops_wb_other_sources$wb_id,
- accounted == 0)
-
- # Remaining WB that don't have outlet COMID defined and aren't in ResOPsUS
- WB_outlet <- filter(wb_table, is.na(resops_flowlcomid),
- !wb_id %in% tmp_POIs$Type_WBOut)
-
- # Size-based waterbody outlet POIs ot in NHDArea, HR, or
- # manually defined with ResOPS
- wbout_COMIDs <- nhd %>%
- filter(wb_id %in% WB_outlet$wb_id) %>%
- group_by(wb_id) %>%
- slice(which.min(Hydroseq)) %>%
- ungroup() %>%
- switchDiv(., nhd) %>%
- select(COMID, wb_id) %>%
- mutate(wbtype = "size") %>%
- st_drop_geometry() %>%
- distinct()
-
- # ResOpsUS defined outlets
- resops_main_outlets <- select(st_drop_geometry(resops_main), COMID = resops_flowlcomid,
- wb_id) %>%
- mutate(wbtype = "resops")
-
- # combine together
- wbout_coms <- distinct(rbind(wbout_COMIDs, resops_main_outlets))
-
- # Get NHDArea and HI-Res waterbodies for ResOpsUS locaitons
- WBs_VPU_areaHR <- HR_Area_coms(nhd, resops_wb_other_sources, data_paths, crs)
-
- # if there are no HR/Area waterbodies generated
- if(all(is.na(WBs_VPU_areaHR))){
- tmp_POIs <- POI_creation(wbout_coms,
- nhd, "wb_id") %>%
- addType(., tmp_POIs, "wb_id", nexus = TRUE) %>%
- mutate(Type_WBOut = ifelse(is.na(Type_WBOut), Type_wb_id, Type_WBOut)) %>%
- select(-Type_wb_id) %>%
- left_join(select(st_drop_geometry(resops_main), GRAND_ID, resops_FlowLcomid),
- by = c("COMID" = "resops_FlowLcomid")) %>%
- mutate(Type_resops = ifelse(!is.na(GRAND_ID), GRAND_ID, Type_resops)) %>%
- select(-GRAND_ID)
- } else {
- # WBarea out comids
- wbareaout_coms <- select(st_drop_geometry(nhd), COMID, Hydroseq) %>%
- inner_join(filter(WBs_VPU_areaHR$nhd_wb_table, !is.na(outlet)),
- by = c("COMID" = "comid")) %>%
- group_by(wb_id) %>%
- filter(Hydroseq == min(Hydroseq)) %>% #, row_number() == 1
- ungroup() %>%
- inner_join(select(WBs_VPU_areaHR$wb, -source), by = "wb_id") %>%
- filter(source != "NHDHR") %>%
- select(COMID, wb_id) %>%
- mutate(wbtype = "resops")
-
- if(nrow(wbareaout_coms) > 0){
- wbout_coms <- rbind(wbout_coms, wbareaout_coms)
- }
-
- tmp_POIs <- POI_creation(select(wbout_coms, COMID, wb_id),
- nhd, "wb_id") %>%
- addType(., tmp_POIs, "wb_id", nexus = TRUE) %>%
- mutate(Type_WBOut = ifelse(is.na(Type_WBOut) & (!nexus | is.na(nexus)),
- Type_wb_id, Type_WBOut)) %>%
- select(-Type_wb_id) %>%
- left_join(select(st_drop_geometry(resops_main), grand_id, resops_flowlcomid),
- by = c("COMID" = "resops_flowlcomid")) %>%
- mutate(Type_resops = ifelse(!is.na(grand_id) & (!nexus | is.na(nexus)),
- grand_id, Type_resops)) %>%
- select(-grand_id)
- }
-
- # If the return is not NA
- if(inherits(WBs_VPU_areaHR, "list")){
- # Attribute flowlines that intersect NHDARea and/or NHD HR waterbodies with
- # waterbody COMIDs
- # For NHDHR, two waterbodies might be back to back and share an intersecting
- # flowline, this may cause problems below, so deal with this case
- WBs_VPU_areaHR$nhd_wb_table <- WBs_VPU_areaHR$nhd_wb_table %>%
- left_join(select(st_drop_geometry(resops_wb_other_sources),
- resops_flowlcomid, member_comid, grand_id),
- by = c("comid" = "resops_flowlcomid",
- "wb_comid" = "member_comid")) %>%
- group_by(comid) %>%
- filter(case_when(n() == 2 ~ !is.na(grand_id),
- TRUE ~ n() == 1)) %>%
- ungroup()
-
- nhd <- nhd %>%
- left_join(distinct(WBs_VPU_areaHR$nhd_wb_table, comid, wb_id_hr = wb_id),
- by = c("COMID" = "comid")) %>%
- mutate(wb_id = ifelse(!is.na(wb_id_hr), wb_id_hr, wb_id)) %>%
- select(-wb_id_hr)
-
- # NHD flowlines within the HR waterbodies
- nhd_WBs <- filter(nhd, wb_id %in%
- filter(WBs_VPU_areaHR$nhd_wb_table, source == "NHDHR")$wb_id)
-
- if(nrow(nhd_WBs) > 0){
- # Create outlet events for splitting HR and NHDArea waterbodies
- wb_outlet_events <- WB_event(WBs_VPU_areaHR, nhd_WBs, "outlet") %>%
- st_compatibalize(., tmp_POIs) %>%
- mutate(nexus = TRUE, Type_resops = 2)
-
- tmp_POIs <- data.table::rbindlist(list(tmp_POIs,
- select(wb_outlet_events, COMID, Type_WBOut = wb_id, Type_resops, nexus)), fill = TRUE) %>%
- st_as_sf()
- } else {
- wb_outlet_events <- NA
- }
-
- WBs_VPU_areaHR$wb <- WBs_VPU_areaHR$wb %>%
- sfheaders::sf_remove_holes(.) %>%
- st_compatibalize(wb_poi_lst)
-
- WBs_VPU_areaHR$wb <- nhdplusTools::st_compatibalize(
- WBs_VPU_areaHR$wb, wb_table)
-
- wb_table_filtered <- wb_table %>%
- filter(!st_is_empty(wb_table))
-
- # Append NHDArea and NHDHR waterbodies to waterbody layer
- WBs_VPU <- data.table::rbindlist(list(wb_table_filtered, WBs_VPU_areaHR$wb),
- fill = TRUE) %>%
- st_as_sf()
-
- st_geometry(WBs_VPU) <- st_make_valid(st_geometry(WBs_VPU))
-
- return(list(POIs = tmp_POIs, events = wb_outlet_events, WBs = WBs_VPU,
- nhd_wb_table = nhd))
-
- } else {
- return(list(POIs = tmp_POIs, events = NA, WBs = wb_table,
- nhd_wb_table = nhd))
- }
-}
#' Creates Waterbody POIs, calls a few other functions
#' @param nhd (sf data.frame) flowline data.frame
@@ -1286,147 +816,6 @@ wb_inlet_collapse <- function(tmp_POIs, nhd, events){
}
}
-### used as part of waterbody poi vreation in the workflow
-#' Collapses POIs us/ds of waterbody POIs
-#' @param (sf data.frame) rolling POI data frame
-#' @param nhd (sf data.frame) nhd flowlines
-#' @param events (sf data.frame) waterbody events
-#'
-#' @return (sf data.frame) dataframe of wb inlet POIs collapsed
-wb_poi_collapse <- function(tmp_POIs, nhd, events){
- gage_dist_node <- function(x, wb_ds_ds, gage_add, events){
- print (x)
- wb_out_fl <- distinct(filter(wb_ds_ds, COMID == x))
- gage_ds <- filter(wb_ds_ds, Hydroseq %in% wb_out_fl$Hydroseq |
- Hydroseq %in% wb_out_fl$DnHydroseq)
-
- gage_reach <- gage_ds %>%
- group_by(REACHCODE) %>%
- summarize(do_union = TRUE,
- total_length = sum(LENGTHKM))
-
- #print(nrow(gage_reach))
- gage_event <- filter(gage_add, COMID %in% gage_ds$COMID) %>%
- filter(reach_meas == min(reach_meas))
- wb_event <- filter(events, COMID == x, event_type == "outlet") %>%
- unique()
-
- if(nrow(gage_reach) == 0){
- print("no gages")
- }
-
- if(nrow(gage_event) == 0){
- return("no events")
- } else if(gage_event$COMID != wb_out_fl$COMID) {
- gage_reach <- gage_reach %>%
- filter(REACHCODE == unique(gage_event$reachcode)) %>%
- mutate(gage_dist = ifelse(gage_event$nexus == TRUE,
- total_length * (1 - (gage_event$reach_meas/100)),
- total_length)) %>%
- mutate(gage_comid = gage_event$COMID,
- wbout_comid = x)
- } else if(gage_event$COMID == wb_out_fl$COMID){
- if(nrow(wb_event) >0){
- wb_out_meas <- min(wb_event$REACH_meas)
- wb_RC <- wb_event$REACHCODE
- } else {
- wb_out_meas <- min(wb_out_fl$FromMeas)
- wb_RC <- wb_out_fl$REACHCODE
- }
-
- # gage reach
- gage_reach <- gage_reach %>%
- filter(REACHCODE == gage_event$reachcode) %>%
- mutate(gage_dist = total_length * (1 - (gage_event$reach_meas/100)))
-
- # wb info
- wb_reach <- gage_reach %>%
- filter(REACHCODE == unique(wb_RC)) %>%
- mutate(wb_dist = total_length * (1 - (wb_out_meas /100)))
-
- gage_reach <- gage_reach %>%
- mutate(gage_dist = abs(wb_reach$wb_dist - gage_dist),
- gage_comid = gage_event$COMID,
- wbout_comid = x)
- }
- }
-
- # Previously identified streamgages within Gage_Selection.Rmd
- streamgages_VPU <- gages %>%
- rename(COMID = comid) %>%
- filter(COMID %in% nhd$COMID) %>%
- switchDiv(., nhd)
-
- # get waterbody outlets
- wb_out <- filter(tmp_POIs, !is.na(Type_WBOut), is.na(Type_Gages))
-
- wb_out_nexus <- filter(wb_out, nexus)
- wb_out_nhd <- filter(nhd, COMID %in% wb_out$COMID)
- wb_ds_ds <- wb_out_nhd %>%
- rbind(filter(nhd, LevelPathI %in% .$LevelPathI & Hydroseq %in% .$DnHydroseq))
-
- # get gages
- gage_wb <- filter(tmp_POIs, !is.na(Type_Gages) & is.na(Type_WBOut)) #%>%
- # filter(COMID %in% wb_ds_ds$COMID)
- #gage_nhd <- filter(nhd, COMID %in% gage_wb$COMID)
- #gage_node <- filter(gage_wb, nexus == FALSE)
- #gage_nexus <- filter(gage_wb, nexus == TRUE)
-
- gage_add <- filter(streamgages_VPU, site_no %in% gage_wb$Type_Gages) %>%
- select(COMID, reachcode, reach_meas, site_no) %>%
- inner_join(select(st_drop_geometry(gage_wb), site_no = Type_Gages, nexus),
- by = "site_no") %>%
- distinct()
-
- # 8693865
- output <- lapply(wb_out$COMID,
- gage_dist_node, wb_ds_ds, gage_add, events)
- output_length <- output[lengths(output) > 1]
-
- if(length(output_length) == 0){
- return(list(POIs = tmp_POIs, events_ret = NA))
- }
-
- output_full <- do.call("rbind", output[lengths(output) > 1]) %>%
- filter(gage_dist < 1) %>%
- left_join(select(st_drop_geometry(nhd), COMID, TotDASqKM_WB = TotDASqKM),
- by = c("wbout_comid" = "COMID")) %>%
- left_join(select(st_drop_geometry(nhd), COMID, TotDASqKM_gage = TotDASqKM),
- by = c("gage_comid" = "COMID")) %>%
- mutate(DAR = TotDASqKM_WB / TotDASqKM_gage) %>%
- filter(gage_dist < 1, DAR > .975) %>%
- st_drop_geometry()
-
- gage_POI <- filter(tmp_POIs, COMID %in% output_full$gage_comid) %>%
- select(COMID, Type_HUC12_ds = Type_HUC12, Type_Gages_ds = Type_Gages,
- Type_TE_ds = Type_TE, nexus) %>%
- st_drop_geometry() %>%
- group_by(COMID) %>%
- summarise(Type_HUC12_ds = last(na.omit(Type_HUC12_ds)),
- Type_Gages_ds = last(na.omit(Type_Gages_ds)),
- Type_TE_ds = last(na.omit(Type_TE_ds)),
- nexus = last(na.omit(nexus)))
-
- WB_POI <- filter(tmp_POIs, COMID %in% output_full$wbout_comid, !is.na(Type_WBOut)) %>%
- inner_join(select(output_full, gage_comid, wbout_comid), by = c("COMID" = "wbout_comid")) %>%
- inner_join(select(gage_POI, -nexus), by = c("gage_comid" = "COMID")) %>%
- mutate(Type_HUC12 = ifelse(!is.na(Type_HUC12_ds), Type_HUC12_ds, Type_HUC12),
- Type_Gages = ifelse(!is.na(Type_Gages_ds), Type_Gages_ds, Type_Gages),
- Type_TE = ifelse(!is.na(Type_TE_ds), Type_TE_ds, Type_TE)) %>%
- select(-c(Type_HUC12_ds, Type_Gages_ds, Type_TE_ds))
-
- tmp_POIs_fin <- filter(tmp_POIs, !COMID %in% c(WB_POI$COMID, WB_POI$gage_comid)) %>%
- rbind(select(WB_POI, -gage_comid))
-
- if(any(gage_POI$nexus == TRUE)){
- gage_events <- filter(gage_POI, nexus == TRUE)
- events <- filter(events, !POI_identifier %in% gage_events$Type_Gages_ds)
- }
-
-
- return(list(POIs = tmp_POIs_fin, events_ret = events))
-}
-
### used when poi_move.R is sourced in POI collapse phase of workflow
### consider calling at the end of each poi creation process?
#' Collapses POIs together based on criteria
diff --git a/workspace/_targets/02_POI_creation_targets.R b/workspace/_targets/02_POI_creation_targets.R
index e7d511f..8a1d8f5 100644
--- a/workspace/_targets/02_POI_creation_targets.R
+++ b/workspace/_targets/02_POI_creation_targets.R
@@ -1,6 +1,6 @@
library(targets)
-tar_option_set(packages = c("dplyr", "sf", "hyfabric", "hydroloom"),
+tar_option_set(packages = c("dplyr", "sf", "hyfabric", "hydroloom", "nhdplusTools"),
memory = "transient", garbage_collection = TRUE, error = "null",
storage = "worker", retrieval = "worker", deployment = "main")
@@ -63,5 +63,22 @@ list(tar_target(data_paths_file, "cache/data_paths.json", format = "file"),
hilarri_data,
resops_pois$wb_pois,
flowline, "hilarri"),
- pattern = map(resops_pois, flowline), deployment = "worker")
+ pattern = map(resops_pois, flowline), deployment = "worker"),
+ # TODO: work on cleaning up parameters of this function
+ tar_target(wb_outlet_pois, create_wb_outlet_pois(hilarri_pois$tmp_POIs,
+ filter(combined_waterbodies, layer == "NHDArea"),
+ hilarri_pois$wb_pois,
+ resops_pois$wb_resops,
+ flowline,
+ gage_pois$events,
+ nav_gpkg,
+ data_paths,
+ proj_crs,
+ read_sf(gage_info_gpkg_file, "Gages"),
+ "wb_id"),
+ pattern = map(hilarri_pois, resops_pois, gage_pois,flowline, nav_gpkg),
+ deployment = "worker"),
+ tar_target(write_wb_flowline_mod, write_sf(wb_outlet_pois$updated_flowline,
+ nav_gpkg, nhd_flowline),
+ pattern = map(wb_outlet_pois, nav_gpkg))
)
\ No newline at end of file
diff --git a/workspace/cache/data_paths.json b/workspace/cache/data_paths.json
index f4eeb12..bfdbba1 100644
--- a/workspace/cache/data_paths.json
+++ b/workspace/cache/data_paths.json
@@ -1,5 +1,6 @@
{
"data_dir": "data",
+ "nhdplushr_dir": "data/nhdplushr",
"hu12_points_path": "data/102020wbd_outlets.gpkg",
"SWIM_points_path": "data/SWIM_gage_loc",
"SWIM_dbf": "data/SWIMGFinfo.dbf",
--
GitLab