diff --git a/hyfabric/DESCRIPTION b/hyfabric/DESCRIPTION
index 24b87d051d84d0f753fb9b75bddc276f8c517661..df518e925e3636203778cce277db074c64dca10a 100644
--- a/hyfabric/DESCRIPTION
+++ b/hyfabric/DESCRIPTION
@@ -1,7 +1,7 @@
Package: hyfabric
Type: Package
Title: Utility functions for creating the reference geospatial fabric.
-Version: 0.5.7
+Version: 0.5.8
Authors@R: c(person(given = "David",
family = "Blodgett",
role = c("aut", "cre"),
diff --git a/hyfabric/NAMESPACE b/hyfabric/NAMESPACE
index b05f4ee65a2b615569b2980a62485b478b866fde..6be8415326a466d5b893a374924abe07f9ed9f45 100644
--- a/hyfabric/NAMESPACE
+++ b/hyfabric/NAMESPACE
@@ -17,3 +17,4 @@ importFrom(nhdplusTools,get_node)
importFrom(nhdplusTools,st_compatibalize)
importFrom(sf,st_drop_geometry)
importFrom(sf,st_layers)
+importFrom(sf, st_as_sf)
diff --git a/hyfabric/R/poi_creation.R b/hyfabric/R/poi_creation.R
index 2cc96a0a5cc42864c0ed0a01e021e7d08b9d6a10..f4be00021fd202a08977209dea8f45a87d936ba9 100644
--- a/hyfabric/R/poi_creation.R
+++ b/hyfabric/R/poi_creation.R
@@ -19,22 +19,15 @@ POI_creation<-function(srcData, nhdDF, IDfield){
POIs <- sub_segs %>%
get_node(., position = "end") %>%
mutate(COMID = sub_segs$COMID) %>%
- mutate(Type_HUC12 = NA, Type_WBIn = NA, Type_WBOut = NA, Type_Gages = NA, Type_TE = NA, Type_NID = NA, Type_Conf = NA,
- Type_Term = NA, Type_Elev = NA, Type_Travel = NA) %>%
inner_join(srcData %>% select(COMID, ID), by = "COMID") %>%
- mutate(!!(paste0("Type_", IDfield)) := ID)
-
- if(!(paste0("Type_", IDfield)) %in% colnames(POIs)){
- POIs <- POIs %>% select(COMID, Type_HUC12, Type_Gages, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf,
- Type_Term, Type_Elev, Type_Travel)
- } else {
- POIs <- POIs %>% select(COMID, Type_HUC12, Type_Gages, Type_TE,
- Type_NID, Type_WBIn, Type_WBOut, Type_Conf, Type_Term, Type_Elev, Type_Travel, !!(paste0("Type_", IDfield)))
- }
+ mutate(!!(paste0("Type_", IDfield)) := ID) %>%
+ select(COMID, !!(paste0("Type_", IDfield))) %>%
+ st_as_sf()
return(POIs)
}
+
#' Adds the type attribute for co-located POIs of multiple themes
#' @param new_POIs (sf data.frame) new POIs to be tested against existing
#' @param POIs (sf data.frame) Existing POIs
@@ -52,14 +45,19 @@ POI_creation<-function(srcData, nhdDF, IDfield){
#'
addType <- function(new_POIs, POIs, IDfield, nexus = TRUE, bind = TRUE){
+ # Compatibalize new points to add
new_POIs <- st_compatibalize(new_POIs, POIs)
- # subset Nexus POIs from Type
- if(nexus){
- nexus_POIs <- filter(POIs, nexus == TRUE)
- POIs <- filter(POIs, nexus %in% c(FALSE, NA))
+ # If nexus not a field, add it to existing POIs
+ if(!"nexus" %in% colnames(POIs)){
+ POIs <- mutate(POIs, nexus = FALSE)
+ }
+
+ if(!"nexus" %in% colnames(new_POIs)){
+ new_POIs <- mutate(new_POIs, nexus = FALSE)
}
+ # Check if ID field exists in POIs, add if not
if(paste0("Type_", IDfield) %in% colnames(POIs)){
POIs_exist <- POIs %>%
rename(ID = !!(paste0("Type_", IDfield)))
@@ -68,31 +66,53 @@ addType <- function(new_POIs, POIs, IDfield, nexus = TRUE, bind = TRUE){
mutate(ID = NA)
}
- # Add columns missing between master POI and new set
- missing_cols <- colnames(POIs)[!colnames(POIs) %in% colnames(new_POIs)]
- for(col in missing_cols){
- new_POIs <- new_POIs %>%
- mutate(!!col := NA)
+ # subset Nexus POIs from old POIs (if they exist)
+ if(nexus){
+ # subset existing nexus/split locations to keep unique
+ nexus_POIs <- filter(POIs, nexus)
+ # Non-nexus POIs
+ POIs <- filter(POIs, !nexus)
}
- POIs_newAtt <- filter(new_POIs, COMID %in% POIs$COMID) %>%
+ # Non-nexus POIs to bind attributes to existing POI
+ POIs_newAtt <- filter(new_POIs, COMID %in% POIs$COMID,
+ !nexus) %>%
rename(ID2 = !!(paste0("Type_", IDfield)))
- POIs_fin <- POIs_exist %>%
+ # Bind attributes
+ POIs_orig <- POIs_exist %>%
left_join(st_drop_geometry(POIs_newAtt) %>%
- select(COMID, ID2), by = "COMID") %>% #, all.x = TRUE) %>%
+ select(COMID, ID2), by = "COMID") %>%
mutate(ID = ifelse(!is.na(ID2), ID2, ID)) %>%
rename(!!(paste0("Type_", IDfield)) := ID) %>%
select(-ID2)
- # Bind unless indicated not to
+ # Non-nex POI new to POI set
+ POIs_new <- filter(new_POIs, !COMID %in% POIs$COMID)
+
+ # If binding and non-nexus HL values exist to bind
if(bind){
- POIs_fin <- rbind(POIs_fin, filter(new_POIs, !COMID %in% POIs_fin$COMID))
+ # If POIs exist for new HL reference
+ if (nrow(POIs_new) > 0){
+ POIs_fin <- data.table::rbindlist(list(POIs_orig,
+ POIs_new),
+ fill = TRUE) %>%
+ st_as_sf()
+ } else {
+ POIs_fin <- POIs_orig
+ }
+ } else {
+ POIs_fin <- POIs_orig
}
- # Add nexus back in if excluded
- if(nexus){
- POIs_fin <- rbind(POIs_fin, nexus_POIs)
+ # Split events that are added
+ new_nexus <- filter(new_POIs, nexus)
+
+ # Bind and is nexus
+ if(nrow(new_nexus) > 0){
+ POIs_fin <- data.table::rbindlist(list(POIs_orig, new_nexus),
+ fill = TRUE) %>%
+ st_as_sf()
}
return(POIs_fin)
diff --git a/hyfabric/tests/testthat/test_poi_creation.R b/hyfabric/tests/testthat/test_poi_creation.R
index 2c337faf3755c23333a7383c505995339d2d3ccc..d6cbb3d644e7977412654c6f2cb1e575182a4644 100644
--- a/hyfabric/tests/testthat/test_poi_creation.R
+++ b/hyfabric/tests/testthat/test_poi_creation.R
@@ -16,18 +16,18 @@ test_that("poi_creation", {
expect_equal(poi_out2$Type_CA_SO, "54321")
- poi_out <- addType(poi_out2, poi_out, "CA_SO", nexus = FALSE)
+ poi_out_at1 <- addType(poi_out2, poi_out, "CA_SO", nexus = FALSE)
+ expect_equal(poi_out_at1$Type_CA_SO, c(NA, "54321"))
- expect_equal(poi_out$Type_CA_SO, c(NA, "54321"))
+ poi_out_at2 <- addType(poi_out2, poi_out, "CA_SO", nexus = FALSE, bind = FALSE)
+ expect_equal(poi_out_at2$Type_CA_SO, c(NA))
poi <- data.frame(COMID = in_data$COMID[1], ID = "57681")
+ poi_out3 <- POI_creation(poi, in_data, "TE")
+ poi_out_at3 <- addType(poi_out3, poi_out, "TE", nexus = FALSE)
+ expect_equal(poi_out_at3$Type_TE, c("57681"))
- poi_out3 <- POI_creation(poi, in_data, "TE") %>%
- dplyr::mutate(nexus = FALSE)
-
- poi_out <- dplyr::mutate(poi_out, nexus = TRUE)
-
- poi_out_fin <- addType(poi_out3, poi_out, "TE", nexus = TRUE)
-
- expect_equal(poi_out_fin$nexus, c(FALSE, TRUE, TRUE))
+ poi_out_at4 <- addType(mutate(poi_out3, nexus = TRUE), poi_out_at1, "TE", nexus = TRUE)
+ expect_equal(poi_out_at4$Type_TE, c(NA, NA, "57681"))
+ expect_equal(poi_out_at4$nexus, c(FALSE, FALSE, TRUE))
})
diff --git a/workspace/00_get_data.Rmd b/workspace/00_get_data.Rmd
index 328abe76df71eb751cd4101830298d1395301e30..8fa7b4b6220d8c36c85a1451a370066d562161ff 100644
--- a/workspace/00_get_data.Rmd
+++ b/workspace/00_get_data.Rmd
@@ -20,6 +20,7 @@ library(dplyr)
library(sf)
source("R/utils.R")
+source("R/config.R")
if(!dir.exists("data")) {dir.create("data")}
if(!dir.exists("bin")) {dir.create("bin")}
@@ -68,11 +69,6 @@ if(!file.exists(SWIM_points_path)) {
item <- retrieve_sb("5ebe92af82ce476925e44b8f", data_dir)
- g3gf <- "gages3GFinfo.dbf"
-
- sbtools::item_file_download("5dcd5f96e4b069579760aedb", names = g3gf,
- destinations = file.path(data_dir, g3gf),
- overwrite_file = TRUE)
}
out_list <- c(out_list, list(SWIM_points_path = SWIM_points_path))
@@ -222,7 +218,8 @@ if(!file.exists(nhdplus_dir)) {
suppressWarnings(staged_nhd <- stage_national_data())
x <- tryCatch(
- download_nhdplusv2(islands_dir, "https://s3.amazonaws.com/edap-nhdplus/NHDPlusV21/Data/NationalData/NHDPlusV21_NationalData_Seamless_Geodatabase_HI_PR_VI_PI_03.7z"),
+ download_nhdplusv2(islands_dir, url = paste0("https://edap-ow-data-commons.s3.amazonaws.com/NHDPlusV21/",
+ "Data/NationalData/NHDPlusV21_NationalData_Seamless", "_Geodatabase_HI_PR_VI_PI_03.7z")),
# Quiet the download, overwrite existing files
error = function(e)
@@ -262,6 +259,44 @@ out_list <- c(out_list, list(nhdplus_dir = nhdplus_dir, nhdplus_gdb = nhdplus_gd
nhdplus_rpu = rpu))
```
+```{r Reference Fabric}
+ref_fab_path <- file.path(data_dir, "reference_fabric")
+ref_cat <- file.path(ref_fab_path, "reference_catchments.gpkg")
+ref_fl <- file.path(ref_fab_path, "reference_flowline.gpkg")
+nwm_fl <- file.path(ref_fab_path, "nwm_network.gpkg")
+
+if(!dir.exists(ref_fab_path)) {
+ dir.create(ref_fab_path, recursive = TRUE, showWarnings = FALSE)
+}
+
+if(!file.exists(ref_cat)){
+ possible_vpu <- c("01", "08", "10L", "15", "02", "04", "05", "06", "07", "09",
+ "03S", "03W", "03N", "10U", "11", "12", "13", "14", "16",
+ "17", "18")
+ if(is.null(sbtools::current_session()))
+ authenticate_sb()
+
+ for (vpu in possible_vpu){
+ vpu_gpkg <- paste0(vpu, "_reference_features.gpkg")
+ sbtools::item_file_download("61295190d34e40dd9c06bcd7", names = vpu_gpkg,
+ destinations = file.path(ref_fab_path, vpu_gpkg))
+ }
+
+ sbtools::item_file_download("61295190d34e40dd9c06bcd7", names = "reference_catchments.gpkg",
+ destinations = ref_cat)
+
+ sbtools::item_file_download("61295190d34e40dd9c06bcd7", names = "reference_flowline.gpkg",
+ destinations = ref_fl)
+
+ sbtools::item_file_download("61295190d34e40dd9c06bcd7", names = "nwm_network.gpkg",
+ destinations = nwm_fl)
+}
+
+
+out_list <- c(out_list, list(ref_fab_path = ref_fab_path, ref_cat = ref_cat, ref_fl = ref_fl, nwm_fl = nwm_fl))
+```
+
+
```{r NHDPlusV2 Waterbodies}
# Waterbodies - derived after downloading and post-processing NHDPlus Seamless National Geodatabase
waterbodies_path <- file.path(nhdplus_dir, "nhdplus_waterbodies.rds")
@@ -273,7 +308,7 @@ if(!file.exists(waterbodies_path)) {
# Read the feature class
wbSF <- read_sf(nhdplus_gdb, "NHDWaterbody") %>%
- st_transform(crs)
+ st_transform(5070)
# Convert to simple feature and save out
saveRDS(wbSF, waterbodies_path)
@@ -326,7 +361,6 @@ out$fac <- as.list(setNames(fac, paste0("rpu_", rpu)))
out_list<- c(out_list, out)
```
-
```{r NHDPlusV2 elev}
# NHDPlus FDR/FAC grids available by raster processing unit
@@ -355,15 +389,18 @@ out_list<- c(out_list, out)
wbd_dir <- file.path(data_dir, "wbd")
+wbd_file <- "WBD_National_GDB"
if(!dir.exists(wbd_dir)) {
dir.create(wbd_dir, recursive = TRUE)
wbd <- download_wbd(wbd_dir, "https://prd-tnm.s3.amazonaws.com/StagedProducts/Hydrography/WBD/National/GDB/WBD_National_GDB.zip")
}
-out_gdb <- file.path(wbd_dir, "WBD_National_GDB.gdb")
+out_gdb <- file.path(wbd_dir, paste0(wbd_file, ".gdb"))
out <- list(latest_wbd = file.path(wbd_dir, "WBD.rds"))
if(!file.exists(out$latest_wbd)) {
+ system(paste0(sevenz, " x ", file.path(wbd_dir, paste0(wbd_file, ".zip")), " -o", wbd_dir), ignore.stdout = TRUE, intern = TRUE)
+
# Read the feature class
wbdfc <- sf::read_sf(out_gdb,
"WBDHU12") %>%
@@ -416,16 +453,17 @@ get_sbfile <- function(f, itm) {
sbtools::item_file_download(itm, names = basename(f),
destinations = f)
-
- zip::unzip(f)
}
return(o)
}
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"))
+ 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)
```
@@ -470,6 +508,7 @@ if(!file.exists(file.path(islands_dir, "hi.gpkg"))) {
out_list <- c(out_list, out_hi)
```
+
```{r nwm_topology}
nwm_targz_url <- "https://www.nohrsc.noaa.gov/pub/staff/keicher/NWM_live/NWM_parameters/NWM_parameter_files_v2.1.tar.gz"
nwm_parm_url <- "https://www.nohrsc.noaa.gov/pub/staff/keicher/NWM_live/web/data_tools/NWM_v2.1_channel_hydrofabric.tar.gz"
@@ -508,7 +547,7 @@ out_list <- c(out_list, out)
```{r nhdplus_attributes}
-out <- list(new_nhdp_atts = file.path("data", (sb_f <- "enhd_nhdplusatts.csv")))
+out <- list(new_nhdp_atts = file.path("cache", (sb_f <- "enhd_nhdplusatts.csv")))
if(!file.exists(out$new_nhdp_atts)) {
if(is.null(sbtools::current_session()))
@@ -581,60 +620,6 @@ if(!file.exists(g2_out$gagesii_lyr)) {
out_list <- c(out_list, g2_out)
```
-```{r updated_flowlines}
-
-out <- list(ref_flowline = file.path(out_list$nhdplus_dir, (sb_f <- "reference_flowline.gpkg")))
-
-if(!file.exists(out$ref_flowline)) {
- if(is.null(sbtools::current_session()))
- authenticate_sb()
-
- fs <- sbtools::item_list_files("61295190d34e40dd9c06bcd7")
-
- if(!sb_f %in% fs$fname) {
-
- stop("you have to run this part interactively")
-
- staged_nhd <- stage_national_data(nhdplus_data = out_list$nhdplus_gdb,
- output_path = out_list$nhdplus_dir)
-
- future::plan(future::multisession(workers = 16))
- # This will generate
- fix_headwaters(staged_nhd$flowline,
- out$ref_flowline,
- new_atts = out_list$new_nhdp_atts,
- nhdpath = out_list$nhdplus_gdb)
-
- sbtools::authenticate_sb()
- sbtools::item_replace_files("61295190d34e40dd9c06bcd7", out$ref_flowline)
-
- }
-
- sbtools::item_file_download("61295190d34e40dd9c06bcd7",
- names = sb_f,
- destinations = out$ref_flowline)
-}
-
-out_list <- c(out_list, out)
-```
-
-```{r reference_catchments}
-
-out <- list(ref_catchment = file.path(out_list$nhdplus_dir, (sb_c <- "reference_catchments.gpkg")))
-
-if(!file.exists(out$ref_catchment)) {
- if(is.null(sbtools::current_session()))
- authenticate_sb()
-
- sbtools::item_file_download("61295190d34e40dd9c06bcd7",
- names = sb_c,
- destinations = out$ref_catchment)
-}
-
-out_list <- c(out_list, out)
-```
-
-
```{r res}
# Download Reservoir and dams data: ResOps, ISTARF and GRanDd
#
diff --git a/workspace/01_Gage_Selection.Rmd b/workspace/01_Gage_Selection.Rmd
index e8638f21e97759124bf7c5eae90888c2757cdabb..5b9827fd9981a54ca9939fd9b71a7f989774fbfa 100644
--- a/workspace/01_Gage_Selection.Rmd
+++ b/workspace/01_Gage_Selection.Rmd
@@ -48,13 +48,20 @@ if(file.exists("cache/dropped_gages.csv")){
file.remove("cache/dropped_gages.csv")
}
-# Update location index information from other sources where necessary from NLDI reference gage set
-ref_gages <- jsonlite::fromJSON("https://github.com/internetofwater/ref_gages/releases/download/v0.5/usgs_nldi_gages.geojson",
- flatten = TRUE)$features %>%
- select(-geometry.type, -geometry.coordinates, ref_site_no = properties.id,
- ref_RC = properties.nhdpv2_REACHCODE, ref_RM = properties.nhdpv2_REACH_measure,
- ref_COM = properties.nhdpv2_COMID) %>%
- filter(!is.na(ref_COM) & ref_COM > 0)
+if(!file.exists("data/ref_gages.gpkg")){
+ # Update location index information from other sources where necessary from NLDI reference gage set
+ ref_gages <- jsonlite::fromJSON("https://github.com/internetofwater/ref_gages/releases/download/v0.5/usgs_nldi_gages.geojson",
+ flatten = TRUE)$features %>%
+ select(-geometry.type, -geometry.coordinates, ref_site_no = properties.id,
+ ref_RC = properties.nhdpv2_REACHCODE, ref_RM = properties.nhdpv2_REACH_measure,
+ ref_COM = properties.nhdpv2_COMID) %>%
+ filter(!is.na(ref_COM) & ref_COM > 0)
+
+ write_sf(ref_gages, "data/ref_gages.gpkg")
+} else {
+ ref_gages <- read_sf("data/ref_gages.gpkg")
+}
+
```
```{r GAGESIII-SWIM}
@@ -66,7 +73,7 @@ gagesIII_db_full <- read.dbf(file = file.path(data_paths$data_dir, "gages3GFinfo
gagesIII_sf_full <- read_sf(data_paths$SWIM_points_path, "SWIM_gage_loc") %>%
st_zm(drop = TRUE) %>%
select(site_no = Gage_no, comid = COMID, reachcode = REACHCODE, reach_meas = REACH_meas) %>%
- left_join(select(ref_gages, ref_site_no, ref_COM, ref_RC, ref_RM), by = c("site_no" = "ref_site_no")) %>%
+ left_join(select(st_drop_geometry(ref_gages), ref_site_no, ref_COM, ref_RC, ref_RM), by = c("site_no" = "ref_site_no")) %>%
# Update index information with NLDI gages
mutate(comid = ref_COM, reachcode = ref_RC, reach_meas = ref_RM) %>%
select(-c(ref_COM, ref_RC, ref_RM))
@@ -166,7 +173,7 @@ gageloc_sf_full <- read_sf(data_paths$nhdplus_gdb, "Gage") %>%
# Disregard locations already removed in GAGESIII
filter(!site_no %in% gagesIII_sf_full$site_no) %>%
# Update index information with NLDI gages
- left_join(select(ref_gages, ref_site_no, ref_COM, ref_RC, ref_RM), by = c("site_no" = "ref_site_no")) %>%
+ left_join(select(st_drop_geometry(ref_gages), ref_site_no, ref_COM, ref_RC, ref_RM), by = c("site_no" = "ref_site_no")) %>%
mutate(comid = ref_COM, reachcode = ref_RC, reach_meas = ref_RM) %>%
select(-c(ref_COM, ref_RC, ref_RM))
@@ -269,7 +276,7 @@ gfv11_sf <- gfv11_sf_full %>%
reach_meas = ifelse(is.na(reach_meas), FromMeas, reach_meas)) %>%
select(-c(NHDPlusID, REACHCODE, FromMeas)) %>%
# Update index information with NLDI gages
- left_join(select(ref_gages, ref_site_no, ref_COM, ref_RC, ref_RM), by = c("site_no" = "ref_site_no")) %>%
+ left_join(select(st_drop_geometry(ref_gages), ref_site_no, ref_COM, ref_RC, ref_RM), by = c("site_no" = "ref_site_no")) %>%
mutate(comid = ref_COM, reachcode = ref_RC, reach_meas = ref_RM) %>%
select(-c(ref_COM, ref_RC, ref_RM)) %>%
st_sf()
diff --git a/workspace/01_HI_prep.Rmd b/workspace/01_HI_prep.Rmd
index 2792b57789b9b545dd892442df6f2b16dd34d2b7..7af7a919a6c4861713fdc0bd1289008271e9b8dc 100644
--- a/workspace/01_HI_prep.Rmd
+++ b/workspace/01_HI_prep.Rmd
@@ -18,7 +18,7 @@ library(readxl)
library(rgdal)
# Numeric code for HI VPU
-vpu <- "20"
+vpu_code <- "20"
# directory of data
data_dir <- "data"
@@ -39,7 +39,7 @@ if(needs_layer(data_paths$hi_source, nhd_flowline)) {
# subset by VPU
nhd <- readRDS(file.path(data_paths$islands_dir, "NHDPlusNationalData", "nhdplus_flowline.rds")) %>%
- filter(VPUID == vpu) %>%
+ filter(VPUID == vpu_code) %>%
# transform to UTM Zone 4N
st_transform(crs)
@@ -60,7 +60,7 @@ if(needs_layer(data_paths$hi_source, raw_wbd)) {
# Use HUC12 layer in the HI_PR GDB
huc12 <- read_sf(data_paths$islands_gdb, layer = "HUC12") %>%
- filter(grepl(paste0("^", vpu, ".*"), .data$HUC_12)) %>%
+ filter(grepl(paste0("^", vpu_code, ".*"), .data$HUC_12)) %>%
st_transform(., crs = st_crs(nhd))
# write to gpkg
diff --git a/workspace/01_NHD_prep.Rmd b/workspace/01_NHD_prep.Rmd
index b48bcc6d45cb7498d1ecccd464160af410231beb..6841eaef02a0f1b7706725be4b5c6ed7ea192d8d 100644
--- a/workspace/01_NHD_prep.Rmd
+++ b/workspace/01_NHD_prep.Rmd
@@ -23,6 +23,8 @@ source("R/NHD_navigate.R")
# Load Configuration of environment
source("R/config.R")
+if(exists("rpu_code")){rm(rpu_code)}
+
if(!exists("vpu_codes")) {
vpu_codes <- unique(rpu_vpu$vpuid)
} else {
@@ -38,13 +40,10 @@ names(out_vpu) <- vpu_codes
out_rpu <- rep(list(list()), length(all_rpu_codes))
names(out_rpu) <- all_rpu_codes
-fline <- sf::read_sf(data_paths$ref_flowline)
-
-fline <- sf::st_cast(fline, "LINESTRING")
-
-catchment <- sf::read_sf(file.path(data_paths$nhdplus_dir, "reference_catchments.gpkg"))
-
-catchment <- sf::st_transform(catchment, crs)
+fline <- sf::read_sf(ref_gpkg, "flowlines")
+catchment <- sf::read_sf(ref_gpkg, "catchments")
+waterbodies <- sf::read_sf(ref_gpkg, "waterbodies")
+full_cat_table <- readRDS(data_paths$fullcats_table)
# we can remove truely degenerate COMIDs
# for 0 upstream area and no catchment area
@@ -58,111 +57,86 @@ if(length(degen_comid[degen_comid %in% keep_tocomid]) > 0) stop("this will break
fline <- fline[!fline$COMID %in% degen_comid, ]
-for(VPU in vpu_codes) {
-
- rm(rpu_code)
+# there are some duplicates with enhd
+nhd <- fline %>%
+ group_by(COMID) %>%
+ filter(row_number() == 1) %>%
+ ungroup() %>%
+ make_standalone()
+
+# 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)
+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))
+
+cat_network <- sf::st_drop_geometry(nhd)
+names(cat_network) <- tolower(names(cat_network))
+cat_network <- select(cat_network, comid, tocomid, lengthkm, areasqkm, totdasqkm,
+ hydroseq, levelpathi, terminalpa, dnlevelpat, dnhydroseq,
+ reachcode, frommeas, tomeas, pathlength, arbolatesu,
+ ftype, fcode, vpuid, rpuid, wbareacomi)
+
+write_sf(nhd, nav_gpkg, nhd_flowline)
+write_sf(cat_network, nav_gpkg, nhd_network)
+write_sf(waterbodies, nav_gpkg, nhd_waterbody)
+
+cats_rpu <- full_cat_table %>%
+ filter(RPUID %in% rpu_codes$rpuid)
+
+cats <- catchment %>%
+ rename(FEATUREID = featureid,
+ AREASQKM = areasqkm) %>%
+ mutate(VPUID = vpu_codes) %>%
+ left_join(select(cats_rpu, FEATUREID, RPUID), by = c("FEATUREID")) %>%
+ filter(FEATUREID %in%
+ unique(c(nhd$COMID,
+ full_cat_table$FEATUREID[full_cat_table$RPUID %in% rpu_codes$rpuid]))) %>%
+ mutate(hy_cats = ifelse(FEATUREID %in% nhd$COMID, 1, 0),
+ full_cats = ifelse(FEATUREID %in% cats_rpu$FEATUREID, 1, 0)) %>%
+ filter(full_cats == 1 | hy_cats == 1) %>%
+ st_sf()
+
+write_sf(cats, nav_gpkg, nhd_catchment)
+
+out_vpu[[vpu_codes]] <- nhd %>%
+ sf::st_drop_geometry() %>%
+ select(COMID, toCOMID) %>%
+ filter(!toCOMID %in% COMID & !toCOMID == 0)
+for(rpu_code in rpu_codes$rpuid) {
source("R/config.R")
- if (needs_layer(nav_gpkg, nhd_flowline)) {
-
- # Subset NHD by VPU
- nhd <- subset_vpu(fline, vpu, run_make_standalone = FALSE)
-
- # there are some duplicates with enhd
- nhd <- nhd %>%
- group_by(COMID) %>%
- filter(row_number() == 1) %>%
- ungroup()
-
- nhd <- make_standalone(nhd)
-
- # 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)
-
- non_dend <- unique(unlist(lapply(zero_order$COMID, NetworkNav, nhdDF = st_drop_geometry(nhd))))
-
- # Add fields to note dendritic and POI flowlines
- nhd <- nhd %>%
- mutate(dend = ifelse(!COMID %in% non_dend, 1, 0),
- poi = ifelse(!COMID %in% non_dend & TotDASqKM >= min_da_km, 1, 0))
-
- cat_network <- sf::st_drop_geometry(nhd)
-
- names(cat_network) <- tolower(names(cat_network))
-
- cat_network <- select(cat_network, comid, tocomid, lengthkm, areasqkm, totdasqkm,
- hydroseq, levelpathi, terminalpa, dnlevelpat, dnhydroseq,
- reachcode, frommeas, tomeas, pathlength, arbolatesu,
- ftype, fcode, vpuid, rpuid, wbareacomi)
+ if(needs_layer(out_refac_gpkg, nhd_flowline)) {
+ nhd_sub <- subset_rpu(nhd, rpu_code, run_make_standalone = TRUE) %>%
+ st_sf()
- write_sf(sf::st_transform(nhd, crs), nav_gpkg, nhd_flowline)
- write_sf(cat_network, nav_gpkg, nhd_network)
+ write_sf(nhd_sub, out_refac_gpkg, nhd_flowline)
cats_rpu <- full_cat_table %>%
- filter(RPUID %in% rpu_codes$rpuid)
+ filter(RPUID == rpu_code)
- cats <- catchment %>%
- rename(FEATUREID = featureid,
- AREASQKM = areasqkm,
- RPUID = rpuid,
- VPUID = vpuid) %>%
- filter(FEATUREID %in%
- unique(c(nhd$COMID,
- full_cat_table$FEATUREID[full_cat_table$RPUID %in% rpu_codes$rpuid]))) %>%
- mutate(hy_cats = ifelse(FEATUREID %in% nhd$COMID, 1, 0),
+ cats %>%
+ mutate(hy_cats = ifelse(FEATUREID %in% nhd_sub$COMID, 1, 0),
full_cats = ifelse(FEATUREID %in% cats_rpu$FEATUREID, 1, 0)) %>%
filter(full_cats == 1 | hy_cats == 1) %>%
- st_sf()
-
- write_sf(cats, nav_gpkg, nhd_catchment)
+ st_sf() %>%
+ write_sf(out_refac_gpkg, nhd_catchment)
} else {
-
- nhd <- read_sf(nav_gpkg, nhd_flowline)
- cats <- read_sf(nav_gpkg, nhd_catchment)
-
+ nhd_sub <- read_sf(out_refac_gpkg, nhd_flowline)
}
- out_vpu[[VPU]] <- nhd %>%
+ out_rpu[[rpu_code]] <- nhd_sub %>%
sf::st_drop_geometry() %>%
select(COMID, toCOMID) %>%
filter(!toCOMID %in% COMID & !toCOMID == 0)
-
- for(rpu_code in rpu_codes$rpuid) {
- source("R/config.R")
-
- if(needs_layer(out_refac_gpkg, nhd_flowline)) {
- nhd_sub <- subset_rpu(nhd, rpu_code, run_make_standalone = TRUE) %>%
- st_sf()
-
- write_sf(nhd_sub, out_refac_gpkg, nhd_flowline)
-
- cats_rpu <- full_cat_table %>%
- filter(RPUID == rpu_code)
-
- cats %>%
- mutate(hy_cats = ifelse(FEATUREID %in% nhd_sub$COMID, 1, 0),
- full_cats = ifelse(FEATUREID %in% cats_rpu$FEATUREID, 1, 0)) %>%
- filter(full_cats == 1 | hy_cats == 1) %>%
- st_sf() %>%
- write_sf(out_refac_gpkg, nhd_catchment)
-
- } else {
- nhd_sub <- read_sf(out_refac_gpkg, nhd_flowline)
- }
-
- out_rpu[[rpu_code]] <- nhd_sub %>%
- sf::st_drop_geometry() %>%
- select(COMID, toCOMID) %>%
- filter(!toCOMID %in% COMID & !toCOMID == 0)
- }
}
```
```{r}
-
if(!file.exists("cache/rpu_vpu_out.csv")) {
make_df <- function(x, d, n) {
y <- d[[x]]
diff --git a/workspace/02_AK_navigate.Rmd b/workspace/02_AK_navigate.Rmd
index e954d00e8708eff04c70cd8af640b6c43536fb48..001cddc6e38c6338f08b2a617dc194a54dde197b 100644
--- a/workspace/02_AK_navigate.Rmd
+++ b/workspace/02_AK_navigate.Rmd
@@ -25,7 +25,7 @@ library(hyRefactor)
library(tidyr)
library(hyfabric)
-vpu <- "19"
+vpu_codes <- "19"
# Load custom functions
source("R/utils.R")
diff --git a/workspace/02_NHD_navigate.Rmd b/workspace/02_NHD_navigate.Rmd
index a03d5026dc43d68bafa2292984e807822be40aee..e2d9ec9f61a8a7ff87ed5f3b41e55b74aae62fb4 100644
--- a/workspace/02_NHD_navigate.Rmd
+++ b/workspace/02_NHD_navigate.Rmd
@@ -22,8 +22,6 @@ knitr::opts_chunk$set(
source("R/utils.R")
source("R/NHD_navigate.R")
source("R/hyrefactor_funs.R")
-#source("R/wb_poi_collapse.R")
-#source("R/wb_inlet_collapse.R")
# increase timeout for data downloads
options(timeout=600)
@@ -34,17 +32,6 @@ source("R/config.R")
# Gages output from Gage_selection
gages <- read_sf(gage_info_gpkg, "Gages")
-# NWM network
-nc <- RNetCDF::open.nc(data_paths$nwm_network)
-# NWM gages
-nwm_gages <- data.frame(
- comid = RNetCDF::var.get.nc(nc, "link"),
- gage_id = RNetCDF::var.get.nc(nc, "gages")) %>%
- mutate(gage_id = gsub(" ", "", gage_id)) %>%
- mutate(gage_id = ifelse(gage_id == "", NA, gage_id))
-
-RNetCDF::close.nc(nc)
-
# need some extra attributes for a few POI analyses
atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.rds"))
```
@@ -55,14 +42,14 @@ if(needs_layer(temp_gpkg, nav_poi_layer)) {
nhd <- read_sf(nav_gpkg, nhd_flowline)
try(nhd <- select(nhd, -c(minNet, WB, struct_POI, struct_Net, POI_ID, dend, poi)), silent = TRUE)
-
+
# Some NHDPlus VPUs include HUCs from other VPUs
- if(vpu == "02"){
+ if(vpu_codes == "02"){
grep_exp <-"^02|^04"
- } else if (vpu == "08") {
+ } else if (vpu_codes == "08") {
grep_exp <- "^03|^08"
} else {
- grep_exp <- paste0("^", substr(vpu, start = 1, stop = 2))
+ grep_exp <- paste0("^", substr(vpu_codes, start = 1, stop = 2))
}
# Join HUC12 outlets with NHD
@@ -90,13 +77,12 @@ mapview(filter(tmp_POIs, Type_HUC12 != ""), layer.name = "HUC12 POIs", col.regio
```
```{r streamgage POIs}
-if(all(is.na(tmp_POIs$Type_Gages))) {
+if(!"Type_Gages" %in% names(tmp_POIs)) {
# Previously identified streamgages within Gage_Selection.Rmd
streamgages_VPU <- gages %>%
rename(COMID = comid) %>%
filter(COMID %in% nhd$COMID) %>%
- #st_drop_geometry() %>%
switchDiv(., nhd)
streamgages <- streamgages_VPU %>%
@@ -106,8 +92,9 @@ if(all(is.na(tmp_POIs$Type_Gages))) {
ungroup()
# Derive GAGE POIs; use NHD as we've already filtered by NWIS DA in the Gage selection step
- gages_POIs <- gage_POI_creation(tmp_POIs, streamgages, nhd, combine_meters, reach_meas_thresh)
+ gages_POIs <- gage_POI_creation(tmp_POIs, streamgages, filter(nhd, poi == 1), combine_meters, reach_meas_thresh)
tmp_POIs <- gages_POIs$tmp_POIs
+ events <- rename(gages_POIs$events, POI_identifier = Type_Gages)
# 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) {
@@ -117,12 +104,12 @@ if(all(is.na(tmp_POIs$Type_Gages))) {
# Start documenting gages that are dropped out; these gages have no mean daily Q
gage_document(filter(streamgages_VPU, !site_no %in% tmp_POIs$Type_Gages) %>%
select(site_no),
- "Gages_info", paste0("VPU ", vpu, "; low gage score"))
+ "Gages_info", paste0("VPU ", vpu_codes, "; low gage score"))
}
# Write out events and outelts
- write_sf(rename(gages_POIs$events, POI_identifier = Type_Gages), temp_gpkg, split_layer)
+ write_sf(events, temp_gpkg, split_layer)
# Write out geopackage layer representing POIs for given theme
write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
} else {
@@ -134,9 +121,9 @@ mapview(filter(tmp_POIs, !is.na(Type_Gages)), layer.name = "Streamgage POIs", co
```
```{r TE POIs}
-if(all(is.na(tmp_POIs$Type_TE))) {
+if(!"Type_TE" %in% names(tmp_POIs)) {
- if(vpu == "08"){
+ if(vpu_codes == "08"){
nhd$VPUID <- "08"
} else {
nhd$VPUID <- substr(nhd$RPUID, 1, 2)
@@ -146,15 +133,15 @@ if(all(is.na(tmp_POIs$Type_TE))) {
TE_COMIDs <- read_sf(data_paths$TE_points_path, "2015_TE_Model_Estimates_lat.long_COMIDs") %>%
mutate(COMID = as.integer(COMID)) %>%
inner_join(., select(st_drop_geometry(nhd), COMID, VPUID), by = "COMID") %>%
- filter(grepl(paste0("^", substr(vpu, 1, 2), ".*"), .data$VPUID), COMID > 0) %>%
+ 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()) %>%
ungroup()
# Derive TE POIs
- tmp_POIs <- POI_creation(st_drop_geometry(TE_COMIDs), filter(nhd, dend == 1), "TE") %>%
- addType(., tmp_POIs, "TE")
+ 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) {
@@ -172,38 +159,9 @@ if(all(is.na(tmp_POIs$Type_TE))) {
mapview(filter(tmp_POIs, Type_TE != ""), layer.name = "TE Plant POIs", col.regions = "blue")
```
-```{r Terminal POIs}
-# Derive or load Terminal POIs ----------------------
-if(all(is.na(tmp_POIs$Type_Term))) {
-
- # Terminal POIs on levelpaths with upstream POIs
- term_paths <- filter(st_drop_geometry(nhd), Hydroseq %in% filter(nhd, COMID %in% tmp_POIs$COMID)$TerminalPa)
-
- # Non-POI levelpath terminal pois, but meet size criteria
- terminal_POIs <- st_drop_geometry(nhd) %>%
- filter(Hydroseq == TerminalPa | (toCOMID == 0 | is.na(toCOMID) | !toCOMID %in% COMID)) %>%
- filter(!COMID %in% term_paths$COMID, TotDASqKM >= min_da_km) %>%
- bind_rows(term_paths) %>%
- # Use level path identifier as Type ID
- select(COMID, LevelPathI)
-
- tmp_POIs <- POI_creation(terminal_POIs, nhd, "Term") %>%
- addType(., tmp_POIs, "Term")
-
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-} else {
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
- nhd <- read_sf(nav_gpkg, nhd_flowline)
-}
-
-mapview(filter(tmp_POIs, !is.na(Type_Term)), layer.name = "Terminal POIs", col.regions = "blue")
-```
-
```{r waterbody outlet POIs}
# Derive or load Waterbody POIs ----------------------
-if(all(is.na(tmp_POIs$Type_WBOut))) {
-
- events <- read_sf(temp_gpkg, split_layer)
+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"),
@@ -211,6 +169,10 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
mutate(FTYPE = as.character(FTYPE)) %>%
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) %>%
@@ -239,14 +201,15 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
# Attribute flowlines that are in waterbodies
nhd <- mutate(nhd, WB = ifelse(WBAREACOMI > 0 & WBAREACOMI %in% WBs_VPU$COMID, 1, 0))
- wb_layers <- wbout_POI_creaton(nhd, WBs_VPU, data_paths, crs)
+ 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
- #tmp_POIs <- wb_layers$POIs
- wb_out_col <- wb_poi_collapse(wb_layers$POIs, nhd, events)
- ho <- filter(wb_layers$POIs, !COMID %in% wb_out_col$POIs$COMID)
- write_sf(wb_out_col$events_ret, temp_gpkg, split_layer)
+ wb_out_col <- wb_poi_collapse(wb_layers$POIs, filter(nhd, WB == 1), events)
tmp_POIs <- wb_out_col$POIs
+ if(!is.na(wb_out_col$events_ret)){
+ write_sf(wb_out_col$events_ret, temp_gpkg, split_layer)
+ }
if(!all(is.na(wb_layers$events))){
wb_events <- select(wb_layers$events, -c(id, offset)) %>%
@@ -263,8 +226,22 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
}
}
+ ref_WB <- filter(ref_WB, COMID %in% WBs_VPU$COMID)
+
+ 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)
+
+ ref_WB <- rbind(ref_WB, WBs_VPU)
+
# Write out waterbodies
- write_sf(WBs_VPU, nav_gpkg, WBs_layer)
+ write_sf(ref_WB, temp_gpkg, WBs_layer)
# Write specific ResOpsUS data to a separate sheet
if(nrow(wb_lst) > 0){
@@ -273,7 +250,8 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
dplyr::filter(!is.na(Type_WBOut)) %>%
dplyr::inner_join(select(resops_wb_df, -FlowLcomid), by = c("Type_WBOut" = "COMID"))
- write.csv(resops_POIs_df, file.path("data/reservoir_Data", paste0("resops_POIs_",vpu,".csv")))
+ write.csv(resops_POIs_df, file.path("data/reservoir_Data",
+ paste0("resops_POIs_", vpu_codes, ".csv")))
}
write_sf(nhd, nav_gpkg, nhd_flowline)
@@ -281,16 +259,43 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
} else {
tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
nhd <- read_sf(nav_gpkg, nhd_flowline)
- WBs_VPU <- read_sf(nav_gpkg, WBs_layer)
+ ref_WB <- read_sf(temp_gpkg, WBs_layer)
resops_POIs_df <- read.csv(file.path("data/reservoir_Data", paste0("resops_POIs_",vpu,".csv")))
}
mapview(filter(tmp_POIs, !is.na(Type_WBOut)), layer.name = "Waterbody outlet POIs", col.regions = "red")
```
+```{r 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) %>%
+ 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")
+
+ write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+} else {
+ tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+ nhd <- read_sf(nav_gpkg, nhd_flowline)
+}
+
+mapview(filter(tmp_POIs, !is.na(Type_Term)), layer.name = "Terminal POIs", col.regions = "blue")
+```
+
```{r Confluence POIs}
# Derive POIs at confluences where they are absent ----------------------
-if(all(is.na(tmp_POIs$Type_Conf))) {
+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)))
@@ -325,16 +330,51 @@ if(all(is.na(tmp_POIs$Type_Conf))) {
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_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)
+
+ # Write out events and outlets
+ if(!needs_layer(temp_gpkg, split_layer)){
+ events <- read_sf(temp_gpkg, split_layer) %>%
+ rbind(st_compatibalize(wb_inlet_events,.))
+ write_sf(events, temp_gpkg, split_layer)
+ } else {
+ write_sf(wb_inlet_events, nav_gpkg, split_layer)
+ }
+ }
+
+ tmp_POIs$nexus <- ifelse(is.na(tmp_POIs$nexus), FALSE, tmp_POIs$nexus)
+ write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+} else {
+ tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
+ nhd <- read_sf(nav_gpkg, nhd_flowline)
+}
+
+mapview(filter(tmp_POIs, !is.na(Type_WBIn)), layer.name = "Waterbody inlet POIs", col.regions = "blue") +
+ mapview(filter(tmp_POIs, !is.na(Type_WBOut)), layer.name = "Waterbody outlet POIs", col.regions = "red")
+```
+
+
```{r NID POIs}
# Derive or load NID POIs ----------------------
-if(all(is.na(tmp_POIs$Type_NID))) {
+if(!"Type_NID" %in% names(tmp_POIs)) {
# Read in NID shapefile
NID_COMIDs <- read_sf(data_paths$NID_points_path, "Final_NID_2018") %>%
st_drop_geometry() %>%
filter(EROM != 0, FlowLcomid %in% filter(nhd, dend ==1)$COMID) %>%
+ rename(COMID = FlowLcomid) %>%
switchDiv(., nhd) %>%
- group_by(FlowLcomid) %>%
+ group_by(COMID) %>%
summarize(Type_NID = paste0(unique(NIDID), collapse = " "))
# Determine which NID facilitites have been co-located with ResOps
@@ -349,14 +389,14 @@ if(all(is.na(tmp_POIs$Type_NID))) {
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, Type_NID)) %>%
+ 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)
}
# Derive other NID POIs
- tmp_POIs <- POI_creation(NID_COMIDs, filter(nhd, dend ==1), "NID") %>%
+ tmp_POIs <- POI_creation(NID_COMIDs, filter(nhd, poi == 1), "NID") %>%
addType(., tmp_POIs, "NID", nexus = TRUE, bind = FALSE)
# Write out geopackage layer representing POIs for given theme
@@ -369,43 +409,6 @@ if(all(is.na(tmp_POIs$Type_NID))) {
mapview(filter(tmp_POIs, Type_NID != ""), layer.name = "NID POIs", col.regions = "blue")
```
-```{r waterbody inlet POIs}
-# Derive or load Waterbody POIs ----------------------
-if(all(is.na(tmp_POIs$Type_WBIn))) {
-
- wb_layers <- wbin_POIcreation(nhd, WBs_VPU, data_paths, crs, split_layer, tmp_POIs)
- #tmp_POIs <- wb_layers$POIs
- wb_in_col <- wb_inlet_collapse(wb_layers$POIs, nhd, events)
- #ho <- filter(wb_layers$POIs, !COMID %in% wb_in_col$POIs$COMID)
- #write_sf(wb_in_col$events_ret, temp_gpkg, split_layer)
- 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)
-
- # Write out events and outlets
- if(!needs_layer(temp_gpkg, split_layer)){
- events <- read_sf(temp_gpkg, split_layer) %>%
- rbind(st_compatibalize(wb_inlet_events,.))
- write_sf(events, temp_gpkg, split_layer)
- } else {
- write_sf(wb_inlet_events, nav_gpkg, split_layer)
- }
- }
-
- tmp_POIs$nexus <- ifelse(is.na(tmp_POIs$nexus), FALSE, tmp_POIs$nexus)
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-} else {
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
- nhd <- read_sf(nav_gpkg, nhd_flowline)
-}
-
-mapview(filter(tmp_POIs, !is.na(Type_WBIn)), layer.name = "Waterbody inlet POIs", col.regions = "blue") +
- mapview(filter(tmp_POIs, !is.na(Type_WBOut)), layer.name = "Waterbody outlet POIs", col.regions = "red")
-```
-
-
This chunk breaks up potential aggregated segments where the elevation range of the are contributing to the segment exceeds 500-m
```{r Elevation Break POIs}
# derive incremental segments from POIs
@@ -424,7 +427,7 @@ elev_pois_split <- inc_segs %>%
mutate(MAXELEVSMO = na_if(MAXELEVSMO, -9998), MINELEVSMO = na_if(MINELEVSMO, -9998),
elev_diff_seg = max(MAXELEVSMO) - min(MINELEVSMO)) %>%
# Filter out to those incremental segments that need splitting above the defined elevation threshold
- filter((max(MAXELEVSMO) - min(MINELEVSMO)) > (elev_diff * 100)) %>%
+ 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
@@ -506,18 +509,18 @@ if(all(is.na(tmp_POIs$Type_Travel))) {
mutate(orig_iter = iter) %>%
ungroup()
- # For reach iteration
- tt_POIs <- do.call("rbind", lapply(unique(tt_pois_split$POI_ID), split_tt,
- tt_pois_split, travt_diff, max_elev_TT_DA)) %>%
- filter(!tt_POI_ID %in% tmp_POIs$COMID, COMID == tt_POI_ID) %>%
- mutate(Type_Travel = 1) %>%
- select(COMID, Type_Travel) %>%
- unique()
-
- tmp_POIs <- POI_creation(select(tt_POIs, COMID, Type_Travel), nhd, "Travel") %>%
- addType(., tmp_POIs, "Travel")
+ # 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()
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
+ tmp_POIs <- POI_creation(select(tt_POIs, COMID, Type_Travel), nhd, "Travel") %>%
+ addType(., tmp_POIs, "Travel")
+
+ 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)
@@ -589,6 +592,7 @@ if(needs_layer(temp_gpkg, nsegments_layer)) {
# create and write out final dissolved segments
nsegments_fin <- segment_creation(nhd_Final, xWalk)
+
nhd_Final <- select(nhd_Final, -POI_ID) %>%
left_join(select(st_drop_geometry(nsegments_fin$raw_segs), COMID, POI_ID), by = "COMID")
nsegments <- nsegments_fin$diss_segs
@@ -611,48 +615,89 @@ if(needs_layer(temp_gpkg, nsegments_layer)) {
# Ensure that all the problem POIs have been taken care of
sub <- nhd_Final %>% filter(COMID %in% final_POIs$COMID)
print (paste0(nrow(sub[sub$AreaSqKM == 0,]), " POIs on flowlines with no local drainage contributions"))
+noDA_pois <- filter(final_POIs, COMID %in% filter(sub, AreaSqKM == 0)$COMID)
+write_sf(noDA_pois, temp_gpkg, "noDA_pois")
```
```{r POI Collapse}
-# Load data
-if(needs_layer(nav_gpkg, final_poi_layer)) {
+# number POIs
+final_POIs <- mutate(final_POIs, id = row_number(), moved = NA) %>%
+ write_sf(temp_gpkg, pois_all_layer)
+collapse <- TRUE
- if(!"Type_Elev" %in% names(final_POIs)) final_POIs$Type_Elev <- rep(NA_real_, nrow(final_POIs))
- # Convert nexus to 0,1
- final_POIs <- mutate(final_POIs, nexus = ifelse(is.na(nexus), 0, nexus))
+# 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"))
+ if(!is.data.frame(moved_pois)){
+ final_POIs <- moved_pois$final_pois
+ moved_pois_table <- moved_pois$moved_points %>%
+ mutate(move_type = "huc12 to other")
+ } else {
+ final_POIs <- moved_POIs
+ }
+
+ # Gages to confluences
+ moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_Conf") # 47
+ if(!is.data.frame(moved_pois)){
+ final_POIs <- moved_pois$final_pois
+ moved_pois_table <- moved_pois_table %>%
+ rbind(moved_pois$moved_points %>%
+ mutate(move_type = "gages to conf"))
+ } else {
+ final_POIs <- moved_POIs
+ }
- #1 Move POIs downstream by category, nexus POIs get removed
- out_HUC12 <- POI_move_down(temp_gpkg, final_POIs, nsegments, filter(nhd_Final, !is.na(POI_ID)), "Type_HUC12",
- poi_dar_move *2)
- out_HUC12$allPOIs$nexus <- as.numeric(out_HUC12$allPOIs$nexus)
- out_gages <- POI_move_down(temp_gpkg, out_HUC12$allPOIs, out_HUC12$segs, out_HUC12$FL,
- "Type_Gages", poi_dar_move)
- # Convert empty strings to NA for handling within R
- out_gages$allPOIs <- mutate_all(out_gages$allPOIs, list(~na_if(.,"")))
+ # Gages to waterbody inlets
+ moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_WBIn") # 2
+ if(!is.data.frame(moved_pois)){
+ final_POIs <- moved_pois$final_pois
+ moved_pois_table <- moved_pois_table %>%
+ rbind(moved_pois$moved_points %>%
+ mutate(move_type = "gages to wbin"))
+ } else {
+ final_POIs <- moved_pois
+ }
- out_gages$allPOIs$snapped <- out_gages$allPOIs$COMID %in% new_POIs$COMID
-
- out_gages$allPOIs$identifier <- seq(1, nrow(out_gages$allPOIs))
+ # Gages to waterbody outlets
+ moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_WBOut") # 6
+ if(!is.data.frame(moved_pois)){
+ final_POIs <- moved_pois$final_pois
+ moved_pois_table <- moved_pois_table %>%
+ rbind(moved_pois$moved_points %>%
+ mutate(move_type = "gages to wbout"))
+ } else {
+ final_POIs <- moved_pois
+ }
- out_gages$allPOIs <- select(out_gages$allPOIs, identifier, everything())
- nhd_Final <- select(nhd_Final, -POI_ID) %>%
- left_join(st_drop_geometry(out_HUC12$FL) %>%
- select(COMID, POI_ID), by = "COMID")
-
- # Write out geopackage layer representing POIs for given theme
- write_sf(out_gages$allPOIs, nav_gpkg, final_poi_layer)
- write_sf(out_gages$segs, temp_gpkg, nsegments_layer)
- write_sf(nhd_Final, temp_gpkg, nhd_flowline)
+ # Streamgage to term outlet
+ moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_Term")
+ if(!is.data.frame(moved_pois)){
+ final_POIs <- moved_pois$final_pois
+ moved_pois_table <- moved_pois_table %>%
+ rbind(moved_pois$moved_points %>%
+ mutate(move_type = "gages to term"))
+ } else {
+ final_POIs <- moved_pois
+ }
- nsegments <- out_gages$segs
- final_POIs <- out_gages$allPOIs
-} else {
- final_POIs <- read_sf(nav_gpkg, final_poi_layer)
- nhd_Final <- read_sf(nav_gpkg, nhd_flowline)
- nsegments <- read_sf(temp_gpkg, nsegments_layer)
-}
+ # NID to waterbody outlet
+ moved_pois <- poi_move(final_POIs, "Type_NID", .025, 1, "Type_WBOut")
+ if(!is.data.frame(moved_pois)){
+ final_POIs <- moved_pois$final_pois
+ moved_pois_table <- moved_pois_table %>%
+ rbind(moved_pois$moved_points %>%
+ mutate(move_type = "nid to wb_out"))
+ } else {
+ final_POIs <- moved_pois
+ }
+
+ write_sf(final_POIs, nav_gpkg, pois_all_layer)
+ write_sf(moved_pois_table, temp_gpkg, "pois_collapsed")
+}
check_dups <- final_POIs %>%
group_by(COMID) %>%
@@ -661,7 +706,7 @@ check_dups <- final_POIs %>%
if(nrow(filter(check_dups, all(c(0,1) %in% nexus))) != nrow(check_dups)){
print("Multiple POI ids at same geometric location")
no_dups <- filter(check_dups, all(c(0,1) %in% nexus))
- dups <- filter(check_dups, !identifier %in% no_dups$identifier)
+ dups <- filter(check_dups, !id %in% no_dups$id)
write_sf(dups, temp_gpkg, dup_pois)
} else {
print("All double COMIDs nexus for gage or WB splitting")
@@ -672,19 +717,89 @@ 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
+final_POI_geom <- POIs %>%
+ select(identifier) %>%
+ cbind(st_coordinates(.)) %>%
+ group_by(X, Y) %>%
+ mutate(geom_id = cur_group_id()) %>%
+ ungroup()
+
+final_POIs_table <- POIs %>%
+ inner_join(select(st_drop_geometry(final_POI_geom), -X, -Y), by = "identifier") %>%
+ select(-identifier)
+
+# POI data theme table
+pois_data_orig <- reshape2::melt(st_drop_geometry(select(final_POIs_table,
+ -nexus)),
+ id.vars = c("COMID", "geom_id")) %>%
+ filter(!is.na(value)) %>%
+ group_by(COMID, geom_id) %>%
+ mutate(identifier = cur_group_id()) %>%
+ rename(hy_id = COMID, poi_id = identifier, hl_reference = variable, hl_link = value) %>%
+ distinct()
+
+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"))
+
+
+# POI Geometry table
+poi_geometry <- select(final_POIs_table, hy_id = COMID, geom_id) %>%
+ inner_join(distinct(pois_data, hy_id, geom_id, poi_id),
+ by = c("hy_id" = "hy_id", "geom_id" = "geom_id")) %>%
+ distinct() %>%
+ st_as_sf()
+
+write_sf(pois_data, nav_gpkg, poi_data_table)
+write_sf(poi_geometry, nav_gpkg, poi_geometry_table)
+
+poi_geom_xy <- cbind(poi_geometry, st_coordinates(poi_geometry)) %>%
+ st_drop_geometry()
+
+events_data <- events %>%
+ arrange(COMID) %>%
+ cbind(st_coordinates(.)) %>%
+ st_drop_geometry() %>%
+ group_by(COMID, REACHCODE, REACH_meas) %>%
+ mutate(event_id = cur_group_id()) %>%
+ rename(hy_id = COMID) %>%
+ ungroup()
+
+nexi <- filter(final_POIs_table, nexus == 1) %>%
+ cbind(st_coordinates(.)) %>%
+ select(hy_id = COMID, X, Y) %>%
+ inner_join(poi_geom_xy, by = c("hy_id" = "hy_id", "X" = "X", "Y" = "Y")) %>%
+ inner_join(events_data, by = c("hy_id" = "hy_id", "X" = "X", "Y" = "Y"), multiple = "all") %>%
+ select(hy_id, REACHCODE, REACH_meas, event_id, poi_id) %>%
+ group_by(hy_id, REACHCODE) %>%
+ filter(REACH_meas == min(REACH_meas)) %>%
+ ungroup()
+ #distinct(hy_id, REACHCODE, REACH_meas, event_id, poi_id)
+
+write_sf(select(events_data, -c(nexus, X, Y)), nav_gpkg, event_table)
+write_sf(nexi, nav_gpkg, event_geometry_table)
+
# Load data
if(needs_layer(nav_gpkg, lookup_table_refactor)) {
full_cats <- readRDS(data_paths$fullcats_table)
-
- lookup <- dplyr::select(sf::st_drop_geometry(nhd_Final),
- NHDPlusV2_COMID = COMID,
+
+ lookup <- dplyr::select(sf::st_drop_geometry(nhd_Final),
+ NHDPlusV2_COMID = COMID,
realized_catchmentID = COMID,
mainstem = LevelPathI) %>%
dplyr::mutate(realized_catchmentID = ifelse(realized_catchmentID %in% full_cats$FEATUREID,
realized_catchmentID, NA)) %>%
- left_join(nexus_from_poi(dplyr::rename(final_POIs, ID = COMID)), by = c("NHDPlusV2_COMID" = "ID"))
-
+ left_join(select(st_drop_geometry(poi_geometry), hy_id, poi_geom_id = geom_id),
+ by = c("NHDPlusV2_COMID" = "hy_id"))
+
sf::write_sf(lookup, nav_gpkg, lookup_table_refactor)
}
-
```
diff --git a/workspace/05_hyRefactor_flines.Rmd b/workspace/05_hyRefactor_flines.Rmd
index 4049aac67e94599faeed24f0ab9b60702f22d889..ccdb4274a4dc8233f31860becd63b3769f340893 100644
--- a/workspace/05_hyRefactor_flines.Rmd
+++ b/workspace/05_hyRefactor_flines.Rmd
@@ -28,36 +28,33 @@ source("R/hyRefactor_funs.R")
source("R/config.R")
```
-Load and filter source NHD Flowline network.
+Load and filter source NHD Flowline network. Bring in POIs
```{r flowlines}
nhd <- read_sf(out_refac_gpkg, nhd_flowline)
cats <- read_sf(out_refac_gpkg, nhd_catchment)
-POIs <- read_sf(nav_gpkg, final_poi_layer)
+POIs_data <- read_sf(nav_gpkg, poi_data_table)
-POI_id <- st_drop_geometry(POIs)
+POIs <- read_sf(nav_gpkg, poi_geometry_table)
-events <- read_sf(temp_gpkg, split_layer)
+POI_id <- st_drop_geometry(POIs)
```
Read POIs and add nhd outlets. Save to a non-spatial table.
```{r refactor}
-POIs <- POIs %>%
- inner_join(select(st_drop_geometry(nhd), TotDASqKM, COMID, DnHydroseq), by = "COMID")
-
-# Gages that were collapsed during navigate; might not correspond to their Gage_info COMIDs;
-# We know DAR for those are acceptable, so will keep those out of the events generation.
-POIs_mv <- filter(POIs, snapped)
+POIs_ref <- POIs %>%
+ inner_join(select(st_drop_geometry(nhd), TotDASqKM, COMID, DnHydroseq), by = c("hy_id" = "COMID"))
# Also need to avoid modification to flowlines immediately downstream of POIs
# This can cause some hydrologically-incorrect catchment aggregation
-POI_downstream <- filter(nhd, Hydroseq %in% POIs$DnHydroseq, AreaSqKM > 0)
+POI_downstream <- filter(nhd, Hydroseq %in% POIs_ref$DnHydroseq, AreaSqKM > 0)
# build final outlets set
-outlets <- POIs %>%
- mutate(type = ifelse(!is.na(Type_Term), "terminal", "outlet"))
-event_POIs <- filter(POIs, nexus == "1")
+term_poi <- filter(POIs_data, hl_reference == "Type_Term")
+
+outlets <- POIs_ref %>%
+ mutate(type = ifelse(poi_id %in% term_poi$poi_id, "terminal", "outlet"))
# derive list of unique terminal paths
TerminalPaths <- unique(nhd$TerminalPa)
@@ -74,7 +71,19 @@ write_sf(nhd, out_refac_gpkg, nhd_flowline)
nhdplus_flines <- sf::st_zm(filter(nhd, refactor == 1)) %>%
st_as_sf()
-events <- filter(events, COMID %in% nhdplus_flines$COMID)
+if(!needs_layer(nav_gpkg, event_geometry_table)){
+ events <- read_sf(nav_gpkg, event_geometry_table) %>%
+ mutate(event_identifier = as.character(row_number()))
+
+ events_ref <- filter(events, hy_id %in% nhdplus_flines$COMID) %>%
+ rename(COMID = hy_id) %>%
+ distinct(COMID, REACHCODE, REACH_meas, event_identifier)
+
+ outlets <- filter(outlets, !poi_id %in% events$poi_id) %>%
+ rename(COMID = hy_id)
+} else {
+ events_ref <- NULL
+}
if(needs_layer(out_refac_gpkg, refactored_layer)) {
@@ -93,7 +102,7 @@ if(needs_layer(out_refac_gpkg, refactored_layer)) {
out_reconciled = tr,
three_pass = TRUE,
purge_non_dendritic = FALSE,
- events = events,
+ events = events_ref,
exclude_cats = c(outlets$COMID, avoid$COMID, POI_downstream$COMID),
warn = TRUE)
@@ -129,49 +138,56 @@ if(needs_layer(out_refac_gpkg, outlets_layer)) {
select(member_COMID = COMID, Hydroseq, event_identifier, event_REACHCODE) %>%
inner_join(select(st_drop_geometry(nhd), orig_COMID = COMID, Hydroseq), by = "Hydroseq")
- # Subset for events
- refactored_events <- refactored %>%
- filter(!is.na(event_REACHCODE), !is.na(event_identifier))
+ if(!is.null(events_ref)){
+ # Subset for events
+ refactored_events <- refactored %>%
+ filter(!is.na(event_REACHCODE), !is.na(event_identifier))
- # Get ref_COMID for events
- # events_ref_COMID <- mutate(events, event_identifier = as.character(row_number())) %>%
- # left_join(select(st_drop_geometry(refactored_events), ref_COMID, event_identifier, orig_COMID),
- # by = "event_identifier")
-
- outlet_events <- filter(outlets, nexus == "1") %>%
- left_join(select(st_drop_geometry(refactored_events), member_COMID, event_identifier, orig_COMID),
- by = c("COMID" = "orig_COMID")) %>%
- filter(!is.na(member_COMID)) %>%
- select(-event_identifier)
-
- # subset for refactored outlets (non-events)
- refactored_outlets <- filter(refactored, !member_COMID %in% outlet_events$member_COMID)
-
- # get ref_COMId for other outlets
- outlets_ref_COMID <- filter(outlets, !identifier %in% outlet_events$identifier) %>%
- left_join(select(st_drop_geometry(refactored_outlets), member_COMID, orig_COMID),
- by = c("COMID" = "orig_COMID")) %>%
- group_by(COMID) %>%
- filter(member_COMID == max(member_COMID)) %>%
- rbind(outlet_events) %>%
+ event_outlets <- events %>%
+ inner_join(st_drop_geometry(refactored_events), by = "event_identifier") %>%
+ select(hy_id, event_id, poi_id, member_COMID)
+
+ # subset for refactored outlets (non-events)
+ refactored_outlets <- filter(refactored, !member_COMID %in% event_outlets$member_COMID)
+
+ # get ref_COMId for other outlets
+ outlets_ref <- outlets %>%
+ left_join(select(st_drop_geometry(refactored_outlets), member_COMID, orig_COMID),
+ by = c("COMID" = "orig_COMID")) %>%
+ group_by(COMID) %>%
+ filter(member_COMID == max(member_COMID)) %>%
+ select(hy_id = COMID, poi_id, member_COMID, type)
+
+ outlets_ref_COMID <- data.table::rbindlist(list(outlets_ref, event_outlets), fill = TRUE) %>%
+ st_as_sf()
+ } else {
+ # get ref_COMId for other outlets
+ outlets_ref_COMID <- outlets %>%
+ left_join(select(st_drop_geometry(refactored), member_COMID, orig_COMID),
+ by = c("COMID" = "orig_COMID")) %>%
+ group_by(COMID) %>%
+ filter(member_COMID == max(member_COMID)) %>%
+ select(hy_id = COMID, poi_id, member_COMID, type)
+ }
+
+ final_outlets <- outlets_ref_COMID %>%
+ st_as_sf() %>%
inner_join(select(lookup_table, member_COMID, reconciled_ID),
- by = "member_COMID")
+ by = "member_COMID")
- write_sf(outlets_ref_COMID, out_refac_gpkg, outlets_layer)
-
-
+ write_sf(final_outlets, out_refac_gpkg, outlets_layer)
} else {
- outlets_ref_COMID <- read_sf(out_refac_gpkg, outlets_layer)
+ final_outlets <- read_sf(out_refac_gpkg, outlets_layer)
}
-check_dups <- outlets_ref_COMID %>%
+check_dups_poi <- final_outlets %>%
group_by(reconciled_ID) %>%
filter(n() > 1) %>%
- ungroup()
+ ungroup
-if(nrow(check_dups) > 1){
+if(nrow(check_dups_poi) > 1){
print("Double-check for double POIs")
- write_sf(check_dups, temp_gpkg, paste0(dup_pois, "_", rpu_code))
+ write_sf(check_dups_poi, out_refac_gpkg, paste0(dup_pois, "_", rpu_code))
} else {
print("no double POIs detected")
}
diff --git a/workspace/06-1_hyRefactor_cats.Rmd b/workspace/06-1_hyRefactor_cats.Rmd
index 58c02af5c1c4fd64bc15b017c7809d9151fa2ec8..3071f1e398c59594f8f0e01ab2ce2c38aa0b04e4 100644
--- a/workspace/06-1_hyRefactor_cats.Rmd
+++ b/workspace/06-1_hyRefactor_cats.Rmd
@@ -24,6 +24,7 @@ source("R/utils.R")
source("R/config.R")
source("R/NHD_navigate.R")
source("R/hyRefactor_funs.R")
+
```
This reconcile-catchments step splits and combines catchment boundaries according to the process that was completed in refactor flines.
@@ -33,13 +34,12 @@ if(needs_layer(out_refac_gpkg, divide_layer)) {
rpu_string <- paste0("rpu_", rpu_code)
fdr <- data_paths$fdr[[rpu_string]]
fac <- data_paths$fac[[rpu_string]]
-
+
fdr_temp <- terra::rast(fdr)
raster_crs <- terra::crs(fdr_temp)
# reconciled have no ID applied in refactoring flowlines
- reconciled <- st_transform(read_sf(out_refac_gpkg, reconciled_layer), raster_crs)
-
+ reconciled <- st_transform(read_sf(out_refac_gpkg, reconciled_layer), raster_crs)
# refactored has the original COMIDs and other attributes
refactored <- st_transform(read_sf(out_refac_gpkg, refactored_layer), raster_crs) %>%
distinct(.keep_all = T)
@@ -51,25 +51,23 @@ if(needs_layer(out_refac_gpkg, divide_layer)) {
fline_rec = reconciled,
fdr = fdr,
fac = fac,
- para = para_reconcile,
+ para = para_reconcile,
cache = cache_split,
- keep = NULL)
-
+ keep = NULL)
if(exists("divides")) {
-
+
load(cache_split)
-
- split_cats <- bind_rows(split_cats)
-
+
+ split_cats <- bind_rows(split_cats[lengths(split_cats) > 0])
+
split_cats$areasqkm <- as.numeric(units::set_units(sf::st_area(split_cats), "km^2"))
-
+
write_sf(sf::st_transform(split_cats, crs), out_refac_gpkg, split_divide_layer)
-
- unlink(cache_split)
-
+
}
+ unlink(cache_split)
write_sf(sf::st_transform(divides, crs), out_refac_gpkg, divide_layer)
}
diff --git a/workspace/06-2_aggregate_cats.Rmd b/workspace/06-2_aggregate_cats.Rmd
index d3cbaf52ef61bd6e25c3eb29daa7f85cd8ed2183..ba68264e01bd1957653f18e7f0f082bb10da000f 100644
--- a/workspace/06-2_aggregate_cats.Rmd
+++ b/workspace/06-2_aggregate_cats.Rmd
@@ -24,6 +24,7 @@ source("R/utils.R")
source("R/config.R")
source("R/NHD_navigate.R")
source("R/hyRefactor_funs.R")
+
```
Derive set of outlets to aggregate catchments for. There are three sets of
@@ -37,39 +38,47 @@ nhd <- read_sf(out_refac_gpkg, nhd_flowline)
# Load outlets
outlets_POI <- read_sf(out_refac_gpkg, outlets_layer) %>%
- distinct()
+ distinct() %>%
+ filter(hy_id %in% nhd$COMID)
+
+POIs_data <- read_sf(nav_gpkg, poi_data_table)
# reconciled have no ID applied in refactoring flowlines
reconciled <- st_transform(read_sf(out_refac_gpkg, reconciled_layer), crs)
divides <- st_transform(read_sf(out_refac_gpkg, divide_layer), crs)
+
+lookup <- read_sf(out_refac_gpkg, lookup_table_refactor) %>%
+ inner_join(select(st_drop_geometry(nhd), COMID, TotDASqKM), by = c("NHDPlusV2_COMID" = "COMID"))
+
+highest_DA <- lookup %>%
+ group_by(reconciled_ID) %>%
+ filter(TotDASqKM == max(TotDASqKM))
```
```{r agg_catchments I}
if(needs_layer(out_agg_gpkg, agg_cats_layer)){
# https://github.com/r-spatial/sf/issues/1094#issuecomment-988933580 for why this is here.
- sf::st_geometry(divides[sf::st_is_empty(divides),]) <-
+ sf::st_geometry(divides[sf::st_is_empty(divides),]) <-
sf::st_cast(sf::st_geometry(divides[sf::st_is_empty(divides),]), "MULTIPOLYGON")
-
- #divides <- hyRefactor::clean_geometry(divides)
# POI outlets
outlets <- select(st_drop_geometry(outlets_POI), ID = reconciled_ID, type)
# Identify reconciled flowpaths that are terminal and above DA threshold
- reconciled_DA <- filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm, is.na(toID))
-
+ reconciled_DA <- filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm, is.na(toID))
+
# Ones not in existing pois
missing_outlets <- filter(reconciled_DA, !ID %in% outlets$ID) %>%
pull(ID)
-
+
# bind missing outlets, re-assign mis-attributed terminals
if(length(missing_outlets) > 0){
outlets <- bind_rows(outlets,
data.frame(ID = missing_outlets,
type = rep("terminal", length(missing_outlets))))
write_sf(filter(reconciled, ID %in% missing_outlets), out_agg_gpkg, "missing_outlets")
-
+
outlets <- outlets %>%
left_join(select(st_drop_geometry(reconciled), ID, toID), by = "ID") %>%
mutate(type = ifelse(type == "terminal" & !is.na(toID), "outlet", type)) %>%
@@ -83,10 +92,10 @@ if(needs_layer(out_agg_gpkg, agg_cats_layer)){
da_thresh = aggregate_da_thresh_sqkm,
only_larger = TRUE,
post_mortem_file = cache_split)
-
+
agg_cats$cat_sets$set <- sapply(agg_cats$cat_sets$set, paste, collapse = ",")
agg_cats$fline_sets$set <- sapply(agg_cats$fline_sets$set, paste, collapse = ",")
-
+
write_sf(agg_cats$cat_sets, out_agg_gpkg, agg_cats_layer)
write_sf(agg_cats$fline_sets, out_agg_gpkg, agg_fline_layer)
@@ -96,18 +105,34 @@ if(needs_layer(out_agg_gpkg, agg_cats_layer)){
st_as_sf()
# Build final POI set
- final_POIs <- st_drop_geometry(agg_cats$fline_set) %>%
+ POIs_att <- select(rec_outlets, ID, toID) %>%
# Bring over POI information
- left_join(st_drop_geometry(outlets_POI), by = c("ID" = "reconciled_ID")) %>%
- # Bring over physical end node XY
- left_join(select(rec_outlets, ID), by = "ID") %>%
- # Mark additional POI/outlets created during aggregate cats
- mutate(Type_Con = ifelse(!ID %in% outlets$ID, 1, 0)) %>%
+ inner_join(st_drop_geometry(outlets_POI), by = c("ID" = "reconciled_ID")) %>%
+ select(-c(type, member_COMID))
+
+ POIs_missing <- filter(select(rec_outlets, -toID), !ID %in% POIs_att$ID) %>%
+ inner_join(st_drop_geometry(agg_cats$fline_sets), by = "ID") %>%
+ arrange(ID) %>%
+ mutate(poi_id = max(POIs_att$poi_id) + row_number()) %>%
+ select(ID, toID, poi_id, event_id = event_identifier, Type_Term = orig_levelpathID) %>%
+ inner_join(select(highest_DA, NHDPlusV2_COMID, reconciled_ID), by = c("ID" = "reconciled_ID")) %>%
+ select(ID, toID, hy_id = NHDPlusV2_COMID, poi_id, event_id, Type_Term) %>%
+ mutate(Type_Term = ifelse(is.na(toID), Type_Term, NA),
+ Type_Con = ifelse(!is.na(toID), 1, NA)) %>%
+ st_as_sf()
+
+ final_POIs <- data.table::rbindlist(list(POIs_att, select(POIs_missing,
+ -c(Type_Term, Type_Con))), fill = TRUE) %>%
st_as_sf()
+
+ unaggregated_outlets <- filter(outlets_POI, !reconciled_ID %in% final_POIs$ID)
+ double_outlets <- final_POIs %>% group_by(ID) %>% filter(n() > 1)
+
# Write out
write_sf(final_POIs, out_agg_gpkg, mapped_outlets_layer)
-
+ write_sf(unaggregated_outlets, out_agg_gpkg, "unmapped_outlets")
+ write_sf(double_outlets, out_agg_gpkg, "double_outlets")
} else {
agg_cats <- list(cat_sets = read_sf(out_agg_gpkg, agg_cats_layer),
fline_sets = read_sf(out_agg_gpkg, agg_fline_layer))
@@ -117,6 +142,25 @@ if(needs_layer(out_agg_gpkg, agg_cats_layer)){
```
+```{r Long form POIs}
+POIs <- final_POIs %>%
+ arrange(ID)
+
+final_POI_geom <- POIs %>%
+ select(ID) %>%
+ cbind(st_coordinates(.)) %>%
+ group_by(ID, X, Y) %>%
+ mutate(geom_id = cur_group_id()) %>%
+ ungroup()
+
+final_POIs_table <- POIs %>%
+ inner_join(select(st_drop_geometry(final_POI_geom), -X, -Y), by = "ID") %>%
+ distinct()
+
+write_sf(final_POIs_table, out_agg_gpkg, poi_data_table)
+```
+
+
```{r lookup table}
refactor_lookup <- dplyr::select(st_drop_geometry(reconciled), ID, member_COMID) %>%
dplyr::mutate(member_COMID = strsplit(member_COMID, ",")) %>%
diff --git a/workspace/07-1_merge.Rmd b/workspace/07-1_merge.Rmd
index 827742da43ba6612b6b827464fa4942c25ece4b2..5a8202306c61905921e72b1ad8c762a7120ee402 100644
--- a/workspace/07-1_merge.Rmd
+++ b/workspace/07-1_merge.Rmd
@@ -28,18 +28,17 @@ source("R/config.R")
For a given VPU bind RPUs into a single geopackage for an entire VPU. This step is intended to be run after Refactor flowlines/cats. For a given VPU, this binds the results of aggregated from the associated RPUs into a single layer and writes them out to a single geopackage for the non-dendritic work.
-```{r Bind RPUs}
-if(needs_layer(gf_gpkg, agg_cats_layer)) {
+```{r refactor}
+if(needs_layer(rfc_gpkg, divide_layer)) {
# Thematic POIs
- POIs <- read_sf(nav_gpkg, final_poi_layer) %>%
- st_transform(crs)
+ POIs <- read_sf(nav_gpkg, poi_data_table) %>%
+ select(-geom_id)
- # write_sf(POIs, rfc_gpkg, final_poi_layer)
- #
- # write_sf(POIs, gf_gpkg, final_poi_layer)
- #
- merged_layers <- merge_refactor(rpu_codes$rpuid, rpu_vpu_out,
+ sf::write_sf(POIs, gf_gpkg, poi_data_table)
+
+ merged_layers <- merge_refactor(rpu_codes$rpuid,
+ rpu_vpu_out,
lookup_table_refactor,
reconciled_layer,
divide_layer,
@@ -61,9 +60,7 @@ if(needs_layer(gf_gpkg, agg_cats_layer)) {
sf::write_sf(merged_layers[[split_divide_layer]],
rfc_gpkg, split_divide_layer)
-
- sf::write_sf(merged_layers[[mapped_outlets_layer]], rfc_gpkg, mapped_outlets_layer)
-
+
sf::write_sf(merged_layers[[mapped_outlets_layer]], gf_gpkg, mapped_outlets_layer)
merged_layers[[agg_cats_layer]] <- merged_layers[[agg_cats_layer]] %>%
@@ -93,6 +90,5 @@ if(needs_layer(gf_gpkg, agg_cats_layer)) {
id = aggregated_flowpath_ID,
levelpathid = LevelPathID))) %>%
sf::write_sf(gf_gpkg, catchment_network_table)
-
}
```
diff --git a/workspace/R/NHD_navigate.R b/workspace/R/NHD_navigate.R
index dfeaba0823f9f1b3395cc7e41ea602783e8b7025..02f8f8c4da6816590772ca63d32d9dcaff64f5d5 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){
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) == FALSE){
+ if (!missing(routing_fix)){
routing_fix <- routing_fix %>%
rename(COMID = oldPOI, new_COMID = COMID)
@@ -207,19 +207,22 @@ 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)
+ nhdDF <- distinct(nhdDF)
# Closest POI/US/DS
- up_segs <- nhdplusTools::get_UM(nhdDF, poi_fix, sort=T)
- seg2fix <- filter(nhdDF, COMID == poi_fix)
+ up_segs <- unique(nhdplusTools::get_UM(nhdDF, poi_fix, sort=T))
+ seg2fix <- filter(nhdDF, COMID == poi_fix) %>%
+ distinct()
# Sorted results and filter out all flowlines w/o catchments
upstuff <- filter(nhdDF, COMID %in% unlist(up_segs)) %>%
- arrange(factor(COMID, levels = up_segs)) %>%
+ arrange(Hydroseq) %>%
filter(AreaSqKM > 0)
- down_segs <- nhdplusTools::get_DM(nhdDF, poi_fix, sort=T)
+ down_segs <- unique(nhdplusTools::get_DM(nhdDF, poi_fix, sort=T))
downstuff <- filter(nhdDF, COMID %in% unlist(down_segs)) %>%
- arrange(factor(COMID, levels = down_segs)) %>%
+ arrange(Hydroseq)%>%
filter(AreaSqKM > 0)
# combine into one dataframe, select up/downstream seg with least change in total drainage area
@@ -966,11 +969,16 @@ gage_POI_creation <- function(tmp_POIs, gages_info, nhd, combine_meters, reach_m
gage_POIs_nonevent <- filter(gage_POIs, !Type_Gages %in% events$Type_Gages) %>%
addType(., tmp_POIs, "Gages", nexus = FALSE, bind = TRUE)
- tmp_POIs <- data.table::rbindlist(list(gage_POIs_nonevent,
- select(events, COMID, Type_Gages, nexus)), fill = TRUE) %>%
- mutate(nexus = ifelse(is.na(nexus), FALSE, nexus)) %>%
- st_as_sf()
-
+ if(nrow(events) > 0){
+ tmp_POIs <- data.table::rbindlist(list(gage_POIs_nonevent,
+ select(events, COMID, Type_Gages, nexus)), fill = TRUE) %>%
+ mutate(nexus = ifelse(is.na(nexus), FALSE, nexus)) %>%
+ st_as_sf()
+ } else {
+ events <- NA
+ tmp_POIs <- mutate(tmp_POIs, nexus = FALSE)
+ }
+
return(list(events = events, tmp_POIs = tmp_POIs))
}
@@ -985,7 +993,7 @@ gage_POI_creation <- function(tmp_POIs, gages_info, nhd, combine_meters, reach_m
#'
wbout_POI_creaton <- function(nhd, WBs_VPU, data_paths, crs){
# Create waterbody outlet POIs for waterbodies that are in NHDv2 waterbody set
- wbout_COMIDs <- filter(nhd, dend == 1 & WB == 1) %>%
+ wbout_COMIDs <- nhd %>%
group_by(WBAREACOMI) %>%
slice(which.min(Hydroseq)) %>%
switchDiv(., nhd) %>%
@@ -1074,7 +1082,7 @@ wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs,
filter(minNet == 1)
# Headwater Waterbodies that may need the network extended through the inlet
- need_wbin <- st_drop_geometry(filter(WBs_VPU, source == "NHDv2WB")) %>%
+ need_wbin <- st_drop_geometry(filter(WBs_VPU, source %in% c("ref_WB", "NHDv2WB"))) %>%
dplyr::select(COMID)%>%
dplyr::filter(!COMID %in% wbin_COMIDs$WBAREACOMI)
@@ -1208,7 +1216,7 @@ wb_inlet_collapse <- function(tmp_POIs, nhd, events){
unique()
if(nrow(gage_reach) == 0){
- print("Wilton")
+ print("no gage reaches")
}
if(nrow(gage_event) == 0){
@@ -1327,17 +1335,15 @@ wb_inlet_collapse <- function(tmp_POIs, nhd, events){
}
}
-#' Creates Waterbody POIs, calls a few other functions
+#' Collapses POIs us/ds of waterbody POIs
#' @param (sf data.frame) rolling POI data frame
#' @param nhd (sf data.frame) nhd flowlines
-#' @param events (sf data.frame) waterbody inlet events
-#'
-#' @return (sf data.frame) dataframe of wb inlet POIs collapse
+#' @param events (sf data.frame) waterbody events
#'
-#' wb_poi_collapse <- function(tmp_POIs, nhd, events){
+#' @return (sf data.frame) dataframe of wb inlet POIs collapsed
wb_poi_collapse <- function(tmp_POIs, nhd, events){
gage_dist_node <- function(x, wb_ds_ds, gage_add, events){
- print (x)
+ print (x)
wb_out_fl <- filter(wb_ds_ds, COMID == x)
gage_ds <- filter(wb_ds_ds, Hydroseq %in% wb_out_fl$Hydroseq |
Hydroseq %in% wb_out_fl$DnHydroseq)
@@ -1353,15 +1359,14 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
unique()
if(nrow(gage_reach) == 0){
- print("Wilton")
+ print("no gages")
}
if(nrow(gage_event) == 0){
- #print("akermayun")
- return("Akermayun")
+ return("no events")
} else if(gage_event$COMID != wb_out_fl$COMID) {
gage_reach <- gage_reach %>%
- filter(REACHCODE == gage_event$reachcode) %>%
+ filter(REACHCODE == unique(gage_event$reachcode)) %>%
mutate(gage_dist = ifelse(gage_event$nexus == TRUE,
total_length * (1 - (gage_event$reach_meas/100)),
total_length)) %>%
@@ -1369,10 +1374,10 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
wbout_comid = x)
} else if(gage_event$COMID == wb_out_fl$COMID){
if(nrow(wb_event) >0){
- wb_out_meas <- wb_event$REACH_meas
+ wb_out_meas <- min(wb_event$REACH_meas)
wb_RC <- wb_event$REACHCODE
} else {
- wb_out_meas <- wb_out_fl$FromMeas
+ wb_out_meas <- min(wb_out_fl$FromMeas)
wb_RC <- wb_out_fl$REACHCODE
}
@@ -1383,7 +1388,7 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
# wb info
wb_reach <- gage_reach %>%
- filter(REACHCODE == wb_RC) %>%
+ filter(REACHCODE == unique(wb_RC)) %>%
mutate(wb_dist = total_length * (1 - (wb_out_meas /100)))
gage_reach <- gage_reach %>%
@@ -1393,14 +1398,10 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
}
}
- #events <- read_sf(temp_gpkg, split_layer) %>%
- # rbind(st_compatibalize(wb_,.))
-
# Previously identified streamgages within Gage_Selection.Rmd
streamgages_VPU <- gages %>%
rename(COMID = comid) %>%
filter(COMID %in% nhd$COMID) %>%
- #st_drop_geometry() %>%
switchDiv(., nhd)
# get waterbody outlets
@@ -1421,10 +1422,16 @@ wb_poi_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") %>%
+ distinct()
output <- lapply(wb_out$COMID,
gage_dist_node, wb_ds_ds, gage_add, events)
+ output_length <- output[lengths(output) > 1]
+
+ if(length(output_length) == 0){
+ return(list(POIs = tmp_POIs, events_ret = NA))
+ }
output_full <- do.call("rbind", output[lengths(output) > 1]) %>%
filter(gage_dist < 1) %>%
@@ -1438,12 +1445,11 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
gage_POI <- filter(tmp_POIs, COMID %in% output_full$gage_comid) %>%
select(COMID, Type_HUC12_ds = Type_HUC12, Type_Gages_ds = Type_Gages,
- Type_TE_ds = Type_TE, Type_Term_ds = Type_Term, nexus) %>%
+ Type_TE_ds = Type_TE, nexus) %>%
st_drop_geometry() %>%
group_by(COMID) %>%
summarise(Type_HUC12_ds = last(na.omit(Type_HUC12_ds)),
Type_Gages_ds = last(na.omit(Type_Gages_ds)),
- Type_Term_ds = last(na.omit(Type_Term_ds)),
Type_TE_ds = last(na.omit(Type_TE_ds)),
nexus = last(na.omit(nexus)))
@@ -1452,9 +1458,8 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
inner_join(select(gage_POI, -nexus), by = c("gage_comid" = "COMID")) %>%
mutate(Type_HUC12 = ifelse(!is.na(Type_HUC12_ds), Type_HUC12_ds, Type_HUC12),
Type_Gages = ifelse(!is.na(Type_Gages_ds), Type_Gages_ds, Type_Gages),
- Type_TE = ifelse(!is.na(Type_TE_ds), Type_TE_ds, Type_TE),
- Type_Term = ifelse(!is.na(Type_Term_ds), Type_Term_ds, Type_Term)) %>%
- select(-c(Type_HUC12_ds, Type_Gages_ds, Type_TE_ds, Type_Term_ds))
+ Type_TE = ifelse(!is.na(Type_TE_ds), Type_TE_ds, Type_TE)) %>%
+ select(-c(Type_HUC12_ds, Type_Gages_ds, Type_TE_ds))
tmp_POIs_fin <- filter(tmp_POIs, !COMID %in% c(WB_POI$COMID, WB_POI$gage_comid)) %>%
rbind(select(WB_POI, -gage_comid))
@@ -1466,4 +1471,135 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
return(list(POIs = tmp_POIs_fin, events_ret = events))
+}
+
+#' Collapses POIs together based on criteria
+#' @param (pois) sf data frame of POIs
+#' @param move_category (character) POI data theme to move
+#' @param DAR (numeric) drainage area threshold to move within
+#' @param dist (numeric) maximum river distance between two points to move within
+#' @param keep_category (character) POI data themes to keep static
+#'
+#' @return (sf data.frame, table) dataframe of pois, table of points that have moved
+poi_move <- function(pois, move_category, DAR, dist, keep_category) {
+ # filter out features with identical geometry
+
+ # Add row_number
+ pois_edit <- pois %>%
+ mutate(nexus = ifelse(is.na(nexus), 0, nexus))
+
+ # Don't consider points already moved
+ if("moved" %in% colnames(pois_edit)){
+ pois_tomove <- filter(pois_edit, is.na(moved)) # change from poi_edit
+ pois_moved_pre <- filter(pois_edit, !is.na(moved))}
+
+ # If 'keep' category included
+ if(!missing(keep_category)){
+ poi2move <- filter(pois_tomove, !is.na(.data[[move_category]]) & nexus == FALSE) %>%
+ filter(if_all(!!as.symbol(keep_category), function(x) is.na(x))) %>%
+ # Never move these
+ filter_at(vars(Type_WBOut, Type_WBIn, Type_Conf, Type_Term), all_vars(is.na(.)))
+
+ pois2keep <- filter(pois_tomove, !id %in% poi2move$id)
+ #is.na(.data[[move_category]]) & nexus == FALSE) #%>%
+ #filter(if_all(!!as.symbol(keep_category), function(x) is.na(x)))
+ } else {
+ # POIs to move
+ poi2move <- pois_tomove %>%
+ filter_at(vars(Type_WBOut, Type_WBIn, Type_Conf, Type_Term), all_vars(is.na(.))) %>%
+ filter(nexus == 0) %>%
+ filter(!is.na(.data[[move_category]]))
+
+ pois2keep <- filter(pois_tomove, !id %in% poi2move$id)
+ }
+
+ # Get relevant NHD data
+ nhd_poi1 <- filter(nhd, COMID %in% pois2keep$COMID)
+ nhd_poi2 <- filter(nhd, COMID %in% poi2move$COMID)
+ # Ensure they are on same level path
+ nhd_poi2 <- filter(nhd_poi2, LevelPathI %in% nhd_poi1$LevelPathI)
+
+ # Join NHD data
+ pois2keep_nhd <- pois2keep %>%
+ inner_join(select(st_drop_geometry(nhd_poi1), COMID, LevelPathI, Hydroseq,
+ DA_keep = TotDASqKM, Pathlength_keep = Pathlength), by = "COMID") %>%
+ rename(COMID_keep = COMID)
+
+ # Join NHD data
+ pois2move_nhd <- select(poi2move, COMID, !!as.symbol(move_category), id_move = id) %>%
+ inner_join(select(st_drop_geometry(nhd_poi2), COMID, LevelPathI, Hydroseq, TotDASqKM, Pathlength),
+ by = "COMID")
+
+ # Candidates to move
+ pois2move_cand <-pois2move_nhd %>%
+ inner_join(select(st_drop_geometry(pois2keep_nhd), COMID_keep, DA_keep, LevelPathI,
+ Pathlength_keep, id_keep = id, nexus),
+ by = "LevelPathI") %>%
+ mutate(river_dist = abs(Pathlength - Pathlength_keep), DAR_poi = abs(DA_keep/TotDASqKM),
+ move_dir = ifelse(Pathlength < Pathlength_keep, "Up", "Down")) %>%
+ group_by(id_move, move_dir) %>%
+ ungroup() %>%
+ filter((river_dist < dist) & (DAR_poi > (1 - DAR)) & (DAR_poi < (1 + DAR))) %>%
+ select(!!as.symbol(move_category), id_move, COMID, id_keep, COMID_keep, river_dist, DAR_poi, move_dir, nexus) %>%
+ st_drop_geometry()
+
+ move_distinct <- pois2move_cand %>%
+ group_by(id_keep) %>%
+ filter(row_number() == 1) %>%
+ ungroup() %>%
+ distinct(id_move, COMID_move = COMID, id_keep, COMID_keep, river_dist, DAR_poi, move_dir, nexus) %>%
+ group_by(id_move) %>%
+ slice(which.min(abs(1 - DAR_poi)))
+
+ if(nrow(move_distinct) == 0){
+ print("no POIs to move")
+ return(pois)
+ }
+
+ pois2_move <- filter(st_drop_geometry(pois_tomove), id %in% move_distinct$id_move) %>%
+ select_if(~sum(!is.na(.)) > 0) %>%
+ select(-c(COMID, nexus)) %>%
+ inner_join(select(move_distinct, id_move, id_keep), by = c("id" = "id_move"))
+
+ move_fields <- colnames(select(pois2_move, -c(id, id_keep)))
+
+ if(length(move_fields) == 1){
+ pois2_keep <- filter(pois_tomove, id %in% pois2_move$id_keep, !id %in% pois2_move$id) %>%
+ inner_join(select(pois2_move, id_move = id, id_keep,
+ new_val = !!as.symbol(move_category)), by = c("id" = "id_keep")) %>%
+ mutate(moved := ifelse(is.na(!!as.symbol(move_category)),
+ id_move, moved),
+ !!as.symbol(move_category) := ifelse(is.na(!!as.symbol(move_category)),
+ new_val, !!as.symbol(move_category)))
+
+ moved_points <- filter(pois2_keep, !is.na(new_val), !is.na(moved)) %>%
+ mutate(moved_value = move_category)
+ } else {
+ for (field in move_fields){
+ pois2_keep <- filter(pois_tomove, id %in% pois2_move$id_keep, !id %in% pois2_move$id) %>%
+ inner_join(select(pois2_move, id_move = id, id_keep, new_val = !!as.symbol(field)),
+ by = c("id" = "id_keep")) %>%
+ mutate(moved := ifelse(is.na(!!as.symbol(field)),
+ id_move, moved),
+ !!as.symbol(field) := ifelse(is.na(!!as.symbol(field)),
+ new_val, !!as.symbol(field)))
+
+ pois_moved <- filter(pois2_keep, !is.na(new_val), !is.na(moved)) %>%
+ mutate(moved_value = field)
+
+ if(!exists("moved_points")){
+ moved_points <- pois_moved
+ } else {
+ moved_points <- rbind(moved_points, pois_moved)
+ }
+ }
+ }
+
+
+ pois_final <- data.table::rbindlist(list(filter(pois_edit, !id %in% c(moved_points$id_move, pois2_keep$id)),
+ select(pois2_keep, -c(new_val, id_move, new_val))), fill = TRUE) %>%
+ st_as_sf()
+
+ return(list(final_pois = pois_final, moved_points = moved_points))
+
}
\ No newline at end of file
diff --git a/workspace/R/config.R b/workspace/R/config.R
index 6ebf0344aeecc62f42ca75cb8adf0e321e06f2d0..4d06f0ce56da7096418a4a79e1b858c7f42a900a 100644
--- a/workspace/R/config.R
+++ b/workspace/R/config.R
@@ -13,24 +13,23 @@ rpu_vpu_out <- readr::read_csv("cache/rpu_vpu_out.csv",
col_types = c("c", "c", "i" , "i"), show_col_types = FALSE)
if(nchar(Sys.getenv("HYDREG")) > 1) {
- VPU <- Sys.getenv("HYDREG")
+ vpu_codes <- Sys.getenv("HYDREG")
}
if(!exists("rpu_code")) {
- if(exists("VPU")) {
- vpu <- VPU
- rpu_codes <- rpu_vpu[rpu_vpu$vpuid %in% vpu, ]
+ if(exists("vpu_codes")) {
+ rpu_codes <- rpu_vpu[rpu_vpu$vpuid %in% vpu_codes, ]
rpu_code <- rpu_codes$rpuid[1]
} else {
- vpu <- get_rpu_dependent_vars()
- rpu_code <- vpu$r
- vpu <- vpu$v
- rpu_codes <- rpu_vpu[rpu_vpu$vpuid %in% vpu, ]
+ vpu_codes <- get_rpu_dependent_vars()
+ rpu_code <- vpu_codes$r
+ vpu_codes <- vpu_codes$v
+ rpu_codes <- rpu_vpu[rpu_vpu$vpuid %in% vpu_codes, ]
}
}
-if(!exists("vpu")) {
- vpu <- get_rpu_dependent_vars(rpu_code)$v
+if(!exists("vpu_codes")) {
+ vpu_codes <- get_rpu_dependent_vars(rpu_code)$v
}
if(!nchar(out_agg_gpkg <- Sys.getenv("OUT_GPKG")) > 5)
@@ -46,7 +45,7 @@ if(!nchar(lookup_table_file <- Sys.getenv("LOOKUP_TABLE")) > 5)
lookup_table_file <- file.path("cache", paste0(rpu_code, "_lookup.csv"))
if(!nchar(vpu_lookup_table_file <- Sys.getenv("VPU_LOOKUP_TABLE")) > 5)
- vpu_lookup_table_file <- file.path("cache", paste0("lookup_", vpu, ".csv"))
+ vpu_lookup_table_file <- file.path("cache", paste0("lookup_", vpu_codes, ".csv"))
options(scipen = 9999)
@@ -61,6 +60,7 @@ suppressWarnings(staged_nhd <- stage_national_data())
# Defined and used broadly
nhd_flowline <- "reference_flowline"
nhd_catchment <- "reference_catchment"
+nhd_waterbody <- "reference_waterbody"
nhd_outlet <- "nhd_outlet"
nhd_network <- "reference_network"
@@ -86,16 +86,22 @@ poi_dar_move <- .05 # maximum proportional drainage area difference to collapse
# Intermediate layers created during 02_Navigate
xwalk_layer <- paste0("HUC12_nhd") # HUC12 - nhdcat crosswalk, built in Nav for VPU 20
-nav_poi_layer <- paste0("POIs_tmp_", vpu) # Rolling Nav POI layer added to/modified througout nav workflow
-WBs_layer <- paste0("WB_", vpu) # Waterbodies within VPU
-WB_events_layer <- paste0("WB_events_", vpu) # inlet and outlet events for NHDArea and HR waterbodies
-poi_moved_layer <- paste0("POIs_mv_", vpu) # POIs moved from original COMID assignment
-nsegments_layer <- paste0("nsegment_", vpu) # Minimally-sufficient network dissolved by POI_ID
-pois_all_layer <- paste0("POIs_", vpu) # All POIs binded together
-poi_xwalk_layer <- paste0("poi_xwalk_layer_", vpu) # POIs that changed COMIDS during the navigate part of the workflow
+nav_poi_layer <- paste0("POIs_tmp_", vpu_codes) # Rolling Nav POI layer added to/modified througout nav workflow
+WBs_layer <- paste0("WB_", vpu_codes) # Waterbodies within VPU
+WB_events_layer <- paste0("WB_events_", vpu_codes) # inlet and outlet events for NHDArea and HR waterbodies
+poi_moved_layer <- paste0("POIs_mv_", vpu_codes) # POIs moved from original COMID assignment
+nsegments_layer <- paste0("nsegment_", vpu_codes) # Minimally-sufficient network dissolved by POI_ID
+pois_all_layer <- paste0("POIs_", vpu_codes) # All POIs binded together
+poi_xwalk_layer <- paste0("poi_xwalk_layer_", vpu_codes) # POIs that changed COMIDS during the navigate part of the workflow
final_poi_layer <- "POIs"
dup_pois <- "dup_POIs"
+poi_data_table <- "poi_data"
+poi_geometry_table <- "poi_geometry_table"
+
+event_table <- "event_data"
+event_geometry_table <- "event_geometry_table"
+
# Settings for refactor workflow
split_meters <- 10000
combine_meters <- 1000
@@ -126,12 +132,13 @@ catchment_network_table <- "catchment_network"
CAC_thresh <- 0.66 # Coefficient of areal correspondence threshold
# output geopackage file names
-nav_gpkg <- file.path("cache", paste0("reference_", vpu,".gpkg"))
-temp_gpkg <- file.path("cache", paste0("nav_", vpu,".gpkg"))
+ref_gpkg <- file.path(data_paths$ref_fab_path, paste0(vpu_codes, "_reference_features.gpkg"))
+nav_gpkg <- file.path("cache", paste0("reference_", vpu_codes,".gpkg"))
+temp_gpkg <- file.path("cache", paste0("nav_", vpu_codes,".gpkg"))
-rfc_gpkg <- file.path("cache", paste0("refactor_", vpu, ".gpkg"))
+rfc_gpkg <- file.path("cache", paste0("refactor_", vpu_codes, ".gpkg"))
-gf_gpkg <- file.path("cache", paste0("GF_", vpu, ".gpkg"))
+gf_gpkg <- file.path("cache", paste0("GF_", vpu_codes, ".gpkg"))
gf_gpkg_conus <- "cache/reference_CONUS.gpkg"
gage_info_gpkg <- "cache/Gages_info.gpkg"
@@ -139,7 +146,7 @@ gage_info_csv <- "cache/Gages_info.csv"
gage_info_csvt <- "cache/Gages_info.csvt"
# Defined during NonDend.Rmd
-ND_gpkg <- file.path("temp", paste0("ND_", vpu,".gpkg"))
+ND_gpkg <- file.path("temp", paste0("ND_", vpu_codes,".gpkg"))
divides_xwalk <- "divides_nhd"
HRU_layer <- "all_cats"
diff --git a/workspace/R/utils.R b/workspace/R/utils.R
index d57ff5a6b17f058fbad10f17dccdc7403e3d0c8b..8c0d6da739b84cce76d787a0bdc7d7345856ea63 100644
--- a/workspace/R/utils.R
+++ b/workspace/R/utils.R
@@ -161,6 +161,7 @@ gage_document <- function(data, source, drop_reason){
}
}
+
#' Merges geopackages together to create CONUs geopackage of features
#' @param feat (character) Type of feature to merge (POIs, segments) character
#' @param in_gpkg (character) geopackage to merge (GF, ND_, etc.)
@@ -177,7 +178,8 @@ Merge_VPU <- function(feat, in_gpkg, out_gpkg){
as.character(11:18), "10U", "10L")
featCON <- do.call("rbind", lapply(VPUs, function(x) {
tryCatch({
- layer <- ifelse(feat %in% c("nhd_flowline", "unassigned_gages", "unassigned_TE"),
+ layer <- ifelse(feat %in% c("POIs", "nhd_flowline", "unassigned_gages", "unassigned_TE",
+ "reference_flowline", "reference_catchment"),
feat, paste0(feat, "_", x))
read_sf(file.path("cache", paste0(in_gpkg, x, ".gpkg")),
layer)},
@@ -248,14 +250,15 @@ merge_conus <- function(rpu_vpu, rpu_vpu_out,
vpus <- unique(rpu_vpu$vpuid)
+ #gf_gpkgs <- paste0("cache/refactor_", vpus, ".gpkg")
gf_gpkgs <- paste0("cache/GF_", vpus, ".gpkg")
hydrofab::assign_global_identifiers(
gpkgs = gf_gpkgs,
outfiles = gsub("cache", "temp", gf_gpkgs),
- flowpath_layer = agg_fline_layer,
+ flowpath_layer = agg_fline_layer, #reconciled_layer,
mapped_POI_layer = mapped_outlets_layer,
- divide_layer = agg_cats_layer,
+ divide_layer = agg_cats_layer, #divide_layer,
lookup_table_layer = lookup_table_refactor,
flowpath_edge_list = catchment_network_table,
update_terminals = TRUE,
@@ -264,7 +267,10 @@ merge_conus <- function(rpu_vpu, rpu_vpu_out,
}
-merge_refactor <- function(rpus, rpu_vpu_out,
+
+
+merge_refactor <- function(rpus,
+ rpu_vpu_out,
lookup_table_refactor,
reconciled_layer,
divide_layer,
@@ -280,6 +286,7 @@ merge_refactor <- function(rpus, rpu_vpu_out,
s <- 10000000
for(rpu in rpus) {
+ print(rpu)
refactor_gpkg <- file.path("cache", paste0(rpu, "_refactor.gpkg"))
aggregate_gpkg <- file.path("cache", paste0(rpu, "_aggregate.gpkg"))
@@ -290,14 +297,14 @@ merge_refactor <- function(rpus, rpu_vpu_out,
sf::read_sf(refactor_gpkg, split_divide_layer),
sf::read_sf(aggregate_gpkg, agg_fline_layer),
sf::read_sf(aggregate_gpkg, agg_cats_layer),
- sf::read_sf(aggregate_gpkg, mapped_outlets_layer)),
+ sf::read_sf(aggregate_gpkg, mapped_outlets_layer)),
c(lookup_table_refactor,
reconciled_layer,
divide_layer,
split_divide_layer,
agg_fline_layer,
agg_cats_layer,
- mapped_outlets_layer))
+ mapped_outlets_layer))#,
### RECONCILED ###
@@ -362,7 +369,7 @@ merge_refactor <- function(rpus, rpu_vpu_out,
#### aggregate layers ####
- for(l in c(agg_cats_layer, agg_fline_layer, mapped_outlets_layer)) {
+ for(l in c(agg_cats_layer, agg_fline_layer)) { #, mapped_outlets_layer
out[[rpu]][[l]] <- left_join(rename(out[[rpu]][[l]], old_ID = ID),
update_id,
diff --git a/workspace/cache/data_paths.json b/workspace/cache/data_paths.json
index 766abadf53d6f89ce77f9e57aae2fa2bbb3d1aab..94643cd131df32ca9e4f99df0ac9f3b7475ee5c8 100644
--- a/workspace/cache/data_paths.json
+++ b/workspace/cache/data_paths.json
@@ -10,9 +10,25 @@
"islands_dir": "data/islands",
"islands_gdb": "data/islands/NHDPlusNationalData/NHDPlusV21_National_Seamless_Flattened_HI_PR_VI_PI.gdb",
"nhdplus_rpu": "data/NHDPlusNationalData/NHDPlusGlobalData/BoundaryUnit.shp",
+ "ref_fab_path": "data/reference_fabric",
+ "ref_cat": "data/reference_fabric/reference_catchments.gpkg",
+ "ref_fl": "data/reference_fabric/reference_flowline.gpkg",
+ "nwm_fl": "data/reference_fabric/nwm_network.gpkg",
"waterbodies_path": "data/NHDPlusNationalData/nhdplus_waterbodies.rds",
"fullcats_table": "data/NHDPlusNationalData/nhdcat_full.rds",
"islandcats_table": "data/islands/nhdcat_full.rds",
+ "latest_wbd_rds": "data/wbd/WBD.rds",
+ "merit_catchments": "data/merged_AK_MERIT_Hydro/cat_pfaf_78_81_82_MERIT_Hydro_v07_Basins_v01.shp",
+ "merit_rivers": "data/merged_AK_MERIT_Hydro/riv_pfaf_78_81_82_MERIT_Hydro_v07_Basins_v01.shp",
+ "ak_source": "data/AK/ak.gpkg",
+ "hi_source": "data/islands/hi.gpkg",
+ "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",
+ "istarf": "data/reservoir_data/ISTARF-CONUS.csv",
+ "resops_NID_CW": "data/reservoir_data/cw_ResOpsUS_NID.csv",
"fdr": {
"rpu_18a": "data/fdrfac/NHDPlusCA/NHDPlus18/NHDPlusFdrFac18a/fdr",
"rpu_18b": "data/fdrfac/NHDPlusCA/NHDPlus18/NHDPlusFdrFac18b/fdr",
@@ -237,26 +253,5 @@
"rpu_12b": "data/nhdplusv2_elev/NHDPlusTX/NHDPlus12/NEDSnapshot/NED12b/elev_cm",
"rpu_12c": "data/nhdplusv2_elev/NHDPlusTX/NHDPlus12/NEDSnapshot/NED12c/elev_cm",
"rpu_12d": "data/nhdplusv2_elev/NHDPlusTX/NHDPlus12/NEDSnapshot/NED12d/elev_cm"
- },
- "latest_wbd_rds": "data/wbd/WBD.rds",
- "merit_catchments": "data/merged_AK_MERIT_Hydro/cat_pfaf_78_81_82_MERIT_Hydro_v07_Basins_v01.shp",
- "merit_rivers": "data/merged_AK_MERIT_Hydro/riv_pfaf_78_81_82_MERIT_Hydro_v07_Basins_v01.shp",
- "aster_dem": "data/merged_AK_MERIT_Hydro/dem.tif",
- "merit_dem": "data/merged_AK_MERIT_Hydro/ak_merit_dem.tif",
- "merit_fdr": "data/merged_AK_MERIT_Hydro/ak_merit_fdr.tif",
- "merit_fac": "data/merged_AK_MERIT_Hydro/ak_merit_fac.tif",
- "ak_source": "data/AK/ak.gpkg",
- "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": "data/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",
- "ref_flowline": "data/NHDPlusNationalData/reference_flowline.gpkg",
- "ref_catchment": "data/NHDPlusNationalData/reference_catchments.gpkg",
- "res_attributes": "data/reservoir_data/ResOpsUS/attributes/reservoir_attributes.csv",
- "istarf": "data/reservoir_data/ISTARF-CONUS.csv",
- "resops_NID_CW": "data/reservoir_data/cw_ResOpsUS_NID.csv"
+ }
}
diff --git a/workspace/hyfabric_0.5.7.tar.gz b/workspace/hyfabric_0.5.7.tar.gz
deleted file mode 100644
index 1a5c37d0d10521ed10cfb76e4b42102f1b6fd52e..0000000000000000000000000000000000000000
Binary files a/workspace/hyfabric_0.5.7.tar.gz and /dev/null differ
diff --git a/workspace/hyfabric_0.5.8.tar.gz b/workspace/hyfabric_0.5.8.tar.gz
new file mode 100644
index 0000000000000000000000000000000000000000..60247d74b598541bbe56e0b25f1a7536a4f89948
Binary files /dev/null and b/workspace/hyfabric_0.5.8.tar.gz differ