diff --git a/workspace/02_NHD_navigate.Rmd b/workspace/02_NHD_navigate.Rmd
index 167a1073e64cc708c590638bf13357c39ae2d377..e2d9ec9f61a8a7ff87ed5f3b41e55b74aae62fb4 100644
--- a/workspace/02_NHD_navigate.Rmd
+++ b/workspace/02_NHD_navigate.Rmd
@@ -32,17 +32,6 @@ source("R/config.R")
 # Gages output from Gage_selection
 gages <- read_sf(gage_info_gpkg, "Gages")
 
-# 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))
-
-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"))
 ```
@@ -52,15 +41,15 @@ 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)), silent = TRUE)
-  
+   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 == "02"){
+  if(vpu_codes == "02"){
     grep_exp <-"^02|^04"
-  } else if (vpu == "08") {
+  } else if (vpu_codes == "08") {
     grep_exp <- "^03|^08"
   } else {
-    grep_exp <- paste0("^", substr(vpu, start = 1, stop = 2))
+    grep_exp <- paste0("^", substr(vpu_codes, start = 1, stop = 2))
   }
   
   # Join HUC12 outlets with NHD
@@ -115,7 +104,7 @@ if(!"Type_Gages" %in% names(tmp_POIs)) {
     # 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"))
+              "Gages_info", paste0("VPU ", vpu_codes, "; low gage score"))
     
   }
   
@@ -134,7 +123,7 @@ mapview(filter(tmp_POIs, !is.na(Type_Gages)), layer.name = "Streamgage POIs", co
 ```{r TE POIs}
 if(!"Type_TE" %in% names(tmp_POIs)) {
   
-  if(vpu == "08"){
+  if(vpu_codes == "08"){
     nhd$VPUID <- "08"
   } else {
     nhd$VPUID <- substr(nhd$RPUID, 1, 2)
@@ -144,7 +133,7 @@ if(!"Type_TE" %in% names(tmp_POIs)) {
   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) %>%
+    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()) %>%
@@ -152,7 +141,7 @@ if(!"Type_TE" %in% names(tmp_POIs)) {
   
    # Derive TE POIs
   tmp_POIs <- POI_creation(st_drop_geometry(TE_COMIDs), filter(nhd, poi == 1), "TE") %>%
-    addType(., tmp_POIs, "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) {
@@ -261,7 +250,8 @@ if(!"Type_WBOut" %in% names(tmp_POIs)) {
       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,".csv")))
+    write.csv(resops_POIs_df, file.path("data/reservoir_Data", 
+                                        paste0("resops_POIs_", vpu_codes, ".csv")))
   }
     
   write_sf(nhd, nav_gpkg, nhd_flowline)
@@ -602,6 +592,7 @@ if(needs_layer(temp_gpkg, nsegments_layer)) {
 
   # 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
@@ -630,7 +621,8 @@ write_sf(noDA_pois, temp_gpkg, "noDA_pois")
 
 ```{r POI Collapse}
 # number POIs
-final_POIs <- mutate(final_POIs, id = row_number(), moved = NA)
+final_POIs <- mutate(final_POIs, id = row_number(), moved = NA) %>%
+  write_sf(temp_gpkg, pois_all_layer)
 collapse <- TRUE
 
 #  Load data
@@ -698,16 +690,15 @@ if(collapse){
     final_POIs <- moved_pois$final_pois
     moved_pois_table <- moved_pois_table %>%
       rbind(moved_pois$moved_points %>%
-              mutate(move_type = "nid to wnpit"))
+              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")
 } 
 
-write_sf(final_POIs, nav_gpkg, pois_all_layer)
-
 check_dups <- final_POIs %>%
   group_by(COMID) %>%
   filter(n() > 1) 
@@ -744,14 +735,21 @@ final_POIs_table <- POIs %>%
   select(-identifier) 
 
 # POI data theme table
-pois_data <- reshape2::melt(st_drop_geometry(select(final_POIs_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()
+  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) %>%