02_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 setup_rmd, echo=FALSE, cache=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=6, 
  fig.height=4,
  cache=FALSE)
```

```{r setup}
# Load custom functions
source("R/utils.R")
source("R/NHD_navigate.R")
source("R/hyrefactor_funs.R")

# Load Configuration of environment
source("R/config.R")

gages <- read_sf("cache/Gages_info.gpkg", "Gages")

nc <- RNetCDF::open.nc(data_paths$nwm_network)

nwm_gages <- data.frame(
  comid = RNetCDF::var.get.nc(nc, "link"), 
  gage_id = RNetCDF::var.get.nc(nc, "gages")) %>%
  mutate(gage_id = gsub(" ", "", gage_id)) %>%
  mutate(gage_id = ifelse(gage_id == "", NA, gage_id))

RNetCDF::close.nc(nc)

# need some extra attributes
atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.rds"))
# atts <- select(atts, COMID, MAXELEVSMO, MINELEVSMO, VA_MA, TOTMA, WBAreaType)

```

```{r huc12 POIs}
#  Derive or load HUC12 POIs
if(needs_layer(temp_gpkg, nav_poi_layer)) {
  
   nhd <- read_sf(nav_gpkg, nhd_flowline)
   try(nhd <- select(nhd, -c(minNet, WB, struct_POI, struct_Net, POI_ID)))
  
  # HUC02 includes some 
  if(vpu == "02"){
    grep_exp <-"^02|^04"
  } else if (vpu == "08") {
    grep_exp <- "^03|^08"
  } else {
    grep_exp <- paste0("^", substr(vpu, start = 1, stop = 2))
  }
  
  # Join HUC12 outlets with NHD
  HUC12_COMIDs <- read_sf(data_paths$hu12_points_path, "hu_points") %>% 
    filter(grepl(grep_exp, .data$HUC12)) %>%
    select(COMID, HUC12) %>%
    # Remove this when HUC12 outlets finished
    group_by(COMID) %>% 
    filter(row_number() == 1) %>%
    ungroup()

  # Create POIs - some r05 HUC12 POIs not in R05 NHD
  huc12_POIs <- POI_creation(st_drop_geometry(HUC12_COMIDs), filter(nhd, poi == 1), "HUC12")

  # Write out geopackage layer representing POIs for given theme
  write_sf(huc12_POIs, temp_gpkg, nav_poi_layer)
  tmp_POIs <- huc12_POIs
} else {
  # Load HUC12 POIs as the tmpPOIs if they already exist
  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer) 
  nhd <- read_sf(nav_gpkg, nhd_flowline)
}

mapview(filter(tmp_POIs, Type_HUC12 != ""), layer.name = "HUC12 POIs", col.regions = "red") 
```

```{r streamgage POIs}
if(all(is.na(tmp_POIs$Type_Gages))) { 

  # Previously identified streamgages within Gage_Selection.Rmd
  streamgages_VPU <- gages %>%
    rename(COMID = comid) %>%
    filter(COMID %in% nhd$COMID) %>%
    st_drop_geometry() %>%
    switchDiv(., nhd) 
  
  streamgages <- streamgages_VPU %>% 
    group_by(COMID) %>%
    # If multiple gages per COMID, pick one with highest rank as rep POI_ID
    filter(gage_score == max(gage_score), !is.na(drain_area)) %>%
    ungroup() %>%
    select(COMID, site_no)
    
  # Derive GAGE POIs; use NHD as we've alredy filtered by NWIS DA in the Gage selection step
  tmp_POIs <- POI_creation(streamgages, nhd, "Gages") %>%
    addType(., tmp_POIs, "Gages")
  
  # As a fail-safe, write out list of gages not assigned a POI
  if(nrow(filter(streamgages_VPU, !site_no %in% tmp_POIs$Type_Gages)) > 0) {
    write_sf(filter(streamgages_VPU, !site_no %in% tmp_POIs$Type_Gages),
             temp_gpkg, "unassigned_gages")
    
    # Start documenting gages that are dropped out; these gages have no mean daily Q
    gage_document(filter(streamgages_VPU, !site_no %in% tmp_POIs$Type_Gages) %>%
                    select(site_no),
              "Gages_info", paste0("VPU ", vpu, "; low gage score"))
    
  }
  
  # Write out geopackage layer representing POIs for given theme
  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
} else {
  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
}

mapview(filter(tmp_POIs, !is.na(Type_Gages)), layer.name = "Streamgage POIs", col.regions = "blue") 
```

```{r TE POIs}
if(all(is.na(tmp_POIs$Type_TE))) {
  
  if(vpu == "08"){
    nhd$VPUID <- "08"
  } else {
    nhd$VPUID <- substr(nhd$RPUID, 1, 2)
  }
  
  # Read in Thermoelectric shapefile
  TE_COMIDs <- read_sf(data_paths$TE_points_path, "2015_TE_Model_Estimates_lat.long_COMIDs") %>%
    mutate(COMID = as.integer(COMID)) %>%
    inner_join(., select(st_drop_geometry(nhd), COMID, VPUID), by = "COMID") %>%
    filter(grepl(paste0("^", substr(vpu, 1, 2), ".*"), .data$VPUID), COMID > 0) %>%
    switchDiv(., nhd) %>%
    group_by(COMID) %>%
    summarize(EIA_PLANT = paste0(unique(EIA_PLANT_), collapse = " "), count = n()) %>%
    ungroup()
  
   # Derive TE POIs
  tmp_POIs <- POI_creation(st_drop_geometry(TE_COMIDs), filter(nhd, dend == 1), "TE") %>%
    addType(., tmp_POIs, "TE")

  # As a fail-safe, write out list of TE plants not assigned a POI
  if(nrow(filter(TE_COMIDs, !COMID %in% tmp_POIs$COMID)) > 0) {
    write_sf(filter(TE_COMIDs, !COMID %in% tmp_POIs$COMID),
             temp_gpkg, "unassigned_TE")
  }
  
  # Write out geopackage layer representing POIs for given theme
  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
} else {
  # Load TE POIs if they already exist
  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
}

mapview(filter(tmp_POIs, Type_TE != ""), layer.name = "TE Plant POIs", col.regions = "blue") 
```

```{r Terminal POIs}
#  Derive or load Terminal POIs ----------------------
if(all(is.na(tmp_POIs$Type_Term))) {
  
  # Terminal POIs on levelpaths with upstream POIs
  term_paths <- filter(st_drop_geometry(nhd), Hydroseq %in% filter(nhd, COMID %in% tmp_POIs$COMID)$TerminalPa)
  
  # Non-POI levelpath terminal pois, but meet size criteria
  terminal_POIs <- st_drop_geometry(nhd) %>%
    filter(Hydroseq == TerminalPa | (toCOMID == 0 | is.na(toCOMID) | !toCOMID %in% COMID)) %>%
    filter(!COMID %in% term_paths$COMID, TotDASqKM >= 20) %>%
    bind_rows(term_paths) %>%
    # Use level path identifier as Type ID
    select(COMID, LevelPathI)

  tmp_POIs <- POI_creation(terminal_POIs, nhd, "Term") %>%
    addType(., tmp_POIs, "Term")

  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
} else {
  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
  nhd <- read_sf(nav_gpkg, nhd_flowline)
}

mapview(filter(tmp_POIs, !is.na(Type_Term)), layer.name = "Terminal POIs", col.regions = "blue") 
```

```{r Confluence POIs}
# Derive POIs at confluences where they are absent ----------------------
if(all(is.na(tmp_POIs$Type_Conf))) {

  # Navigate upstream from each POI and determine minimally-sufficient network between current POI sets
  up_net <- unique(unlist(lapply(unique(tmp_POIs$COMID), NetworkNav, nhd))) 
  finalNet <- unique(NetworkConnection(up_net, st_drop_geometry(nhd))) 

  # Subset NHDPlusV2 flowlines to navigation results and write to shapefile
  nhd <- mutate(nhd, minNet = ifelse(COMID %in% finalNet, 1, 0))
  
  # Create new confluence POIs
  conf_COMIDs <- st_drop_geometry(filter(nhd, minNet == 1)) %>%
    # Downstream hydrosequence of 0 indicates
    #   the flowline is terminating or
    #   leaving the domain, so they
    #   are excluded from this process
    filter(DnHydroseq > 0) %>%
    group_by(DnHydroseq) %>%
    filter(n()> 1) %>%
    mutate(Type_Conf = LevelPathI) %>%
    ungroup() %>% 
    select(COMID, Type_Conf)

  tmp_POIs <- POI_creation(conf_COMIDs, filter(nhd, minNet == 1), "Conf") %>%
    addType(., tmp_POIs, "Conf")
  
  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
} else {
  nhd <- read_sf(nav_gpkg, nhd_flowline)
  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
}

mapview(filter(tmp_POIs, !is.na(Type_Conf)), layer.name = "Confluence POIs", col.regions = "blue") 
```

```{r NID POIs}
#  Derive or load NID POIs ----------------------
if(all(is.na(tmp_POIs$Type_NID))) {

  # Read in NID shapefile
  NID_COMIDs <- read.csv(file.path(data_paths$NID_points_path, "NID_attributes_20170612.txt"), 
                         stringsAsFactors = FALSE) %>%
    filter(ONNHDPLUS == 1, FlowLcomid %in% filter(nhd, dend ==1)$COMID) %>%
    switchDiv(., nhd) %>%
    group_by(FlowLcomid) %>%
    summarize(Type_NID = paste0(unique(NIDID), collapse = " ")) 
  
  # Derive NID POIs
  tmp_POIs <- POI_creation(NID_COMIDs, filter(nhd, dend ==1), "NID") %>%
    addType(., tmp_POIs, "NID", bind = FALSE)

  # Write out geopackage layer representing POIs for given theme
  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
} else {
  # Load NID POIs if they already exist
  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
}

mapview(filter(tmp_POIs, Type_NID != ""), layer.name = "NID POIs", col.regions = "blue") 
```

```{r waterbody POIs}
#  Derive or load Waterbody POIs ----------------------
if(all(is.na(tmp_POIs$Type_WBOut))) {

  # Read in waterbodies
  WBs_VPU_sel <- readRDS("data/NHDPlusNationalData/nhdplus_waterbodies.rds") %>%
    mutate(FTYPE = as.character(FTYPE)) %>%
    dplyr::filter(COMID %in% filter(nhd, minNet == 1)$WBAREACOMI) %>%
    st_sf()
  
  wb_lst <- st_drop_geometry(WBs_VPU_sel) %>%
    dplyr::filter(FTYPE %in% c("LakePond", "Reservoir") & AREASQKM >= (min_da_km/2)) %>%
    dplyr::select(COMID)
  
  # if RESOPS cw exists will add WBs (that are nhd with minNet == 1) to the WBPOIS  
  resops <- data_paths$resops_NID_CW
  if (file.exists(resops)){
    resops_wb_df <- read.csv(resops) %>%
      dplyr::filter(!is.na(WBAREACOMI)) %>%
      dplyr::filter(WBAREACOMI %in% WBs_VPU_sel$COMID) %>%
      dplyr::select(DAM_ID, DAM_NAME, WBAREACOMI)
      
    resops_wb_lst <- resops_wb_df %>%
      dplyr::select(COMID = WBAREACOMI) %>%
      distinct()
      
    if(nrow(resops_wb_lst) > 0){
      wb_lst <- rbind(wb_lst, resops_wb_lst) %>%
        distinct()
    }
    wb_lst <- rbind(wb_lst, resops_wb_lst) %>%
      distinct()
  }
  
  WBs_VPU <- WBs_VPU_sel %>%
    dplyr::filter(COMID %in% wb_lst$COMID)
  
  # Write out waterbodies
  write_sf(st_transform(WBs_VPU, sf::st_crs(nhd)), nav_gpkg, WBs_layer)

  # Segments that are in waterbodies
  nhd <- mutate(nhd, WB = ifelse(WBAREACOMI > 0 & WBAREACOMI %in% WBs_VPU$COMID, 1, 0))

  # Create waterbody outlet POIs
  wbout_COMIDs <- filter(nhd, dend == 1 & WB == 1) %>%
    group_by(WBAREACOMI) %>%
    slice(which.min(Hydroseq)) %>%
    switchDiv(., nhd) %>%
    select(COMID, WBAREACOMI) 

  # Create waterbody inlet POIs
  wbin_COMIDs <- filter(nhd, WB == 0, 
                        DnHydroseq %in% filter(nhd, WB == 1)$Hydroseq) %>%
    select(-WBAREACOMI) %>%
    switchDiv(., nhd) %>%
    inner_join(st_drop_geometry(filter(nhd, minNet == 1)) %>%
                 select(Hydroseq, WBAREACOMI), by = c("DnHydroseq" = "Hydroseq")) %>%
    select(COMID, WBAREACOMI) %>%
    group_by(COMID) %>%
    filter(row_number() == 1) %>%
    ungroup()
  
  tmp_POIs <- POI_creation(filter(st_drop_geometry(wbout_COMIDs), !COMID %in% wbin_COMIDs$COMID), 
                           filter(nhd, poi == 1), "WBOut") %>%
    addType(., tmp_POIs, "WBOut")
  
  tmp_POIs <- POI_creation(filter(st_drop_geometry(wbin_COMIDs), !COMID %in% wbout_COMIDs$COMID), 
                           filter(nhd, poi == 1), "WBIn") %>%
    addType(., tmp_POIs, "WBIn")
  
  # sel the resops POIS, write to csv
  if(wb_lst > 0){
    resops_POIs_df <- st_drop_geometry(tmp_POIs) %>%
      dplyr::select(COMID, Type_WBIn, Type_WBOut) %>%
      dplyr::filter(!is.na(Type_WBIn) | !is.na(Type_WBOut)) %>%
      dplyr::mutate(WBAREACOMI = if_else(is.na(Type_WBIn), as.integer(Type_WBOut), as.integer(Type_WBIn))) %>%
      dplyr::filter(WBAREACOMI %in% resops_wb_lst$COMID) %>%
      dplyr::left_join(resops_wb_df, by = "WBAREACOMI")
    write.csv(resops_POIs_df, file.path("data/reservoir_Data", paste0("resops_POIs_",vpu,".csv")))
  }
    
  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
} else {
  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
  nhd <- read_sf(nav_gpkg, nhd_flowline)
}

mapview(filter(tmp_POIs, !is.na(Type_WBIn)), layer.name = "Waterbody inlet POIs", col.regions = "blue") +
mapview(filter(tmp_POIs, !is.na(Type_WBOut)), layer.name = "Waterbody outlet POIs", col.regions = "red")
```


This chunk breaks up potential aggregated segments where the elevation range of the are contributing to the segment exceeds 500-m
```{r Elevation Break POIs}
# derive incremental segments from POIs
inc_segs <- segment_increment(filter(nhd, minNet == 1), 
                              filter(st_drop_geometry(nhd),
                                     COMID %in% tmp_POIs$COMID, COMID %in% filter(nhd, minNet == 1)$COMID)) %>%
  # bring over VAA data
  inner_join(select(atts, COMID, DnHydroseq, VA_MA, TOTMA, LENGTHKM, MAXELEVSMO, MINELEVSMO, WBAREACOMI, WBAreaType, FTYPE,
                    AreaSqKM, TotDASqKM), by = "COMID") 

# Build dataframe for creation of points based on breaks
elev_pois_split <- inc_segs %>%
  group_by(POI_ID) %>%
  # Get elevation info
  mutate(MAXELEVSMO = na_if(MAXELEVSMO, -9998), MINELEVSMO = na_if(MINELEVSMO, -9998), 
         elev_diff_seg = max(MAXELEVSMO) - min(MINELEVSMO)) %>%
  # Filter out to those incremental segments that need splitting above the defined elevation threshold
  filter((max(MAXELEVSMO) - min(MINELEVSMO)) > (elev_diff * 100)) %>%
  # Do not aggregate on fake lfowlines
  filter(AreaSqKM > 0, TotDASqKM > max_elev_TT_DA) %>% 
  # Order from upstream to downstream
  arrange(desc(Hydroseq)) %>%
  # Get attributes used for splitting
  # csum_length = cumulative length US to DS along segment, cumsum_elev = cumulative US to DS elev diff along segment
  # inc_elev = elevation range of each segment
  # iter = estimated number of times we'd need to split existing agg flowlines, or number of rows in set
  mutate(csum_length = cumsum(LENGTHKM),
         inc_elev = cumsum(MAXELEVSMO - MINELEVSMO),
         #nc_elev_diff = c(inc_elev[1], (inc_elev - lag(inc_elev))[-1]),
         iter = min(round(max(inc_elev) / (elev_diff * 100)), n()),
         elev_POI_ID = NA) %>%
  # Remove segments we can't break
  filter(iter > 0, n() > 1) %>%
  # Track original iteration number
  mutate(orig_iter = iter) %>%
  ungroup() 

if(nrow(elev_pois_split) > 0) {
  
  # For reach iteration
  elev_POIs <- do.call("rbind", lapply(unique(elev_pois_split$POI_ID), split_elev, 
                                      elev_pois_split, elev_diff*100, max_elev_TT_DA)) %>%
    filter(!elev_POI_ID %in% tmp_POIs$COMID, COMID == elev_POI_ID) %>%
    mutate(Type_Elev = 1) %>%
    select(COMID, Type_Elev) %>%
    unique()
  
  if(nrow(elev_POIs) > 0){
    tmp_POIs <- POI_creation(select(elev_POIs, COMID, Type_Elev), nhd, "Elev") %>%
      addType(., tmp_POIs, "Elev")
  } 
  

} else {
  
  tmp_POIs$Type_Elev <- rep(NA, nrow(tmp_POIs))

  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
    
  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer) 
  nhd <- read_sf(nav_gpkg, nhd_flowline)
}

if(!all(is.na(tmp_POIs$Type_Elev)))
  mapview(filter(tmp_POIs, Type_Elev == 1), layer.name = "Elevation break POIs", col.regions = "blue") 
```

```{r Time of travel POIs}
if(all(is.na(tmp_POIs$Type_Travel))) {
  
  # derive incremental segments from POIs
  inc_segs <- segment_increment(filter(nhd, minNet == 1), 
                                filter(st_drop_geometry(nhd),
                                       COMID %in% tmp_POIs$COMID, COMID %in% filter(nhd, minNet == 1)$COMID)) %>%
    # bring over VAA data
    inner_join(select(atts, COMID, DnHydroseq, VA_MA, TOTMA, LENGTHKM, MAXELEVSMO, MINELEVSMO, WBAREACOMI, WBAreaType, FTYPE,
                      AreaSqKM, TotDASqKM), by = "COMID") 

  # TT POIs
  tt_pois_split <- inc_segs %>%
    # Should we substitute with a very small value
    mutate(VA_MA = ifelse(VA_MA < 0, NA, VA_MA)) %>%
    mutate(FL_tt_hrs = (LENGTHKM * 3280.84)/ VA_MA / 3600 ) %>%
    group_by(POI_ID) %>%
    filter(sum(FL_tt_hrs) > travt_diff, TotDASqKM > max_elev_TT_DA) %>%
    # Sort upstream to downsream
    arrange(desc(Hydroseq)) %>%
    # Get cumulative sums of median elevation
    mutate(csum_length = cumsum(LENGTHKM), 
           cumsum_tt = cumsum(FL_tt_hrs), 
           iter = min(round(sum(FL_tt_hrs) / travt_diff), n()),
           tt_POI_ID = NA) %>%
    # Only look to split segments based on travel time longer than 10 km
    #filter(sum(LENGTHKM) > (split_meters/1000)) %>%
    filter(iter > 0, n() > 1) %>%
    # Track original iteration number
    mutate(orig_iter = iter) %>%
    ungroup() 

  # For reach iteration
  tt_POIs <- do.call("rbind", lapply(unique(tt_pois_split$POI_ID), split_tt, 
                                      tt_pois_split, travt_diff, max_elev_TT_DA)) %>%
    filter(!tt_POI_ID %in% tmp_POIs$COMID, COMID == tt_POI_ID) %>%
    mutate(Type_Travel = 1) %>%
    select(COMID, Type_Travel) %>%
    unique()
    
  tmp_POIs <- POI_creation(select(tt_POIs, COMID, Type_Travel), nhd, "Travel") %>%
    addType(., tmp_POIs, "Travel")
  
  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
}else{
  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
  nhd <- read_sf(nav_gpkg, nhd_flowline)
}

mapview(filter(tmp_POIs, Type_Travel == 1), layer.name = "Travel time break POIs", col.regions = "blue") 
```
  
```{r Final POIs}
# Derive final POI set ----------------------
if(needs_layer(temp_gpkg, pois_all_layer)) {

  unCon_POIs <- filter(st_drop_geometry(filter(nhd, minNet == 1)), COMID %in% tmp_POIs$COMID, AreaSqKM == 0)

  # If any POIs happened to fall on flowlines w/o catchment
  if (nrow(unCon_POIs) >0){
    # For confluence POIs falling on Flowlines w/o catchments, derive upstream valid flowline,
    poi_fix <- DS_poiFix(tmp_POIs, filter(nhd, dend == 1))
    new_POIs <- st_compatibalize(poi_fix$new_POIs, tmp_POIs)
    xWalk <- poi_fix$xWalk

    # POIs that didn't need to be moved
    tmp_POIs_fixed <- filter(tmp_POIs, !COMID %in% c(poi_fix$xWalk$oldPOI, poi_fix$xWalk$COMID))
    # bind together
    final_POIs <- bind_rows(tmp_POIs_fixed, select(poi_fix$new_POIs, -c(oldPOI, to_com)))
    # Write out fixes
    write_sf(new_POIs, temp_gpkg, poi_moved_layer)
    write_sf(xWalk, temp_gpkg, poi_xwalk_layer)
    # write out final POIs
    write_sf(final_POIs, temp_gpkg, pois_all_layer)
  } else {
    # If no fixes designate as NA
    poi_fix <- NA
    # write out final POIs
    write_sf(tmp_POIs, temp_gpkg, pois_all_layer)

  }

} else {
  # Need all three sets for attribution below
  final_POIs <- read_sf(temp_gpkg, pois_all_layer)
  new_POIs <- if(!needs_layer(temp_gpkg, poi_moved_layer)) read_sf(temp_gpkg, poi_moved_layer) else (NA)
  xWalk <- if(!needs_layer(temp_gpkg, poi_xwalk_layer)) read_sf(temp_gpkg, poi_xwalk_layer) else (NA)
  unCon_POIs <- filter(st_drop_geometry(filter(nhd, minNet == 1)), COMID %in% tmp_POIs$COMID, AreaSqKM == 0)
}
``` 
  
```{r Draft nsegment generation}
if(needs_layer(temp_gpkg, nsegments_layer)) {

  # Sort POIs by Levelpath and Hydrosequence in upstream to downstream order
  seg_POIs <-  filter(st_drop_geometry(nhd), COMID %in% tmp_POIs$COMID, COMID %in% filter(nhd, minNet == 1)$COMID)
  # derive incremental segments from POIs
  inc_segs <- segment_increment(filter(nhd, minNet == 1), seg_POIs)

  nhd_Final <- nhd %>%
    left_join(select(inc_segs, COMID, POI_ID), by = "COMID")

  # create and write out final dissolved segments
  nsegments_fin <- segment_creation(nhd_Final, xWalk)
  nhd_Final <- select(nhd_Final, -POI_ID) %>%
    left_join(select(st_drop_geometry(nsegments_fin$raw_segs), COMID, POI_ID), by = "COMID")
  nsegments <- nsegments_fin$diss_segs

  # Produce the minimal POIs needed to derive the network based on LP and terminal outlets
  strucFeat <- structPOIsNet(nhd_Final, final_POIs)
  nhd_Final <- nhd_Final %>%
    mutate(struct_POI = ifelse(COMID %in% strucFeat$struc_POIs$COMID, 1, 0),
           struct_Net = ifelse(COMID %in% strucFeat$structnet$COMID, 1, 0))

  write_sf(nsegments, temp_gpkg, nsegments_layer)
} else {
  # Read in NHDPlusV2 flowline simple features and filter by vector processing unit (VPU)
  final_POIs <- read_sf(temp_gpkg, pois_all_layer)
  nhd_Final <- read_sf(nav_gpkg, nhd_flowline)
  nsegments <- read_sf(temp_gpkg, nsegments_layer)
}

# Ensure that all the problem POIs have been taken care of
sub <- nhd_Final %>% filter(COMID %in% final_POIs$COMID)
print (paste0(nrow(sub[sub$AreaSqKM == 0,]), " POIs on flowlines with no local drainage contributions"))
```

```{r POI Collapse}
#  Load data
if(needs_layer(nav_gpkg, final_poi_layer)) {

  if(!"Type_Elev" %in% names(final_POIs)) final_POIs$Type_Elev <- rep(NA_real_, nrow(final_POIs))
  
  #1 Move POIs downstream by category
  out_HUC12 <- POI_move_down(temp_gpkg, final_POIs, nsegments, filter(nhd_Final, !is.na(POI_ID)), "Type_HUC12", .10)
  out_gages <- POI_move_down(temp_gpkg, out_HUC12$allPOIs, out_HUC12$segs, out_HUC12$FL, "Type_Gages", .05)
  
  # Convert empty strings to NA for handling within R
  out_gages$allPOIs <- mutate_all(out_gages$allPOIs, list(~na_if(.,"")))
  
  out_gages$allPOIs$snapped <- out_gages$allPOIs$COMID %in% new_POIs$COMID

  nhd_Final <- select(nhd_Final, -POI_ID) %>%
    left_join(st_drop_geometry(out_HUC12$FL) %>%
                select(COMID, POI_ID), by = "COMID")

  # Write out geopackage layer representing POIs for given theme
  write_sf(out_gages$allPOIs, nav_gpkg, final_poi_layer)
  write_sf(out_gages$segs, temp_gpkg, nsegments_layer)
  
  nsegments <- out_gages$segs
  final_POIs <- out_gages$allPOIs
} else {
  final_POIs <- read_sf(nav_gpkg, final_poi_layer)
  nhd_Final <- read_sf(nav_gpkg, nhd_flowline)
  nsegments <- read_sf(temp_gpkg, nsegments_layer)
}

mapview(final_POIs, layer.name = "All POIs", col.regions = "red") + 
  mapview(nsegments, layer.name = "Nsegments", col.regions = "blue")
```