NHD_navigate.R

# NHDPlus V2 Region
if(!exists("hydReg") && !nchar(hydReg <- Sys.getenv("HYDREG")) > 1) 
  hydReg <- "01"

# Layers
huc12_pois <- paste0("HUC12_", hydReg) # HUC12 POIs
dup_comids <- paste0("cache/dupCOMIDs_",hydReg, ".rds") # data frame of COMIDs with multiple HUC12 assignments
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
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
hucgage_segs <- paste0("hucgagesegs_", hydReg) # inimally-sufficient network
conf_pois <- paste0("Conf_", hydReg) # Confluence POIs
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

# 
NetworkNav <- function(inCom, nhdDF, navType){
  ##########################################
  # 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)
}


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
  while (length(upNet_DF$COMID) > 0){
    # create item for number of dangles
    count <- dim(upNet_DF)[1]
    print (dim(upNet_DF))
    # find out which level paths are downstream of dangling huc12 POIs
    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 <- nhd$COMID[nhd$Hydroseq %in% DSLP]
    # Run the upstream navigation code
    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<-nhd %>% dplyr::filter(COMID %in% inCom) %>% 
      filter(!DnHydroseq %in% Hydroseq)
    # Get the count
    count2 <- dim(upNet_DF)[1]
    # if the count has remained the same we are done and return the flowline list
    if (count == count2){
      return (inCom)
    }
  }
  # Not sure this other return is needed
  return(inCom)
}


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")
    # Look Upstream
    Upstream <- st_drop_geometry(nhd) %>%  inner_join(st_drop_geometry(divSegs) %>% filter(Divergence == 2) %>% 
                                                        select(COMID, Hydroseq), by = c("DnMinorHyd" = "Hydroseq"))
    
    # From Upstream Segment switch POI to the downstream major/main path
    inSegs_Maj <- st_drop_geometry(nhd) %>% inner_join(Upstream %>% select(COMID.y, DnHydroseq), by = c("Hydroseq" = "DnHydroseq")) %>% 
      select(COMID, COMID.y)
    
    inSegs <- inSegs %>% left_join(inSegs_Maj, by = c("COMID" = "COMID.y")) %>%
      mutate(COMID = ifelse(!is.na(COMID.y), COMID.y, COMID)) %>% select(-COMID.y)

  } else {
    print ("no divergences present in POI set")
  }
  return(inSegs)
}

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
  ##########################################
  
  # Give generic ID to Identity field
  colnames(srcData) <- c("COMID", "ID")
  
  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)
}

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)))
  
  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(POIs_fin)
}

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(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_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_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) %>%
      select (POI_ID, To_POI_ID, TotalLength.y, DAR, IncDA, geom)
  } else {
    SegSub <- segs_Down %>% filter_at(Types, all_vars(is.na(.))) %>%
      select (POI_ID, To_POI_ID, TotalLength.y, DAR, IncDA, geom)
  }
  
  #1 - POIs that need to be moved downstream
  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)
  # 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) %>%
    mutate(merged_COMID = NA)

  #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))

  # Join SegSub mod to downstream POIs
  MergedPOIs <- stationaryPOI %>%
    inner_join(st_drop_geometry(MoveDownPOI_withAtt),
               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_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_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_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)), 
                                  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))

  # Write to regional geopackage
  write_sf(nseg_raw, out_gpkg, paste0("nsegment_raw_", hydReg))
  
  # Aggregate flowlines per POI_ID to a single segment, carrying over useful information
  nsegments<-nseg_raw %>% group_by(POI_ID) %>% mutate(PathTimeMA = na_if(PathTimeMA, -9999)) %>%
    summarize(TotalLength = sum(LENGTHKM), TotalDA = max(TotDASqKM), HW = max(StartFlag),
              TT = sum(PathTimeMA)) %>%
    inner_join(st_drop_geometry(Seg2) %>% select(COMID, Hydroseq, DnHydroseq),   by=c("POI_ID"="COMID"))

  # produce a short data frame for populating TO_POI for downstream segment
  to_from<-st_drop_geometry(Seg2) %>% inner_join(st_drop_geometry(nseg_raw), by=c("DnHydroseq" = "Hydroseq")) %>%
    select(COMID.x, Hydroseq, DnHydroseq, POI_ID.x, POI_ID.y) %>% rename(POI_ID = POI_ID.x, To_POI_ID = POI_ID.y)
  
  # 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)

  # Write out dissolved segments
  write_sf(nsegments, out_gpkg, paste0("nsegment_", hydReg))
  # Write out new merged POIs
  write_sf(rbind(MergedPOIs, FinalPOI), out_gpkg, paste0("allPOIs_", hydReg))
  
  # Return POIs, segments, and raw segments
  return (list(allPOIs = rbind(MergedPOIs, FinalPOI), segs = nsegments, segsRaw = nseg_raw))
}

# Adjust confluence POIs based on if they are associated with flowline with no
#        catchment/IncDA == 0
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
  ##########################################
  
  # 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% 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 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 <- POIs %>% filter(COMID %in% repPOIs$COMID) %>% rbind(repPOIs %>% select(-old_COMID)) %>% 
    group_by(COMID) %>% summarise_each(funs(max(., na.rm=T))) 
  
  # 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)
  repPOIs_final <- repPOIs_unique %>% inner_join(poiFix %>% select(COMID, New_COMID), by = c("COMID" = "New_COMID")) %>%
    inner_join(nhd %>% select(COMID, DnHydroseq), by = c("COMID.y" = "COMID")) %>% 
    rename(new_COMID = COMID, old_COMID = COMID.y)
  
  # Get the COMID of downstream flowlines of repPOIs_final still on flowline w/o catchment 
  more_Replace <- nhdDF %>% filter(Hydroseq %in% repPOIs_final$DnHydroseq, AreaSqKM == 0) 
  
  # Subset to produce replacement data frame of POIs that still need valid downstream COMID
  moveDown <- more_Replace %>% select (COMID, DnHydroseq) %>% inner_join(repPOIs_final %>% select(new_COMID, old_COMID), 
                                                                         by = c("COMID" = "old_COMID"))
  # Column to hold the final toCOMID values
  moveDown$ToCom <- NA
  
  # If the DS replacement flowlines with Inc Area > 0 equals the number of replacement comIDs  
  while(sum(is.na(moveDown$ToCom)) > 0){
    print(sum(is.na(moveDown$ToCom)))
    
    # DS replacement COMIDs
    replacementComs <- st_drop_geometry(nhd) %>% filter(Hydroseq %in% moveDown$DnHydroseq)#, AreaSqKM > 0)
    
    # Add downstream COMIDS to the newCOMID
    moveDown <- moveDown %>% left_join(replacementComs %>% select(COMID, Hydroseq, DnHydroseq, AreaSqKM), 
                                       by = c("DnHydroseq" = "Hydroseq")) %>% 
      mutate(ToCom = ifelse(AreaSqKM > 0 & is.na(ToCom), COMID.y, ToCom), DnHydroseq = DnHydroseq.y) %>% rename(COMID = COMID.x) %>% 
      select(COMID, DnHydroseq, new_COMID, ToCom, -c(COMID.y, AreaSqKM, DnHydroseq.y))
  }
  
  # Fold this information back in to the master data frame
  repPOIs_unique <- repPOIs_unique %>% left_join(moveDown %>% select(new_COMID, ToCom), by = c("COMID" = "new_COMID")) %>%
    right_join(poiFix %>% select(New_COMID, old_COMID), by = c("COMID" = "New_COMID")) %>% 
    left_join(ToComDF, by = c("old_COMID" = "COMID.x")) %>%
    mutate(ToCom = ifelse(!is.na(COMID) & is.na(ToCom), COMID.y, ToCom)) %>% select(-c(DnHydroseq, COMID.y))
  
  return (repPOIs_unique)
}


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(nhd, poiFix, sort=T)
  seg2fix <-nhd %>% filter(COMID == poiFix)
  
  # Sorted results and filter out all flowlines w/o catchments
  upStuff <- nhd %>% filter(COMID %in% unlist(upSegs)) %>% arrange(factor(COMID, levels=upSegs)) %>%
    filter(AreaSqKM > 0)
  
  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, 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))
  
  # If 1 or other adjacent flowlines are coupled with a catchment
  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
  }
  return(newPOI)
}


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
    POInames <- lyrs$name[!is.na(lyrs$geomtype) & lyrs$geomtype== "Point"]
    print (POInames)
    for (poi in POInames){
      print (poi)
      poiType <- unlist(strsplit(poi, "_"))[1]
      subPOIs <- POIs %>%
        filter(!is.na(!!as.name(paste0("Type_", poiType))))
      
      write_sf(subPOIs, out_gpkg, poi)
    }
  } else {
    subPOIs <- POIs %>%
      filter(!is.na(!!as.name(paste0("Type_", Type))))
    write_sf(subPOIs, out_gpkg, poi)
  }
}


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") 
}