From 4c68c7c40c7f124fa34ba8b050d72818552f6d18 Mon Sep 17 00:00:00 2001
From: Bock <abock@usgs.gov>
Date: Thu, 18 May 2023 08:04:04 -0600
Subject: [PATCH] moved some code around
---
workspace/02_NHD_navigate_Sp.Rmd | 823 ++++++++++++++++++++++++++++++
workspace/06-2_aggregate_cats.Rmd | 19 +-
2 files changed, 836 insertions(+), 6 deletions(-)
create mode 100644 workspace/02_NHD_navigate_Sp.Rmd
diff --git a/workspace/02_NHD_navigate_Sp.Rmd b/workspace/02_NHD_navigate_Sp.Rmd
new file mode 100644
index 0000000..daebbea
--- /dev/null
+++ b/workspace/02_NHD_navigate_Sp.Rmd
@@ -0,0 +1,823 @@
+---
+title: "NHD Navigate"
+output: html_document
+editor_options:
+ chunk_output_type: console
+---
+This notebook Generate Segments and POIS from POI data sources and builds
+a minimally-sufficient stream network connecting all POis, as well as generating first-cut of
+national segment network.
+
+```{r setup_rmd, echo=FALSE, cache=FALSE}
+knitr::opts_chunk$set(
+ collapse = TRUE,
+ comment = "#>",
+ fig.width=6,
+ fig.height=4,
+ cache=FALSE)
+```
+
+```{r setup}
+# Load custom functions
+source("R/utils.R")
+source("R/NHD_navigate.R")
+source("R/hyrefactor_funs.R")
+
+# increase timeout for data downloads
+options(timeout=600)
+
+# Load Configuration of environment
+source("R/config.R")
+
+# Gages output from Gage_selection
+gages <- read_sf(gage_info_gpkg, "Gages")
+# need some extra attributes for a few POI analyses
+atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.rds"))
+
+sparrow_dir <- "data/sparrow"
+
+collapse <- TRUE
+```
+
+```{r huc12 POIs}
+# Derive or load HUC12 POIs
+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, dend, poi)), silent = TRUE)
+
+ # Some NHDPlus VPUs include HUCs from other VPUs
+ if(vpu == "02"){
+ grep_exp <-"^02|^04"
+ } else if (vpu == "08") {
+ grep_exp <- "^03|^08"
+ } else {
+ grep_exp <- paste0("^", substr(vpu, start = 1, stop = 2))
+ }
+
+ # Join HUC12 outlets with NHD
+ HUC12_COMIDs <- read_sf(data_paths$hu12_points_path, "hu_points") %>%
+ filter(grepl(grep_exp, .data$HUC12)) %>%
+ select(COMID, HUC12) %>%
+ # Remove this when HUC12 outlets finished
+ group_by(COMID) %>%
+ filter(row_number() == 1) %>%
+ ungroup()
+
+ # Create POIs - some r05 HUC12 POIs not in R05 NHD
+ huc12_POIs <- POI_creation(st_drop_geometry(HUC12_COMIDs), nhd, "HUC12")
+ tmp_POIs <- huc12_POIs
+
+ # Write out geopackage layer representing POIs for given theme
+ write_sf(huc12_POIs, temp_gpkg, nav_poi_layer)
+ tmp_POIs <- huc12_POIs
+} else {
+ # Load HUC12 POIs as the tmpPOIs if they already exist
+ tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+ nhd <- read_sf(nav_gpkg, nhd_flowline)
+}
+
+mapview(filter(tmp_POIs, Type_HUC12 != ""), layer.name = "HUC12 POIs", col.regions = "red")
+```
+
+```{r streamgage POIs}
+if(!"Type_Gages" %in% names(tmp_POIs)) {
+
+ ref_gages <- read_sf("data/ref_gages.gpkg")
+
+ # Previously identified streamgages within Gage_Selection.Rmd
+ streamgages_VPU <- gages %>%
+ rename(COMID = comid) %>%
+ filter(COMID %in% nhd$COMID) %>%
+ #st_drop_geometry() %>%
+ switchDiv(., nhd)
+
+ streamgages <- streamgages_VPU %>%
+ group_by(COMID) %>%
+ # If multiple gages per COMID, pick one with highest rank as rep POI_ID
+ filter(gage_score == max(gage_score), !is.na(drain_area)) %>%
+ ungroup() %>%
+ select(COMID, site_no)
+
+ # Derive TE POIs
+ tmp_POIs <- POI_creation(st_drop_geometry(streamgages), nhd, "Gages") %>%
+ addType(., tmp_POIs, "Gages", nexus = FALSE)
+
+ # Write out geopackage layer representing POIs for given theme
+ write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+} else {
+ tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+}
+
+mapview(filter(tmp_POIs, !is.na(Type_Gages)), layer.name = "Streamgage POIs", col.regions = "blue")
+```
+
+# R01 - missing one
+```{r Sparrow Cal POIs}
+if(!"Type_Sparrow" %in% names(tmp_POIs)) {
+
+ # Load Sparrow calibration POIs
+ sparrow_cal_pois <- read_sf("data/sparrow/sparrow_source.gpkg", "calibration_poi") %>%
+ mutate(COMID = as.integer(COMID)) %>%
+ filter(COMID %in% nhd$COMID) %>%
+ select(COMID, MonitoringLocationIdentifier, site_no, station_nm,
+ usgs_hy_id, distance_m, longitude, latitude, nat_id) %>%
+ st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>%
+ select(COMID, nat_id)
+
+ # Derive TE POIs
+ tmp_POIs <- POI_creation(st_drop_geometry(sparrow_cal_pois), nhd, "Sparrow") %>%
+ addType(., tmp_POIs, "Sparrow", nexus = FALSE)
+
+ missing_sparrow <- filter(sparrow_cal_pois, !nat_id %in% tmp_POIs$Type_Sparrow)
+ if(nrow(missing_sparrow) > 0){write_sf(missing_sparrow, temp_gpkg, "missing_sparrow")}
+ write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+} else {
+ tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+}
+
+mapview(filter(tmp_POIs, !is.na(Type_Gages)), layer.name = "Streamgage POIs", col.regions = "blue")
+```
+
+```{r TE POIs}
+if(!"Type_TE" %in% names(tmp_POIs)) {
+
+ if(vpu == "08"){
+ nhd$VPUID <- "08"
+ } else {
+ nhd$VPUID <- substr(nhd$RPUID, 1, 2)
+ }
+
+ # Read in Thermoelectric shapefile
+ TE_COMIDs <- read_sf(data_paths$TE_points_path, "2015_TE_Model_Estimates_lat.long_COMIDs") %>%
+ mutate(COMID = as.integer(COMID)) %>%
+ inner_join(., select(st_drop_geometry(nhd), COMID, VPUID), by = "COMID") %>%
+ filter(grepl(paste0("^", substr(vpu, 1, 2), ".*"), .data$VPUID), COMID > 0) %>%
+ switchDiv(., nhd) %>%
+ group_by(COMID) %>%
+ summarize(EIA_PLANT = paste0(unique(EIA_PLANT_), collapse = " "), count = n()) %>%
+ ungroup() %>%
+ select(COMID, EIA_PLANT)
+
+ # Derive TE POIs
+ tmp_POIs <- POI_creation(st_drop_geometry(TE_COMIDs), nhd, "TE") %>%
+ addType(., tmp_POIs, "TE", nexus = FALSE, bind = TRUE)
+
+ # As a fail-safe, write out list of TE plants not assigned a POI
+ if(nrow(filter(TE_COMIDs, !COMID %in% tmp_POIs$COMID)) > 0) {
+ write_sf(filter(TE_COMIDs, !COMID %in% tmp_POIs$COMID),
+ temp_gpkg, "unassigned_TE")
+ }
+
+ # Write out geopackage layer representing POIs for given theme
+ write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+} else {
+ # Load TE POIs if they already exist
+ tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+}
+
+mapview(filter(tmp_POIs, Type_TE != ""), layer.name = "TE Plant POIs", col.regions = "blue")
+```
+
+```{r waterbody outlet POIs}
+# Derive or load Waterbody POIs ----------------------
+if(!"Type_WBOut" %in% names(tmp_POIs)) {
+
+ sparrow_wb_outlet <- read_sf("data/sparrow/sparrow_source.gpkg", "wb_outlets") %>%
+ filter(COMID %in% nhd$COMID)
+
+ sparrow_wb <- read_sf("data/sparrow/sparrow_source.gpkg", "waterbodies") %>%
+ filter(COMID %in% nhd$WBAREACOMI)
+
+ # Waterbodies sourced from NHD waterbody layer for entire VPU
+ WBs_VPU_all <- filter(readRDS("data/NHDPlusNationalData/nhdplus_waterbodies.rds"),
+ COMID %in% nhd$WBAREACOMI) %>%
+ mutate(FTYPE = as.character(FTYPE)) %>%
+ mutate(source = "NHDv2WB") %>%
+ st_transform(crs)
+
+ # Waterbody list that strictly meet the size criteria
+ wb_lst <- st_drop_geometry(WBs_VPU_all) %>%
+ dplyr::filter(FTYPE %in% c("LakePond", "Reservoir") & AREASQKM >= (min_da_km/2)) %>%
+ dplyr::select(COMID) %>%
+ filter(!COMID == 166410600) %>%
+ rbind(dplyr::select(st_drop_geometry(sparrow_wb), COMID))
+
+ if (file.exists(data_paths$resops_NID_CW)){
+ # ResOpsUS locations in the VPU waterbody set
+ resops_wb_df <- read.csv(data_paths$resops_NID_CW) %>%
+ dplyr::filter(WBAREACOMI %in% WBs_VPU_all$COMID) %>%
+ dplyr::select(DAM_ID, DAM_NAME, COMID = WBAREACOMI, FlowLcomid) %>%
+ distinct()
+
+ # Add ResOPsUS locations to waterbody list
+ if(nrow(resops_wb_df) > 0){
+ wb_lst <- rbind(wb_lst, select(resops_wb_df, COMID)) %>%
+ distinct()
+ }
+ }
+
+ # Waterbodies that meet the size criteria and/or are in ResOpsUS datasets
+ WBs_VPU <- WBs_VPU_all %>%
+ dplyr::filter(COMID %in% wb_lst$COMID)
+
+ WBs_VPU <- filter(WBs_VPU, COMID != 666170) #1R16
+
+ # Attribute flowlines that are in waterbodies
+ nhd <- mutate(nhd, WB = ifelse(WBAREACOMI > 0 & WBAREACOMI %in% WBs_VPU$COMID, 1, 0))
+
+ wbout_COMIDs <- nhd %>% #filter(nhd, dend == 1 & WB == 1)
+ filter(WBAREACOMI %in% WBs_VPU$COMID) %>%
+ filter(!WBAREACOMI %in% sparrow_wb_outlet$WB_COMID) %>%
+ group_by(WBAREACOMI) %>%
+ slice(which.min(Hydroseq)) %>%
+ switchDiv(., nhd) %>%
+ select(COMID, WBAREACOMI, LevelPathI) %>%
+ st_drop_geometry() %>%
+ rbind(select(st_drop_geometry(sparrow_wb_outlet), COMID, WBAREACOMI = WB_COMID))
+
+ tmp_POIs <- POI_creation(wbout_COMIDs, nhd, "WBOut") %>%
+ addType(., tmp_POIs, "WBOut", nexus = FALSE, bind = TRUE)
+
+ # Write out waterbodies
+ write_sf(WBs_VPU, nav_gpkg, WBs_layer)
+
+ # Write specific ResOpsUS data to a separate sheet
+ if(nrow(wb_lst) > 0){
+ resops_POIs_df <- st_drop_geometry(tmp_POIs) %>%
+ dplyr::select(COMID, Type_WBOut) %>%
+ dplyr::filter(!is.na(Type_WBOut)) %>%
+ dplyr::inner_join(select(resops_wb_df, -FlowLcomid), by = c("Type_WBOut" = "COMID"))
+
+ write.csv(resops_POIs_df, file.path("data/reservoir_Data", paste0("resops_POIs_",vpu,".csv")))
+ }
+
+
+ wb_out <- filter(tmp_POIs, !is.na(Type_WBOut))
+ missing_wb_out <- filter(sparrow_wb_outlet, !WB_COMID %in% wb_out$Type_WBOut)
+ missing_wb_out_com <- filter(sparrow_wb_outlet,
+ !COMID %in% wb_out$COMID)
+ if(nrow(missing_wb_out) > 0){write_sf(missing_wb_out, temp_gpkg, "missing_sparrow_WB")}
+ if(nrow(missing_wb_out_com) > 0){write_sf(missing_wb_out_com, temp_gpkg, "missing_sparrow_WBout")}
+
+ sparrow_wb_ids <- select(sparrow_wb_outlet, COMID, WB_COMID, nat_id)
+
+ tmp_POIs <- tmp_POIs %>%
+ left_join(select(st_drop_geometry(sparrow_wb_ids), -COMID), by = c("Type_WBOut" = "WB_COMID")) %>%
+ mutate(Type_Sparrow = ifelse(is.na(Type_Sparrow), nat_id, Type_Sparrow)) %>%
+ select(-nat_id)
+
+ write_sf(nhd, nav_gpkg, nhd_flowline)
+ write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+} else {
+ tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+ nhd <- read_sf(nav_gpkg, nhd_flowline)
+ WBs_VPU <- read_sf(nav_gpkg, WBs_layer)
+ resops_POIs_df <- read.csv(file.path("data/reservoir_Data", paste0("resops_POIs_",vpu,".csv")))
+}
+
+mapview(filter(tmp_POIs, !is.na(Type_WBOut)), layer.name = "Waterbody outlet POIs", col.regions = "red")
+```
+
+```{r NID POIs}
+# Derive or load NID POIs ----------------------
+if(all(is.na(tmp_POIs$Type_NID))) {
+
+ sparrow_nid <- read_sf("data/sparrow/sparrow_source.gpkg", "nid_pois")
+
+ sparrow_wb_nid <- sparrow_nid %>%
+ #rename(COMID = ComID) %>%
+ filter(COMID %in% nhd$COMID) %>%
+ select(COMID, NIDID, nat_id)
+
+ if(nrow(sparrow_wb_nid) > 0){
+ # Determine which NID facilitites have been co-located with ResOps
+ res_ops_table <- read.csv(data_paths$resops_NID_CW) %>%
+ mutate(new_WBARECOMI = ifelse(is.na(WBAREACOMI) & is.na(NHDAREA), HR_NHDPLUSID, WBAREACOMI)) %>%
+ mutate(new_WBARECOMI = ifelse(is.na(new_WBARECOMI) & is.na(HR_NHDPLUSID), NHDAREA, new_WBARECOMI)) %>%
+ filter(new_WBARECOMI %in% tmp_POIs$Type_WBOut, !is.na(new_WBARECOMI)) %>%
+ select(NID_ID = NID_Source_FeatureID, Type_WBOut = new_WBARECOMI) %>%
+ filter()
+
+ # Attribute those NID facilities precisely at the outlet
+ if(nrow(res_ops_table) > 0){
+ tmp_POIs <- tmp_POIs %>%
+ left_join(res_ops_table, by = "Type_WBOut") %>%
+ mutate(Type_NID = ifelse(!is.na(NID_ID), NID_ID, NA)) %>%
+ select(-NID_ID)
+
+ sparrow_wb_nid <- filter(sparrow_wb_nid, !NIDID %in% tmp_POIs$Type_NID)
+ }
+
+ # Derive other NID POIs
+ tmp_POIs <- POI_creation(st_drop_geometry(sparrow_wb_nid), nhd, "NID") %>%
+ addType(., tmp_POIs, "NID", nexus = FALSE, bind = TRUE) %>%
+ left_join(select(st_drop_geometry(sparrow_wb_nid), NIDID, nat_id), by = c("Type_NID" = "NIDID")) %>%
+ mutate(Type_Sparrow = ifelse(is.na(Type_Sparrow), nat_id, Type_Sparrow)) %>%
+ select(-nat_id)
+
+ sparrow_wb_nid_missing <- filter(sparrow_wb_nid, !nat_id %in% tmp_POIs$Type_Sparrow) %>%
+ mutate(WB_Comid = paste("Sparrow_nid_", row_number()))
+
+ if(nrow(sparrow_wb_nid_missing) > 0){write_sf(sparrow_wb_nid_missing, temp_gpkg, "missing_wb_nid")}
+ # Write out geopackage layer representing POIs for given theme
+ write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+ mapview(filter(tmp_POIs, Type_NID != ""), layer.name = "NID POIs", col.regions = "blue")
+ }
+
+} else {
+ # Load NID POIs if they already exist
+ tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+}
+```
+
+
+```{r Diversion POIs}
+# # Derive POIs at confluences where they are absent ----------------------
+# if(needs_layer(temp_gpkg, "diversions")) {
+#
+# div_pois_source <- read_sf("data/sparrow/sparrow_source.gpkg", "div_poi") %>%
+# filter(COMID %in% nhd$COMID) %>%
+# select(COMID, nat_id)
+#
+# div_POIs <- POI_creation(st_drop_geometry(div_pois_source), nhd, "Div") %>%
+# mutate(Type_Sparrow = NA, nexus = NA) %>%
+# st_compatibalize(., tmp_POIs)
+#
+# tmp_POIs <- data.table::rbindlist(list(tmp_POIs,
+# select(div_POIs, -c(Type_Div))), fill = TRUE) %>%
+# st_as_sf() %>%
+# left_join(select(st_drop_geometry(div_pois_source), COMID, nat_id), by = "COMID") %>%
+# mutate(Type_Sparrow = ifelse(is.na(Type_Sparrow), nat_id, Type_Sparrow)) %>%
+# select(-nat_id)
+#
+# missing_div <- filter(div_pois_source, !nat_id %in% tmp_POIs$Type_Sparrow)
+# if(nrow(missing_div) > 0){write_sf(missing_div, temp_gpkg, "missing_div")}
+#
+# write_sf(div_pois_source, temp_gpkg, "diversions")
+# # Write out geopackage layer representing POIs for given theme
+# write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+# } else {
+# # Load NID POIs if they already exist
+# tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+# }
+```
+
+```{r Confluence POIs}
+# # Derive POIs at confluences where they are absent ----------------------
+# if(all(is.na(tmp_POIs$Type_Conf)) | !"Type_Conf" %in% colnames(tmp_POIs)) {
+#
+# # Navigate upstream from each POI and determine minimally-sufficient network between current POI sets
+# up_net <- unique(unlist(lapply(unique(tmp_POIs$COMID), NetworkNav, nhd)))
+# finalNet <- unique(NetworkConnection(up_net, st_drop_geometry(nhd)))
+#
+# # Subset NHDPlusV2 flowlines to navigation results and write to shapefile
+# nhd <- mutate(nhd, minNet = ifelse(COMID %in% finalNet, 1, 0))
+#
+# # Create new confluence POIs
+# conf_COMIDs <- st_drop_geometry(filter(nhd, minNet == 1)) %>%
+# # Downstream hydrosequence of 0 indicates
+# # the flowline is terminating or
+# # leaving the domain, so they
+# # are excluded from this process
+# filter(DnHydroseq > 0) %>%
+# group_by(DnHydroseq) %>%
+# filter(n()> 1) %>%
+# mutate(Type_Conf = LevelPathI) %>%
+# ungroup() %>%
+# select(COMID, Type_Conf)
+#
+# tmp_POIs <- POI_creation2(conf_COMIDs, filter(nhd, minNet == 1), "Conf") %>%
+# addType2(., tmp_POIs, "Conf", nexus = FALSE, bind = TRUE)
+#
+# write_sf(nhd, nav_gpkg, nhd_flowline)
+# write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+# } else {
+# nhd <- read_sf(nav_gpkg, nhd_flowline)
+# tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+# }
+#
+# mapview(filter(tmp_POIs, !is.na(Type_Conf)), layer.nameble = "Confluence POIs", col.regions = "blue")
+```
+
+```{r Terminal POIs}
+# # Derive or load Terminal POIs ----------------------
+# if(all(is.na(tmp_POIs$Type_Term))) {
+#
+# # Terminal POIs on levelpaths with upstream POIs
+# term_paths <- filter(st_drop_geometry(nhd), Hydroseq %in% filter(nhd, COMID %in% tmp_POIs$COMID)$TerminalPa)
+#
+# # Non-POI levelpath terminal pois, but meet size criteria
+# terminal_POIs <- st_drop_geometry(nhd) %>%
+# filter(Hydroseq == TerminalPa | (toCOMID == 0 | is.na(toCOMID) | !toCOMID %in% COMID)) %>%
+# filter(!COMID %in% term_paths$COMID, TotDASqKM >= min_da_km) %>%
+# bind_rows(term_paths) %>%
+# # Use level path identifier as Type ID
+# select(COMID, LevelPathI)
+#
+# tmp_POIs <- POI_creation2(terminal_POIs, nhd, "Term") %>%
+# addType2(., tmp_POIs, "Term", nexus = FALSE, bind = TRUE)
+#
+# write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+# } else {
+# tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+# nhd <- read_sf(nav_gpkg, nhd_flowline)
+# }
+#
+# mapview(filter(tmp_POIs, !is.na(Type_Term)), layer.name = "Terminal POIs", col.regions = "blue")
+```
+
+```{r waterbody inlet POIs}
+# # Derive or load Waterbody POIs ----------------------
+# if(all(is.na(tmp_POIs$Type_WBIn))) {
+#
+# sparrow_inlets <- read_sf("data/sparrow/sparrow_source.gpkg", "wb_inlets") %>%
+# filter(COMID %in% nhd$COMID)
+#
+# WBs_VPU <- filter(WBs_VPU, COMID != 666170) #1R16
+#
+# wbout_COMIDs <- filter(tmp_POIs, !is.na(Type_WBOut))
+# nhd_out <- filter(nhd, COMID %in% wbout_COMIDs$COMID)
+#
+# # Create waterbody inlet POIs
+# wbin_COMIDs <- filter(nhd, WB == 0,
+# DnHydroseq %in% filter(nhd, WB == 1)$Hydroseq) %>%
+# select(-WBAREACOMI) %>%
+# switchDiv(., nhd) %>%
+# inner_join(st_drop_geometry(filter(nhd, minNet == 1)) %>%
+# select(Hydroseq, WBAREACOMI), by = c("DnHydroseq" = "Hydroseq")) %>%
+# select(COMID, WBAREACOMI, minNet) %>%
+# group_by(COMID) %>%
+# # Ensure the inlets are on the network defined by confluence POIs
+# filter(minNet == 1) %>%
+# select(COMID, WBAREACOMI)
+#
+# # Headwater Waterbodies that may need the network extended through the inlet
+# need_wbin <- st_drop_geometry(filter(WBs_VPU, source == "NHDv2WB")) %>%
+# dplyr::select(COMID)%>%
+# dplyr::filter(!COMID %in% wbin_COMIDs$WBAREACOMI)
+#
+# up_lp <- filter(nhd, LevelPathI %in% nhd_out$LevelPathI,
+# WBAREACOMI %in% need_wbin$COMID) %>%
+# group_by(WBAREACOMI) %>%
+# filter(n() > 1, StartFlag != 1) %>%
+# slice(which.min(Hydroseq))
+#
+# wbin_pois <- select(nhd, COMID, DnHydroseq, LevelPathI) %>%
+# inner_join(select(st_drop_geometry(up_lp), Hydroseq, WBAREACOMI, LevelPathI),
+# by = c("DnHydroseq" = "Hydroseq", "LevelPathI" = "LevelPathI")) %>%
+# select(COMID, WBAREACOMI) %>%
+# rbind(wbin_COMIDs)
+#
+# tmp_POIs <- POI_creation2(st_drop_geometry(wbin_pois), nhd, "WBIn") %>%
+# addType2(., tmp_POIs, "WBIn", nexus = FALSE, bind = TRUE)
+#
+# tmp_POIs <- tmp_POIs %>%
+# left_join(select(st_drop_geometry(sparrow_inlets), WB_COMID, nat_id),
+# by = c("Type_WBIn" = "WB_COMID"), multiple = "all") %>%
+# mutate(Type_Sparrow = ifelse(is.na(Type_Sparrow), nat_id, Type_Sparrow)) %>%
+# select(-nat_id)
+#
+# missing_wb_in <- filter(sparrow_inlets, !nat_id %in% tmp_POIs$Type_Sparrow |
+# !WB_COMID %in% tmp_POIs$Type_WBIn)
+#
+# if(nrow(missing_wb_in) > 0){write_sf(missing_wb_in, temp_gpkg, "missing_wbin")}
+# #tmp_POIs$nexus <- ifelse(is.na(tmp_POIs$nexus), FALSE, tmp_POIs$nexus)
+# write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+# } else {
+# tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+# nhd <- read_sf(nav_gpkg, nhd_flowline)
+# }
+#
+# mapview(filter(tmp_POIs, !is.na(Type_WBIn)), layer.name = "Waterbody inlet POIs", col.regions = "blue") +
+# mapview(filter(tmp_POIs, !is.na(Type_WBOut)), layer.name = "Waterbody outlet POIs", col.regions = "red")
+```
+
+
+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, 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,
+# AreaSqKM, TotDASqKM), by = "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(MINELEVSMO) - min(MINELEVSMO)) %>%
+# # Filter out to those incremental segments that need splitting above the defined elevation threshold
+# filter((max(MINELEVSMO) - 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(MINELEVSMO - 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")
+# }
+#
+#
+# } else {
+#
+# tmp_POIs$Type_Elev <- rep(NA, nrow(tmp_POIs))
+#
+# write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+#
+# tmp_POIs <- read_sf(temp_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")
+```
+
+```{r Time of travel POIs}
+# if(all(is.na(tmp_POIs$Type_Travel))) {
+#
+# # derive incremental segments from POIs
+# 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,
+# AreaSqKM, TotDASqKM), by = "COMID")
+#
+# # TT POIs
+# tt_pois_split <- inc_segs %>%
+# # Should we substitute with a very small value
+# mutate(VA_MA = ifelse(VA_MA < 0, NA, VA_MA)) %>%
+# mutate(FL_tt_hrs = (LENGTHKM * ft_per_km)/ VA_MA / s_per_hr ) %>%
+# group_by(POI_ID) %>%
+# filter(sum(FL_tt_hrs) > travt_diff, TotDASqKM > max_elev_TT_DA) %>%
+# # Sort upstream to downsream
+# arrange(desc(Hydroseq)) %>%
+# # Get cumulative sums of median elevation
+# mutate(csum_length = cumsum(LENGTHKM),
+# cumsum_tt = cumsum(FL_tt_hrs),
+# iter = min(round(sum(FL_tt_hrs) / travt_diff), n()),
+# tt_POI_ID = NA) %>%
+# # Only look to split segments based on travel time longer than 10 km
+# #filter(sum(LENGTHKM) > (split_meters/1000)) %>%
+# filter(iter > 0, n() > 1) %>%
+# # Track original iteration number
+# mutate(orig_iter = iter) %>%
+# ungroup()
+#
+# # For reach iteration
+# tt_POIs <- do.call("rbind", lapply(unique(tt_pois_split$POI_ID), split_tt,
+# tt_pois_split, travt_diff, max_elev_TT_DA)) %>%
+# filter(!tt_POI_ID %in% tmp_POIs$COMID, COMID == tt_POI_ID) %>%
+# mutate(Type_Travel = 1) %>%
+# select(COMID, Type_Travel) %>%
+# unique()
+#
+# tmp_POIs <- POI_creation2(select(tt_POIs, COMID, Type_Travel), nhd, "Travel") %>%
+# addType2(., tmp_POIs, "Travel", nexus = FALSE, bind = TRUE)
+#
+# write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+# }else{
+# tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+# nhd <- read_sf(nav_gpkg, nhd_flowline)
+# }
+#
+# mapview(filter(tmp_POIs, Type_Travel == 1), layer.name = "Travel time break POIs", col.regions = "blue")
+```
+
+```{r Final POIs}
+# # Derive final POI set ----------------------
+# if(needs_layer(nav_gpkg, pois_all_layer)) {
+#
+# unCon_POIs <- filter(nhd, COMID %in% tmp_POIs$COMID, AreaSqKM == 0) %>%
+# distinct()
+#
+# # If any POIs happened to fall on flowlines w/o catchment
+# if (nrow(unCon_POIs) >0){
+# # For confluence POIs falling on Flowlines w/o catchments, derive upstream valid flowline,
+# poi_fix <- DS_poiFix2(tmp_POIs, nhd %>% distinct())
+# new_POIs <- st_compatibalize(poi_fix$new_POIs, tmp_POIs)
+# xWalk <- distinct(poi_fix$xWalk)
+# not_matched <- distinct(poi_fix$unmatched)
+#
+# # POIs that didn't need to be moved
+# tmp_POIs_fixed <- filter(tmp_POIs, nexus == TRUE |
+# !COMID %in% c(poi_fix$xWalk$oldPOI, poi_fix$xWalk$COMID))
+# # bind together
+# final_POIs <- bind_rows(tmp_POIs_fixed, select(poi_fix$new_POIs, -c(oldPOI, to_com)))
+#
+# # if a POI will be a nexus, it can not be terminal.
+# final_POIs <- mutate(final_POIs, Type_Term = ifelse(nexus == 1, NA, Type_Term))
+#
+# # Write out fixes
+# write_sf(new_POIs, temp_gpkg, poi_moved_layer)
+# write_sf(xWalk, temp_gpkg, poi_xwalk_layer)
+# # write out final POIs
+# write_sf(final_POIs, nav_gpkg, pois_all_layer)
+# write_sf(not_matched, temp_gpkg, "on_zeroDA")
+# } else {
+# # If no fixes designate as NA
+# xWalk <- NA
+#
+# tmp_POIs$nexus <- as.integer(tmp_POIs$nexus)
+#
+# # if a POI will be a nexus, it can not be terminal.
+# tmp_POIs <- mutate(tmp_POIs, Type_Term = ifelse(nexus == 1, NA, Type_Term))
+#
+# # write out final POIs
+# write_sf(tmp_POIs, temp_gpkg, pois_all_layer)
+# final_POIs <- tmp_POIs
+# }
+#
+# } else {
+# # Need all three sets for attribution below
+# final_POIs <- read_sf(nav_gpkg, pois_all_layer)
+# new_POIs <- if(!needs_layer(temp_gpkg, poi_moved_layer)) read_sf(temp_gpkg, poi_moved_layer) else (NA)
+# xWalk <- if(!needs_layer(temp_gpkg, poi_xwalk_layer)) read_sf(temp_gpkg, poi_xwalk_layer) else (NA)
+# unCon_POIs <- filter(st_drop_geometry(filter(nhd, minNet == 1)),
+# COMID %in% tmp_POIs$COMID, AreaSqKM == 0)
+# }
+```
+
+```{r Draft nsegment generation}
+# if(needs_layer(temp_gpkg, nsegments_layer)) {
+#
+# # Sort POIs by Levelpath and Hydrosequence in upstream to downstream order
+# seg_POIs <- filter(st_drop_geometry(nhd),
+# COMID %in% tmp_POIs$COMID, COMID %in% filter(nhd, minNet == 1)$COMID)
+# # derive incremental segments from POIs
+# inc_segs <- segment_increment(filter(nhd, minNet == 1), seg_POIs)
+#
+# nhd_Final <- nhd %>%
+# left_join(select(inc_segs, COMID, POI_ID), by = "COMID")
+#
+# # create and write out final dissolved segments
+# nsegments_fin <- segment_creation(nhd_Final, xWalk)
+#
+# nhd_Final <- select(nhd_Final, -POI_ID) %>%
+# left_join(select(st_drop_geometry(nsegments_fin$raw_segs), COMID, POI_ID), by = "COMID")
+# nsegments <- nsegments_fin$diss_segs
+#
+# # Produce the minimal POIs needed to derive the network based on LP and terminal outlets
+# strucFeat <- structPOIsNet(nhd_Final, final_POIs)
+# nhd_Final <- nhd_Final %>%
+# mutate(struct_POI = ifelse(COMID %in% strucFeat$struc_POIs$COMID, 1, 0),
+# struct_Net = ifelse(COMID %in% strucFeat$structnet$COMID, 1, 0))
+#
+# write_sf(nhd_Final, nav_gpkg, nhd_flowline)
+# write_sf(nsegments, temp_gpkg, nsegments_layer)
+# } else {
+# # Read in NHDPlusV2 flowline simple features and filter by vector processing unit (VPU)
+# final_POIs <- read_sf(nav_gpkg, pois_all_layer)
+# nhd_Final <- read_sf(nav_gpkg, nhd_flowline)
+# nsegments <- read_sf(temp_gpkg, nsegments_layer)
+# }
+#
+# # Ensure that all the problem POIs have been taken care of
+# sub <- nhd_Final %>% filter(COMID %in% final_POIs$COMID)
+# print (paste0(nrow(sub[sub$AreaSqKM == 0,]), " POIs on flowlines with no local drainage contributions"))
+# noDA_pois <- filter(final_POIs, COMID %in% filter(sub, AreaSqKM == 0)$COMID)
+# write_sf(noDA_pois, temp_gpkg, "noDA_pois")
+```
+
+# Start from here next time
+```{r POI Collapse}
+# # number POIs
+# final_POIs <- mutate(final_POIs, id = row_number(), moved = NA)
+# write_sf(final_POIs, nav_gpkg, pois_all_layer)
+#
+# # Load data
+# if(collapse){
+#
+# # Move HUC12 to other POIs
+# huc_to_poi <- poi_move(final_POIs, "Type_HUC12", .10, 2.5)
+# move_table <- huc_to_poi$moved_points %>%
+# mutate(move_cat = "HUC12 to other")
+#
+# # Gages to confluences
+# gage_to_conf <- poi_move(huc_to_poi$final_pois, "Type_Gages", .05, 2.5, "Type_Conf") # 47
+# move_table <- gage_to_conf$moved_points %>%
+# mutate(move_cat = "Gage to conf") %>%
+# rbind(move_table)
+#
+# # Gages to waterbody inlets
+# gage_to_wbin <- poi_move(gage_to_conf$final_pois, "Type_Gages", .05, 2.5, "Type_WBIn") # 2
+# move_table <- gage_to_wbin$moved_points %>%
+# mutate(move_cat = "Gage to WBin") %>%
+# rbind(move_table)
+#
+# # Gages to waterbody outlets
+# gage_to_wbout <- poi_move(gage_to_wbin$final_pois, "Type_Gages", .05, 2.5, "Type_WBOut") # 6
+# move_table <- gage_to_wbout$moved_points %>%
+# mutate(move_cat = "Gage to WBout") %>%
+# rbind(move_table)
+#
+# # Streamgage to term outlet
+# gage_to_term <- poi_move(gage_to_wbout$final_pois, "Type_Gages", .05, 2.5, "Type_Term")
+# # NID to waterbody outlet
+# nid_to_wbout <- poi_move(gage_to_wbout$final_pois, "Type_NID", .025, 1, "Type_WBOut")
+#
+# # Sparrow to confluence
+# sparrow_to_conf <- poi_move(gage_to_wbout$final_pois, "Type_Sparrow", .05, 2.5, "Type_Conf")
+# move_table <- sparrow_to_conf$moved_points %>%
+# mutate(move_cat = "Sparrow to Conf") %>%
+# rbind(move_table)
+#
+# # Sparrow to WBoutlet
+# sparrow_to_wbout <- poi_move(sparrow_to_conf$final_pois, "Type_Sparrow", .05, 2.5, "Type_WBOut")
+# # Sparrow to WBInlet
+# sparrow_to_wbin <- poi_move(gage_to_wbout$final_pois, "Type_Sparrow", .05, 2.5, "Type_WBIn")
+#
+# POIs <- sparrow_to_conf$final_pois
+# write_sf(POIs, nav_gpkg, final_poi_layer)
+# write_sf(move_table, temp_gpkg, "pois_collapsed")
+# } else {
+# write_sf(final_POIs, nav_gpkg, final_poi_layer)
+# POIs <- final_POIs
+# }
+```
+
+# ```{r lookup}
+# # Final POI layer
+# POIs <- POIs %>%
+# arrange(COMID) %>%
+# mutate(identifier = row_number())
+#
+# # Unique POI geometry
+# final_POI_geom <- POIs %>%
+# select(identifier) %>%
+# cbind(st_coordinates(.)) %>%
+# group_by(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 = "identifier") %>%
+# select(-identifier)
+#
+# # POI data theme table
+# pois_data_table <- reshape2::melt(st_drop_geometry(select(final_POIs_table,
+# -nexus)),
+# id.vars = c("COMID", "geom_id")) %>%
+# filter(!is.na(value)) %>%
+# group_by(COMID, geom_id) %>%
+# mutate(identifier = cur_group_id()) %>%
+# rename(hy_id = COMID, poi_id = identifier, hl_reference = variable, hl_link = value) %>%
+# distinct()
+#
+# # POI Geometry table
+# poi_geometry_table <- select(final_POIs_table, hy_id = COMID, geom_id) %>%
+# inner_join(distinct(pois_data_table, hy_id, geom_id, poi_id),
+# by = c("hy_id" = "hy_id", "geom_id" = "geom_id")) %>%
+# distinct()
+#
+# write_sf(pois_data_table, nav_gpkg, "poi_data")
+# write_sf(poi_geometry_table, nav_gpkg, "poi_geometry")
+#
+# # Create xwalk table
+# if(needs_layer(nav_gpkg, lookup_table_refactor)) {
+# full_cats <- readRDS(data_paths$fullcats_table)
+#
+# lookup <- dplyr::select(sf::st_drop_geometry(nhd_Final),
+# NHDPlusV2_COMID = COMID,
+# realized_catchmentID = COMID,
+# mainstem = LevelPathI) %>%
+# dplyr::mutate(realized_catchmentID = ifelse(realized_catchmentID %in% full_cats$FEATUREID,
+# realized_catchmentID, NA)) %>%
+# left_join(select(st_drop_geometry(poi_geometry_table), hy_id, poi_geom_id = geom_id),
+# by = c("NHDPlusV2_COMID" = "hy_id"))
+#
+# sf::write_sf(lookup, nav_gpkg, lookup_table_refactor)
+# }
+#
+# ```
diff --git a/workspace/06-2_aggregate_cats.Rmd b/workspace/06-2_aggregate_cats.Rmd
index ba68264..1b1584c 100644
--- a/workspace/06-2_aggregate_cats.Rmd
+++ b/workspace/06-2_aggregate_cats.Rmd
@@ -44,15 +44,20 @@ outlets_POI <- read_sf(out_refac_gpkg, outlets_layer) %>%
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)
+reconciled <- st_transform(read_sf(out_refac_gpkg, reconciled_layer), crs)
+refactored <- st_drop_geometry(read_sf(out_refac_gpkg, refactored_layer)) %>%
+ select(COMID, Hydroseq) %>%
+ mutate(COMID = as.integer(COMID))
+
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 %>%
+ inner_join(refactored, by = c("NHDPlusV2_COMID" = "COMID")) %>%
group_by(reconciled_ID) %>%
- filter(TotDASqKM == max(TotDASqKM))
+ filter(Hydroseq == min(Hydroseq))
```
```{r agg_catchments I}
@@ -79,12 +84,14 @@ if(needs_layer(out_agg_gpkg, agg_cats_layer)){
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)) %>%
- select(-toID)
}
+ outlets <- outlets %>%
+ left_join(select(st_drop_geometry(reconciled), ID, toID), by = "ID") %>%
+ mutate(type = ifelse(is.na(toID), "terminal", "outlet")) %>%
+ mutate(type = ifelse(type == "terminal" & !is.na(toID), "outlet", type)) %>%
+ select(-toID)
+
outlets <- unique(outlets)
agg_cats <- aggregate_catchments(flowpath = reconciled,
divide = divides,
--
GitLab