diff --git a/workspace/02_HI_navigate.Rmd b/workspace/02_HI_navigate.Rmd
index 1af53d2a52d523a5a8b1b3fcf9de0b988c0d0fac..8e17e7066e36e45e87db049cab71ce2127b6b730 100644
--- a/workspace/02_HI_navigate.Rmd
+++ b/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}
diff --git a/workspace/05_hyRefactor_HI.Rmd b/workspace/05_hyRefactor_HI.Rmd
index fd6c7ac2e39985488f190b1388ddb53afa0294a5..87d5a95bdd93c4279831c58068584093a0a0e84d 100644
--- a/workspace/05_hyRefactor_HI.Rmd
+++ b/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)
diff --git a/workspace/R/ExCONUS_config.R b/workspace/R/ExCONUS_config.R
index 37ea5fc72f7334406584634890ab8b4a9a39ba21..9888dbf9c8383530ad4782b9e48b3bf2d36c54c7 100644
--- a/workspace/R/ExCONUS_config.R
+++ b/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