diff --git a/.gitignore b/.gitignore
index a7f31b22ec56a5cff613aa38aec0e56a93e79e53..719b9046a64288fc47cf19f3b9da9630a39c4590 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,6 @@
 .DS_Store
 data/
+hyfabric/
 NID.R
 *.ipynb_checkpoints
 workspace/data
diff --git a/hyfabric/DESCRIPTION b/hyfabric/DESCRIPTION
index 98ead0626063a8ba032a69e37dff8c03eec20851..24b87d051d84d0f753fb9b75bddc276f8c517661 100644
--- a/hyfabric/DESCRIPTION
+++ b/hyfabric/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: hyfabric
 Type: Package
 Title: Utility functions for creating the reference geospatial fabric.
-Version: 0.5.6
+Version: 0.5.7
 Authors@R: c(person(given = "David", 
                     family = "Blodgett", 
                     role = c("aut", "cre"), 
diff --git a/hyfabric/R/poi_creation.R b/hyfabric/R/poi_creation.R
index 63f070c5f18d1fb40d9349972984d2248555a2a1..2cc96a0a5cc42864c0ed0a01e021e7d08b9d6a10 100644
--- a/hyfabric/R/poi_creation.R
+++ b/hyfabric/R/poi_creation.R
@@ -53,7 +53,7 @@ POI_creation<-function(srcData, nhdDF, IDfield){
 addType <- function(new_POIs, POIs, IDfield, nexus = TRUE, bind = TRUE){
 
   new_POIs <- st_compatibalize(new_POIs, POIs)
-  
+
   # subset Nexus POIs from Type
   if(nexus){
     nexus_POIs <- filter(POIs, nexus == TRUE)
@@ -80,7 +80,7 @@ addType <- function(new_POIs, POIs, IDfield, nexus = TRUE, bind = TRUE){
 
   POIs_fin <- POIs_exist %>%
     left_join(st_drop_geometry(POIs_newAtt) %>%
-                select(COMID, ID2), by = "COMID", all.x = TRUE) %>%
+                select(COMID, ID2), by = "COMID") %>% #, all.x = TRUE) %>%
     mutate(ID = ifelse(!is.na(ID2), ID2, ID)) %>%
     rename(!!(paste0("Type_", IDfield)) := ID) %>%
     select(-ID2)
@@ -89,7 +89,7 @@ addType <- function(new_POIs, POIs, IDfield, nexus = TRUE, bind = TRUE){
   if(bind){
     POIs_fin <- rbind(POIs_fin, filter(new_POIs, !COMID %in% POIs_fin$COMID))
   }
-  
+
   # Add nexus back in if excluded
   if(nexus){
     POIs_fin <- rbind(POIs_fin, nexus_POIs)
diff --git a/workspace/02_NHD_navigate.Rmd b/workspace/02_NHD_navigate.Rmd
index 985cc60b9fb8a7bb4993a2c9096a9a4c90fd094a..a03d5026dc43d68bafa2292984e807822be40aee 100644
--- a/workspace/02_NHD_navigate.Rmd
+++ b/workspace/02_NHD_navigate.Rmd
@@ -22,6 +22,8 @@ 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)
@@ -29,10 +31,12 @@ options(timeout=600)
 # Load Configuration of environment
 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")) %>%
@@ -41,9 +45,8 @@ nwm_gages <- data.frame(
 
 RNetCDF::close.nc(nc)
 
-# need some extra attributes
+# need some extra attributes for a few POI analyses
 atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.rds"))
-
 ```
 
 ```{r huc12 POIs}
@@ -51,9 +54,9 @@ 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)), silent = TRUE)
+   try(nhd <- select(nhd, -c(minNet, WB, struct_POI, struct_Net, POI_ID, dend, poi)), silent = TRUE)
   
-  # HUC02 includes some 
+  # Some NHDPlus VPUs include HUCs from other VPUs
   if(vpu == "02"){
     grep_exp <-"^02|^04"
   } else if (vpu == "08") {
@@ -124,6 +127,7 @@ if(all(is.na(tmp_POIs$Type_Gages))) {
   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") 
@@ -199,6 +203,8 @@ mapview(filter(tmp_POIs, !is.na(Type_Term)), layer.name = "Terminal POIs", col.r
 #  Derive or load Waterbody POIs ----------------------
 if(all(is.na(tmp_POIs$Type_WBOut))) {
   
+  events <- read_sf(temp_gpkg, split_layer)
+  
   # Waterbodies sourced from NHD waterbody layer for entire VPU
   WBs_VPU_all <- filter(readRDS("data/NHDPlusNationalData/nhdplus_waterbodies.rds"), 
                         COMID %in% nhd$WBAREACOMI) %>%
@@ -235,7 +241,12 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
   
   wb_layers <- wbout_POI_creaton(nhd, WBs_VPU, data_paths, crs)
   WBs_VPU <- wb_layers$WBs
-  tmp_POIs <- wb_layers$POIs
+  #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)
+  tmp_POIs <- wb_out_col$POIs
   
   if(!all(is.na(wb_layers$events))){
     wb_events <- select(wb_layers$events, -c(id, offset)) %>%
@@ -244,7 +255,8 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
     # Write out events and outlets
     if(!needs_layer(temp_gpkg, split_layer)){
       events <- read_sf(temp_gpkg, split_layer) %>%
-        rbind(st_compatibalize(wb_events,.))
+        rbind(st_compatibalize(wb_events,.)) %>%
+        unique()
       write_sf(events, temp_gpkg, split_layer)
     } else {
       write_sf(wb_events, nav_gpkg, split_layer)
@@ -362,7 +374,11 @@ mapview(filter(tmp_POIs, Type_NID != ""), layer.name = "NID POIs", col.regions =
 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 
+  #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)) %>%
@@ -397,7 +413,8 @@ 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,
+  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
@@ -609,7 +626,7 @@ if(needs_layer(nav_gpkg, final_poi_layer)) {
                              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)
+                             "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(.,""))) 
diff --git a/workspace/R/NHD_navigate.R b/workspace/R/NHD_navigate.R
index 9deee9b020b0a075d48be9cf1241b82595d715cf..dfeaba0823f9f1b3395cc7e41ea602783e8b7025 100644
--- a/workspace/R/NHD_navigate.R
+++ b/workspace/R/NHD_navigate.R
@@ -15,162 +15,6 @@ NetworkNav <- function(inCom, nhdDF, withTrib){
   return(seg)
 }
 
-#' ### DEPRECATED ###
-#' #' Identifies and connects dangles in network generated by Network Nav function
-#' #' @param inCOM  (list) list of input COMIDs 
-#' #' @param nhdDF (sf data.frame) (data frame) valid data frame of NHD flowlines
-#' #' @param withTrib (logical) flag for if the upstream navigation should include tributaries
-#' #              or stick to mainstem level path
-#' #
-#' #' @return (list) list of COMIDs connecting dangle to existing network
-#' NetworkConnection <- function(incom, nhd){
-#' 
-#'   upnet_DF <- filter(nhd, COMID %in% incom) %>%
-#'     filter(!DnHydroseq %in% Hydroseq)
-#'   
-#'   # while the number of dangles is greater than 0
-#'   while (length(upnet_DF$COMID) > 0){
-#'     # create item for number of dangles
-#'     count <- dim(upnet_DF)[1]
-#'     print (dim(upnet_DF))
-#'     # find out which level paths are downstream of dangling huc12 POIs
-#'     DSLP <- upnet_DF %>% pull(DnLevelPat)#[upnet_DF$COMID %in% incom]
-#'     # Get the COMID of the hydroseq with level path value
-#'     # the lowest downstream flowline within the levelpath
-#'     inCom2 <- nhd$COMID[nhd$Hydroseq %in% DSLP]
-#'     # Run the upstream navigation code
-#'     upNet <- unique(unlist(lapply(inCom2, NetworkNav, nhd)))
-#'     # Append result to existing segment list
-#'     incom <- append(incom, upNet)
-#'     # Get the same variable as above
-#'     upnet_DF <- filter(nhd, COMID %in% incom, !DnHydroseq %in% Hydroseq)
-#'     # Get the count
-#'     count2 <- dim(upnet_DF)[1]
-#'     # if the count has remained the same we are done and return the flowline list
-#'     if (count == count2){
-#'       return (incom)
-#'     }
-#'   }
-#'   # Not sure this other return is needed
-#'   return(incom)
-#' }
-
-## Deprecated -- see hyfabric package
-#' #' Switches valid POIs from minor to major path divergences
-#' #' @param inSegs  (list) list of input COMIDs representing POIs
-#' #' @param nhdDF (sf data.frame) (data frame) valid data frame of NHD flowlines
-#' #' 
-#' #' @return (sf data.frame) Corresponding major path COMID for POI
-#' switchDiv <- function(insegs, nhdDF){
-#' 
-#'   div_segs <- filter(nhdDF, COMID %in% insegs$COMID)
-#'   if (2 %in% div_segs$Divergence){
-#'     print ("Switching divergence to other fork")
-#'     
-#'     # Look Upstream
-#'     upstream <- st_drop_geometry(nhdDF) %>% 
-#'       inner_join(st_drop_geometry(div_segs) %>% 
-#'                    filter(Divergence == 2) %>% 
-#'                    select(COMID, Hydroseq), 
-#'                  by = c("DnMinorHyd" = "Hydroseq"))
-#'     
-#'     # From Upstream Segment switch POI to the downstream major/main path
-#'     insegs_maj <- st_drop_geometry(nhdDF) %>% 
-#'       inner_join(upstream %>% select(COMID.y, DnHydroseq), 
-#'                  by = c("Hydroseq" = "DnHydroseq")) %>% 
-#'       select(COMID, COMID.y)
-#'     
-#'     insegs <- insegs %>% 
-#'       left_join(insegs_maj, by = c("COMID" = "COMID.y")) %>%
-#'       mutate(COMID = ifelse(!is.na(COMID.y), COMID.y, COMID)) %>% select(-COMID.y)
-#' 
-#'   } else {
-#'     print ("no divergences present in POI set")
-#'   }
-#'   return(insegs)
-#' }
-
-
-## Deprecated -- see hyfabric package
-#' #' Creates POIs for a given data theme
-#' #' @param srcData  (data.frame) (data frame) DF with two columns: 
-#' #             1 - COMID
-#' #             2 - Unique ID value for POI (Streamgage ID, etc.)
-#' #' @param nhdDF (sf data.frame) valid data frame of NHD flowlines
-#' #' @param IDfield character) Name of 'Type' field to be modified in POI 
-#' #' 
-#' #' @return (sf data.frame) OIs for the specific data theme
-#' POI_creation<-function(srcData, nhdDF, IDfield){
-#' 
-#'   # Give generic ID to Identity field
-#'   colnames(srcData) <- c("COMID", "ID")
-#'   
-#'   sub_segs <- filter(nhdDF, COMID %in% srcData$COMID) 
-#' 
-#'   POIs <- sub_segs %>% 
-#'     get_node(., position = "end") %>% 
-#'     mutate(COMID = sub_segs$COMID) %>%
-#'     mutate(Type_HUC12 = NA, Type_WBIn = NA, Type_WBOut = NA, Type_Gages = NA, Type_TE = NA, Type_NID = NA, Type_Conf = NA) %>%
-#'     inner_join(srcData %>% select(COMID, ID), by = "COMID") %>% 
-#'     mutate(!!(paste0("Type_", IDfield)) := ID)
-#'   
-#'   if(!(paste0("Type_", IDfield)) %in% colnames(POIs)){
-#'     POIs <- POIs %>% select(COMID, Type_HUC12, Type_Gages, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf)
-#'   } else {
-#'     POIs <- POIs %>% select(COMID, Type_HUC12, Type_Gages, Type_TE, 
-#'                             Type_NID, Type_WBIn, Type_WBOut, Type_Conf, !!(paste0("Type_", IDfield)))
-#'   }
-#' 
-#'   return(POIs)
-#' }
-
-# Deprecated -- see hyfabric package.
-#' #' Adds the type attribute for co-located POIs of multiple themes
-#' #' @param new_POIs (sf data.frame) new POIs to be tested against existing
-#' #' @param POIs  (sf data.frame) Existing POIs
-#' #' @param IDfield (character) Name of 'Type' field to be modified in existing POI 
-#' #' @param bind (logical) whether to bind non co-located POIs to data frame or just
-#' #'             attribute existing POIs 
-#' #' 
-#' #' @return (sf data.frame) xisting POIs with Type field modified for
-#' #            duplicate POIs for given data theme
-#' addType<-function(new_POIs, POIs, IDfield, bind){
-#' 
-#'   new_POIs <- st_compatibalize(new_POIs, POIs) 
-#'   
-#'   if(paste0("Type_", IDfield) %in% colnames(POIs)){
-#'     POIs_exist <- POIs %>% 
-#'       rename(ID = !!(paste0("Type_", IDfield)))
-#'   } else {
-#'     POIs_exist <- POIs %>%
-#'       mutate(ID = NA)
-#'   }
-#' 
-#'   # Add columns missing between master POI and new set
-#'   missing_cols <- colnames(POIs)[!colnames(POIs) %in% colnames(new_POIs)]
-#'   for(col in missing_cols){
-#'     new_POIs <-  new_POIs %>%
-#'       mutate(!!col := NA)
-#'   }
-#'   
-#'   POIs_newAtt <- filter(new_POIs, COMID %in% POIs$COMID) %>%
-#'     rename(ID2 = !!(paste0("Type_", IDfield)))
-#' 
-#'   POIs_fin <- POIs_exist %>% 
-#'     left_join(st_drop_geometry(POIs_newAtt) %>% 
-#'                                  select(COMID, ID2), by = "COMID", all.x = TRUE) %>%
-#'     mutate(ID = ifelse(!is.na(ID2), ID2, ID)) %>% 
-#'     rename(!!(paste0("Type_", IDfield)) := ID) %>% 
-#'     select(-ID2) 
-#'  
-#'   # Bind unless indicated not to
-#'   if(missing(bind)){
-#'     POIs_fin <- rbind(POIs_fin, filter(new_POIs, !COMID %in% POIs_fin$COMID))
-#'   }
-#'   
-#'   return(POIs_fin)
-#' }
-
 #' Creates raw and dissolved segment layers with necessaary
 #         upstream/downstream routing attribution
 #'  @param nhdDF (sf data.frame) valid data frame of NHD flowlines
@@ -325,10 +169,12 @@ DS_poiFix <- function(POIs_wgeom, nhd){
     filter(n() > 1) %>%
     # Spurs warnings where there are NAs for a column 
     #       for a given grouped COMID, but output is what I expect
-    summarise_all(funs(min(unique(.[!is.na(.)]), na.rm = T))) %>%
+    summarise_all(funs(min(unique(.[!is.na(.)]), na.rm = T))) #%>%
     # Replace -Inf with NA where applicaable
-    mutate_all(~na_if(., -Inf))
+    #mutate_all(~na_if(., -Inf))
 
+  poi_fix_unique[mapply(is.infinite, poi_fix_unique)] <- NA  
+    
   # Join new POI COMIDs and geometry with the old Type fields
   new_POIs <- bind_rows(filter(poi_fix, !COMID %in% poi_fix_unique$COMID), poi_fix_unique) %>%
     arrange(COMID) %>%
@@ -505,7 +351,9 @@ POI_move_down<-function(out_gpkg, POIs, Seg, Seg2, exp, DA_diff){
   merged_POIs <- stationary_POI %>%
     inner_join(movedownPOI_withatt,
                by = c("COMID" = "Move"), suffix = c(".x", ".y")) %>%
-    mutate_all(., list(~na_if(.,""))) %>%
+  #merged_POIs[mapply(is.infinite, merged_POIs)] <- NA
+  #merged_POIs <- merged_POIs %>%
+    #mutate_all(., list(~na_if(.,""))) %>%
     # Don't overwrite existing gages or HUC12s
     filter(is.na(!!as.symbol(paste0(exp, ".x")))) %>%
     # Bring over relevant attributes
@@ -789,8 +637,8 @@ split_tt <- function(in_POI_ID, split_DF, tt_diff, max_DA){
 
 #'  Retrieves waterbodies from NHDArea layer or NHD Hi-Res for ResOpsUs 
 #'  locations if absent from NHD waterbody
-#'  @param nhd (sf data.frame) flowlines for given VPU
-#'  @param WBs  (sf data.frame) waterbodys for discretization within VPU
+#'  @param nhd (sf data.frame) Fowlines for given VPU
+#'  @param WBs  (sf data.frame) waterbodiess for discretization within VPU
 #'  @param data_paths (list) data paths to data layers
 #' 
 #'  @return (list) wbs - sf data frame for NHDArea and HR waterbodies
@@ -1081,6 +929,16 @@ WB_event <- function(WBs, nhd_wb, type){
   return(wb_events)
 } 
 
+
+#'  Creates POIs for gages using refactor criteria
+#'  @param temp_POIs  (sf data.frame) current data frame of POIs
+#'  @param gages_info (sf data.frame) VPU specific streamgages from 01_gage_selection
+#'  @param nhd (sf data.frame) flowline data.frame 
+#'  @param combine_meters (integer) refactor threshold (m) for if if two adjacent fl should be combined
+#'  @reach_meas_thresh (integer) threshold added or substracted from reach_meas to determine if a gage will split fl
+#' 
+#'  @return (sf data.frame) dataframe of gage POIs
+#'  
 gage_POI_creation <- function(tmp_POIs, gages_info, nhd, combine_meters, reach_meas_thresh){
   # Create streamgage POIs
   gage_POIs <- POI_creation(select(st_drop_geometry(gages_info), COMID, site_no), nhd, "Gages") %>%
@@ -1117,6 +975,14 @@ gage_POI_creation <- function(tmp_POIs, gages_info, nhd, combine_meters, reach_m
   return(list(events = events, tmp_POIs = tmp_POIs))
 }
 
+#'  Creates Waterbody POIs, calls a few other functions
+#'  @param nhd (sf data.frame) flowline data.frame 
+#'  @param WBs_VPU (sf data.frame) 
+#'  @param data_paths (list) data paths to data layers
+#'  @param crs (integer) CRS to use (epsg code) 
+#' 
+#'  @return (sf data.frame) dataframe of waterbody outlet POIs
+#'  
 wbout_POI_creaton <- function(nhd, WBs_VPU, data_paths, crs){
   # Create waterbody outlet POIs for waterbodies that are in NHDv2 waterbody set
   wbout_COMIDs <- filter(nhd, dend == 1 & WB == 1) %>%
@@ -1180,6 +1046,16 @@ wbout_POI_creaton <- function(nhd, WBs_VPU, data_paths, crs){
   
 }
 
+#'  Creates Waterbody POIs, calls a few other functions
+#'  @param nhd (sf data.frame) flowline data.frame 
+#'  @param WBs_VPU (sf data.frame) waterbodies to create inlet pois from
+#'  @param data_paths (list) list of data paths 
+#'  @param crs (integer) CRS to use (epsg code)
+#'  @param split_layer (sf data.frame) events to split flowlines with
+#'  @param tmp_POIs (sf data.frame) rolling POI data frame
+#' 
+#'  @return (sf data.frame) dataframe of WB inlet POIs
+#'  
 wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs, 
                              split_layer, tmp_POIs){
   
@@ -1305,3 +1181,289 @@ wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs,
   }
   
 }
+
+
+#'  collapses close upstream gages with wb inlet events.
+#'  @param temp_POIs (sf data.frame) rolling POI data frame
+#'  @param nhd (sf data.frame) nhd flowlines
+#'  @param events (list) waterbody inlet events
+#' 
+#'  @return (sf data.frame) dataframe of gage POIs
+#' 
+wb_inlet_collapse <- function(tmp_POIs, nhd, events){
+  gage_dist_node <- function(x, wb_ds_ds, gage_add, events){
+    print (x)
+    wb_in_fl <- filter(wb_ds_ds, COMID == x)
+    gage_ds <- filter(wb_ds_ds, Hydroseq %in% wb_in_fl$Hydroseq |
+                        Hydroseq %in% wb_in_fl$DnHydroseq) 
+    
+    gage_reach <- gage_ds %>%
+      group_by(REACHCODE) %>%
+      summarize(do_union = T,
+                total_length = sum(LENGTHKM))
+    
+    #print(nrow(gage_reach))
+    gage_event <- filter(gage_add, COMID %in% gage_ds$COMID)
+    wb_event <- filter(events, COMID == x, event_type == "inlet") %>%
+      unique()
+    
+    if(nrow(gage_reach) == 0){
+      print("Wilton")
+    }
+    
+    if(nrow(gage_event) == 0){
+      return("No gage events")
+    } else if(gage_event$COMID != wb_in_fl$COMID) {
+      gage_reach <- gage_reach %>%
+        filter(REACHCODE == gage_event$reachcode) %>%
+        mutate(gage_dist = ifelse(gage_event$nexus == TRUE,
+                                  total_length * (1 - (gage_event$reach_meas/100)),
+                                  total_length)) %>%
+        mutate(gage_comid = gage_event$COMID,
+               wbin_comid = x)
+    } else if(gage_event$COMID == wb_in_fl$COMID){
+      if(nrow(wb_event) >0){
+        wb_in_meas <- wb_event$REACH_meas
+        wb_RC <- wb_event$REACHCODE
+      } else {
+        wb_out_meas <- wb_in_fl$FromMeas
+        wb_RC <- wb_in_fl$REACHCODE
+      }
+      
+      # gage reach
+      gage_reach <- gage_reach %>%
+        filter(REACHCODE == gage_event$reachcode) %>%
+        mutate(gage_dist = total_length * (1 - (gage_event$reach_meas/100))) 
+      
+      # wb info
+      wb_reach <- gage_reach %>%
+        filter(REACHCODE == wb_RC) %>%
+        mutate(wb_dist = total_length * (1 - (wb_out_meas /100)))
+      
+      gage_reach <- gage_reach %>%
+        mutate(gage_dist = abs(wb_reach$wb_dist - gage_dist),
+               gage_comid = gage_event$COMID,
+               wbin_comid = x)
+    }
+  }
+  
+  #events <- read_sf(temp_gpkg, split_layer) %>%
+  #  rbind(st_compatibalize(wb_,.))
+  
+  # Previously identified streamgages within Gage_Selection.Rmd
+  streamgages_VPU <- gages %>%
+    rename(COMID = comid) %>%
+    filter(COMID %in% nhd$COMID) %>%
+    #st_drop_geometry() %>%
+    switchDiv(., nhd) 
+  
+  # get waterbody outlets
+  wb_in <- filter(tmp_POIs, !is.na(Type_WBIn), is.na(Type_Gages))
+  
+  wb_in_nexus <- filter(wb_in, nexus == TRUE)
+  wb_in_nhd <- filter(nhd, COMID %in% wb_in$COMID) 
+  wb_ds_ds <- wb_in_nhd %>%
+    rbind(filter(nhd, LevelPathI %in% .$LevelPathI & DnHydroseq %in% .$Hydroseq))
+  
+  # get gages
+  gage_wb <- filter(tmp_POIs, !is.na(Type_Gages) & is.na(Type_WBIn)) #%>%
+  #  filter(COMID %in% wb_ds_ds$COMID) 
+  #gage_nhd <- filter(nhd, COMID %in% gage_wb$COMID)
+  #gage_node <- filter(gage_wb, nexus == FALSE)
+  #gage_nexus <- filter(gage_wb, nexus == TRUE)
+  
+  gage_add <- filter(streamgages_VPU, site_no %in% gage_wb$Type_Gages) %>%
+    select(COMID, reachcode, reach_meas, site_no) %>%
+    inner_join(select(st_drop_geometry(gage_wb), site_no = Type_Gages, nexus), 
+               by = "site_no")
+  
+  output <- lapply(wb_in$COMID, 
+                   gage_dist_node, wb_ds_ds, gage_add, events)
+  
+  output_full <- do.call("rbind", output[lengths(output) > 1]) 
+  
+  if(!is.null(output_full)){
+    output_full <- output_full %>%
+      filter(gage_dist < 1) %>%
+      left_join(select(st_drop_geometry(nhd), COMID, TotDASqKM_WB = TotDASqKM),
+                by = c("wbin_comid" = "COMID")) %>%
+      left_join(select(st_drop_geometry(nhd), COMID, TotDASqKM_gage = TotDASqKM),
+                by = c("gage_comid" = "COMID")) %>%
+      mutate(DAR = TotDASqKM_gage / TotDASqKM_WB) %>%
+      filter(gage_dist < 1, DAR > .975) %>%
+      st_drop_geometry()
+    
+    gage_POI <- filter(tmp_POIs, COMID %in% output_full$gage_comid) %>%
+      select(COMID, Type_HUC12_ds = Type_HUC12, Type_Gages_ds = Type_Gages, 
+             Type_TE_ds = Type_TE, Type_Term_ds = Type_Term, nexus) %>%
+      st_drop_geometry() %>%
+      group_by(COMID) %>%
+      summarise(Type_HUC12_ds = last(na.omit(Type_HUC12_ds)), 
+                Type_Gages_ds = last(na.omit(Type_Gages_ds)),
+                Type_Term_ds = last(na.omit(Type_Term_ds)),
+                Type_TE_ds = last(na.omit(Type_TE_ds)),
+                nexus = last(na.omit(nexus)))
+    
+    WB_POI <- filter(tmp_POIs, COMID %in% output_full$wbin_comid, !is.na(Type_WBIn)) %>%
+      left_join(select(output_full, gage_comid, wbin_comid), by = c("COMID" = "wbin_comid")) %>%
+      left_join(select(gage_POI, -nexus), by = c("gage_comid" = "COMID")) %>%
+      mutate(Type_HUC12 = ifelse(!is.na(Type_HUC12_ds), Type_HUC12_ds, Type_HUC12),
+             Type_Gages = ifelse(!is.na(Type_Gages_ds), Type_Gages_ds, Type_Gages),
+             Type_TE = ifelse(!is.na(Type_TE_ds), Type_TE_ds, Type_TE),
+             Type_Term = ifelse(!is.na(Type_Term_ds), Type_Term_ds, Type_Term)) %>%
+      select(-c(Type_HUC12_ds, Type_Gages_ds, Type_TE_ds, Type_Term_ds))
+    
+    tmp_POIs_fin <- filter(tmp_POIs, !COMID %in% c(WB_POI$COMID, WB_POI$gage_comid)) %>%
+      rbind(select(WB_POI, -gage_comid))
+    
+    if(any(gage_POI$nexus == TRUE)){
+      gage_events <- filter(gage_POI, nexus == TRUE)
+      events <- filter(events, !POI_identifier %in% gage_events$Type_Gages_ds)
+    }
+    return(list(POIs = tmp_POIs_fin, events_ret = events))
+  } else {
+    print ("no points collapse")
+    return(list(POIs = tmp_POIs, events_ret = NA))
+  }
+}
+
+#'  Creates Waterbody POIs, calls a few other functions
+#'  @param (sf data.frame) rolling POI data frame
+#'  @param nhd (sf data.frame) nhd flowlines 
+#'  @param events (sf data.frame) waterbody inlet events
+#' 
+#'  @return (sf data.frame) dataframe of wb inlet POIs collapse
+#' 
+#' wb_poi_collapse <- function(tmp_POIs, nhd, events){
+wb_poi_collapse <- function(tmp_POIs, nhd, events){
+  gage_dist_node <- function(x, wb_ds_ds, gage_add, events){
+    print (x)
+    wb_out_fl <- filter(wb_ds_ds, COMID == x)
+    gage_ds <- filter(wb_ds_ds, Hydroseq %in% wb_out_fl$Hydroseq |
+                        Hydroseq %in% wb_out_fl$DnHydroseq) 
+    
+    gage_reach <- gage_ds %>%
+      group_by(REACHCODE) %>%
+      summarize(do_union = T,
+                total_length = sum(LENGTHKM))
+    
+    #print(nrow(gage_reach))
+    gage_event <- filter(gage_add, COMID %in% gage_ds$COMID)
+    wb_event <- filter(events, COMID == x, event_type == "outlet") %>%
+      unique()
+    
+    if(nrow(gage_reach) == 0){
+      print("Wilton")
+    }
+    
+    if(nrow(gage_event) == 0){
+      #print("akermayun")
+      return("Akermayun")
+    } else if(gage_event$COMID != wb_out_fl$COMID) {
+      gage_reach <- gage_reach %>%
+        filter(REACHCODE == gage_event$reachcode) %>%
+        mutate(gage_dist = ifelse(gage_event$nexus == TRUE,
+                                  total_length * (1 - (gage_event$reach_meas/100)),
+                                  total_length)) %>%
+        mutate(gage_comid = gage_event$COMID,
+               wbout_comid = x)
+    } else if(gage_event$COMID == wb_out_fl$COMID){
+      if(nrow(wb_event) >0){
+        wb_out_meas <- wb_event$REACH_meas
+        wb_RC <- wb_event$REACHCODE
+      } else {
+        wb_out_meas <- wb_out_fl$FromMeas
+        wb_RC <- wb_out_fl$REACHCODE
+      }
+      
+      # gage reach
+      gage_reach <- gage_reach %>%
+        filter(REACHCODE == gage_event$reachcode) %>%
+        mutate(gage_dist = total_length * (1 - (gage_event$reach_meas/100))) 
+      
+      # wb info
+      wb_reach <- gage_reach %>%
+        filter(REACHCODE == wb_RC) %>%
+        mutate(wb_dist = total_length * (1 - (wb_out_meas /100)))
+      
+      gage_reach <- gage_reach %>%
+        mutate(gage_dist = abs(wb_reach$wb_dist - gage_dist),
+               gage_comid = gage_event$COMID,
+               wbout_comid = x)
+    }
+  }
+  
+  #events <- read_sf(temp_gpkg, split_layer) %>%
+  #  rbind(st_compatibalize(wb_,.))
+  
+  # Previously identified streamgages within Gage_Selection.Rmd
+  streamgages_VPU <- gages %>%
+    rename(COMID = comid) %>%
+    filter(COMID %in% nhd$COMID) %>%
+    #st_drop_geometry() %>%
+    switchDiv(., nhd) 
+  
+  # get waterbody outlets
+  wb_out <- filter(tmp_POIs, !is.na(Type_WBOut), is.na(Type_Gages))
+  
+  wb_out_nexus <- filter(wb_out, nexus == TRUE)
+  wb_out_nhd <- filter(nhd, COMID %in% wb_out$COMID) 
+  wb_ds_ds <- wb_out_nhd %>%
+    rbind(filter(nhd, LevelPathI %in% .$LevelPathI & Hydroseq %in% .$DnHydroseq))
+  
+  # get gages
+  gage_wb <- filter(tmp_POIs, !is.na(Type_Gages) & is.na(Type_WBOut)) #%>%
+  #  filter(COMID %in% wb_ds_ds$COMID) 
+  #gage_nhd <- filter(nhd, COMID %in% gage_wb$COMID)
+  #gage_node <- filter(gage_wb, nexus == FALSE)
+  #gage_nexus <- filter(gage_wb, nexus == TRUE)
+  
+  gage_add <- filter(streamgages_VPU, site_no %in% gage_wb$Type_Gages) %>%
+    select(COMID, reachcode, reach_meas, site_no) %>%
+    inner_join(select(st_drop_geometry(gage_wb), site_no = Type_Gages, nexus), 
+               by = "site_no")
+  
+  output <- lapply(wb_out$COMID, 
+                   gage_dist_node, wb_ds_ds, gage_add, events)
+  
+  output_full <- do.call("rbind", output[lengths(output) > 1]) %>%
+    filter(gage_dist < 1) %>%
+    left_join(select(st_drop_geometry(nhd), COMID, TotDASqKM_WB = TotDASqKM),
+              by = c("wbout_comid" = "COMID")) %>%
+    left_join(select(st_drop_geometry(nhd), COMID, TotDASqKM_gage = TotDASqKM),
+              by = c("gage_comid" = "COMID")) %>%
+    mutate(DAR = TotDASqKM_WB / TotDASqKM_gage) %>%
+    filter(gage_dist < 1, DAR > .975) %>%
+    st_drop_geometry()
+  
+  gage_POI <- filter(tmp_POIs, COMID %in% output_full$gage_comid) %>%
+    select(COMID, Type_HUC12_ds = Type_HUC12, Type_Gages_ds = Type_Gages, 
+           Type_TE_ds = Type_TE, Type_Term_ds = Type_Term, nexus) %>%
+    st_drop_geometry() %>%
+    group_by(COMID) %>%
+    summarise(Type_HUC12_ds = last(na.omit(Type_HUC12_ds)), 
+              Type_Gages_ds = last(na.omit(Type_Gages_ds)),
+              Type_Term_ds = last(na.omit(Type_Term_ds)),
+              Type_TE_ds = last(na.omit(Type_TE_ds)),
+              nexus = last(na.omit(nexus)))
+  
+  WB_POI <- filter(tmp_POIs, COMID %in% output_full$wbout_comid, !is.na(Type_WBOut)) %>%
+    inner_join(select(output_full, gage_comid, wbout_comid), by = c("COMID" = "wbout_comid")) %>%
+    inner_join(select(gage_POI, -nexus), by = c("gage_comid" = "COMID")) %>%
+    mutate(Type_HUC12 = ifelse(!is.na(Type_HUC12_ds), Type_HUC12_ds, Type_HUC12),
+           Type_Gages = ifelse(!is.na(Type_Gages_ds), Type_Gages_ds, Type_Gages),
+           Type_TE = ifelse(!is.na(Type_TE_ds), Type_TE_ds, Type_TE),
+           Type_Term = ifelse(!is.na(Type_Term_ds), Type_Term_ds, Type_Term)) %>%
+    select(-c(Type_HUC12_ds, Type_Gages_ds, Type_TE_ds, Type_Term_ds))
+  
+  tmp_POIs_fin <- filter(tmp_POIs, !COMID %in% c(WB_POI$COMID, WB_POI$gage_comid)) %>%
+    rbind(select(WB_POI, -gage_comid))
+  
+  if(any(gage_POI$nexus == TRUE)){
+    gage_events <- filter(gage_POI, nexus == TRUE)
+    events <- filter(events, !POI_identifier %in% gage_events$Type_Gages_ds)
+  }
+  
+  
+  return(list(POIs = tmp_POIs_fin, events_ret = events))
+}
\ No newline at end of file
diff --git a/workspace/hyfabric_0.5.6.tar.gz b/workspace/hyfabric_0.5.6.tar.gz
deleted file mode 100644
index 2d025292b5b132001c49e665b072f01e97b884d5..0000000000000000000000000000000000000000
Binary files a/workspace/hyfabric_0.5.6.tar.gz and /dev/null differ
diff --git a/workspace/hyfabric_0.5.7.tar.gz b/workspace/hyfabric_0.5.7.tar.gz
new file mode 100644
index 0000000000000000000000000000000000000000..1a5c37d0d10521ed10cfb76e4b42102f1b6fd52e
Binary files /dev/null and b/workspace/hyfabric_0.5.7.tar.gz differ