diff --git a/workspace/02_NHD_navigate.Rmd b/workspace/02_NHD_navigate.Rmd
index faa5ab54503095e0654b2efca622085fd395579f..238b5aafd99e3961530901a43e5673fe0c0f697c 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
 if(needs_layer(temp_gpkg, nav_poi_layer)) {
   
    nhd <- read_sf(nav_gpkg, nhd_flowline)
-   try(nhd <- select(nhd, -c(minNet, WB, struct_POI, struct_Net, POI_ID)))
+   try(nhd <- select(nhd, -c(minNet, WB, struct_POI, struct_Net, POI_ID)), silent = TRUE)
   
   # HUC02 includes some 
   if(vpu == "02"){
@@ -237,7 +237,7 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
   # Get NHDArea and HI-Res waterbodies for ResOpsUS locaitons
   WBs_VPU_areaHR <- HR_Area_coms(nhd, WBs_VPU, data_paths)
   # If the return is not NA
-  if(!is.na(WBs_VPU_areaHR)){
+  if(inherits(WBs_VPU_areaHR, "list")){
     # Attribute flowlines that intersect NHDARea and/or NHD HR waterbodies with 
     #           waterbody COMIDs
     nhd <- nhd %>%
@@ -264,10 +264,14 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
                       select(wb_outlet_events, COMID, WBAREACOMI))
     }
     
+    WBs_VPU_areaHR$wb <- nhdplusTools::st_compatibalize(
+      WBs_VPU_areaHR$wb, WBs_VPU)
+    
     # Append NHDArea and NHDHR waterbodies to waterbody layer
     WBs_VPU <- data.table::rbindlist(list(WBs_VPU, WBs_VPU_areaHR$wb), fill = TRUE) %>%
-      st_as_sf() %>%
-      st_make_valid()
+      st_as_sf()
+    
+    st_geometry(WBs_VPU) <- st_make_valid(st_geometry(WBs_VPU))
   }
     
   # Write out waterbodies
@@ -288,9 +292,9 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
   }
     
   write_sf(nhd, nav_gpkg, nhd_flowline)
-  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
   nhd <- read_sf(nav_gpkg, nhd_flowline)
   WBs_VPU <- read_sf(nav_gpkg, WBs_layer)
   resops_POIs_df <- read.csv(file.path("data/reservoir_Data", paste0("resops_POIs_",vpu,".csv")))
@@ -637,6 +641,7 @@ if(needs_layer(temp_gpkg, nsegments_layer)) {
     mutate(struct_POI = ifelse(COMID %in% strucFeat$struc_POIs$COMID, 1, 0),
            struct_Net = ifelse(COMID %in% strucFeat$structnet$COMID, 1, 0))
 
+  write_sf(nhd_Final, nav_gpkg, nhd_flowline)
   write_sf(nsegments, temp_gpkg, nsegments_layer)
 } else {
   # Read in NHDPlusV2 flowline simple features and filter by vector processing unit (VPU)