Added functions for WB events

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