diff --git a/workspace/01_NHD_prep.Rmd b/workspace/01_NHD_prep.Rmd
index 1f977dd5338bfe315652b9d3c2934a627d934435..fa257fcb1b5cd8910e06bd54f3da6e6ccefd887b 100644
--- a/workspace/01_NHD_prep.Rmd
+++ b/workspace/01_NHD_prep.Rmd
@@ -44,6 +44,8 @@ fline <- sf::st_cast(fline, "LINESTRING")
 
 catchment <- sf::read_sf(file.path(data_paths$nhdplus_dir, "reference_catchments.gpkg"))
 
+catchment <- sf::st_transform(catchment, 5070)
+
 # we can remove truely degenerate COMIDs 
 # for 0 upstream area and no catchment area
 degen_comid <- fline[fline$TotDASqKM == 0 & 
diff --git a/workspace/02_NHD_navigate.Rmd b/workspace/02_NHD_navigate.Rmd
index 94d74a4bf7035980fe0053fd6890743bd2648f12..93ed09c5d7c9b57674746d07749f2e88bf7a0c3c 100644
--- a/workspace/02_NHD_navigate.Rmd
+++ b/workspace/02_NHD_navigate.Rmd
@@ -41,11 +41,12 @@ RNetCDF::close.nc(nc)
 # need some extra attributes
 atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.rds"))
 # atts <- select(atts, COMID, MAXELEVSMO, MINELEVSMO, VA_MA, TOTMA, WBAreaType)
+
 ```
 
 ```{r huc12 POIs}
 #  Derive or load HUC12 POIs
-if(needs_layer(nav_gpkg, nav_poi_layer)) {
+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)))
@@ -72,11 +73,11 @@ if(needs_layer(nav_gpkg, nav_poi_layer)) {
   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)
+  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(nav_gpkg, nav_poi_layer) 
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer) 
   nhd <- read_sf(nav_gpkg, nhd_flowline)
 }
 
@@ -107,7 +108,7 @@ if(all(is.na(tmp_POIs$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),
-             nav_gpkg, "unassigned_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) %>%
@@ -117,9 +118,9 @@ if(all(is.na(tmp_POIs$Type_Gages))) {
   }
   
   # Write out geopackage layer representing POIs for given theme
-  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
 }
 
 mapview(filter(tmp_POIs, !is.na(Type_Gages)), layer.name = "Streamgage POIs", col.regions = "blue") 
@@ -151,14 +152,14 @@ if(all(is.na(tmp_POIs$Type_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")
+             temp_gpkg, "unassigned_TE")
   }
   
   # Write out geopackage layer representing POIs for given theme
-  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
   # Load TE POIs if they already exist
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
 }
 
 mapview(filter(tmp_POIs, Type_TE != ""), layer.name = "TE Plant POIs", col.regions = "blue") 
@@ -182,9 +183,9 @@ if(all(is.na(tmp_POIs$Type_Term))) {
   tmp_POIs <- POI_creation(terminal_POIs, nhd, "Term") %>%
     addType(., tmp_POIs, "Term")
 
-  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
   nhd <- read_sf(nav_gpkg, nhd_flowline)
 }
 
@@ -218,11 +219,10 @@ if(all(is.na(tmp_POIs$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)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
   nhd <- read_sf(nav_gpkg, nhd_flowline)
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
 }
 
 mapview(filter(tmp_POIs, !is.na(Type_Conf)), layer.name = "Confluence POIs", col.regions = "blue") 
@@ -245,10 +245,10 @@ if(all(is.na(tmp_POIs$Type_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)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
   # Load NID POIs if they already exist
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
 }
 
 mapview(filter(tmp_POIs, Type_NID != ""), layer.name = "NID POIs", col.regions = "blue") 
@@ -292,7 +292,7 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
     dplyr::filter(COMID %in% wb_lst$COMID)
   
   # Write out waterbodies
-  write_sf(st_transform(WBs_VPU, 4269), nav_gpkg, WBs_layer)
+  write_sf(st_transform(WBs_VPU, sf::st_crs(nhd)), nav_gpkg, WBs_layer)
 
   # Segments that are in waterbodies
   nhd <- mutate(nhd, WB = ifelse(WBAREACOMI > 0 & WBAREACOMI %in% WBs_VPU$COMID, 1, 0))
@@ -335,10 +335,9 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
     write.csv(resops_POIs_df, file.path("data/reservoir_Data", paste0("resops_POIs_",vpu,".csv")))
   }
     
-  write_sf(nhd, nav_gpkg, nhd_flowline)
-  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
   nhd <- read_sf(nav_gpkg, nhd_flowline)
 }
 
@@ -404,9 +403,9 @@ if(nrow(elev_pois_split) > 0) {
   
   tmp_POIs$Type_Elev <- rep(NA, nrow(tmp_POIs))
 
-  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
     
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer) 
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer) 
   nhd <- read_sf(nav_gpkg, nhd_flowline)
 }
 
@@ -457,9 +456,9 @@ if(all(is.na(tmp_POIs$Type_Travel))) {
   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)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 }else{
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
   nhd <- read_sf(nav_gpkg, nhd_flowline)
 }
 
@@ -468,7 +467,7 @@ mapview(filter(tmp_POIs, Type_Travel == 1), layer.name = "Travel time break POIs
   
 ```{r Final POIs}
 # Derive final POI set ----------------------
-if(needs_layer(nav_gpkg, pois_all_layer)) {
+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)
 
@@ -484,29 +483,29 @@ if(needs_layer(nav_gpkg, pois_all_layer)) {
     # 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_sf(new_POIs, temp_gpkg, poi_moved_layer)
+    write_sf(xWalk, temp_gpkg, poi_xwalk_layer)
     # write out final POIs
-    write_sf(final_POIs, nav_gpkg, pois_all_layer)
+    write_sf(final_POIs, temp_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)
+    write_sf(tmp_POIs, temp_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)
+  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(nav_gpkg, nsegments_layer)) {
+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)
@@ -528,13 +527,12 @@ if(needs_layer(nav_gpkg, nsegments_layer)) {
     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)
+  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(nav_gpkg, pois_all_layer)
+  final_POIs <- read_sf(temp_gpkg, pois_all_layer)
   nhd_Final <- read_sf(nav_gpkg, nhd_flowline)
-  nsegments <- read_sf(nav_gpkg, nsegments_layer)
+  nsegments <- read_sf(temp_gpkg, nsegments_layer)
 }
 
 # Ensure that all the problem POIs have been taken care of
@@ -549,11 +547,13 @@ if(needs_layer(nav_gpkg, final_poi_layer)) {
   if(!"Type_Elev" %in% names(final_POIs)) final_POIs$Type_Elev <- rep(NA_real_, nrow(final_POIs))
   
   #1 Move POIs downstream by category
-  out_HUC12 <- POI_move_down(nav_gpkg, final_POIs, nsegments, filter(nhd_Final, !is.na(POI_ID)), "Type_HUC12", .10)
-  out_gages <- POI_move_down(nav_gpkg, out_HUC12$allPOIs, out_HUC12$segs, out_HUC12$FL, "Type_Gages", .05)
+  out_HUC12 <- POI_move_down(temp_gpkg, final_POIs, nsegments, filter(nhd_Final, !is.na(POI_ID)), "Type_HUC12", .10)
+  out_gages <- POI_move_down(temp_gpkg, out_HUC12$allPOIs, out_HUC12$segs, out_HUC12$FL, "Type_Gages", .05)
   
   # Convert empty strings to NA for handling within R
   out_gages$allPOIs <- mutate_all(out_gages$allPOIs, list(~na_if(.,"")))
+  
+  out_gages$allPOIs$snapped <- out_gages$allPOIs$COMID %in% new_POIs$COMID
 
   nhd_Final <- select(nhd_Final, -POI_ID) %>%
     left_join(st_drop_geometry(out_HUC12$FL) %>%
@@ -561,15 +561,14 @@ if(needs_layer(nav_gpkg, final_poi_layer)) {
 
   # Write out geopackage layer representing POIs for given theme
   write_sf(out_gages$allPOIs, nav_gpkg, final_poi_layer)
-  write_sf(nhd_Final, nav_gpkg, nhd_flowline)
-  write_sf(out_gages$segs, nav_gpkg, nsegments_layer)
+  write_sf(out_gages$segs, temp_gpkg, nsegments_layer)
   
   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(nav_gpkg, nsegments_layer)
+  nsegments <- read_sf(temp_gpkg, nsegments_layer)
 }
 
 mapview(final_POIs, layer.name = "All POIs", col.regions = "red") + 
diff --git a/workspace/05_hyRefactor_AK.Rmd b/workspace/05_hyRefactor_AK.Rmd
index 705a4973e991282241f0a0ab6ef94cc59606e047..aae92ca8d4cd8496e5d571074d185b0650854a23 100644
--- a/workspace/05_hyRefactor_AK.Rmd
+++ b/workspace/05_hyRefactor_AK.Rmd
@@ -375,11 +375,11 @@ refactor_lookup <- dplyr::select(st_drop_geometry(reconciled), ID, member_COMID)
 
 aggregate_lookup_fline <- dplyr::select(st_drop_geometry(agg_cats$fline_sets), ID, set) %>%
   tidyr::unnest(cols = set) %>%
-  dplyr::rename(aggregated_flowline_ID = ID, reconciled_ID = set)
+  dplyr::rename(aggregated_flowpath_ID = ID, reconciled_ID = set)
 
 aggregate_lookup_catchment <- dplyr::select(st_drop_geometry(agg_cats$cat_sets), ID, set) %>%
   tidyr::unnest(cols = set) %>%
-  dplyr::rename(aggregated_catchment_ID = ID, reconciled_ID = set)
+  dplyr::rename(aggregated_divide_ID = ID, reconciled_ID = set)
 
 if(is.character(aggregate_lookup_catchment$reconciled_ID)) 
   aggregate_lookup_catchment$reconciled_ID <- as.integer(aggregate_lookup_catchment$reconciled_ID)
diff --git a/workspace/05_hyRefactor_flines.Rmd b/workspace/05_hyRefactor_flines.Rmd
index 1015377e582a83726265601eee26674f5aa5bfc8..6efa5e6dd93a9006e59cdadbd05fed216e913f49 100644
--- a/workspace/05_hyRefactor_flines.Rmd
+++ b/workspace/05_hyRefactor_flines.Rmd
@@ -43,8 +43,7 @@ POIs <- read_sf(nav_gpkg, final_poi_layer) %>%
 
 # Gages that were collapsed during navigate; might not correspond to their Gage_info COMIDs;
 #   We know DAR for those are acceptable, so will keep those out of the events generation.
-POIs_mv <- read_sf(nav_gpkg, poi_moved_layer) %>% 
-  filter(!is.na(Type_Gages))
+POIs_mv <- filter(POIs, snapped)
 
 # Also need to avoid modification to flowlines immediately downstream of POIs
 #      This can cause some hydrologically-incorrect catchment aggregation
diff --git a/workspace/06-2_aggregate_cats.Rmd b/workspace/06-2_aggregate_cats.Rmd
index a138ee666878e2688befa3d1a437d5d5c6df395e..0e239b8c574188e898b9a0895eab3d28d6d47af5 100644
--- a/workspace/06-2_aggregate_cats.Rmd
+++ b/workspace/06-2_aggregate_cats.Rmd
@@ -51,54 +51,50 @@ reconciled <- st_transform(read_sf(out_refac_gpkg, reconciled_layer), crs)
 # refactored has the original COMIDs and other attributes
 refactored <- st_transform(read_sf(out_refac_gpkg, refactored_layer), crs) 
 
-# Summarize total drainage area for reconciled flowlines
-reconciled_DA <- select(reconciled, ID, toID) %>%
-  left_join(select(st_drop_geometry(divides), ID, area = areasqkm), by = "ID") %>%
-  mutate(ID = as.integer(ID), toID = as.integer(toID),
-         area = ifelse(is.na(area), 0, area)) %>%
-  mutate(rec_DA = calculate_total_drainage_area(.)) %>%
-  left_join(select(lookup, reconciled_ID, NHDPlusV2_COMID), by = c("ID" = "reconciled_ID")) %>%
-  mutate(NHDPlusV2_COMID = as.integer(NHDPlusV2_COMID))
-  
 # Create unified POI/outlet layer
 if(needs_layer(out_agg_gpkg, mapped_outlets_layer)) {
   
-  # subset reconciled flowlines by minimum drainage area for aggcats
-  reconciled_sub <- filter(reconciled, TotDASqKM > 1) 
+  # paths that will for sure be "in network" after aggregate
+  dominant_paths <- filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm)
   
   # Get the end of the reconciled flowlines and bind with data frame
-  reconciled_end <- get_node(filter(reconciled, TotDASqKM > 1) , "end") %>%
-    cbind(st_drop_geometry(filter(reconciled, TotDASqKM > 1) ))
+  reconciled_end <- get_node(dominant_paths, "end") %>%
+    cbind(st_drop_geometry(dominant_paths))
   
   # If there are events on split flowlines
   if(nrow(read_sf(out_refac_gpkg, split_layer)) > 0){
     events <- read_sf(out_refac_gpkg, split_layer) 
-    
-    # Get rows of lookup for split events
-    lookup_events <- filter(lookup, NHDPlusV2_COMID %in% events$COMID)
-    reconciled_events <- filter(reconciled_end, ID %in% lookup_events$reconciled_ID)
 
     # Outlets that are not associated with flowlines with events
-    poi_outlets <- filter(outlets_all, !Type_Gages %in% events$site_no) %>%
+    poi_outlets_mapped <- filter(outlets_all, !Type_Gages %in% events$site_no) %>%
       # Reset Type_Gages attrbute to NA if event exists
       distinct(.keep_all = TRUE) %>%
-      filter(TotDASqKM > 1)
-    
-    # Map non-event outlets ID to reconciled ID 
-    poi_outlets_mapped <- map_outlet_ids(select(st_drop_geometry(poi_outlets), COMID, type), 
-                                     filter(reconciled, TotDASqKM > 1)) %>%
+      filter(TotDASqKM > aggregate_da_thresh_sqkm) %>%
+      st_drop_geometry() %>%
+      select(COMID, type) %>%
+      map_outlet_ids(dominant_paths) %>%
       mutate(COMID = as.integer(COMID)) %>%
       inner_join(select(outlets_all, -type), by = "COMID") %>%
       select(-TotDASqKM) %>%
       st_as_sf() 
-
+    
+    # Summarize total drainage area for reconciled flowlines
+    reconciled_DA <- left_join(select(st_drop_geometry(reconciled), ID, toID),
+                               select(st_drop_geometry(divides), ID, area = areasqkm), by = "ID") %>%
+      mutate(ID = as.integer(ID), toID = as.integer(toID), area = ifelse(is.na(area), 0, area)) %>%
+      mutate(rec_DA = calculate_total_drainage_area(.)) %>%
+      left_join(select(lookup, reconciled_ID, NHDPlusV2_COMID), by = c("ID" = "reconciled_ID")) %>%
+      mutate(NHDPlusV2_COMID = as.integer(NHDPlusV2_COMID)) %>%
+      select(ID, rec_DA, NHDPlusV2_COMID)
+    
     # Pois which share same COMID as events, reset Type_Gages to NA 
     poi_events_mapped <- filter(outlets_all, COMID %in% events$COMID) %>%
       mutate(Type_Gages = NA) %>%
       distinct() %>%
       # Subset to POIs that have other POI flags than Type_Gage
-      filter_at(vars(-c(COMID, TotDASqKM, type, geom)), any_vars(!is.na(.)), TotDASqKM > 1) %>%
-      inner_join(select(st_drop_geometry(reconciled_DA), ID, rec_DA, NHDPlusV2_COMID), 
+      filter_at(vars(-c(COMID, TotDASqKM, type, geom)), any_vars(!is.na(.)), 
+                TotDASqKM > aggregate_da_thresh_sqkm) %>%
+      inner_join(reconciled_DA, 
                  by = c("COMID" = "NHDPlusV2_COMID")) %>% 
       group_by(COMID) %>% 
       # Max DA associated with each orig_COMID will be the proper reconciled ID value
@@ -110,7 +106,7 @@ if(needs_layer(out_agg_gpkg, mapped_outlets_layer)) {
     events_mapped <- select(st_drop_geometry(events), COMID, Type_Gages = site_no) %>%
       #IDs that events are associated with
       mutate(ID = reconciled_end[st_nearest_feature(st_transform(events, st_crs(reconciled_end)), 
-                                                 reconciled_end),]$ID) %>%
+                                                    reconciled_end),]$ID) %>%
       mutate(type = "outlet", COMID = as.character(COMID)) %>%
       # Replace event geometry (off-flowline) with POI geometry (on-flowline)
       inner_join(reconciled_end, by = "ID") %>%
@@ -141,11 +137,11 @@ if(needs_layer(out_agg_gpkg, mapped_outlets_layer)) {
   } else {
     # Map non-event outlets ID into one dataset here with reconciled ID 
     poi_outlets <- distinct(outlets_all, .keep_all = TRUE) %>%
-      filter(TotDASqKM > 1)
+      filter(TotDASqKM > aggregate_da_thresh_sqkm)
     
     # Map non-event outlets ID to reconciled ID 
     outlets_POI<- map_outlet_ids(select(st_drop_geometry(poi_outlets), COMID, type), 
-                                     filter(reconciled, TotDASqKM > 1)) %>%
+                                 dominant_paths) %>%
       inner_join(select(outlets_all, -type), by = "COMID") %>%
       select(-TotDASqKM) %>%
       st_as_sf()
@@ -185,7 +181,7 @@ if(needs_layer(out_agg_gpkg, agg_cats_layer)){
   agg_cats <- aggregate_catchments(flowpath = reconciled,
                                    divide = divides,
                                    outlets = outlets,
-                                   da_thresh = 1,
+                                   da_thresh = aggregate_da_thresh_sqkm,
                                    only_larger = TRUE,
                                    post_mortem_file = cache_split)
     
@@ -232,12 +228,12 @@ refactor_lookup <- dplyr::select(st_drop_geometry(reconciled), ID, member_COMID)
 aggregate_lookup_fline <- dplyr::select(st_drop_geometry(agg_cats$fline_sets), ID, set) %>%
   dplyr::mutate(set = strsplit(set, ",")) %>%
   tidyr::unnest(cols = set) %>%
-  dplyr::rename(aggregated_flowline_ID = ID, reconciled_ID = set)
+  dplyr::rename(aggregated_flowpath_ID = ID, reconciled_ID = set)
 
 aggregate_lookup_catchment <- dplyr::select(st_drop_geometry(agg_cats$cat_sets), ID, set) %>%
   dplyr::mutate(set = strsplit(set, ",")) %>%
   tidyr::unnest(cols = set) %>%
-  dplyr::rename(aggregated_catchment_ID = ID, reconciled_ID = set)
+  dplyr::rename(aggregated_divide_ID = ID, reconciled_ID = set)
 
 if(is.character(aggregate_lookup_catchment$reconciled_ID)) aggregate_lookup_catchment$reconciled_ID <- as.integer(aggregate_lookup_catchment$reconciled_ID)
 if(is.character(aggregate_lookup_fline$reconciled_ID)) aggregate_lookup_fline$reconciled_ID <- as.integer(aggregate_lookup_fline$reconciled_ID)
diff --git a/workspace/07-1_merge.Rmd b/workspace/07-1_merge.Rmd
index bb91278db1708bbcdaf7eaef73296a905712a7fd..f9f4cb63687c906de9b9d447e3835eba1d569cbc 100644
--- a/workspace/07-1_merge.Rmd
+++ b/workspace/07-1_merge.Rmd
@@ -35,10 +35,10 @@ if(needs_layer(gf_gpkg, agg_cats_layer)) {
   POIs <- read_sf(nav_gpkg,  final_poi_layer) %>% 
     st_transform(crs) 
   
-  write_sf(rfc_gpkg, final_poi_layer)
-  
-  write_sf(gf_gpkg, final_poi_layer)
-  
+  # write_sf(POIs, rfc_gpkg, final_poi_layer)
+  # 
+  # write_sf(POIs, gf_gpkg, final_poi_layer)
+  # 
   merged_layers <- merge_refactor(rpu_codes$rpuid, rpu_vpu_out, 
                                   lookup_table_refactor, 
                                   reconciled_layer, 
@@ -47,26 +47,27 @@ if(needs_layer(gf_gpkg, agg_cats_layer)) {
                                   agg_cats_layer,
                                   mapped_outlets_layer)
   
-  sf::write_sf(merged_layers[[lookup_table_refactor]], rfc_gpkg, lookup_table_refactor)
+  sf::write_sf(select(merged_layers[[lookup_table_refactor]], 
+                      -aggregated_flowpath_ID, 
+                      -aggregated_divide_ID), 
+               rfc_gpkg, lookup_table_refactor)
   
-  sf::write_sf(merged_layers[[reconciled_layer]], rfc_gpkg, reconciled_layer)
+  sf::write_sf(merged_layers[[reconciled_layer]], 
+               rfc_gpkg, reconciled_layer)
   
-  sf::write_sf(merged_layers[[divide_layer]], rfc_gpkg, divide_layer)
+  sf::write_sf(merged_layers[[divide_layer]], 
+               rfc_gpkg, divide_layer)
   
   sf::write_sf(merged_layers[[mapped_outlets_layer]], rfc_gpkg, mapped_outlets_layer)
   
-    sf::write_sf(merged_layers[[mapped_outlets_layer]], gf_gpkg, mapped_outlets_layer)
+  sf::write_sf(merged_layers[[mapped_outlets_layer]], gf_gpkg, mapped_outlets_layer)
 
   sf::write_sf(merged_layers[[agg_cats_layer]], gf_gpkg, agg_cats_layer)
     
   sf::write_sf(merged_layers[[agg_fline_layer]], gf_gpkg, agg_fline_layer)  
   
-} else {
-  # Load layers described above 
-  nhd <- read_sf(gf_gpkg, nhd_flowline)
-  cats <- read_sf(gf_gpkg, nhd_catchment)
-  reconciled <- read_sf(gf_gpkg, reconciled_layer)
-  divides <- read_sf(gf_gpkg, divide_layer)
-  lookup <- read_sf(gf_gpkg, lookup_table_refactor)
+  sf::write_sf(merged_layers[[lookup_table_refactor]], 
+               gf_gpkg, lookup_table_refactor)
+  
 }
 ```
diff --git a/workspace/07-2_NonDend.Rmd b/workspace/07-2_NonDend.Rmd
index 59548818df57dbc6377b56548522bbcc7b54d454..ea46205c3cf8bcf957dac5ed6077213c51525b8f 100644
--- a/workspace/07-2_NonDend.Rmd
+++ b/workspace/07-2_NonDend.Rmd
@@ -139,16 +139,16 @@ if(needs_layer(ND_gpkg, divides_nd)){
   
   # Bind divide with lookup, identify aggregation step
   divides_lu <- divides %>%
-    left_join(distinct(select(lookup, reconciled_ID, aggregated_catchment_ID)), 
+    left_join(distinct(select(lookup, reconciled_ID, aggregated_divide_ID)), 
               by = c("ID" = "reconciled_ID")) %>%
-    rename(POI_ID = aggregated_catchment_ID) %>%
+    rename(POI_ID = aggregated_divide_ID) %>%
     # attribute that hyrefactor was the step used to aggregate these catchments
     mutate(aggStep = ifelse(!is.na(POI_ID), "hyRef", NA))
 
   # Identify cats witin the full catchment set not refactored or aggregated, add
   #          to divides data frame
   cats_miss <- cats %>%
-    left_join(select(lookup, NHDPlusV2_COMID, POI_ID = aggregated_catchment_ID),
+    left_join(select(lookup, NHDPlusV2_COMID, POI_ID = aggregated_divide_ID),
               by = c("FEATUREID" =  "NHDPlusV2_COMID")) %>%
     filter(is.na(POI_ID)) %>%
     select(FEATUREID, HUC_12_int, intArea, AreaHUC12, POI_ID) %>%
diff --git a/workspace/R/config.R b/workspace/R/config.R
index 36a3dc2fbde159266a27cc5bafb8d84e1365b5c1..1bc57dfaedc4000affe509bcab9d6b83250bd89b 100644
--- a/workspace/R/config.R
+++ b/workspace/R/config.R
@@ -25,6 +25,7 @@ if(!exists("rpu_code")) {
     vpu <- get_rpu_dependent_vars()
     rpu_code <- vpu$r
     vpu <- vpu$v
+    rpu_codes <- rpu_vpu[rpu_vpu$vpuid %in% vpu, ]
   }
 }
 
@@ -69,7 +70,7 @@ travt_diff <- 24 # Max number of hours between adjacent POIS
 
 xwalk_layer <- paste0("HUC12_nhd") # HUC12 - nhdcat crosswalk, built in Nav for VPU 20
 nav_poi_layer <- paste0("POIs_tmp_", vpu) # Rolling Nav POI layer added to/modified througout nav workflow
-WBs_layer <-  paste0("WB_", vpu) # Waterbodies within VPU
+WBs_layer <-  paste0("WB") # Waterbodies within VPU
 poi_moved_layer <- paste0("POIs_mv_", vpu) # POIs moved from original COMID assignment
 nsegments_layer <- paste0("nsegment_", vpu) # Minimally-sufficient network dissolved by POI_ID
 pois_all_layer <- paste0("POIs_", vpu) # All POIs binded together
@@ -81,6 +82,8 @@ split_meters <- 10000
 combine_meters <- 1000
 min_da_km <- 20
 
+aggregate_da_thresh_sqkm <- 1
+
 # parallel factor for catchment reconciliation.
 para_reconcile <- 2
 para_split_flines <- 2
@@ -101,6 +104,7 @@ lookup_table_refactor <- "lookup_table"
 
 # output geopackage file names
 nav_gpkg <- file.path("cache", paste0("reference_", vpu,".gpkg"))
+temp_gpkg <- file.path("cache", paste0("nav_", vpu,".gpkg"))
 
 rfc_gpkg <- file.path("cache", paste0("refactor_", vpu, ".gpkg"))
 
@@ -108,7 +112,7 @@ gf_gpkg <- file.path("cache", paste0("GF_", vpu, ".gpkg"))
 gf_gpkg_conus <- "cache/reference_CONUS.gpkg"
 
 # Defined during NonDend.Rmd
-ND_gpkg <- file.path("cache", paste0("ND_", vpu,".gpkg"))
+ND_gpkg <- file.path("temp", paste0("ND_", vpu,".gpkg"))
 
 divides_xwalk <- paste0("divides_nhd")
 HRU_layer <- paste0("all_cats")
diff --git a/workspace/R/utils.R b/workspace/R/utils.R
index 9e318a8688c68096e72931b8fe576efea7e802db..61cf195a5c0190c9a0fe04f8dd6d26ae8061c882 100644
--- a/workspace/R/utils.R
+++ b/workspace/R/utils.R
@@ -366,10 +366,10 @@ merge_refactor <- function(rpus, rpu_vpu_out,
     out[[rpu]][[lookup_table_refactor]] <- out[[rpu]][[lookup_table_refactor]] %>%
       left_join(update_id, by = c("reconciled_ID" = "old_ID")) %>%
       select(-reconciled_ID, reconciled_ID = ID) %>%
-      left_join(update_id, by = c("aggregated_flowline_ID" = "old_ID")) %>%
-      select(-aggregated_flowline_ID, aggregated_flowline_ID = ID) %>%
-      left_join(update_id, by = c("aggregated_catchment_ID" = "old_ID")) %>%
-      select(-aggregated_catchment_ID, aggregated_catchment_ID = ID) %>%
+      left_join(update_id, by = c("aggregated_flowpath_ID" = "old_ID")) %>%
+      select(-aggregated_flowpath_ID, aggregated_flowpath_ID = ID) %>%
+      left_join(update_id, by = c("aggregated_divide_ID" = "old_ID")) %>%
+      select(-aggregated_divide_ID, aggregated_divide_ID = ID) %>%
       left_join(sf::st_drop_geometry(select(out[[rpu]][[reconciled_layer]], 
                                             ID = newID, LevelPathID)),
                 by = c("reconciled_ID" = "ID"))
@@ -393,29 +393,45 @@ merge_refactor <- function(rpus, rpu_vpu_out,
                                    update_id, 
                                    by = c("old_ID"))
       
-      
-      out[[rpu]][[l]] <- left_join(rename(out[[rpu]][[l]], old_toID = toID),
-                                   rename(update_id, toID = ID), 
-                                   by = c("old_toID" = "old_ID"))
-      
-     sets <- select(st_drop_geometry(out[[rpu]][[l]]), ID, old_ID, old_set = set) %>%
-       mutate(old_set = strsplit(old_set, ",")) %>%
-       tidyr::unnest(cols = old_set) %>%
-       mutate(old_set = as.integer(old_set)) %>%
-       left_join(rename(update_id, set = ID), by  = c("old_set" = "old_ID")) %>%
-       select(-old_set, -old_ID)
-       
-     sets <- split(sets$set, sets$ID)
-     
-     sets <- data.frame(ID = as.integer(names(sets))) %>%
-       mutate(set = unname(sets))
-     
-     out[[rpu]][[l]] <- left_join(select(out[[rpu]][[l]], -set),
-                                  sets, by = "ID")
-     
-     out[[rpu]][[l]] <- select(out[[rpu]][[l]], ID, toID, set)
-     
-     out[[rpu]][[l]]$set <- sapply(out[[rpu]][[l]]$set, paste, collapse = ",")
+      if("toID" %in% names(out[[rpu]][[l]])) {
+        out[[rpu]][[l]] <- left_join(rename(out[[rpu]][[l]], old_toID = toID),
+                                     rename(update_id, toID = ID), 
+                                     by = c("old_toID" = "old_ID"))
+        
+        sets <- select(st_drop_geometry(out[[rpu]][[l]]), ID, old_ID, old_set = set) %>%
+          mutate(old_set = strsplit(old_set, ",")) %>%
+          tidyr::unnest(cols = old_set) %>%
+          mutate(old_set = as.integer(old_set)) %>%
+          left_join(rename(update_id, set = ID), by  = c("old_set" = "old_ID")) %>%
+          select(-old_set, -old_ID)
+        
+        sets <- split(sets$set, sets$ID)
+        
+        sets <- data.frame(ID = as.integer(names(sets))) %>%
+          mutate(set = unname(sets))
+        
+        out[[rpu]][[l]] <- left_join(select(out[[rpu]][[l]], -set),
+                                     sets, by = "ID")
+        
+        out[[rpu]][[l]] <- select(out[[rpu]][[l]], -old_ID, -old_toID)
+        
+        out[[rpu]][[l]] <- out[[rpu]][[l]][c("ID", "toID", "set", 
+                                             names(out[[rpu]][[l]])[!names(out[[rpu]][[l]]) %in% 
+                                                                      c("ID", "toID", "set", 
+                                                                        attr(out[[rpu]][[l]], "sf_column"))],
+                                             attr(out[[rpu]][[l]], "sf_column"))]
+        
+        out[[rpu]][[l]]$set <- sapply(out[[rpu]][[l]]$set, paste, collapse = ",")
+      } else {
+        
+        out[[rpu]][[l]] <- select(out[[rpu]][[l]], -old_ID)
+        
+        out[[rpu]][[l]] <- out[[rpu]][[l]][c("ID", 
+                                             names(out[[rpu]][[l]])[!names(out[[rpu]][[l]]) %in% 
+                                                                      c("ID", 
+                                                                        attr(out[[rpu]][[l]], "sf_column"))],
+                                             attr(out[[rpu]][[l]], "sf_column"))]
+      }
     }
     
   }
diff --git a/workspace/README.md b/workspace/README.md
index e52b4af95ab3e814ce7e4ec98b35638e4e9cfb8e..314a40a8567e959080708931740fe34912dedea5 100644
--- a/workspace/README.md
+++ b/workspace/README.md
@@ -21,32 +21,32 @@ This is where scripts and data live.
 - nwm_enhd__nhdplusatts.csv: updated network attributes for NWM ids.
 - nwm_network.gpkg: NWM flowlines with updated network attributes.
 - reference*: by VPU files containing all reference fabric content at the native source resolution.
-  - nhd_flowline
-  - nhd_catchment
-  - final_POIs*
-  - poi_xwalk
+  - reference_flowlines
+  - reference_catchments
+  - POIs
   - WB
 - national_reference_POI.gpkg: National set of POIs determined by `02_NHD_Navigate.Rmd.`
-- refactor*: by RPU files containing a set of refactored network artifacts preserving all network junctions from the source data.
-  - nhd_flowline
-  - nhd_catchment
+- refactor*: by VPU refactored outputs.
+  - refactored_flowpaths
+  - refactored_divides
+  - mapped_POIs
+  - lookup_table
+- GF*: by VPU files containing aggregated outputs.
+  - aggregated_flowpaths
+  - aggregated_divides
+  - mapped_POIs
+  - lookup_table
+- *refactor: by RPU files containing a set of refactored network artifacts preserving all network junctions from the source data.
+  - reconciled_flowpaths
+  - reconciled_divides
   - split_events
   - outlets
-  - reconciled
   - lookup_table
-- aggregate*: by RPU files containing a set of aggregated network artifacts resolving, at a minimum the POIs in national_reference_POI.gpkg.
+- *aggregate: by RPU files containing a set of aggregated network artifacts resolving, at a minimum the POIs in national_reference_POI.gpkg.
   - agg_cats
   - agg_flines
   - mapped_outlets
   - lookup_table
-- GF*: by VPU files containing contents from previous by RPU steps.
-  - agg_cats
-  - agg_files
-  - mapped_outlets
-  - lookup_table
-  - reconciled
-  - nhd_flowline
-  - nhd_catchment
 - ND*: by VPU files containing contents "non-dendritic" workflow that seeks to identify areas of non-networked catchments and how to connect them to the network appropriately.
   - all_cats
   - divides_nd
diff --git a/workspace/hyfabric_0.5.4.tar.gz b/workspace/hyfabric_0.5.4.tar.gz
deleted file mode 100644
index 86456c90fedfb77623c1daaddada4a34b52c167a..0000000000000000000000000000000000000000
Binary files a/workspace/hyfabric_0.5.4.tar.gz and /dev/null differ
diff --git a/workspace/runners/run_one_R.R b/workspace/runners/run_one_R.R
index a720a1e00f4d90fb94d139f00d763f7133dbf187..ab7ab168798275b730b3a39fb54ec02722bdf6b0 100644
--- a/workspace/runners/run_one_R.R
+++ b/workspace/runners/run_one_R.R
@@ -2,7 +2,9 @@
 ## workflows. These should be run prior to checking in any changes to these 
 ## workflows.
 
-VPU <- vpu_codes <- "12"
+VPU <- vpu_codes <- "01"
+
+# try 15a 17a
 
 source("R/config.R")
 source("R/utils.R")
@@ -48,8 +50,10 @@ rmarkdown::render("07-2_NonDend.Rmd",
 
 sb_reference <- "61295190d34e40dd9c06bcd7"
 sb_refactor <- "61fbfdced34e622189cb1b0a"
+sb_gf <- "60be1502d34e86b9389102cc"
 
 sbtools::authenticate_sb()
 
 sbtools::item_replace_files(sb_reference, nav_gpkg)
-sbtools::item_replace_files(sb_refactor, gf_gpkg)
+sbtools::item_replace_files(sb_refactor, rfc_gpkg)
+sbtools::item_replace_files(sb_gf, gf_gpkg)