Merge branch 'main' into 'main'

Updates to HI to sync up

See merge request wma/nhgf/gfv2.0!130

workspace/02_HI_navigate.Rmd

@@ -23,7 +23,7 @@ library(hyRefactor)
 library(tidyr)
 library(hyfabric)
 
-vpu <- "20"
+vpu <- VPU <- "20"
 
 # Load data paths
 data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
@@ -34,6 +34,8 @@ source("R/NHD_navigate.R")
 source("R/hyRefactor_funs.R")
 source("R/config.R")
 source("R/ExCONUS_config.R")
+
+atts <- readRDS(file.path(data_paths$islands_dir, "nhdplus_flowline_attributes.rds"))
 ```
 
 Define NHD dataset to use for refactoring and discretization.  RPUs 20f and 20g had NHD attribution that behaved differently
@@ -101,7 +103,7 @@ if(needs_layer(nav_gpkg, nhd_flowline)) {
     mutate(dend = ifelse(COMID %in% dend_COM$COMID, 1, 0))
 
   # Bring in catchments
-  nhd_cats <- readRDS(staged_nhdplus$catchment) %>%
+  nhd_cats <- readRDS(staged_nhdplus_islands$catchment) %>%
     st_transform(crs) %>%
     st_make_valid()
 
@@ -423,6 +425,67 @@ if(all(is.na(nav_POIs$Type_NID))) {
 mapview(filter(nav_POIs, !is.na(Type_NID)))
 ```
 
+
+This chunk breaks up potential aggregated segments where the elevation range of the are contributing to the segment exceeds 500-m
+```{r Elevation Break POIs}
+# derive incremental segments from POIs
+inc_segs <- segment_increment(filter(nhd_prep, minNet == 1), filter(st_drop_geometry(nhd),
+                    COMID %in% nav_POIs$COMID, COMID %in% filter(nhd_prep, minNet == 1)$COMID)) %>%
+  # bring over VAA data
+  inner_join(select(atts, COMID, LENGTHKM, MAXELEVSMO, MINELEVSMO, AreaSqKM, TotDASqKM), by = c("COMID" = "ComID")) 
+
+  # Build dataframe for creation of points based on breaks
+elev_pois_split <- inc_segs %>%
+  group_by(POI_ID) %>%
+  # Get elevation info
+  mutate(MAXELEVSMO = na_if(MAXELEVSMO, -9998), MINELEVSMO = na_if(MINELEVSMO, -9998), 
+         elev_diff_seg = max(MAXELEVSMO) - min(MINELEVSMO)) %>%
+  # Filter out to those incremental segments that need splitting above the defined elevation threshold
+  filter((max(MAXELEVSMO) - min(MINELEVSMO)) > (elev_diff * 100)) #%>%
+  # Do not aggregate on fake lfowlines
+  filter(AreaSqKM > 0, TotDASqKM > max_elev_TT_DA) #%>% 
+  # Order from upstream to downstream
+  arrange(desc(Hydroseq)) #%>%
+  # Get attributes used for splitting
+  # csum_length = cumulative length US to DS along segment, cumsum_elev = cumulative US to DS elev diff along segment
+  # inc_elev = elevation range of each segment
+  # iter = estimated number of times we'd need to split existing agg flowlines, or number of rows in set
+  mutate(csum_length = cumsum(LENGTHKM),
+         inc_elev = cumsum(MAXELEVSMO - MINELEVSMO),
+         #nc_elev_diff = c(inc_elev[1], (inc_elev - lag(inc_elev))[-1]),
+         iter = min(round(max(inc_elev) / (elev_diff * 100)), n()),
+         elev_POI_ID = NA) %>%
+  # Remove segments we can't break
+  filter(iter > 0, n() > 1) %>%
+  # Track original iteration number
+  mutate(orig_iter = iter) %>%
+  ungroup() 
+
+if(nrow(elev_pois_split) > 0) {
+  
+  # For reach iteration
+  elev_POIs <- do.call("rbind", lapply(unique(elev_pois_split$POI_ID), split_elev, 
+                                      elev_pois_split, elev_diff*100, max_elev_TT_DA)) %>%
+    filter(!elev_POI_ID %in% tmp_POIs$COMID, COMID == elev_POI_ID) %>%
+    mutate(Type_Elev = 1) %>%
+    select(COMID, Type_Elev) %>%
+    unique()
+  
+  if(nrow(elev_POIs) > 0){
+    tmp_POIs <- POI_creation(select(elev_POIs, COMID, Type_Elev), nhd, "Elev") %>%
+      addType(., tmp_POIs, "Elev")
+  } 
+  
+  write_sf(tmp_POIs, nav_gpkg, nav_poi_layer)
+} else {
+  tmp_POIs <- read_sf(nav_gpkg, nav_poi_layer) 
+  nhd <- read_sf(nav_gpkg, nhd_flowline)
+}
+
+if(!all(is.na(tmp_POIs$Type_Elev)))
+  mapview(filter(tmp_POIs, Type_Elev == 1), layer.name = "Elevation break POIs", col.regions = "blue") 
+```
+
 Clean up and map final POIs
 
 ```{r}

workspace/05_hyRefactor_HI.Rmd

@@ -86,7 +86,7 @@ POIs <- read_sf(nav_gpkg, final_poi_layer) %>%
 if (rpu_code == "20f"){
   new_outlet_POI_atts <- filter(st_drop_geometry(nhd), COMID == 800001119) %>%
     select(COMID, Type_Conf_new = LevelPathI)
-  
+
   POIs <- POIs %>%
     left_join(new_outlet_POI_atts, by = "COMID") %>%
     mutate(Type_Conf = ifelse(!is.na(Type_Conf_new), Type_Conf_new, Type_Conf)) %>%
@@ -126,12 +126,13 @@ avoid <- dplyr::filter(nhd, (sqrt(AreaSqKM) / LENGTHKM) > 3 & AreaSqKM > 1)
 # Verified that it works on all RPUS
 events <- read_sf(nav_gpkg, gage_events_layer) %>%
   #filter(Split == 1) %>%
-  inner_join(select(st_drop_geometry(filter(nhd, refactor == 1)), AreaSqKM, LENGTHKM, COMID, FromMeas, ToMeas), by = "COMID") %>%
+  inner_join(select(st_drop_geometry(nhd), AreaSqKM, LENGTHKM, COMID, FromMeas, ToMeas), by = "COMID") %>%
   #filter(AreaSqKM > 1 & LENGTHKM > (combine_meters / 250)) %>%
-  filter(site_no == "16265600" | REACH_meas - FromMeas > 5 & AreaSqKM > 0.5 & ToMeas - REACH_meas > 5) %>% 
+  filter(site_no %in% c("16265600", "16759600", "16700600", "16682000", "16681000") | 
+           REACH_meas - FromMeas > 5 & AreaSqKM > 0.5 & ToMeas - REACH_meas > 5) %>% 
   select(site_no, COMID, REACHCODE, REACH_meas) %>%
   # Events cannot be in terminal POIs, code seems to ignore the command not to split/combine those flowlines
-  filter(!COMID %in% filter(outlets, type == "terminal")$COMID) %>%
+  #filter(!COMID %in% filter(outlets, type == "terminal")$COMID) %>%
   filter(!COMID %in% avoid$COMID)
 
 # write out events
@@ -221,7 +222,7 @@ if(needs_layer(out_refac_gpkg, divide_layer)) {
                                          fline_rec = reconciled,
                                          fdr = fdr,
                                          fac = fac,
-                                         para = para_reconcile, 
+                                         para = 6, 
                                          cache = cache_split,
                                          keep = NULL,
                                          vector_crs = crs) 

workspace/R/ExCONUS_config.R

@@ -25,6 +25,8 @@ if(vpu == "20"){
   min_da_km <- 5
   split_meters <- 2500
   combine_meters <- 500
+  elev_diff <- 300
+  max_elev_DA <- 1.5
 } else {
   # Cusomizations for AK
   crs <- 3338