Improved event handling and documentation

workspace/02_NHD_navigate.Rmd

@@ -22,6 +22,8 @@ knitr::opts_chunk$set(
 source("R/utils.R")
 source("R/NHD_navigate.R")
 source("R/hyrefactor_funs.R")
+#source("R/wb_poi_collapse.R")
+#source("R/wb_inlet_collapse.R")
 
 # increase timeout for data downloads
 options(timeout=600)
@@ -29,10 +31,12 @@ options(timeout=600)
 # Load Configuration of environment
 source("R/config.R")
 
+# Gages output from Gage_selection
 gages <- read_sf(gage_info_gpkg, "Gages")
 
+# NWM network
 nc <- RNetCDF::open.nc(data_paths$nwm_network)
-
+# NWM gages
 nwm_gages <- data.frame(
   comid = RNetCDF::var.get.nc(nc, "link"), 
   gage_id = RNetCDF::var.get.nc(nc, "gages")) %>%
@@ -41,9 +45,8 @@ nwm_gages <- data.frame(
 
 RNetCDF::close.nc(nc)
 
-# need some extra attributes
+# need some extra attributes for a few POI analyses
 atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.rds"))
-
 ```
 
 ```{r huc12 POIs}
@@ -51,9 +54,9 @@ 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)), silent = TRUE)
+   try(nhd <- select(nhd, -c(minNet, WB, struct_POI, struct_Net, POI_ID, dend, poi)), silent = TRUE)
   
-  # HUC02 includes some 
+  # Some NHDPlus VPUs include HUCs from other VPUs
   if(vpu == "02"){
     grep_exp <-"^02|^04"
   } else if (vpu == "08") {
@@ -124,6 +127,7 @@ if(all(is.na(tmp_POIs$Type_Gages))) {
   write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
 } else {
   tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+  events <- read_sf(temp_gpkg, split_layer)
 }
 
 mapview(filter(tmp_POIs, !is.na(Type_Gages)), layer.name = "Streamgage POIs", col.regions = "blue") 
@@ -199,6 +203,8 @@ mapview(filter(tmp_POIs, !is.na(Type_Term)), layer.name = "Terminal POIs", col.r
 #  Derive or load Waterbody POIs ----------------------
 if(all(is.na(tmp_POIs$Type_WBOut))) {
   
+  events <- read_sf(temp_gpkg, split_layer)
+  
   # Waterbodies sourced from NHD waterbody layer for entire VPU
   WBs_VPU_all <- filter(readRDS("data/NHDPlusNationalData/nhdplus_waterbodies.rds"), 
                         COMID %in% nhd$WBAREACOMI) %>%
@@ -235,7 +241,12 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
   
   wb_layers <- wbout_POI_creaton(nhd, WBs_VPU, data_paths, crs)
   WBs_VPU <- wb_layers$WBs
-  tmp_POIs <- wb_layers$POIs
+  #tmp_POIs <- wb_layers$POIs
+  
+  wb_out_col <- wb_poi_collapse(wb_layers$POIs, nhd, events)
+  ho <- filter(wb_layers$POIs, !COMID %in% wb_out_col$POIs$COMID)
+  write_sf(wb_out_col$events_ret, temp_gpkg, split_layer)
+  tmp_POIs <- wb_out_col$POIs
   
   if(!all(is.na(wb_layers$events))){
     wb_events <- select(wb_layers$events, -c(id, offset)) %>%
@@ -244,7 +255,8 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
     # Write out events and outlets
     if(!needs_layer(temp_gpkg, split_layer)){
       events <- read_sf(temp_gpkg, split_layer) %>%
-        rbind(st_compatibalize(wb_events,.))
+        rbind(st_compatibalize(wb_events,.)) %>%
+        unique()
       write_sf(events, temp_gpkg, split_layer)
     } else {
       write_sf(wb_events, nav_gpkg, split_layer)
@@ -362,7 +374,11 @@ mapview(filter(tmp_POIs, Type_NID != ""), layer.name = "NID POIs", col.regions =
 if(all(is.na(tmp_POIs$Type_WBIn))) {
 
   wb_layers <- wbin_POIcreation(nhd, WBs_VPU, data_paths, crs, split_layer, tmp_POIs)
-  tmp_POIs <- wb_layers$POIs 
+  #tmp_POIs <- wb_layers$POIs 
+  wb_in_col <- wb_inlet_collapse(wb_layers$POIs, nhd, events)
+  #ho <- filter(wb_layers$POIs, !COMID %in% wb_in_col$POIs$COMID)
+  #write_sf(wb_in_col$events_ret, temp_gpkg, split_layer)
+  tmp_POIs <- wb_in_col$POIs
   
   if(!all(is.na(wb_layers$events))) {
     wb_inlet_events <- select(wb_layers$events, -c(id, offset, Hydroseq, ToMeas)) %>%
@@ -397,7 +413,8 @@ inc_segs <- segment_increment(filter(nhd, minNet == 1),
                               filter(st_drop_geometry(nhd),
                                      COMID %in% tmp_POIs$COMID, COMID %in% filter(nhd, minNet == 1)$COMID)) %>%
   # bring over VAA data
-  inner_join(select(atts, COMID, DnHydroseq, VA_MA, TOTMA, LENGTHKM, MAXELEVSMO, MINELEVSMO, WBAREACOMI, WBAreaType, FTYPE,
+  inner_join(select(atts, COMID, DnHydroseq, VA_MA, TOTMA, LENGTHKM, MAXELEVSMO, 
+                    MINELEVSMO, WBAREACOMI, WBAreaType, FTYPE,
                     AreaSqKM, TotDASqKM), by = "COMID") 
 
 # Build dataframe for creation of points based on breaks
@@ -609,7 +626,7 @@ if(needs_layer(nav_gpkg, final_poi_layer)) {
                              poi_dar_move *2)
   out_HUC12$allPOIs$nexus <- as.numeric(out_HUC12$allPOIs$nexus)
   out_gages <- POI_move_down(temp_gpkg, out_HUC12$allPOIs, out_HUC12$segs, out_HUC12$FL, 
-                             "Type_Gages", poi_dar_move)
+                             "Type_Gages", poi_dar_move) 
   
   # Convert empty strings to NA for handling within R
   out_gages$allPOIs <- mutate_all(out_gages$allPOIs, list(~na_if(.,"")))