diff --git a/.gitignore b/.gitignore
index 3efeadd482ddba04b638393d0eac437d9f608755..b42611b0a09b2b69f2c3f43c3caec72bb054f5ce 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,3 +8,4 @@ workspace/rosm.cache/
 *.rds
 *.csv
 *.gpkg
+*.html
diff --git a/workspace/NHD_navigate.Rmd b/workspace/NHD_navigate.Rmd
index f26aa2e3c2246144e60852bc3b5e7e34d08b9526..eb24a044589927a4abdd193d0af96b75349692d7 100644
--- a/workspace/NHD_navigate.Rmd
+++ b/workspace/NHD_navigate.Rmd
@@ -1,151 +1,91 @@
-##################################################
-## Project: GFv2.0
-## Script purpose: Generate Segments and POIS from HUC12 outlets
-## Date: 12/4/2019
-## Author: abock@usgs.gov
-## Notes:  First cut, includes both HUC12 outlets as POIS and 
-##         connecting confluences; Run over Hydrologic Region 06
-##         in this example
-## NEEDED MODS
-#  0 - a) create setup script for downloading and formatting data 
-#      b) ensuring the proper libraries are installed and loaded.
-#  1 - Functionalize script with input arguments to cover workspace, results, etc.
-#  2 - Add Gages to workflow
-#  3 - Add NID/Waterbodies COMIDs to workflow
-#  3 - Collapse Points close together
-#  4 - Spit out segment length, travel time
-#  5 - identify incremental catchments associated with each downstream POI
-##################################################
-
-# Network navigation for upstream/downstream from a COMID of interest
-NetworkNav<-function(inCom,nhdDF,navType){
-  # Upstream network navigation is easier
-  Seg<-nhdplusTools::get_UM(nhdDF,inCom,include=T)
-  # Assign POI_ID to flowlines associated with each POI/Incremental Segment
-  if ("POI_ID" %in% colnames(nhdDF)){
-    nhdDF$POI_ID[nhdDF$COMID %in% Seg & nhdDF$POI_ID==0]<<-inCom
-  }
-  return(Seg)
-}
-
-# Network navigation to connect dangles in the network
-NetworkConnection<-function(inCom,nhdDF){
-  # Subset by dangling segments that need downstream navigation to be connected
-  #        reduces the number of features we need to connect
-  upNet_DF<-nhdDF %>% dplyr::filter(COMID %in% inCom) %>% 
-    filter(!DnHydroseq %in% Hydroseq)
-  
-  # while the number of dangles is greater than 0
-  while (length(upNet_DF$COMID)>0){
-    # create item for number of dangles
-    count=dim(upNet_DF)[1]
-    print (dim(upNet_DF))
-    # find out which level paths are downstream of dangling huc12 POIs
-    DSLP<-upNet_DF$DnLevelPat[upNet_DF$COMID %in% inCom]
-    # Get the COMID of the hydroseq with level path value
-    # the lowest downstream flowline within the levelpath
-    inCom2<-nhd$COMID[nhd$Hydroseq %in% DSLP]
-    # Run the upstream navigation code
-    upNet<-unique(unlist(lapply(inCom2,NetworkNav,nhdDF,"up")))
-    # Append result to existing segment list
-    inCom<-append(inCom,upNet)
-    # Get the same variable as above
-    upNet_DF<-nhdDF %>% dplyr::filter(COMID %in% inCom) %>% 
-      filter(!DnHydroseq %in% Hydroseq)
-    # Get the count
-    count2=dim(upNet_DF)[1]
-    # if the count has remained the same we are done and return the flowline list
-    if (count == count2){
-      return (inCom)
-    }
-  }
-  # Not sure this other return is needed
-  return(inCom)
-}
-
-# POI Creation
-POI_creation<-function(inSegs){
-  # break segs into points
-  Segpts <- st_geometry(inSegs) %>% 
-    lapply(., function(x) {sf::st_sfc(x) %>% sf::st_cast(., 'POINT')})
-  # subset the last point from each geometry, make a POINT sf object
-  Seg_ends <- sapply(Segpts, function(p) {
-    p[c(length(p))]}) %>% sf::st_sfc() %>% sf::st_sf('geom' = .)
-  # Add the COMID, otherwise comes out as feature with no relating attributes
-  Seg_ends$COMID<-inSegs$COMID
-  return(Seg_ends)
-}
-
-library(nhdplusTools) #https://cran.r-project.org/web/packages/nhdplusTools/nhdplusTools.pdf
+---
+title: "NHD Navigate"
+output: html_document
+---
+
+Project: GFv2.0  
+Script purpose: Generate Segments and POIS from HUC12 outlets  
+Date: 12/4/2019  
+Author: (Andy Bock)[mailto:abock@usgs.gov]  
+Notes:  First cut, includes both HUC12 outlets as POIS and 
+        connecting confluences; Run over Hydrologic Region 06
+        in this example
+
+NEEDED MODS
+1. a) create setup script for downloading and formatting data  
+   b) ensuring the proper libraries are installed and loaded.  
+2. Functionalize script with input arguments to cover workspace, results, etc.  
+3. Add Gages to workflow  
+4. Add NID/Waterbodies COMIDs to workflow  
+5. Collapse Points close together  
+6. Spit out segment length, travel time  
+7. identify incremental catchments associated with each downstream POI  
+
+```{r setup}
+# https://cran.r-project.org/web/packages/nhdplusTools/nhdplusTools.pdf
+library(nhdplusTools)
 library(tidyverse) # a collection of libraries that work with eacher other (dplyr)
 library(sf) # simple features - a library for handling spatial data
 library(rgdal) # rgdal - an open source GIS library that's more akin to ESRI and handles shapefile
 
+source("R/NHD_navigate.R")
+data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
+nhdplus_path(data_paths$nhdplus_gdb)
+staged_nhd <- stage_national_data(nhdplus_data = data_paths$nhdplus_gdb, 
+                                  output_path = data_paths$nhdplus_dir)
+
+huc12_pois_geo <- file.path("cache", "huc12pois.geojson")
+huc12_segs <- file.path("cache", "huc12segs.gpkg")
+conf_pois <- file.path("cache", "confpois.geojson")
+huc12segs_all <- file.path("cache", "huc12segs_all.gpkg")
+```
+
+```{r, eval = FALSE}
 #######################
-# Read Region 06 HUC12 outlets, get this data at https://www.sciencebase.gov/catalog/item/5dbc53d4e4b06957974eddae
-HUC12<-rgdal::readOGR("data/hu_points.gpkg","hu_points")
-# Add new region field
-HUC12@data$Region<-substr(HUC12@data$HUC12,1,2)
-# Get the COMIDS of HUC12 outlets within Region 06
-comIDs<-HUC12@data %>% dplyr::filter(Region == "06") %>% dplyr::select (COMID)
-
-# Set path to national NHDPlus v2 seamless data; 
-print ("This will take awhile, go take a walk")
-if (!file.exists("data/nhdplus_flowline_attributes.rds")){ 
-  # Download seamless data to this location
-  nhdplusTools::download_nhdplusv2("data")
-  # Set the path
-  nhdplusTools::nhdplus_path("data/NHDPlusNationalData/NHDPlusV21_National_Seamless_Flattened_Lower48.gdb")
-  # Stage the data as RDS for easy access
-  nhdplusTools::stage_national_data(include=c("attribute","flowline","catchment"),output_path="data") 
-}else{
-  print ("NHDPlusv2 data already downloaded and staged, reading in flowlines information; will do catchment work in another step")
-  # Read flowline attribute data frames for entire seamless database
-  nhdDF<-readRDS("data/nhdplus_flowline_attributes.rds") # %>% select(COMID, Pathlength, LENGTHKM, LevelPathI, UpHydroseq, Hydroseq, DnHydroseq)
-  # Read flowline simple features
-  nhd<-readRDS("data/nhdplus_flowline.rds") 
-}
+HUC12 <- read_sf(data_paths$hu12_points_path, "hu_points") %>%
+  filter(grepl("^06.*", .data$HUC12))
+
+nhdDF<-readRDS(staged_nhd$attributes) # %>% select(COMID, Pathlength, LENGTHKM, LevelPathI, UpHydroseq, Hydroseq, DnHydroseq)
+# Read flowline simple features
+nhd<-readRDS(staged_nhd$flowline) 
 
 # Return flowlinest that are HUC12 POIs
-huc12POIsegs<-nhd %>% dplyr::filter(COMID %in% comIDs[,1]) # 15 missing?
+huc12POIsegs<-nhd %>% filter(COMID %in% HUC12$COMID) # 15 missing?
 # Sort data frame by highest-> lowest LP, then subsequently highest->lowest HS within each LP
-huc12POIsegs<-huc12POIsegs %>% dplyr::arrange(desc(LevelPathI),desc(Hydroseq))
+huc12POIsegs<-huc12POIsegs %>% arrange(desc(LevelPathI),desc(Hydroseq))
 
 #########################
 # Create POIs from HUC12 segments
 huc12POIs<-POI_creation(huc12POIsegs)
 # Write out shapefile representing POIs for given theme
-sf::st_write(huc12POIs,dsn="workspace/huc12POIs.shp",delete_dsn = T)
-write_rds(huc12POIs,path="workspace/huc12POIS.rds")
+write_sf(huc12POIs, dsn = huc12_pois_geo, delete_dsn = TRUE)
 
 ##########################
 # Navigate upstream along the mainstems from HUC12 POIS
 
 # Navigate upstream from each POI
-upNet<-unique(unlist(lapply(comIDs[,1],NetworkNav,nhdDF,"up")))
+upNet<-unique(unlist(lapply(HUC12$COMID,NetworkNav,nhdDF,"up")))
 # Connect POIs that don't have connecting levelpath looking downstream and 
 #         produce unique set of flowlines linking POIs
 finalNet<-unique(NetworkConnection(upNet,nhdDF))
 
 # Subset nhd segs to navigation results and write to shapefile
 huc12Segs<-nhd %>% filter(COMID %in% finalNet)
-sf::st_write(huc12Segs,dsn="workspace/huc12segs_all.shp",delete_dsn = T)
-write_rds(huc12Segs,path="workspace/huc12segs_US.rds")
+write_sf(huc12Segs, dsn = huc12_segs)
 
 ##########################
 # Derive confluence POIs where they are absent
-nhdConf<-nhd %>% dplyr::filter(COMID %in% huc12Segs$COMID) # check on redundancy here
+nhdConf<-nhd %>% filter(COMID %in% huc12Segs$COMID) # check on redundancy here
 
 # find confluences 
-confluences<-nhdConf %>% dplyr::group_by(DnHydroseq) %>% dplyr::filter(n()>1) %>% 
-  dplyr::filter(!COMID %in% huc12POIsegs$COMID)
+confluences<-nhdConf %>% group_by(DnHydroseq) %>% filter(n()>1) %>% 
+  filter(!COMID %in% huc12POIsegs$COMID)
   
 # Create POIs from HUC12 segments
 confPOIs<-POI_creation(confluences)
 
 # Write out shapefile representing POIs for given theme
-sf::st_write(confPOIs,dsn="workspace/confPOIs.shp",delete_dsn = T)
-write_rds(confPOIs,path="workspace/confPOIs.rds")
+write_sf(confPOIs, dsn = conf_pois, delete_dsn = TRUE)
 
 #################################################################
 # Assign POI_ID to each flowline related to POI/Incremental catchment
@@ -158,16 +98,18 @@ nhdDF$POI_ID<-0
 # Combine POI Types
 allPOIs<-rbind(huc12POIs,confPOIs)
 # sort data frame by
-POIs_ordered<-nhdDF %>% filter(COMID %in% allPOIs$COMID) %>% arrange(desc(LevelPathI,Hydroseq))
+POIs_ordered<-nhdDF %>% filter(COMID %in% allPOIs$COMID) %>% 
+  arrange(desc(LevelPathI,Hydroseq))
 # Get flowlines (proto-segments) associated with each COMID
 upNet<-unique(unlist(lapply(POIs_ordered$COMID,NetworkNav,nhdDF,"up")))
 
 # Subset nhd segs to navigation results and write to shapefile
 allSegs<-nhd %>% filter(COMID %in% upNet)
 
-allSegs<-allSegs %>% inner_join(nhdDF, by=c("COMID","COMID")) %>% select(COMID,LevelPathI.x,Hydroseq.x,POI_ID,Shape)
-sf::st_write(allSegs,dsn="workspace/huc12segs_all.shp",delete_dsn = T)
-write_rds(allSegs,path="workspace/huc12segs_all.rds")
+allSegs<-allSegs %>% 
+  inner_join(nhdDF, by=c("COMID","COMID")) %>%
+  select(COMID,LevelPathI.x,Hydroseq.x,POI_ID,Shape)
+write_sf(allSegs, dsn = huc12segs_all, delete_dsn = TRUE)
 
 
 ##############DEPRECATED/In Development SECTION##################
@@ -187,7 +129,7 @@ write_rds(allSegs,path="workspace/huc12segs_all.rds")
 #     
 #     # Subset by dangling segments that need downstream navigation to be connected
 #     #        reduces the number of features we need to connect
-#     upNet_DF<-nhd %>% dplyr::filter(COMID %in% upNet) %>% 
+#     upNet_DF<-nhd %>% filter(COMID %in% upNet) %>% 
 #       filter(!DnHydroseq %in% Hydroseq)
 #     
 #   }else{
@@ -201,4 +143,5 @@ write_rds(allSegs,path="workspace/huc12segs_all.rds")
 #     #nhdDF$POI_ID[nhdDF$COMID %in% Seg & nhdDF$POI_ID==0]<<-comID
 #   }
 #   return(Seg)
-# }
\ No newline at end of file
+# }
+```
\ No newline at end of file
diff --git a/workspace/R/NHD_navigate.R b/workspace/R/NHD_navigate.R
new file mode 100644
index 0000000000000000000000000000000000000000..5e59f8f53fce32756c3bd4ab30a7854ac4da21a6
--- /dev/null
+++ b/workspace/R/NHD_navigate.R
@@ -0,0 +1,60 @@
+# Network navigation for upstream/downstream from a COMID of interest
+NetworkNav<-function(inCom,nhdDF,navType){
+  # Upstream network navigation is easier
+  Seg<-nhdplusTools::get_UM(nhdDF,inCom,include=T)
+  # Assign POI_ID to flowlines associated with each POI/Incremental Segment
+  if ("POI_ID" %in% colnames(nhdDF)){
+    nhdDF$POI_ID[nhdDF$COMID %in% Seg & nhdDF$POI_ID==0]<<-inCom
+  }
+  return(Seg)
+}
+
+# Network navigation to connect dangles in the network
+NetworkConnection<-function(inCom,nhdDF){
+  # Subset by dangling segments that need downstream navigation to be connected
+  #        reduces the number of features we need to connect
+  upNet_DF<-nhdDF %>% dplyr::filter(COMID %in% inCom) %>% 
+    filter(!DnHydroseq %in% Hydroseq)
+  
+  # while the number of dangles is greater than 0
+  while (length(upNet_DF$COMID)>0){
+    # create item for number of dangles
+    count=dim(upNet_DF)[1]
+    print (dim(upNet_DF))
+    # find out which level paths are downstream of dangling huc12 POIs
+    DSLP<-upNet_DF$DnLevelPat[upNet_DF$COMID %in% inCom]
+    # Get the COMID of the hydroseq with level path value
+    # the lowest downstream flowline within the levelpath
+    inCom2<-nhd$COMID[nhd$Hydroseq %in% DSLP]
+    # Run the upstream navigation code
+    upNet<-unique(unlist(lapply(inCom2,NetworkNav,nhdDF,"up")))
+    # Append result to existing segment list
+    inCom<-append(inCom,upNet)
+    # Get the same variable as above
+    upNet_DF<-nhdDF %>% dplyr::filter(COMID %in% inCom) %>% 
+      filter(!DnHydroseq %in% Hydroseq)
+    # Get the count
+    count2=dim(upNet_DF)[1]
+    # if the count has remained the same we are done and return the flowline list
+    if (count == count2){
+      return (inCom)
+    }
+  }
+  # Not sure this other return is needed
+  return(inCom)
+}
+
+# POI Creation
+POI_creation<-function(inSegs){
+  # break segs into points
+  Segpts <- st_geometry(inSegs) %>% 
+    lapply(., function(x) {sf::st_sfc(x) %>% sf::st_cast(., 'POINT')})
+  # subset the last point from each geometry, make a POINT sf object
+  Seg_ends <- sapply(Segpts, function(p) {
+    p[c(length(p))]}) %>% 
+    sf::st_sfc(crs = st_crs(inSegs)) %>% 
+    sf::st_sf('geom' = .)
+  # Add the COMID, otherwise comes out as feature with no relating attributes
+  Seg_ends$COMID<-inSegs$COMID
+  return(Seg_ends)
+}
\ No newline at end of file
diff --git a/workspace/cache/data_paths.json b/workspace/cache/data_paths.json
index acb6d8adcedb231268cd02968962a4b966a5d8a7..e33bcab6042253c0458a8a33b7095b69bcbaf523 100644
--- a/workspace/cache/data_paths.json
+++ b/workspace/cache/data_paths.json
@@ -2,5 +2,6 @@
   "data_dir": "data",
   "hu12_points_path": "data/hu_points.gpkg",
   "gagesiii_points_path": "data/GAGESIII_gages",
-  "nhdplus_dir": "data/NHDPlusNationalData"
+  "nhdplus_dir": "data/NHDPlusNationalData",
+  "nhdplus_gdb": "data/NHDPlusNationalData/NHDPlusV21_National_Seamless.gdb"
 }
diff --git a/workspace/get_data.Rmd b/workspace/get_data.Rmd
index 7b27f5e52ad1db86eb8bbd17939d99ae78863a84..3ceba36aba6274a018eec0d6e8aa17a3ad669d75 100644
--- a/workspace/get_data.Rmd
+++ b/workspace/get_data.Rmd
@@ -58,19 +58,21 @@ out_list <- c(out_list, list(gagesiii_points_path = gagesiii_points_path))
 
 ```{r}
 nhdplus_dir <- file.path(data_dir, "NHDPlusNationalData")
-if(!file.exists(nhdplus_dir)) {
-  nhdplus_gdb <- download_nhdplusv2(data_dir)
 
-  nhdplus_path(nhdplus_gdb)
+if(!file.exists(nhdplus_dir)) message("downloading NHDPlus...")
 
-  staged_nhdplus <- stage_national_data()
-}
+suppressMessages(nhdplus_gdb <- download_nhdplusv2(data_dir))
+
+nhdplus_path(nhdplus_gdb)
+
+suppressWarnings(staged_nhdplus <- stage_national_data())
 
-out_list <- c(out_list, list(nhdplus_dir = nhdplus_dir))
+out_list <- c(out_list, list(nhdplus_dir = nhdplus_dir, nhdplus_gdb = nhdplus_gdb))
 ```
 
 Don't edit the following chunk.
 
 ```{r}
 write_json(out_list, path = out_file, pretty = TRUE, auto_unbox = TRUE)
+rm(out_list)
 ```