NHD_navigate.Rmd

---
title: "NHD Navigate"
output: html_document
editor_options:
  chunk_output_type: console
---
This notebook Generate Segments and POIS from POI data sources and builds
a minimally-sufficient stream network connecting all POis, as well as generating first-cut of
national segment network.

```{r}
knitr::opts_chunk$set(error = TRUE)
```

```{r setup}
#  Load data and libraries ----------------------

ptm<-proc.time()

library(nhdplusTools)
library(dplyr)
library(sf)
library(hyRefactor)
library(tidyr)

hydReg <- "01"

# Load custom functions
source("R/utils.R")
source("R/NHD_navigate.R")

# Load user-defined data paths
data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
# Load NHDPlus Data if precprocessed to RDS format
nhdplus_path(data_paths$nhdplus_gdb)
staged_nhd <- stage_national_data(nhdplus_data = data_paths$nhdplus_gdb, 
                                  output_path  = data_paths$nhdplus_dir)

# Output Geopackage where we are writing output layers
out_gpkg <- file.path("cache", paste0("GF_",hydReg,".gpkg"))
```

```{r huc12_pois}
#  Derive or load HUC12 POIs ----------------------
if(needs_layer(out_gpkg, huc12_pois) | needs_layer(out_gpkg, nhdflowline)) {

  # Subset NHD by VPU
  nhd <- VPU_Subset(staged_nhd$flowline, hydReg)
  write_sf(nhd, out_gpkg, nhdflowline)
  
  # Read in HUC12 outlets, get rid of duplicate COMID/HUC12 rows
  hydReg2 <- ifelse(hydReg %in% c("10U", "10L"), substr(hydReg, start = 1, stop = 2), hydReg)
  HUC12 <- read_sf(data_paths$hu12_points_path, "hu_outlets") %>% 
    filter(grepl(paste0("^", hydReg2,".*"), .data$HUC12)) %>% st_drop_geometry() %>%
    select(COMID, HUC12) %>% distinct() %>% switchDiv(., nhd)

  # Identify multiple COMIDs for HUC12
  dupCOMIDs <- HUC12 %>% group_by(COMID) %>% select(COMID, HUC12) %>% filter(n()>1)
  saveRDS(dupCOMIDs, dup_comids)
  
  # Craete physical POIs from segments
  huc12POIs <- nhd %>% filter(COMID %in% HUC12$COMID) %>% POI_creation(., nhd) %>%  
    #merge(st_drop_geometry(HUC12), by = "COMID") %>% #Need HUC12 mods before uncommenting
    mutate(Type_HUC12 = 1, Type_WBIn = NA, Type_WBOut = NA, Type_GagesIII = NA, Type_TE = NA, Type_NID = NA, Type_Conf = NA) %>% 
    select(COMID, Type_HUC12, Type_GagesIII, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf, geometry)
  
  # Write out geopackage layer representing POIs for given theme
  write_sf(huc12POIs, out_gpkg, huc12_pois)
  tmpPOIs <- st_drop_geometry(huc12POIs)
} else {
  # Load HUC12 POIs as the tmpPOIs if they already exist
  tmpPOIs <- st_drop_geometry(read_sf(out_gpkg, huc12_pois))   
  nhd <- read_sf(out_gpkg, nhdflowline)
}
```

```{r gagesIII_pois}
#  Derive or load GAGESIII POIs ----------------------
# Note several POIs have multiple gages indexed
if(needs_layer(out_gpkg, gagesIII_pois)) {
  
  # Read in GAGESIII shapefile
  gagesIII <- read_sf(data_paths$gagesiii_points_path, "GagesIII_gages") %>% 
    filter(grepl(paste0("^", hydReg, ".*"), .data$NHDPlusReg)) %>% 
    mutate(COMID = as.integer(COMID)) %>% inner_join(st_drop_geometry(nhd), on = "COMID") %>%
    switchDiv(., nhd) %>% filter(TotDASqKM > 1)
  
  # Derive GAGESIII POIs
  gagesIIIPOIs <- nhd %>% filter(COMID %in% gagesIII$COMID) %>% POI_creation(., nhd) %>% 
    left_join(st_drop_geometry(gagesIII), by = "COMID") %>% 
    mutate(Type_HUC12 = NA, Type_GagesIII = Gage_no, Type_TE = NA, Type_NID = NA, Type_WBIn = NA, Type_WBOut = NA, Type_Conf = NA) %>%
    select(COMID, Type_HUC12, Type_GagesIII, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf, geometry) 
  
  # Assign HUC12 POIs streamgage ID if they are co-located
  tmpPOIs <- tmpPOIs %>% left_join(filter(st_drop_geometry(gagesIIIPOIs) %>% select(COMID, Type_GagesIII), COMID %in% tmpPOIs$COMID), by = "COMID", all.x = TRUE) %>%
    mutate(Type_GagesIII = ifelse(!is.na(Type_GagesIII.y), Type_GagesIII.y, NA)) %>% 
    select(COMID, Type_HUC12, Type_GagesIII, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf) 
  
  # Write out geopackage layer representing POIs for given theme
  write_sf(gagesIIIPOIs, out_gpkg, gagesIII_pois)
  tmpPOIs <- do.call("rbind", list(tmpPOIs, st_drop_geometry(gagesIIIPOIs) %>% filter(!COMID %in% tmpPOIs$COMID)))
} else {
  # Load GAGESIII POIs if they already exist and append to tmpPOIs
  gagesIIIPOIs <- read_sf(out_gpkg, gagesIII_pois)
  tmpPOIs <- do.call("rbind", list(tmpPOIs, st_drop_geometry(gagesIIIPOIs) %>% filter(!COMID %in% tmpPOIs$COMID)))
}
```

```{r TE_pois}
#  Derive or load Thermoelectric POIs ----------------------
if(needs_layer(out_gpkg, TE_pois)) {
  
  # Read in Thermoelectric shapefile
  TE_fac <- read_sf(data_paths$TE_points_path, "2015_TE_Model_Estimates_lat.long_COMIDs") %>% 
    mutate(COMID=as.integer(COMID)) %>% filter(COMID > 0) %>%
    inner_join(., select(st_drop_geometry(nhd), c(COMID, VPUID)), by = "COMID") %>%
    filter(grepl(paste0("^", hydReg, ".*"), .data$VPUID)) 
  
  # Get unique list of COMIDs
  TE_fac_unique <- TE_fac %>% group_by(COMID) %>% summarize(EIA_PLANT = paste0(unique(EIA_PLANT_), collapse = " ")) %>%
    switchDiv(., nhd) %>% st_drop_geometry
  
  # Derive TE POIs
  TE_POIs <- nhd %>% filter(COMID %in% TE_fac_unique$COMID) %>% POI_creation(., nhd) %>% 
    inner_join(TE_fac_unique %>% select(COMID, EIA_PLANT), by = "COMID") %>% 
    mutate(Type_HUC12 = NA, Type_GagesIII = NA, Type_TE = EIA_PLANT, Type_NID = NA, Type_WBIn = NA, Type_WBOut = NA, Type_Conf = NA) %>%
    select(-EIA_PLANT)
    
  # Assign existing POIs as a Type_TE if they are co-located
  tmpPOIs <- tmpPOIs %>% left_join(filter(st_drop_geometry(TE_POIs) %>% select(COMID, Type_TE)), by = "COMID", all.x = TRUE) %>%
    mutate(Type_TE = ifelse(!is.na(Type_TE.y), Type_TE.y, NA)) %>%
    select(COMID, Type_HUC12, Type_GagesIII, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf)
  
  # Write out geopackage layer representing POIs for given theme
  write_sf(TE_POIs, out_gpkg, TE_pois)
  tmpPOIs <- do.call("rbind", list(tmpPOIs, st_drop_geometry(TE_POIs) %>% filter(!COMID %in% tmpPOIs$COMID)))
} else {
  # Load TE POIs if they already exist
  TE_POIs <- read_sf(out_gpkg, TE_pois) 
  tmpPOIs <- do.call("rbind", list(tmpPOIs, st_drop_geometry(TE_POIs) %>% filter(!COMID %in% tmpPOIs$COMID))) 
}
```

```{r NID_pois}
#  Derive or load Thermoelectric POIs ----------------------
if(needs_layer(out_gpkg, NID_pois)) {
  
  # Read in Thermoelectric shapefile
  NID_fac <- read.csv(file.path(data_paths$NID_points_path, "/NID_attributes_20170612.txt")) %>% 
    filter(ONNHDPLUS == 1, grepl(paste0("^", hydReg, ".*"), .data$VPU)) %>% rename(COMID = FlowLcomid) %>% 
    switchDiv(., nhd)
  
  # Get unique list of COMIDs and NID facilities
  NID_fac_unique <- NID_fac %>% group_by(COMID) %>% summarize(NID_ID = paste0(unique(NIDID), collapse = " ")) 
  
  # Derive TE POIs
  NID_POIs <- nhd %>% filter(COMID %in% NID_fac$COMID) %>% POI_creation(., nhd) %>% 
    inner_join(NID_fac_unique %>% select(COMID, NID_ID), by = "COMID") %>% 
    mutate(Type_HUC12 = NA, Type_GagesIII = NA, Type_Conf = NA, Type_TE = NA, Type_WBIn = NA, Type_WBOut = NA, Type_NID = NID_ID) %>%
    select(-NID_ID)
  
  # Assign existing POIs as a Type_TE if they are co-located
  tmpPOIs <- tmpPOIs %>% left_join(filter(st_drop_geometry(NID_POIs) %>% select(COMID, Type_NID)), by = "COMID", all.x=TRUE) %>%
    mutate(Type_NID = ifelse(!is.na(Type_NID.y), Type_NID.y, NA)) %>% 
    select(COMID, Type_HUC12, Type_GagesIII, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf) 
  
  # Write out geopackage layer representing POIs for given theme
  write_sf(NID_POIs, out_gpkg, NID_pois)
  tmpPOIs <- do.call("rbind", list(tmpPOIs, st_drop_geometry(NID_POIs) %>% filter(!COMID %in% tmpPOIs$COMID)))
} else {
  # Load TE POIs if they already exist
  NID_POIs <- read_sf(out_gpkg, NID_pois)   
  tmpPOIs <- do.call("rbind", list(tmpPOIs, st_drop_geometry(NID_POIs) %>% filter(!COMID %in% tmpPOIs$COMID)))  
}
```

```{r hucgages_segs}
# Build minimally-sufficient network to connect POIs ----------------------
if(needs_layer(out_gpkg, hucgage_segs)) {
  
  # Navigate upstream from each POI
  upNet <- unique(unlist(lapply(unique(tmpPOIs$COMID), NetworkNav, st_drop_geometry(nhd), "up")))
  
  # Connect POIs that don't have connecting levelpath looking downstream and 
  #         produce unique set of flowlines linking POIs
  finalNet <- unique(NetworkConnection(upNet, st_drop_geometry(nhd)))
   
  # Subset NHDPlusV2 flowlines to navigation results and write to shapefile
  hucgagessegs <- nhd %>% filter(COMID %in% finalNet & grepl(paste0("^", hydReg, ".*"), .data$VPUID))
  
  # Write out minimally-sufficient network
  write_sf(hucgagessegs, out_gpkg, hucgage_segs) 
} else {
  # Load minimally-sufficient network
  hucgagessegs <- read_sf(out_gpkg, hucgage_segs)
}
```

```{r waterbody_pois}
#  Derive or load GAGESIII POIs ----------------------
if(needs_layer(out_gpkg, WBin_POIs)) {
  
  # Read in waterobdies
  wb <- readRDS("data/NHDPlusNationalData/nhdplus_waterbodies.rds") %>% 
    mutate(FTYPE = as.character(FTYPE)) %>% filter(FTYPE %in% c("LakePond", "Reservoir")) %>%
    filter(AREASQKM > 1)
  # Get Waterbodies that intersect with Region flowlines
  WBs_VPU <- wb %>% filter(COMID %in% hucgagessegs$WBAREACOMI)
  # Write out waterbodies
  write_sf(WBs_VPU, out_gpkg, WBs)
  
  # Segments that are in waterbodies
  wbSegs <- nhd %>% filter(WBAREACOMI > 0 & WBAREACOMI %in% WBs_VPU$COMID)
  wbOut <- hucgagessegs %>% filter(WBAREACOMI %in% WBs_VPU$COMID) %>% group_by(WBAREACOMI) %>% slice(which.min(Hydroseq)) %>%
    switchDiv(., nhd) 
  
  wbIn <- hucgagessegs %>% filter(!WBAREACOMI %in% WBs_VPU$COMID) %>% filter(DnHydroseq %in% wbSegs$Hydroseq) %>%
    switchDiv(., nhd) %>% inner_join(st_drop_geometry(hucgagessegs) %>% select(Hydroseq, WBAREACOMI), by = c("DnHydroseq" = "Hydroseq"))
  
  # Craete physical POIs from segments
  WBoutPOIs <- POI_creation(wbOut, nhd) %>% 
    inner_join(st_drop_geometry(wbOut) %>% select(COMID, WBAREACOMI), by = "COMID") %>%
    mutate(Type_HUC12 = NA, Type_GagesIII = NA, Type_TE = NA, Type_NID = NA, Type_WBIn = NA, Type_WBOut = WBAREACOMI, Type_Conf = NA) %>% 
    select(COMID, Type_HUC12, Type_GagesIII, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf, geometry) %>% 
    filter(!COMID %in% wbIn$COMID) 
  
  WBinPOIs <- POI_creation(wbIn, nhd) %>% 
    inner_join(st_drop_geometry(wbIn) %>% select(COMID, WBAREACOMI.y), by = "COMID") %>%
    mutate(Type_HUC12 = NA, Type_GagesIII = NA, Type_TE = NA, Type_NID = NA, Type_WBIn = WBAREACOMI.y, Type_WBOut = NA, Type_Conf = NA) %>% 
    select(COMID, Type_HUC12, Type_GagesIII, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf, geometry, -WBAREACOMI.y) %>%
    filter(!COMID %in% wbOut$COMID)
  
    # Assign HUC12 POIs WB in/out if they are co-located
  tmpPOIs <- tmpPOIs %>% left_join(filter(st_drop_geometry(WBinPOIs) %>% select(COMID, Type_WBIn)), by = "COMID", all.x = TRUE) %>%
    mutate(Type_WBIn = ifelse(!is.na(Type_WBIn.y), Type_WBIn.y, NA)) %>%
    select(COMID, Type_HUC12, Type_GagesIII, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf) %>%
    left_join(filter(st_drop_geometry(WBoutPOIs) %>% select(COMID, Type_WBOut)), by = "COMID", all.x = TRUE) %>%
    mutate(Type_WBOut = ifelse(!is.na(Type_WBOut.y), Type_WBOut.y, NA)) %>%
    select(COMID, Type_HUC12, Type_GagesIII, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf)
  
  write_sf(WBoutPOIs, out_gpkg, WBout_POIs)
  write_sf(WBinPOIs, out_gpkg, WBin_POIs)
  
  tmpPOIs <- do.call("rbind",list(tmpPOIs, st_drop_geometry(WBinPOIs) %>% filter(!COMID %in% tmpPOIs$COMID), 
                                  st_drop_geometry(WBoutPOIs) %>% filter(!COMID %in% tmpPOIs$COMID)))
  
  # Some of the WBAReaCOMIDs grab small parts of headwater lakes that don't have segments upstream
  missingSegs <- tmpPOIs %>% filter(!COMID %in% hucgagessegs$COMID)
  hucgagessegs <- hucgagessegs %>% rbind(., nhd %>% filter(COMID %in% NetworkConnection(missingSegs$COMID, st_drop_geometry(nhd)))) 
  
} else {
  # Load HUC12 POIs if they already exist
  WBinPOIs <- read_sf(out_gpkg, WBin_POIs)   
  WBoutPOIs <- read_sf(out_gpkg, WBout_POIs) 
  tmpPOIs <- do.call("rbind", list(tmpPOIs, st_drop_geometry(WBinPOIs) %>% filter(!COMID %in% tmpPOIs$COMID),
                                   st_drop_geometry(WBoutPOIs) %>% filter(!COMID %in% tmpPOIs$COMID)))
}
```

```{r}
# Derive POIs at confluences where they are absent ----------------------
if(needs_layer(out_gpkg, conf_pois)) {

  # Create new confluence POIs
  confPOIs <- hucgagessegs %>% filter(DnHydroseq > 0) %>% group_by(DnHydroseq) %>% filter(n()>1) %>%
    switchDiv(., nhd)  
  
  confPOIs <- POI_creation(confPOIs, nhd) %>% 
    left_join(st_drop_geometry(.), by = "COMID") %>% 
    mutate(Type_HUC12 = NA, Type_GagesIII = NA, Type_TE = NA, Type_NID = NA, Type_WBIn = NA, Type_WBOut = NA, Type_Conf = 1) %>%
    select(COMID, Type_HUC12, Type_GagesIII, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf, geometry) 
  
  # Mark existing POIs if they are a confluence
  tmpPOIs <- tmpPOIs %>% 
    left_join(filter(st_drop_geometry(confPOIs) %>% select(COMID, Type_Conf)), by = "COMID", all.x = TRUE) %>%
    mutate(Type_Conf=ifelse(!is.na(Type_Conf.y), 1, NA)) %>%
    select(COMID, Type_HUC12, Type_GagesIII, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf)
  
  # Write out shapefile representing POIs for given theme
  write_sf(confPOIs, out_gpkg, conf_pois)
  tmpPOIs<-rbind(tmpPOIs, st_drop_geometry(confPOIs) %>% filter(!COMID %in% tmpPOIs$COMID))
} else {
  confPOIs <- read_sf(out_gpkg, conf_pois)
  tmpPOIs <- do.call("rbind", list(tmpPOIs, st_drop_geometry(confPOIs) %>% filter(!COMID %in% tmpPOIs$COMID)))
}
```

```{r Final POIs}
# Derive final POI set ----------------------
if(needs_layer(out_gpkg, pois_all)) {
  
  # Find segments with POIs where there is no corresponding catchment
  unConPOIs <- st_drop_geometry(nhd) %>% filter(COMID %in% tmpPOIs$COMID, AreaSqKM == 0) 
  
  # For confluence POIs falling on Flowlines w/o catchments, derive upstream valid flowline, 
  #     and retrieve original downstream COMID if necessary
  poiFix <- DS_poiFix(tmpPOIs, nhd) 
  
  poiFix_final <- poiFix %>% select(-c(ToCom, old_COMID)) %>% filter(!is.na(COMID)) %>% distinct()
  
  # Replace the POI_ID values, first with existing POIs
  finalPOIs_att <- tmpPOIs %>% filter(!COMID %in% c(poiFix$old_COMID, poiFix$COMID)) %>% rbind(poiFix_final)
  finalPOI_geom <- nhd %>% filter(COMID %in% finalPOIs_att$COMID) %>% POI_creation(., nhd)
  finalPOIs <- finalPOIs_att %>% left_join(finalPOI_geom, by = "COMID") %>% st_as_sf() 

  poiFix_sf <- nhd %>% filter(COMID %in% poiFix$old_COMID) %>% select(COMID) %>% POI_creation(., nhd) %>%
    inner_join(poiFix, by = c("COMID" = "old_COMID")) %>% rename(new_COMID = COMID.y)

  write_sf(finalPOIs, out_gpkg, pois_all) 
  writePOIs(finalPOIs, out_gpkg)
  write_sf(poiFix_sf, out_gpkg, "poi_Fix")
} else {
  finalPOIs <- read_sf(out_gpkg, pois_all)
  poiFix_sf <- read_sf(out_gpkg, "poi_Fix")
}
```

```{r Final segments}
# Derive first cut of segments ----------------------
if(needs_layer(out_gpkg, nsegment_raw)) {
  
  # Sort POIs by Levelpath and Hydrosequence in upstream to downstream order
  segPOIs <- st_drop_geometry(nhd) %>% filter(COMID %in% finalPOIs$COMID, COMID %in% hucgagessegs$COMID) %>%
    arrange(desc(LevelPathI), desc(Hydroseq))
  hucgagessegs_ns <- st_drop_geometry(hucgagessegs)
  
  # Create column to hold  POI_ID values (most downstream COMID of each segment)
  hucgagessegs_ns$POI_ID<-0
  # 1733450 1733448 1732360 1734154 1734156
  # Begin with upstream most Levelpath/Hydrosequenc
  for (LP in unique(segPOIs$LevelPathI)){
    # For each level path
    segPOIs_LP <- segPOIs %>% filter(LevelPathI == LP)
    # Go through sorted list of US -> DS comids
    for (POI in segPOIs_LP$COMID){

      # Re-iterate only on flowlines that have no POI_ID to increase speed
      hucgagessegs_ns<-hucgagessegs_ns %>% filter(POI_ID == 0)
      
      subPOI <- hucgagessegs_ns %>% filter (COMID == POI)
      
      # Get upstream segments
      Seg<-nhdplusTools::get_UM(hucgagessegs_ns, POI, include = T)
      
      # Assign POI_ID
      hucgagessegs_ns <- hucgagessegs_ns %>%  
        mutate(POI_ID = ifelse(COMID %in% Seg & POI_ID == 0, POI, POI_ID))
      
      # Derive incremental segment composed of flowlines
      incSeg<-hucgagessegs_ns %>% filter(POI_ID == POI)
  
      #print (dim(incSeg))
      if(!exists("IncSegs")){
        IncSegs<-incSeg
      }else{
        IncSegs<-rbind(IncSegs, incSeg)
      }
    }
  }
  print(proc.time()-ptm)
  
  # Above we generated the network using the initial set of POIs; here we crosswalk over the old COMIDs to the new
  IncSegs_fin<- IncSegs %>% left_join(st_drop_geometry(poiFix_sf) %>% select(COMID, new_COMID), by = c("POI_ID" = "COMID")) %>% 
    mutate(POI_ID = ifelse(is.na(new_COMID), POI_ID, new_COMID)) %>% filter(!POI_ID %in% unConPOIs$COMID) %>%
    select(-new_COMID)
  
  #poiFix_DS <- poiFix_sf %>% select(new_COMID, ToCom) %>% filter(!is.na(ToCom)) %>% st_drop_geometry() %>% distinct() %>%
  #  inner_join(IncSegs_fin %>% select(COMID, POI_ID), by = c("ToCom" = "COMID"))
  
  # Bring over NHDv2 attributes
  ncombined <- nhd %>% inner_join(select(IncSegs_fin, c(COMID, POI_ID)), by = "COMID")
  
  # Write out nsegments composed of nhdflowlines
  write_sf(ncombined, out_gpkg, nsegment_raw)

  # Dissolve flowlines to aggregated segments
  nsegments <- ncombined %>% group_by(POI_ID) %>% arrange(desc(LevelPathI), desc(Hydroseq)) %>%
    mutate(TT_Hours = (LENGTHKM * 3280.84) / (VA_MA * 3600)) %>% 
    summarize(TotalLength = sum(LENGTHKM),TotalDA = max(TotDASqKM), HW = max(StartFlag), TT = sum(PathTimeMA),
            do_union=FALSE) %>% st_cast("MULTILINESTRING")  %>%
    inner_join(st_drop_geometry(hucgagessegs) %>% select(COMID, Hydroseq, DnHydroseq),   by = c("POI_ID" = "COMID"))
  
  # produce a short data frame for populating TO_POI for downstream segment
  to_from <- segPOIs %>% inner_join(st_drop_geometry(ncombined), by = c("DnHydroseq" = "Hydroseq")) %>%
    select(COMID.x, Hydroseq, DnHydroseq, POI_ID) %>% distinct()

  # Add To_POI_ID to dissolved segments
  nsegments_fin <- nsegments %>% left_join(to_from, by = c("POI_ID" = "COMID.x")) %>%
                                    mutate(To_POI_ID = POI_ID.y) %>%
    select(POI_ID, TotalLength, TotalDA, HW, TT, Hydroseq.x, To_POI_ID) %>%
    left_join(st_drop_geometry(poiFix_sf) %>% select(new_COMID, ToCom) %>% distinct(), by = c("POI_ID" = "new_COMID")) %>%
    left_join(st_drop_geometry(ncombined) %>% select(COMID, POI_ID), by = c("ToCom" = "COMID")) %>%
    mutate(To_POI_ID = ifelse(!is.na(ToCom), POI_ID.y, To_POI_ID)) %>% select (-c(ToCom, POI_ID.y)) %>%
    rename(Hydroseq = Hydroseq.x, POI_ID = POI_ID.x)

  # Write out dissolved segments
  write_sf(nsegments_fin, out_gpkg, n_segments)
} else {
  nsegmentraw <- read_sf(out_gpkg, nsegment_raw)
  nsegments <- read_sf(out_gpkg, n_segments)
  tmpPOIs <- read_sf(out_gpkg, pois_all)
}

sub <- nhd %>% filter(COMID %in% finalPOIs$COMID)
print (dim(sub[sub$AreaSqKM == 0,]))

siteAtt(hydReg)
#knitr::purl("NHD_navigate.Rmd")
```