diff --git a/hyfabric/DESCRIPTION b/hyfabric/DESCRIPTION
index 4a5916087a9fb348f937c09555cb739df6d4cbc3..54a6cef14a392446b78faaaa86b634662959e2d6 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.2
+Version: 0.5.3
 Authors@R: c(person(given = "David", 
                     family = "Blodgett", 
                     role = c("aut", "cre"), 
diff --git a/hyfabric/inst/extdata/rpu_vpu.csv b/hyfabric/inst/extdata/rpu_vpu.csv
new file mode 100644
index 0000000000000000000000000000000000000000..2c8101ec92c9142b7b0667a18544b7402db8e721
--- /dev/null
+++ b/hyfabric/inst/extdata/rpu_vpu.csv
@@ -0,0 +1,60 @@
+rpuid,vpuid
+18c,18
+18b,18
+18a,18
+14a,14
+14b,14
+15b,15
+15a,15
+16b,16
+16a,16
+04b,04
+04d,04
+04c,04
+04a,04
+02b,02
+02a,02
+05d,05
+05c,05
+05a,05
+05b,05
+06a,06
+07c,07
+07b,07
+07a,07
+03g,08
+08b,08
+08a,08
+10c,10L
+10d,10L
+10b,10L
+10a,10L
+10i,10U
+10h,10U
+10g,10U
+10f,10U
+10e,10U
+11c,11
+11d,11
+11b,11
+11a,11
+01a,01
+17b,17
+17a,17
+17c,17
+17d,17
+13a,13
+13b,13
+13c,13
+13d,13
+03a,03N
+03b,03N
+03c,03S
+03d,03S
+03e,03W
+03f,03W
+09a,09
+12a,12
+12b,12
+12c,12
+12d,12
diff --git a/workspace/00_get_data.Rmd b/workspace/00_get_data.Rmd
index 555397477e7f2436e5921d1f56285526b946cae0..93a4e2e1590dbbd9baf2b4b1f406edb58a228efb 100644
--- a/workspace/00_get_data.Rmd
+++ b/workspace/00_get_data.Rmd
@@ -561,10 +561,11 @@ if(!file.exists(out$new_nwm_nhdp_atts)) {
 out_list <- c(out_list, out)
 ```
 
-```{r Gages_GFv1.1}
+```{r GFv1.1}
 GFv11_dir <- file.path(data_dir, "GFv11")
 out <- list(GFv11_gages_lyr = file.path(data_dir, "GFv11/GFv11_gages.rds"),
-            GFv11_gdb = file.path(GFv11_dir, "GFv1.1.gdb"))
+            GFv11_gdb = file.path(GFv11_dir, "GFv1.1.gdb"),
+            GFv11_tgf = file.path(GFv11_dir, "TGF.gdb"))
 
 # Download the GFv1.1 geodatabase
 if(!dir.exists(GFv11_dir)) {
@@ -574,10 +575,18 @@ if(!dir.exists(GFv11_dir)) {
     sb <- sbtools::authenticate_sb()
   
   sbtools::item_file_download("5e29d1a0e4b0a79317cf7f63", names = "GFv1.1.gdb.zip", 
-                              destinations = file.path(GFv11_dir, "GFv1.1.gdb.zip"), session = sb)
+                              destinations = file.path(GFv11_dir, "GFv1.1.gdb.zip"), 
+                              session = sb)
+  
+  tgf_f <- file.path(GFv11_dir, "TGF.gdb.zip")
+  
+  sbtools::item_file_download("5d967365e4b0c4f70d113923", names = basename(tgf_f), 
+                              destinations = tgf_f)
   
   unzip(file.path(GFv11_dir, "GFv1.1.gdb.zip"), exdir = GFv11_dir)
+  unzip(tgf_f, exdir = GFv11_dir)
   
+  file.remove(tgf_f)
   # Extract gages
   GFv11_gages <- read_sf(out$GFv11_gdb, "POIs_v1_1") %>%
     filter(Type_Gage != 0) 
@@ -613,15 +622,35 @@ out_list <- c(out_list, g2_out)
 
 ```{r updated_flowlines}
 
-out <- list(new_nhdp_rds = file.path(out_list$nhdplus_dir, (sb_f <- "nhdplus_flowline_update.rds")))
+out <- list(ref_flowline = file.path(out_list$nhdplus_dir, (sb_f <- "reference_flowline.gpkg")))
 
-if(!file.exists(out$new_nhdp_rds)) {
+if(!file.exists(out$ref_flowline)) {
   if(is.null(sbtools::current_session()))
     authenticate_sb()
   
+  fs <- sbtools::item_list_files("61295190d34e40dd9c06bcd7")
+  
+  if(!sb_f %in% fs$fname) {
+    
+    stop("you have to run this part interactively")
+    
+    staged_nhd <- stage_national_data(nhdplus_data = out_list$nhdplus_gdb,
+                                      output_path  = out_list$nhdplus_dir)
+    
+    # This will generate 
+    fix_headwaters(staged_nhd$flowline, 
+                   out$ref_flowline,
+                   new_atts = out_list$new_nhdp_atts, 
+                   nhdpath = out_list$nhdplus_gdb)
+    
+    sbtools::authenticate_sb()
+    sbtools::item_append_files("5dcd5f96e4b069579760aedb", out$ref_flowline)
+    
+  }
+    
   sbtools::item_file_download("5dcd5f96e4b069579760aedb",
                               names = sb_f,
-                              destinations = out$new_nhdp_rds)
+                              destinations = out$ref_flowline)
 }
 
 out_list <- c(out_list, out)
diff --git a/workspace/01_Gage_Selection.Rmd b/workspace/01_Gage_Selection.Rmd
index c380f56a4ae1e81ae0894bff33e9a4a751fbf770..50a29f9ce3e280aaf1b88eb3dd7ec1065ef13861 100644
--- a/workspace/01_Gage_Selection.Rmd
+++ b/workspace/01_Gage_Selection.Rmd
@@ -4,11 +4,6 @@ output: html_document
 ---
 This notebook builds a final set of gages for GFv2.0 from data sources assembled in 00_gage_info.rmd
 
-Retrieves and pre-processes streamgage location information for GFv2.0 from
-1. [Geospatial Fabric Version 1.1](https://www.sciencebase.gov/catalog/item/5e29b87fe4b0a79317cf7df5)
-2. GAGES-III (in review)
-3. [NHDPlus GageLoc](https://www.sciencebase.gov/catalog/item/577445bee4b07657d1a991b6)
-
 ```{r setup_0, echo=FALSE, cache=FALSE}
 knitr::opts_chunk$set(
   collapse = TRUE,
@@ -61,7 +56,7 @@ if(file.exists("cache/dropped_gages.csv")){
 }
 
 # Update location index information from other sources where necessary from NLDI reference gage set
-ref_gages <- jsonlite::fromJSON("https://github.com/dblodgett-usgs/ref_gages/releases/download/v0.5/usgs_nldi_gages.geojson", 
+ref_gages <- jsonlite::fromJSON("https://github.com/internetofwater/ref_gages/releases/download/v0.5/usgs_nldi_gages.geojson", 
                                 flatten = TRUE)$features %>%
   select(-geometry.type, -geometry.coordinates, ref_site_no = properties.id, 
          ref_RC = properties.nhdpv2_REACHCODE, ref_RM = properties.nhdpv2_REACH_measure, 
@@ -327,12 +322,12 @@ gage_final_ranked <- gage_final %>%
   # Arrange by number of days
   arrange(total_obs, .by_group = TRUE) %>%
   # Populate by rank (1 - lowest number of days to n(highest))
-  mutate(gage_score = ifelse(n > 1, seq_along(row_number()), 1)) %>%
+  mutate(gage_score = ifelse(n > 1, row_number(), 1)) %>%
   # If a reference gage is one of gages, bump that to the top
   #mutate(gage_score = ifelse(n > 1,
   #                        ifelse(is.na(gagesII_cl), gage_score, max(gage_score) + 1), gage_score)) %>%
   mutate(gage_score = ifelse(n > 1,
-                          ifelse(gagesII_cl == "Ref", max(gage_score) + 1, gage_score), gage_score)) %>%
+                          ifelse(!is.na(gagesII_cl) && gagesII_cl == "Ref", max(gage_score) + 1, gage_score), gage_score)) %>%
   ungroup() %>%
   # Some gages have multiple gages with same days (inlets/outlets around dam)
   # This further prioritizes to the gage with the smallest reach measure (furthest downstream)
diff --git a/workspace/01_NHD_prep.Rmd b/workspace/01_NHD_prep.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..4d85989798689d52d4dc2c0ca8c3ada0b2f5e498
--- /dev/null
+++ b/workspace/01_NHD_prep.Rmd
@@ -0,0 +1,159 @@
+---
+title: "NHD Prep"
+output: html_document
+editor_options:
+  chunk_output_type: console
+---
+Prepares the NHD network for navigate and refactor. Adds base layers to the reference fabric.
+
+```{r setup_rmd, echo=FALSE, cache=FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  fig.width=6, 
+  fig.height=4,
+  cache=FALSE)
+```
+
+```{r setup}
+# Load custom functions
+source("R/utils.R")
+source("R/NHD_navigate.R")
+
+# Load Configuration of environment
+source("R/config.R")
+
+if(!exists("vpu_codes")) {
+  vpu_codes <- unique(rpu_vpu$vpuid)
+} else {
+  rpu_vpu <- filter(rpu_vpu, vpuid %in% vpu_codes)
+}
+
+all_rpu_codes <- unique(rpu_vpu$rpuid)
+full_cat_table <- readRDS(data_paths$fullcats_table)
+
+out_vpu <- rep(list(list()), length(vpu_codes))
+names(out_vpu) <- vpu_codes
+
+out_rpu <- rep(list(list()), length(all_rpu_codes))
+names(out_rpu) <- all_rpu_codes
+
+fline <- sf::read_sf(data_paths$ref_flowline)
+catchment <- readRDS(staged_nhd$catchment)
+
+for(VPU in vpu_codes) {
+  
+  rm(rpu_code)
+  
+  source("R/config.R")
+  
+  if (needs_layer(nav_gpkg, nhd_flowline)) {
+    
+    # Subset NHD by VPU
+    nhd <- subset_vpu(fline, vpu) 
+    
+    # there are some duplicates with enhd
+    nhd <- nhd %>%
+      group_by(COMID) %>%
+      filter(row_number() == 1) %>%
+      ungroup()
+    
+    # Filter and write dendritic/non-coastal subset to gpkg
+    # This will be iterated over to produce the final network after POIs identified
+    zero_order <- filter(nhd, TerminalFl == 1 & TotDASqKM < min_da_km)
+    
+    non_dend <- unique(unlist(lapply(zero_order$COMID, NetworkNav, nhdDF = st_drop_geometry(nhd))))
+    
+    # Add fields to note dendritic and POI flowlines
+    nhd <- nhd %>% 
+      mutate(dend = ifelse(!COMID %in% non_dend, 1, 0),
+             poi = ifelse(!COMID %in% non_dend & TotDASqKM >= min_da_km, 1, 0)) 
+    
+    write_sf(nhd, nav_gpkg, nhd_flowline)
+    
+    cats_rpu <- full_cat_table %>%
+      filter(RPUID %in% rpu_codes$rpuid)
+    
+    cats <- catchment %>% 
+      filter(FEATUREID %in% 
+               unique(c(nhd$COMID, 
+                        full_cat_table$FEATUREID[full_cat_table$RPUID %in% rpu_codes$rpuid]))) %>%
+      mutate(hy_cats = ifelse(FEATUREID %in% nhd$COMID, 1, 0),
+             full_cats = ifelse(FEATUREID %in% cats_rpu$FEATUREID, 1, 0)) %>%
+      filter(full_cats == 1 | hy_cats == 1) %>%
+      st_sf()
+    
+    write_sf(cats, nav_gpkg, nhd_catchment)
+    
+  } else {
+    
+    nhd <- read_sf(nav_gpkg, nhd_flowline)
+    cats <- read_sf(nav_gpkg, nhd_catchment)
+    
+  }
+  
+  out_vpu[[VPU]] <- nhd %>%
+    sf::st_drop_geometry() %>%
+    select(COMID, toCOMID) %>%
+    filter(!toCOMID %in% COMID & !toCOMID == 0)
+  
+  for(rpu_code in rpu_codes$rpuid) {
+    source("R/config.R")
+    
+    if(needs_layer(out_refac_gpkg, nhd_flowline)) {
+      nhd_sub <- subset_rpu(nhd, rpu_code, run_make_standalone = TRUE) %>%
+        st_sf()
+      
+      write_sf(nhd_sub, out_refac_gpkg, nhd_flowline)
+      
+      cats_rpu <- full_cat_table %>%
+        filter(RPUID == rpu_code)
+      
+      cats %>%
+        mutate(hy_cats = ifelse(FEATUREID %in% nhd_sub$COMID, 1, 0),
+               full_cats = ifelse(FEATUREID %in% cats_rpu$FEATUREID, 1, 0)) %>%
+        filter(full_cats == 1 | hy_cats == 1) %>%
+        st_sf() %>%
+        write_sf(out_refac_gpkg, nhd_catchment)
+      
+    } else {
+      nhd_sub <- read_sf(out_refac_gpkg, nhd_flowline)
+    }
+    
+    out_rpu[[rpu_code]] <- nhd_sub %>%
+      sf::st_drop_geometry() %>%
+      select(COMID, toCOMID) %>%
+      filter(!toCOMID %in% COMID & !toCOMID == 0)
+  }
+}
+```
+
+```{r}
+
+if(!file.exists("cache/rpu_vpu_out.csv")) {
+  make_df <- function(x, d, n) {
+    y <- d[[x]]
+    nr <- nrow(y)
+    na <- names(d)[x]
+    o <- data.frame(d = rep(na, nr), 
+                    COMID = d[[x]]$COMID, 
+                    toCOMID = d[[x]]$toCOMID)
+    names(o) <- c(n, "COMID", "toCOMID") 
+    o
+  }
+  
+  rpu <- do.call(rbind, lapply(1:length(out_rpu), make_df, d = out_rpu, n = "rpu"))
+  
+  vpu <- do.call(rbind, lapply(1:length(out_vpu), make_df, d = out_vpu, n = "vpu"))
+  
+  out_rpu_vpu <- left_join(rpu, vpu, by = "COMID") 
+  
+  out_rpu_vpu <- select(out_rpu_vpu, RPUID = rpu, VPUID = vpu, 
+                        COMID = COMID, toCOMID = toCOMID.x)
+  
+  out_rpu_vpu <- left_join(out_rpu_vpu, select(sf::st_drop_geometry(fline), COMID, toRPUID = RPUID, toVPUID = VPUID), by = c("toCOMID" = "COMID"))
+  
+  readr::write_csv(out_rpu_vpu, "cache/rpu_vpu_out.csv")
+}
+
+```
diff --git a/workspace/01_nwm_prep.Rmd b/workspace/01_nwm_prep.Rmd
index 21cd0396ba4c332ddaf96e03b8759e3e58beb1cc..c9787e24b6497049406e6500dddf4417c0bc1ee1 100644
--- a/workspace/01_nwm_prep.Rmd
+++ b/workspace/01_nwm_prep.Rmd
@@ -137,5 +137,6 @@ sf::write_sf(nhd, "cache/nwm_network.gpkg")
 # sbtools::authenticate_sb()
 # sbtools::item_replace_files("61295190d34e40dd9c06bcd7", "cache/nwm_network.gpkg")
 
+```
 
-```
\ No newline at end of file
+TODO: create nwm-appropriate flowline and catchment sets by VPU.
\ No newline at end of file
diff --git a/workspace/02_NHD_navigate.Rmd b/workspace/02_NHD_navigate.Rmd
index 4aed686b71b073ba3cf9ad332c02e862dc3952c9..cf67a7101d8d191173c0b493fba8371161a0cdb0 100644
--- a/workspace/02_NHD_navigate.Rmd
+++ b/workspace/02_NHD_navigate.Rmd
@@ -38,50 +38,17 @@ nwm_gages <- data.frame(
 
 RNetCDF::close.nc(nc)
 
-# Load NHDPlus Data if precprocessed to RDS format
-nhdplus_path(data_paths$nhdplus_gdb)
-staged_nhd <- stage_national_data(nhdplus_data = data_paths$nhdplus_gdb,
-                                  output_path  = data_paths$nhdplus_dir)
-
-nhd_rds <- fix_headwaters(staged_nhd$flowline, gsub("flowline.rds", "flowline_update.rds", staged_nhd$flowline),
-                          new_atts = data_paths$new_nhdp_atts, nhdpath = data_paths$nhdplus_gdb)
-
-nhdplusTools_data_dir(dir = data_paths$nhdplus_dir)
+# 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(nav_gpkg, nhd_flowline)){
-
-    nhd <- readRDS(nhd_rds) 
-    
-    # Subset NHD by VPU
-    nhd <- subset_vpu(nhd, vpu) 
-    
-    # there are some duplicates with enhd
-    nhd <- nhd %>%
-      group_by(COMID) %>%
-      filter(row_number() == 1) %>%
-      ungroup()
-
-    # Filter and write dendritic/non-coastal subset to gpkg
-    # This will be iterated over to produce the final network after POIs identified
-    non_dend <- unique(unlist(lapply(filter(nhd, TerminalFl == 1 & TotDASqKM < min_da_km)
-                                     %>% pull(COMID), NetworkNav, st_drop_geometry(nhd))))
-      
-    # Add fields to note dendritic and POI flowlines
-    nhd <- nhd %>% 
-      mutate(dend = ifelse(!COMID %in% non_dend, 1, 0),
-             poi = ifelse(!COMID %in% non_dend & TotDASqKM >= min_da_km, 1, 0)) 
-    
-    write_sf(nhd, nav_gpkg, nhd_flowline)
-  } else {
-     nhd <- read_sf(nav_gpkg, nhd_flowline)
-     try(nhd <- select(nhd, -c(minNet, WB, struct_POI, struct_Net, POI_ID)))
-  }
+   nhd <- read_sf(nav_gpkg, nhd_flowline)
+   try(nhd <- select(nhd, -c(minNet, WB, struct_POI, struct_Net, POI_ID)))
   
   # HUC02 includes some 
   if(vpu == "02"){
diff --git a/workspace/04_merge.Rmd b/workspace/04_merge.Rmd
index dfd0a77dc63c93b194e2dcf0dc94a0de8939d1cd..9d6c273fa8b2ab0da997055ace969eab3566599a 100644
--- a/workspace/04_merge.Rmd
+++ b/workspace/04_merge.Rmd
@@ -9,14 +9,9 @@ library(hyfabric)
 source("R/utils.R")
 source("R/hyRefactor_funs.R")
 
-# 1 - Feat type
-# 2 - Name of geopackage (""- Navigate; "Collapse", "ND")
 in_gpkg <- "reference_"
 
-gf_gpkg <- "cache/reference_Nav_CONUS.gpkg"
+gf_gpkg <- "cache/reference_CONUS.gpkg"
 
-Merge_VPU("WB", in_gpkg, gf_gpkg)
 Merge_VPU("final_POIs", in_gpkg, gf_gpkg)
-Merge_VPU("nsegment", in_gpkg, gf_gpkg)
-Merge_VPU("nhd_flowline", in_gpkg, gf_gpkg)
 ```
diff --git a/workspace/05_hyRefactor_AK.Rmd b/workspace/05_hyRefactor_AK.Rmd
index 5faf9a1f2d997b1ee375fbc61965174184a4b4aa..705a4973e991282241f0a0ab6ef94cc59606e047 100644
--- a/workspace/05_hyRefactor_AK.Rmd
+++ b/workspace/05_hyRefactor_AK.Rmd
@@ -146,7 +146,7 @@ lookup_table <- tibble::tibble(merit_COMID = unique(as.integer(lookup_table$memb
 
 
 readr::write_csv(lookup_table, lookup_table_file)
-write_sf(lookup_table, out_refac_gpkg, lookup_table_gpkg)
+write_sf(lookup_table, out_refac_gpkg, lookup_table_refactor)
 ```
 
 
@@ -200,7 +200,7 @@ If no events, the section skips to line 282
 
 ```{r POIs/outlets}
 # Read lookup for original COMID
-lookup_table <- read_sf(out_refac_gpkg, lookup_table_gpkg)
+lookup_table <- read_sf(out_refac_gpkg, lookup_table_refactor)
 
 # Create unified POI/outlet layer
 if(needs_layer(out_agg_gpkg, mapped_outlets_layer)) {
@@ -393,6 +393,6 @@ lookup_table <- tibble::tibble(merit_COMID = unique(as.integer(refactor_lookup$m
 
 if(!file.exists(lookup_table_file)) {
   readr::write_csv(lookup_table, lookup_table_file)
-  write_sf(lookup_table, out_refac_gpkg, lookup_table_gpkg)
+  write_sf(lookup_table, out_refac_gpkg, lookup_table_refactor)
 }
 ```
diff --git a/workspace/05_hyRefactor_HI.Rmd b/workspace/05_hyRefactor_HI.Rmd
index 199ead3ee7378f72c8162a9694e9b16832457e13..fd6c7ac2e39985488f190b1388ddb53afa0294a5 100644
--- a/workspace/05_hyRefactor_HI.Rmd
+++ b/workspace/05_hyRefactor_HI.Rmd
@@ -191,7 +191,7 @@ lookup_table <- tibble::tibble(NHDPlusV2_COMID = unique(as.integer(refactor_look
 
 if(!file.exists(lookup_table_file)) {
   readr::write_csv(lookup_table, lookup_table_file)
-  write_sf(lookup_table, out_refac_gpkg, lookup_table_gpkg)
+  write_sf(lookup_table, out_refac_gpkg, lookup_table_refactor)
 }
 ```
 
diff --git a/workspace/05_hyRefactor_flines.Rmd b/workspace/05_hyRefactor_flines.Rmd
index af5d0f69ba926464f0f159c780daa17187782d02..2bf60151c2b90e982cbe8496991f819876d6f5cd 100644
--- a/workspace/05_hyRefactor_flines.Rmd
+++ b/workspace/05_hyRefactor_flines.Rmd
@@ -26,54 +26,15 @@ See `hyRefactory_config.R` for all constants.
 source("R/utils.R")
 source("R/hyRefactor_funs.R")
 source("R/config.R")
-
-staged_nhd <- stage_national_data(nhdplus_data = data_paths$nhdplus_gdb,
-                                  output_path  = data_paths$nhdplus_dir)
-
-# Source modifications if iterating for Alaska
-if (vpu %in% c("19", "20")) source("R/ExCONUS_config.R")
 ```
 
 Load and filter source NHD Flowline network. 
 ```{r flowlines}
-if(needs_layer(out_refac_gpkg, nhd_flowline)) {
-  
-  # Determine which Flowlines are in scope, previously implemented in NHD_Navigate
-  nhd <- read_sf(nav_gpkg, nhd_flowline) 
-    
-  nhd <- nhd %>%
-    subset_rpu(rpu_code, run_make_standalone = TRUE) %>%
-    st_sf()
 
-  write_sf(nhd, out_refac_gpkg, nhd_flowline)
-} else {
   nhd <- read_sf(out_refac_gpkg, nhd_flowline)
-}
-```
 
-Read and filter catchments.
-```{r cats}
-if(needs_layer(out_refac_gpkg, nhd_catchment)) {
-  # Fullset of catchments from RPU
-  cats_rpu <- readRDS(data_paths$fullcats_table) %>%
-    filter(RPUID == rpu_code)
-  
-  # Full cats with geom
-  cat_rds <- readRDS(staged_nhd$catchment)
-  
-  sf::st_crs(cat_rds) <- sf::st_crs(cat_rds)
-  
-  # Attribute cat sets
-  cats <- cat_rds %>%
-    mutate(hy_cats = ifelse(FEATUREID %in% nhd$COMID, 1, 0),
-           full_cats = ifelse(FEATUREID %in% cats_rpu$FEATUREID, 1, 0)) %>%
-    filter(full_cats == 1 | hy_cats == 1) %>%
-    st_sf()
-
-  write_sf(cats, out_refac_gpkg, nhd_catchment)
-} else {
   cats <- read_sf(out_refac_gpkg, nhd_catchment)
-}
+
 ```
 
 
@@ -115,9 +76,6 @@ events <- read_sf("cache/gages_info.gpkg", "Gages") %>%
   filter(reach_meas - FromMeas > 10 & AreaSqKM > 2 & 
            ToMeas - reach_meas > 10 & LENGTHKM > (combine_meters / 1000)) %>%
   select(COMID = comid, REACHCODE = reachcode, REACH_meas = reach_meas, site_no) %>%
-  # Events cannot be in terminal POIs, 
-  #  code seems to ignore the command not to split/combine those flowlines
-  filter(!COMID %in% filter(outlets, type == "terminal")$COMID) %>%
   filter(!COMID %in% avoid$COMID)
 
 # Write out events and outelts
@@ -179,6 +137,6 @@ lookup_table <- tibble::tibble(NHDPlusV2_COMID = unique(as.integer(refactor_look
 
 
 readr::write_csv(lookup_table, lookup_table_file)
-write_sf(lookup_table, out_refac_gpkg, lookup_table_gpkg)
+write_sf(lookup_table, out_refac_gpkg, lookup_table_refactor)
 ```
 
diff --git a/workspace/06-2_aggregate_cats.Rmd b/workspace/06-2_aggregate_cats.Rmd
index d4c8e5d692b72980a8b467cfd754af7e81c7215f..d18ae6e360dc6724666c64a6d052b0e6c7d8bfbc 100644
--- a/workspace/06-2_aggregate_cats.Rmd
+++ b/workspace/06-2_aggregate_cats.Rmd
@@ -44,7 +44,7 @@ divides <- read_sf(out_refac_gpkg, divide_layer) %>%
   mutate(areasqkm = as.numeric(st_area(.)))
 
 # Read lookup for original COMID
-lookup <- read_sf(out_refac_gpkg, lookup_table_gpkg)
+lookup <- read_sf(out_refac_gpkg, lookup_table_refactor)
 
 # reconciled have no ID applied in refactoring flowlines
 reconciled <- st_transform(read_sf(out_refac_gpkg, reconciled_layer), crs)
@@ -249,7 +249,7 @@ lookup_table <- tibble::tibble(NHDPlusV2_COMID = unique(as.integer(refactor_look
     dplyr::left_join(aggregate_lookup_catchment, by = "reconciled_ID")
 
 readr::write_csv(lookup_table, lookup_table_file)
-write_sf(lookup_table, out_agg_gpkg, lookup_table_gpkg)
+write_sf(lookup_table, out_agg_gpkg, lookup_table_refactor)
 
 sf::st_layers(out_agg_gpkg)
 
diff --git a/workspace/07-1_merge.Rmd b/workspace/07-1_merge.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..7ba7e88da21712983b6b6dc87363222a2ba83d75
--- /dev/null
+++ b/workspace/07-1_merge.Rmd
@@ -0,0 +1,74 @@
+---
+title: "07-1_merge.Rmd"
+output: html_document
+editor_options:
+  chunk_output_type: console
+---
+title: 07-1_merge.Rmd
+Project: GFv2.0
+Date: 4/2022 
+Author: [David Blodgett](dblodgett@usgs.gov)
+Script purpose: Merge refactored outputs into VPU basis.
+
+          
+```{r rmd, echo=F, message=F, warning=FALSE, fig.width=4}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  fig.width=10, 
+  fig.height=8)
+```
+
+```{r setup}
+# Load custom functions
+source("R/config.R")
+source("R/utils.R")
+
+```
+
+For a given VPU bind RPUs into a single geopackage for an entire VPU. This step is intended to be run after Refactor flowlines/cats.  For a given VPU, this binds the results of aggregated from the associated RPUs into a single layer and writes them out to a single geopackage for the non-dendritic work. 
+
+```{r Bind RPUs}
+if(needs_layer(gf_gpkg, agg_cats_layer)) {
+  
+  # Thematic POIs
+  read_sf(nav_gpkg,  final_poi_layer) %>% 
+    st_transform(crs) %>%
+    write_sf(gf_gpkg, final_poi_layer)
+  
+  fline <- sf::read_sf(nav_gpkg, nhd_flowline)
+  
+  sf::write_sf(fline, gf_gpkg, nhd_flowline)
+
+  cats <- sf::read_sf(nav_gpkg, nhd_catchment)
+  
+  sf::write_sf(cats, gf_gpkg, nhd_catchment)
+  
+  merged_layers <- merge_refactor(rpu_codes$rpuid, rpu_vpu_out, 
+                                  lookup_table_refactor, 
+                                  reconciled_layer, 
+                                  divide_layer, 
+                                  agg_fline_layer,
+                                  agg_cats_layer,
+                                  mapped_outlets_layer)
+  
+  sf::write_sf(merged_layers[[lookup_table_refactor]], gf_gpkg, lookup_table_refactor)
+  
+  sf::write_sf(merged_layers[[reconciled_layer]], gf_gpkg, reconciled_layer)
+  
+  sf::write_sf(merged_layers[[divide_layer]], gf_gpkg, divide_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)
+  divides <- read_sf(gf_gpkg, divide_layer)
+  lookup <- read_sf(gf_gpkg, lookup_table_refactor)
+}
+```
diff --git a/workspace/07_NonDend.Rmd b/workspace/07-2_NonDend.Rmd
similarity index 61%
rename from workspace/07_NonDend.Rmd
rename to workspace/07-2_NonDend.Rmd
index f47bc0955ec4bb3deb26a10d3adbb8b7baf27409..0c2c2257bf78c49d30ec0d9b2351cfe34ceffc16 100644
--- a/workspace/07_NonDend.Rmd
+++ b/workspace/07-2_NonDend.Rmd
@@ -1,5 +1,5 @@
 ---
-title: "06_NonDend.Rmd"
+title: "07-2_NonDend.Rmd"
 output: html_document
 editor_options:
   chunk_output_type: console
@@ -22,43 +22,13 @@ knitr::opts_chunk$set(
 ```
 
 ```{r setup}
+# set to true to output extra layers.
+debug <- FALSE
+
 # Load custom functions
 source("R/config.R")
-# source("R/utils.R")
-# source("R/NHD_navigate.R")
-# source("R/hyRefactor_funs.R")
 source("R/non_dend.R")
 
-# Source modifications if iterating for Alaska
-if (vpu %in% c("19", "20")) source("R/ExCONUS_config.R")
-
-# FIXME: this shouldn't be in here... create a stand alone merge step?
-# Bind together layers for R03 and R10
-if (vpu %in% c("03", "10")){
- if(!file.exists(nav_gpkg)){
-   # filter geopackages to those of particular VPU
-   nav_gpkgs <- grep(paste0("reference_", vpu),
-                     list.files("cache", pattern = "gpkg", 
-                                full.names = TRUE), value = TRUE)
-   
-   # Binding function
-   featCON <- function(x, y){
-     vpu <- substr(basename(x), 11, 13)
-     print(vpu)
-     layername <- paste0(y, "_", vpu)
-     print(layername)
-     outlayer <- read_sf(x, layername)
-     print(dim(outlayer))
-     return(outlayer)
-     }
-   # Write out to new VPU gpkg
-   nsegment <- data.table::rbindlist(lapply(nav_gpkgs, featCON, "nsegment"), fill = TRUE)
-   write_sf(nsegment, nav_gpkg, paste0("nsegment_", vpu))
-   pois <- data.table::rbindlist(lapply(nav_gpkgs, featCON, "final_POIS"), fill = TRUE)
-   write_sf(pois, nav_gpkg, paste0("final_POIS_", vpu))
-  }
-}
-
 # Read in full VPU NHD set
 if (vpu == "20"){
    full_nhd <- readRDS(file.path(data_paths$islands_dir, "NHDPlusNationalData/nhdplus_flowline.rds")) %>%
@@ -77,6 +47,8 @@ if(vpu == "02"){
   grep_exp <- paste0("^", substr(vpu, start = 1, stop = 2))
 }
 
+cat_rpu_table <- readRDS(data_paths$fullcats_table)
+
 # get full flowline set including coastlines
 full_nhd <- full_nhd %>%
   filter(grepl(paste0("^", grep_exp, ".*"), .data$VPUID)) %>%
@@ -87,106 +59,11 @@ vpu_WBD <- readRDS(file.path(data_paths$nhdplus_dir, "HUC12.rds")) %>%
 
 elev <- data_paths$elev_cm[grepl(paste0("Ned", vpu), data_paths$elev_cm)]
 
-```
-
-For a given VPU bind RPUs into a single geopackage for an entire VPU in __ND_VPU.gpkg__. This step is intended to be run after Refactor flowlines/cats.  For a given VPU, this binds the results of aggregated from the associated RPUs into a single layer and writes them out to a single geopackage for the non-dendritic work. 
-
-Note that below I am combining the aggregated ID and the RPUID to create a distinct ID __*grID*__ for each aggregated feature within a VPU to avoid duplicate IDs at this stage.
+nhd <- st_transform(read_sf(gf_gpkg, nhd_flowline), crs)
+cats <- st_transform(read_sf(gf_gpkg, nhd_catchment), crs)
+divides <- st_transform(read_sf(gf_gpkg, divide_layer), crs)
+lookup <- read_sf(gf_gpkg, lookup_table_refactor)
 
-```{r Bind RPUs}
-if(needs_layer(ND_gpkg, agg_cats_layer)) {
-  
-  # FIXME: move this to output creation step?
-  # Thematic POIs
-  read_sf(nav_gpkg,  final_poi_layer) %>% 
-    st_transform(crs) %>%
-    write_sf(ND_gpkg, final_poi_layer)
-  
-  # FIXME: Should this be handled elsewhere? We have hyfabric::get_rpu_dependent_vars()
-  # Should we use something like that for lookup an rpu_codes?
-  # Maybe just move to a stand alone merge step?
-  # List of lookup files for this VPU
-  lookup <- list.files("cache", paste0("^", vpu, "+[a-z]+_lookup"))
-  
-  if(vpu == "03") lookup <- lookup[lookup != "03g_lookup.csv"]
-  if(vpu == "08") lookup <- c(lookup, "03g_lookup.csv")
-  
-  rpu_codes <- unlist(lapply(lookup, function(x) {substr(x, 1, 3)}))
-  
-  # FIXME move to a stand alone merge step?
-  # Combind layers from multiple RPUs 
-  layer_list <- c(nhd_flowline, nhd_catchment, reconciled_layer, divide_layer, 
-                  mapped_outlets_layer, split_layer)
-  
-  for (layer in layer_list) {
-    print(layer)
-    # List of NHDflowlines used to aggregate (note they are NOT the full flowlines)
-    do.call("rbind", lapply(rpu_codes, function(x) {
-      print(x)
-      if(layer == mapped_outlets_layer)
-        {read_sf(file.path("cache", paste0(x, "_aggregate.gpkg")), layer)}
-      else
-        {read_sf(file.path("cache", paste0(x, "_refactor.gpkg")), layer)}} %>%
-        # Bind together and create VPU-centric ID
-        mutate(rpu = x))) %>%
-        mutate(vpu_ID = row_number()) %>%
-        st_make_valid() %>%
-        st_transform(crs) %>%
-        write_sf(ND_gpkg, layer)
-  }
-  
-  # Read in to environment full files
-  cats <- read_sf(ND_gpkg, nhd_catchment)
-  divides <- read_sf(ND_gpkg, divide_layer)
-  nhd <- read_sf(ND_gpkg, nhd_flowline)
-  
-  # Lookup tables, the grID (group ID) field is essential to avoid duplicated 
-  #        aggregated IDs for VPUs with multiple RPUs
-  lookup <- do.call("rbind", lapply(lookup, function(x) 
-    read.csv(file.path("cache", x), header = T) %>%
-                                      mutate(rpu = substr(x, 1, 3)))) %>%
-    mutate(intCOM = as.integer(NHDPlusV2_COMID)) %>%
-    # Create new VPU-unique aggregated ID
-    group_by(intCOM) %>%
-    group_by(aggregated_catchment_ID, rpu) %>%
-    mutate(grID = cur_group_id(), intCOM = as.integer(NHDPlusV2_COMID)) %>%
-    ungroup() %>%
-    # Create new VPU-unique reconciled ID
-    group_by(reconciled_ID, rpu) %>%
-    mutate(recID = cur_group_id()) %>%
-    ungroup()
-  
-  write.csv(lookup, vpu_lookup_table_file, row.names = F)
-  
-  # Get distinct lookup of aggregated ID
-  lookup <- lookup %>%
-    distinct(aggregated_catchment_ID, rpu, grID)
-    
-  # Aggregated Flowlines, add group ID
-  do.call("rbind", lapply(rpu_codes, function(x) 
-    read_sf(file.path("cache", paste0(x, "_aggregate.gpkg")), agg_fline_layer) %>%
-      mutate(rpu = substr(x, 1, 3)))) %>%
-    inner_join(lookup, by = c("ID" = "aggregated_catchment_ID", "rpu" = "rpu")) %>%
-    write_sf(ND_gpkg, agg_fline_layer)
-  
-  read_agg_cats <- function(x) {
-    mutate(read_sf(file.path("cache", paste0(x, "_aggregate.gpkg")), 
-                   agg_cats_layer), 
-           rpu = substr(x, 1, 3))
-  }
-  
-  # Aggregated Catchments, add group ID
-  do.call("rbind", lapply(rpu_codes, read_agg_cats)) %>%
-    inner_join(lookup, by = c("ID" = "aggregated_catchment_ID", "rpu" = "rpu")) %>%
-    write_sf(ND_gpkg, agg_cats_layer)
-  
-} else {
-  # Load layers described above 
-  nhd <- read_sf(ND_gpkg, nhd_flowline)
-  cats <- read_sf(ND_gpkg, nhd_catchment)
-  divides <- read_sf(ND_gpkg, divide_layer)
-  lookup <- read.csv(file.path("cache", paste0("lookup_", vpu, ".csv")))
-}
 ```
 
 This next section intersects HUC12 for a given VPU with NHD+ catchments and catchment divides.
@@ -198,6 +75,9 @@ if(needs_layer(ND_gpkg, xwalk_layer)){
   # Read in full NHD cats
   cats <- filter(cats, full_cats == 1)
 
+  cats <- sf::st_make_valid(cats)
+  divides <- sf::st_make_valid(divides)
+  
   # Intersect NHD catchments and divides with HUC12
   nhd_wbd_int <- get_HUC12_Xwalk(vpu, cats, divides,
                                 file.path(data_paths$nhdplus_dir,
@@ -210,11 +90,11 @@ if(needs_layer(ND_gpkg, xwalk_layer)){
   
   divides <- divides %>%
     left_join(xwalk_divides_wbd %>%
-                group_by(vpu_ID) %>%
+                group_by(ID) %>%
                 filter(intArea == max(intArea)) %>%
                 ungroup() %>%
-                select(vpu_ID, HUC_12_int, intArea, AreaHUC12), 
-              by = "vpu_ID") 
+                select(ID, HUC_12_int, intArea, AreaHUC12), 
+              by = "ID") 
   
   # Bring over divides/HUC12 intersection information into divides layer
   xwalk_nhd_wbd <- st_drop_geometry(nhd_wbd_int$cats_HUC12) %>%
@@ -249,56 +129,60 @@ if(needs_layer(ND_gpkg, xwalk_layer)){
 
 This next step takes catchment divides and:
 
-1. Populates them with the unique per-VPU identifier (*__grID__*) from the lookup we added in the __bind RPUs__ chunk
-2. Adds a field that indicates how a divide without an aggregated catchment ID assigned in hyRefactor::aggregate_catchments was aggregated (*__aggStep__*)
-3. Identifies internal sinks and non-dendritic areas entirely surrounded by a single aggregated catchment, and assigns those the catchment ID (see __dissolve_holes__ in __non_dend.R__ function)
+1. Adds a field that indicates how a divide without an aggregated catchment ID assigned in hyRefactor::aggregate_catchments was aggregated (*__aggStep__*)
+2. Identifies internal sinks and non-dendritic areas entirely surrounded by a single aggregated catchment, and assigns those the catchment ID (see __dissolve_holes__ in __non_dend.R__ function)
   
 ```{r First-cut}
 if(needs_layer(ND_gpkg, divides_nd)){
-
+  
   # Bind divide with lookup, identify aggregation step
   divides_lu <- divides %>%
-    left_join(select(lookup, reconciled_ID, rpu, grID) %>%
-                distinct(), by = c("ID" = "reconciled_ID", "rpu")) %>%
-    rename(POI_ID = grID) %>%
+    left_join(distinct(select(lookup, reconciled_ID, aggregated_catchment_ID)), 
+              by = c("ID" = "reconciled_ID")) %>%
+    rename(POI_ID = aggregated_catchment_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, intCOM, POI_ID = aggregated_catchment_ID),
-              by = c("FEATUREID" =  "intCOM")) %>%
+    left_join(select(lookup, NHDPlusV2_COMID, POI_ID = aggregated_catchment_ID),
+              by = c("FEATUREID" =  "NHDPlusV2_COMID")) %>%
     filter(is.na(POI_ID)) %>%
-    select(member_COMID = FEATUREID, rpu, vpu_ID, HUC_12_int, intArea, AreaHUC12, POI_ID) %>%
-    mutate(ID = NA, vpu_ID = vpu, aggStep = NA) %>%
-    distinct()
+    select(FEATUREID, HUC_12_int, intArea, AreaHUC12, POI_ID) %>%
+    mutate(ID = NA, aggStep = NA) %>%
+    distinct() %>%
+    left_join(select(cat_rpu_table, FEATUREID, rpu = RPUID)) %>%
+    mutate(member_COMID = as.character(FEATUREID)) %>%
+    select(-FEATUREID)
 
   divides_lu <- divides_lu %>%
     select(-areasqkm) %>%
     rbind(cats_miss) %>%
     # resolve terminal non-POI outlets where we can
-    nhdplusTools:::check_valid() %>%
+    nhdplusTools:::check_valid()
   # clean_geometry(divides_lu, ID = "POI_ID", keep = 1.0, crs = st_crs(divides_lu))
-    dissolve_holes()
+  divides_lu <- dissolve_holes(divides_lu)
   
   rm(cats_miss)
 
   # DEBUG:
   # HRU layer
-  protoHRUs <- union_polygons_geos(divides_lu, "POI_ID") %>%
-    sfheaders::sf_remove_holes(.) %>%
-    st_make_valid()
-  write_sf(protoHRUs, ND_gpkg, HRU_layer)
-  rm(protoHRUs) 
-   
+  if(debug) {
+    protoHRUs <- union_polygons_geos(divides_lu, "POI_ID") %>%
+      sfheaders::sf_remove_holes(.) %>%
+      st_make_valid()
+    write_sf(protoHRUs, ND_gpkg, HRU_layer)
+    rm(protoHRUs) 
+  }
+  
   write_sf(divides_lu, ND_gpkg, divides_nd)
 } else {
   divides_lu <- read_sf(ND_gpkg, divides_nd)
 }
 ```
 
-This next section takes divides not assigned an aggregated ID (*__grID__*) and determines if they can be grouped on the basis of their HUC12/HUC10 ID.  See __assign_HUC10/12__ in __non_dend.R__. This assigns a temporary aggregated ID that reflects the HUC10/12 the divides overlap with.
+This next section takes divides not assigned an aggregated ID (*__POI_ID__*) and determines if they can be grouped on the basis of their HUC12/HUC10 ID.  See __assign_HUC10/12__ in __non_dend.R__. This assigns a temporary aggregated ID that reflects the HUC10/12 the divides overlap with.
 
 ```{r HUC12 cats}
 # Check for HUC12 aggregation assignments
@@ -315,13 +199,15 @@ if(!"HUC_12" %in% unique(divides_lu$aggStep)){
                              filter(full_nhd, FTYPE != "Coastline"), 0.66) %>%
     select(-comid_char)
   
-  # DEBUG:
-  # HRU layer
-  protoHRUs <- union_polygons_geos(divides_lu, "POI_ID") %>%
-    sfheaders::sf_remove_holes(.) %>%
-    st_make_valid()
-  write_sf(protoHRUs, ND_gpkg, HRU_layer)
-  rm(protoHRUs)
+  if(debug) {
+    # DEBUG:
+    # HRU layer
+    protoHRUs <- union_polygons_geos(divides_lu, "POI_ID") %>%
+      sfheaders::sf_remove_holes(.) %>%
+      st_make_valid()
+    write_sf(protoHRUs, ND_gpkg, HRU_layer)
+    rm(protoHRUs)
+  }
   
   # Update missing_cats
   write_sf(divides_lu, ND_gpkg, divides_nd)
@@ -347,14 +233,16 @@ if("Coastline" %in% unique(full_nhd$FTYPE)){
     # Function to create coastal catchments by HUC12 and catch inflows.
     divides_lu <- coastal_cats(divides_lu, full_nhd, vpu_WBD)
 
-    # DEBUG:
-    # HRU layer
-    protoHRUs <- union_polygons_geos(divides_lu, "POI_ID") %>%
-      sfheaders::sf_remove_holes(.) %>%
-      st_make_valid()
-    # Write out updates
-    write_sf(protoHRUs, ND_gpkg, HRU_layer)
-    rm(protoHRUs)
+    if(debug) {
+      # DEBUG:
+      # HRU layer
+      protoHRUs <- union_polygons_geos(divides_lu, "POI_ID") %>%
+        sfheaders::sf_remove_holes(.) %>%
+        st_make_valid()
+      # Write out updates
+      write_sf(protoHRUs, ND_gpkg, HRU_layer)
+      rm(protoHRUs)
+    }
     
     # Update missing_cats
     write_sf(divides_lu, ND_gpkg, divides_nd)
@@ -377,7 +265,7 @@ This is where alot of the more complex steps happen. For most VPUs, there are no
 ```{r FixingtheHoles}
 if(needs_layer(ND_gpkg,  "missing_cats")){
   print(paste0("Currently there are ", sum(is.na(divides_lu$POI_ID)), " cats without a POI_ID"))
-
+  
   # Handle nhd sinks
   # arrange divides_lu by member_COMID and populate with RowID for easy iteration within next two steps
   divides_lu <- divides_lu %>%
@@ -445,8 +333,11 @@ if(needs_layer(ND_gpkg,  "missing_cats")){
              aggStep = ifelse(!is.na(new_POI_ID), "boundary DEM", aggStep)) %>%
       select(-new_POI_ID)
     
-    write_sf(filter(divides_lu, is.na(POI_ID)), ND_gpkg, "missing_cats")
-    write_sf(term_outlets_sol$term_outlets, ND_gpkg, "unassigned_term")
+    write_sf(filter(divides_lu, is.na(POI_ID)), ND_gpkg, missing_cats)
+    
+    if(debug)
+      write_sf(term_outlets_sol$term_outlets, ND_gpkg, "unassigned_term")
+    
   } else {
     print ("all unaggregated catchments assigned")
     divides_lu <- term_outlets_sol$divides_lu
@@ -472,5 +363,3 @@ if(needs_layer(ND_gpkg,  "missing_cats")){
   all_hrus <- read_sf(ND_gpkg, HRU_layer)
 }
 ```
-
-
diff --git a/workspace/08_hyRefactor_HI.Rmd b/workspace/08_hyRefactor_HI.Rmd
deleted file mode 100644
index 796a719e800ac229c4fcd459ee4c05f6d8b26b98..0000000000000000000000000000000000000000
--- a/workspace/08_hyRefactor_HI.Rmd
+++ /dev/null
@@ -1,238 +0,0 @@
----
-title: "08_hyRefactor_HI"
-output: html_document
-editor_options:
-  chunk_output_type: console
----
-
-Project: GFv2.0  
-Script purpose: Refactor HI flowlines and catchments
-Author: [Andy Bock](abock@usgs.gov) 
-Author: [Dave Blodgett](dblodgett@usgs.gov)
-Notes: 
-
-```{r setup_0, echo=FALSE, cache=FALSE}
-knitr::opts_chunk$set(
-  collapse = TRUE,
-  comment = "#>",
-  fig.width=6, 
-  fig.height=4,
-  cache=FALSE)
-```
-
-Load functions and configuration. 
-See `hyRefactory_config.R` for all constants.
-```{r setup}
-vpu <- "20"
-
-# Hawaii RPUs - a, b, c, d, e ,f ,g, h
-if(!exists(rpu_code))(rpu_code <- "20a")
-
-source("R/utils.R")
-source("R/hyRefactor_funs.R")
-source("R/config.R")
-source("R/ExCONUS_config.R")
-```
-
-Load and filter source NHD Flowline network. 
-Load and filter source NHD Flowline network. 
-```{r flowlines}
-if(needs_layer(out_refac_gpkg, nhd_flowline)) {
-  
-  # Determine which Flowlines are in scope, previously implemented in NHD_Navigate
-  nhd <- read_sf(nav_gpkg, nhd_flowline) 
-    
-  # Ensure flowlines play nice within refactor
-  nhd <- nhd %>%
-    subset_rpu(rpu_code, run_make_standalone = TRUE) %>%
-    st_sf()
-
-  write_sf(nhd, out_refac_gpkg, nhd_flowline)
-} else {
-  nhd <- read_sf(out_refac_gpkg, nhd_flowline)
-}
-```
-
-Read and filter catchments.
-```{r cats}
-if(needs_layer(out_refac_gpkg, nhd_catchment)) {
-  # Fullset of catchments from RPU
-  cats_rpu <- readRDS(data_paths$fullcats_table) %>%
-    filter(RPUID == rpu_code)
-  
-  # Full cats with geom
-  cat_rds <- readRDS(staged_nhd$catchment)
-  
-  # Attribute cat sets for dendritic and full sets.
-  cats <- cat_rds %>%
-    mutate(hy_cats = ifelse(FEATUREID %in% nhd$COMID, 1, 0),
-           full_cats = ifelse(FEATUREID %in% cats_rpu$FEATUREID, 1, 0)) %>%
-    filter(full_cats == 1 | hy_cats == 1) %>%
-    st_sf()
-
-  write_sf(cats, out_refac_gpkg, nhd_catchment)
-} else {
-  cats <- read_sf(out_refac_gpkg, nhd_catchment)
-}
-```
-
-Add new POIs
-```{r new_POIs}
-# POIs, join with nhd to bring over attributes
-POIs <- read_sf(nav_gpkg, final_poi_layer) %>%
-  inner_join(select(st_drop_geometry(nhd), TotDASqKM, COMID, DnHydroseq), by = "COMID")
-
-# A minor modification to 20f
-if (rpu_code == "20f"){
-  new_outlet_POI_atts <- filter(st_drop_geometry(nhd), COMID == 800001119) %>%
-    select(COMID, Type_Conf_new = LevelPathI)
-  
-  POIs <- POIs %>%
-    left_join(new_outlet_POI_atts, by = "COMID") %>%
-    mutate(Type_Conf = ifelse(!is.na(Type_Conf_new), Type_Conf_new, Type_Conf)) %>%
-    select(-Type_Conf_new)
-}
-```
-
-Refactor flowlines, determine outlets and events
-
-```{r refactor}
-# 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))
-
-# Also need to avoid modification to flowlines immediately downstream of POIs
-#      This can cause some hydrologically-incorrect catchment aggregation
-POI_downstream <- filter(nhd, Hydroseq %in% POIs$DnHydroseq, AreaSqKM > 0)
-
-# build final outlets set
-outlets <- POIs %>%
-  mutate(type = ifelse(!is.na(Type_Term), "terminal", "outlet")) %>%
-  filter(COMID %in% cats$FEATUREID)
-
-# derive list of unique terminal paths
-TerminalPaths <- unique(filter(nhd, COMID %in% outlets$COMID)$TerminalPa)
-
-# attribute flowline network used for refactor
-nhd <- mutate(nhd, refactor = ifelse(TerminalPa %in% TerminalPaths, 1, 0))
-write_sf(nhd, out_refac_gpkg, nhd_flowline)
-
-# Need to avoid refactoring catchments that are long and thing and 
-# reasonably large. They cause problems in a downstream step.
-avoid <- dplyr::filter(nhd, (sqrt(AreaSqKM) / LENGTHKM) > 3 & AreaSqKM > 1)
-
-# Some minor customization to split thresholds due to HI gage customizations
-# Verified that it works on all RPUS
-events <- read_sf(nav_gpkg, gage_events_layer) %>%
-  #filter(Split == 1) %>%
-  inner_join(select(st_drop_geometry(filter(nhd, refactor == 1)), AreaSqKM, LENGTHKM, COMID, FromMeas, ToMeas), by = "COMID") %>%
-  #filter(AreaSqKM > 1 & LENGTHKM > (combine_meters / 250)) %>%
-  filter(site_no == "16265600" | REACH_meas - FromMeas > 5 & AreaSqKM > 0.5 & ToMeas - REACH_meas > 5) %>% 
-  select(site_no, COMID, REACHCODE, REACH_meas) %>%
-  # Events cannot be in terminal POIs, code seems to ignore the command not to split/combine those flowlines
-  filter(!COMID %in% filter(outlets, type == "terminal")$COMID) %>%
-  filter(!COMID %in% avoid$COMID)
-
-# write out events
-if(nrow(events) > 0) {write_sf(events, out_refac_gpkg, "split_events")}
-
-# write out outlets
-if(needs_layer(out_refac_gpkg, outlets_layer)) {
-  write_sf(outlets, out_refac_gpkg, outlets_layer)
-}
-
-# Prep flowlines for refactor
-nhdplus_flines <- sf::st_zm(filter(nhd, refactor == 1)) %>%
-  st_as_sf() 
-
-if(needs_layer(out_refac_gpkg, refactored_layer)) {
-  
-  tf <- paste0("temp/refactored_", rpu_code,".gpkg")
-  tr <- paste0("temp/reconciled_", rpu_code, ".gpkg")
-  # see hyRefactor_funs.R for these parameters.
-  # Original Refactor
-  refactor_nhdplus(nhdplus_flines = dplyr::select(nhdplus_flines, -FTYPE), # Select controls whether prepare_nhdplus is used. 
-                   split_flines_meters = split_meters, 
-                   split_flines_cores = 6,#para_split_flines, 
-                   collapse_flines_meters = combine_meters,
-                   collapse_flines_main_meters = combine_meters,
-                   out_refactored = tf, 
-                   out_reconciled = tr, 
-                   three_pass = TRUE, 
-                   purge_non_dendritic = FALSE, 
-                   events = events,
-                   exclude_cats = c(outlets$COMID, 
-                                    avoid$COMID, POI_downstream$COMID),
-                   warn = TRUE)
-  
-  write_sf(read_sf(tf), out_refac_gpkg, refactored_layer)
-  write_sf(read_sf(tr), out_refac_gpkg, reconciled_layer)
-  
-  unlink(list(tf, tr))
-  rm(tf, tr)
-}
-```
-
-Write NHD to refactored flowlines look up table
-
-```{r lookup}
-# create lookup for ref flowlines to use in the non-dendritic steps
-refactor_lookup <- dplyr::select(st_drop_geometry(read_sf(out_refac_gpkg, reconciled_layer)), ID, member_COMID) %>%
-  dplyr::mutate(member_COMID = strsplit(member_COMID, ",")) %>%
-  tidyr::unnest(cols = member_COMID) %>%
-  dplyr::mutate(NHDPlusV2_COMID = as.integer(member_COMID)) %>% # note as.integer truncates
-  dplyr::rename(reconciled_ID = ID)
-
-if(is.character(refactor_lookup$reconciled_ID)) refactor_lookup$reconciled_ID <- as.integer(refactor_lookup$reconciled_ID)
-
-lookup_table <- tibble::tibble(NHDPlusV2_COMID = unique(as.integer(refactor_lookup$member_COMID))) %>%
-    dplyr::left_join(refactor_lookup, by = "NHDPlusV2_COMID") 
-
-if(!file.exists(lookup_table_file)) {
-  readr::write_csv(lookup_table, lookup_table_file)
-  write_sf(lookup_table, out_refac_gpkg, lookup_table_gpkg)
-}
-```
-
-Reconcile catchments.  Note that aggregate catchments can be run with the 06-2 workbook.
-
-```{r reconcile_catchments}
-rpu_string <- paste0("rpu_", rpu_code)
-
-fac <- raster::raster(data_paths$fac[[rpu_string]])
-fdr <- raster::raster(data_paths$fdr[[rpu_string]])
-rast_crs <- raster::crs(fac)
-# reconciled have no ID applied in refactoring flowlines
-reconciled <- st_transform(read_sf(out_refac_gpkg, reconciled_layer), rast_crs) 
-# refactored has the original COMIDs and other attributes
-refactored <- st_transform(read_sf(out_refac_gpkg, refactored_layer), rast_crs) 
-nhd <- read_sf(out_refac_gpkg, nhd_flowline) #%>%
-  #filter(refactor == 1)
-cats <- st_transform(read_sf(out_refac_gpkg, nhd_catchment), rast_crs) 
-st_precision(cats) <- 30 # meters
-cats <- sf::st_simplify(cats, dTolerance = 0)
-cats <- sf::st_make_valid(cats)
-
-if(needs_layer(out_refac_gpkg, divide_layer)) {
-  
-  divides <- reconcile_catchment_divides(catchment = cats,
-                                         fline_ref = refactored,
-                                         fline_rec = reconciled,
-                                         fdr = fdr,
-                                         fac = fac,
-                                         para = para_reconcile, 
-                                         cache = cache_split,
-                                         keep = NULL,
-                                         vector_crs = crs) 
-  
-  if(exists("divides")) unlink(cache_split)
-  
-  write_sf(divides, out_refac_gpkg, divide_layer)
-} else {
-  divides <- read_sf(out_refac_gpkg, divide_layer) %>%
-    st_make_valid()
-}
-
-sf::st_layers(out_refac_gpkg)
-```
diff --git a/workspace/R/NHD_navigate.R b/workspace/R/NHD_navigate.R
index 149d9516f1b9d9f2490da110b182a3630188f880..55ff4a56aef08a7372ca2e653f361c9a4d02ce98 100644
--- a/workspace/R/NHD_navigate.R
+++ b/workspace/R/NHD_navigate.R
@@ -736,7 +736,7 @@ split_tt <- function(in_POI_ID, split_DF, tt_diff, max_DA){
     if(nrow() == 0){
       #print("done")
       # Done iterating on possible splits
-      return(split_tt_DF)}tt_pois_init_iter
+      return(split_tt_DF)}
     
     # Get the flowline POI
     tt_pois_init_pre <- tt_pois_init_iter %>%
diff --git a/workspace/R/config.R b/workspace/R/config.R
index 1dff406cba1207c9ee77b5f19f92255d1d9938e2..4a003f15d15c3a1c908b401cb1dcab1cfeb66372 100644
--- a/workspace/R/config.R
+++ b/workspace/R/config.R
@@ -7,7 +7,17 @@ library(hyfabric)
 library(mapview)
 library(terra)
 
+rpu_vpu <- readr::read_csv(system.file("extdata/rpu_vpu.csv", package = "hyfabric"), 
+                           show_col_types = FALSE)
+rpu_vpu_out <- readr::read_csv("cache/rpu_vpu_out.csv", 
+                               col_types = c("c", "c", "i" , "i"), show_col_types = FALSE)
+
 if(!exists("rpu_code")) {
+  if(exists("VPU")) {
+    vpu <- VPU
+    rpu_codes <- rpu_vpu[rpu_vpu$vpuid %in% vpu, ]
+    rpu_code <- rpu_codes$rpuid[1]
+  }
   vpu <- get_rpu_dependent_vars()
   rpu_code <- vpu$r
   vpu <- vpu$v
@@ -29,7 +39,7 @@ if(!nchar(cache_split <- Sys.getenv("CACHE_SPLIT")) > 5)
 if(!nchar(lookup_table_file <- Sys.getenv("LOOKUP_TABLE")) > 5) 
   lookup_table_file <- file.path("cache", paste0(rpu_code, "_lookup.csv"))
 
-if(!nchar(lookup_table_file <- Sys.getenv("VPU_LOOKUP_TABLE")) > 5) 
+if(!nchar(vpu_lookup_table_file <- Sys.getenv("VPU_LOOKUP_TABLE")) > 5) 
   vpu_lookup_table_file <- file.path("cache", paste0("lookup_", vpu, ".csv"))
 
 
@@ -39,6 +49,9 @@ options(scipen = 9999)
 crs <- 5070
 data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
 
+nhdplus_path(data_paths$nhdplus_gdb)
+suppressWarnings(staged_nhd <- stage_national_data())
+
 # Defined and used broadly
 nhd_flowline <- "nhd_flowline"
 nhd_catchment <- "nhd_catchment"
@@ -48,14 +61,14 @@ elev_diff <- 500 # Max difference in elevation within a single segment/POI
 max_elev_TT_DA <- 10 # Min drainage area to consider for elevation/travel_time breaks
 tt_diff <- 24 # Max number of hours between adjacent POIS
 nav_gpkg <- file.path("cache", paste0("reference_", vpu,".gpkg"))
-xwalk_layer <- paste0("HUC12_nhd_", vpu) # HUC12 - nhdcat crosswalk, built in Nav for VPU 20
+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
 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
 poi_xwalk_layer <- paste0("poi_xwalk_layer_", vpu) # POIs that changed COMIDS during the navigate part of the workflow
-final_poi_layer <- paste0("final_POIS_", vpu)
+final_poi_layer <- paste0("final_POIS")
 
 # Settings for refactor workflow
 split_meters <- 10000
@@ -79,16 +92,17 @@ agg_fline_layer <- "agg_fline"
 agg_cats_layer <- "agg_cats"
 outlets_layer <- "outlets"
 mapped_outlets_layer <- "mapped_outlets"
-lookup_table_gpkg <- paste0(rpu_code, "_lookup")
+lookup_table_refactor <- "lookup_table"
+
+gf_gpkg <- file.path("cache", paste0("GF_", vpu, ".gpkg"))
 
 # Defined during NonDend.Rmd
 ND_gpkg <- file.path("cache", paste0("ND_", vpu,".gpkg"))
-divides_xwalk <- paste0("divides_nhd_", vpu)
-HRU_layer <- paste0("poi_cats_", vpu)
-divides_nd <- paste0("divides_nd_", vpu)
-missing_terms <- paste0("miss_terms_", vpu)
-cat_boundary <- paste0("boundaries_", vpu)
-missing_cats <- paste0("miss_cats_", vpu)
+divides_xwalk <- paste0("divides_nhd")
+HRU_layer <- paste0("all_cats")
+divides_nd <- paste0("divides_nd")
+missing_terms <- paste0("miss_terms")
+missing_cats <- paste0("miss_cats")
 
 
 # only used in HI?
diff --git a/workspace/R/utils.R b/workspace/R/utils.R
index fa8a581e8eb75c83982da8688888507a2a444f98..8bef62d09cdc7f3a41655a20e83f799463027bb1 100644
--- a/workspace/R/utils.R
+++ b/workspace/R/utils.R
@@ -10,11 +10,11 @@ if(!(require(hyfabric, quietly = TRUE) && packageVersion("hyfabric") == hyfabric
                     repos = NULL, type = "source")
 }
 
-fix_headwaters <- function(nhd_rds, out_rds, new_atts = NULL, 
+fix_headwaters <- function(nhd_rds, out_gpkg, new_atts = NULL, 
                            nhdpath = nhdplusTools::nhdplus_path(), 
                            par = 6) {
   
-  if(!file.exists(out_rds)) {
+  if(!file.exists(out_gpkg)) {
   nhd <- readRDS(nhd_rds)
   
   # Need to replace headwater flowlines with the trimmed burn line event.
@@ -25,7 +25,7 @@ fix_headwaters <- function(nhd_rds, out_rds, new_atts = NULL,
                           ble, by = "COMID")
   
   ble <- sf::st_zm(sf::st_as_sf(ble))
-  fline <- sf::st_zm(nhd)
+  nhd <- sf::st_zm(nhd)
   
   sf::st_geometry(nhd)[nhd$StartFlag == 1] <- sf::st_geometry(ble)[nhd$StartFlag == 1]
   
@@ -38,7 +38,7 @@ fix_headwaters <- function(nhd_rds, out_rds, new_atts = NULL,
                                                         ArbolateSu = weight)
     
     nhd <- select(nhd, COMID, VPUID, RPUID, FTYPE, FromNode, ToNode, 
-                  StartFlag, StreamOrde, StreamCalc, Divergence, WBAREACOMI, 
+                  StartFlag, StreamCalc, Divergence, WBAREACOMI, 
                   VPUID, RPUID, DnMinorHyd) %>%
       right_join(new_atts, by = c("COMID" = "comid")) %>%
       nhdplusTools::align_nhdplus_names()
@@ -75,13 +75,10 @@ fix_headwaters <- function(nhd_rds, out_rds, new_atts = NULL,
   new_geom <- do.call(c, new_geom)
   st_geometry(nhd)[check] <- new_geom
   
-  saveRDS(nhd, out_rds)
+  sf::write_sf(nhd, out_gpkg, "reference_flowlines")
   
-  ### ONLY RUN IF CHANGED ###
-  # sbtools::authenticate_sb()
-  # sbtools::item_replace_files("5dcd5f96e4b069579760aedb", "data/NHDPlusNationalData/nhdplus_flowline_update.rds")
   }
-  return(out_rds)
+  return(out_gpkg)
 }
 
 ## Is available in sbtools now, but will leave till on cran.
@@ -273,3 +270,149 @@ cat_rpu <- function(fcats, nhd_gdb){
   lkup_rpu2 <- rbind(lkup_rpu, missrec_df)
   return(lkup_rpu2)
 }
+
+merge_refactor <- function(rpus, rpu_vpu_out, 
+                           lookup_table_refactor, 
+                           reconciled_layer,
+                           divide_layer, 
+                           agg_fline_layer,
+                           agg_cats_layer,
+                           mapped_outlets_layer) {
+  
+  
+  out <- rep(list(list()), length(rpus))
+  names(out) <- rpus
+  
+  s <- 10000000
+  
+  for(rpu in rpus) {
+
+    refactor_gpkg <- file.path("cache", paste0(rpu, "_refactor.gpkg"))
+    aggregate_gpkg <- file.path("cache", paste0(rpu, "_aggregate.gpkg"))
+    
+    out[[rpu]] <- setNames(list(sf::read_sf(aggregate_gpkg, lookup_table_refactor), 
+                                 sf::read_sf(refactor_gpkg, reconciled_layer), 
+                                 sf::read_sf(refactor_gpkg, divide_layer), 
+                                 sf::read_sf(aggregate_gpkg, agg_fline_layer), 
+                                 sf::read_sf(aggregate_gpkg, agg_cats_layer), 
+                                 sf::read_sf(aggregate_gpkg, mapped_outlets_layer)), 
+                            c(lookup_table_refactor, 
+                              reconciled_layer, 
+                              divide_layer, 
+                              agg_fline_layer, 
+                              agg_cats_layer, 
+                              mapped_outlets_layer))
+    
+    ### RECONCILED ###
+    
+    out[[rpu]][[reconciled_layer]] <- select(out[[rpu]][[reconciled_layer]], 
+                                              ID, toID, LENGTHKM, TotDASqKM, member_COMID)
+    
+    out[[rpu]][[reconciled_layer]]$newID <- 
+      seq(s, by = 1, length.out = nrow(out[[rpu]][[reconciled_layer]]))
+    
+    s <- max(out[[rpu]][[reconciled_layer]]$newID) + 1
+    
+    # get a toID that matches the new IDs
+    out[[rpu]][[reconciled_layer]] <- left_join(out[[rpu]][[reconciled_layer]], select(st_drop_geometry(out[[rpu]][[reconciled_layer]]), 
+                             ID, newtoID = newID), 
+                   by = c("toID" = "ID"))
+    
+    # find updates that we need to apply
+    update <- rpu_vpu_out[!rpu_vpu_out$toCOMID %in% out[[rpu]][[lookup_table_refactor]]$NHDPlusV2_COMID &
+                            rpu_vpu_out$COMID %in% out[[rpu]][[lookup_table_refactor]]$NHDPlusV2_COMID, ]
+    
+    # apply these updates for follow up if we got any
+    if(nrow(update) > 0) {
+      update <- left_join(update, out[[rpu]][[lookup_table_refactor]], by = c("COMID" = "NHDPlusV2_COMID")) %>%
+        left_join(select(out[[rpu]][[lookup_table_refactor]], NHDPlusV2_COMID, to_member_COMID = member_COMID), 
+                  by = c("toCOMID" = "NHDPlusV2_COMID"))
+      
+      out[[rpu]][[reconciled_layer]] <- left_join(out[[rpu]][[reconciled_layer]], select(update, reconciled_ID, toCOMID),
+                     by = c("ID" = "reconciled_ID"))
+    }
+    
+    #### LOOKUP TABLE ####
+    
+    update_id <- sf::st_drop_geometry(select(out[[rpu]][[reconciled_layer]], 
+                                             old_ID = ID, ID = newID))
+    
+    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)
+    
+    #### DIVIDE LAYER ####
+    
+    out[[rpu]][[divide_layer]] <- left_join(rename(out[[rpu]][[divide_layer]], old_ID = ID),
+                                            update_id, 
+                                            by = c("old_ID"))
+    
+    out[[rpu]][[divide_layer]] <- select(out[[rpu]][[divide_layer]], ID, areasqkm, member_COMID)
+    
+    out[[rpu]][[divide_layer]]$rpu <- rep(rpu, nrow(out[[rpu]][[divide_layer]]))
+    
+    #### aggregate layers ####
+    
+    for(l in c(agg_cats_layer, agg_fline_layer, mapped_outlets_layer)) {
+      
+      out[[rpu]][[l]] <- left_join(rename(out[[rpu]][[l]], old_ID = ID),
+                                   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 = ",")
+    }
+    
+  }
+  
+  out <- setNames(lapply(names(out[[1]]), function(x, out) {
+    dplyr::bind_rows(lapply(out, function(df, n) df[[n]], n = x))
+  }, out = out), names(out[[1]]))
+  
+  # blow up so we have unique COMIDs to join on.
+  long_form <- st_drop_geometry(out[[reconciled_layer]]) %>%
+    select(newID, member_COMID) %>%
+    mutate(member_COMID = strsplit(member_COMID, ",")) %>%
+    tidyr::unnest(cols = member_COMID) %>%
+    mutate(NHDPlusV2_COMID = as.integer(member_COMID)) %>%
+    select(-member_COMID, update_newtoID = newID)
+  
+  if("toCOMID" %in% names(out[[reconciled_layer]])) {
+    # NOTE: if the missing tocomid is in the next VPU they will still be NA
+    out[[reconciled_layer]] <- left_join(out[[reconciled_layer]], 
+                                         long_form, by = c("toCOMID" = "NHDPlusV2_COMID"))
+    out[[reconciled_layer]]$newtoID[!is.na(out[[reconciled_layer]]$toCOMID)] <- 
+      out[[reconciled_layer]]$update_newtoID[!is.na(out[[reconciled_layer]]$toCOMID)]
+  }
+  
+  out[[reconciled_layer]] <- select(out[[reconciled_layer]], ID = newID, 
+                                    toID = newtoID, LENGTHKM, TotDASqKM, 
+                                    member_COMID, refactor_ID = ID)
+  
+  out
+}
diff --git a/workspace/README.md b/workspace/README.md
index 93d16d57dced160e0b15c9340a1b95cf869a9e55..e52b4af95ab3e814ce7e4ec98b35638e4e9cfb8e 100644
--- a/workspace/README.md
+++ b/workspace/README.md
@@ -1,16 +1,103 @@
 # gfv2 workspace
 
-This is where you scripts and data live.
+This is where scripts and data live. 
 
-## 00_ workflow files
+- Rmd files: contain all workflow execution logic.
+- bin: folder for executable files (7z, etc).
+- cache: folder for cached process files.
+- data: folder for static downloaded data.
+- R: folder for R source files used in workflow.
+- reports: folder for process result reports.
+- runners: folder for orchestration files.
+- temp: folder for temporary data.
 
-These files are pre-processors that should not be run by all workflow users. They all generate files that are cached and downloaded elsewhere for general use.
+## cache files
 
-The exception to this is `00_get_data.Rmd` -- this should be run interactively to download required data.
+- data_paths.json: all data used by the workflow is referenced here. (checked in)
+- dropped_gages.csv: list of gages NOT included and why.
+- Gages_info.*: gages included with useful attribues.
+- enhd_nhdplusatts.csv: updated network attributes for NHDPlus ids.
+- enhd_atts_fline.gpkg: flowlines with updated network attributes for NHDPlus ids.
+- 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
+  - 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
+  - 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.
+  - 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
+  - divides_nhd
+  - HUC12_nhd
+  - miss_terms
+  - miss_cats
 
-## 01_ workflow files 
+The following summarizes the role of each workflow file.
 
-These files are pre-processors that work quickly enough that not all users must run them.
+## 00_ workflow files: preprocess caches
+
+- `00_get_data.Rmd` -- should be run interactively to download required data.
+- `00_AK_prep.Rmd` -- prepares AK domain data. All outputs are cached and can be downloaded with `00_get_data.Rmd`.
+- `00_enhd_nwm_network.Rmd` -- processes the core NHDPlus network along with updates from the enhd and nwm. All outputs are cached and can be downloaded with `00_get_data.Rmd`. 
+- `00_gage_info.Rmd` -- retrieves gage information. All outputs are cached and can be downloaded with `00_get_data.Rmd`. 
+
+## 01_ workflow files: preparation
+
+- `01_NHD_prep.Rmd` -- sets up basic data in VPU-based reference fabric files.
+- `01_Gage_selection.Rmd` -- refines gage list for inclusion in the rest of the workflow.
+- `01_HI_prep.Rmd` -- prepares the Hawaii domain for processing.
+- `01_nwm_prep.Rmd` -- specific processing for NWM network with handling for the extended domain NWM (non-NHDPlus) attributes.
+
+## 02_ workflow files: navigate / POI selection
+
+- `02_AK_navigate.Rmd`
+- `02_HI_navigate.Rmd` 
+- `02_NHD_navigate.Rmd`
+
+These files navigate the network determining a set of required points of interest, refining the points of interest for later use (snapping/merging, etc.), and generating initial output artifacts.
+
+## 04_ workflow files: merge navigate
+
+- `04_merge.Rmd` -- merges VPU-based processing from 02 step into a national summary.
+
+## 05_ workflow files: refactor
+
+- `05_hyrefactor_flines.Rmd` -- refactors NHDPlus network for CONUS domain, eliminating very small catchments and splitting very large catchments and those with points of interest that need to be enforced in the network. Runs on an RPU basis.
+- `05_hyrefactor_HI.Rmd` -- complete refactor and aggregate workflow for the HI domain.
+- `05_hyrefactor_AK.Rmd` -- complete refactor and aggregate workflow for the AK domain.
+
+## 06_ workflow files: catchment processing
+
+- `06-1_hyrefactor_cats.Rmd` -- processes catchments (splits and merges) according to 05 refactored flowlines.
+- `06-2_aggregate_cats.Rmd` -- aggregates network to minimal set of catchments for points of interest. Adds some additional points of interest to ensure network continuity.
+
+## 07_ workflow files: non-network handling
+
+- `07-1_merge.Rmd` -- Merges data from multiple RPUs into VPU-based files with internally consistent identifiers.
+- `07-2_NonDend.Rmd` -- post processing of catchment processing aiming to identify and fill in gaps in the catchment coverage along coast lines and in closed basins that aren't otherwise included.
 
 ## Running Markdown Programmatically
 
diff --git a/workspace/cache/data_paths.json b/workspace/cache/data_paths.json
index 21d1a451abfc25d567f108bf226ee518fd9a74d9..6dd6047ce7d81226d5c71a2b22d35c4b24a7221d 100644
--- a/workspace/cache/data_paths.json
+++ b/workspace/cache/data_paths.json
@@ -255,8 +255,8 @@
   "GFv11_gages_lyr": "data/GFv11/GFv11_gages.rds",
   "GFv11_gdb": "data/GFv11/GFv1.1.gdb",
   "gagesii_lyr": "data/SWIM_gage_loc/gagesII_9322_point_shapefile.shp",
-  "new_nhdp_rds": "data/NHDPlusNationalData/nhdplus_flowline_update.rds",
-  "resOPs_path" : "data/reservoir_data/ResOPsUS/attributes/reservoir_attributes.csv",
-  "istarf_url" : "data/reservoir_data/ISTARF-CONUS.csv",
-  "resops_NID_CW" : "data/reservoir_data/cw_ResOpsUS_NID.csv"
+  "ref_flowline": "data/NHDPlusNationalData/reference_flowline.gpkg",
+  "res_attributes": "data/reservoir_data/ResOpsUS/attributes/reservoir_attributes.csv",
+  "istarf": "data/reservoir_data/ISTARF-CONUS.csv",
+  "resops_NID_CW": "data/reservoir_data/cw_ResOpsUS_NID.csv"
 }
diff --git a/workspace/cache/rpu_vpu_out.csv b/workspace/cache/rpu_vpu_out.csv
new file mode 100644
index 0000000000000000000000000000000000000000..5c0f61a2ce495ba90217f22c64de594db6d2f90d
--- /dev/null
+++ b/workspace/cache/rpu_vpu_out.csv
@@ -0,0 +1,73 @@
+RPUID,VPUID,COMID,toCOMID,toRPUID,toVPUID
+14a,14,20734037,20734041,15b,15
+14a,14,18267749,20734041,15b,15
+14b,NA,1356530,1230231,14a,14
+14b,NA,10037848,16977472,14a,14
+14b,NA,4932636,4884933,14a,14
+15a,NA,20380457,21413723,15b,15
+16b,NA,10782860,11225406,16a,16
+04d,NA,13032045,13032047,04c,04
+02b,NA,4652274,4763680,02a,02
+05d,NA,19443509,3369818,05c,05
+05d,NA,19314608,3369818,05c,05
+05c,NA,3408307,10161930,05a,05
+05c,NA,1824414,10161930,05a,05
+05a,05,1844789,7469418,08b,08
+05b,NA,10384401,11866814,05a,05
+06a,06,1862004,1840191,05a,05
+06a,06,1861888,1840095,05a,05
+07c,NA,13208194,13208654,07b,07
+07c,NA,13208570,13208604,07b,07
+07c,NA,13208612,13208610,07b,07
+07c,NA,13208574,13208604,07b,07
+07c,NA,2456378,13208610,07b,07
+07c,NA,2463151,13208654,07b,07
+07b,NA,6966883,5802988,07a,07
+07b,NA,6966909,5803242,07a,07
+07b,NA,5007665,5802988,07a,07
+07a,NA,5802986,5002739,07b,07
+07a,07,5093446,7469418,08b,08
+03g,NA,933180378,938090363,08a,08
+03g,NA,15714825,22789074,08a,08
+03g,NA,22763604,938090362,08a,08
+03g,NA,22762964,167578936,08a,08
+08b,NA,17948698,19266232,08a,08
+08b,NA,18061868,19266232,08a,08
+08a,NA,167578937,933180292,03g,08
+10c,NA,940190279,19058746,10d,10L
+10c,NA,19058874,19058890,10d,10L
+10d,NA,7266957,18760000,10b,10L
+10d,NA,17219530,7310885,10b,10L
+10d,NA,18758588,18760000,10b,10L
+10d,NA,17416474,7310885,10b,10L
+10b,10L,6618232,6619612,07b,07
+10b,NA,2255489,4391417,10a,10L
+10b,NA,3727783,4391417,10a,10L
+10a,10L,6018266,3624763,07a,07
+10i,NA,12736390,12736376,10h,10U
+10i,NA,12736378,12736376,10h,10U
+10h,NA,3014997,21538028,10f,10U
+10g,NA,21873545,21538028,10f,10U
+10f,NA,19251225,11542462,10e,10U
+10f,NA,940130432,11543264,10e,10U
+10e,10U,11764402,17244390,10d,10L
+10e,10U,7227390,17244390,10d,10L
+11c,11,941140164,19408128,08a,08
+11d,NA,20973380,369129,11c,11
+11b,NA,512825,512821,11c,11
+11b,NA,254998,255644,11d,11
+11b,NA,512823,512821,11c,11
+11a,11,941140164,19408128,08a,08
+17b,NA,947020332,24233052,17c,17
+17b,NA,23475172,24233052,17c,17
+17a,NA,24192562,24219835,17b,17
+17d,NA,23894396,23894394,NA,17
+13a,NA,943030230,266815,13b,13
+13b,NA,285570,317247,13d,13
+13c,NA,334016,285370,13b,13
+03b,NA,11562080,8835948,03a,03N
+03c,NA,18254165,16660931,03d,03S
+12b,NA,1572640,1572612,12a,12
+12b,NA,1572614,1572612,12a,12
+12c,NA,1609042,3125214,12b,12
+12d,NA,9357548,3766322,12c,12
diff --git a/workspace/hyfabric_0.5.2.tar.gz b/workspace/hyfabric_0.5.2.tar.gz
deleted file mode 100644
index b75dcfc6ae8052544026771b1f728aea792e1419..0000000000000000000000000000000000000000
Binary files a/workspace/hyfabric_0.5.2.tar.gz and /dev/null differ
diff --git a/workspace/hyfabric_0.5.3.tar.gz b/workspace/hyfabric_0.5.3.tar.gz
new file mode 100644
index 0000000000000000000000000000000000000000..a6ea142351a796276f4ec5a793bf595e54b56b94
Binary files /dev/null and b/workspace/hyfabric_0.5.3.tar.gz differ
diff --git a/workspace/runners/run_one_R.R b/workspace/runners/run_one_R.R
index a3ed68c70f532ea6ccbbdb65b0429ffec9cce2c6..8d8d28256924e7dc0dfc4eb290df4cc7083a1487 100644
--- a/workspace/runners/run_one_R.R
+++ b/workspace/runners/run_one_R.R
@@ -2,24 +2,31 @@
 ## workflows. These should be run prior to checking in any changes to these 
 ## workflows.
 
-VPU <- "01"
-rpu_code <- "01a"
+VPU <- vpu_codes <- "14"
+
+source("R/config.R")
+source("R/utils.R")
 
 ref <- paste0("cache/reference_", VPU, ".gpkg")
 
 unlink(ref)
 
+for(r in rpu_codes$rpuid) {
+  unlink(paste0("cache/", r, "_refactor.gpkg"))
+}
+
+rmarkdown::render('01_NHD_prep.Rmd', 
+                  output_file='temp/01_NHD_prep_test.html')
+
 rmarkdown::render('02_NHD_navigate.Rmd', 
                   output_file='temp/02_NHD_navigate_test.html')
 
-rpus <- unique(sf::read_sf(ref, "nhd_flowline")$RPUID)
-rpus <- rpus[!is.na(rpus)]
-
-for(rpu_code in rpus) {
+for(rpu_code in rpu_codes$rpuid) {
   unlink(paste0("cache/", rpu_code, "_aggregate.gpkg"))
-  unlink(paste0("cache/", rpu_code, "_refactor.gpkg"))
   unlink(paste0("cache/", rpu_code, "_lookup.csv"))
   
+  source("R/config.R")
+  
   rmarkdown::render("05_hyRefactor_flines.Rmd",
                     output_file = "temp/05_hyRefactor_flines_test.html")
   
@@ -30,6 +37,10 @@ for(rpu_code in rpus) {
                     output_file = "temp/06-2_aggregate_cats_test.html")
 }
 
+unlink(paste0("cache/GF_", VPU, ".gpkg"))
+rmarkdown::render("07-1_merge.Rmd",
+                  output_file = "temp/07-1_merge_test.html")
+
 unlink(paste0("cache/ND_", VPU, ".gpkg"))
-rmarkdown::render("07_NonDend.Rmd",
-                  output_file = "temp/07_NonDend_test.html")
+rmarkdown::render("07-2_NonDend.Rmd",
+                  output_file = "temp/07-2_NonDend_test.html")