diff --git a/workspace/06-2_aggregate_cats.Rmd b/workspace/06-2_aggregate_cats.Rmd
index d3cbaf52ef61bd6e25c3eb29daa7f85cd8ed2183..7180fe2d3e3d3cc5695b5754921457a2e0ea2063 100644
--- a/workspace/06-2_aggregate_cats.Rmd
+++ b/workspace/06-2_aggregate_cats.Rmd
@@ -24,6 +24,9 @@ source("R/utils.R")
source("R/config.R")
source("R/NHD_navigate.R")
source("R/hyRefactor_funs.R")
+
+old_path <- Sys.getenv("PATH")
+Sys.setenv(PATH = paste(old_path, "C:\\Users\\abock\\AppData\\Roaming\\npm", sep = ";"))
```
Derive set of outlets to aggregate catchments for. There are three sets of
@@ -37,39 +40,49 @@ nhd <- read_sf(out_refac_gpkg, nhd_flowline)
# Load outlets
outlets_POI <- read_sf(out_refac_gpkg, outlets_layer) %>%
- distinct()
+ distinct() %>%
+ filter(hy_id %in% nhd$COMID)
+
+POIs_data <- read_sf(nav_gpkg, poi_data_table)
# reconciled have no ID applied in refactoring flowlines
reconciled <- st_transform(read_sf(out_refac_gpkg, reconciled_layer), crs)
divides <- st_transform(read_sf(out_refac_gpkg, divide_layer), crs)
+
+lookup <- read_sf(out_refac_gpkg, lookup_table_refactor) %>%
+ inner_join(select(st_drop_geometry(nhd), COMID, TotDASqKM), by = c("NHDPlusV2_COMID" = "COMID"))
+
+highest_DA <- lookup %>%
+ group_by(reconciled_ID) %>%
+ filter(TotDASqKM == max(TotDASqKM))
```
```{r agg_catchments I}
if(needs_layer(out_agg_gpkg, agg_cats_layer)){
# https://github.com/r-spatial/sf/issues/1094#issuecomment-988933580 for why this is here.
- sf::st_geometry(divides[sf::st_is_empty(divides),]) <-
+ sf::st_geometry(divides[sf::st_is_empty(divides),]) <-
sf::st_cast(sf::st_geometry(divides[sf::st_is_empty(divides),]), "MULTIPOLYGON")
-
- #divides <- hyRefactor::clean_geometry(divides)
+
+ #divides <- hyRefactor::clean_geometry(divides)
# POI outlets
outlets <- select(st_drop_geometry(outlets_POI), ID = reconciled_ID, type)
# Identify reconciled flowpaths that are terminal and above DA threshold
- reconciled_DA <- filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm, is.na(toID))
-
+ reconciled_DA <- filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm, is.na(toID))
+
# Ones not in existing pois
missing_outlets <- filter(reconciled_DA, !ID %in% outlets$ID) %>%
pull(ID)
-
+
# bind missing outlets, re-assign mis-attributed terminals
if(length(missing_outlets) > 0){
outlets <- bind_rows(outlets,
data.frame(ID = missing_outlets,
type = rep("terminal", length(missing_outlets))))
write_sf(filter(reconciled, ID %in% missing_outlets), out_agg_gpkg, "missing_outlets")
-
+
outlets <- outlets %>%
left_join(select(st_drop_geometry(reconciled), ID, toID), by = "ID") %>%
mutate(type = ifelse(type == "terminal" & !is.na(toID), "outlet", type)) %>%
@@ -83,10 +96,10 @@ if(needs_layer(out_agg_gpkg, agg_cats_layer)){
da_thresh = aggregate_da_thresh_sqkm,
only_larger = TRUE,
post_mortem_file = cache_split)
-
+
agg_cats$cat_sets$set <- sapply(agg_cats$cat_sets$set, paste, collapse = ",")
agg_cats$fline_sets$set <- sapply(agg_cats$fline_sets$set, paste, collapse = ",")
-
+
write_sf(agg_cats$cat_sets, out_agg_gpkg, agg_cats_layer)
write_sf(agg_cats$fline_sets, out_agg_gpkg, agg_fline_layer)
@@ -96,18 +109,34 @@ if(needs_layer(out_agg_gpkg, agg_cats_layer)){
st_as_sf()
# Build final POI set
- final_POIs <- st_drop_geometry(agg_cats$fline_set) %>%
+ POIs_att <- select(rec_outlets, ID, toID) %>%
# Bring over POI information
- left_join(st_drop_geometry(outlets_POI), by = c("ID" = "reconciled_ID")) %>%
- # Bring over physical end node XY
- left_join(select(rec_outlets, ID), by = "ID") %>%
- # Mark additional POI/outlets created during aggregate cats
- mutate(Type_Con = ifelse(!ID %in% outlets$ID, 1, 0)) %>%
+ inner_join(st_drop_geometry(outlets_POI), by = c("ID" = "reconciled_ID")) %>%
+ select(-c(type, member_COMID))
+
+ POIs_missing <- filter(select(rec_outlets, -toID), !ID %in% POIs_att$ID) %>%
+ inner_join(st_drop_geometry(agg_cats$fline_sets), by = "ID") %>%
+ arrange(ID) %>%
+ mutate(poi_id = max(POIs_att$poi_id) + row_number()) %>%
+ select(ID, toID, poi_id, event_id = event_identifier, Type_Term = orig_levelpathID) %>%
+ inner_join(select(highest_DA, NHDPlusV2_COMID, reconciled_ID), by = c("ID" = "reconciled_ID")) %>%
+ select(ID, toID, hy_id = NHDPlusV2_COMID, poi_id, event_id, Type_Term) %>%
+ mutate(Type_Term = ifelse(is.na(toID), Type_Term, NA),
+ Type_Con = ifelse(!is.na(toID), 1, NA)) %>%
+ st_as_sf()
+
+ final_POIs <- data.table::rbindlist(list(POIs_att, select(POIs_missing,
+ -c(Type_Term, Type_Con))), fill = TRUE) %>%
st_as_sf()
+
+ unaggregated_outlets <- filter(outlets_POI, !reconciled_ID %in% final_POIs$ID)
+ double_outlets <- final_POIs %>% group_by(ID) %>% filter(n() > 1)
+
# Write out
write_sf(final_POIs, out_agg_gpkg, mapped_outlets_layer)
-
+ write_sf(unaggregated_outlets, out_agg_gpkg, "unmapped_outlets")
+ write_sf(double_outlets, out_agg_gpkg, "double_outlets")
} else {
agg_cats <- list(cat_sets = read_sf(out_agg_gpkg, agg_cats_layer),
fline_sets = read_sf(out_agg_gpkg, agg_fline_layer))
@@ -117,6 +146,25 @@ if(needs_layer(out_agg_gpkg, agg_cats_layer)){
```
+```{r Long form POIs}
+POIs <- final_POIs %>%
+ arrange(ID)
+
+final_POI_geom <- POIs %>%
+ select(ID) %>%
+ cbind(st_coordinates(.)) %>%
+ group_by(ID, X, Y) %>%
+ mutate(geom_id = cur_group_id()) %>%
+ ungroup()
+
+final_POIs_table <- POIs %>%
+ inner_join(select(st_drop_geometry(final_POI_geom), -X, -Y), by = "ID") %>%
+ distinct()
+
+write_sf(final_POIs_table, out_agg_gpkg, poi_data_table)
+```
+
+
```{r lookup table}
refactor_lookup <- dplyr::select(st_drop_geometry(reconciled), ID, member_COMID) %>%
dplyr::mutate(member_COMID = strsplit(member_COMID, ",")) %>%