From 975155e8058fc1ae194010405deb3b32d62c5989 Mon Sep 17 00:00:00 2001
From: Bock <abock@usgs.gov>
Date: Wed, 8 Jul 2020 14:49:09 -0600
Subject: [PATCH] Added structural POIs func., cleaned up others and added
 comments on input/output

---
 workspace/R/NHD_navigate.R | 316 +++++++++++++++++++++++++++----------
 1 file changed, 236 insertions(+), 80 deletions(-)

diff --git a/workspace/R/NHD_navigate.R b/workspace/R/NHD_navigate.R
index bba4fc8..72c961f 100644
--- a/workspace/R/NHD_navigate.R
+++ b/workspace/R/NHD_navigate.R
@@ -9,7 +9,7 @@ nhdflowline <- "NHDFlowline" # NHD flowlines pers VPU
 WBs <-  paste0("WB_", hydReg) # Waterbodies within VPU
 WBout_POIs <- paste0("WBOut_", hydReg) # Waterbody Outlet POIs
 WBin_POIs <- paste0("WBIn_", hydReg) # Waterbody Inlet POIs
-gagesIII_pois <- paste0("GagesIII_", hydReg) # GAGESIII POIs
+gages_pois <- paste0("Gages_", hydReg) # GAGESIII POIs
 gageLoc <- paste0("gageLoc_", hydReg) # GageLoc Files
 TE_pois <- paste0("TE_", hydReg) # Thermoelectric POIs
 NID_pois <- paste0("NID_", hydReg) # NID POIs
@@ -19,20 +19,38 @@ pois_all <- paste0("POIs_", hydReg) # All POIs binded together
 nsegment_raw <- paste0("nsegment_raw_", hydReg) # Minimally-sufficient network attributed with POI_ID
 n_segments <- paste0("Nsegment_", hydReg) # Minimally-sufficient network dissolved by POI_ID
 
-# Network navigation for upstream/downstream from a COMID of interest
+# 
 NetworkNav <- function(inCom, nhdDF, navType){
-  #print (inCom)
-  # Upstream network navigation is easier
+  ##########################################
+  # Network navigation for upstream/downstream from a COMID of interest
+  # 
+  # Arguments:
+  #   inCOM : (list) list of input COMIDs that are 'dangles'
+  #   nhdDF : (data frame) valid data frame of NHD flowlines
+  #
+  # Returns:
+  #   mydata : (list) list of COMIDs to connect dangle to network
+  ##########################################
   Seg <- nhdplusTools::get_UM(nhdDF, inCom, include = TRUE)
 
   return(Seg)
 }
 
-# Network navigation to connect dangles in the network
-NetworkConnection <- function(inCom, nhdDF){
-  # Subset by dangling segments that need downstream navigation to be connected
-  #        reduces the number of features we need to connect
-  upNet_DF <- nhdDF %>% dplyr::filter(COMID %in% inCom) %>% 
+
+NetworkConnection <- function(inCom, nhd){
+  ##########################################
+  # Connects dangles in the network that are not
+  #          terminal outlets
+  # 
+  # Arguments:
+  #   inCOM : (list) list of input COMIDs that are 'dangles'
+  #   nhdDF : (data frame) valid data frame of NHD flowlines
+  #
+  # Returns:
+  #   mydata : (list) list of COMIDs to connect dangle to network
+  ##########################################
+  # data frame of connections that need to be made
+  upNet_DF <- nhd %>% dplyr::filter(COMID %in% inCom) %>% 
     filter(!DnHydroseq %in% Hydroseq)
   
   # while the number of dangles is greater than 0
@@ -44,13 +62,13 @@ NetworkConnection <- function(inCom, nhdDF){
     DSLP <- upNet_DF$DnLevelPat[upNet_DF$COMID %in% inCom]
     # Get the COMID of the hydroseq with level path value
     # the lowest downstream flowline within the levelpath
-    inCom2 <- nhdDF$COMID[nhdDF$Hydroseq %in% DSLP]
+    inCom2 <- nhd$COMID[nhd$Hydroseq %in% DSLP]
     # Run the upstream navigation code
-    upNet <- unique(unlist(lapply(inCom2, NetworkNav, nhdDF, "up")))
+    upNet <- unique(unlist(lapply(inCom2, NetworkNav, nhd, "up")))
     # Append result to existing segment list
     inCom <- append(inCom, upNet)
     # Get the same variable as above
-    upNet_DF<-nhdDF %>% dplyr::filter(COMID %in% inCom) %>% 
+    upNet_DF<-nhd %>% dplyr::filter(COMID %in% inCom) %>% 
       filter(!DnHydroseq %in% Hydroseq)
     # Get the count
     count2 <- dim(upNet_DF)[1]
@@ -63,9 +81,19 @@ NetworkConnection <- function(inCom, nhdDF){
   return(inCom)
 }
 
-# Swith divergence path for POIs from minor (2) to major (1)
-# reduces the number of short, spurious segments and POIS as locations such as waterbody outlets
+
 switchDiv <- function(inSegs, nhd){
+  ##########################################
+  # Swith divergence path for POIs from minor (2) to major (1)
+  # reduces the number of short, spurious segments and POIS as locations such as waterbody outlets
+  # 
+  # Arguments:
+  #   inSegs : (list) list of input COMIDs that are 'dangles'
+  #   nhd : (data frame) valid data frame of NHD flowlines
+  #
+  # Returns:
+  #   mydata : (list) list of COMIDs for major diversions
+  ##########################################
   divSegs <- nhd %>% filter(COMID %in% inSegs$COMID)
   if (2 %in% divSegs$Divergence){
     print ("Switching divergence to other fork")
@@ -86,54 +114,84 @@ switchDiv <- function(inSegs, nhd){
   return(inSegs)
 }
 
-
-# POI Creation
-POI_creation<-function(inSegs, nhd){
-
-  # break segs into points
-  Segpts <- st_geometry(inSegs) %>% 
-    lapply(., function(x) {sf::st_sfc(x) %>% sf::st_cast(., 'POINT')})
+POI_creation<-function(srcData, nhd, IDfield){
+  ##########################################
+  # Creates POIs for a given data theme
+  # 
+  # Arguments:
+  #   srcData : (data frame) DF with two columns: 
+  #             1 - COMID
+  #             2 - Unique ID value for POI (Streamgage ID, etc.)
+  #   nhd : (data frame) valid data frame of NHD flowlines
+  #   IDfield : (character) Name of 'Type' field to be modified in POI 
+  #
+  # Returns:
+  #   mydata : (simple features) POIs for the specific data theme
+  ##########################################
   
-  # subset the last point from each geometry, make a POINT sf object
-  Seg_ends <- sapply(Segpts, function(p) {
-    p[c(length(p))]}) %>% 
-    sf::st_sfc(crs = st_crs(inSegs)) %>% 
-    sf::st_sf('geom' = .)
+  # Give generic ID to Identity field
+  colnames(srcData) <- c("COMID", "ID")
   
-  # Add the COMID, otherwise comes out as feature with no relating attributes
-  Seg_ends$COMID<-inSegs$COMID
-  
-  if("geom" %in% colnames(Seg_ends)) Seg_ends <- Seg_ends %>% rename(geometry=geom)
-  
-  return(Seg_ends)
+  subSegs <- nhd %>% filter(COMID %in% srcData$COMID) 
+  #merge(st_drop_geometry(subSegs %>% select(COMID))
+  POIs <-subSegs %>% get_node(., position = "end") %>% mutate(COMID = subSegs$COMID) %>%
+    mutate(Type_HUC12 = NA, Type_WBIn = NA, Type_WBOut = NA, Type_Gages = NA, Type_TE = NA, Type_NID = NA, Type_Conf = NA) %>%
+    inner_join(srcData %>% select(COMID, ID), by = "COMID") %>% mutate(!!(paste0("Type_", IDfield)) := ID) %>% 
+    select(COMID, Type_HUC12, Type_Gages, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf, geometry)
+
+  #if("geom" %in% colnames(Seg_ends)) Seg_ends <- Seg_ends %>% rename(geometry=geom)
+
+  return(POIs)
 }
 
-WB_area_width <- function(polygon){
+addType<-function(POIs, newPOIs, IDfield){
+  ##########################################
+  # Checks if existing POIs are co-located with new POI set
+  #        Adds the type attribute for co-located POIs of multiple themes
+  # 
+  # Arguments:
+  #   POIs : (data frame) Existing POIs
+  #   newPOIs: (data frame) newly-derived POIs of a given data theme
+  #   IDfield : (character) Name of 'Type' field to be modified in POI 
+  #
+  # Returns:
+  #   mydata : (data frame) Existing POIs with Type field modified for
+  #            duplicate POIs for given data theme
+  ##########################################
+  
+  POIs <- POIs %>% rename(ID = !!(paste0("Type_", IDfield)))
   
-  medWidth <- median(st_distance(st_cast(st_geometry(polygon),"POINT"))) 
+  POIs_fin <- POIs %>% left_join(filter(st_drop_geometry(newPOIs) %>% 
+                                 select(COMID, !!paste0("Type_", IDfield)), COMID %in% POIs$COMID), by = "COMID", all.x = TRUE) %>%
+    rename(ID2 = !!(paste0("Type_", IDfield))) %>%
+    mutate(ID = ifelse(!is.na(ID2), ID2, NA)) %>% rename(!!(paste0("Type_", IDfield)) := ID) %>% select(-c(ID2)) 
   
-  return (medWidth)
+  return(POIs_fin)
 }
 
-POI_move_down<-function(out_gpkg, pois, Seg, Seg2, exp){
+POI_move_down<-function(out_gpkg, POIs, Seg, Seg2, exp){
+  ##########################################
+  # Need to update & review this;
+  #      seems really long and could be shortened
+  ##########################################
   
   # Segs to collapse downstream (fold POIs together downstream if valid catchment downstream)
-  segs_Down <- Seg %>%  inner_join(select(st_drop_geometry(.),c(POI_ID, TotalDA, TotalLength)), 
-                                    by=c("To_POI_ID" = "POI_ID")) %>%
-    inner_join(select(st_drop_geometry(poisAll), 
-                      c(COMID, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID)),by=c("POI_ID"="COMID")) %>%
+  segs_Down <- Seg %>%  inner_join(select(st_drop_geometry(.), c(POI_ID, TotalDA, TotalLength)), 
+                                    by = c("To_POI_ID" = "POI_ID")) %>%
+    inner_join(select(st_drop_geometry(POIs), 
+                      c(COMID, Type_HUC12, Type_Gages, Type_Conf, Type_TE, Type_NID)), by = c("POI_ID" = "COMID")) %>%
     # Select POIs within the correct drainage area ratio and fit the size criteria
     mutate(DAR = TotalDA.x/TotalDA.y, IncDA = TotalDA.y - TotalDA.x) %>% 
     # If the drainage area ratio is within acceptable range
     filter(DAR < 1 & DAR >= 0.95 | TotalLength.y < 1, IncDA > 0)
 
   # Filter by POI Type
-  Types <- c("Type_HUC12", "Type_GagesIII", "Type_TE", "Type_Conf", "Type_NID", "Type_TE")
+  Types <- c("Type_HUC12", "Type_Gages", "Type_TE", "Type_Conf", "Type_NID", "Type_TE")
   Types <- Types[Types != exp]
   
   # Divide to segments of interest based on POI Type
   # Only move gages downstream if they are upstream of confluence
-  if (exp %in% c("Type_GagesIII", "Type_TE")){
+  if (exp %in% c("Type_Gages", "Type_TE")){
     ConfPOIs <- pois %>% filter(Type_Conf == 1)
     SegSub <- segs_Down %>% filter_at(Types, all_vars(is.na(.))) %>%
       filter(To_POI_ID %in% ConfPOIs$COMID) %>%
@@ -144,14 +202,14 @@ POI_move_down<-function(out_gpkg, pois, Seg, Seg2, exp){
   }
   
   #1 - POIs that need to be moved downstream
-  MoveDownPOI <- pois %>% filter(COMID %in% SegSub$POI_ID)
+  MoveDownPOI <- POIs %>% filter(COMID %in% SegSub$POI_ID)
   # POIs that are stationary and their infomration will be merged with upstream  POIs
-  stationaryPOI <- pois %>% filter(!COMID %in% SegSub$POI_ID)
+  stationaryPOI <- POIs %>% filter(!COMID %in% SegSub$POI_ID)
   # Final Set to be merged with that don't involve either sets above
-  FinalPOI <- pois %>% filter(!COMID %in% SegSub$POI_ID & !COMID %in% SegSub$To_POI_ID) %>%
+  FinalPOI <- POIs %>% filter(!COMID %in% SegSub$POI_ID & !COMID %in% SegSub$To_POI_ID) %>%
     mutate(merged_COMID = NA)
 
-  #2 - Join SegSub assignment to pois to bring over POI attributes
+  #2 - Join SegSub assignment to POIs to bring over POI attributes
   MoveDownPOI_withAtt <- MoveDownPOI %>%
     inner_join(st_drop_geometry(SegSub), by = c("COMID" = "POI_ID"), suffix = c(".x", ".y"))  %>%
     select(-c(geom, TotalLength.y, DAR, IncDA))
@@ -162,22 +220,24 @@ POI_move_down<-function(out_gpkg, pois, Seg, Seg2, exp){
                by = c("COMID" = "To_POI_ID"), suffix = c(".x", ".y")) %>%
     # Bring over relevant attributes
     mutate(Type_HUC12 = ifelse(!is.na(Type_HUC12.y), 1, Type_HUC12.x)) %>%
-    mutate(Type_GagesIII = ifelse(is.na(Type_GagesIII.x) & !is.na(Type_GagesIII.y), Type_GagesIII.y, NA)) %>%
+    mutate(Type_Gages = ifelse(is.na(Type_Gages.x) & !is.na(Type_Gages.y), Type_Gages.y, NA)) %>%
     # Gages_B is incase there are two gages being merged together, not writing out for now
-    mutate(Type_GagesIII_B = ifelse(!is.na(Type_GagesIII.y), Type_GagesIII.y, NA)) %>%
+    mutate(Type_Gages_B = ifelse(!is.na(Type_Gages.y), Type_Gages.y, NA)) %>%
     mutate(Type_Conf = ifelse(!is.na(Type_Conf.y), Type_Conf.y, Type_Conf.x)) %>%
     mutate(Type_TE = ifelse(!is.na(Type_TE.y), Type_TE.y, Type_TE.x)) %>%
     mutate(Type_NID = ifelse(!is.na(Type_NID.y), Type_NID.y, Type_NID.x)) %>%
     mutate(merged_COMID = COMID.y) %>%
-    select(COMID, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID, merged_COMID, geom)
+    select(COMID, Type_HUC12, Type_Gages, Type_Conf, Type_TE, Type_NID, merged_COMID, geom)
   
   # Write to regional geopackage
   write_sf(MergedPOIs, out_gpkg, paste0("MergedPOIs_", exp, "_", hydReg))
 
+  #***********************************************************************************************
+  # Maybe we don't need to this now?
   # Raw nsegments (original resolution of NHDPlus flowlines)
-  nseg_raw <- Seg2 %>% left_join(select(st_drop_geometry(MergedPOIs),c(COMID,  merged_COMID)), 
+  nseg_raw <- Seg2 %>% left_join(select(st_drop_geometry(MergedPOIs),c(COMID, merged_COMID)), 
                                   by = c("POI_ID" = "merged_COMID")) %>%
-    mutate(POI_ID = ifelse(!is.na(COMID.y), COMID.y, POI_ID)) %>% rename(COMID = COMID.x) %>% select (-c(COMID.y))
+    mutate(POI_ID = ifelse(!is.na(COMID.y), COMID.y, POI_ID)) %>% rename(COMID = COMID.x) %>% select(-c(COMID.y))
 
   # Write to regional geopackage
   write_sf(nseg_raw, out_gpkg, paste0("nsegment_raw_", hydReg))
@@ -194,8 +254,8 @@ POI_move_down<-function(out_gpkg, pois, Seg, Seg2, exp){
   
   # Add To_POI_ID to dissolved segments
   nsegments<-nsegments %>% left_join(to_from, by=c("POI_ID" = "COMID.x")) %>%
-    select(POI_ID,TotalLength,TotalDA,HW,TT,Hydroseq.x,To_POI_ID) %>%
-    rename(Hydroseq=Hydroseq.x)
+    select(POI_ID, TotalLength, TotalDA, HW, TT, Hydroseq.x, To_POI_ID) %>%
+    rename(Hydroseq = Hydroseq.x)
 
   # Write out dissolved segments
   write_sf(nsegments, out_gpkg, paste0("nsegment_", hydReg))
@@ -208,28 +268,42 @@ POI_move_down<-function(out_gpkg, pois, Seg, Seg2, exp){
 
 # Adjust confluence POIs based on if they are associated with flowline with no
 #        catchment/IncDA == 0
-DS_poiFix <- function(tmpPOIs, nhd){
+DS_poiFix <- function(POIs, nhd){
+  ##########################################
+  # Moves POI Upstream or downstream if it falls on COMID
+  #       of flowline with no corresponding catchment
+  # 
+  # Arguments:
+  #   POIs : (data frame) POI data set
+  #   nhd : (data frame) valid data frame of NHD flowlines
+  #
+  # Returns:
+  #   repPOIs_unique : (data frame) DF of POIs with new COMID associated
+  ##########################################
   
-  # Have subset nhd with no geometry
+  # Remove geom from NHD
   nhdDF <- st_drop_geometry(nhd)
+  # Include downstream HS into data frame
   ToComDF <- nhdDF %>% select(COMID, DnHydroseq) %>% 
     inner_join(st_drop_geometry(nhd) %>% select(COMID, Hydroseq), by = c("DnHydroseq" = "Hydroseq"))
   
   # Find segments with POIs where there is no corresponding catchment
-  unConPOIs <- nhdDF %>% filter(COMID %in% tmpPOIs$COMID, AreaSqKM == 0) 
-  poiFix <- tmpPOIs %>% filter(COMID %in% unConPOIs$COMID) %>%
-    mutate(New_COMID = unlist(lapply(.$COMID, movePOI_NA_DA, nhdDF)), old_COMID = COMID) 
+  unConPOIs <- nhdDF %>% filter(COMID %in% POIs$COMID, AreaSqKM == 0) 
+  # apply initial move function to the whole POI Set
+  poiFix <- POIs %>% filter(COMID %in% unConPOIs$COMID) %>%
+    mutate(New_COMID = unlist(lapply(.$COMID, movePOI_NA_DA, nhdDF)), old_COMID = COMID)
   
-  # Convert to POI data frame format with appropriate fields
-  repPOIs <- nhd %>% filter(COMID %in% poiFix$New_COMID) %>% POI_creation(., nhd) %>%
+  # Convert to POI data frame format with appropriate fields to match input POIs
+  repPOIs <- nhd %>% filter(COMID %in% poiFix$New_COMID) %>% get_node(., position = "end") %>% 
+    mutate(COMID = (nhd %>% filter(COMID %in% poiFix$New_COMID) %>% pull(COMID))) %>% 
     left_join(poiFix %>% filter(!is.na(New_COMID)), by = c("COMID" = "New_COMID")) %>% 
     select(-COMID.y) %>% st_drop_geometry()
   
   # Fold in new POIs with existing POIs so all the "Type" attribution will carry over
-  repPOIs_unique <- tmpPOIs %>% filter(COMID %in% repPOIs$COMID) %>% rbind(repPOIs %>% select(-old_COMID)) %>% 
+  repPOIs_unique <- POIs %>% filter(COMID %in% repPOIs$COMID) %>% rbind(repPOIs %>% select(-old_COMID)) %>% 
     group_by(COMID) %>% summarise_each(funs(max(., na.rm=T))) 
   
-  # Replace -Inf values
+  # Replace -Inf values that is an output of the 'summarise-each' application
   repPOIs_unique[mapply(is.infinite, repPOIs_unique)] <- NA
   
   # Get final list of POIs being moved, and the downstream POI they are linking to (many : 1 join)
@@ -269,23 +343,34 @@ DS_poiFix <- function(tmpPOIs, nhd){
   return (repPOIs_unique)
 }
 
-# Move POIs that fall on flowlines with no catchment upstream/downstream
-#  to adjacent flowline with most similar total drainage area
-movePOI_NA_DA <- function(poiFix, nhdDF){
-  # by most reandom to least random
+
+movePOI_NA_DA <- function(poiFix, nhd){
+  ##########################################
+  # Move POIs that fall on flowlines with no catchment upstream/downstream
+  #  to adjacent flowline with most similar total drainage area. Called from 
+  #  DS_poi_fix function above
+  # 
+  # Arguments:
+  #   poiFix : (data frame) POI data set of COMIDs to be changed
+  #   nhd : (data frame) valid data frame of NHD flowlines
+  #
+  # Returns:
+  #   newPOI : (data frame) DF of POIs with new COMID associated
+  ##########################################
   
   # Closest POI/US/DS
-  upSegs <- nhdplusTools::get_UM(nhdDF, poiFix, sort=T)
-  seg2fix <-nhdDF %>% filter(COMID == poiFix)
+  upSegs <- nhdplusTools::get_UM(nhd, poiFix, sort=T)
+  seg2fix <-nhd %>% filter(COMID == poiFix)
+  
   # Sorted results and filter out all flowlines w/o catchments
-  upStuff <- nhdDF %>% filter(COMID %in% unlist(upSegs)) %>% arrange(factor(COMID, levels=upSegs)) %>%
+  upStuff <- nhd %>% filter(COMID %in% unlist(upSegs)) %>% arrange(factor(COMID, levels=upSegs)) %>%
     filter(AreaSqKM > 0)
   
-  downSegs <- nhdplusTools::get_DM(nhdDF, poiFix, sort=T)
-  downStuff <- nhdDF %>% filter(COMID %in% unlist(downSegs)) %>% arrange(factor(COMID, levels=downSegs)) %>%
+  downSegs <- nhdplusTools::get_DM(nhd, poiFix, sort=T)
+  downStuff <- nhd %>% filter(COMID %in% unlist(downSegs)) %>% arrange(factor(COMID, levels=downSegs)) %>%
     filter(AreaSqKM > 0)
   
-  # combine into one dataframe
+  # combine into one dataframe, select up/downstream seg with least change in total drainage area
   nearFL <- rbind(upStuff %>% select(COMID, TotDASqKM, AreaSqKM) %>% slice(1), downStuff %>% filter(AreaSqKM > 0) %>%
                     select(COMID, TotDASqKM, AreaSqKM) %>% slice(1))
   
@@ -293,6 +378,7 @@ movePOI_NA_DA <- function(poiFix, nhdDF){
   if (sum(nearFL$AreaSqKM) > 0){
     newPOI <- nearFL$COMID[(which.min(abs(seg2fix$TotDASqKM - nearFL$TotDASqKM)))]
   } else {
+    # Remove POI if not catchment associated with flowlines upstream or downstream
     print (poiFix)
     print ("US and DS flowlines also have no catchment, removing POI")
     newPOI <- NA
@@ -300,26 +386,96 @@ movePOI_NA_DA <- function(poiFix, nhdDF){
   return(newPOI)
 }
 
-# Write out final POI datasets with information
-#       used to generate segments
-writePOIs <- function(inDF, out_gpkg, Type){
+
+writePOIs <- function(POIs, out_gpkg, Type){
+  ##########################################
+  # Write out final POI datasets with information
+  # 
+  # Arguments:
+  #   POIs : (data frame) POI data set 
+  #   out_gkpg : (geopackage) Geopackage where final POI layers written
+  #   Type : (character) Type of POI being written; default is write features for all types
+  #
+  # Returns:
+  #   finPOIs : (data frame) DF of final POIs 
+  ##########################################
+  
   print ("Writing out final POIs")
+  # If type is missing write out all flowlines
   if (missing(Type)){
+    print ("Writing out all POI Types")
     lyrs <- st_layers(out_gpkg)
     # get subcategory of POIs
-    POIs <- lyrs$name[lyrs$geomtype == "Point"]
-    for (poi in POIs[-(grep("POIs_*", POIs))]){
+    POInames <- lyrs$name[!is.na(lyrs$geomtype) & lyrs$geomtype== "Point"]
+    print (POInames)
+    for (poi in POInames){
       print (poi)
       poiType <- unlist(strsplit(poi, "_"))[1]
-      subPOIs <- inDF %>%
+      subPOIs <- POIs %>%
         filter(!is.na(!!as.name(paste0("Type_", poiType))))
       
       write_sf(subPOIs, out_gpkg, poi)
     }
   } else {
-    subPOIs <- inDF %>%
+    subPOIs <- POIs %>%
       filter(!is.na(!!as.name(paste0("Type_", Type))))
     write_sf(subPOIs, out_gpkg, poi)
   }
-  return(inDF)
+}
+
+
+structPOIsNet <- function(ncombined, nhd, finalPOIs, out_gpkg){
+  ##########################################
+  # Produce the structural POIs
+  # 
+  # Arguments:
+  #   ncombined : (data frame) final Segments 
+  #   nhd : (data frame) valid data frame of NHD flowlines
+  #   finalPOIs :  (data frame) final POIs
+  #   out_gkpg : (geopackage) Geopackage where final POI layers written
+  #   Type : (character) Type of POI being written; default is write features for all types
+  #
+  # Returns:
+  #   writes Structural POIs and segments to geopackage 
+  ##########################################
+  # Terminal outlets from the initial network
+  termOut <- nhd %>% filter(!Hydroseq %in% DnHydroseq & TerminalFl == 1 & !COMID %in% finalPOIs$COMID)
+  
+  # POIs that are also terminal outlets
+  outPOIs <- nhd %>% filter(COMID %in% finalPOIs$COMID & TerminalFl == 1)
+  
+  # Confluence POIs
+  confPOIs <- ncombined %>% filter(COMID %in% finalPOIs$COMID[finalPOIs$Type_Conf == 1])
+  
+  # Waterbody outlet POIs
+  WBPOIs <- ncombined %>% filter(COMID %in% finalPOIs$COMID[!is.na(finalPOIs$Type_WBOut) | !is.na(finalPOIs$Type_WBIn)])
+  
+  # Waterbody inlets POIs
+  strucFlines <- termOut %>% bind_rows(outPOIs, confPOIs, WBPOIs) %>% arrange(COMID)
+  strucPOIs <- get_node(strucFlines, position = "end") %>% mutate(COMID = strucFlines$COMID, LevelPathI = strucFlines$LevelPathI) %>% 
+    st_drop_geometry()
+  
+  # Add back in any levelpaths missing (shouldn't be much)
+  miss_LP <- ncombined %>% filter(COMID %in% finalPOIs$COMID) %>% filter(!LevelPathI %in% strucPOIs$LevelPathI)
+  
+  # Only bind if there are rows present
+  if (nrow(miss_LP) > 0){
+    LP_pois <- c(miss_LP$LevelPathI, strucPOIs$LevelPathI)
+  } else {
+    LP_pois <- strucPOIs$LevelPathI
+  }
+  
+  # final stuctural POIs, order by COMID to match with strucPOIs
+  finalStruct <- nhd %>% filter(Hydroseq %in% LP_pois) %>% arrange(COMID)
+  structPOIs <- get_node(finalStruct, position = "end") %>% mutate(COMID = finalStruct$COMID) 
+  
+  #  produce unique set of flowlines linking POIs
+  finalNet <- unique(unlist(lapply(unique(finalStruct$COMID), NetworkNav, st_drop_geometry(nhd), "up")))
+  
+  # Subset NHDPlusV2 flowlines to navigation results and write to shapefile
+  StructNet <- nhd %>% filter(COMID %in% finalNet & grepl(paste0("^", hydReg, ".*"), .data$VPUID)) 
+  
+  # Write out minimally-sufficient network
+  write_sf(structPOIs, out_gpkg, "struct_POIs2")
+  write_sf(StructNet, out_gpkg, "struct_Net2") 
 }
-- 
GitLab