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

# increase timeout for data downloads
options(timeout=600)

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

# Gages output from Gage_selection
gages <- read_sf(gage_info_gpkg, "Gages")

# need some extra attributes for a few POI analyses
atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.rds"))
```

```{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, dend, poi)), silent = TRUE)

  # Some NHDPlus VPUs include HUCs from other VPUs
  if(vpu_codes == "02"){
    grep_exp <-"^02|^04"
  } else if (vpu_codes == "08") {
    grep_exp <- "^03|^08"
  } else {
    grep_exp <- paste0("^", substr(vpu_codes, 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(!"Type_Gages" %in% names(tmp_POIs)) { 

  # Previously identified streamgages within Gage_Selection.Rmd
  streamgages_VPU <- gages %>%
    rename(COMID = comid) %>%
    filter(COMID %in% nhd$COMID) %>%
    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() 
    
  # Derive GAGE POIs; use NHD as we've already filtered by NWIS DA in the Gage selection step
  gages_POIs <- gage_POI_creation(tmp_POIs, streamgages, filter(nhd, poi == 1), combine_meters, reach_meas_thresh)
  tmp_POIs <- gages_POIs$tmp_POIs
  events <- rename(gages_POIs$events, POI_identifier = Type_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_codes, "; low gage score"))
    
  }
  
  # Write out events and outelts
  write_sf(events, temp_gpkg, split_layer)
  # 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)
  events <- read_sf(temp_gpkg, split_layer)
}

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

```{r TE POIs}
if(!"Type_TE" %in% names(tmp_POIs)) {
  
  if(vpu_codes == "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_codes, 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, poi == 1), "TE") %>%
    addType(., tmp_POIs, "TE", nexus = FALSE) 

  # 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 waterbody outlet POIs}
#  Derive or load Waterbody POIs ----------------------
if(!"Type_WBOut" %in% names(tmp_POIs)) {
  
  # Waterbodies sourced from NHD waterbody layer for entire VPU
  WBs_VPU_all <- filter(readRDS("data/NHDPlusNationalData/nhdplus_waterbodies.rds"), 
                        COMID %in% nhd$WBAREACOMI) %>%
    mutate(FTYPE = as.character(FTYPE)) %>%
    mutate(source = "NHDv2WB") %>%
    st_transform(crs)
  
  ref_WB <- read_sf(nav_gpkg, nhd_waterbody) %>%
    filter(COMID %in% nhd$WBAREACOMI) %>%
    mutate(source = "Ref_WB")

  # Waterbody list that strictly meet the size criteria
  wb_lst <- st_drop_geometry(WBs_VPU_all) %>%
    dplyr::filter(FTYPE %in% c("LakePond", "Reservoir") & AREASQKM >= (min_da_km/2)) %>%
    dplyr::select(COMID) %>%
    filter(!COMID == 166410600)

  if (file.exists(data_paths$resops_NID_CW)){
    # ResOpsUS locations in the VPU waterbody set
    resops_wb_df <- read.csv(data_paths$resops_NID_CW) %>%
      dplyr::filter(WBAREACOMI %in% WBs_VPU_all$COMID) %>%
      dplyr::select(DAM_ID, DAM_NAME, COMID = WBAREACOMI, FlowLcomid) %>%
      distinct()

    # Add ResOPsUS locations to waterbody list 
    if(nrow(resops_wb_df) > 0){
      wb_lst <- rbind(wb_lst, select(resops_wb_df, COMID)) %>%
        distinct()
    }
  }
  
  # Waterbodies that meet the size criteria and/or are in ResOpsUS datasets
  WBs_VPU <- WBs_VPU_all %>%
    dplyr::filter(COMID %in% wb_lst$COMID) 
  
  # Attribute flowlines that are in waterbodies
  nhd <- mutate(nhd, WB = ifelse(WBAREACOMI > 0 & WBAREACOMI %in% WBs_VPU$COMID, 1, 0))
  
  wb_layers <- wbout_POI_creaton(filter(nhd, WB == 1), WBs_VPU, data_paths, crs)
  if(!is.na(wb_layers$events) > 0) {events <- rbind(events, wb_layers$events)}
  WBs_VPU <- wb_layers$WBs
  
  wb_out_col <- wb_poi_collapse(wb_layers$POIs, filter(nhd, WB == 1), events)
  tmp_POIs <- wb_out_col$POIs
  if(!is.na(wb_out_col$events_ret)){
    write_sf(wb_out_col$events_ret, temp_gpkg, split_layer)
  }
  
  if(!all(is.na(wb_layers$events))){
    wb_events <- select(wb_layers$events, -c(id, offset)) %>%
      rename(POI_identifier = WBAREACOMI)
      
    # Write out events and outlets
    if(!needs_layer(temp_gpkg, split_layer)){
      events <- read_sf(temp_gpkg, split_layer) %>%
        rbind(st_compatibalize(wb_events,.)) %>%
        unique()
      write_sf(events, temp_gpkg, split_layer)
    } else {
      write_sf(wb_events, nav_gpkg, split_layer)
    }
  }
  
  ref_WB <- filter(ref_WB, COMID %in% WBs_VPU$COMID)
  
  WBs_VPU <- filter(WBs_VPU, !COMID %in% ref_WB$COMID) %>%
    group_by(GNIS_ID, GNIS_NAME, COMID, FTYPE, source) %>%
    summarize(do_union = T) %>%
    ungroup() %>% 
    st_make_valid() %>%
    sfheaders::sf_remove_holes(.) %>%
    mutate(member_comid = NA, 
           area_sqkm = as.numeric(st_area(.))/1000000) %>%
    st_compatibalize(., ref_WB) 
  
  ref_WB <- rbind(ref_WB, WBs_VPU)
    
  # Write out waterbodies
  write_sf(ref_WB, temp_gpkg, WBs_layer)

  # Write specific ResOpsUS data to a separate sheet
  if(nrow(wb_lst) > 0){
    resops_POIs_df <- st_drop_geometry(tmp_POIs) %>%
      dplyr::select(COMID, Type_WBOut) %>%
      dplyr::filter(!is.na(Type_WBOut)) %>%
      dplyr::inner_join(select(resops_wb_df, -FlowLcomid), by = c("Type_WBOut" = "COMID")) 
    
    write.csv(resops_POIs_df, file.path("data/reservoir_Data", 
                                        paste0("resops_POIs_", vpu_codes, ".csv")))
  }
    
  write_sf(nhd, nav_gpkg, nhd_flowline)
  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)
  ref_WB <- read_sf(temp_gpkg, WBs_layer)
  resops_POIs_df <- read.csv(file.path("data/reservoir_Data", paste0("resops_POIs_", VPU,".csv")))
}

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

```{r Terminal POIs}
#  Derive or load Terminal POIs ----------------------
if(!"Type_Term" %in% names(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 >= min_da_km) %>%
    bind_rows(term_paths) %>%
    # Use level path identifier as Type ID
    select(COMID, LevelPathI)

  tmp_POIs <- POI_creation(terminal_POIs, filter(nhd, poi == 1), "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(!"Type_Conf" %in% names(tmp_POIs)) {

  # 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, 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 waterbody inlet POIs}
#  Derive or load Waterbody POIs ----------------------
if(!"Type_WBIn" %in% names(tmp_POIs)) {

  wb_layers <- wbin_POIcreation(nhd, ref_WB, data_paths, crs, split_layer, tmp_POIs)
  wb_in_col <- wb_inlet_collapse(wb_layers$POIs, nhd, events)
  tmp_POIs <- wb_in_col$POIs
  
  if(!all(is.na(wb_layers$events))) {
    wb_inlet_events <- select(wb_layers$events, -c(id, offset, Hydroseq, ToMeas)) %>%
      rename(POI_identifier = WBAREACOMI)
      
    # Write out events and outlets
    if(!needs_layer(temp_gpkg, split_layer)){
      events <- read_sf(temp_gpkg, split_layer) %>%
        rbind(st_compatibalize(wb_inlet_events,.))
      write_sf(events, temp_gpkg, split_layer)
    } else {
      write_sf(wb_inlet_events, nav_gpkg, split_layer)
    }
  }
  
  tmp_POIs$nexus <- ifelse(is.na(tmp_POIs$nexus), FALSE, tmp_POIs$nexus)
  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") 
```


```{r NID POIs}
#  Derive or load NID POIs ----------------------
if(!"Type_NID" %in% names(tmp_POIs)) {

  # Read in NID shapefile
  NID_COMIDs <- read_sf(data_paths$NID_points_path, "Final_NID_2018") %>%
    st_drop_geometry() %>%
    filter(EROM != 0, FlowLcomid %in% filter(nhd, dend ==1)$COMID) %>%
    rename(COMID = FlowLcomid) %>%
    switchDiv(., nhd) %>%
    group_by(COMID) %>%
    summarize(Type_NID = paste0(unique(NIDID), collapse = " ")) 
  
  # Determine which NID facilitites have been co-located with ResOps
  res_ops_table <- read.csv(data_paths$resops_NID_CW) %>%
    mutate(new_WBARECOMI = ifelse(is.na(WBAREACOMI) & is.na(NHDAREA), HR_NHDPLUSID, WBAREACOMI)) %>%
    mutate(new_WBARECOMI = ifelse(is.na(new_WBARECOMI) & is.na(HR_NHDPLUSID), NHDAREA, new_WBARECOMI)) %>%
    filter(new_WBARECOMI %in% tmp_POIs$Type_WBOut, !is.na(new_WBARECOMI)) %>%
    select(NID_ID = NID_Source_FeatureID, Type_WBOut = new_WBARECOMI) %>%
    filter()
  
  # Attribute those NID facilities precisely at the outlet
  if(nrow(res_ops_table) > 0){
    tmp_POIs <- tmp_POIs %>%
      left_join(res_ops_table, by = "Type_WBOut") %>%
      mutate(Type_NID = ifelse(!is.na(NID_ID), NID_ID, NA)) %>%
      select(-NID_ID)
    
    NID_COMIDs <- filter(NID_COMIDs, !Type_NID %in% res_ops_table$NID_ID)
  }

  # Derive other NID POIs
  tmp_POIs <- POI_creation(NID_COMIDs, filter(nhd, poi == 1), "NID") %>%
    addType(., tmp_POIs, "NID", nexus = TRUE, 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") 
```

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(MINELEVSMO) - 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 * ft_per_km)/ VA_MA / s_per_hr ) %>%
    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, nexus == TRUE | !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)))
    
    # if a POI will be a nexus, it can not be terminal.
    final_POIs <- mutate(final_POIs, Type_Term = ifelse(nexus == 1, NA, Type_Term))
    
    # 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
    
    tmp_POIs$nexus <- as.integer(tmp_POIs$nexus)

    # if a POI will be a nexus, it can not be terminal.
    tmp_POIs <- mutate(tmp_POIs, Type_Term = ifelse(nexus == 1, NA, Type_Term))
        
    # 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(nhd_Final, nav_gpkg, nhd_flowline)
  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"))
noDA_pois <- filter(final_POIs, COMID %in% filter(sub, AreaSqKM == 0)$COMID)
write_sf(noDA_pois, temp_gpkg, "noDA_pois")
```

```{r POI Collapse}
# number POIs
final_POIs <- mutate(final_POIs, id = row_number(), moved = NA) %>%
  write_sf(temp_gpkg, pois_all_layer)
collapse <- TRUE

#  Load data
if(collapse){
  
  # Move HUC12 to other POIs
  moved_pois <- poi_move(final_POIs, "Type_HUC12", .05, 2.5, c("Type_Gages", "Type_TE", "Type_NID")) 
  if(!is.data.frame(moved_pois)){
    final_POIs <- moved_pois$final_pois
    moved_pois_table <- moved_pois$moved_points %>%
      mutate(move_type = "huc12 to other")
  } else {
    final_POIs <- moved_POIs
  }
  
  # Gages to confluences
  moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_Conf") # 47
  if(!is.data.frame(moved_pois)){
    final_POIs <- moved_pois$final_pois
    moved_pois_table <- moved_pois_table %>%
      rbind(moved_pois$moved_points %>%
              mutate(move_type = "gages to conf"))
  } else {
    final_POIs <- moved_POIs
  }
  
  
  # Gages to waterbody inlets
  moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_WBIn") # 2
  if(!is.data.frame(moved_pois)){
    final_POIs <- moved_pois$final_pois
    moved_pois_table <- moved_pois_table %>%
      rbind(moved_pois$moved_points %>%
              mutate(move_type = "gages to wbin"))
  } else {
    final_POIs <- moved_pois
  }
  
  # Gages to waterbody outlets
  moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_WBOut") # 6
  if(!is.data.frame(moved_pois)){
    final_POIs <- moved_pois$final_pois
    moved_pois_table <- moved_pois_table %>%
      rbind(moved_pois$moved_points %>%
              mutate(move_type = "gages to wbout"))
  } else {
    final_POIs <- moved_pois
  }
  
  
  # Streamgage to term outlet
  moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_Term") 
  if(!is.data.frame(moved_pois)){
    final_POIs <- moved_pois$final_pois
    moved_pois_table <- moved_pois_table %>%
      rbind(moved_pois$moved_points %>%
              mutate(move_type = "gages to term"))
  } else {
    final_POIs <- moved_pois
  }
  
  # NID to waterbody outlet
  moved_pois <- poi_move(final_POIs, "Type_NID", .025, 1, "Type_WBOut") 
  if(!is.data.frame(moved_pois)){
    final_POIs <- moved_pois$final_pois
    moved_pois_table <- moved_pois_table %>%
      rbind(moved_pois$moved_points %>%
              mutate(move_type = "nid to wb_out"))
  } else {
    final_POIs <- moved_pois
  }

  write_sf(final_POIs, nav_gpkg, pois_all_layer)
  write_sf(moved_pois_table, temp_gpkg, "pois_collapsed")
} 

check_dups <- final_POIs %>%
  group_by(COMID) %>%
  filter(n() > 1) 

if(nrow(filter(check_dups, all(c(0,1) %in% nexus))) != nrow(check_dups)){
  print("Multiple POI ids at same geometric location")
  no_dups <- filter(check_dups, all(c(0,1) %in% nexus))
  dups <- filter(check_dups, !id %in% no_dups$id)
  write_sf(dups, temp_gpkg, dup_pois)
} else {
  print("All double COMIDs nexus for gage or WB splitting")
}

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

```{r lookup}
# Final POI layer
POIs <- final_POIs %>%
  arrange(COMID) %>%
  mutate(identifier = row_number())

# Unique POI geometry
final_POI_geom <- POIs %>%
  select(identifier) %>%
  cbind(st_coordinates(.)) %>%
  group_by(X, Y) %>%
  mutate(geom_id = cur_group_id()) %>%
  ungroup()

final_POIs_table <- POIs %>%
  inner_join(select(st_drop_geometry(final_POI_geom), -X, -Y), by = "identifier")  %>%
  select(-identifier) 

# POI data theme table
pois_data_orig <- reshape2::melt(st_drop_geometry(select(final_POIs_table,
                                             -nexus)),
                            id.vars = c("COMID", "geom_id")) %>%
  filter(!is.na(value)) %>%
  group_by(COMID, geom_id) %>%
  mutate(identifier = cur_group_id()) %>%
  rename(hy_id = COMID, poi_id = identifier, hl_reference = variable, hl_link = value) %>%
  distinct() 

pois_data_moved <- select(st_drop_geometry(moved_pois_table), hy_id = COMID, hl_link = new_val, hl_reference = moved_value) %>%
  inner_join(distinct(pois_data_orig, hy_id, geom_id, poi_id), by = "hy_id") 

pois_data <- data.table::rbindlist(list(pois_data_moved, pois_data_orig), use.names = TRUE) %>%
  filter(!hl_reference %in% c("id", "moved"))


# POI Geometry table
poi_geometry <- select(final_POIs_table, hy_id = COMID, geom_id) %>%
  inner_join(distinct(pois_data, hy_id, geom_id, poi_id),
             by = c("hy_id" = "hy_id", "geom_id" = "geom_id")) %>%
  distinct() %>%
  st_as_sf()

write_sf(pois_data, nav_gpkg, poi_data_table)
write_sf(poi_geometry, nav_gpkg, poi_geometry_table)

poi_geom_xy <- cbind(poi_geometry, st_coordinates(poi_geometry)) %>%
 st_drop_geometry()

events_data <- events %>%
  arrange(COMID) %>%
  cbind(st_coordinates(.)) %>%
  st_drop_geometry() %>%
  group_by(COMID, REACHCODE, REACH_meas) %>%
  mutate(event_id = cur_group_id()) %>%
  rename(hy_id = COMID) %>%
  ungroup()

nexi <- filter(final_POIs_table, nexus == 1) %>%
  cbind(st_coordinates(.)) %>%
  select(hy_id = COMID, X, Y) %>%
  inner_join(poi_geom_xy, by = c("hy_id" = "hy_id", "X" = "X", "Y" = "Y")) %>%
  inner_join(events_data, by = c("hy_id" = "hy_id", "X" = "X", "Y" = "Y"), multiple = "all") %>%
  select(hy_id, REACHCODE, REACH_meas, event_id, poi_id) %>%
  group_by(hy_id, REACHCODE) %>%
  filter(REACH_meas == min(REACH_meas)) %>%
  ungroup()
  #distinct(hy_id, REACHCODE, REACH_meas, event_id, poi_id)

write_sf(select(events_data, -c(nexus, X, Y)), nav_gpkg, event_table)
write_sf(nexi, nav_gpkg, event_geometry_table)

#  Load data
if(needs_layer(nav_gpkg, lookup_table_refactor)) {
  full_cats <- readRDS(data_paths$fullcats_table)

  lookup <- dplyr::select(sf::st_drop_geometry(nhd_Final),
                          NHDPlusV2_COMID = COMID,
                          realized_catchmentID = COMID,
                          mainstem = LevelPathI) %>%
    dplyr::mutate(realized_catchmentID = ifelse(realized_catchmentID %in% full_cats$FEATUREID,
                                                realized_catchmentID, NA)) %>%
    left_join(select(st_drop_geometry(poi_geometry), hy_id, poi_geom_id = geom_id), 
              by = c("NHDPlusV2_COMID" = "hy_id"))

  sf::write_sf(lookup, nav_gpkg, lookup_table_refactor)
}
```