diff --git a/workspace/02_NHD_navigate.Rmd b/workspace/02_NHD_navigate.Rmd
index c86fd93f59833a1658710d97910da488f985c7e1..93ed09c5d7c9b57674746d07749f2e88bf7a0c3c 100644
--- a/workspace/02_NHD_navigate.Rmd
+++ b/workspace/02_NHD_navigate.Rmd
@@ -48,7 +48,7 @@ atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.r
 #  Derive or load HUC12 POIs
 if(needs_layer(temp_gpkg, nav_poi_layer)) {
   
-   nhd <- read_sf(ref_gpkg, nhd_flowline)
+   nhd <- read_sf(nav_gpkg, nhd_flowline)
    try(nhd <- select(nhd, -c(minNet, WB, struct_POI, struct_Net, POI_ID)))
   
   # HUC02 includes some 
@@ -78,7 +78,7 @@ if(needs_layer(temp_gpkg, nav_poi_layer)) {
 } else {
   # Load HUC12 POIs as the tmpPOIs if they already exist
   tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer) 
-  nhd <- read_sf(ref_gpkg, nhd_flowline)
+  nhd <- read_sf(nav_gpkg, nhd_flowline)
 }
 
 mapview(filter(tmp_POIs, Type_HUC12 != ""), layer.name = "HUC12 POIs", col.regions = "red") 
@@ -186,7 +186,7 @@ if(all(is.na(tmp_POIs$Type_Term))) {
   write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
   tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-  nhd <- read_sf(ref_gpkg, nhd_flowline)
+  nhd <- read_sf(nav_gpkg, nhd_flowline)
 }
 
 mapview(filter(tmp_POIs, !is.na(Type_Term)), layer.name = "Terminal POIs", col.regions = "blue") 
@@ -221,7 +221,7 @@ if(all(is.na(tmp_POIs$Type_Conf))) {
   
   write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
-  nhd <- read_sf(ref_gpkg, nhd_flowline)
+  nhd <- read_sf(nav_gpkg, nhd_flowline)
   tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
 }
 
@@ -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, sf::st_crs(nhd)), ref_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))
@@ -338,7 +338,7 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
   write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
   tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-  nhd <- read_sf(ref_gpkg, nhd_flowline)
+  nhd <- read_sf(nav_gpkg, nhd_flowline)
 }
 
 mapview(filter(tmp_POIs, !is.na(Type_WBIn)), layer.name = "Waterbody inlet POIs", col.regions = "blue") +
@@ -406,7 +406,7 @@ if(nrow(elev_pois_split) > 0) {
   write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
     
   tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer) 
-  nhd <- read_sf(ref_gpkg, nhd_flowline)
+  nhd <- read_sf(nav_gpkg, nhd_flowline)
 }
 
 if(!all(is.na(tmp_POIs$Type_Elev)))
@@ -459,7 +459,7 @@ if(all(is.na(tmp_POIs$Type_Travel))) {
   write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 }else{
   tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-  nhd <- read_sf(ref_gpkg, nhd_flowline)
+  nhd <- read_sf(nav_gpkg, nhd_flowline)
 }
 
 mapview(filter(tmp_POIs, Type_Travel == 1), layer.name = "Travel time break POIs", col.regions = "blue") 
@@ -531,7 +531,7 @@ if(needs_layer(temp_gpkg, nsegments_layer)) {
 } else {
   # Read in NHDPlusV2 flowline simple features and filter by vector processing unit (VPU)
   final_POIs <- read_sf(temp_gpkg, pois_all_layer)
-  nhd_Final <- read_sf(ref_gpkg, nhd_flowline)
+  nhd_Final <- read_sf(nav_gpkg, nhd_flowline)
   nsegments <- read_sf(temp_gpkg, nsegments_layer)
 }
 
@@ -542,7 +542,7 @@ print (paste0(nrow(sub[sub$AreaSqKM == 0,]), " POIs on flowlines with no local d
 
 ```{r POI Collapse}
 #  Load data
-if(needs_layer(ref_gpkg, final_poi_layer)) {
+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))
   
@@ -560,14 +560,14 @@ if(needs_layer(ref_gpkg, final_poi_layer)) {
                 select(COMID, POI_ID), by = "COMID")
 
   # Write out geopackage layer representing POIs for given theme
-  write_sf(out_gages$allPOIs, ref_gpkg, final_poi_layer)
+  write_sf(out_gages$allPOIs, nav_gpkg, final_poi_layer)
   write_sf(out_gages$segs, temp_gpkg, nsegments_layer)
   
   nsegments <- out_gages$segs
   final_POIs <- out_gages$allPOIs
 } else {
-  final_POIs <- read_sf(ref_gpkg, final_poi_layer)
-  nhd_Final <- read_sf(ref_gpkg, nhd_flowline)
+  final_POIs <- read_sf(nav_gpkg, final_poi_layer)
+  nhd_Final <- read_sf(nav_gpkg, nhd_flowline)
   nsegments <- read_sf(temp_gpkg, nsegments_layer)
 }
 
diff --git a/workspace/07-1_merge.Rmd b/workspace/07-1_merge.Rmd
index 84b8d9559c36a28e65c572970be9e6644088460c..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(POIs, rfc_gpkg, final_poi_layer)
-  
-  write_sf(POIs, 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, 
@@ -66,12 +66,8 @@ if(needs_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/R/utils.R b/workspace/R/utils.R
index 810dfb237f3d3e556adda2421aad04c209a948e0..61cf195a5c0190c9a0fe04f8dd6d26ae8061c882 100644
--- a/workspace/R/utils.R
+++ b/workspace/R/utils.R
@@ -393,35 +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]], -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 = ",")
+      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 437a13e5ab4cad7a14220bb2a40284aadce6db6b..314a40a8567e959080708931740fe34912dedea5 100644
--- a/workspace/README.md
+++ b/workspace/README.md
@@ -31,6 +31,11 @@ This is where scripts and data live.
   - 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
@@ -42,14 +47,6 @@ This is where scripts and data live.
   - 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/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)