diff --git a/workspace/00_get_data.Rmd b/workspace/00_get_data.Rmd
index 10cfa0f2746d51726e6162115e8e07343d52f5f7..ab03de56e22722f17bca79b65e0e0101ea6f434d 100644
--- a/workspace/00_get_data.Rmd
+++ b/workspace/00_get_data.Rmd
@@ -1,9 +1,9 @@
---
-title: "GFv2 Get Data"
+title: "Get Hydro Location and Network Data"
output: html_document
---
-This notebook pulls data from a number of sources and populates the GFv2 data
+This notebook pulls data from a number of sources and populates the data
directory. Any new data requirements should be added as code chunks here.
Each code chunk should create a path to the file you want to use in a process
@@ -41,7 +41,9 @@ initialize_sciencebase_session(username = Sys.getenv("sb_user"))
# Enable mapview rendering if desired
mapview <- FALSE
```
-
+HUC12 (Hydrologic Unit Code, Level 12) outlets derived from the Watershed
+Boundary Dataset indexed to the reference fabricform the baseline and extent of
+national modeling fabrics.
```{r HUC12 outlets}
# Blodgett, D.L., 2022, Mainstem Rivers of the Conterminous United States:
# U.S. Geological Survey data release, https://doi.org/10.5066/P9BTKP3T.
@@ -59,6 +61,13 @@ out_list <- c(out_list, list(hu12_points_path = hu12_points_path))
if(mapview)(mapview(read_sf(out_list$hu12_points_path)))
```
+SWIM (Streamgage Watershed InforMation(SWIM) includes locations for 12,422 USGS
+streamgages as indexed along the network of streams (flowlines) in NHDPlus
+Version 2.1 (NHDPlus v2, Moore and Dewald, 2016). The dataset is one of two
+datasets developed for the Streamgage Watershed InforMation (SWIM) project. This
+dataset, which is referred to as “SWIM streamgage locations,†was created in
+support of the second dataset of basin characteristics and disturbance indexes.
+
```{r SWIM}
# Hayes, L., Chase, K.J., Wieczorek, M.E., and Jackson, S.E., 2021,
# USGS streamgages in the conterminous United States indexed to NHDPlus v2.1
@@ -78,6 +87,11 @@ out_list <- c(out_list, list(SWIM_points_path = SWIM_points_path))
if(mapview)(mapview(read_sf(out_list$SWIM_points_path)))
```
+Sites associated with Work by the U.S. Geological Survey (USGS) to estimate
+the amount of water that is withdrawn and consumed by thermoelectric power
+plants (Diehl and others, 2013; Diehl and Harris, 2014; Harris and Diehl, 2019
+Galanter and othes, 2023).
+
```{r Thermoelectric Facilities}
# Harris, Melissa A. and Diehl, Timothy H., 2017. A Comparison of Three
# Federal Datasets for Thermoelectric Water Withdrawals in the United States
@@ -123,6 +137,12 @@ out_list <- c(out_list, list(TE_points_path = TE_points_path))
if(mapview)(mapview(read_sf(out_list$TE_points_path)))
```
+Network locations made to improve the routing capabilities
+and ancillary hydrologic attributes of NHDPlusV2 to support modeling and other
+hydrologic analyses. The resulting enhanced network is named E2NHDPlusV2_us.
+This includes the network locations associated with some diversions and
+water use withdrawals.
+
```{r e2nhd supplemental data - USGS}
# Schwarz, G.E., 2019, E2NHDPlusV2_us: Database of Ancillary Hydrologic
# Attributes and Modified Routing for NHDPlus Version 2.1 Flowlines: U.S.
@@ -146,6 +166,11 @@ if(!file.exists(USGS_IT_path)) {
out_list <- c(out_list, list(USGS_IT_path = USGS_IT_path))
```
+Two datasets relate hydro location information from the National Inventory of
+Dams to the NHDPlus network. One effort is related to the SPARROW work
+(Wieczorek and others, 2018), the other related to work quantifying impacts on
+natural flow (Wieczorek and others, 2021).
+
```{r National Inventory of Dams}
# Wieczorek, M.E., Jackson, S.E., and Schwarz, G.E., 2018, Select Attributes
@@ -186,6 +211,13 @@ out_list <- c(out_list, list(NID_points_path = NID_points_path))
if(mapview)(mapview(read_sf(out_list$NID_points_path)))
```
+This next section retrieves NHDPlus datasets related to national modeling
+efforts. These include:
+1. National Geodatbase
+2. Hawaii, Puerto Rico, and Islands Geodatabase
+3. GageLoc file of streamgages indexed to NHDPlusv2 flowlines
+4. NHDPlusv2 catchment - HUC12 crosswalk.
+
```{r NHDPlusV2}
# NHDPlus Seamless National Data - pulled from NHDPlus national data server;
# post-processed to RDS files by NHDPlusTools GageLoc - Gages snapped to
@@ -238,11 +270,9 @@ if(!file.exists(nhdplus_dir)) {
, " -o", data_dir), ignore.stdout = T)}
)
- suppressWarnings(staged_nhd <- stage_national_data(nhdplus_data = nhdplus_gdb))
-
x <- tryCatch(
download_nhdplusv2(islands_dir,
- url = paste0("https://edap-ow-data-commons.s3.amazonaws.com/NHDPlusV21/",
+ url = paste0("https://dmap-data-commons-ow.s3.amazonaws.com/NHDPlusV21/",
"Data/NationalData/NHDPlusV21_NationalData_Seamless",
"_Geodatabase_HI_PR_VI_PI_03.7z")),
@@ -252,8 +282,6 @@ if(!file.exists(nhdplus_dir)) {
, " -o", islands_dir), ignore.stdout = T)}
)
- suppressWarnings(staged_nhdplus_islands <- stage_national_data(nhdplus_data = islands_gdb))
-
rm(gLz, xWalk)
HUC12 <- read_sf(nhdplus_gdb, layer = "HUC12") %>%
@@ -285,6 +313,15 @@ out_list <- c(out_list, list(nhdplus_dir = nhdplus_dir,
nhdplus_rpu = rpu))
```
+Reference catchments and flowlines are hydrographic products that are derived
+from the USGS National Hydrologic Geospatial Fabric and the National Oceanic
+and Atmospheric Administration. The Reference Flowlines include modifications
+from the National Water Model and e2nhd networks integrated into NHDPlusv2.1,
+and the reference catchments are geometrically-simplified to POLYGON geometry to
+improve rasterre-gridding efficiency, and have a large number of DEM artifacts
+removed.
+
+
```{r Reference Fabric}
# Reference Fabric flowpaths and catchments derived by Mike Johnson (NOAA)
@@ -305,7 +342,7 @@ if(!file.exists(ref_cat)){
for (vpu in possible_vpu){
vpu_gpkg <- paste0(vpu, "_reference_features.gpkg")
check_auth()
- sbtools::item_file_download("61295190d34e40dd9c06bcd7", names = vpu_gpkg,
+ sbtools::item_file_download("6317a72cd34e36012efa4e8a", names = vpu_gpkg,
destinations = file.path(ref_fab_path, vpu_gpkg))
}
check_auth()
@@ -327,6 +364,9 @@ out_list <- c(out_list, list(ref_fab_path = ref_fab_path, ref_cat = ref_cat,
ref_fl = ref_fl, nwm_fl = nwm_fl))
```
+NHDPlus Value-added attributes used to characterize the network and get
+information such as velocity and mean streamflow.
+
```{r NHDPlusV2 VAA}
# Waterbodies - derived after downloading and post-processing
# NHDPlus Seamless National Geodatabase
@@ -335,7 +375,7 @@ VAA_fst_path <- file.path(nhdplus_dir, "nhdplusVAA.fst")
VAA_rds_path <- file.path(nhdplus_dir, "nhd_vaa.rds")
#TODO pull attributes needed for threshold POIs
-if(!file.exists(VAA_fst_path)) {
+if(!file.exists(VAA_rds_path)) {
message("downloading NHDPlusVAA...")
get_vaa(path = nhdplus_dir)
@@ -345,9 +385,11 @@ if(!file.exists(VAA_fst_path)) {
saveRDS(nhdpv2_vaa, file.path(nhdplus_dir, "nhd_vaa.rds"))
}
-out_list <- c(out_list, list(VAA_fst = VAA_path, VAA_rds = VAA_rds_path))
+out_list <- c(out_list, list(VAA_fst = VAA_fst_path, VAA_rds = VAA_rds_path))
```
+NHDPlus Waterbody and Area Polygons converted to an RDS file for easier
+loading within R.
```{r NHDPlusV2 Waterbodies}
# Waterbodies - derived after downloading and post-processing
@@ -377,6 +419,10 @@ if(!file.exists(waterbodies_path)) {
out_list <- c(out_list, list(waterbodies_path = waterbodies_path))
```
+Formatting a full list of network and non-network catchments for the NHDPlus
+domains. This more easily tracks catchments were are off and on the network
+when aggregating at points of interest.
+
```{r full cats}
# Modification to NHDPlus catchments
fullcat_path <- file.path(nhdplus_dir, "nhdcat_full.rds")
@@ -385,12 +431,10 @@ islandcat_path <- file.path(islands_dir, "nhdcat_full.rds")
# Create full cat dataset
if(!file.exists(fullcat_path)){
- cat_tab <- cat_rpu(staged_nhd$catchment, nhdplus_gdb)
+ cat_tab <- cat_rpu(out_list$ref_cat, nhdplus_gdb)
saveRDS(cat_tab, fullcat_path)
- island_cats <- file.path(islands_dir,
- "NHDPlusNationalData/nhdplus_catchment.rds")
- island_tab <- cat_rpu(island_cats, islands_gdb)
+ island_tab <- cat_rpu(out_list$islands_gdb, islands_gdb)
saveRDS(island_tab, islandcat_path)
}
@@ -398,6 +442,8 @@ out_list <- c(out_list, list(fullcats_table = fullcat_path))
out_list <- c(out_list, list(islandcats_table = islandcat_path))
```
+Download NHDPlusV2 FDR and FAC grids for refactoring and catcment splitting.
+
```{r NHDPlusV2 FDR_FAC}
# NHDPlus FDR/FAC grids available by raster processing unit
@@ -424,6 +470,9 @@ out$fac <- as.list(setNames(fac, paste0("rpu_", rpu)))
out_list<- c(out_list, out)
```
+Download NHDPlusV2 elevation grids for headwater extensions and splitting
+catchments into left and right banks.
+
```{r NHDPlusV2 elev}
# NHDPlus FDR/FAC grids available by raster processing unit
@@ -448,6 +497,8 @@ out$elev_cm <- as.list(setNames(elev_cm, paste0("rpu_", rpu)))
out_list<- c(out_list, out)
```
+Download the current WBD snapshot.
+
```{r WBD}
# Snapshot of National WBD
@@ -455,8 +506,6 @@ wbd_dir <- file.path(data_dir, "wbd")
wbd_file <- "WBD_National_GDB"
if(!dir.exists(wbd_dir)) {
- if(is.null(sbtools::current_session()))
- sb <- authenticate_sb()
dir.create(wbd_dir, recursive = TRUE)
wbd <- download.file(
@@ -486,6 +535,10 @@ out_rds <- list(latest_wbd_rds = out$latest_wbd)
out_list <- c(out_list, out_rds)
```
+
+Merrit Topographic and Hydrographic data for deriving GIS Features of the
+National Hydrologic Modeling, Alaska Domain
+
```{r MERIT HydroDEM}
# MERIT HydroDEM - used for AK Geospatial Fabric, and potentially
# Mexico portion of R13
@@ -534,9 +587,7 @@ if(!dir.exists(merit_dir)) {
merit_fac = get_sbfile(file.path(merit_dir,
"ak_merit_fac.zip"),
"64ff628ed34ed30c2057b430"))
-}
-
-out <- list(merit_catchments = file.path(merit_dir,
+ out <- list(merit_catchments = file.path(merit_dir,
"cat_pfaf_78_81_82_MERIT_Hydro_v07_Basins_v01.shp"),
merit_rivers = file.path(merit_dir,
"riv_pfaf_78_81_82_MERIT_Hydro_v07_Basins_v01.shp"),
@@ -546,17 +597,24 @@ out <- list(merit_catchments = file.path(merit_dir,
merit_fac = file.path(merit_dir, "ak_merit_fac.tif"))
-out <- list(aster_dem = get_sbfile(file.path(merit_dir, "dem.zip"), "5fbbc6b6d34eb413d5e21378"),
- merit_dem = get_sbfile(file.path(merit_dir,
- "ak_merit_dem.zip"), "5fc51e65d34e4b9faad8877b"),
- merit_fdr = get_sbfile(file.path(merit_dir,
- "ak_merit_fdr.zip"), "5fc51e65d34e4b9faad8877b"),
- merit_fac = get_sbfile(file.path(merit_dir,
- "ak_merit_fac.zip"), "5fc51e65d34e4b9faad8877b"))
+ out <- list(aster_dem = get_sbfile(file.path(merit_dir, "dem.zip"), "5fbbc6b6d34eb413d5e21378"),
+ merit_dem = get_sbfile(file.path(merit_dir,
+ "ak_merit_dem.zip"), "5fc51e65d34e4b9faad8877b"),
+ merit_fdr = get_sbfile(file.path(merit_dir,
+ "ak_merit_fdr.zip"), "5fc51e65d34e4b9faad8877b"),
+ merit_fac = get_sbfile(file.path(merit_dir,
+ "ak_merit_fac.zip"), "5fc51e65d34e4b9faad8877b"))
+
+}
+
+
out_list <- c(out_list, out)
```
+Source data for deriving GIS Featurs of the National Hydrologic Modeling,
+Alaska Domain
+
```{r AK GF Source data}
# Bock, A.R., Rosa, S.N., McDonald, R.R., Wieczorek, M.E., Santiago, M.,
# Blodgett, D.L., and Norton, P.A., 2024, Geospatial Fabric for National
@@ -580,6 +638,9 @@ if(!dir.exists(AK_dir)) {
}
```
+Source data for deriving GIS Featurs of the National Hydrologic Modeling,
+Hawaii Domain
+
```{r HI GF Source data}
# Bock, A.R., Rosa, S.N., McDonald, R.R., Wieczorek, M.E., Santiago, M.,
# Blodgett, D.L., and Norton, P.A., 2024, Geospatial Fabric for National
@@ -602,6 +663,7 @@ if(!file.exists(file.path(islands_dir, "hi.gpkg"))) {
out_list <- c(out_list, out_hi)
```
+National Water Model Network Topology
```{r nwm_topology}
nwm_targz_url <-
@@ -641,6 +703,8 @@ out_list <- c(out_list, out)
```
+e2nhd Network Attributes
+
```{r nhdplus_attributes}
# Blodgett, D.L., 2023, Updated CONUS river network attributes based on the
# E2NHDPlusV2 and NWMv2.1 networks (ver. 2.0, February 2023): U.S. Geological
@@ -659,6 +723,9 @@ if(!file.exists(out$new_nhdp_atts)) {
out_list <- c(out_list, out)
```
+GIS Features of the Geospatial Fabric for National Hydrologic Modeling,
+version 1.1, Transboundary Geospatial Fabric
+
```{r GFv1.1}
# Bock, A.E, Santiago,M., Wieczorek, M.E., Foks, S.S., Norton, P.A., and
# Lombard, M.A., 2020, Geospatial Fabric for National Hydrologic Modeling,
@@ -702,6 +769,8 @@ out_list <- c(out_list, out)
if(mapview)(mapview(readRDS(out_list$GFv11_gages_lyr)))
```
+GAGESII dataset
+
```{r Gages_II}
# https://doi.org/10.3133/70046617
@@ -726,6 +795,9 @@ out_list <- c(out_list, g2_out)
if(mapview)(mapview(read_sf(out_list$gagesii_lyr)))
```
+HILARRI dataset of Network-indexed Hydropower structures, reservoirs, and
+locations
+
```{r HILARRI}
# Carly H. Hansen and Paul G. Matson. 2023. Hydropower Infrastructure - LAkes,
# Reservoirs, and RIvers (HILARRI), # Version 2. HydroSource. Oak Ridge
@@ -746,7 +818,7 @@ if(!file.exists(file.path(hilarri_dir, "HILARRI_v2.csv"))){
dest = file.path(hilarri_dir, basename(hilarri_url)),
mode = "wb")
- unzip(file.path(hillari_dir, basename(hilarri_url)), exdir = hilarri_dir)
+ unzip(file.path(hilarri_dir, basename(hilarri_url)), exdir = hilarri_dir)
}
out_list <- c(out_list, hilarri_out)
@@ -757,6 +829,8 @@ if(mapview){
```
+ResOpsUS dataset and indexed locations
+
```{r Reservoir datasets}
# ResOpsUS
# Steyaert, J.C., Condon, L.E., W.D. Turner, S. et al. ResOpsUS, a dataset of
@@ -818,6 +892,8 @@ tab_out <- list(resops_NID_CW = resops_to_nid_path)
out_list <- c(out_list, tab_out)
```
+All Hydro-linked Network Data Index (NLDI) datasets
+
```{r nldi}
# NLDI feature data sources
# https://www.sciencebase.gov/catalog/item/60c7b895d34e86b9389b2a6c
@@ -835,7 +911,7 @@ for (fname in nldi_list$fname){
print(fname)
floc <- file.path(nldi_dir, fname)
if(!file.exists(floc)){
- check_auth()
+ #check_auth()
sbtools::item_file_download("60c7b895d34e86b9389b2a6c",
names = fname,
destinations = floc)
diff --git a/workspace/01_Gage_Selection.Rmd b/workspace/01_Gage_Selection.Rmd
index 8230103bf5826058fc58587962607a42f7a81ff3..c39ee2d293d9c984177529ab2416bb343301012b 100644
--- a/workspace/01_Gage_Selection.Rmd
+++ b/workspace/01_Gage_Selection.Rmd
@@ -347,7 +347,7 @@ gage_final <- data.table::rbindlist(list(gagesIII_final, gageloc_final,
# which gagesii_ref gages are missing from final set
missing_gagesii_ref <- filter(gages_ii_ref, !id %in% gage_final$site_no) %>%
rename(site_no = id, comid = nhdpv2_COMID, reachcode = nhdpv2_REACHCODE,
- reach_meas = nhdpv2_REACH_measure) #%>%
+ reach_meas = nhdpv2_REACH_measure) %>%
inner_join(gageloc_db_full, by = "site_no") %>%
select(-c(agency_cd, site_tp_cd, coord_acy_, huc_cd, parm_cd, stat_cd, ts_id, missing_da)) %>%
st_compatibalize(gage_final)
diff --git a/workspace/01_NHD_prep.Rmd b/workspace/01_NHD_prep.Rmd
index 6841eaef02a0f1b7706725be4b5c6ed7ea192d8d..c21d5108cc81e0e55c6daebb16406d7887e12ae3 100644
--- a/workspace/01_NHD_prep.Rmd
+++ b/workspace/01_NHD_prep.Rmd
@@ -22,6 +22,7 @@ source("R/NHD_navigate.R")
# Load Configuration of environment
source("R/config.R")
+source("R/user_vars.R")
if(exists("rpu_code")){rm(rpu_code)}
@@ -66,11 +67,11 @@ nhd <- fline %>%
# Filter and write dendritic/non-coastal subset to gpkg
# This will be iterated over to produce the final network after POIs identified
-zero_order <- filter(nhd, TerminalFl == 1 & TotDASqKM < min_da_km)
+zero_order <- filter(nhd, TerminalFl == 1 & TotDASqKM < min_da_km_gages)
non_dend <- unique(unlist(lapply(zero_order$COMID, NetworkNav, nhdDF = st_drop_geometry(nhd))))
nhd <- nhd %>%
mutate(dend = ifelse(!COMID %in% non_dend, 1, 0),
- poi = ifelse(!COMID %in% non_dend & TotDASqKM >= min_da_km, 1, 0))
+ poi = ifelse(!COMID %in% non_dend & TotDASqKM >= min_da_km_gages, 1, 0))
cat_network <- sf::st_drop_geometry(nhd)
names(cat_network) <- tolower(names(cat_network))
diff --git a/workspace/02_NHD_navigate.Rmd b/workspace/02_NHD_navigate.Rmd
index 9bf3eaa3d6c78ff8d16a9c3b710942e47af57dac..7dfde7228a9466264b3ba4e62817313552f2d0fd 100644
--- a/workspace/02_NHD_navigate.Rmd
+++ b/workspace/02_NHD_navigate.Rmd
@@ -1,12 +1,14 @@
---
-title: "NHD Navigate"
+title: "MAPNAT POI Geneation"
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.
+This notebook generates a network of Points of Interest harvested from Public
+source that are used to generate a hydrofabric for use within the USGS
+MAPNAT Project.
+
+
```{r setup_rmd, echo=FALSE, cache=FALSE}
knitr::opts_chunk$set(
@@ -22,18 +24,21 @@ knitr::opts_chunk$set(
source("R/utils.R")
source("R/NHD_navigate.R")
source("R/hyrefactor_funs.R")
+#source("R/test.R")
# increase timeout for data downloads
options(timeout=600)
# Load Configuration of environment
source("R/config.R")
+source("R/user_vars.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"))
+#atts <- read_fst(data_paths$VAA_fst)
+atts_rds <- readRDS(data_paths$VAA_rds)
```
```{r huc12 POIs}
@@ -41,7 +46,10 @@ atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.r
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)
+ if("POI_ID" %in% colnames(nhd)){
+ try(nhd <- select(nhd, -c(wb_id, POI_ID)), silent = TRUE) %>%
+ distinct()
+ }
# Some NHDPlus VPUs include HUCs from other VPUs
if(vpu_codes == "02"){
@@ -83,6 +91,7 @@ if(!"Type_Gages" %in% names(tmp_POIs)) {
streamgages_VPU <- gages %>%
rename(COMID = comid) %>%
filter(COMID %in% nhd$COMID) %>%
+ inner_join(st_drop_geometry(st_drop_geometry(nhd), COMID, Divergence), by = "COMID") %>%
switchDiv(., nhd)
streamgages <- streamgages_VPU %>%
@@ -98,7 +107,8 @@ if(!"Type_Gages" %in% names(tmp_POIs)) {
# As a fail-safe, write out list of gages not assigned a POI
if(nrow(filter(streamgages_VPU, !site_no %in% tmp_POIs$Type_Gages)) > 0) {
- write_sf(filter(streamgages_VPU, !site_no %in% tmp_POIs$Type_Gages),
+ write_sf(filter(streamgages_VPU, !site_no %in% tmp_POIs$Type_Gages) %>%
+ rename(nhd_reachcode = REACHCODE),
temp_gpkg, "unassigned_gages")
# Start documenting gages that are dropped out; these gages have no mean daily Q
@@ -117,7 +127,7 @@ if(!"Type_Gages" %in% names(tmp_POIs)) {
events <- read_sf(temp_gpkg, split_layer)
}
-mapview(filter(tmp_POIs, !is.na(Type_Gages)), layer.name = "Streamgage POIs", col.regions = "blue")
+#mapview(filter(tmp_POIs, !is.na(Type_Gages)), layer.name = "Streamgage POIs", col.regions = "blue")
```
```{r TE POIs}
@@ -130,159 +140,410 @@ if(!"Type_TE" %in% names(tmp_POIs)) {
}
# 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)) %>%
+ TE_COMIDs <- readxl::read_xlsx(file.path(data_paths$TE_points_path,
+ "nhdv1_2_comids_all_plants_all_years.xlsx")) %>%
+ rename(COMID = comid) %>%
+ #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_codes, 1, 2), ".*"), .data$VPUID), COMID > 0) %>%
switchDiv(., nhd) %>%
group_by(COMID) %>%
- summarize(EIA_PLANT = paste0(unique(EIA_PLANT_), collapse = " "), count = n()) %>%
+ summarize(eia_id = paste0(unique(eia_id), collapse = " "), count = n()) %>%
ungroup()
- # Derive TE POIs
- tmp_POIs <- POI_creation(st_drop_geometry(TE_COMIDs), filter(nhd, poi == 1), "TE") %>%
- addType(., tmp_POIs, "TE", nexus = FALSE)
-
- # 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")
+ if(nrow(TE_COMIDs) > 0){
+ # Derive TE POIs
+ tmp_POIs <- POI_creation(st_drop_geometry(TE_COMIDs), filter(nhd, poi == 1), "TE") %>%
+ addType(., tmp_POIs, "TE", nexus = 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 {
+ print ("no TEs in VPU")
}
- # 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")
+#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)) {
-
- # 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)) %>%
+```{r WB mods}
+
+# Waterbodies sourced from NHD waterbody layer for entire VPU
+WBs_VPU_all <- filter(readRDS(data_paths$waterbodies_path),
+ COMID %in% nhd$WBAREACOMI) %>%
+ filter(FTYPE != "SwampMarsh") %>%
+ mutate(FTYPE = as.character(FTYPE),
+ source = "NHDv2WB",
+ wb_id = ifelse(GNIS_ID == " ", COMID, GNIS_ID)) %>%
+ st_transform(geo_crs) %>%
+ st_make_valid()
+
+saveRDS(WBs_VPU_all, data_paths$waterbodies_path)
+
+if(needs_layer(temp_gpkg, "WBs_layer_orig")){
+ sf_use_s2(FALSE)
+ ref_WB <- WBs_VPU_all %>%
+ group_by(wb_id, GNIS_NAME) %>%
+ sfheaders::sf_remove_holes() %>%
+ summarize(do_union = T, member_comid = paste(COMID, collapse = ",")) %>%
+ st_make_valid() %>%
+ ungroup() %>%
+ mutate(area_sqkm = as.numeric(st_area(.)) / 1000000) %>%
mutate(source = "NHDv2WB") %>%
- st_transform(crs)
-
- ref_WB <- read_sf(nav_gpkg, nhd_waterbody) %>%
- filter(COMID %in% nhd$WBAREACOMI) %>%
- mutate(source = "Ref_WB")
-
- # 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)
-
- 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()
+ st_cast("MULTIPOLYGON")
+ #sf_use_s2(TRUE)
+
+ write_sf(ref_WB, temp_gpkg, "WBs_layer_orig")
+} else {
+ ref_WB <- read_sf(temp_gpkg, "WBs_layer_orig")
+}
+```
+
+```{r resops}
+# # 1: Many based on original GNIS_ID
+if(!"Type_resops" %in% names(tmp_POIs)){
+
+ # Unnest wb_poi_lst
+ wb_table <- st_drop_geometry(ref_WB) %>%
+ dplyr::mutate(member_comid = strsplit(member_comid, ",")) %>%
+ tidyr::unnest(cols = member_comid) %>%
+ mutate(member_comid = as.character(member_comid)) %>%
+ distinct()
+
+ # ResOpsUS locations with attributes from the VPU NHDv2 wb set
+ resops_wb_df <- read.csv("data/reservoir_data/resops_crosswalk_AB.csv") %>%
+ # Subset to VPU, only one DAMID per waterbody
+ filter(flowlcomid %in% nhd$COMID |
+ wbareacomi %in% wb_table$member_comid) %>%
+ dplyr::select(grand_id, nid_source_featureid, source = comi_srclyr,
+ flowlcomid, wbareacomi, hr_permid, onoffnet) %>%
+ mutate(wbareacomi = as.character(wbareacomi)) %>%
+ # Link to waterbody table
+ left_join(distinct(wb_table, member_comid, wb_id),
+ by = c("wbareacomi" = "member_comid")) %>%
+ left_join(select(st_drop_geometry(WBs_VPU_all), wbareacomi = COMID,
+ wb_id2 = wb_id) %>%
+ mutate(wbareacomi = as.character(wbareacomi)),
+ by = "wbareacomi") %>%
+ mutate(wb_id = ifelse(is.na(wb_id), wb_id2, wb_id),
+ wbareacomi = ifelse(wbareacomi == -2,
+ hr_permid, wbareacomi)) %>%
+ select(-c(wb_id2))
+
+ # POIs with waterbody IDs in reference waterbodies
+ resops_wb_pois <- filter(ref_WB, wb_id %in% resops_wb_df$wb_id) %>%
+ inner_join(select(st_drop_geometry(resops_wb_df), grand_id,
+ NID_ID = nid_source_featureid,
+ resops_flowlcomid = flowlcomid, wb_id),
+ by = "wb_id") %>%
+ #mutate(source = "ref_WB") %>%
+ group_by(wb_id) %>%
+ filter(n() == 1) %>%
+ st_as_sf()
+
+ # Add ResOPsUS locations to waterbody list with attributed NID and resops data
+ if(nrow(resops_wb_pois) > 0){
+ wb_poi_lst_filtered <- filter(ref_WB, !wb_id %in% resops_wb_pois$wb_id,
+ !is.na(wb_id))
+
+ wb_poi_lst <- data.table::rbindlist(list(wb_poi_lst_filtered, resops_wb_pois),
+ fill = TRUE) %>%
+ mutate(accounted = 0) %>%
+ st_as_sf()
- # 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)
+ # Reach resopsus
+ reach_resops <- filter(resops_wb_df, !wb_id %in% wb_poi_lst$wb_id,
+ !source %in% c("NHDAREA", "HR ID AVAILABLE") |
+ onoffnet == 0) %>%
+ mutate(Type_WBOut = NA) %>%
+ select(COMID = flowlcomid, resops = grand_id, Type_WBOut)
- # Attribute flowlines that are in waterbodies
- nhd <- mutate(nhd, WB = ifelse(WBAREACOMI > 0 & WBAREACOMI %in% WBs_VPU$COMID, 1, 0))
+ # Existing reservoir-waterbody outlet POIs
+ exist_POIs_WB <- filter(resops_wb_df, flowlcomid %in% tmp_POIs$COMID) %>%
+ mutate(wb_id = ifelse(source == "NHDAREA", wbareacomi, wb_id)) %>%
+ filter(!is.na(wb_id)) %>%
+ select(COMID = flowlcomid, resops = grand_id, Type_WBOut = wb_id)
- wb_layers <- wbout_POI_creaton(filter(nhd, WB == 1), WBs_VPU, data_paths, crs)
- if(!is.na(wb_layers$events) > 0) {events <- rbind(events, wb_layers$events)}
- WBs_VPU <- wb_layers$WBs
+ # Resops POIs
+ resops_pois <- rbind(reach_resops, exist_POIs_WB)
- wb_out_col <- wb_poi_collapse(wb_layers$POIs, filter(nhd, WB == 1), events)
- tmp_POIs <- wb_out_col$POIs
- if(!all(is.na(wb_out_col$events_ret))) {
- write_sf(wb_out_col$events_ret, temp_gpkg, split_layer)
+ # Resops defined by reach
+ if(nrow(resops_pois) > 0){
+ # Resops POIs with no reservoirs, defined by reach
+ tmp_POIs <- POI_creation(resops_pois, filter(nhd, poi == 1), "resops") %>%
+ addType(., tmp_POIs, "resops", nexus = TRUE)
+
+ # Add waterbody attribute
+ tmp_POIs <- tmp_POIs %>%
+ left_join(select(resops_pois, COMID, Type_WBOut),
+ by = "COMID") %>%
+ mutate(Type_WBOut = ifelse(!nexus,
+ Type_WBOut, NA))
+
+ # Resops with reservoirs
+ resops_wb_df <- resops_wb_df %>%
+ mutate(accounted = ifelse(grand_id %in% tmp_POIs$Type_resops, 1, 0),
+ source = ifelse(grand_id %in% reach_resops$resops, "REACH", source))
}
+
+ write_sf(wb_poi_lst, temp_gpkg, WBs_layer)
+ write_sf(resops_wb_df, temp_gpkg, "wb_resops")
+ write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+} else {
+ tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+ wb_poi_lst <- read_sf(temp_gpkg, WBs_layer)
+ resops_wb_df <- read_sf(temp_gpkg, "wb_resops")
+}
+```
+
+```{r Hilari}
+if(!"Type_hilarri" %in% names(tmp_POIs)){
- if(!all(is.na(wb_layers$events))){
- wb_events <- select(wb_layers$events, -c(id, offset)) %>%
- rename(POI_identifier = WBAREACOMI)
-
- # Write out events and outlets
- if(!needs_layer(temp_gpkg, split_layer)){
- events <- read_sf(temp_gpkg, split_layer) %>%
- rbind(st_compatibalize(wb_events,.)) %>%
- unique()
- write_sf(events, temp_gpkg, split_layer)
- } else {
- write_sf(wb_events, nav_gpkg, split_layer)
- }
+ # 1: Many based on original GNIS_ID
+ wb_table <- st_drop_geometry(wb_poi_lst) %>%
+ dplyr::mutate(member_comid = strsplit(member_comid, ",")) %>%
+ tidyr::unnest(cols = member_comid) %>%
+ mutate(member_comid = as.integer(member_comid)) %>%
+ distinct()
+
+ # Hilarri POI Points
+ hilarri_points <- read.csv("data/HILARRI/HILARRI_v2.csv",
+ colClasses = "character") %>%
+ mutate(nhdwbcomid = as.integer(nhdwbcomid)) %>%
+ filter(!dataset %in% c('Power plant only; no reservoir or inventoried dam'),
+ nhdv2comid %in% nhd$COMID) %>%
+ left_join(select(wb_table, member_comid, wb_id),
+ by = c("nhdwbcomid" = "member_comid"))
+
+ # Waterbodies linked to hilarri information
+ hil_wb_pois_rf <- wb_poi_lst %>%
+ inner_join(select(st_drop_geometry(hilarri_points), hilarriid, nhdwbcomid,
+ nidid_hill = nidid, wb_id),
+ by = "wb_id") %>%
+ group_by(wb_id) %>%
+ mutate(hilarriid = paste0(hilarriid, collapse = ","),
+ nidid_hill = paste0(nidid_hill, collapse = ",")) %>%
+ ungroup() %>%
+ distinct() %>%
+ st_as_sf() %>%
+ st_compatibalize(wb_poi_lst)
+
+ # Add ResOPsUS locations to waterbody list
+ if(nrow(hil_wb_pois_rf) > 0){
+ wb_poi_lst <- wb_poi_lst %>%
+ #filter(!is.na(resops_FlowLcomid)) %>%
+ left_join(select(st_drop_geometry(hil_wb_pois_rf), wb_id, hilarriid,
+ nidid_hill), by = "wb_id") %>%
+ mutate(NID_ID = ifelse(is.na(NID_ID), nidid_hill, NID_ID)) %>%
+ select(-nidid_hill) %>%
+ distinct()
}
- ref_WB <- filter(ref_WB, COMID %in% WBs_VPU$COMID)
+ # Reach POIs
+ reach_pois <- filter(hilarri_points, !hilarriid %in% wb_poi_lst$hilarriid) %>%
+ select(COMID = nhdv2comid, hilarriid) %>%
+ mutate(COMID = as.integer(COMID)) %>%
+ group_by(COMID) %>%
+ filter(row_number() == 1) %>%
+ st_drop_geometry()
- WBs_VPU <- filter(WBs_VPU, !COMID %in% ref_WB$COMID) %>%
- group_by(GNIS_ID, GNIS_NAME, COMID, FTYPE, source) %>%
- summarize(do_union = T) %>%
- ungroup() %>%
- st_make_valid() %>%
- sfheaders::sf_remove_holes(.) %>%
- mutate(member_comid = NA,
- area_sqkm = as.numeric(st_area(.))/1000000) %>%
- st_compatibalize(., ref_WB)
+ # Make POIs for reach POIs
+ tmp_POIs <- POI_creation(reach_pois, filter(nhd, poi == 1), "hilarri") %>%
+ addType(., tmp_POIs, "hilarri", nexus = TRUE)
- ref_WB <- rbind(ref_WB, WBs_VPU)
+ # Write out geopackage layer representing POIs for given theme
+ write_sf(wb_poi_lst, temp_gpkg, WBs_layer)
+ write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+} else {
+ tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+ wb_poi_lst <- read_sf(temp_gpkg, WBs_layer)
+}
+```
+
+```{r waterbody outlet POIs}
+# Derive or load Waterbody POIs ----------------------
+if("Type_WBOut" %in% names(tmp_POIs)) {
+
+ # Hillari IDs
+ hilarri_pois <- filter(tmp_POIs, !is.na(Type_hilarri))
+
+ # Get NHDARea polygons already define
+ nhd_area_vpu <- read_sf(data_paths$nhdplus_gdb, "NHDArea") %>%
+ select(COMID, AREASQKM) %>%
+ mutate(source = "NHDv2Area") %>%
+ st_compatibalize(wb_poi_lst)
+
+ # Non-nhdv2 waterbodies from resopsus
+ other_wbs <- filter(resops_wb_df, !wb_id %in% wb_poi_lst$wb_id,
+ source != "REACH") %>%
+ select(grand_id, member_comid = wbareacomi, resops_flowlcomid = flowlcomid,
+ source, accounted) %>%
+ left_join(select(st_drop_geometry(hilarri_pois),
+ COMID, hilarriid = Type_hilarri),
+ by = c("resops_flowlcomid" = "COMID")) %>%
+ group_by(resops_flowlcomid) %>%
+ filter(row_number() == 1) %>%
+ ungroup() %>%
+ left_join(select(nhd_area_vpu, COMID) %>%
+ mutate(COMID = as.character(COMID)),
+ by = c("member_comid" = "COMID")) %>%
+ st_as_sf() %>%
+ mutate(wb_id = ifelse(source == "NHDAREA", member_comid, NA)) %>%
+ st_drop_geometry()
+
+ # copy empty geometries to wb_list for use in wb POI function
+ if(nrow(other_wbs) > 0){
+ wb_poi_list_fn <- data.table::rbindlist(list(wb_poi_lst,
+ other_wbs),
+ fill = TRUE) %>%
+ distinct() %>%
+ st_as_sf()
+ } else {
+ wb_poi_list_fn <- wb_poi_lst
+ }
+
+ # Final waterbody list with all requirements
+ final_wb_list <- filter(wb_poi_list_fn,
+ area_sqkm > wb_area_thresh |
+ !is.na(resops_flowlcomid) |
+ !is.na(hilarriid))
+
+ new_path <- "?prefix=StagedProducts/Hydrography/NHDPlusHR/VPU/Current/GDB/"
+ assign("nhdhr_file_list", new_path, envir = nhdplusTools:::nhdplusTools_env)
+
+ # fix something here
+ wb_layers <- wbout_POI_creaton(filter(nhd, poi == 1), final_wb_list, data_paths, tmp_POIs, proj_crs)
+
+ if(!is.data.frame(wb_layers)){
+ # write_sf(wb_layers$nhd_wb_table, temp_gpkg, "wb_flowpaths")
+ #
+ wbs <- wb_layers$WBs
+ wb_pois <- wb_layers$POIs %>% filter(!is.na(Type_WBOut))
+ missing_wb <- filter(wbs, !wb_id %in% wb_pois$Type_WBOut)
+ print(nrow(missing_wb))
+
+ double_wb <- wbs %>% group_by(wb_id) %>% filter(n() > 1)
+ write_sf(double_wb, temp_gpkg, "double_wb")
- # Write out waterbodies
- write_sf(ref_WB, temp_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"))
+ double_wb_pois <- wb_pois %>%
+ group_by(Type_WBOut) %>%
+ filter(n() > 1)
+ write_sf(double_wb_pois, temp_gpkg, "double_wb_pois")
- write.csv(resops_POIs_df, file.path("data/reservoir_Data",
- paste0("resops_POIs_", vpu_codes, ".csv")))
- }
+ wb_missing_poi <- filter(wbs,
+ !wb_id %in% wb_pois$Type_WBOut)
+ write_sf(wb_missing_poi, temp_gpkg, "wb_missing_poi")
+ poi_missing_wb <- filter(wb_pois, !Type_WBOut %in% wb_layers$WBs$wb_id)
+ write_sf(poi_missing_wb, temp_gpkg, "poi_missing_wb")
+
+ #*************************************
+ if(all(!is.na(wb_layers$events))){
+ wb_events <- wb_layers$events %>%
+ select(COMID, REACHCODE, REACH_meas, POI_identifier = wb_id,
+ event_type, nexus) %>%
+ st_compatibalize(., events)
+
+ events <- rbind(events, wb_events)
+ write_sf(events, temp_gpkg, split_layer)
+ }
+
+ if(all(!is.na(wb_layers$events))){
+ # Write out events and outlets
+ if(!needs_layer(temp_gpkg, split_layer)){
+ events <- read_sf(temp_gpkg, split_layer) %>%
+ rbind(st_compatibalize(wb_events,.)) %>%
+ unique()
+
+ write_sf(events, temp_gpkg, split_layer)
+ } else {
+ write_sf(wb_events, nav_gpkg, split_layer)
+ }
+ }
+
+ tmp_POIs <- wb_layers$POIs
+ # Write out waterbodies
+ WB_VPU_pois <- wb_layers$WBs
+ } else {
+ tmp_POIs <- wb_layers
+ WB_VPU_pois <- final_wb_list
+ }
+
+ #12047928
+ wb_out_col <- wb_poi_collapse(tmp_POIs, nhd, events)
+ tmp_POIs <- wb_out_col$POIs
+
+ if(all(!is.na(wb_out_col$events_ret))){
+ write_sf(wb_out_col$events_ret, temp_gpkg, split_layer)
+ }
+
+ nhd <- wb_layers$nhd_wb_table
write_sf(nhd, nav_gpkg, nhd_flowline)
+ write_sf(WB_VPU_pois, temp_gpkg, WBs_layer)
write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-} else {
+ } else {
tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
nhd <- read_sf(nav_gpkg, nhd_flowline)
- ref_WB <- read_sf(temp_gpkg, WBs_layer)
- resops_POIs_df <- read.csv(file.path("data/reservoir_Data", paste0("resops_POIs_", VPU,".csv")))
+ WB_VPU_pois <- read_sf(temp_gpkg, WBs_layer)
}
-mapview(filter(tmp_POIs, !is.na(Type_WBOut)), layer.name = "Waterbody outlet POIs", col.regions = "red")
+#mapview(filter(tmp_POIs, !is.na(Type_WBOut)), layer.name = "Waterbody outlet POIs", col.regions = "red")
+```
+
+```{r AR POIs}
+# Derive or load Waterbody POIs ----------------------
+if(!"Type_AR" %in% names(tmp_POIs)) {
+ # Sparrow AR Events
+ ar_events <- readxl::read_xlsx(file.path(data_paths$USGS_IT_path, "ARevents_final7.xlsx")) %>%
+ filter(fromComid > 0 | ToComid > 0) %>%
+ select(COMID = NHDComid, AREventID) %>%
+ filter(COMID %in% nhd$COMID)
+
+ if(nrow(ar_events) > 0){
+ # AR POIs
+ tmp_POIs <- POI_creation(ar_events, nhd, "AR") %>%
+ addType(., tmp_POIs, "AR", nexus = TRUE)
+
+ write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+ } else {
+ print(paste0("No AR events in ", vpu_codes))
+ }
+
+} else {
+ tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+}
```
```{r Terminal POIs}
# Derive or load Terminal POIs ----------------------
if(!"Type_Term" %in% names(tmp_POIs)) {
-
+
# 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) %>%
+ filter(!COMID %in% term_paths$COMID, TotDASqKM >= min_da_km_terminal) %>%
bind_rows(term_paths) %>%
# Use level path identifier as Type ID
select(COMID, LevelPathI)
- tmp_POIs <- POI_creation(terminal_POIs, filter(nhd, poi == 1), "Term") %>%
- addType(., tmp_POIs, "Term")
+ tmp_POIs <- POI_creation(terminal_POIs, nhd, "Term") %>%
+ addType(., tmp_POIs, "Term", nexus = TRUE)
write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
} else {
@@ -290,7 +551,7 @@ if(!"Type_Term" %in% names(tmp_POIs)) {
nhd <- read_sf(nav_gpkg, nhd_flowline)
}
-mapview(filter(tmp_POIs, !is.na(Type_Term)), layer.name = "Terminal POIs", col.regions = "blue")
+mapview(filter(tmp_POIs, !is.na(Type_Term)), layer.name = "Terminal POIs", col.regions = "blue")
```
```{r Confluence POIs}
@@ -298,12 +559,12 @@ mapview(filter(tmp_POIs, !is.na(Type_Term)), layer.name = "Terminal POIs", col.r
if(!"Type_Conf" %in% names(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)))
+ 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))
-
+ 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
@@ -314,12 +575,12 @@ if(!"Type_Conf" %in% names(tmp_POIs)) {
group_by(DnHydroseq) %>%
filter(n()> 1) %>%
mutate(Type_Conf = LevelPathI) %>%
- ungroup() %>%
+ ungroup() %>%
select(COMID, Type_Conf)
tmp_POIs <- POI_creation(conf_COMIDs, filter(nhd, minNet == 1), "Conf") %>%
- addType(., tmp_POIs, "Conf")
-
+ addType(., tmp_POIs, "Conf", nexus = TRUE)
+
write_sf(nhd, nav_gpkg, nhd_flowline)
write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
} else {
@@ -327,21 +588,22 @@ if(!"Type_Conf" %in% names(tmp_POIs)) {
tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
}
-mapview(filter(tmp_POIs, !is.na(Type_Conf)), layer.name = "Confluence POIs", col.regions = "blue")
+mapview(filter(tmp_POIs, !is.na(Type_Conf)), layer.name = "Confluence POIs", col.regions = "blue")
```
```{r waterbody inlet POIs}
# Derive or load Waterbody POIs ----------------------
if(!"Type_WBIn" %in% names(tmp_POIs)) {
- wb_layers <- wbin_POIcreation(nhd, ref_WB, data_paths, crs, split_layer, tmp_POIs)
+ wb_layers <- wbin_POIcreation(filter(nhd, minNet == 1), WB_VPU_pois,
+ data_paths, proj_crs, split_layer, tmp_POIs)
wb_in_col <- wb_inlet_collapse(wb_layers$POIs, nhd, events)
tmp_POIs <- wb_in_col$POIs
if(!all(is.na(wb_layers$events))) {
- wb_inlet_events <- select(wb_layers$events, -c(id, offset, Hydroseq, ToMeas)) %>%
- rename(POI_identifier = WBAREACOMI)
-
+ wb_inlet_events <- select(wb_layers$events,
+ -c(LevelPathI, Hydroseq, ToMeas))
+
# Write out events and outlets
if(!needs_layer(temp_gpkg, split_layer)){
events <- read_sf(temp_gpkg, split_layer) %>%
@@ -351,18 +613,19 @@ if(!"Type_WBIn" %in% names(tmp_POIs)) {
write_sf(wb_inlet_events, nav_gpkg, split_layer)
}
}
-
- tmp_POIs$nexus <- ifelse(is.na(tmp_POIs$nexus), FALSE, tmp_POIs$nexus)
+
+ #tmp_POIs$nexus <- ifelse(is.na(tmp_POIs$nexus), FALSE, tmp_POIs$nexus)
+ #tmp_POIs <- wb_layers$POIs
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")
-```
-
+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")
+```
+
```{r NID POIs}
# Derive or load NID POIs ----------------------
@@ -375,29 +638,11 @@ if(!"Type_NID" %in% names(tmp_POIs)) {
rename(COMID = FlowLcomid) %>%
switchDiv(., nhd) %>%
group_by(COMID) %>%
- summarize(Type_NID = paste0(unique(NIDID), collapse = " "))
-
- # 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)
-
- NID_COMIDs <- filter(NID_COMIDs, !Type_NID %in% res_ops_table$NID_ID)
- }
+ summarize(Type_NID = paste0(unique(NIDID), collapse = " "))
# Derive other NID POIs
- tmp_POIs <- POI_creation(NID_COMIDs, filter(nhd, poi == 1), "NID") %>%
- addType(., tmp_POIs, "NID", nexus = TRUE, bind = FALSE)
+ tmp_POIs <- POI_creation(NID_COMIDs, nhd, "NID") %>%
+ addType(., tmp_POIs, "NID", nexus = FALSE, bind = FALSE)
# Write out geopackage layer representing POIs for given theme
write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
@@ -406,150 +651,190 @@ if(!"Type_NID" %in% names(tmp_POIs)) {
tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
}
-mapview(filter(tmp_POIs, Type_NID != ""), layer.name = "NID POIs", col.regions = "blue")
+mapview(filter(tmp_POIs, Type_NID != ""), layer.name = "NID POIs", col.regions = "blue")
```
-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(MAXELEVSMO) - 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(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")
- }
+```{r HW DA POIs}
+if(!"Type_DA" %in% names(tmp_POIs)){
+
+ # derive incremental segments from POIs
+ inc_segs <- segment_increment(nhd,
+ filter(st_drop_geometry(nhd),
+ COMID %in% tmp_POIs$COMID)) %>%
+ # bring over VAA data
+ inner_join(select(atts_rds, COMID,
+ DnHydroseq, VA_MA, TOTMA, LENGTHKM, MAXELEVSMO,
+ MINELEVSMO, WBAREACOMI, WBAreaType, FTYPE, StartFlag,
+ AreaSqKM, TotDASqKM), by = "COMID")
+ hw_segs <- inc_segs %>%
+ group_by(POI_ID) %>%
+ filter(any(StartFlag == 1)) %>%
+ filter(any(TotDASqKM > max_da_km_hw)) %>%
+ ungroup()
+ #TODO Add magic numbers to config file
+ hw_DA_splits <- hw_segs %>%
+ st_drop_geometry() %>%
+ #filter(.data$ID %in% cat$ID) %>%
+ group_by(POI_ID) %>%
+ arrange(-Hydroseq) %>%
+ mutate(ind = att_group(TotDASqKM, min_da_km_hw)) %>%
+ ungroup() %>%
+ group_by(POI_ID, ind) %>%
+ mutate(total_length = cumsum(LENGTHKM)) %>%
+ ungroup() %>%
+ group_by(POI_ID, ind) %>%
+ mutate(set = cur_group_id()) %>%
+ filter(TotDASqKM == max(TotDASqKM) &
+ total_length > 5) %>%
+ mutate(da_class = paste(POI_ID, ind, sep = "_")) %>%
+ ungroup() %>%
+ select(-ind) %>%
+ filter(!COMID %in% tmp_POIs$COMID)
+
+ if(nrow(hw_DA_splits) > 0) {
+ tmp_POIs <- POI_creation(select(hw_DA_splits, COMID, da_class),
+ filter(nhd, poi == 1), "DA") %>%
+ addType(., tmp_POIs, "DA", nexus = TRUE)
+ }
+ write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
} else {
+ tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+}
+```
+
+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}
+if(!"Type_Elev" %in% names(tmp_POIs)){
+
+ inc_segs <- segment_increment(filter(nhd, minNet == 1),
+ filter(st_drop_geometry(nhd),
+ COMID %in% tmp_POIs$COMID)) %>%
+ # bring over VAA data
+ inner_join(select(atts_rds, COMID,
+ DnHydroseq, VA_MA, TOTMA, LENGTHKM, MAXELEVSMO,
+ MINELEVSMO, WBAREACOMI, WBAreaType, FTYPE, StartFlag,
+ AreaSqKM, TotDASqKM), by = "COMID")
- tmp_POIs$Type_Elev <- rep(NA, nrow(tmp_POIs))
+ elev_fp <- inc_segs %>%
+ group_by(POI_ID) %>%
+ arrange(Hydroseq) %>%
+ # Get elevation info
+ mutate(MAXELEVSMO = na_if(MAXELEVSMO, -9998), MINELEVSMO = na_if(MINELEVSMO, -9998),
+ elev_diff_seg = max(MAXELEVSMO) - min(MINELEVSMO),
+ total_length = cumsum(LENGTHKM)) %>%
+ filter((max(MINELEVSMO) - min(MINELEVSMO)) > elev_diff) %>%
+ mutate(inc_elev_diff = c(MINELEVSMO[1], (MINELEVSMO - lag(MINELEVSMO))[-1])) %>%
+ mutate(inc_elev_diff = ifelse(inc_elev_diff == MINELEVSMO, 0, inc_elev_diff)) %>%
+ ungroup()
+
+ if(nrow(elev_fp) > 0){
+ #TODO Add magic numbers to config file
+ hw_elev_splits <- elev_fp %>%
+ st_drop_geometry() %>%
+ #filter(.data$ID %in% cat$ID) %>%
+ group_by(POI_ID) %>%
+ arrange(-Hydroseq) %>%
+ mutate(ind = cs_group(inc_elev_diff, elev_diff/2)) %>%
+ ungroup() %>%
+ group_by(POI_ID, ind) %>%
+ filter(TotDASqKM == max(TotDASqKM) &
+ total_length > 5) %>%
+ mutate(elev_class = paste(POI_ID, ind, sep = "_")) %>%
+ ungroup() %>%
+ select(-ind) %>%
+ filter(!COMID %in% tmp_POIs$COMID)
+
+ if(nrow(hw_elev_splits) > 0) {
+ tmp_POIs <- POI_creation(select(hw_elev_splits, COMID, elev_class),
+ filter(nhd, poi == 1),
+ "elev") %>%
+ addType(., tmp_POIs, "elev", nexus = TRUE)
+ write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+ }
+ }
- 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)
+
+} else {
+ tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
}
-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),
+ inc_segs <- segment_increment(nhd,
filter(st_drop_geometry(nhd),
- COMID %in% tmp_POIs$COMID, COMID %in% filter(nhd, minNet == 1)$COMID)) %>%
+ COMID %in% tmp_POIs$COMID,
+ COMID %in% nhd$COMID)) %>%
# bring over VAA data
- inner_join(select(atts, COMID, DnHydroseq, VA_MA, TOTMA, LENGTHKM, MAXELEVSMO, MINELEVSMO, WBAREACOMI, WBAreaType, FTYPE,
- AreaSqKM, TotDASqKM), by = "COMID")
+ inner_join(select(atts_rds, COMID, DnHydroseq, VA_MA, TOTMA, LENGTHKM,
+ MAXELEVSMO, MINELEVSMO, WBAREACOMI, WBAreaType, FTYPE,
+ AreaSqKM, TotDASqKM), by = "COMID")
# TT POIs
tt_pois_split <- inc_segs %>%
+ arrange(-Hydroseq) %>%
# 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()
+ filter(sum(FL_tt_hrs) > travt_diff, max(TotDASqKM) > max_elev_TT_DA) %>%
+ mutate(cum_tt = cumsum(FL_tt_hrs),
+ hrs = sum(FL_tt_hrs))
- # 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_creation(select(tt_POIs, COMID, Type_Travel), nhd, "Travel") %>%
- addType(., tmp_POIs, "Travel")
+ #TODO Add magic numbers to config file
+ tt_splits <- tt_pois_split %>%
+ st_drop_geometry() %>%
+ group_by(POI_ID) %>%
+ arrange(-Hydroseq) %>%
+ mutate(ind = cs_group(FL_tt_hrs, travt_diff * 0.75)) %>%
+ ungroup() %>%
+ group_by(POI_ID, ind) %>%
+ filter(TotDASqKM == max(TotDASqKM)) %>%
+ mutate(tt_class = paste(POI_ID, ind, sep = "_")) %>%
+ ungroup() %>%
+ select(-ind) %>%
+ filter(!COMID %in% tmp_POIs$COMID)
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+ if(nrow(tt_splits) > 0) {
+ tmp_POIs <- POI_creation(select(tt_splits, COMID, tt_class),
+ filter(nhd, poi == 1), "tt") %>%
+ addType(., tmp_POIs, "tt", nexus = 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(temp_gpkg, pois_all_layer)) {
- unCon_POIs <- filter(st_drop_geometry(filter(nhd, minNet == 1)), COMID %in% tmp_POIs$COMID, AreaSqKM == 0)
-
+ unCon_POIs <- filter(tmp_POIs, COMID %in% filter(nhd, AreaSqKM == 0)$COMID)
+
# 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_poiFix(tmp_POIs, filter(nhd, dend == 1))
+ poi_fix <- DS_poiFix(tmp_POIs, filter(nhd, minNet == 1))
new_POIs <- st_compatibalize(poi_fix$new_POIs, tmp_POIs)
xWalk <- poi_fix$xWalk
# 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))
+ 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))
-
+ final_POIs <- bind_rows(tmp_POIs_fixed, new_POIs) %>%
+ mutate(Type_Term = ifelse(nexus == 1, NA, Type_Term)) %>%
+ select(-ID)
+
# Write out fixes
write_sf(new_POIs, temp_gpkg, poi_moved_layer)
write_sf(xWalk, temp_gpkg, poi_xwalk_layer)
@@ -558,12 +843,12 @@ if(needs_layer(temp_gpkg, pois_all_layer)) {
} else {
# If no fixes designate as NA
poi_fix <- 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)
@@ -574,16 +859,20 @@ if(needs_layer(temp_gpkg, pois_all_layer)) {
final_POIs <- read_sf(temp_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)),
+ 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)) {
-
+
+ if("POI_ID" %in% colnames(nhd)){
+ nhd <- select(nhd, -POI_ID)
+ }
# 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)
+ 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)
@@ -594,15 +883,9 @@ if(needs_layer(temp_gpkg, nsegments_layer)) {
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")
+ left_join(distinct(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 {
@@ -621,15 +904,19 @@ write_sf(noDA_pois, temp_gpkg, "noDA_pois")
```{r POI Collapse}
# number POIs
-final_POIs <- mutate(final_POIs, id = row_number(), moved = NA) %>%
- write_sf(temp_gpkg, pois_all_layer)
+final_POIs_prec <- mutate(final_POIs, id = row_number(), moved = NA) %>%
+ write_sf(temp_gpkg, pois_all_layer) %>%
+ write_sf(nav_gpkg, nav_poi_layer)
collapse <- TRUE
# Load data
if(collapse){
# Move HUC12 to other POIs
- moved_pois <- poi_move(final_POIs, "Type_HUC12", .05, 2.5, c("Type_Gages", "Type_TE", "Type_NID"))
+ moved_pois <- poi_move(final_POIs_prec, "Type_HUC12", poi_dar_move,
+ poi_distance_move,
+ c("Type_Gages", "Type_TE", "Type_NID"))
+
if(!is.data.frame(moved_pois)){
final_POIs <- moved_pois$final_pois
moved_pois_table <- moved_pois$moved_points %>%
@@ -639,7 +926,8 @@ if(collapse){
}
# Gages to confluences
- moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_Conf") # 47
+ moved_pois <- poi_move(final_POIs, "Type_Gages", poi_dar_move,
+ poi_distance_move, "Type_Conf")
if(!is.data.frame(moved_pois)){
final_POIs <- moved_pois$final_pois
moved_pois_table <- moved_pois_table %>%
@@ -651,7 +939,8 @@ if(collapse){
# Gages to waterbody inlets
- moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_WBIn") # 2
+ moved_pois <- poi_move(final_POIs, "Type_Gages", poi_dar_move,
+ poi_distance_move, "Type_WBIn")
if(!is.data.frame(moved_pois)){
final_POIs <- moved_pois$final_pois
moved_pois_table <- moved_pois_table %>%
@@ -662,7 +951,8 @@ if(collapse){
}
# Gages to waterbody outlets
- moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_WBOut") # 6
+ moved_pois <- poi_move(final_POIs, "Type_Gages", poi_dar_move,
+ poi_distance_move, "Type_WBOut") # 6
if(!is.data.frame(moved_pois)){
final_POIs <- moved_pois$final_pois
moved_pois_table <- moved_pois_table %>%
@@ -672,9 +962,9 @@ if(collapse){
final_POIs <- moved_pois
}
-
# Streamgage to term outlet
- moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_Term")
+ moved_pois <- poi_move(final_POIs, "Type_Gages", poi_dar_move,
+ poi_distance_move, "Type_Term")
if(!is.data.frame(moved_pois)){
final_POIs <- moved_pois$final_pois
moved_pois_table <- moved_pois_table %>%
@@ -683,9 +973,10 @@ if(collapse){
} else {
final_POIs <- moved_pois
}
-
+
# NID to waterbody outlet
- moved_pois <- poi_move(final_POIs, "Type_NID", .025, 1, "Type_WBOut")
+ moved_pois <- poi_move(final_POIs, "Type_hilarri", poi_dar_move/10,
+ poi_distance_move * 0.4, c("Type_WBOut", "Type_TE"))
if(!is.data.frame(moved_pois)){
final_POIs <- moved_pois$final_pois
moved_pois_table <- moved_pois_table %>%
@@ -719,7 +1010,7 @@ mapview(final_POIs, layer.name = "All POIs", col.regions = "red") +
```{r lookup}
# Final POI layer
POIs <- final_POIs %>%
- arrange(COMID) %>%
+
mutate(identifier = row_number())
# Unique POI geometry
@@ -736,7 +1027,7 @@ final_POIs_table <- POIs %>%
# POI data theme table
pois_data_orig <- reshape2::melt(st_drop_geometry(select(final_POIs_table,
- -nexus)),
+ -c(nexus, nonrefactor, id, moved))),
id.vars = c("COMID", "geom_id")) %>%
filter(!is.na(value)) %>%
group_by(COMID, geom_id) %>%
@@ -744,12 +1035,16 @@ pois_data_orig <- reshape2::melt(st_drop_geometry(select(final_POIs_table,
rename(hy_id = COMID, poi_id = identifier, hl_reference = variable, hl_link = value) %>%
distinct()
-pois_data_moved <- select(st_drop_geometry(moved_pois_table), hy_id = COMID, hl_link = new_val, hl_reference = moved_value) %>%
- inner_join(distinct(pois_data_orig, hy_id, geom_id, poi_id), by = "hy_id")
-
-pois_data <- data.table::rbindlist(list(pois_data_moved, pois_data_orig), use.names = TRUE) %>%
- filter(!hl_reference %in% c("id", "moved"))
-
+if(exists("moved_pois_table")){
+ pois_data_moved <- select(st_drop_geometry(moved_pois_table),
+ hy_id = COMID, hl_link = new_val, hl_reference = moved_value) %>%
+ inner_join(distinct(pois_data_orig, hy_id, geom_id, poi_id), by = "hy_id")
+
+ pois_data <- data.table::rbindlist(list(pois_data_moved, pois_data_orig), use.names = TRUE) %>%
+ filter(!hl_reference %in% c("id", "moved"))
+} else {
+ pois_data <- pois_data_orig
+}
# POI Geometry table
poi_geometry <- select(final_POIs_table, hy_id = COMID, geom_id) %>%
diff --git a/workspace/R/NHD_navigate.R b/workspace/R/NHD_navigate.R
index c40dafe2b53936a46694dfe3a11057d5a0043ce0..06fa3c76d7d06525a0dfa0a8273288a8b84420f5 100644
--- a/workspace/R/NHD_navigate.R
+++ b/workspace/R/NHD_navigate.R
@@ -89,7 +89,7 @@ segment_increment <- function(nhdDF, POIs){
#' @param routing_fix (sf data.frame) any additional routing fixes
#'
#' @return (sf data.frame) data.frame of segments
-segment_creation <- function(nhdDF, routing_fix){
+segment_creation <- function(nhdDF, routing_fix = NULL){
if(!"StartFlag" %in% names(nhdDF)) {
nhdDF$StartFlag <- ifelse(nhdDF$Hydroseq %in% nhdDF$DnHydroseq, 0, 1)
@@ -98,7 +98,7 @@ segment_creation <- function(nhdDF, routing_fix){
in_segs <- filter(nhdDF, !is.na(POI_ID))
# If there are routing fixes to account for if a POI with a DA of 0 is moved upsream or downstream
- if (!missing(routing_fix)){
+ if (is.data.frame(routing_fix)){
routing_fix <- routing_fix %>%
rename(COMID = oldPOI, new_COMID = COMID)
@@ -145,57 +145,94 @@ segment_creation <- function(nhdDF, routing_fix){
#'
#' @return (sf data.frame) data.frame of POIs with new COMID associated
DS_poiFix <- function(POIs_wgeom, nhd){
-
+ nhd <- distinct(nhd)
POIs <- st_drop_geometry(POIs_wgeom) %>%
arrange(COMID) %>%
filter(nexus == FALSE)
# DF of downstream segment
- tocomDF <- select(nhd, COMID, DnHydroseq) %>%
- inner_join(select(st_drop_geometry(nhd), COMID, Hydroseq), by = c("DnHydroseq" = "Hydroseq"))
-
+ tocomDF <- select(st_drop_geometry(nhd), COMID, Hydroseq, TotDASqKM,
+ DnHydroseq, WBAREACOMI) %>%
+ inner_join(select(st_drop_geometry(nhd), COMID_ds = COMID, Hydroseq,
+ WBAREACOMI_down = WBAREACOMI, totda_ds = TotDASqKM),
+ by = c("DnHydroseq" = "Hydroseq")) %>%
+ inner_join(select(st_drop_geometry(nhd), COMID_us = COMID, DnHydroseq,
+ WBAREACOMI_up = WBAREACOMI, totda_us = TotDASqKM),
+ by = c("Hydroseq" = "DnHydroseq"))
+
# Find segments with POIs where there is no corresponding catchment that are not terminal
- unCon_POIs <- filter(nhd, COMID %in% POIs$COMID, AreaSqKM == 0)# & Hydroseq != TerminalPa)
+ unCon_fl <- filter(nhd, COMID %in% POIs$COMID, AreaSqKM == 0)# & Hydroseq != TerminalPa)
+ unCon_POIs <- filter(POIs, COMID %in% unCon_fl$COMID)
+
+ # Get specific fixes for waterbody inlets and outlets
+ wbout <- filter(unCon_POIs, !is.na(Type_WBOut)) %>%
+ inner_join(tocomDF, by = "COMID") %>%
+ mutate(nonrefactor = ifelse(WBAREACOMI %in% WBAREACOMI_up, COMID_ds, 0),
+ new_POI = COMID_us)
+
+ wb_pois <- filter(unCon_POIs, !is.na(Type_WBIn)) %>%
+ inner_join(tocomDF, by = "COMID") %>%
+ mutate(nonrefactor = ifelse(WBAREACOMI < 0, COMID_ds, 0),
+ new_POI = COMID_us) %>%
+ rbind(wbout) %>%
+ select(-c(nexus, Hydroseq, TotDASqKM, DnHydroseq, WBAREACOMI, COMID_ds,
+ WBAREACOMI_down, totda_ds, COMID_us, WBAREACOMI_up, totda_us)) %>%
+ rename(COMID = new_POI, oldPOI = COMID)
+
+ # The rest can be resolved with drainage are ratio
+ unCon_POIs <- filter(unCon_POIs, !COMID %in% wb_pois$oldPOI)
poi_fix <- as.data.frame(do.call("rbind", lapply(unCon_POIs$COMID, movePOI_NA_DA, st_drop_geometry(nhd)))) %>%
inner_join(POIs, by = c("oldPOI" = "COMID")) %>%
- inner_join(select(st_drop_geometry(nhd), COMID), by = c("oldPOI" = "COMID"))
+ inner_join(select(st_drop_geometry(nhd), COMID), by = c("oldPOI" = "COMID")) %>%
+ select(-c(AreaSqKM, DnHydroseq, nexus, TotDASqKM)) %>%
+ distinct() %>%
+ bind_rows(wb_pois)
# Fold in new POIs with existing POIs so all the "Type" attribution will carry over
# using the minimum will ensure correct downstream hydrosequence gets carried over
- poi_fix_unique <- filter(POIs, COMID %in% poi_fix$COMID) %>%
- bind_rows(poi_fix) %>%
- group_by(COMID) %>%
- filter(n() > 1) %>%
- # Spurs warnings where there are NAs for a column
- # for a given grouped COMID, but output is what I expect
- summarise_all(funs(min(unique(.[!is.na(.)]), na.rm = T))) #%>%
- # Replace -Inf with NA where applicaable
- #mutate_all(~na_if(., -Inf))
-
- poi_fix_unique[mapply(is.infinite, poi_fix_unique)] <- NA
-
+ poi_orig <- filter(POIs, COMID %in% poi_fix$COMID) %>%
+ bind_rows(poi_fix) %>%
+ select(-oldPOI)
+
+ list_df <- dplyr::group_by(poi_orig, COMID) |>
+ group_split()
+
+ compact <- function(l) {
+ if(nrow(l) == 1) return(as.data.frame(l))
+ lapply(names(l), \(x) {
+ out <- unique(l[[x]][!is.na(l[[x]])])
+ if(!length(out)) NA
+ else if(length(out) == 1) out
+ else {
+ cat(paste("duplicate id for", unique(l$COMID),
+ "column", x,
+ "values", paste(out, collapse = ", "),
+ "using", out[1]), file = "POI-issues.log")
+ out[1]
+ }
+ }) |>
+ setNames(names(l)) |>
+ as.data.frame()
+ }
+
+ poi_merged <- bind_rows(lapply(list_df, compact))
+
+ # # Combine POI information together for redundant pois
+ # poi_merged <- poi_orig %>%
+ # select(-c(nexus, AreaSqKM, oldPOI, DnHydroseq, TotDASqKM)) %>%
+ # group_by(COMID) %>%
+ # summarise_each(funs(toString(na.omit(.))))
+ # is.na(poi_merged) <- poi_merged == ""
+
# Join new POI COMIDs and geometry with the old Type fields
- new_POIs <- bind_rows(filter(poi_fix, !COMID %in% poi_fix_unique$COMID), poi_fix_unique) %>%
+ fin_POIs <- poi_merged %>%
arrange(COMID) %>%
bind_cols(get_node(filter(nhd, COMID %in% .$COMID) %>% arrange(COMID), position = "end")) %>%
- mutate(to_com = NA) %>%
st_sf() %>%
st_compatibalize(., POIs_wgeom)
-
- # If the DS replacement flowlines with Inc Area > 0 equals the number of replacement comIDs
- replacement_tocoms <- st_drop_geometry(new_POIs) %>%
- inner_join(select(st_drop_geometry(nhd), COMID, Hydroseq), by = c("DnHydroseq" = "Hydroseq")) %>%
- mutate(to_com = COMID.y) %>%
- select(oldPOI, COMID.x, COMID.y, AreaSqKM)
- # Add downstream COMIDS to the newCOMID
- fin_POIs <- new_POIs %>%
- left_join(replacement_tocoms, by = "oldPOI") %>%
- mutate(to_com = ifelse(!is.na(COMID.y), COMID.y, NA)) %>%
- select(-c(TotDASqKM, AreaSqKM.x, DnHydroseq, COMID.x, COMID.y, AreaSqKM.y))
-
- return (list(xWalk = select(poi_fix, COMID, oldPOI), new_POIs = fin_POIs))
+ return (list(xWalk = poi_fix, new_POIs = fin_POIs))
}
@@ -207,7 +244,7 @@ DS_poiFix <- function(POIs_wgeom, nhd){
#'
#' @return (data frame) DF of POIs with new COMID associated
movePOI_NA_DA <- function(poi_fix, nhdDF){
- print(poi_fix)
+ #print(poi_fix)
nhdDF <- distinct(nhdDF)
# Closest POI/US/DS
@@ -648,21 +685,20 @@ split_tt <- function(in_POI_ID, split_DF, tt_diff, max_DA){
#' wb_table - table of flowlines and outlet info for each
#' feature in wb
HR_Area_coms <- function(nhd, WBs, data_paths, crs){
- # Read in resops table
- resops_wb_df <- read.csv(data_paths$resops_NID_CW)
# Pull out rows for VPU that are NHDArea
- nhd_area_resops <- resops_wb_df %>%
- filter(FlowLcomid %in% nhd$COMID, COMI_srclyr == "NHDAREA")
+ nhd_area_resops <- WBs %>%
+ filter(resops_flowlcomid %in% nhd$COMID, source == "NHDAREA")
# Pull out rows for VPU that are NHD HR
- nhd_hr_wb_resops <- resops_wb_df %>%
- filter(FlowLcomid %in% nhd$COMID, COMI_srclyr == "HR ID AVAILABLE")
+ nhd_hr_wb_resops <- WBs %>%
+ filter(member_comid != 65000300139130) %>% # R09, in Canada
+ filter(resops_flowlcomid %in% nhd$COMID, source == "HR ID AVAILABLE") %>%
+ st_drop_geometry()
# Get reachcodes for waterbody outlets, so we have an idea of which
# NHD HR 4-digit geodatabase we may need to retrieve
- RC <- filter(nhd, COMID %in% c(nhd_area_resops$FlowLcomid,
- nhd_hr_wb_resops$FlowLcomid))$REACHCODE
+ RC <- filter(nhd, COMID %in% nhd_hr_wb_resops$resops_flowlcomid)$REACHCODE
# If no NHDArea or HR waterbodies needed return NULL
if (nrow(nhd_area_resops) == 0 & nrow(nhd_hr_wb_resops) == 0){
@@ -672,9 +708,9 @@ HR_Area_coms <- function(nhd, WBs, data_paths, crs){
# If NHDArea feature needed retrieve from National GDB
if (nrow(nhd_area_resops) > 0){
nhd_area_vpu <- read_sf(data_paths$nhdplus_gdb, "NHDArea") %>%
- filter(COMID %in% nhd_area_resops$NHDAREA) %>%
+ filter(COMID %in% nhd_area_resops$member_comid) %>%
mutate(source = "NHDv2Area")
-
+
wb <- st_transform(nhd_area_vpu, crs)
}
@@ -682,57 +718,87 @@ HR_Area_coms <- function(nhd, WBs, data_paths, crs){
if (nrow(nhd_hr_wb_resops) > 0){
# HUC04 we need to download
huc04 <- substr(RC, 1, 4)
-
+
# Download NHD HR HUC04 if we dont' have it, other wise load and
# Bind NHDHR waterbodies into one layer
hr_wb <- do.call("rbind", lapply(unique(huc04), function(x){
print(x)
-
- if(!file.exists(file.path(data_paths$nhdplus_dir, vpu,
+ vpu <- unique(substr(x, 1, 2))
+
+ if(!file.exists(file.path(data_paths$nhdplus_dir, vpu,
paste0("NHDPLUS_H_", x, "_HU4_GDB.gdb")))){
download_nhdplushr(data_paths$nhdplus_dir, unique(huc04))
}
+ gdb <- list.files(file.path(data_paths$nhdplus_dir, substr(vpu, 1, 2)),
+ pattern = paste0("_", x, "_.+gdb"), full.names = TRUE)
+
# Format to similar to NHDArea/Waterbody layers
- read_sf(file.path(data_paths$nhdplus_dir, substr(vpu, 1, 2),
- paste0("NHDPLUS_H_", x, "_HU4_GDB.gdb")), "NHDWaterbody")})) %>%
- filter(NHDPlusID %in% nhd_hr_wb_resops$HR_NHDPLUSID) %>%
- rename(COMID = NHDPlusID, GNIS_NAME = GNIS_Name, RESOLUTION = Resolution, AREASQKM = AreaSqKm, ELEVATION = Elevation,
- FTYPE = FType, FCODE = FCode, FDATE = FDate, REACHCODE = ReachCode) %>%
- select(-c(Permanent_Identifier, VisibilityFilter, VPUID)) %>%
+ read_sf(gdb, "NHDWaterbody")})) #%>%
+
+ # convert to lower case and designate new shape field
+ names(hr_wb) <- tolower(names(hr_wb))
+ st_geometry(hr_wb) <- "shape"
+
+ hr_wb <- filter(hr_wb, permanent_identifier %in% WBs$member_comid) %>%
+ rename(GNIS_NAME = gnis_name, GNIS_ID = gnis_id,
+ Permanent_Identifier = permanent_identifier, AREASQKM = areasqkm,
+ FTYPE = ftype, FCODE = fcode, FDATE = fdate,
+ REACHCODE = reachcode) %>%
+ #select(-c(permanent_identifier, visibilityfilter, vpuid)) %>%
st_zm(.) %>%
st_as_sf() %>%
- mutate(source = "NHDHR") %>%
- st_transform(crs) #%>%
- #st_compatibalize(wb)
-
- # Bind or create new object
- if(exists("wb")){
- wb <- data.table::rbindlist(list(wb, hr_wb), fill = TRUE) %>%
- st_as_sf()
- } else {
- wb <- hr_wb
- }
+ mutate(source = "NHDHR",
+ wb_id = ifelse(!is.na(GNIS_ID), GNIS_ID, Permanent_Identifier)) %>%
+ st_transform(crs)
+
+ # Bind or create new object
+ if(exists("wb")){
+ hr_wb <- st_compatibalize(hr_wb, wb)
+
+ wb <- data.table::rbindlist(list(wb, hr_wb), fill = TRUE) %>%
+ st_as_sf()
+ } else {
+ wb <- hr_wb
+ }
+
}
+
+ # # get the outlt rows from the table
+ # resops_outlet <- read.csv("data/reservoir_data/ISTARF-CONUS.csv") %>%
+ # filter(WBAREACOMI %in% wb$COMID | HR_NHDPLUSID %in% wb$COMID)
- # get the outlt rows from the table
- resops_outlet <- filter(resops_wb_df, NHDAREA %in% wb$COMID | HR_NHDPLUSID %in% wb$COMID)
+ # Clear out columns not needed
+ WBs_sub <- WBs %>%
+ select(-any_of(c("COMID", "GNIS_ID")))
+
+ WBs_fin <- st_drop_geometry(WBs_sub) %>%
+ mutate(member_comid = member_comid) %>%
+ inner_join(select(wb, Permanent_Identifier, GNIS_ID, GNIS_NAME2 = GNIS_NAME),
+ by = c("member_comid" = "Permanent_Identifier")) %>%
+ mutate(wb_id = ifelse(GNIS_ID %in% c(NA, " "), member_comid, GNIS_ID),
+ GNIS_NAME = ifelse(is.na(GNIS_NAME), GNIS_NAME2, GNIS_NAME)) %>%
+ select(-c(GNIS_ID, GNIS_NAME2)) %>%
+ st_as_sf()
# Create table of all flowlines that intersect the waterbody
- nhd_wb <- st_intersects(nhd, wb)
+ nhd_wb <- st_intersects(st_transform(nhd, st_crs(wb)), wb)
comid <- nhd[lengths(nhd_wb) > 0,]$COMID
nhd_wb_all <- nhd_wb[lengths(nhd_wb) > 0] %>%
purrr::set_names(comid) %>%
stack() %>%
# Ind is the default name for the set_names
rename(comid = ind, nhd_wb = values) %>%
- mutate(wb_comid = wb[nhd_wb,]$COMID,
- outlet = ifelse(comid %in% resops_outlet$FlowLcomid, "outlet", NA),
+ mutate(wb_comid = wb[nhd_wb,]$Permanent_Identifier,
+ outlet = ifelse(comid %in% WBs_fin$resops_flowlcomid, "outlet", NA),
comid = as.integer(as.character(comid))) %>%
- left_join(select(resops_wb_df, DAM_ID, DAM_NAME, FlowLcomid), by = c("comid" = "FlowLcomid")) %>%
- left_join(select(st_drop_geometry(wb), COMID, source), by = c("wb_comid" = "COMID"))
+ #left_join(select(wb_table, DAM_ID, DAM_NAME, resops_FlowLcomid),
+ # by = c("comid" = "resops_FlowLcomid")) %>%
+ left_join(select(st_drop_geometry(wb), wb_id, Permanent_Identifier, GNIS_ID,
+ source),
+ by = c("wb_comid" = "Permanent_Identifier"))
- return(list(wb_table = nhd_wb_all, wb = wb))
+ return(list(nhd_wb_table = nhd_wb_all, wb = WBs_fin))
}
#' Creates wb inlet and outlet events for splitting in hyRefactor
@@ -745,9 +811,9 @@ HR_Area_coms <- function(nhd, WBs, data_paths, crs){
#'
WB_event <- function(WBs, nhd_wb, type){
# split into features and table
- WBs_table <- WBs$wb_table
+ WBs_table <- WBs$nhd_wb_table
WBs_layer <- WBs$wb
-
+
if (type == "outlet"){
# get the outlet comid from the ResOps Table
outlet_fl <- filter(nhd_wb, COMID %in% filter(WBs_table, outlet == "outlet")$comid)
@@ -762,111 +828,91 @@ WB_event <- function(WBs, nhd_wb, type){
summarize(do_union = T) %>%
st_cast("LINESTRING")
- WBs_area <- filter(WBs_layer, source == "NHDv2Area")
- WBs_HR <- filter(WBs_layer, source == "NHDHR")
-
- # For NHDArea waterbodies (better congruity with th flowline st)
- if (nrow(WBs_area) > 0){
-
- # Intersect unioned FL with waterbody and get intersecting point
- outlet_pnts_int <- sf::st_intersection(ds_fl, WBs_area) #%>%
- # # Cast to point
- # st_cast("POINT")
- # group_by(LevelPathI) %>%
- # # Get furthest downstream point of intersection
- # filter(row_number() == max(row_number(), na.rm = T)) %>%
- # ungroup()
-
- outlet_point <- outlet_pnts_int[st_geometry_type(outlet_pnts_int) %in% c("POINT"),]
- if(nrow(outlet_point) > 0){
- outlet_fl <- outlet_pnts_int[st_geometry_type(outlet_pnts_int) %in% c("LINESTRING", "MULTILINESTRING"),] %>%
- st_cast("LINESTRING") %>%
- st_cast("POINT")
- outlet_pnts_int <- rbind(outlet_point, outlet_fl)
- } else {
- outlet_pnts_int <- st_cast(outlet_pnts_int, "MULTILINESTRING")
- }
-
- outlet_pnts <- outlet_pnts_int %>%
- group_by(LevelPathI) %>%
- st_cast("POINT") %>%
- mutate(id = row_number()) %>%# Was comid
- filter(id == max(id))
-
- # Derive event from point to use for splitting within hyRefactor
- wb_events <- get_flowline_index(nhd_wb, outlet_pnts) %>%
- inner_join(select(st_drop_geometry(nhd_wb), COMID, WBAREACOMI), by = "COMID") %>%
- mutate(event_type = type) %>%
- cbind(select(outlet_pnts, geom)) %>%
- st_as_sf()
- }
+ WBs_HR <- filter(WBs_layer, source == "HR ID AVAILABLE")
# For NHD HR waterbody outlets its a bit more complicated
if (nrow(WBs_HR) > 0){
# Get the flowlines intersecting the HR waterbody and find one with the
# max DA
- outlet_wb_int <- nhd_wb[lengths(st_intersects(nhd_wb, WBs_HR)) > 0,] %>%
- group_by(WBAREACOMI) %>%
+ outlet_wb_int <- nhd_wb[lengths(st_intersects(st_transform(nhd_wb, st_crs(WBs_HR)),
+ WBs_HR)) > 0,] %>%
+ group_by(wb_id) %>%
filter(TotDASqKM == max(TotDASqKM)) %>%
- ungroup()
-
+ ungroup()
+
# get the ds flo with the same levepath (JIC)
- ds_fl <- filter(nhd_wb, DnHydroseq %in% outlet_wb_int$Hydroseq,
+ ds_fl <- filter(nhd_wb, DnHydroseq %in% outlet_wb_int$Hydroseq,
LevelPathI %in% outlet_wb_int$LevelPathI)
-
+
outlet_fl <- rbind(outlet_wb_int, ds_fl)
-
+
# Cast flowlines within NHDHR waterbody to point
WB_FL_pnts <- outlet_wb_int %>%
st_cast("POINT") %>%
- group_by(WBAREACOMI) %>%
+ group_by(wb_id) %>%
mutate(pnt_id = row_number()) %>%
ungroup()
-
- # Determine which points intersect waterbody
- WB_outlet_pnts <- WB_FL_pnts[lengths(st_intersects(WB_FL_pnts, WBs_HR)) > 0,] %>%
- st_drop_geometry() %>%
- group_by(WBAREACOMI) %>%
- mutate(within_wb_id = row_number()) %>%
- filter(within_wb_id >= max(within_wb_id)) %>%
- ungroup() %>%
- select(WBAREACOMI, orig_pnt_id = pnt_id, within_wb_id)
-
- # Deriv new linestring by concating points from most upstream point
- # within waterbody to downstream so we can split at FL/waterbody
- # nexus
- outlet_FL <- WB_FL_pnts %>%
- inner_join(WB_outlet_pnts, by = "WBAREACOMI") %>%
- select(WBAREACOMI, pnt_id, orig_pnt_id, within_wb_id) %>%
- filter(pnt_id >= orig_pnt_id) %>%
- group_by(WBAREACOMI) %>%
- summarize(do_union = F) %>%
- st_cast("LINESTRING")
-
- # now run the intersection with the modified flowline
- outlet_pnts <- sf::st_intersection(outlet_FL, WBs_HR) %>%
- st_cast("POINT") %>%
- group_by(WBAREACOMI) %>%
- filter(row_number() == max(row_number(), na.rm = T)) %>%
- ungroup()
+
+ outlet_pnts <- do.call("rbind", lapply(unique(WB_FL_pnts$wb_id), function(x){
+ fl <- filter(WB_FL_pnts, wb_id == x)
+ wb <- filter(WBs_HR, wb_id == x)
+
+ # Determine which points intersect waterbody
+ WB_outlet_pnts <- fl[lengths(st_intersects(fl, wb)) > 0,] %>%
+ st_drop_geometry() %>%
+ group_by(wb_id) %>%
+ mutate(within_wb_id = row_number()) %>%
+ filter(within_wb_id >= max(within_wb_id)) %>%
+ ungroup() %>%
+ select(wb_id, orig_pnt_id = pnt_id, within_wb_id)
+
+ # Deriv new linestring by concating points from most upstream point
+ # within waterbody to downstream so we can split at FL/waterbody
+ # nexus
+ outlet_FL <- fl %>%
+ inner_join(WB_outlet_pnts, by = "wb_id") %>%
+ select(wb_id, pnt_id, orig_pnt_id, within_wb_id) %>%
+ filter(pnt_id >= orig_pnt_id) %>%
+ group_by(wb_id) %>%
+ summarize(do_union = F) %>%
+ st_cast("LINESTRING") %>%
+ filter(wb_id %in% wb$wb_id)
+
+ outlet_pnt <- sf::st_intersection(outlet_FL, WBs_HR) %>%
+ st_cast("MULTIPOINT") %>%
+ st_cast("POINT") %>% # was point
+ group_by(wb_id) %>%
+ filter(row_number() == max(row_number(), na.rm = T)) %>%
+ ungroup() %>%
+ mutate(id = row_number()) %>%
+ filter(wb_id == wb_id.1)
+ }))
# Derive the events
if(exists("wb_events")){
- hr_events <- get_flowline_index(nhd_wb, outlet_pnts) %>%
- inner_join(select(st_drop_geometry(nhd_wb), COMID, WBAREACOMI), by = "COMID") %>%
+ hr_events <- get_flowline_index(st_transform(nhd_wb, st_crs(outlet_pnts)),
+ outlet_pnts) %>%
+ inner_join(select(st_drop_geometry(nhd_wb), COMID, wb_id), by = "COMID") %>%
+ filter(wb_id %in% WBs_HR$wb_id) %>%
mutate(event_type = type) %>%
cbind(select(outlet_pnts, geom)) %>%
- st_as_sf()
-
- wb_events <- rbind(wb_events, hr_events)
+ st_as_sf() %>%
+ st_transform(st_crs(wb_events)) %>%
+ st_compatibalize(., wb_events)
+
+ wb_events <- rbind(wb_events, hr_events) %>%
+ distinct()
} else {
- wb_events <- get_flowline_index(nhd_wb, outlet_pnts) %>%
- inner_join(select(st_drop_geometry(nhd_wb), COMID, WBAREACOMI), by = "COMID") %>%
+ wb_events <- do.call("rbind", lapply(unique(outlet_pnts$wb_id), function(x){
+ outlet_pnt <- filter(outlet_pnts, wb_id == x)
+ get_flowline_index(nhd_wb, outlet_pnt)})) %>%
+ left_join(distinct(st_drop_geometry(outlet_pnts), resops_flowlcomid, wb_id),
+ by = c("COMID" = "resops_flowlcomid")) %>%
mutate(event_type = type) %>%
cbind(select(outlet_pnts, geom)) %>%
st_as_sf()
}
-
+
}
# For inlet points its alot easier for both NHDARea and NHDHR
@@ -874,7 +920,8 @@ WB_event <- function(WBs, nhd_wb, type){
start_pts <- get_node(nhd_wb, position = "start") %>%
cbind(st_drop_geometry(nhd_wb))
- inlet_FL <- nhd_wb[lengths(st_intersects(start_pts, WBs_layer)) == 0,] %>%
+ inlet_FL <- nhd_wb[lengths(st_intersects(
+ st_transform(start_pts, st_crs(WBs_layer)), WBs_layer)) == 0,] %>%
rbind(filter(nhd_wb, Hydroseq %in% .$DnHydroseq,
LevelPathI %in% .$LevelPathI)) %>%
arrange(desc(Hydroseq)) %>%
@@ -887,10 +934,10 @@ WB_event <- function(WBs, nhd_wb, type){
st_cast("LINESTRING") %>%
ungroup()
- #inlet_pnts <- sf::st_intersection(inlet_ls, st_cast(WBs_layer, "MULTILINESTRING", group_or_split = FALSE)) #%>%
-
- inlet_pnts_int <- sf::st_intersection(inlet_ls, WBs_layer) #%>%
- #st_cast("LINESTRING")
+ inlet_pnts_int <- sf::st_intersection(
+ st_transform(inlet_ls, st_crs(WBs_layer)), WBs_layer) #%>%
+ #st_cast("LINESTRING")
+
int_point <- inlet_pnts_int[st_geometry_type(inlet_pnts_int) %in% c("POINT"),]
if(nrow(int_point) > 0){
inlet_fl <- inlet_pnts_int[st_geometry_type(inlet_pnts_int) %in% c("LINESTRING", "MULTILINESTRING"),] %>%
@@ -899,34 +946,21 @@ WB_event <- function(WBs, nhd_wb, type){
}
inlet_pnts <- inlet_pnts_int %>%
- group_by(LevelPathI) %>%
+ group_by(LevelPathI) %>%
st_cast("POINT") %>%
- mutate(id = row_number()) %>%# Was comid
- filter(id == min(id))
-
- # # the intersection can return linestring features.
- # if(sf::st_geometry_type(inlet_pnts, by_geometry = FALSE) %in% c("LINESTRING", "MULTILINESTRING")) {
- # inlet_pnts2 <- inlet_pnts %>%
- # #st_cast("LINESTRING") %>%
- # st_cast("POINT") %>%
- # group_by(COMID) %>%
- # mutate(id = row_number()) %>%# Was comid
- # filter(id == 1)
- # #filter(row_number() == 1) %>%
- # #filter(row_number() == min(row_number())) %>%
- # ungroup()
- # } else {
- # inlet_pnts <- inlet_pnts %>%
- # st_cast("POINT") %>%
- # st_as_sf()
- # }
-
- wb_events <- get_flowline_index(nhd_wb, inlet_pnts) %>%
- inner_join(select(st_drop_geometry(nhd_wb), COMID, WBAREACOMI, LevelPathI), by = "COMID") %>%
+ mutate(id = row_number()) %>%
+ filter(id == min(id), !is.na(wb_id)) %>%
+ ungroup()
+
+ wb_events <- get_flowline_index(st_transform(nhd_wb, st_crs(inlet_pnts)),
+ inlet_pnts) %>%
+ inner_join(select(st_drop_geometry(nhd_wb), COMID, wb_id, LevelPathI), by = "COMID") %>%
mutate(event_type = type) %>%
- inner_join(select(inlet_pnts, LevelPathI, WB_COMID = COMID), by = "LevelPathI") %>%
- mutate(WBAREACOMI = WB_COMID) %>%
- select(-c(LevelPathI, WB_COMID)) %>%
+ inner_join(select(inlet_pnts, LevelPathI), by = "LevelPathI") %>%
+ #group_by(COMID, LEVELPATHI) %>%
+ mutate(nexus = T) %>%
+ rename(POI_identifier = wb_id) %>%
+ select(-c(id, offset)) %>%
st_as_sf()
}
return(wb_events)
@@ -946,6 +980,7 @@ gage_POI_creation <- function(tmp_POIs, gages_info, nhd, combine_meters, reach_m
# Create streamgage POIs
gage_POIs <- POI_creation(select(st_drop_geometry(gages_info), COMID, site_no), nhd, "Gages") %>%
st_compatibalize(., tmp_POIs)
+
# Avoid these for refactoring
avoid <- dplyr::filter(nhd, (sqrt(AreaSqKM) / LENGTHKM) > 3 & AreaSqKM > 1)
@@ -954,8 +989,8 @@ gage_POI_creation <- function(tmp_POIs, gages_info, nhd, combine_meters, reach_m
# bring over NHD information
inner_join(select(st_drop_geometry(gage_POIs), Type_Gages), by = c("site_no" = "Type_Gages")) %>%
# Get over the NHD attributes
- inner_join(select(st_drop_geometry(nhd), AreaSqKM, REACHCODE, LENGTHKM, COMID, FromMeas, ToMeas),
- by = "COMID") %>%
+ #inner_join(select(st_drop_geometry(nhd), AreaSqKM, REACHCODE, LENGTHKM, COMID, FromMeas, ToMeas),
+ # by = "COMID") %>%
# Apply reach measure thresholds
filter(reach_meas - FromMeas > reach_meas_thresh & AreaSqKM > 2 &
ToMeas - reach_meas > reach_meas_thresh & LENGTHKM > (combine_meters / m_per_km)) %>%
@@ -991,67 +1026,175 @@ gage_POI_creation <- function(tmp_POIs, gages_info, nhd, combine_meters, reach_m
#'
#' @return (sf data.frame) dataframe of waterbody outlet POIs
#'
-wbout_POI_creaton <- function(nhd, WBs_VPU, data_paths, crs){
- # Create waterbody outlet POIs for waterbodies that are in NHDv2 waterbody set
- wbout_COMIDs <- nhd %>%
- group_by(WBAREACOMI) %>%
+wbout_POI_creaton <- function(nhd, wb_table, data_paths, tmp_POIs, crs){
+
+ wb_components_table <- st_drop_geometry(wb_table) %>%
+ dplyr::mutate(member_comid = strsplit(member_comid, ",")) %>%
+ tidyr::unnest(cols = member_comid) %>%
+ mutate(member_comid = as.integer(member_comid)) %>%
+ distinct()
+
+ # R17, some wierd stuff going on in the waterbodies
+ if("947070203" %in% wb_components_table$member_comid){
+ wb_components_table <- wb_components_table %>%
+ filter(!wb_id %in% c(1503847, 1513298))
+ }
+
+ # Bring wb_id to NHD data frame
+ if(("wb_id") %in% colnames(nhd)){
+ nhd <- select(nhd, -wb_id)
+ }
+
+ nhd <- nhd %>%
+ left_join(distinct(st_drop_geometry(wb_components_table), wb_id, member_comid) %>%
+ mutate(member_comid = as.numeric(member_comid)),
+ by = c("WBAREACOMI" = "member_comid"))
+
+
+ # Resops that are either NHDArea or MHD HR source and are missing in the
+ # nhd attributiong table
+ resops_wb_other_sources <- filter(wb_table, source ==
+ "HR ID AVAILABLE" | (source == "NHDAREA" &
+ !wb_id %in% nhd$wb_id) &
+ accounted == 0)
+
+ # Resops that are not accounted for and NHDv2 waterbody or NHDArea sourced
+ resops_main <- filter(wb_table, !is.na(resops_flowlcomid) &
+ !wb_id %in% resops_wb_other_sources$wb_id,
+ accounted == 0)
+
+ # Remaining WB that don't have outlet COMID defined and aren't in ResOPsUS
+ WB_outlet <- filter(wb_table, is.na(resops_flowlcomid),
+ !wb_id %in% tmp_POIs$Type_WBOut)
+
+ # Size-based waterbody outlet POIs ot in NHDArea, HR, or
+ # manually defined with ResOPS
+ wbout_COMIDs <- nhd %>%
+ filter(wb_id %in% WB_outlet$wb_id) %>%
+ group_by(wb_id) %>%
slice(which.min(Hydroseq)) %>%
+ ungroup() %>%
switchDiv(., nhd) %>%
- select(COMID, WBAREACOMI, LevelPathI)
-
- tmp_POIs <- POI_creation(wbout_COMIDs, nhd, "WBOut") %>%
- addType(., tmp_POIs, "WBOut")
-
- # Get NHDArea and HI-Res waterbodies for ResOpsUS locaitons
- WBs_VPU_areaHR <- HR_Area_coms(nhd, WBs_VPU, data_paths, crs)
- print(length(WBs_VPU_areaHR))
-
- # If the return is not NA
- if(inherits(WBs_VPU_areaHR, "list")){
- # Attribute flowlines that intersect NHDARea and/or NHD HR waterbodies with
- # waterbody COMIDs
- nhd <- nhd %>%
- left_join(select(WBs_VPU_areaHR$wb_table, comid, wb_comid), by = c("COMID" = "comid")) %>%
- mutate(WBAREACOMI = ifelse(!is.na(wb_comid), wb_comid, WBAREACOMI)) %>%
- select(-wb_comid)
-
- # Full table of HR and NHDArea waterbodies and their flowlines
- # Subset down to specific outlets
- resops_wb_hr <- WBs_VPU_areaHR$wb_table %>%
- filter(!wb_comid %in% WBs_VPU$COMID, !is.na(outlet)) %>%
- rename(WBAREACOMI = wb_comid) %>%
- select(DAM_ID, DAM_NAME, COMID = WBAREACOMI, FlowLcomid = comid) %>%
- rbind(resops_wb_df)
-
- # NHD flowlines within the HR/NHDArea waterbodies
- nhd_WBs <- filter(nhd, WBAREACOMI %in% WBs_VPU_areaHR$wb_table$wb_comid)
- if(nrow(nhd_WBs) > 0){
-
- # Create outlet events for splitting HR and NHDArea waterbodies
- wb_outlet_events <- WB_event(WBs_VPU_areaHR, nhd_WBs, "outlet") %>%
- st_compatibalize(., tmp_POIs) %>%
- mutate(nexus = TRUE)
+ select(COMID, wb_id) %>%
+ mutate(wbtype = "size") %>%
+ st_drop_geometry() %>%
+ distinct()
+
+ # ResOpsUS defined outlets
+ resops_main_outlets <- select(st_drop_geometry(resops_main), COMID = resops_flowlcomid,
+ wb_id) %>%
+ mutate(wbtype = "resops")
+
+ # combine together
+ wbout_coms <- distinct(rbind(wbout_COMIDs, resops_main_outlets))
- tmp_POIs <- data.table::rbindlist(list(tmp_POIs,
- select(wb_outlet_events, COMID, Type_WBOut = WBAREACOMI, nexus)), fill = TRUE) %>%
- st_as_sf()
+ # Get NHDArea and HI-Res waterbodies for ResOpsUS locaitons
+ WBs_VPU_areaHR <- HR_Area_coms(nhd, resops_wb_other_sources, data_paths, crs)
+
+ # if there are no HR/Area waterbodies generated
+ if(all(is.na(WBs_VPU_areaHR))){
+ tmp_POIs <- POI_creation(wbout_coms,
+ nhd, "wb_id") %>%
+ addType(., tmp_POIs, "wb_id", nexus = TRUE) %>%
+ mutate(Type_WBOut = ifelse(is.na(Type_WBOut), Type_wb_id, Type_WBOut)) %>%
+ select(-Type_wb_id) %>%
+ left_join(select(st_drop_geometry(resops_main), GRAND_ID, resops_FlowLcomid),
+ by = c("COMID" = "resops_FlowLcomid")) %>%
+ mutate(Type_resops = ifelse(!is.na(GRAND_ID), GRAND_ID, Type_resops)) %>%
+ select(-GRAND_ID)
+ } else {
+ # WBarea out comids
+ wbareaout_coms <- select(st_drop_geometry(nhd), COMID, Hydroseq) %>%
+ inner_join(filter(WBs_VPU_areaHR$nhd_wb_table, !is.na(outlet)),
+ by = c("COMID" = "comid")) %>%
+ group_by(wb_id) %>%
+ filter(Hydroseq == min(Hydroseq)) %>% #, row_number() == 1
+ ungroup() %>%
+ inner_join(select(WBs_VPU_areaHR$wb, -source), by = "wb_id") %>%
+ filter(source != "NHDHR") %>%
+ select(COMID, wb_id) %>%
+ mutate(wbtype = "resops")
+
+ if(nrow(wbareaout_coms) > 0){
+ wbout_coms <- rbind(wbout_coms, wbareaout_coms)
}
-
- WBs_VPU_areaHR$wb <- nhdplusTools::st_compatibalize(
- WBs_VPU_areaHR$wb, WBs_VPU)
-
- # Append NHDArea and NHDHR waterbodies to waterbody layer
- WBs_VPU <- data.table::rbindlist(list(WBs_VPU, WBs_VPU_areaHR$wb), fill = TRUE) %>%
+
+ tmp_POIs <- POI_creation(select(wbout_coms, COMID, wb_id),
+ nhd, "wb_id") %>%
+ addType(., tmp_POIs, "wb_id", nexus = TRUE) %>%
+ mutate(Type_WBOut = ifelse(is.na(Type_WBOut) & (!nexus | is.na(nexus)),
+ Type_wb_id, Type_WBOut)) %>%
+ select(-Type_wb_id) %>%
+ left_join(select(st_drop_geometry(resops_main), grand_id, resops_flowlcomid),
+ by = c("COMID" = "resops_flowlcomid")) %>%
+ mutate(Type_resops = ifelse(!is.na(grand_id) & (!nexus | is.na(nexus)),
+ grand_id, Type_resops)) %>%
+ select(-grand_id)
+ }
+
+ # If the return is not NA
+ if(inherits(WBs_VPU_areaHR, "list")){
+ # Attribute flowlines that intersect NHDARea and/or NHD HR waterbodies with
+ # waterbody COMIDs
+ # For NHDHR, two waterbodies might be back to back and share an intersecting
+ # flowline, this may cause problems below, so deal with this case
+ WBs_VPU_areaHR$nhd_wb_table <- WBs_VPU_areaHR$nhd_wb_table %>%
+ left_join(select(st_drop_geometry(resops_wb_other_sources),
+ resops_flowlcomid, member_comid, grand_id),
+ by = c("comid" = "resops_flowlcomid",
+ "wb_comid" = "member_comid")) %>%
+ group_by(comid) %>%
+ filter(case_when(n() == 2 ~ !is.na(grand_id),
+ TRUE ~ n() == 1)) %>%
+ ungroup()
+
+ nhd <- nhd %>%
+ left_join(distinct(WBs_VPU_areaHR$nhd_wb_table, comid, wb_id_hr = wb_id),
+ by = c("COMID" = "comid")) %>%
+ mutate(wb_id = ifelse(!is.na(wb_id_hr), wb_id_hr, wb_id)) %>%
+ select(-wb_id_hr)
+
+ # NHD flowlines within the HR waterbodies
+ nhd_WBs <- filter(nhd, wb_id %in%
+ filter(WBs_VPU_areaHR$nhd_wb_table, source == "NHDHR")$wb_id)
+
+ if(nrow(nhd_WBs) > 0){
+ # Create outlet events for splitting HR and NHDArea waterbodies
+ wb_outlet_events <- WB_event(WBs_VPU_areaHR, nhd_WBs, "outlet") %>%
+ st_compatibalize(., tmp_POIs) %>%
+ mutate(nexus = TRUE, Type_resops = 2)
+
+ tmp_POIs <- data.table::rbindlist(list(tmp_POIs,
+ select(wb_outlet_events, COMID, Type_WBOut = wb_id, Type_resops, nexus)), fill = TRUE) %>%
st_as_sf()
-
- st_geometry(WBs_VPU) <- st_make_valid(st_geometry(WBs_VPU))
-
- return(list(WBs = WBs_VPU, POIs = tmp_POIs, events = wb_outlet_events))
-
} else {
- return(list(WBs = WBs_VPU, POIs = tmp_POIs, events = NA))
+ wb_outlet_events <- NA
+ }
+
+ WBs_VPU_areaHR$wb <- WBs_VPU_areaHR$wb %>%
+ sfheaders::sf_remove_holes(.) %>%
+ st_compatibalize(wb_poi_lst)
+
+ WBs_VPU_areaHR$wb <- nhdplusTools::st_compatibalize(
+ WBs_VPU_areaHR$wb, wb_table)
+
+ wb_table_filtered <- wb_table %>%
+ filter(!st_is_empty(wb_table))
+
+ # Append NHDArea and NHDHR waterbodies to waterbody layer
+ WBs_VPU <- data.table::rbindlist(list(wb_table_filtered, WBs_VPU_areaHR$wb),
+ fill = TRUE) %>%
+ st_as_sf()
+
+ st_geometry(WBs_VPU) <- st_make_valid(st_geometry(WBs_VPU))
+
+ return(list(POIs = tmp_POIs, events = wb_outlet_events, WBs = WBs_VPU,
+ nhd_wb_table = nhd))
+
+ } else {
+ return(list(POIs = tmp_POIs, events = NA, WBs = wb_table,
+ nhd_wb_table = nhd))
}
-
}
#' Creates Waterbody POIs, calls a few other functions
@@ -1064,69 +1207,90 @@ wbout_POI_creaton <- function(nhd, WBs_VPU, data_paths, crs){
#'
#' @return (sf data.frame) dataframe of WB inlet POIs
#'
-wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs,
+wbin_POIcreation <- function(nhd, wb_lyr, data_paths, crs,
split_layer, tmp_POIs){
wbout_COMIDs <- filter(tmp_POIs, !is.na(Type_WBOut))
+ WBs_in_POIs <- filter(wb_lyr, wb_id %in% wbout_COMIDs$Type_WBOut)
+
+ # Unnest to get list of waterbody COMIDS (member_comid)
+ wb_table <- st_drop_geometry(wb_lyr) %>%
+ dplyr::mutate(member_comid = strsplit(member_comid, ",")) %>%
+ tidyr::unnest(cols = member_comid) %>%
+ mutate(member_comid = as.integer(member_comid)) %>%
+ distinct()
+
+ # # Attribute flowline if waterbody exists for POIs
+ # nhd_wb <- nhd %>%
+ # left_join(distinct(wb_table, wb_id, member_comid),
+ # by = c("WBAREACOMI" = "member_comid"))
# Create waterbody inlet POIs
- wbin_COMIDs <- filter(nhd, WB == 0,
- DnHydroseq %in% filter(nhd, WB == 1)$Hydroseq) %>%
- select(-WBAREACOMI) %>%
+ wbin_COMIDs <- filter(nhd, is.na(wb_id) &
+ DnHydroseq %in% filter(nhd, !is.na(wb_id))$Hydroseq) %>%
+ select(-wb_id) %>%
switchDiv(., nhd) %>%
inner_join(st_drop_geometry(filter(nhd, minNet == 1)) %>%
- select(Hydroseq, WBAREACOMI), by = c("DnHydroseq" = "Hydroseq")) %>%
- select(COMID, WBAREACOMI, minNet) %>%
+ select(wb_id, Hydroseq), by = c("DnHydroseq" = "Hydroseq")) %>%
+ select(COMID, wb_id, minNet) %>%
group_by(COMID) %>%
# Ensure the inlets are on the network defined by confluence POIs
filter(minNet == 1)
# Headwater Waterbodies that may need the network extended through the inlet
- need_wbin <- st_drop_geometry(filter(WBs_VPU, source %in% c("ref_WB", "NHDv2WB"))) %>%
- dplyr::select(COMID)%>%
- dplyr::filter(!COMID %in% wbin_COMIDs$WBAREACOMI)
-
+ need_wbin <- st_drop_geometry(wb_lyr) %>%
+ #dplyr::select(COMID)%>%
+ dplyr::filter(!wb_id %in% wbin_COMIDs$wb_id)
+
if(nrow(need_wbin)>0){
-
- nhd_inlet <- mutate(nhd, WB = ifelse(WBAREACOMI > 0 & WBAREACOMI %in% need_wbin$COMID, 1, 0))
-
+
+ nhd_inlet <- mutate(nhd, WB = ifelse(wb_id > 0 & wb_id %in% need_wbin$wb_id, 1, 0))
+
missing_wbin_COMIDs <- filter(nhd_inlet, WB == 0,
DnHydroseq %in% filter(nhd_inlet, WB == 1)$Hydroseq) %>%
- select(-WBAREACOMI) %>%
+ select(-wb_id) %>%
switchDiv(., nhd_inlet) %>%
inner_join(st_drop_geometry(filter(nhd_inlet, minNet == 1)) %>%
- select(Hydroseq, WBAREACOMI), by = c("DnHydroseq" = "Hydroseq")) %>%
- select(COMID, WBAREACOMI) %>%
+ select(Hydroseq, wb_id),
+ by = c("DnHydroseq" = "Hydroseq")) %>%
+ select(COMID, wb_id) %>%
group_by(COMID) %>%
filter(row_number() == 1) %>%
ungroup()
-
+
missing_wbin_COMIDs <- missing_wbin_COMIDs %>%
- left_join(select(st_drop_geometry(nhd_inlet), COMID, TotDASqKM, minNet), by = c("COMID")) %>%
- group_by(COMID) %>%
+ left_join(select(st_drop_geometry(nhd_inlet), COMID, TotDASqKM)
+ , by = c("COMID")) %>%
+ group_by(COMID) %>%
slice(which.max(TotDASqKM))%>%
+ ungroup() %>%
select(-TotDASqKM)
-
+
if(nrow(missing_wbin_COMIDs) > 0){
- wbin_COMIDs <- rbind(wbin_COMIDs, missing_wbin_COMIDs)
+ wbin_COMIDs <- data.table::rbindlist(list(wbin_COMIDs,
+ missing_wbin_COMIDs),
+ fill = TRUE) %>%
+ select(COMID, wb_id)
}
}
wbin_POIs <- POI_creation(filter(st_drop_geometry(wbin_COMIDs), !COMID %in% wbout_COMIDs$COMID),
nhd, "WBIn")
-
+
# Get NHDArea and HR waterbodies
- WBs_VPU_areaHR <- HR_Area_coms(nhd, WBs_VPU, data_paths, crs)
- # Get flowlines that intersect waterbodies and are on the minimal network
- # (outlets should already be defined)
- nhd_WBs <- data.frame()
+ WBs_VPU_areaHR <- HR_Area_coms(nhd,
+ filter(wb_lyr, source == "HR ID AVAILABLE"),
+ data_paths, crs)
+
if(is.list(WBs_VPU_areaHR)) {
nhd_WBs <- filter(nhd, minNet == 1,
- COMID %in% WBs_VPU_areaHR$wb_table$comid,
+ COMID %in% WBs_VPU_areaHR$nhd_wb_table$comid,
AreaSqKM > 0)
+ } else {
+ nhd_WBs <- NA
}
- if(nrow(nhd_WBs) > 0){
+ if(is.data.frame(nhd_WBs)){
# Get the outlet events again
wb_outlet_events <- read_sf(temp_gpkg, split_layer) %>%
@@ -1159,7 +1323,8 @@ wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs,
mutate(dsCOMID = COMID, COMID = usCOMID)
if(nrow(wb_inlet_POIs) > 0) {
- wbin_POIs <- bind_rows(wbin_POIs, POI_creation(select(st_drop_geometry(wb_inlet_POIs), dsCOMID, Type_WBIn = WBAREACOMI),
+ wbin_POIs <- bind_rows(wbin_POIs, POI_creation(select(st_drop_geometry(wb_inlet_POIs), dsCOMID,
+ Type_WBIn = POI_identifier),
nhd, "WBIn"))
wb_inlet_events <- filter(wb_inlet_events, !COMID %in% wb_inlet_POIs$dsCOMID)
@@ -1168,11 +1333,11 @@ wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs,
if(nrow(wb_inlet_events) > 0){
- wbin_POIs <- addType(wbin_POIs, tmp_POIs, "WBIn")
+ wbin_POIs <- addType(wbin_POIs, tmp_POIs, "WBIn", nexus = TRUE)
wb_inlet_events <- st_compatibalize(wb_inlet_events, wbin_POIs)
wbin_POIs_fin <- data.table::rbindlist(list(wbin_POIs,
select(wb_inlet_events, COMID,
- Type_WBIn = WBAREACOMI, nexus)),
+ Type_WBIn = POI_identifier, nexus)),
fill = TRUE) %>%
st_as_sf()
}
@@ -1203,7 +1368,8 @@ wb_inlet_collapse <- function(tmp_POIs, nhd, events){
print (x)
wb_in_fl <- filter(wb_ds_ds, COMID == x)
gage_ds <- filter(wb_ds_ds, Hydroseq %in% wb_in_fl$Hydroseq |
- Hydroseq %in% wb_in_fl$DnHydroseq)
+ Hydroseq %in% wb_in_fl$DnHydroseq) %>%
+ filter(TotDASqKM == max(TotDASqKM))
gage_reach <- gage_ds %>%
group_by(REACHCODE) %>%
@@ -1271,7 +1437,8 @@ wb_inlet_collapse <- function(tmp_POIs, nhd, events){
wb_in_nexus <- filter(wb_in, nexus == TRUE)
wb_in_nhd <- filter(nhd, COMID %in% wb_in$COMID)
wb_ds_ds <- wb_in_nhd %>%
- rbind(filter(nhd, LevelPathI %in% .$LevelPathI & DnHydroseq %in% .$Hydroseq))
+ rbind(filter(nhd, LevelPathI %in% .$LevelPathI & DnHydroseq %in% .$Hydroseq)) %>%
+ distinct()
# get gages
gage_wb <- filter(tmp_POIs, !is.na(Type_Gages) & is.na(Type_WBIn)) #%>%
@@ -1283,7 +1450,8 @@ wb_inlet_collapse <- function(tmp_POIs, nhd, events){
gage_add <- filter(streamgages_VPU, site_no %in% gage_wb$Type_Gages) %>%
select(COMID, reachcode, reach_meas, site_no) %>%
inner_join(select(st_drop_geometry(gage_wb), site_no = Type_Gages, nexus),
- by = "site_no")
+ by = "site_no") %>%
+ filter(!nexus)
output <- lapply(wb_in$COMID,
gage_dist_node, wb_ds_ds, gage_add, events)
@@ -1344,7 +1512,7 @@ wb_inlet_collapse <- function(tmp_POIs, nhd, events){
wb_poi_collapse <- function(tmp_POIs, nhd, events){
gage_dist_node <- function(x, wb_ds_ds, gage_add, events){
print (x)
- wb_out_fl <- filter(wb_ds_ds, COMID == x)
+ wb_out_fl <- distinct(filter(wb_ds_ds, COMID == x))
gage_ds <- filter(wb_ds_ds, Hydroseq %in% wb_out_fl$Hydroseq |
Hydroseq %in% wb_out_fl$DnHydroseq)
@@ -1354,7 +1522,8 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
total_length = sum(LENGTHKM))
#print(nrow(gage_reach))
- gage_event <- filter(gage_add, COMID %in% gage_ds$COMID)
+ gage_event <- filter(gage_add, COMID %in% gage_ds$COMID) %>%
+ filter(reach_meas == min(reach_meas))
wb_event <- filter(events, COMID == x, event_type == "outlet") %>%
unique()
@@ -1407,7 +1576,7 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
# get waterbody outlets
wb_out <- filter(tmp_POIs, !is.na(Type_WBOut), is.na(Type_Gages))
- wb_out_nexus <- filter(wb_out, nexus == TRUE)
+ wb_out_nexus <- filter(wb_out, nexus)
wb_out_nhd <- filter(nhd, COMID %in% wb_out$COMID)
wb_ds_ds <- wb_out_nhd %>%
rbind(filter(nhd, LevelPathI %in% .$LevelPathI & Hydroseq %in% .$DnHydroseq))
@@ -1425,6 +1594,7 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
by = "site_no") %>%
distinct()
+ # 8693865
output <- lapply(wb_out$COMID,
gage_dist_node, wb_ds_ds, gage_add, events)
output_length <- output[lengths(output) > 1]
@@ -1602,4 +1772,36 @@ poi_move <- function(pois, move_category, DAR, dist, keep_category) {
return(list(final_pois = pois_final, moved_points = moved_points))
+}
+
+
+att_group <- function(a, athres) {
+ #cumsum <- 0
+ group <- 1
+ result <- numeric()
+ for (i in 1:length(a)) {
+ #cumsum <- cumsum + a[i]
+ tot_DA <- a[i]
+ if (tot_DA > athres) {
+ group <- group + 1
+ athres <- athres + athres
+ }
+ result = c(result, group)
+ }
+ return (result)
+}
+
+cs_group <- function(a, athres) {
+ cumsum <- 0
+ group <- 1
+ result <- numeric()
+ for (i in 1:length(a)) {
+ cumsum <- cumsum + a[i]
+ if (cumsum > athres) {
+ group <- group + 1
+ cumsum <- a[i]
+ }
+ result = c(result, group)
+ }
+ return (result)
}
\ No newline at end of file
diff --git a/workspace/R/user_vars.R b/workspace/R/user_vars.R
index af848b4fdd194d50876745734548981d142cb6d6..b6488a7437f2ca05cda0511b4c8aa620edd4e847 100644
--- a/workspace/R/user_vars.R
+++ b/workspace/R/user_vars.R
@@ -22,6 +22,7 @@ min_obs <- 90
# POI thresholds
# Minimum drainage area for gages (square kilometers)
min_da_km_gages <- 5
+min_da_km_terminal <- 10
# Variables used during 02_Navigate
wb_area_thresh <- 1 # Waterbody size for inclusion as POI
# minimum/max drainage area for headwaters(square kilometers)
@@ -41,7 +42,6 @@ poi_distance_move <- 2.5
# Settings for refactor workflow
split_meters <- 10000
combine_meters <- 1000
-min_da_km <- 40
reach_meas_thresh <- 10
aggregate_da_thresh_sqkm <- 3
diff --git a/workspace/R/utils.R b/workspace/R/utils.R
index d80d18cb2122a0050e45f72de9e1c3dbc79471b5..621b42af79db88fbdb54b5cea2e54a6b3432dbaa 100644
--- a/workspace/R/utils.R
+++ b/workspace/R/utils.R
@@ -10,6 +10,14 @@ if(!(require(hyfabric, quietly = TRUE) && packageVersion("hyfabric") == hyfabric
repos = NULL, type = "source")
}
+
+## check or re-initiate authorizatoin
+#' Retrieves files from SB in facet file_structure
+#' @param self
+check_auth <- function(){
+ sbtools::initialize_sciencebase_session()
+}
+
## Is available in sbtools now, but will leave till on cran.
#' Retrieves files from SB in facet file_structure
#' @param sbitem (character) ScienceBase item ID
@@ -96,11 +104,17 @@ cat_rpu <- function(fcats, nhd_gdb){
fl <- read_sf(nhd_gdb, "NHDFlowline_Network") %>%
align_nhdplus_names()
- # read the CatchmentSP
- cats <- readRDS(fcats) %>%
- st_as_sf() %>%
- st_drop_geometry() %>%
- dplyr::select(FEATUREID, SOURCEFC)
+ # Reference catchments vs. island catchments
+ if(basename(fcats) == "reference_catchments.gpkg"){
+ cats <- read_sf(fcats) %>%
+ st_as_sf() %>%
+ st_drop_geometry() %>%
+ dplyr::select(FEATUREID = featureid)
+ } else {
+ cats <- read_sf(fcats, "CatchmentSP") %>%
+ st_as_sf() %>%
+ st_drop_geometry()
+ }
# read the Flowlines
flowln_df <- fl %>%
diff --git a/workspace/cache/data_paths.json b/workspace/cache/data_paths.json
index 8c08252e542e6eb539ae0c530ccb9d9b7bc6b5ab..058366f9008126e96bdf0a5593e5369c458165d0 100644
--- a/workspace/cache/data_paths.json
+++ b/workspace/cache/data_paths.json
@@ -14,6 +14,8 @@
"ref_cat": "data/reference_fabric/reference_catchments.gpkg",
"ref_fl": "data/reference_fabric/reference_flowline.gpkg",
"nwm_fl": "data/reference_fabric/nwm_network.gpkg",
+ "VAA_fst": "data/NHDPlusNationalData/nhdplusVAA.fst",
+ "VAA_rds": "data/NHDPlusNationalData/nhd_vaa.rds",
"waterbodies_path": "data/NHDPlusNationalData/nhdplus_waterbodies.rds",
"fullcats_table": "data/NHDPlusNationalData/nhdcat_full.rds",
"islandcats_table": "data/islands/nhdcat_full.rds",
@@ -253,11 +255,13 @@
"hi_source": "data/islands/hi.gpkg",
"nwm_network": "data/NWM_parameters_v2.1/RouteLink_CONUS.nc",
"nwm_parm": "data/NWM_v2.1_channel_hydrofabric_10262020/nwm_v2_1_hydrofabric.gdb",
+ "new_nhdp_atts": "cache/enhd_nhdplusatts.csv",
"GFv11_gages_lyr": "data/GFv11/GFv11_gages.rds",
"GFv11_gdb": "data/GFv11/GFv1.1.gdb",
"GFv11_tgf": "data/GFv11/TGF.gdb",
- "gagesii_lyr": "data/SWIM_gage_loc/gagesII_9322_point_shapefile.shp",
- "res_attributes": "data/reservoir_data/ResOpsUS/attributes/reservoir_attributes.csv",
+ "gagesii_lyr": "data/SWIM_gage_loc/gagesII_9322_point_shapefile",
+ "hilari_sites": "data/HILARRI/HILARRI_v2.csv",
"istarf": "data/reservoir_data/ISTARF-CONUS.csv",
+ "istarf.1": "data/reservoir_data/ISTARF-CONUS.csv",
"resops_NID_CW": "data/reservoir_data/cw_ResOpsUS_NID.csv"
}