From 5be9b6d00dbf1949a8f15f842aacfb56ae7241b9 Mon Sep 17 00:00:00 2001
From: David Blodgett <dblodgett@usgs.gov>
Date: Mon, 13 Jun 2022 20:29:03 -0500
Subject: [PATCH] add aggregate_da_thresh_sqkm config parameter

---
 workspace/06-2_aggregate_cats.Rmd | 16 ++++++++--------
 workspace/R/config.R              |  8 ++++++--
 2 files changed, 14 insertions(+), 10 deletions(-)

diff --git a/workspace/06-2_aggregate_cats.Rmd b/workspace/06-2_aggregate_cats.Rmd
index a138ee6..5856968 100644
--- a/workspace/06-2_aggregate_cats.Rmd
+++ b/workspace/06-2_aggregate_cats.Rmd
@@ -64,11 +64,11 @@ reconciled_DA <- select(reconciled, ID, toID) %>%
 if(needs_layer(out_agg_gpkg, mapped_outlets_layer)) {
   
   # subset reconciled flowlines by minimum drainage area for aggcats
-  reconciled_sub <- filter(reconciled, TotDASqKM > 1) 
+  reconciled_sub <- filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm) 
   
   # Get the end of the reconciled flowlines and bind with data frame
-  reconciled_end <- get_node(filter(reconciled, TotDASqKM > 1) , "end") %>%
-    cbind(st_drop_geometry(filter(reconciled, TotDASqKM > 1) ))
+  reconciled_end <- get_node(filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm) , "end") %>%
+    cbind(st_drop_geometry(filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm) ))
   
   # If there are events on split flowlines
   if(nrow(read_sf(out_refac_gpkg, split_layer)) > 0){
@@ -82,11 +82,11 @@ if(needs_layer(out_agg_gpkg, mapped_outlets_layer)) {
     poi_outlets <- filter(outlets_all, !Type_Gages %in% events$site_no) %>%
       # Reset Type_Gages attrbute to NA if event exists
       distinct(.keep_all = TRUE) %>%
-      filter(TotDASqKM > 1)
+      filter(TotDASqKM > aggregate_da_thresh_sqkm)
     
     # Map non-event outlets ID to reconciled ID 
     poi_outlets_mapped <- map_outlet_ids(select(st_drop_geometry(poi_outlets), COMID, type), 
-                                     filter(reconciled, TotDASqKM > 1)) %>%
+                                     filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm)) %>%
       mutate(COMID = as.integer(COMID)) %>%
       inner_join(select(outlets_all, -type), by = "COMID") %>%
       select(-TotDASqKM) %>%
@@ -141,11 +141,11 @@ if(needs_layer(out_agg_gpkg, mapped_outlets_layer)) {
   } else {
     # Map non-event outlets ID into one dataset here with reconciled ID 
     poi_outlets <- distinct(outlets_all, .keep_all = TRUE) %>%
-      filter(TotDASqKM > 1)
+      filter(TotDASqKM > aggregate_da_thresh_sqkm)
     
     # Map non-event outlets ID to reconciled ID 
     outlets_POI<- map_outlet_ids(select(st_drop_geometry(poi_outlets), COMID, type), 
-                                     filter(reconciled, TotDASqKM > 1)) %>%
+                                     filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm)) %>%
       inner_join(select(outlets_all, -type), by = "COMID") %>%
       select(-TotDASqKM) %>%
       st_as_sf()
@@ -185,7 +185,7 @@ if(needs_layer(out_agg_gpkg, agg_cats_layer)){
   agg_cats <- aggregate_catchments(flowpath = reconciled,
                                    divide = divides,
                                    outlets = outlets,
-                                   da_thresh = 1,
+                                   da_thresh = aggregate_da_thresh_sqkm,
                                    only_larger = TRUE,
                                    post_mortem_file = cache_split)
     
diff --git a/workspace/R/config.R b/workspace/R/config.R
index 36a3dc2..1bc57df 100644
--- a/workspace/R/config.R
+++ b/workspace/R/config.R
@@ -25,6 +25,7 @@ if(!exists("rpu_code")) {
     vpu <- get_rpu_dependent_vars()
     rpu_code <- vpu$r
     vpu <- vpu$v
+    rpu_codes <- rpu_vpu[rpu_vpu$vpuid %in% vpu, ]
   }
 }
 
@@ -69,7 +70,7 @@ travt_diff <- 24 # Max number of hours between adjacent POIS
 
 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
+WBs_layer <-  paste0("WB") # 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
@@ -81,6 +82,8 @@ split_meters <- 10000
 combine_meters <- 1000
 min_da_km <- 20
 
+aggregate_da_thresh_sqkm <- 1
+
 # parallel factor for catchment reconciliation.
 para_reconcile <- 2
 para_split_flines <- 2
@@ -101,6 +104,7 @@ lookup_table_refactor <- "lookup_table"
 
 # output geopackage file names
 nav_gpkg <- file.path("cache", paste0("reference_", vpu,".gpkg"))
+temp_gpkg <- file.path("cache", paste0("nav_", vpu,".gpkg"))
 
 rfc_gpkg <- file.path("cache", paste0("refactor_", vpu, ".gpkg"))
 
@@ -108,7 +112,7 @@ gf_gpkg <- file.path("cache", paste0("GF_", vpu, ".gpkg"))
 gf_gpkg_conus <- "cache/reference_CONUS.gpkg"
 
 # Defined during NonDend.Rmd
-ND_gpkg <- file.path("cache", paste0("ND_", vpu,".gpkg"))
+ND_gpkg <- file.path("temp", paste0("ND_", vpu,".gpkg"))
 
 divides_xwalk <- paste0("divides_nhd")
 HRU_layer <- paste0("all_cats")
-- 
GitLab