split nav_gpkg into ref and temp, add "snapped" attribute to POIs

workspace/02_NHD_navigate.Rmd

@@ -41,13 +41,14 @@ RNetCDF::close.nc(nc)
 # need some extra attributes
 atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.rds"))
 # atts <- select(atts, COMID, MAXELEVSMO, MINELEVSMO, VA_MA, TOTMA, WBAreaType)
+
 ```
 
 ```{r huc12 POIs}
 #  Derive or load HUC12 POIs
-if(needs_layer(nav_gpkg, nav_poi_layer)) {
+if(needs_layer(temp_gpkg, nav_poi_layer)) {
   
-   nhd <- read_sf(nav_gpkg, nhd_flowline)
+   nhd <- read_sf(ref_gpkg, nhd_flowline)
    try(nhd <- select(nhd, -c(minNet, WB, struct_POI, struct_Net, POI_ID)))
   
   # HUC02 includes some 
@@ -72,12 +73,12 @@ if(needs_layer(nav_gpkg, nav_poi_layer)) {
   huc12_POIs <- POI_creation(st_drop_geometry(HUC12_COMIDs), filter(nhd, poi == 1), "HUC12")
 
   # Write out geopackage layer representing POIs for given theme
-  write_sf(huc12_POIs, nav_gpkg, nav_poi_layer)
+  write_sf(huc12_POIs, temp_gpkg, nav_poi_layer)
   tmp_POIs <- huc12_POIs
 } else {
   # Load HUC12 POIs as the tmpPOIs if they already exist
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer) 
-  nhd <- read_sf(nav_gpkg, nhd_flowline)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer) 
+  nhd <- read_sf(ref_gpkg, nhd_flowline)
 }
 
 mapview(filter(tmp_POIs, Type_HUC12 != ""), layer.name = "HUC12 POIs", col.regions = "red") 
@@ -107,7 +108,7 @@ if(all(is.na(tmp_POIs$Type_Gages))) {
   # As a fail-safe, write out list of gages not assigned a POI
   if(nrow(filter(streamgages_VPU, !site_no %in% tmp_POIs$Type_Gages)) > 0) {
     write_sf(filter(streamgages_VPU, !site_no %in% tmp_POIs$Type_Gages),
-             nav_gpkg, "unassigned_gages")
+             temp_gpkg, "unassigned_gages")
     
     # Start documenting gages that are dropped out; these gages have no mean daily Q
     gage_document(filter(streamgages_VPU, !site_no %in% tmp_POIs$Type_Gages) %>%
@@ -117,9 +118,9 @@ if(all(is.na(tmp_POIs$Type_Gages))) {
   }
   
   # Write out geopackage layer representing POIs for given theme
-  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
 }
 
 mapview(filter(tmp_POIs, !is.na(Type_Gages)), layer.name = "Streamgage POIs", col.regions = "blue") 
@@ -151,14 +152,14 @@ if(all(is.na(tmp_POIs$Type_TE))) {
   # As a fail-safe, write out list of TE plants not assigned a POI
   if(nrow(filter(TE_COMIDs, !COMID %in% tmp_POIs$COMID)) > 0) {
     write_sf(filter(TE_COMIDs, !COMID %in% tmp_POIs$COMID),
-             nav_gpkg, "unassigned_TE")
+             temp_gpkg, "unassigned_TE")
   }
   
   # Write out geopackage layer representing POIs for given theme
-  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
   # Load TE POIs if they already exist
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
 }
 
 mapview(filter(tmp_POIs, Type_TE != ""), layer.name = "TE Plant POIs", col.regions = "blue") 
@@ -182,10 +183,10 @@ if(all(is.na(tmp_POIs$Type_Term))) {
   tmp_POIs <- POI_creation(terminal_POIs, nhd, "Term") %>%
     addType(., tmp_POIs, "Term")
 
-  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
-  nhd <- read_sf(nav_gpkg, nhd_flowline)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+  nhd <- read_sf(ref_gpkg, nhd_flowline)
 }
 
 mapview(filter(tmp_POIs, !is.na(Type_Term)), layer.name = "Terminal POIs", col.regions = "blue") 
@@ -218,11 +219,10 @@ if(all(is.na(tmp_POIs$Type_Conf))) {
   tmp_POIs <- POI_creation(conf_COMIDs, filter(nhd, minNet == 1), "Conf") %>%
     addType(., tmp_POIs, "Conf")
   
-  write_sf(nhd, nav_gpkg, nhd_flowline)
-  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
-  nhd <- read_sf(nav_gpkg, nhd_flowline)
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
+  nhd <- read_sf(ref_gpkg, nhd_flowline)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
 }
 
 mapview(filter(tmp_POIs, !is.na(Type_Conf)), layer.name = "Confluence POIs", col.regions = "blue") 
@@ -245,10 +245,10 @@ if(all(is.na(tmp_POIs$Type_NID))) {
     addType(., tmp_POIs, "NID", bind = FALSE)
 
   # Write out geopackage layer representing POIs for given theme
-  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
   # Load NID POIs if they already exist
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
 }
 
 mapview(filter(tmp_POIs, Type_NID != ""), layer.name = "NID POIs", col.regions = "blue") 
@@ -292,7 +292,7 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
     dplyr::filter(COMID %in% wb_lst$COMID)
   
   # Write out waterbodies
-  write_sf(st_transform(WBs_VPU, 4269), nav_gpkg, WBs_layer)
+  write_sf(st_transform(WBs_VPU, sf::st_crs(nhd)), ref_gpkg, WBs_layer)
 
   # Segments that are in waterbodies
   nhd <- mutate(nhd, WB = ifelse(WBAREACOMI > 0 & WBAREACOMI %in% WBs_VPU$COMID, 1, 0))
@@ -335,11 +335,10 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
     write.csv(resops_POIs_df, file.path("data/reservoir_Data", paste0("resops_POIs_",vpu,".csv")))
   }
     
-  write_sf(nhd, nav_gpkg, nhd_flowline)
-  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
-  nhd <- read_sf(nav_gpkg, nhd_flowline)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+  nhd <- read_sf(ref_gpkg, nhd_flowline)
 }
 
 mapview(filter(tmp_POIs, !is.na(Type_WBIn)), layer.name = "Waterbody inlet POIs", col.regions = "blue") +
@@ -404,10 +403,10 @@ if(nrow(elev_pois_split) > 0) {
   
   tmp_POIs$Type_Elev <- rep(NA, nrow(tmp_POIs))
 
-  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
     
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer) 
-  nhd <- read_sf(nav_gpkg, nhd_flowline)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer) 
+  nhd <- read_sf(ref_gpkg, nhd_flowline)
 }
 
 if(!all(is.na(tmp_POIs$Type_Elev)))
@@ -457,10 +456,10 @@ if(all(is.na(tmp_POIs$Type_Travel))) {
   tmp_POIs <- POI_creation(select(tt_POIs, COMID, Type_Travel), nhd, "Travel") %>%
     addType(., tmp_POIs, "Travel")
   
-  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
+  write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 }else{
-  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer)
-  nhd <- read_sf(nav_gpkg, nhd_flowline)
+  tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+  nhd <- read_sf(ref_gpkg, nhd_flowline)
 }
 
 mapview(filter(tmp_POIs, Type_Travel == 1), layer.name = "Travel time break POIs", col.regions = "blue") 
@@ -468,7 +467,7 @@ mapview(filter(tmp_POIs, Type_Travel == 1), layer.name = "Travel time break POIs
   
 ```{r Final POIs}
 # Derive final POI set ----------------------
-if(needs_layer(nav_gpkg, pois_all_layer)) {
+if(needs_layer(temp_gpkg, pois_all_layer)) {
 
   unCon_POIs <- filter(st_drop_geometry(filter(nhd, minNet == 1)), COMID %in% tmp_POIs$COMID, AreaSqKM == 0)
 
@@ -484,29 +483,29 @@ if(needs_layer(nav_gpkg, pois_all_layer)) {
     # bind together
     final_POIs <- bind_rows(tmp_POIs_fixed, select(poi_fix$new_POIs, -c(oldPOI, to_com)))
     # Write out fixes
-    write_sf(new_POIs, nav_gpkg, poi_moved_layer)
-    write_sf(xWalk, nav_gpkg, poi_xwalk_layer)
+    write_sf(new_POIs, temp_gpkg, poi_moved_layer)
+    write_sf(xWalk, temp_gpkg, poi_xwalk_layer)
     # write out final POIs
-    write_sf(final_POIs, nav_gpkg, pois_all_layer)
+    write_sf(final_POIs, temp_gpkg, pois_all_layer)
   } else {
     # If no fixes designate as NA
     poi_fix <- NA
     # write out final POIs
-    write_sf(tmp_POIs, nav_gpkg, pois_all_layer)
+    write_sf(tmp_POIs, temp_gpkg, pois_all_layer)
 
   }
 
 } else {
   # Need all three sets for attribution below
-  final_POIs <- read_sf(nav_gpkg, pois_all_layer)
-  new_POIs <- if(!needs_layer(nav_gpkg, poi_moved_layer)) read_sf(nav_gpkg, poi_moved_layer) else (NA)
-  xWalk <- if(!needs_layer(nav_gpkg, poi_xwalk_layer)) read_sf(nav_gpkg, poi_xwalk_layer) else (NA)
+  final_POIs <- read_sf(temp_gpkg, pois_all_layer)
+  new_POIs <- if(!needs_layer(temp_gpkg, poi_moved_layer)) read_sf(temp_gpkg, poi_moved_layer) else (NA)
+  xWalk <- if(!needs_layer(temp_gpkg, poi_xwalk_layer)) read_sf(temp_gpkg, poi_xwalk_layer) else (NA)
   unCon_POIs <- filter(st_drop_geometry(filter(nhd, minNet == 1)), COMID %in% tmp_POIs$COMID, AreaSqKM == 0)
 }
 ``` 
   
 ```{r Draft nsegment generation}
-if(needs_layer(nav_gpkg, nsegments_layer)) {
+if(needs_layer(temp_gpkg, nsegments_layer)) {
 
   # Sort POIs by Levelpath and Hydrosequence in upstream to downstream order
   seg_POIs <-  filter(st_drop_geometry(nhd), COMID %in% tmp_POIs$COMID, COMID %in% filter(nhd, minNet == 1)$COMID)
@@ -528,13 +527,12 @@ if(needs_layer(nav_gpkg, nsegments_layer)) {
     mutate(struct_POI = ifelse(COMID %in% strucFeat$struc_POIs$COMID, 1, 0),
            struct_Net = ifelse(COMID %in% strucFeat$structnet$COMID, 1, 0))
 
-  write_sf(nhd_Final, nav_gpkg, nhd_flowline)
-  write_sf(nsegments, nav_gpkg, nsegments_layer)
+  write_sf(nsegments, temp_gpkg, nsegments_layer)
 } else {
   # Read in NHDPlusV2 flowline simple features and filter by vector processing unit (VPU)
-  final_POIs <- read_sf(nav_gpkg, pois_all_layer)
-  nhd_Final <- read_sf(nav_gpkg, nhd_flowline)
-  nsegments <- read_sf(nav_gpkg, nsegments_layer)
+  final_POIs <- read_sf(temp_gpkg, pois_all_layer)
+  nhd_Final <- read_sf(ref_gpkg, nhd_flowline)
+  nsegments <- read_sf(temp_gpkg, nsegments_layer)
 }
 
 # Ensure that all the problem POIs have been taken care of
@@ -544,32 +542,33 @@ print (paste0(nrow(sub[sub$AreaSqKM == 0,]), " POIs on flowlines with no local d
 
 ```{r POI Collapse}
 #  Load data
-if(needs_layer(nav_gpkg, final_poi_layer)) {
+if(needs_layer(ref_gpkg, final_poi_layer)) {
 
   if(!"Type_Elev" %in% names(final_POIs)) final_POIs$Type_Elev <- rep(NA_real_, nrow(final_POIs))
   
   #1 Move POIs downstream by category
-  out_HUC12 <- POI_move_down(nav_gpkg, final_POIs, nsegments, filter(nhd_Final, !is.na(POI_ID)), "Type_HUC12", .10)
-  out_gages <- POI_move_down(nav_gpkg, out_HUC12$allPOIs, out_HUC12$segs, out_HUC12$FL, "Type_Gages", .05)
+  out_HUC12 <- POI_move_down(temp_gpkg, final_POIs, nsegments, filter(nhd_Final, !is.na(POI_ID)), "Type_HUC12", .10)
+  out_gages <- POI_move_down(temp_gpkg, out_HUC12$allPOIs, out_HUC12$segs, out_HUC12$FL, "Type_Gages", .05)
   
   # Convert empty strings to NA for handling within R
   out_gages$allPOIs <- mutate_all(out_gages$allPOIs, list(~na_if(.,"")))
+  
+  out_gages$allPOIs$snapped <- out_gages$allPOIs$COMID %in% new_POIs$COMID
 
   nhd_Final <- select(nhd_Final, -POI_ID) %>%
     left_join(st_drop_geometry(out_HUC12$FL) %>%
                 select(COMID, POI_ID), by = "COMID")
 
   # Write out geopackage layer representing POIs for given theme
-  write_sf(out_gages$allPOIs, nav_gpkg, final_poi_layer)
-  write_sf(nhd_Final, nav_gpkg, nhd_flowline)
-  write_sf(out_gages$segs, nav_gpkg, nsegments_layer)
+  write_sf(out_gages$allPOIs, ref_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(nav_gpkg, final_poi_layer)
-  nhd_Final <- read_sf(nav_gpkg, nhd_flowline)
-  nsegments <- read_sf(nav_gpkg, nsegments_layer)
+  final_POIs <- read_sf(ref_gpkg, final_poi_layer)
+  nhd_Final <- read_sf(ref_gpkg, nhd_flowline)
+  nsegments <- read_sf(temp_gpkg, nsegments_layer)
 }
 
 mapview(final_POIs, layer.name = "All POIs", col.regions = "red") + 

workspace/05_hyRefactor_flines.Rmd

@@ -43,8 +43,7 @@ POIs <- read_sf(nav_gpkg, final_poi_layer) %>%
 
 # Gages that were collapsed during navigate; might not correspond to their Gage_info COMIDs;
 #   We know DAR for those are acceptable, so will keep those out of the events generation.
-POIs_mv <- read_sf(nav_gpkg, poi_moved_layer) %>% 
-  filter(!is.na(Type_Gages))
+POIs_mv <- filter(POIs, snapped)
 
 # Also need to avoid modification to flowlines immediately downstream of POIs
 #      This can cause some hydrologically-incorrect catchment aggregation