NHD_navigate.R

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

# Source Data
data_paths <- data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
huc12_xWalk <-  file.path(data_paths$nhdplus_dir, "CrosswalkTable_NHDplus_HU12.csv")
HUC12 <- file.path(data_paths$nhdplus_dir, "HUC12.rds")

# Defined during NHD_Navigate.Rmd
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
struct_POIs <- paste0("structPOIs_", hydReg) # structual POIs
struc_Net <- paste0("structNet", hydReg) # structural network
coastFlowlines <- "coast_NCA" # Coastal dendritic segments (modify this)

# Defined during POI_Collapse.Rmd
nav_gpkg <- file.path("cache", paste0("GF_",hydReg,".gpkg"))
col_gpkg <- file.path("cache", paste0("GF_",hydReg,"_Collapse.gpkg"))
gage_gpkg <- "cache/Gage_info.gpkg"
gage_cand <- "Potential_Gages"
crs <- 5070
# Output Layers
pois_merge <- paste0("merPOIs_", hydReg) # All POIs binded together
nsegment_raw_merge <- paste0("nsegment_rawMG_", hydReg)
nsegment_merge <- paste0("nsegment_MG_", hydReg)
hydReg_gage_Table <- file.path("cache", paste0("R", hydReg, "_Gages.csv"))
CONUS_gage_Table <- "data/gages_MDA.rds"

# Defined during NonDend.Rmd
nd_gpkg <- file.path("cache", paste0("ND_",hydReg,".gpkg"))
coast_FL <- paste0("coastal_FL_", hydReg)
nhd_aggLine <- paste0("nhd_aggline_", hydReg)
agg_fline <- paste0("agg_fline_", hydReg)
reg_cats_lyr <- paste0("nhd_cats", hydReg)
agg_cats_lyr <- paste0("agg_cats", hydReg)
outlets_lyr <- paste0("outlets", hydReg)
lookup_table <- file.path("cache/", paste0(hydReg, "_lookup.csv"))
cat0_lyr <- paste0("cat0_", hydReg)
hru0_lyr <- paste0("hru0_", hydReg)
nd_table_out <- file.path("cache/", paste0(hydReg, "_ND_lookup.csv"))
nhdCoast_net <- paste0("nhdCoast_network_", hydReg)
catCoast_net <- paste0("catCoast_", hydReg)
cat0_lyr <- "catchment_coast" # modify this
hru0_lyr <- "hru_coast" # modify this
sink_lyr_out <- paste0("sink_", hydReg)
missing_cat_lyr <- paste0("missingCats_", hydReg)
sink_cats_out <- paste0("Terminal_Cats_", hydReg)
HUC12_lyr_out <- paste0("HUC12_", hydReg)
cat_final_out <- paste0("Final_Cats_", hydReg)
sink_net_out <- paste0("sink_Net_", hydReg)
rem_cat_out <- paste0("Rem_Cats_", hydReg)
hru_NCA_out <- paste0("HRU_NCA_", hydReg)


NetworkNav <- function(inCom, nhdDF, withTrib){
  ##########################################
  # 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
  #   withTrib : (logical) flag for if the upstream navigation should include tributaries
  #              or stick to mainstem level path
  #
  # Returns:
  #   mydata : (list) list of COMIDs to connect dangle to network
  ##########################################
  if (missing(withTrib)){
    seg <- nhdplusTools::get_UM(nhdDF, inCom, include = TRUE)
  } else {
    seg <- nhdplusTools::get_UT(nhdDF, inCom)
  }
  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)))
    # Append result to existing segment list
    incom <- append(incom, upNet)
    # Get the same variable as above
    upnet_DF <- filter(nhd, COMID %in% incom, !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
  ##########################################
  div_segs <- filter(nhd, COMID %in% insegs$COMID)
  if (2 %in% div_segs$Divergence){
    print ("Switching divergence to other fork")
    # Look Upstream
    upstream <- st_drop_geometry(nhd) %>% 
      inner_join(st_drop_geometry(div_segs) %>% 
                   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")
  
  sub_segs <- filter(nhd, COMID %in% srcData$COMID) 
  #merge(st_drop_geometry(sub_segs %>% select(COMID))
  POIs <- sub_segs %>% 
    get_node(., position = "end") %>% 
    mutate(COMID = sub_segs$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, DA_diff){
  ##########################################
  # Moves down/up POIs based on data theme
  # 
  # Arguments:
  #   out_gpkg : (geopackage) Output geopackage
  #   POIs: (data frame) Original  POIs
  #   Seg : (data frame) dissolved segments
  #   Seg2 : (data frame) raw segments
  #   exp : (string) Type of thematic POI being moved round
  #   DA_diff : (float) threshold of drainage area difference to
  #             consider when merging two POIs
  #
  # Returns:
  #   list : 1 - New set of POIs
  #          2/3 - correpsonding segments (both raw and dissolved)
  ##########################################
  
  # Segs/POIs to collapse downstream (fold POIs together downstream if valid catchment downstream)
  # Don't want to move down to TE facility location
  segs_down <- Seg %>% 
    inner_join(select(st_drop_geometry(.), c(POI_ID, TotalDA, TotalLength)), 
                                    by = c("To_POI_ID" = "POI_ID")) %>%
    inner_join(st_drop_geometry(POIs), by = c("POI_ID" = "COMID")) %>%
    # Keep WBOut and In to preserve NHD Waterbody geometry, Keep Confluences and TE where they are
    filter(is.na(Type_WBOut), is.na(Type_WBIn), is.na(Type_Conf), is.na(Type_TE)) %>% 
    # Select POIs within the correct drainage area ratio and fit the size criteria
    mutate(DAR = TotalDA.y/TotalDA.x, IncDA = TotalDA.y - TotalDA.x) %>% 
    st_drop_geometry() %>%
    # If the drainage area ratio is within acceptable range (1.05)
    filter(DAR < (1 + DA_diff)) %>% 
    select(POI_ID, HW, To_POI_ID, Type_HUC12, Type_Gages, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf, DAR, IncDA)
  
  # Segs/POIS to collapose upstream if HUC12 or Gage
  up_POIs <- POIs %>% filter(is.na(Type_Conf))
  segs_up <- Seg %>% inner_join(select(filter(st_drop_geometry(.), POI_ID %in% up_POIs$COMID), 
                                       c(POI_ID, To_POI_ID, TotalDA, TotalLength)), 
                                by = c("POI_ID" = "To_POI_ID")) %>%
    inner_join(st_drop_geometry(POIs), by = c("POI_ID" = "COMID")) %>% 
    rename(From_POI_ID = POI_ID.y) %>%
    # Don't want to move outlets upstream, don't move confluences upstream
    filter(is.na(Type_WBOut), is.na(Type_WBIn), is.na(Type_Conf), is.na(Type_TE)) %>% 
    # Select POIs within the correct drainage area ratio and fit the size criteria
    mutate(DAR = TotalDA.y/TotalDA.x, IncDA = TotalDA.x - TotalDA.y) %>% st_drop_geometry() %>% 
    # If the drainage area ratio is within acceptable range (0.95)
    filter(DAR >= (1 - DA_diff)) %>% 
    select(POI_ID, HW, From_POI_ID, Type_HUC12, Type_Gages, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf, DAR, IncDA)
  
  # compiled list of POIs to move up or down
  seg_choices <- POIs %>% 
    left_join(segs_down %>% 
                select(POI_ID, To_POI_ID, DAR, IncDA), by = c("COMID" = "POI_ID")) %>%
    left_join(segs_up %>% 
                select(POI_ID, From_POI_ID, DAR, IncDA), by = c("COMID" = "POI_ID")) #%>%
    #filter(IncDA.x < 10 | IncDA.y < 10)

  # Filter by POI Type
  Types <- c("Type_HUC12", "Type_Gages", "Type_TE", "Type_Conf", "Type_NID", "Type_WBIn", "Type_WBOut", "Type_TE")
  Types <- Types[Types != exp]
  static_POIs <- POIs %>% 
    filter(!is.na(Type_Gages) & (Type_Conf == 1 | !is.na(Type_WBOut) | !is.na(Type_WBIn)))
  
  seg_sub <- seg_choices %>% 
    filter_at(Types, all_vars(is.na(.))) %>%
    filter(!COMID %in% static_POIs$COMID) %>%
    select (COMID, To_POI_ID, DAR.x, IncDA.x, From_POI_ID, DAR.y, DAR.y, IncDA.y, geom) %>%
    st_as_sf()
  
  st_write(seg_sub, out_gpkg, "Seg_choices", delete_layer = F, append = T)
  
  segsub_dir <- seg_sub %>% 
    mutate(DirDA = ifelse(is.na(DAR.y), "Down", 
                          ifelse(is.na(DAR.x), "Up",
                                 ifelse(IncDA.x < IncDA.y, "Down", "Up"))),
                                  Move = ifelse(DirDA == "Down", To_POI_ID, From_POI_ID))
  
  Mult <- select(segsub_dir, COMID, Move) %>% 
    inner_join(st_drop_geometry(segsub_dir) %>% 
                                         select(COMID, Move), by = c("Move" = "COMID"), suffix = c(".A", ".B")) 
  
  segsub_dir <- filter(segsub_dir, !COMID %in% Mult$Move)
  
  #1 - POIs that need to be moved downstream
  move_POI <- filter(POIs, COMID %in% segsub_dir$COMID)
  # POIs that are stationary and their infomration will be merged with upstream  POIs
  stationary_POI <- filter(POIs, !COMID %in% move_POI$COMID)
  # Final Set to be merged with that don't involve either sets above
  #FinalPOI <- POIs %>% filter(!COMID %in% SegSub_Dir$COMID & !COMID %in% SegSub_Dir$Move) %>%
  #  mutate(merged_COMID = NA)

  #2 - Join SegSub assignment to POIs to bring over POI attributes
  movedownPOI_withatt <- select(st_drop_geometry(segsub_dir), COMID, DirDA, Move) %>%
    left_join(st_drop_geometry(move_POI), by = "COMID")  

  # Join SegSub mod to downstream POIs
  merged_POIs <- stationary_POI %>%
    inner_join(movedownPOI_withatt,
               by = c("COMID" = "Move"), suffix = c(".x", ".y")) %>%
    # Bring over relevant attributes
    mutate(Type_HUC12 = ifelse(is.na(Type_HUC12.y), Type_HUC12.x, Type_HUC12.y)) %>%
    mutate(Type_Gages_A = ifelse(is.na(Type_Gages.x) & !is.na(Type_Gages.y), Type_Gages.y, Type_Gages.x)) %>%
    # 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) & !is.na(Type_Gages.x), Type_Gages.y, NA)) %>%
    mutate(Type_WBIn = ifelse(is.na(Type_WBIn.y), Type_WBIn.x, Type_WBIn.y)) %>% 
    mutate(Type_WBOut = ifelse(is.na(Type_WBOut.y), Type_WBOut.x, Type_WBOut.y)) %>%
    mutate(Type_TE = ifelse(is.na(Type_TE.y), Type_TE.x, Type_TE.y)) %>%
    mutate(Type_NID = ifelse(is.na(Type_NID.y), Type_NID.x, Type_NID.y)) %>%
    mutate(Type_Conf = ifelse(is.na(Type_Conf.y), Type_Conf.x, Type_Conf.y)) %>%
    mutate(old_COMID = COMID.y) %>% 
    st_as_sf()

  # Bring new POI data over to new data
  fin_POIs <- filter(POIs, !COMID %in% merged_POIs$old_COMID) %>% 
    left_join(st_drop_geometry(merged_POIs) %>% select(COMID, Type_Gages_A), by = "COMID") %>%
    mutate(Type_Gages = ifelse(!is.na(Type_Gages_A), Type_Gages_A, Type_Gages)) %>% select(-Type_Gages_A)
  
  # Write to regional geopackage (x = new POI, y = old)
  write_sf(segsub_dir, out_gpkg, paste0("POI_Change_", hydReg), append = T)
  write_sf(merged_POIs, 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(merged_POIs),c(COMID, old_COMID)), 
                                  by = c("POI_ID" = "old_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(fin_POIs, out_gpkg, paste0("allPOIs_", hydReg))
  
  # Return POIs, segments, and raw segments
  return (list(allPOIs = fin_POIs, segs = nsegments, segsRaw = nseg_raw))
}

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 <- select(nhdDF, COMID, DnHydroseq) %>% 
    inner_join(select(nhdDF, COMID, Hydroseq), by = c("DnHydroseq" = "Hydroseq"))
  
  # Find segments with POIs where there is no corresponding catchment
  unCon_POIs <- filter(nhdDF, COMID %in% POIs$COMID, AreaSqKM == 0) 
  # apply initial move function to the whole POI Set
  poi_fix <- filter(POIs, COMID %in% unCon_POIs$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
  rep_POIs <- filter(nhd, COMID %in% poi_fix$New_COMID) %>% 
    get_node(., position = "end") %>% 
    mutate(COMID = (filter(nhd, COMID %in% poi_fix$New_COMID) %>% 
                      pull(COMID))) %>% 
    left_join(poi_fix %>% 
                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
  rep_POIs_unique <- filter(POIs, COMID %in% rep_POIs$COMID) %>% 
    rbind(rep_POIs %>% 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
  rep_POIs_unique[mapply(is.infinite, rep_POIs_unique)] <- NA
  
  # Get final list of POIs being moved, and the downstream POI they are linking to (many : 1 join)
  rep_POIs_final <- rep_POIs_unique %>% 
    inner_join(poi_fix %>% 
                 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 <- filter(nhdDF, Hydroseq %in% rep_POIs_final$DnHydroseq, AreaSqKM == 0) 
  
  # Subset to produce replacement data frame of POIs that still need valid downstream COMID
  move_down <- select(more_replace, COMID, DnHydroseq) %>% 
    inner_join(rep_POIs_final %>% select(new_COMID, old_COMID), by = c("COMID" = "old_COMID"))
  
  # Column to hold the final toCOMID values
  move_down$ToCom <- NA
  
  # If the DS replacement flowlines with Inc Area > 0 equals the number of replacement comIDs  
  while(sum(is.na(move_down$ToCom)) > 0){
    print(sum(is.na(move_down$ToCom)))
    
    # DS replacement COMIDs
    replacement_coms <- filter(st_drop_geometry(nhd), Hydroseq %in% move_down$DnHydroseq)#, AreaSqKM > 0)
    
    # Add downstream COMIDS to the newCOMID
    move_down <- move_down %>% 
      left_join(replacement_coms %>% 
                  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
  rep_POIs_unique <- rep_POIs_unique %>% 
    left_join(move_down %>% select(new_COMID, ToCom), by = c("COMID" = "new_COMID")) %>%
    right_join(poi_fix %>% 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 (rep_POIs_unique)
}


movePOI_NA_DA <- function(poi_fix, 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:
  #   poi_fix : (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
  up_segs <- nhdplusTools::get_UM(nhd, poi_fix, sort=T)
  seg2fix <- filter(nhd, COMID == poi_fix)
  
  # Sorted results and filter out all flowlines w/o catchments
  upstuff <- filter(nhd, COMID %in% unlist(up_segs)) %>% 
    arrange(factor(COMID, levels = up_segs)) %>%
    filter(AreaSqKM > 0)
  
  down_segs <- nhdplusTools::get_DM(nhd, poi_fix, sort=T)
  downstuff <- filter(nhd, COMID %in% unlist(down_segs)) %>% 
    arrange(factor(COMID, levels = down_segs)) %>%
    filter(AreaSqKM > 0)
  
  # combine into one dataframe, select up/downstream seg with least change in total drainage area
  near_FL <- rbind(select(upstuff, 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(near_FL$AreaSqKM) > 0){
    new_POI <- near_FL$COMID[(which.min(abs(seg2fix$TotDASqKM - near_FL$TotDASqKM)))]
  } else {
    # Remove POI if not catchment associated with flowlines upstream or downstream
    print (poi_fix)
    print ("US and DS flowlines also have no catchment, removing POI")
    new_POI <- NA
  }
  return(new_POI)
}


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


structPOIsNet <- function(ncombined, nhd, final_POIs, 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
  #
  # Returns:
  #   writes Structural POIs and segments to geopackage 
  ##########################################
  # Terminal outlets from the initial network
  termout <- filter(nhd, !Hydroseq %in% DnHydroseq & TerminalFl == 1 & !COMID %in% final_POIs$COMID)
  
  # POIs that are also terminal outlets
  out_POIs <- filter(nhd, COMID %in% final_POIs$COMID & TerminalFl == 1)
  
  # Confluence POIs
  conf_POIs <- filter(ncombined, COMID %in% final_POIs$COMID[final_POIs$Type_Conf == 1])
  
  # Waterbody outlet POIs
  wb_POIs <- filter(ncombined, COMID %in% final_POIs$COMID[!is.na(final_POIs$Type_WBOut) | !is.na(final_POIs$Type_WBIn)])
  
  # Waterbody inlets POIs
  struc_flines <- termout %>% 
    bind_rows(out_POIs, conf_POIs, wb_POIs) %>% 
    arrange(COMID)
  struc_POIs <- get_node(struc_flines, position = "end") %>% 
    mutate(COMID = struc_flines$COMID, LevelPathI = struc_flines$LevelPathI) %>% 
    st_drop_geometry()
  
  # Add back in any levelpaths missing (shouldn't be much)
  miss_LP <- filter(ncombined, COMID %in% final_POIs$COMID) %>% 
    filter(!LevelPathI %in% struc_POIs$LevelPathI)
  
  # Only bind if there are rows present
  if (nrow(miss_LP) > 0){
    lp_POIs <- c(miss_LP$LevelPathI, struc_POIs$LevelPathI)
  } else {
    lp_POIs <- struc_POIs$LevelPathI
  }
  
  # final stuctural POIs, order by COMID to match with strucPOIs
  final_struct <- filter(nhd, Hydroseq %in% lp_POIs) %>% 
    arrange(COMID)
  struct_POIs <- get_node(final_struct, position = "end") %>% 
    mutate(COMID = final_struct$COMID) 
  
  #  produce unique set of flowlines linking POIs
  final_net <- unique(unlist(lapply(unique(final_struct$COMID), NetworkNav, st_drop_geometry(nhd), "up")))
  
  # Subset NHDPlusV2 flowlines to navigation results and write to shapefile
  structnet <- filter(nhd, COMID %in% final_net & grepl(paste0("^", hydReg, ".*"), .data$VPUID)) 
  
  # Write out minimally-sufficient network and POIs as a list
  return(list(struc_POIs = struc_POIs, structnet = structnet))
}