From 2a830f5d87fe2d1ce8c54685f5662624485feeb4 Mon Sep 17 00:00:00 2001
From: Bock <abock@usgs.gov>
Date: Mon, 8 May 2023 10:00:00 -0600
Subject: [PATCH] Major mods see MR for details

---
 workspace/02_NHD_navigate.Rmd | 411 ++++++++++++++++++++++------------
 1 file changed, 265 insertions(+), 146 deletions(-)

diff --git a/workspace/02_NHD_navigate.Rmd b/workspace/02_NHD_navigate.Rmd
index a03d502..4ad3e90 100644
--- a/workspace/02_NHD_navigate.Rmd
+++ b/workspace/02_NHD_navigate.Rmd
@@ -22,8 +22,6 @@ knitr::opts_chunk$set(
 source("R/utils.R")
 source("R/NHD_navigate.R")
 source("R/hyrefactor_funs.R")
-#source("R/wb_poi_collapse.R")
-#source("R/wb_inlet_collapse.R")
 
 # increase timeout for data downloads
 options(timeout=600)
@@ -35,15 +33,15 @@ source("R/config.R")
 gages <- read_sf(gage_info_gpkg, "Gages")
 
 # NWM network
-nc <- RNetCDF::open.nc(data_paths$nwm_network)
+#nc <- RNetCDF::open.nc(data_paths$nwm_network)
 # NWM gages
-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))
+#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)
+#RNetCDF::close.nc(nc)
 
 # need some extra attributes for a few POI analyses
 atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.rds"))
@@ -54,7 +52,7 @@ atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.r
 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)
+   #try(nhd <- select(nhd, -c(minNet, WB, struct_POI, struct_Net, POI_ID, dend)), silent = TRUE)
   
   # Some NHDPlus VPUs include HUCs from other VPUs
   if(vpu == "02"){
@@ -90,13 +88,12 @@ mapview(filter(tmp_POIs, Type_HUC12 != ""), layer.name = "HUC12 POIs", col.regio
 ```
 
 ```{r streamgage POIs}
-if(all(is.na(tmp_POIs$Type_Gages))) { 
+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) %>%
-    #st_drop_geometry() %>%
     switchDiv(., nhd) 
   
   streamgages <- streamgages_VPU %>% 
@@ -106,8 +103,9 @@ if(all(is.na(tmp_POIs$Type_Gages))) {
     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, nhd, combine_meters, reach_meas_thresh)
+  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) {
@@ -122,7 +120,7 @@ if(all(is.na(tmp_POIs$Type_Gages))) {
   }
   
   # Write out events and outelts
-  write_sf(rename(gages_POIs$events, POI_identifier = Type_Gages), temp_gpkg, split_layer)
+  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 {
@@ -134,7 +132,7 @@ mapview(filter(tmp_POIs, !is.na(Type_Gages)), layer.name = "Streamgage POIs", co
 ```
 
 ```{r TE POIs}
-if(all(is.na(tmp_POIs$Type_TE))) {
+if(!"Type_TE" %in% names(tmp_POIs)) {
   
   if(vpu == "08"){
     nhd$VPUID <- "08"
@@ -153,7 +151,7 @@ if(all(is.na(tmp_POIs$Type_TE))) {
     ungroup()
   
    # Derive TE POIs
-  tmp_POIs <- POI_creation(st_drop_geometry(TE_COMIDs), filter(nhd, dend == 1), "TE") %>%
+  tmp_POIs <- POI_creation(st_drop_geometry(TE_COMIDs), filter(nhd, poi == 1), "TE") %>%
     addType(., tmp_POIs, "TE") 
 
   # As a fail-safe, write out list of TE plants not assigned a POI
@@ -172,38 +170,9 @@ if(all(is.na(tmp_POIs$Type_TE))) {
 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 >= min_da_km) %>%
-    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 waterbody outlet POIs}
 #  Derive or load Waterbody POIs ----------------------
-if(all(is.na(tmp_POIs$Type_WBOut))) {
-  
-  events <- read_sf(temp_gpkg, split_layer)
+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"), 
@@ -211,6 +180,10 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
     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) %>%
@@ -239,14 +212,15 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
   # 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(nhd, WBs_VPU, data_paths, crs)
+  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
-  #tmp_POIs <- wb_layers$POIs
   
-  wb_out_col <- wb_poi_collapse(wb_layers$POIs, nhd, events)
-  ho <- filter(wb_layers$POIs, !COMID %in% wb_out_col$POIs$COMID)
-  write_sf(wb_out_col$events_ret, temp_gpkg, split_layer)
+  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)) %>%
@@ -263,8 +237,22 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
     }
   }
   
+  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(WBs_VPU, nav_gpkg, WBs_layer)
+  write_sf(ref_WB, temp_gpkg, WBs_layer)
 
   # Write specific ResOpsUS data to a separate sheet
   if(nrow(wb_lst) > 0){
@@ -281,16 +269,43 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
 } else {
   tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
   nhd <- read_sf(nav_gpkg, nhd_flowline)
-  WBs_VPU <- read_sf(nav_gpkg, WBs_layer)
+  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(all(is.na(tmp_POIs$Type_Conf))) {
+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))) 
@@ -325,16 +340,53 @@ if(all(is.na(tmp_POIs$Type_Conf))) {
 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)
+  #ho <- filter(wb_layers$POIs, !COMID %in% wb_in_col$POIs$COMID)
+  #write_sf(wb_in_col$events_ret, temp_gpkg, split_layer)
+  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(all(is.na(tmp_POIs$Type_NID))) {
+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(FlowLcomid) %>%
+    group_by(COMID) %>%
     summarize(Type_NID = paste0(unique(NIDID), collapse = " ")) 
   
   # Determine which NID facilitites have been co-located with ResOps
@@ -349,14 +401,14 @@ if(all(is.na(tmp_POIs$Type_NID))) {
   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, Type_NID)) %>%
+      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, dend ==1), "NID") %>%
+  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
@@ -369,43 +421,6 @@ if(all(is.na(tmp_POIs$Type_NID))) {
 mapview(filter(tmp_POIs, Type_NID != ""), layer.name = "NID POIs", col.regions = "blue") 
 ```
 
-```{r waterbody inlet POIs}
-#  Derive or load Waterbody POIs ----------------------
-if(all(is.na(tmp_POIs$Type_WBIn))) {
-
-  wb_layers <- wbin_POIcreation(nhd, WBs_VPU, data_paths, crs, split_layer, tmp_POIs)
-  #tmp_POIs <- wb_layers$POIs 
-  wb_in_col <- wb_inlet_collapse(wb_layers$POIs, nhd, events)
-  #ho <- filter(wb_layers$POIs, !COMID %in% wb_in_col$POIs$COMID)
-  #write_sf(wb_in_col$events_ret, temp_gpkg, split_layer)
-  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") 
-```
-
-
 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
@@ -424,7 +439,7 @@ elev_pois_split <- inc_segs %>%
   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)) %>%
+  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
@@ -506,18 +521,18 @@ if(all(is.na(tmp_POIs$Type_Travel))) {
     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") 
+  # 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()
     
-    write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+  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)
@@ -611,48 +626,89 @@ if(needs_layer(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}
-#  Load data
-if(needs_layer(nav_gpkg, final_poi_layer)) {
+# number POIs
+final_POIs <- mutate(final_POIs, id = row_number(), moved = NA)
+collapse <- TRUE
 
-  if(!"Type_Elev" %in% names(final_POIs)) final_POIs$Type_Elev <- rep(NA_real_, nrow(final_POIs))
-  # Convert nexus to 0,1
-  final_POIs <- mutate(final_POIs, nexus = ifelse(is.na(nexus), 0, nexus))
+#  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
+  }
   
-  #1 Move POIs downstream by category, nexus POIs get removed
-  out_HUC12 <- POI_move_down(temp_gpkg, final_POIs, nsegments, filter(nhd_Final, !is.na(POI_ID)), "Type_HUC12",
-                             poi_dar_move *2)
-  out_HUC12$allPOIs$nexus <- as.numeric(out_HUC12$allPOIs$nexus)
-  out_gages <- POI_move_down(temp_gpkg, out_HUC12$allPOIs, out_HUC12$segs, out_HUC12$FL, 
-                             "Type_Gages", poi_dar_move) 
   
-  # Convert empty strings to NA for handling within R
-  out_gages$allPOIs <- mutate_all(out_gages$allPOIs, list(~na_if(.,""))) 
+  # 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
+  }
   
-  out_gages$allPOIs$snapped <- out_gages$allPOIs$COMID %in% new_POIs$COMID
-
-  out_gages$allPOIs$identifier <- seq(1, nrow(out_gages$allPOIs))
+  # 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
+  }
   
-  out_gages$allPOIs <- select(out_gages$allPOIs, identifier, everything())
   
-  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)
-  write_sf(nhd_Final, temp_gpkg, nhd_flowline)
+  # 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
+  }
   
-  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)
-}
+  # 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 wnpit"))
+  } else {
+    final_POIs <- moved_pois
+  }
+
+  write_sf(moved_pois_table, temp_gpkg, "pois_collapsed")
+} 
+
+write_sf(final_POIs, nav_gpkg, pois_all_layer)
 
 check_dups <- final_POIs %>%
   group_by(COMID) %>%
@@ -661,7 +717,7 @@ check_dups <- final_POIs %>%
 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, !identifier %in% no_dups$identifier)
+  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")
@@ -672,19 +728,82 @@ mapview(final_POIs, layer.name = "All POIs", col.regions = "red") +
 ```
 
 ```{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 <- 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()
+
+# 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, 
+
+  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(nexus_from_poi(dplyr::rename(final_POIs, ID = COMID)), by = c("NHDPlusV2_COMID" = "ID"))
-  
+    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)
 }
-
 ```
-- 
GitLab