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)

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

nhd_rds <- fix_headwaters(staged_nhd$flowline, gsub("flowline.rds", "flowline_update.rds", staged_nhd$flowline),
                          new_atts = data_paths$new_nhdp_atts, nhdpath = data_paths$nhdplus_gdb)

nhdplusTools_data_dir(dir = data_paths$nhdplus_dir)
atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.rds"))
```

```{r huc12 POIs}
#  Derive or load HUC12 POIs
if(needs_layer(nav_gpkg, nav_poi_layer)) {
  
  if (needs_layer(nav_gpkg, nhd_flowline)){

    nhd <- readRDS(nhd_rds)
    
    # Subset NHD by VPU
    nhd <- subset_vpu(nhd, VPU) 
    
    # accidental duplocates?
    nhd <- nhd %>%
      group_by(COMID) %>%
      filter(n() == 1) %>%
      ungroup()

    # Filter and write dendritic/non-coastal subset to gpkg
    # This will be iterated over to produce the final network after POIs identified
    non_dend <- unique(unlist(lapply(filter(nhd, TerminalFl == 1 & TotDASqKM < min_da_km)
                                     %>% pull(COMID), NetworkNav, st_drop_geometry(nhd))))
      
    # Add fields to note dendritic and POI flowlines
    nhd <- nhd %>% 
      mutate(dend = ifelse(!COMID %in% non_dend, 1, 0),
             poi = ifelse(!COMID %in% non_dend & TotDASqKM >= min_da_km, 1, 0)) 
    
    write_sf(nhd, nav_gpkg, nhd_flowline)
  } else {
     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) %>% 
    slice(1)

  # 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, nav_gpkg, nav_poi_layer)
  tmp_POIs <- huc12_POIs
} else {
  # Load HUC12 POIs as the tmpPOIs if they already exist
  tmp_POIs <- read_sf(nav_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),
             nav_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, nav_gpkg, nav_poi_layer)
} else {
  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
}

mapview(filter(tmp_POIs, 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()) 
  
   # 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),
             nav_gpkg, "unassigned_TE")
  }
  
  # Write out geopackage layer representing POIs for given theme
  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
} else {
  # Load TE POIs if they already exist
  tmp_POIs <- read_sf(nav_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(!"Type_Term" %in% colnames(tmp_POIs)) {
  
  # 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, nav_gpkg, nav_poi_layer)
} else {
  tmp_POIs <- read_sf(nav_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(nhd, nav_gpkg, nhd_flowline)
  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
} else {
  nhd <- read_sf(nav_gpkg, nhd_flowline)
  tmp_POIs <- read_sf(nav_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, nav_gpkg, nav_poi_layer)
} else {
  # Load NID POIs if they already exist
  tmp_POIs <- read_sf(nav_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 <- readRDS("data/NHDPlusNationalData/nhdplus_waterbodies.rds") %>%
    mutate(FTYPE = as.character(FTYPE)) %>%
    filter(FTYPE %in% c("LakePond", "Reservoir") &
             AREASQKM >= (min_da_km/2) &
           COMID %in% filter(nhd, minNet == 1)$WBAREACOMI) %>%
    st_sf()

  # Write out waterbodies
  write_sf(st_transform(WBs_VPU, 4269), 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,
                        TotDASqKM >= min_da_km) %>%
    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) %>%
    slice(n = 1)

  tmp_POIs <- POI_creation(filter(wbout_COMIDs, !COMID %in% wbin_COMIDs$COMID), filter(nhd, poi == 1), "WBOut") %>%
    addType(., tmp_POIs, "WBOut")
  
  tmp_POIs <- POI_creation(filter(wbin_COMIDs, !COMID %in% wbout_COMIDs$COMID), filter(nhd, poi == 1), "WBIn") %>%
    addType(., tmp_POIs, "WBIn")

  write_sf(nhd, nav_gpkg, nhd_flowline)
  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
} else {
  tmp_POIs <- read_sf(nav_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 = "blue")
```


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}
if(all(is.na(tmp_POIs$Type_Elev))) {

  # 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, LENGTHKM, MAXELEVSMO, MINELEVSMO, AreaSqKM, TotDASqKM), by = "COMID") 

  # Run after incremental segments
  elev_pois_init <- 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 that need splitting
    filter(max(MAXELEVSMO) - min(MINELEVSMO) > (elev_diff * 100)) %>%
    # Order from upstream to downstream
    arrange(desc(Hydroseq)) %>%
    # Get cumulative sums of median elevation, Estimate number of times segments need splitting
    # csum_length = cumulative length US to DS along segment, cumsum_elev = cumulative US to DS elev diff along segment
    # iter = estimated number of times segment may need to be split (modified below if need be)
    mutate(csum_length = cumsum(LENGTHKM), cumsum_elev = cumsum(MAXELEVSMO - MINELEVSMO), iter = round(elev_diff_seg / (elev_diff * 100))) %>%
    # Only look to split segments based on elevation longer than 10 km
    filter(sum(LENGTHKM) > (split_meters/1000)) %>%
    ungroup() 

  if(nrow(elev_pois_init > 0)){

    # Iterate through segs that need splitting
    for (i in c(1:max(elev_pois_init$iter))){

      # first split
      elev_pois_init_pre <- elev_pois_init %>%
        group_by(POI_ID) %>%
        # Split threshold is half-value between min and max (median)
        mutate(split_elev = (max(MAXELEVSMO) - min(MINELEVSMO)) / 2) %>%
        ungroup()

      # Test if split_elev will be adquate for all cases
      elev_pois_init_post <- elev_pois_init_pre %>%
        group_by(POI_ID) %>%
        filter(split_elev > cumsum_elev) %>%
        mutate(cumsum_elev = cumsum(MAXELEVSMO - MINELEVSMO)) %>%
        # This will determine cases where a second split is needed
        filter(max(cumsum_elev) > (elev_diff * 100)) %>%
        ungroup()

      # Test for if we need to split iter ==2 again
      if(nrow(elev_pois_init_post) > 0){
        # Re-cacl the iter based on results above
        elev_pois_init <- mutate(elev_pois_init_pre,
                                 iter = ifelse(POI_ID %in% elev_pois_init_post$POI_ID, iter, 1))

        # Re-caculate split if  more than one split required
        elev_pois_init <- elev_pois_init %>%
          group_by(POI_ID) %>%
          mutate(split_elev = (max(MAXELEVSMO) - min(MINELEVSMO)) / (unique(iter) + 1)) %>%
          ungroup()

      } else {
        print ("no further divisions needed")
        # If no further division needed, set all to 1
        elev_pois_init <- mutate(elev_pois_init_pre, iter = 1)
      }

      # For reach iteration
      new_POIs <- elev_pois_init %>%
        group_by(POI_ID) %>%
        # select for new POIs that meet criteria; avoid splits in small, steep HW cats
        filter(split_elev < cumsum_elev, TotDASqKM > 1, csum_length > ((split_meters/1000) / 2)) %>% #, iter >= i) #%>%
        # for each poi get the most upstream fl that matches the above conditions
        filter(Hydroseq == max(Hydroseq)) %>%
        mutate(Type_Elev = 1) %>%
        ungroup()

      # Create new POI data frame
      if(exists("elev_POIs")){
        elev_POIs <- rbind(elev_POIs, new_POIs)
      } else {
        elev_POIs <- new_POIs
      }

      # create new set for next iteration
      elev_pois_init <- elev_pois_init_post %>%
        filter(POI_ID %in% new_POIs$POI_ID & iter != i)

      # IF further case exist, re-calc the cum. elev change
      if(nrow(elev_pois_init) > 0){
        elev_pois_init <- elev_pois_init %>%
          group_by(POI_ID) %>%
          mutate(cumsum_elev = cumsum(MAXELEVSMO - MINELEVSMO)) %>%
          filter(max(cumsum_elev) > 50000) %>%
          ungroup()
      }
    }

    print(nrow(elev_POIs))

    tmp_POIs <- POI_creation(select(elev_POIs, COMID, Type_Elev), nhd, "Elev") %>%
      addType(., tmp_POIs, "Elev")
  } else {
    tmp_POIs <- mutate(tmp_POIs, Type_Elev = NA)
  }

  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
} else {
  tmp_POIs <- read_sf(nav_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)) %>%
    inner_join(select(atts, COMID, VA_MA, TOTMA, LENGTHKM, WBAREACOMI, WBAreaType, FTYPE, TotDASqKM, AreaSqKM), by = "COMID") 

  tt_pois_init <- 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) > tt_diff) %>%
    # 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 = round(sum(FL_tt_hrs) / tt_diff)) %>%
    # Only look to split segments based on travel time longer than 10 km
    filter(sum(LENGTHKM) > (split_meters/1000)) %>%
    ungroup() 

  for (i in c(1:max(tt_pois_init$iter))){

    # first split
    tt_pois_init_pre <- tt_pois_init %>%
      group_by(POI_ID) %>%
      mutate(split_tt = (sum(FL_tt_hrs) / 2)) %>%
      ungroup()
    
    # Test if split_elev will be adquate for all cases
    tt_pois_init_post <- tt_pois_init_pre %>%
      group_by(POI_ID) %>%
      filter(split_tt > cumsum_tt) %>%
      # Re-calculate cumulative travel time once flowlines affected by first split
      #    taken out
      mutate(cumsum_tt = cumsum(FL_tt_hrs)) %>%
      # This will determine cases where a second split is needed
      filter(max(cumsum_tt) > tt_diff) %>%
      ungroup()
    
    # Test for if we need to split iter ==2 again
    if(nrow(tt_pois_init_post) > 0){
      # If some of the FL (iter > 1) need a second split
      elev_pois_init <- mutate(tt_pois_init_pre,
                               iter = ifelse(POI_ID %in% tt_pois_init_post$POI_ID, iter, 1))

      # Re-caculate split if  more than than one split required
      tt_pois_init <- tt_pois_init %>%
        group_by(POI_ID) %>%
        mutate(split_tt = (sum(FL_tt_hrs) / (unique(iter) + 1))) %>%
        ungroup()

    } else {
      print ("no further divisions needed")
      # If no further division/splits needed, set all iter values to 1
      tt_pois_init <- mutate(tt_pois_init_pre, iter = 1)
    }
    
    # For reach iteration, generate new P OIs
    new_POIs <- tt_pois_init %>%
      group_by(POI_ID) %>%
      # select for new POIs that meet criteria, put some size criteria within
      filter(split_tt < cumsum_tt, TotDASqKM > 1, csum_length > ((split_meters/1000) / 2)) %>%
      # for each poi get the most upstream fl that matches
      filter(Hydroseq == max(Hydroseq)) %>%
      mutate(Type_Travel = 1) %>%
      ungroup()
    
    # Create new POI data frame
    if(exists("tt_POIs")){
      tt_POIs <- rbind(tt_POIs, new_POIs)
    } else {
      tt_POIs <- new_POIs
    }
    
    # create new set for next iteration
    tt_pois_init <- tt_pois_init_post %>%
      filter(POI_ID %in% new_POIs$POI_ID & iter != i)

    # IF further case exist, re-calc the cum. elev change
    if(nrow(tt_pois_init) > 0){
      tt_pois_init <- tt_pois_init %>%
        group_by(POI_ID) %>%
        mutate(cumsum_tt = cumsum(FL_tt_hrs)) %>%
        filter(max(cumsum_tt) > tt_diff) %>%
        ungroup()
    }

  }

  tmp_POIs <- POI_creation(select(tt_POIs, COMID, Type_Travel), nhd, "Travel") %>%
    addType(., tmp_POIs, "Travel")

  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
}else{
  tmp_POIs <- read_sf(nav_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(nav_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, nav_gpkg, poi_moved_layer)
    write_sf(xWalk, nav_gpkg, poi_xwalk_layer)
    # write out final POIs
    write_sf(final_POIs, nav_gpkg, pois_all_layer)
  } else {
    # If no fixes designate as NA
    poi_fix <- NA
    # write out final POIs
    write_sf(tmp_POIs, nav_gpkg, pois_all_layer)

  }

} else {
  # Need all three sets for attribution below
  final_POIs <- read_sf(nav_gpkg, pois_all_layer)
  new_POIs <- if(!needs_layer(nav_gpkg, poi_moved_layer)) read_sf(nav_gpkg, poi_moved_layer) else (NA)
  xWalk <- if(!needs_layer(nav_gpkg, poi_xwalk_layer)) read_sf(nav_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(nav_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(nhd_Final, nav_gpkg, nhd_flowline)
  write_sf(nsegments, nav_gpkg, nsegments_layer)
} else {
  # Read in NHDPlusV2 flowline simple features and filter by vector processing unit (VPU)
  final_POIs <- read_sf(nav_gpkg, pois_all_layer)
  nhd_Final <- read_sf(nav_gpkg, nhd_flowline)
  nsegments <- read_sf(nav_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)) {

  #1 Move POIs downstream by category
  out_gages <- POI_move_down(nav_gpkg, final_POIs, nsegments, filter(nhd_Final, !is.na(POI_ID)), "Type_Gages", .05)
  out_HUC12 <- POI_move_down(nav_gpkg, out_gages$allPOIs, out_gages$segs, out_gages$FL, "Type_HUC12", .10)

  # Convert empty strings to NA for handling within R
  out_HUC12$allPOIs <- mutate_all(out_HUC12$allPOIs, list(~na_if(.,"")))

  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_HUC12$allPOIs, nav_gpkg, final_poi_layer)
  write_sf(nhd_Final, nav_gpkg, nhd_flowline)
  write_sf(out_HUC12$segs, nav_gpkg, nsegments_layer)
  
  nsegments <- out_HUC12$segs
  final_POIs <- out_HUC12$allPOIs
} else {
  final_POIs <- read_sf(nav_gpkg, final_poi_layer)
  nhd_Final <- read_sf(nav_gpkg, nhd_flowline)
  nsegments <- read_sf(nav_gpkg, nsegments_layer)
}

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