diff --git a/workspace/00_AK_prep.Rmd b/workspace/00_AK_prep.Rmd
index a6bc276231c4ddf328cc9f62d6278d38f7630e17..0dd4430ef5cbebbd7d31ee5e403d463fd8166c1f 100644
--- a/workspace/00_AK_prep.Rmd
+++ b/workspace/00_AK_prep.Rmd
@@ -46,18 +46,24 @@ if(needs_layer(ak_gpkg, ak_flowline)) {
ak_rivers[["arbolate_sum"]] <- st_drop_geometry(ak_rivers) %>%
select(ID = COMID, toID = NextDownID, length = lengthkm) %>%
- nhdplusTools::calculate_arbolate_sum()
+ hydroloom::accumulate_downstream("length")
ak_rivers[["order"]] <- st_drop_geometry(ak_rivers) %>%
select(ID = COMID, toID = NextDownID) %>%
- get_streamorder()
+ hydroloom::add_streamorder() %>%
+ dplyr::pull(stream_order)
lp <- st_drop_geometry(ak_rivers) %>%
select(ID = COMID, toID = NextDownID, weight = arbolate_sum) %>%
mutate(nameID = "constant") %>%
- get_levelpaths()
+ hydroloom::add_levelpaths(weight_attribute = "weight",
+ name_attribute = "nameID",
+ override_factor = NULL) %>%
+ select(dplyr::all_of(c("ID" = "ID", "outletID" = "levelpath_outlet_id",
+ "topo_sort" = "topo_sort", "levelpath" = "levelpath"))) |>
+ as.data.frame()
- lp <- rename(lp, COMID = .data$ID, outletCOMID = .data$outletID)
+ lp <- rename(lp, COMID = ID, outletCOMID = outletID)
ak_rivers <- left_join(ak_rivers, lp, by = "COMID")
@@ -70,8 +76,13 @@ if(needs_layer(ak_gpkg, ak_flowline)) {
ak_outlets <- ak_rivers$COMID[ak_rivers$toCOMID == 0]
- terminal_paths <- get_terminal(select(st_drop_geometry(ak_rivers), COMID, toCOMID), ak_outlets)
-
+ terminal_paths <- dplyr::left_join(data.frame(ID = ak_rivers$COMID),
+ hydroloom::sort_network(
+ select(st_drop_geometry(ak_rivers), COMID, toCOMID),
+ split = TRUE,
+ outlets = ak_outlets) |>
+ dplyr::select(dplyr::all_of(c("ID" = "COMID", "terminalID" = "terminal_id"))))
+
terminal_paths <- left_join(terminal_paths,
select(st_drop_geometry(ak_rivers),
terminalID = COMID,
@@ -86,15 +97,16 @@ if(needs_layer(ak_gpkg, ak_flowline)) {
ak_rivers$toCOMID[ak_rivers$toCOMID == 0] <- NA
- pathlength <- get_pathlength(select(st_drop_geometry(ak_rivers),
- ID = COMID,
- toID = toCOMID,
- length = LENGTHKM))
+ pathlength <- hydroloom::add_pathlength(select(st_drop_geometry(ak_rivers),
+ ID = COMID,
+ toID = toCOMID,
+ length_km = LENGTHKM)) |>
+ dplyr::select(ID, pathlength = pathlength_km)
ak_rivers <- left_join(ak_rivers,
- select(pathlength,
- COMID = ID,
- Pathlength = pathlength),
+ select(pathlength, dplyr::all_of(c(
+ COMID = "ID",
+ Pathlength = "pathlength"))),
by = "COMID")
# Same as mainstems:::prep_net
diff --git a/workspace/00_get_data.Rmd b/workspace/00_get_data.Rmd
index ab03de56e22722f17bca79b65e0e0101ea6f434d..9c8b496866baf5e911b85c2dbf3a936dfb9e1247 100644
--- a/workspace/00_get_data.Rmd
+++ b/workspace/00_get_data.Rmd
@@ -14,12 +14,14 @@ are made to the output of this notebook, they should be checked in.
**If resources from ScienceBase need to be downloaded Rmarkdown document should be run from RStudio so username and password authentication will work**
```{r}
-source("R/config.R")
source("R/utils.R")
+source("R/config.R")
source("R/user_vars.R")
+source("R/1_get_data.R")
if(!dir.exists("data")) {dir.create("data")}
if(!dir.exists("bin")) {dir.create("bin")}
+
data_dir <- "data"
out_list <- list("data_dir" = data_dir)
out_file <- file.path("cache", "data_paths.json")
@@ -45,20 +47,19 @@ 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.
-hu12_points_path <- file.path(data_dir, "102020wbd_outlets.gpkg")
-
-if(!file.exists(hu12_points_path)) {
- check_auth()
- sbtools::item_file_download("63cb38b2d34e06fef14f40ad",
- names = "102020wbd_outlets.gpkg",
- destinations = hu12_points_path)
-}
+out_list <- c(
+ out_list,
+ list(hu12_points_path =
+ get_sb_file(item = "63cb38b2d34e06fef14f40ad",
+ item_files = "102020wbd_outlets.gpkg",
+ out_destination = data_dir)))
-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
@@ -74,16 +75,13 @@ support of the second dataset of basin characteristics and disturbance indexes.
# flowlines to support Streamgage Watershed InforMation (SWIM), 2021: U.S.
# Geological Survey data release, https://doi.org/10.5066/P9J5CK2Y.
-SWIM_points_path <- file.path(data_dir, "SWIM_gage_loc")
-
-if(!file.exists(SWIM_points_path)) {
- dir.create(SWIM_points_path, recursive = TRUE, showWarnings = FALSE)
- check_auth()
- item <- retrieve_sb("5ebe92af82ce476925e44b8f", data_dir)
-
-}
+out_list <- c(
+ out_list,
+ list(SWIM_points_path =
+ get_sb_file(item = "5ebe92af82ce476925e44b8f",
+ item_files = "all",
+ out_destination = file.path(data_dir, "SWIM_gage_loc"))))
-out_list <- c(out_list, list(SWIM_points_path = SWIM_points_path))
if(mapview)(mapview(read_sf(out_list$SWIM_points_path)))
```
@@ -108,32 +106,18 @@ Galanter and othes, 2023).
TE_points_path <- file.path(data_dir, "TE_points")
-if(!file.exists(TE_points_path)) {
- dir.create(TE_points_path, recursive = TRUE, showWarnings = FALSE)
-
- TEz <- "2015_TE_Model_Estimates_lat.long_COMIDs.7z"
- check_auth()
- sbtools::item_file_download("5dbc53d4e4b06957974eddae", names = TEz,
- destinations = file.path(TE_points_path, TEz))
-
- system(paste0(sevenz, " e -o", TE_points_path, " ",
- file.path(TE_points_path, TEz)))
-
- rm(TEz)
-
- TEz_2023 <- "galanter_and_others_2023.zip"
- check_auth()
- sbtools::item_file_download("63adc826d34e92aad3ca5af4", names = TEz_2023,
- destinations =
- file.path(TE_points_path, TEz_2023))
-
- system(paste0(sevenz, " e -o", TE_points_path, " ",
- file.path(TE_points_path, TEz_2023)))
-
- rm(TEz_2023)
-}
+dir.create(TE_points_path, recursive = TRUE, showWarnings = FALSE)
+
+get_sb_file("5dbc53d4e4b06957974eddae",
+ "2015_TE_Model_Estimates_lat.long_COMIDs.7z",
+ TE_points_path, unzip = TRUE)
+
+get_sb_file("63adc826d34e92aad3ca5af4",
+ "galanter_and_others_2023.zip",
+ TE_points_path, unzip = TRUE)
out_list <- c(out_list, list(TE_points_path = TE_points_path))
+
if(mapview)(mapview(read_sf(out_list$TE_points_path)))
```
@@ -378,10 +362,12 @@ VAA_rds_path <- file.path(nhdplus_dir, "nhd_vaa.rds")
if(!file.exists(VAA_rds_path)) {
message("downloading NHDPlusVAA...")
- get_vaa(path = nhdplus_dir)
+ download_vaa(VAA_fst_path)
+
# Extract other attributes
nhdpv2_vaa <- read_sf(nhdplus_gdb, "NHDFlowline_Network") %>%
st_drop_geometry()
+
saveRDS(nhdpv2_vaa, file.path(nhdplus_dir, "nhd_vaa.rds"))
}
@@ -550,6 +536,15 @@ National Hydrologic Modeling, Alaska Domain
merit_dir <- file.path(data_dir, "merged_AK_MERIT_Hydro")
+out <- list(merit_catchments = file.path(merit_dir, "merged_AK_MERIT_Hydro",
+ "cat_pfaf_78_81_82_MERIT_Hydro_v07_Basins_v01.shp"),
+ merit_rivers = file.path(merit_dir, "merged_AK_MERIT_Hydro",
+ "riv_pfaf_78_81_82_MERIT_Hydro_v07_Basins_v01.shp"),
+ aster_dem = file.path(merit_dir, "dem.tif"),
+ merit_dem = file.path(merit_dir, "ak_merit_dem.tif"),
+ merit_fdr = file.path(merit_dir, "ak_merit_fdr.tif"),
+ merit_fac = file.path(merit_dir, "ak_merit_fac.tif"))
+
if(!dir.exists(merit_dir)) {
dir.create(merit_dir, recursive = TRUE, showWarnings = FALSE)
@@ -587,15 +582,6 @@ 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,
- "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"),
- aster_dem = file.path(merit_dir, "dem.tif"),
- merit_dem = file.path(merit_dir, "ak_merit_dem.tif"),
- merit_fdr = file.path(merit_dir, "ak_merit_fdr.tif"),
- 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,
diff --git a/workspace/01_HI_prep.Rmd b/workspace/01_HI_prep.Rmd
index 7af7a919a6c4861713fdc0bd1289008271e9b8dc..48e4d3b8f5eda272f69237a74f3c1dee1c6d40e1 100644
--- a/workspace/01_HI_prep.Rmd
+++ b/workspace/01_HI_prep.Rmd
@@ -15,10 +15,10 @@ library(dplyr)
library(sf)
library(hyRefactor)
library(readxl)
-library(rgdal)
# Numeric code for HI VPU
vpu_code <- "20"
+vpu <- "20"
# directory of data
data_dir <- "data"
@@ -149,7 +149,7 @@ if(needs_layer(data_paths$hi_source, NID_layer)) {
message("formatting NID locations...")
# Download NID file from website for 2019
- NID_file <- download.file(url = "https://nid.usace.army.mil/api/nation/csv",
+ NID_file <- download.file(url = "https://nid.sec.usace.army.mil/api/nation/csv",
destfile = file.path(data_paths$NID_points_path, "nation.csv"), method = "curl", mode = "wb")
# Convert to simple features
diff --git a/workspace/02_AK_navigate.Rmd b/workspace/02_AK_navigate.Rmd
index 001cddc6e38c6338f08b2a617dc194a54dde197b..4a2808622dbb90c5cf8e6162f77a1c1cb8de3669 100644
--- a/workspace/02_AK_navigate.Rmd
+++ b/workspace/02_AK_navigate.Rmd
@@ -31,7 +31,6 @@ vpu_codes <- "19"
source("R/utils.R")
source("R/config.R")
source("R/NHD_navigate.R")
-source("R/hyRefactor_funs.R")
source("R/ExCONUS_config.R")
# Load user-defined data paths
diff --git a/workspace/02_HI_navigate.Rmd b/workspace/02_HI_navigate.Rmd
index 8e17e7066e36e45e87db049cab71ce2127b6b730..6deba6b39ccfe43daf79909506291b1cc8a6821f 100644
--- a/workspace/02_HI_navigate.Rmd
+++ b/workspace/02_HI_navigate.Rmd
@@ -23,7 +23,7 @@ library(hyRefactor)
library(tidyr)
library(hyfabric)
-vpu <- VPU <- "20"
+vpu_code <- vpu <- VPU <- "20"
# Load data paths
data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
@@ -31,10 +31,12 @@ data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
# Load custom functions
source("R/utils.R")
source("R/NHD_navigate.R")
-source("R/hyRefactor_funs.R")
source("R/config.R")
source("R/ExCONUS_config.R")
+# TODO: Fix up configuration
+max_elev_TT_DA <- 10
+
atts <- readRDS(file.path(data_paths$islands_dir, "nhdplus_flowline_attributes.rds"))
```
@@ -87,9 +89,9 @@ if(needs_layer(nav_gpkg, nhd_flowline)) {
rbind(missing_LP_20g)
# Re-calc drainage area, set ND to false
- nhd_prep <- filter(nhd, COMID %in% prepare_nhdplus(nhd, 0, 0, 0, FALSE, skip_toCOMID = TRUE)$COMID) %>%
+ nhd_prep <- filter(nhd, COMID %in% dplyr::filter(sf::st_drop_geometry(nhd), FCODE != 56600)$COMID) %>%
mutate(TotDASqKM =
- nhdplusTools::calculate_total_drainage_area(select(., ID = COMID, toID = toCOMID, area = AreaSqKM)))
+ hydroloom::accumulate_downstream(select(., ID = COMID, toID = toCOMID, area = AreaSqKM), "area"))
# Coastal and non-dendritic flowlines below the size threshold for poi selection and refactoring
non_dend_COM <- unique(unlist(lapply(filter(nhd, TerminalFl == 1 & TotDASqKM < min_da_km)
@@ -126,6 +128,7 @@ if(needs_layer(nav_gpkg, nhd_flowline)) {
write_sf(nhd_prep, nav_gpkg, nhd_flowline)
} else {
nhd_prep <- read_sf(nav_gpkg, nhd_flowline)
+ nhd <- read_sf(nav_gpkg, nhd_flowline)
}
mapview(filter(nhd_prep, dend == 1))
@@ -234,7 +237,8 @@ if(needs_layer(nav_gpkg, nav_poi_layer)) {
mutate(numCom = n()) %>%
# Ensure specific stations listed above are preserved as POIs,
# otherwise use one with the most records
- filter(if(numCom == 1 | Station %in% other_loc_gages) TRUE else nrecords == max(nrecords, na.rm = T)) %>%
+ filter(if(any(numCom == 1 | Station %in% other_loc_gages)) TRUE else
+ nrecords == max(nrecords, na.rm = TRUE)) %>%
select(COMID, Station)
# Write raw gages data out to use as split events
@@ -380,7 +384,7 @@ if(all(is.na(nav_POIs$Type_WBOut))) {
select(Hydroseq, WBAREACOMI), by = c("DnHydroseq" = "Hydroseq")) %>%
select(COMID, WBAREACOMI) %>%
group_by(COMID) %>%
- slice(n = 1)
+ slice(1)
nav_POIs <- POI_creation(filter(wbout_COMIDs, !COMID %in% wbin_COMIDs$COMID), filter(nhd_prep, dend == 1), "WBOut") %>%
addType(., nav_POIs, "WBOut", bind = FALSE)
@@ -432,18 +436,20 @@ This chunk breaks up potential aggregated segments where the elevation range of
inc_segs <- segment_increment(filter(nhd_prep, minNet == 1), filter(st_drop_geometry(nhd),
COMID %in% nav_POIs$COMID, COMID %in% filter(nhd_prep, minNet == 1)$COMID)) %>%
# bring over VAA data
- inner_join(select(atts, COMID, LENGTHKM, MAXELEVSMO, MINELEVSMO, AreaSqKM, TotDASqKM), by = c("COMID" = "ComID"))
+ inner_join(select(atts, ComID, LengthKM, MaxElevSmo, MinElevSmo, AreaSqKm, TotDASqKm), by = c("COMID" = "ComID"))
# Build dataframe for creation of points based on breaks
elev_pois_split <- inc_segs %>%
group_by(POI_ID) %>%
# Get elevation info
- mutate(MAXELEVSMO = na_if(MAXELEVSMO, -9998), MINELEVSMO = na_if(MINELEVSMO, -9998),
- elev_diff_seg = max(MAXELEVSMO) - min(MINELEVSMO)) %>%
+ mutate(MaxElevSmo = na_if(MaxElevSmo, -9998), MinElevSmo = na_if(MinElevSmo, -9998),
+ elev_diff_seg = max(MaxElevSmo) - min(MinElevSmo)) %>%
+
+ # TODO: I have no ide awhy the pipes below here were commented out!?!
# Filter out to those incremental segments that need splitting above the defined elevation threshold
- filter((max(MAXELEVSMO) - min(MINELEVSMO)) > (elev_diff * 100)) #%>%
+ filter((max(MaxElevSmo) - min(MinElevSmo)) > (elev_diff * 100)) #%>%
# Do not aggregate on fake lfowlines
- filter(AreaSqKM > 0, TotDASqKM > max_elev_TT_DA) #%>%
+ filter(AreaSqKm > 0, TotDASqKm > max_elev_TT_DA) #%>%
# Order from upstream to downstream
arrange(desc(Hydroseq)) #%>%
# Get attributes used for splitting
@@ -451,7 +457,7 @@ elev_pois_split <- inc_segs %>%
# 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),
+ 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) %>%
@@ -460,7 +466,12 @@ elev_pois_split <- inc_segs %>%
# Track original iteration number
mutate(orig_iter = iter) %>%
ungroup()
-
+ # TODO END
+
+ elev_pois_split <- dplyr::rename(elev_pois_split, AreaSqKM = AreaSqKm,
+ TotDASqKM = TotDASqKm)
+
+ # TODO: This is failing and I don't know why!
if(nrow(elev_pois_split) > 0) {
# For reach iteration
diff --git a/workspace/02_NHD_navigate.Rmd b/workspace/02_NHD_navigate.Rmd
index 1df172aba6adebe87732e6cd87c81eeb75638b2e..fd0e904059f3369b9321ad6b5b1702d831a851de 100644
--- a/workspace/02_NHD_navigate.Rmd
+++ b/workspace/02_NHD_navigate.Rmd
@@ -23,7 +23,6 @@ knitr::opts_chunk$set(
# Load custom functions
source("R/utils.R")
source("R/NHD_navigate.R")
-source("R/hyrefactor_funs.R")
#source("R/test.R")
# increase timeout for data downloads
@@ -850,7 +849,7 @@ if(needs_layer(temp_gpkg, pois_all_layer)) {
# bind together
final_POIs <- bind_rows(tmp_POIs_fixed, new_POIs) %>%
mutate(Type_Term = ifelse(nexus == 1, NA, Type_Term)) %>%
- select(-ID)
+ select(-dplyr::any_of(c("ID")))
# Write out fixes
write_sf(new_POIs, temp_gpkg, poi_moved_layer)
diff --git a/workspace/02_NHD_navigate_Sp.Rmd b/workspace/02_NHD_navigate_Sp.Rmd
index daebbeaf4e90d4b3cc80a279b442aabd346ae54c..336bc0f4d81f12b06d4c1166294bcf87b6b433c9 100644
--- a/workspace/02_NHD_navigate_Sp.Rmd
+++ b/workspace/02_NHD_navigate_Sp.Rmd
@@ -21,7 +21,6 @@ knitr::opts_chunk$set(
# Load custom functions
source("R/utils.R")
source("R/NHD_navigate.R")
-source("R/hyrefactor_funs.R")
# increase timeout for data downloads
options(timeout=600)
diff --git a/workspace/03_hyRefactor_AK.Rmd b/workspace/03_hyRefactor_AK.Rmd
index 9e2c113e433b4d0015c67b54750f66a994b92352..65477afa9da0cdf138530d5e462c51c8b7f53998 100644
--- a/workspace/03_hyRefactor_AK.Rmd
+++ b/workspace/03_hyRefactor_AK.Rmd
@@ -28,7 +28,6 @@ VPU <- substr(rpu_code, 1, 2)
out_agg_gpkg <- file.path("cache", paste0(rpu_code, ".gpkg"))
source("R/utils.R")
-source("R/hyRefactor_funs.R")
source("R/config.R")
source("R/ExCONUS_config.R")
@@ -104,7 +103,6 @@ if(needs_layer(out_refac_gpkg, reconciled_layer)){
tf <- paste0("temp/refactored_", rpu_code,".gpkg")
tr <- paste0("temp/reconciled_", rpu_code, ".gpkg")
- # see hyRefactor_funs.R for these parameters.
refactor_nhdplus(nhdplus_flines = nhd,
split_flines_meters = 10000000,
split_flines_cores = 6,
@@ -211,7 +209,7 @@ if(needs_layer(out_agg_gpkg, mapped_outlets_layer)) {
# Join with divides to get the total drainage area
left_join(select(st_drop_geometry(divides), ID, area = areasqkm), by = "ID") %>%
mutate(area = ifelse(is.na(area), 0, area)) %>%
- mutate(TotDASqKM = calculate_total_drainage_area(select(., ID, toID, area)))
+ mutate(TotDASqKM = hydroloom::accumulate_downstream(select(., ID, toID, area), "area"))
# If there are events on split flowlines
if(!needs_layer(out_refac_gpkg, split_layer)){
diff --git a/workspace/03_hyRefactor_HI.Rmd b/workspace/03_hyRefactor_HI.Rmd
index 87d5a95bdd93c4279831c58068584093a0a0e84d..9372684bd55d9d0fd20190ea6ca8ce2f97d16367 100644
--- a/workspace/03_hyRefactor_HI.Rmd
+++ b/workspace/03_hyRefactor_HI.Rmd
@@ -29,7 +29,6 @@ vpu <- "20"
if(!exists("rpu_code")) (rpu_code <- "20a")
source("R/utils.R")
-source("R/hyRefactor_funs.R")
source("R/config.R")
source("R/ExCONUS_config.R")
```
@@ -151,7 +150,6 @@ if(needs_layer(out_refac_gpkg, refactored_layer)) {
tf <- paste0("temp/refactored_", rpu_code,".gpkg")
tr <- paste0("temp/reconciled_", rpu_code, ".gpkg")
- # see hyRefactor_funs.R for these parameters.
# Original Refactor
refactor_nhdplus(nhdplus_flines = dplyr::select(nhdplus_flines, -FTYPE), # Select controls whether prepare_nhdplus is used.
split_flines_meters = split_meters,
diff --git a/workspace/03_hyRefactor_flines.Rmd b/workspace/03_hyRefactor_flines.Rmd
index 0d0f5e17d2a77055f0d59991a04cc8db11f61262..d23d0b3321d7d3728ff9ac1d305b315b085b9f9d 100644
--- a/workspace/03_hyRefactor_flines.Rmd
+++ b/workspace/03_hyRefactor_flines.Rmd
@@ -24,7 +24,6 @@ Load functions and configuration.
See `hyRefactory_config.R` for all constants.
```{r setup}
source("R/utils.R")
-source("R/hyRefactor_funs.R")
source("R/config.R")
source("R/user_vars.R")
```
@@ -96,9 +95,7 @@ if(needs_layer(out_refac_gpkg, refactored_layer)) {
tr <- paste0("temp/reconciled_", rpu_code, ".gpkg")
unlink(list(tf, tr))
- # see hyRefactor_funs.R for these parameters.
- # Select controls whether prepare_nhdplus is used.
- refactor_nhdplus(nhdplus_flines = dplyr::select(nhdplus_flines, -FTYPE),
+ refactor_nhdplus(nhdplus_flines = dplyr::select(nhdplus_flines, -FTYPE), # Select controls whether prepare_nhdplus is used.
split_flines_meters = split_meters,
split_flines_cores = para_split_flines,
collapse_flines_meters = combine_meters,
diff --git a/workspace/04-1_hyRefactor_cats.Rmd b/workspace/04-1_hyRefactor_cats.Rmd
index d705f80e823f937525fa69f58fb21296214146ba..e920d7a8b7af7953e907cf6cf0fcb49bf3d3fcca 100644
--- a/workspace/04-1_hyRefactor_cats.Rmd
+++ b/workspace/04-1_hyRefactor_cats.Rmd
@@ -24,8 +24,7 @@ source("R/utils.R")
source("R/config.R")
source("R/user_vars_ww.R")
source("R/NHD_navigate.R")
-source("R/hyRefactor_funs.R")
-
+source("R/user_vars.R")
```
This reconcile-catchments step splits and combines catchment boundaries according to the process that was completed in refactor flines.
diff --git a/workspace/04-2_aggregate_cats.Rmd b/workspace/04-2_aggregate_cats.Rmd
index b6c95fe8f340e5fa3f711a25444c55f2c93add12..bb49aa1f076d253620ad695f2bfbb06b7441571f 100644
--- a/workspace/04-2_aggregate_cats.Rmd
+++ b/workspace/04-2_aggregate_cats.Rmd
@@ -24,8 +24,7 @@ source("R/utils.R")
source("R/config.R")
source("R/user_vars.R")
source("R/NHD_navigate.R")
-source("R/hyRefactor_funs.R")
-
+source("R/user_vars.R")
```
Derive set of outlets to aggregate catchments for. There are three sets of
@@ -94,8 +93,8 @@ if(needs_layer(out_agg_gpkg, agg_cats_layer)){
extra_outlets <- filter(reconciled_DA, !ID %in% outlets$ID) %>%
pull(ID)
- extra_nets <- nhdplusTools::get_sorted(select(sf::st_drop_geometry(reconciled),
- ID, toID), outlets = extra_outlets)
+ extra_nets <- hydroloom::sort_network(select(sf::st_drop_geometry(reconciled),
+ ID, toID), outlets = extra_outlets)
# remove extra network from contention
reconciled <- filter(reconciled, !ID %in% extra_nets$ID)
diff --git a/workspace/04_merge.Rmd b/workspace/04_merge.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..1ca27838dd81d5c3788b4b52f923895692bb3d40
--- /dev/null
+++ b/workspace/04_merge.Rmd
@@ -0,0 +1,15 @@
+This notebook is used for creation of national scale output.
+
+The outputs are not used later in the workflow.
+
+```{r}
+library(sf)
+library(hyfabric)
+
+source("R/utils.R")
+source("R/config.R")
+
+in_gpkg <- "reference_"
+
+Merge_VPU(final_poi_layer, in_gpkg, gf_gpkg_conus)
+```
diff --git a/workspace/05-2_NonDend.Rmd b/workspace/05-2_NonDend.Rmd
index 7822d10b241f2326777cefa1b7234e572a3ead0d..77f64e42263f83af84c0428a225b3eb6a151bb7a 100644
--- a/workspace/05-2_NonDend.Rmd
+++ b/workspace/05-2_NonDend.Rmd
@@ -30,6 +30,7 @@ source("R/utils.R")
source("R/config.R")
source("R/user_vars.R")
source("R/non_dend.R")
+source("R/user_vars.R")
# Read in full VPU NHD set
if (vpu_codes == "20"){
diff --git a/workspace/R/1_get_data.R b/workspace/R/1_get_data.R
new file mode 100644
index 0000000000000000000000000000000000000000..19dd825c09acd180855122093277737325701211
--- /dev/null
+++ b/workspace/R/1_get_data.R
@@ -0,0 +1,37 @@
+#' get_sb_file
+#' @param item sciencebase item id
+#' @param item_files character vector of file names. "all" will download all files.
+#' @param out_destination character path to save files into
+#' @param unzip logical if TRUE, files ending with .7z will be un7zipped
+#' @return path that data was saved to
+get_sb_file <- function(item, item_files, out_destination, unzip = TRUE) {
+
+ check_auth()
+
+ if(length(item_files) == 1 && item_files == "all") {
+ item_files <- sbtools::item_list_files(item)$fname
+ }
+
+ out_files <- file.path(out_destination, item_files)
+
+ missing <- out_files[!file.exists(out_files)]
+
+ if(length(missing) > 0) {
+ sbtools::item_file_download(item,
+ names = basename(missing),
+ destinations = missing)
+ }
+
+ if(unzip) {
+
+ un7zip_fs <- missing[grepl("7z$", missing)]
+
+ for(f in un7zip_fs) {
+ system(paste0(sevenz, " e -o", out_destination, " ", f))
+ }
+
+ }
+
+ out_destination
+
+}
diff --git a/workspace/R/NHD_navigate.R b/workspace/R/NHD_navigate.R
index 06fa3c76d7d06525a0dfa0a8273288a8b84420f5..7a4d3f9eeba2ce1070e23da6c68f2cff17f9fa45 100644
--- a/workspace/R/NHD_navigate.R
+++ b/workspace/R/NHD_navigate.R
@@ -590,91 +590,6 @@ split_elev <- function(in_POI_ID, split_DF, elev_diff, max_DA){
return(split_elev_DF)
}
-#' Splits a network based on estimated travel time and drainage area
-#' @param in_POI_ID (integer) POI_ID of aggregated flowline that needs to be
-#' split based on travel time
-#' @param split_DF (sf data.frame) flowlines attributed with POI_ID
-#' @param tt_diff (numeric) max travel time change threshold within a segment
-#' @param max_DA (numeric) minimum drainage area for resulting split
-#'
-#' @return (sf data.frame) flowlines with new POI_IDs identified (elev_POI_ID)
-split_tt <- function(in_POI_ID, split_DF, tt_diff, max_DA){
-
- # subset to a given POI_ID
- split_tt_DF <- filter(split_DF, POI_ID == in_POI_ID) %>%
- filter(AreaSqKM > 0) %>%
- mutate(inc_seg_area = c(TotDASqKM[1], (TotDASqKM - lag(TotDASqKM))[-1]))
-
- if(nrow(split_tt_DF) == 0) {
- return(NULL)
- }
-
- first_iter <- unique(split_tt_DF$iter)
- # Iterate through segs that need splitting
- for (i in c(first_iter:max(split_DF$iter, na.rm = T))){
- #print(i)
-
- # first split
- tt_pois_init_iter <- split_tt_DF %>%
- # filter by iteration, comes into play when multiple splits
- # are required per single POI_ID
- filter(iter == i) %>%
- # Recalculate inc elev, length, and area based on the last split
- mutate(csum_length = cumsum(LENGTHKM),
- cumsum_tt = cumsum(FL_tt_hrs),
- sum_area = cumsum(inc_seg_area)) %>%
- # Determine if split actually necessary on this iteration
- filter(sum(FL_tt_hrs) > tt_diff, TotDASqKM > max_DA) %>%
- # This is an important step, identify which row to split on based on incremental
- # elevaton value. Some flowlines exceed the elevation difference threshold
- # The 'split_index' values determiens which row to create a new POI on
- mutate(split_index = ifelse(nrow(filter(., cumsum_tt > (tt_diff))) == nrow(.),
- 1, nrow(filter(., cumsum_tt < (tt_diff))))) %>%
- # Take out if at last iteration to avoid duplicating existng POIs
- filter(!COMID == POI_ID)
-
- if(nrow(tt_pois_init_iter) == 0){
- #print("done")
- # Done iterating on possible splits
- return(split_tt_DF)}
-
- # Get the flowline POI
- tt_pois_init_pre <- tt_pois_init_iter %>%
- filter(row_number() == split_index) %>%
- mutate(new_POI_ID = COMID) %>%
- select(COMID, POI_ID, new_POI_ID, Hydroseq_cut = Hydroseq)
-
- # Remaining rows downstream after split
- tt_pois_init_post <- tt_pois_init_iter %>%
- left_join(select(tt_pois_init_pre, -COMID), by = "POI_ID") %>%
- filter(Hydroseq < Hydroseq_cut) %>%
- ungroup()
-
- # Test for if we need to split again
- if(nrow(tt_pois_init_post) > 0){
- #print("yope")
- # Re-cacl the iter based on results above
- split_tt_DF <- split_tt_DF %>%
- left_join(select(tt_pois_init_pre, -COMID), by = "POI_ID") %>%
- mutate(tt_POI_ID = ifelse((is.na(tt_POI_ID) & Hydroseq >= Hydroseq_cut), new_POI_ID, tt_POI_ID),
- iter = ifelse((!is.na(new_POI_ID) & Hydroseq < Hydroseq_cut), iter + 1, orig_iter)) %>%
- select(-c(new_POI_ID, Hydroseq_cut)) %>%
- ungroup()
-
- } else {
- split_tt_DF <- split_tt_DF %>%
- left_join(select(tt_pois_init_pre, -COMID), by = "POI_ID") %>%
- mutate(tt_POI_ID = ifelse((is.na(tt_POI_ID) & Hydroseq >= Hydroseq_cut), new_POI_ID, tt_POI_ID),
- iter = ifelse((!is.na(new_POI_ID) & Hydroseq < Hydroseq_cut), 0, orig_iter)) %>%
- select(-c(new_POI_ID, Hydroseq_cut)) %>%
- ungroup()
-
- }
- }
- return(split_tt_DF)
-}
-
-
#' Retrieves waterbodies from NHDArea layer or NHD Hi-Res for ResOpsUs
#' locations if absent from NHD waterbody
#' @param nhd (sf data.frame) Fowlines for given VPU
diff --git a/workspace/R/hyRefactor_funs.R b/workspace/R/hyRefactor_funs.R
deleted file mode 100644
index b2451fb0c0fdb8d85c60484b2522b2535d62de2f..0000000000000000000000000000000000000000
--- a/workspace/R/hyRefactor_funs.R
+++ /dev/null
@@ -1,82 +0,0 @@
-#' Makes run subsets for RPUs with multiple large outlets for parallelization
-#' @param terminal_path : (Character) Terminal Path ID
-#' @param nhd : (SF, data frame) nhd flowlines formatted for hyRefactor
-#' @param reconciled : (SF, data frame) reconciled flowlines
-#' @param divides : (SF, data frame) output of divides operation
-#' @param outlets : (SF, data frame) outlet points for divides/reconciled
-#' @events : (SF, data frame) outlet points for events generated and used to split flowlines/catchments
-#' @return List of SF layers specific to each terminal path which are are inputs to aggregate_catchments
-make_run_subsets <- function(terminal_path, nhd, reconciled, divides, outlets) {
-
- requireNamespace("sf", quietly = TRUE)
-
- # Filter nhd by terminal path
- nhd_sub <- dplyr::filter(nhd, TerminalPa == terminal_path)
-
- comid <- lapply(strsplit(reconciled$member_COMID, split = ","),
- function(x) as.integer(x)[1])
-
- reconciled_sub <- dplyr::filter(reconciled, comid %in% nhd_sub$COMID)
-
- divides_sub <- dplyr::filter(divides, ID %in% reconciled_sub$ID)
-
- # if no divides return nothing
- if(nrow(divides_sub) == 0) return(NULL)
-
- divides_sub <- nhdplusTools:::check_valid(divides_sub)
-
- outlets_sub <- dplyr::filter(outlets, ID %in% reconciled_sub$ID)
-
- # return output ready to go into agg_Cats
- list(nhd = nhd_sub, reconciled = reconciled_sub,
- divides = divides_sub, outlets_fin = outlets_sub)
-}
-
-run_agg2 <- function(rs) {
- ##########################################
- # Slightly modified version of Aggregate_catchments that writes features of terminal path run subset where there was a
- # processing error, while not interrupting the aggregation of remaining run subsets
- #
- # Arguments:
- # rs : (list) list of objects needed for aggregated catchments and linked to a
- # single terminal path
- # rs$Objects:
- # terminal_path : (Character) Terminal Path ID
- # nhd : (SF, data frame) nhd flowlines formatted for hyRefactor
- # reconciled : (SF, data frame) reconciled flowlines
- # divides : (SF, data frame) output of divides operation
- # gf_outlets : (SF, data frame) outlet points for divides/reconciled
- #
- # Returns:
- # Aggregated catchments for a given terminal path
- #
- ########################################
-
- print(rs$number)
- message(paste(Sys.time(), "running", rs$number, "with", nrow(rs$reconciled), "rows \n"))
-
- tryCatch({
- # Run aggregate catchments
- hyRefactor::aggregate_catchments(flowpath = rs$reconciled,
- divide = rs$divides,
- outlets = rs$outlets,
- da_thresh = 1,
- only_larger = TRUE)
- },
- error = function(e) {
-
- write(toString(e), file.path("temp", paste0(unique(rs$nhd$RPUID),
- "_", rs$number, ".txt")), append=TRUE)
- # write out run subset where there was an issue for inspection
- message("ERROR \n")
- error_gpkg <- file.path("temp", paste0(unique(rs$nhd$RPUID),
- "_", rs$number, ".gpkg"))
- sf::write_sf(rs$nhd, error_gpkg, "nhd")
- sf::write_sf(rs$reconciled, error_gpkg, "reconciled")
- sf::write_sf(rs$divides, error_gpkg, "divides")
- sf::write_sf(rs$outlets, error_gpkg, "outlets")
- return(NA)
- })
-}
-
-
diff --git a/workspace/R/non_dend.R b/workspace/R/non_dend.R
index c6f002b1af16c44d1b70ec0fc405a99b708cb027..008afdd66dbfd488ced3f8bd4c926fe127517319 100644
--- a/workspace/R/non_dend.R
+++ b/workspace/R/non_dend.R
@@ -425,83 +425,6 @@ coastal_cats <- function(divides_poi, full_nhd, vpu_HUC12){
}
}
-#' Find little terminal/non-dendritic outlets and determine closest valid catchment with HRU ID
-#' @param nhd_full (sf data frame) full NHD flowline dataset
-#' @param nhd_sub (sf data frame) refactored nhd flowlines
-#' @param divides (sf data frame) data frame of reconciled divides
-#
-#' @return (sf data frame) terminal outlets whose upstream contributing divides are not assigned
-#' an aggregated ID
-outlets_close <- function(nhd_full, nhd_sub, divides_poi){
- # Non-coastal terminal paths
- nhd_full <- filter(nhd_full, FTYPE != "Coastline")
- nhd_full <- st_compatibalize(nhd_full, nhd)
-
- # Switched from little terminal paths to all just to make sure we catch all
- little_outlets <- filter(nhd_full, TerminalFl == 1,
- #TotDASqKM <= min_da_km &
- COMID %in% filter(divides_poi, is.na(POI_ID))$member_COMID) %>%
- dplyr::select(COMID, Hydroseq, DnHydroseq, TerminalPa, rpu = RPUID)
-
- # Networks that terminate and are smaller than the threshold, previously defined
- nhd_sub_outlets <- filter(nhd_sub, !DnHydroseq %in% Hydroseq &
- COMID %in% filter(divides_poi, is.na(POI_ID))$member_COMID)%>%
- dplyr::select(COMID, Hydroseq, DnHydroseq, TerminalPa, rpu = RPUID)
-
- # little outlets
- little_outlets_fin <- little_outlets %>%
- rbind(filter(nhd_sub_outlets, !COMID %in% little_outlets$COMID)) %>%
- arrange(COMID) %>%
- mutate(COMID = as.character(COMID))
-
- if(nrow(little_outlets_fin) == 0){
- return(NULL)
- }
-
- # some VPUs may not have these outlets
- if (nrow(little_outlets_fin) == 0){
- missing_outlets <- NA
- } else {
- # convert flowlines to end points
- xy <- get_node(st_cast(little_outlets_fin, "LINESTRING")) %>%
- st_transform(crs)
-
- # All cats except those with same COMID as little_outlets
- nonout_cats <- filter(divides_poi, !member_COMID %in% little_outlets_fin$COMID)
-
- # All divide catchments with aggregated or POI IDs
- poi_cats <- filter(divides_poi, !is.na(POI_ID))
-
- # cast to multlinestring sf feature for proximity analysis
- cats_lines <- nonout_cats %>%
- group_by(row_ID) %>%
- summarize(do_union = FALSE) %>%
- st_cast("MULTILINESTRING") %>%
- arrange(row_ID)
-
- # find and measure distance from the missing outlet to closest catchment
- missing_outlets <- bind_cols(st_drop_geometry(little_outlets_fin) %>% dplyr::select(outlet_COMID = COMID, rpu), xy) %>%
- st_sf() %>%
- # nearest neighboring catchment
- mutate(closest_cat = nonout_cats[st_nearest_feature(., nonout_cats),]$member_COMID,
- closest_cat_HUC12 = nonout_cats[st_nearest_feature(., nonout_cats),]$HUC_12_int) %>%
- # nearest catchment that has a valid aggregated catchment ID value
- mutate(closest_cat_wHRU = poi_cats[st_nearest_feature(., poi_cats),]$member_COMID,
- closest_cat_wHRU_HUC12 = poi_cats[st_nearest_feature(., poi_cats),]$HUC_12_int,
- closest_cat_wHRU_POI_ID = poi_cats[st_nearest_feature(., poi_cats),]$POI_ID) %>%
- # Distance to nearest catchments determined above
- mutate(dist2edge = round(as.numeric(st_distance(., nonout_cats[st_nearest_feature(., nonout_cats),],
- by_element = TRUE)), 2)) %>%
- mutate(dist2hru = round(as.numeric(st_distance(., poi_cats[st_nearest_feature(., poi_cats),],
- by_element = TRUE)), 2)) %>%
- left_join(dplyr::select(st_drop_geometry(poi_cats), member_COMID, HUC_12_int), by = c("outlet_COMID" = "member_COMID")) %>%
- inner_join(dplyr::select(st_drop_geometry(nonout_cats), member_COMID, HUC_12_int, POI_ID, row_ID), by = c("closest_cat" = "member_COMID")) %>%
- rename(outlet_HUC12 = HUC_12_int.x)
- }
-
- return(missing_outlets)
-}
-
#' Determines length of shared perimeters with adjacent catchments
#' @param divides (sf data frame) data frame of reconciled divides
#
@@ -528,140 +451,6 @@ perims <- function(divides_poi){
return(divides_nopoi_int)
}
-#' Determines length of shared perimeters with adjacent catchments
-#' @param inData (list) list of input arguments
-# cats_in[1] : (data frame, sf) valid data frame of NHD catchments
-# nhd[2] : (data frame, sf) data frame of nhd flowlines
-# missing outlets[3] : (data frame) outlets of catchments w/o POI_ID
-# xWalk[4] : (data frame) HUC12-NHD catchment intersection
-#
-#' @return (sf data frame) DF of cats with populated POI_ID
-assignable_cats <- function(inData){
- cats_in <- inData[[1]]
- nhd_full <- inData[[2]]
- term_outlets <- inData[[3]]
- xWalk <- inData[[4]]
-
- cat_bound <- perims(cats_in)
-
- cats_miss <- filter(cats_in, is.na(POI_ID))
-
- term_outlets_POI <- filter(term_outlets, !is.na(POI_ID) | dist2hru < 1)
-
- # First case, if the outlet is terminal, unconnected, but physically within
- # a catchment with an assigned aggregated catchment ID, assign that outlet the POI_ID of cat
- #term_outlets_inHRU <- filter(term_outlets, dist2hru < 1, !is.na(closest_cat_wHRU_POI_ID))
-
- if (nrow(term_outlets_POI) > 0) {
- # Move POI_ID values over for close distance cases
- term_outlets_POI <- term_outlets_POI %>%
- mutate(POI_ID = ifelse(is.na(POI_ID), closest_cat_wHRU_POI_ID, POI_ID))
-
- # for each missing outlt
- for (inCOMID in term_outlets_POI$outlet_COMID){
- current_miss <- filter(term_outlets_POI, outlet_COMID == inCOMID)
-
- # get the upstream flowline/cat contributiong if it exists
- up <- get_UT(nhd_full, inCOMID)
-
- if(length(up) == 0){
- next
- }
-
- print(paste0("sol. 1, terminal outlet in agg. cat, ", inCOMID))
- cats_in <- cats_in %>%
- mutate(POI_ID = ifelse(member_COMID %in% up, current_miss$POI_ID, POI_ID),
- aggStep = ifelse(member_COMID %in% up, "term_into_HRU", aggStep))
-
- }
- term_outlets <- filter(term_outlets, !outlet_COMID %in% term_outlets_POI$outlet_COMID)
- cat_bound <- perims(cats_in)
- }
-
- # Resolve rest of outlets
- if(nrow(term_outlets) > 0){
-
- #term_outlets_HUC12 <- filter(term_outlets, HUC_12_outlet == HUC_12_nrHRU)
- # for each missing outlt
- for (inCOMID in unique(term_outlets$outlet_COMID)){
- print(inCOMID)
- current_miss <- filter(term_outlets, outlet_COMID == inCOMID) %>%
- group_by(closest_cat_wHRU_POI_ID) %>%
- slice(n = 1)
-
- # get the upstream flowline/cat contributiong if it exists
- up <- get_UT(nhd_full, inCOMID)
-
- up <- up[up %in% cats_miss$member_COMID]
-
- if(length(up) == 0){
- next
- }
-
- # majority border-length assigned POI IDs for cats that drain to missing outlet
- nhdsub <- filter(cat_bound, member_COMID %in% up) %>%
- mutate(outlet_POI = current_miss$POI_ID, outlet_HUC12 = current_miss$outlet_HUC12,
- outlet_toHUC12 = current_miss$closest_cat_wHRU_HUC12, outlet_toPOI_ID = current_miss$closest_cat_wHRU_POI_ID) %>%
- # POI of bordering catchments are not NA, and the majority border is the same HUC12
- # as that closest to the terminal outlet
- filter(!is.na(bound_cat_POI_ID), bound_cat_HUC12 == outlet_toHUC12)
-
- # Summarize by each catchment and its shared perimeters
- summary_HUC12 <- nhdsub %>%
- filter(member_COMID == inCOMID, bound_cat_HUC12 == HUC_12_int) %>%
- distinct(member_COMID, bound_cat_POI_ID)
-
- summary_HUC12_multi <- nhdsub %>%
- group_by(member_COMID, HUC_12_int, bound_cat_HUC12, outlet_toHUC12, outlet_toPOI_ID) %>%
- filter(member_COMID == inCOMID, bound_cat_HUC12 == outlet_toHUC12) %>%
- ungroup() %>%
- distinct(member_COMID, outlet_toPOI_ID)
-
- # where the closest HRU has the same HUC12 as the shared boundary of missing cats
- if(nrow(summary_HUC12) == 1){
- print(paste0("sol. 2, terminal outlet near agg. cat ", inCOMID, " ", summary_HUC12$bound_cat_POI_ID))
-
- final_df_single <- filter(st_drop_geometry(cats_in), member_COMID %in% up) %>%
- dplyr::select(member_COMID, bound_cat_POI_ID = POI_ID) %>%
- mutate(bound_cat_POI_ID = summary_HUC12$bound_cat_POI_ID)
-
- cats_in <- cats_in %>%
- left_join(final_df_single, by = "member_COMID") %>%
- mutate(POI_ID = ifelse(!is.na(bound_cat_POI_ID), bound_cat_POI_ID, POI_ID),
- aggStep = ifelse(!is.na(bound_cat_POI_ID), "term_nearHRU", aggStep)) %>%
- dplyr::select(-bound_cat_POI_ID)
-
- # filter out successfully matched outlet
- term_outlets <- filter(term_outlets, !outlet_COMID %in% current_miss$outlet_COMID)
- #term_outlets$POI_ID <- unique(final_df_single$bound_cat_POI_ID)
-
- } else if(nrow(summary_HUC12_multi) == 1) {
-
- print(paste0("solution 3, ", inCOMID, " ", summary_HUC12_multi$outlet_toPOI_ID))
-
- final_df_double <- filter(st_drop_geometry(cats_in), member_COMID %in% up) %>%
- dplyr::select(member_COMID, bound_cat_POI_ID = POI_ID) %>%
- mutate(bound_cat_POI_ID = unique(summary_HUC12_multi$outlet_toPOI_ID))
-
- cats_in <- cats_in %>%
- left_join(final_df_double, by = "member_COMID") %>%
- mutate(POI_ID = ifelse(!is.na(bound_cat_POI_ID), bound_cat_POI_ID, POI_ID),
- aggStep = ifelse(!is.na(bound_cat_POI_ID), "term_nearmultiHRU", aggStep)) %>%
- dplyr::select(-bound_cat_POI_ID)
-
- term_outlets <- filter(term_outlets, !outlet_COMID %in% current_miss$outlet_COMID)
- #term_outlets$POI_ID <- unique(final_df_double$bound_cat_POI_ID)
- } else {
- print ("no way to aggregate")
- }
- }
- }
-
- cats_in <- dissolve_holes(cats_in)
-
- return(list(divides_lu = cats_in, term_outlets = term_outlets))
-}
-
#' Find and lump remaining catchment sinks where possible
#' @param divides_lu (sf data frame) full set of divides and catchments with no assigned POI_ID
#' @param area_thresh (numeric) size threshold for considering sink as its own aggregated unit
@@ -981,43 +770,3 @@ get_HUC12_Xwalk <- function(vpu, full_cats = FALSE, divides = FALSE, wbd = FALSE
return(list(cats_HUC12 = cats_HUC12, divides_HUC12 = divides_HUC12))
}
-
-# Think this is superseeded by Mike's new code
-# final_assignment <- function (divides_lu){
-# cat_bound <- perims(divides_lu)
-#
-# HUC12_assignment <- filter(cat_bound, HUC_12_int == bound_cat_HUC12) %>%
-# distinct(member_COMID, bound_cat_POI_ID) %>%
-# group_by(member_COMID) %>%
-# filter(n() == 1) %>%
-# ungroup()
-#
-# # Longest shared perimeter with same HUC12
-# cat_bound2 <- filter(cat_bound, !member_COMID %in% HUC12_assignment$member_COMID) %>%
-# filter(!member_COMID %in% nhd$COMID) %>%
-# group_by(member_COMID) %>%
-# filter(HUC_12_int == bound_cat_HUC12) %>%
-# filter(perim == max(perim)) %>%
-# ungroup() %>%
-# distinct(member_COMID, bound_cat_POI_ID) %>%
-# rbind(HUC12_assignment)
-#
-# # For coastal unaggregated catchments that may be outside HUC12
-# cat_bound3 <- filter(cat_bound, !member_COMID %in% cat_bound2$member_COMID) %>%
-# filter(!member_COMID %in% nhd$COMID, is.na(HUC_12_int)) %>%
-# group_by(member_COMID) %>%
-# filter(perim == max(perim)) %>%
-# ungroup() %>%
-# distinct(member_COMID, bound_cat_POI_ID) %>%
-# rbind(cat_bound2)
-#
-# cats_in <- divides_lu %>%
-# left_join(cat_bound3, by = "member_COMID") %>%
-# mutate(POI_ID = ifelse(!is.na(bound_cat_POI_ID), bound_cat_POI_ID, POI_ID),
-# aggStep = ifelse(!is.na(bound_cat_POI_ID), "leftovers", aggStep)) %>%
-# select(-bound_cat_POI_ID)
-#
-# cats_in <- dissolve_holes(cats_in)
-#
-# return(cats_in)
-# }
diff --git a/workspace/R/utils.R b/workspace/R/utils.R
index 2b90f6aba056189adc0f7006bd3c4d16e7d72f8e..7e495b2f0ae3c464ef06f750ed105876da8c7b3b 100644
--- a/workspace/R/utils.R
+++ b/workspace/R/utils.R
@@ -10,40 +10,15 @@ 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
-#' @param out_dir (directory) directory where files are downloaded
-#'
-# Retrieve files from SB if files not listed under sbtools
-retrieve_sb <- function(sbitem, out_dir){
-
- item <- item_get(sbitem)
-
- linked_files <-
- lapply(item$facets, function(x)
- {
- list(name = x$name,
- files = lapply(x$files, function(y, out_path)
- {
-
- out_file <- file.path(out_path, y$name)
- download.file(y$downloadUri, out_file, mode = "wb")
-
- list(fname = y$name,
- url = y$downloadUri)
-
- },
- out_path = file.path(data_dir, x$name)))
- })
+ if(!sbtools::is_logged_in()) {
+ try(sbtools::initialize_sciencebase_session())
+ if(!sbtools::is_logged_in()) stop("stale sciencebase login")
+ }
+ TRUE
}
@@ -436,72 +411,6 @@ nexus_from_poi = function(mapped_POIs,
nexus_locations
}
-#' Generate Lookup table for Aggregated Nextwork
-#' @param gpkg to the aggregated network. The look up table will be added here as a layer called "lookup_table"
-#' @param flowpath_name the name of the flowpath layer in gpkg
-#' @param refactored_fabric the gpkg for the corresponding VPU reference fabric
-#' @return data.frame
-#' @export
-#' @importFrom sf read_sf st_drop_geometry write_sf
-#' @importFrom dplyr mutate select full_join left_join
-#' @importFrom tidyr unnest
-
-generate_lookup_table = function(nl,
- mapped_POIs) {
-
- # THIS IS THE CODE ABOVE!!!!
- nexus_locations = nexus_from_poi(mapped_POIs, verbose = FALSE)
-
- lu = nl$lookup_table %>%
- select(reconciled_ID, NHDPlusV2_COMID)
-
- names(nl$flowpaths) <- tolower(names(nl$flowpaths))
-
- nl$flowpaths <- rename(nl$flowpaths, member_COMID = set)
-
- lu2 = nl$flowpaths %>%
- st_drop_geometry() %>%
- select(
- aggregated_ID = id,
- toID = toid,
- member_COMID = member_comid,
- divide_ID = id,
- POI_ID = poi_id,
- mainstem_id = levelpathid
- ) %>%
- mutate(member_COMID = strsplit(member_COMID, ","),
- POI_ID = as.integer(POI_ID)) %>%
- unnest(col = 'member_COMID') %>%
- mutate(NHDPlusV2_COMID_part = sapply( strsplit(member_COMID, "[.]"),
- FUN = function(x){ x[2] }),
- NHDPlusV2_COMID_part = ifelse(is.na(NHDPlusV2_COMID_part), 1L,
- as.integer(NHDPlusV2_COMID_part)),
- NHDPlusV2_COMID = sapply( strsplit(member_COMID, "[.]"),
- FUN = function(x){ as.numeric(x[1]) })
- ) %>%
- full_join(lu,
- by = "NHDPlusV2_COMID")
-
-
- group_by(aggregated_ID) %>%
- arrange(hydroseq) %>%
- mutate(POI_ID = c(POI_ID[1], rep(NA, n()-1))) %>%
- ungroup() %>%
- select(-hydroseq, - member_COMID) %>%
- left_join(select(nexus_locations, POI_ID = poi_id, value, type),
- by = "POI_ID") %>%
- select(NHDPlusV2_COMID, NHDPlusV2_COMID_part,
- reconciled_ID,
- aggregated_ID, toID, mainstem = mainstem_id,
- POI_ID, POI_TYPE = type, POI_VALUE = value)
-
-
- write_sf(lu2, gpkg, "lookup_table")
-
- return(gpkg)
-
-}
-
# deprecated function from nhdplusTools -- should remove from package but here for compatibility
stage_national_data <- function(include = c("attribute",
"flowline",
@@ -536,7 +445,7 @@ stage_national_data <- function(include = c("attribute",
if (!(file.exists(out_path_flines) | file.exists(out_path_attributes))) {
fline <- sf::st_zm(sf::read_sf(nhdplus_data,
- get_flowline_layer_name(nhdplus_data)))
+ nhdplusTools:::get_flowline_layer_name(nhdplus_data)))
}
if ("attribute" %in% include) {
@@ -566,7 +475,7 @@ stage_national_data <- function(include = c("attribute",
warning("catchment already exists.")
} else {
- layer_name <- get_catchment_layer_name(simplified, nhdplus_data)
+ layer_name <- nhdplusTools:::get_catchment_layer_name(simplified, nhdplus_data)
saveRDS(sf::st_zm(sf::read_sf(nhdplus_data, layer_name)),
out_path_catchments)
diff --git a/workspace/cache/data_paths.json b/workspace/cache/data_paths.json
index 058366f9008126e96bdf0a5593e5369c458165d0..515f3d201042ce7fc440cddf0fd7f51827535a3f 100644
--- a/workspace/cache/data_paths.json
+++ b/workspace/cache/data_paths.json
@@ -245,8 +245,8 @@
"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",
+ "merit_catchments": "data/merged_AK_MERIT_Hydro/merged_AK_MERIT_Hydro/cat_pfaf_78_81_82_MERIT_Hydro_v07_Basins_v01.shp",
+ "merit_rivers": "data/merged_AK_MERIT_Hydro/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",