From 8fdb231187793e2f730d4f93d06deb09962cedad Mon Sep 17 00:00:00 2001
From: David Blodgett <dblodgett@usgs.gov>
Date: Thu, 19 Dec 2024 11:51:34 -0600
Subject: [PATCH] remove unused files
---
workspace/02_NHD_navigate.Rmd | 1070 ----------------------
workspace/02_NHD_navigate_Sp.Rmd | 822 -----------------
workspace/03_hyRefactor_flines.Rmd | 198 ----
workspace/04-1_hyRefactor_cats.Rmd | 76 --
workspace/04-2_aggregate_cats.Rmd | 222 -----
workspace/04_merge.Rmd | 15 -
workspace/05-1_merge.Rmd | 94 --
workspace/05-2_NonDend.Rmd | 412 ---------
workspace/runners/README.md | 19 -
workspace/runners/batch.sh | 13 -
workspace/runners/cleaner.R | 70 --
workspace/runners/local_batch.sh | 9 -
workspace/runners/make_docker_calls.Rmd | 191 ----
workspace/runners/make_shifter_calls.Rmd | 63 --
workspace/runners/run_all_R.R | 16 -
workspace/runners/run_one_R.R | 74 --
16 files changed, 3364 deletions(-)
delete mode 100644 workspace/02_NHD_navigate.Rmd
delete mode 100644 workspace/02_NHD_navigate_Sp.Rmd
delete mode 100644 workspace/03_hyRefactor_flines.Rmd
delete mode 100644 workspace/04-1_hyRefactor_cats.Rmd
delete mode 100644 workspace/04-2_aggregate_cats.Rmd
delete mode 100644 workspace/04_merge.Rmd
delete mode 100644 workspace/05-1_merge.Rmd
delete mode 100644 workspace/05-2_NonDend.Rmd
delete mode 100644 workspace/runners/README.md
delete mode 100644 workspace/runners/batch.sh
delete mode 100644 workspace/runners/cleaner.R
delete mode 100644 workspace/runners/local_batch.sh
delete mode 100644 workspace/runners/make_docker_calls.Rmd
delete mode 100644 workspace/runners/make_shifter_calls.Rmd
delete mode 100644 workspace/runners/run_all_R.R
delete mode 100644 workspace/runners/run_one_R.R
diff --git a/workspace/02_NHD_navigate.Rmd b/workspace/02_NHD_navigate.Rmd
deleted file mode 100644
index 0fb348f..0000000
--- a/workspace/02_NHD_navigate.Rmd
+++ /dev/null
@@ -1,1070 +0,0 @@
----
-title: "MAPNAT POI Geneation"
-output: html_document
-editor_options:
- chunk_output_type: console
----
-This notebook generates a network of Points of Interest harvested from Public
-source that are used to generate a hydrofabric for use within the USGS
-MAPNAT Project.
-
-
-
-```{r setup_rmd, echo=FALSE, cache=FALSE}
-knitr::opts_chunk$set(
- collapse = TRUE,
- comment = "#>",
- fig.width=6,
- fig.height=4,
- cache=FALSE)
-```
-
-```{r setup}
-# Load custom functions
-source("R/utils.R")
-source("R/NHD_navigate.R")
-#source("R/test.R")
-
-# increase timeout for data downloads
-options(timeout=600)
-
-# Load Configuration of environment
-source("R/config.R")
-source("R/user_vars.R")
-
-# Gages output from Gage_selection
-gages <- read_sf(gage_info_gpkg, "Gages")
-
-# need some extra attributes for a few POI analyses
-#atts <- read_fst(data_paths$VAA_fst)
-atts_rds <- readRDS(data_paths$VAA_rds)
-```
-
-```{r huc12 POIs}
-# Derive or load HUC12 POIs
-if(needs_layer(temp_gpkg, nav_poi_layer)) {
-
- nhd <- read_sf(nav_gpkg, nhd_flowline)
- if("POI_ID" %in% colnames(nhd)){
- try(nhd <- select(nhd, -c(wb_id, POI_ID)), silent = TRUE) %>%
- distinct()
- }
-
- # Some NHDPlus VPUs include HUCs from other VPUs
- if(vpu_codes == "02"){
- grep_exp <-"^02|^04"
- } else if (vpu_codes == "08") {
- grep_exp <- "^03|^08"
- } else {
- grep_exp <- paste0("^", substr(vpu_codes, start = 1, stop = 2))
- }
-
- # Join HUC12 outlets with NHD
- HUC12_COMIDs <- read_sf(data_paths$hu12_points_path, "hu_points") %>%
- filter(grepl(grep_exp, .data$HUC12)) %>%
- select(COMID, HUC12) %>%
- # Remove this when HUC12 outlets finished
- group_by(COMID) %>%
- filter(row_number() == 1) %>%
- ungroup()
-
- # Create POIs - some r05 HUC12 POIs not in R05 NHD
- huc12_POIs <- POI_creation(st_drop_geometry(HUC12_COMIDs), filter(nhd, poi == 1), "HUC12")
-
- # Write out geopackage layer representing POIs for given theme
- write_sf(huc12_POIs, temp_gpkg, nav_poi_layer)
- tmp_POIs <- huc12_POIs
-} else {
- # Load HUC12 POIs as the tmpPOIs if they already exist
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
- nhd <- read_sf(nav_gpkg, nhd_flowline)
-}
-
-mapview(filter(tmp_POIs, Type_HUC12 != ""), layer.name = "HUC12 POIs", col.regions = "red")
-```
-
-```{r streamgage POIs}
-if(!"Type_Gages" %in% names(tmp_POIs)) {
-
- # Previously identified streamgages within Gage_Selection.Rmd
- streamgages_VPU <- gages %>%
- rename(COMID = comid) %>%
- filter(COMID %in% nhd$COMID) %>%
- inner_join(st_drop_geometry(st_drop_geometry(nhd), COMID, Divergence), by = "COMID") %>%
- switchDiv(., nhd) %>%
- filter(drain_area > min_da_km_gages)
-
- streamgages <- streamgages_VPU %>%
- group_by(COMID) %>%
- # If multiple gages per COMID, pick one with highest rank as rep POI_ID
- filter(gage_score == max(gage_score), !is.na(drain_area)) %>%
- ungroup()
-
- # 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, 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) {
- write_sf(filter(streamgages_VPU, !site_no %in% tmp_POIs$Type_Gages) %>%
- rename(nhd_reachcode = REACHCODE),
- temp_gpkg, "unassigned_gages")
-
- # Start documenting gages that are dropped out; these gages have no mean daily Q
- document_dropped_gage(filter(streamgages_VPU, !site_no %in% tmp_POIs$Type_Gages) %>%
- select(site_no),
- "Gages_info", paste0("VPU ", vpu_codes, "; low gage score"))
-
- }
-
- # Write out events and outelts
- 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 {
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
- events <- read_sf(temp_gpkg, split_layer)
-}
-
-#mapview(filter(tmp_POIs, !is.na(Type_Gages)), layer.name = "Streamgage POIs", col.regions = "blue")
-```
-
-```{r TE POIs}
-if(!"Type_TE" %in% names(tmp_POIs)) {
-
- if(vpu_codes == "08"){
- nhd$VPUID <- "08"
- } else {
- nhd$VPUID <- substr(nhd$RPUID, 1, 2)
- }
-
- # Read in Thermoelectric shapefile for 2015 estimates
- TE_v1 <- read_sf(data_paths$TE_points_path,
- "2015_TE_Model_Estimates_lat.long_COMIDs")
- # Read in Thermoelectric Plants csv for updated estimates
- TE_v2 <- read.csv(file.path(data_paths$TE_points_path,
- "te_plants.csv"),
- colClasses = "character") %>%
- mutate(Plant.Code = as.integer(Plant.Code))
-
- # Bind hydrographic address information together and get final COMIDs
- TE_COMIDs_full <- TE_v2 %>%
- filter(water_source %in% c("sw_fresh", "gw_fresh", "gw_sw_mix")) %>%
- left_join(select(st_drop_geometry(TE_v1),
- Plant.Code = EIA_PLANT_, COMID), by = "Plant.Code") %>%
- left_join(rename(st_drop_geometry(HUC12_COMIDs), huc12_comid = COMID),
- by = c("huc_12" = "HUC12")) %>%
- mutate(COMID = ifelse(is.na(COMID), huc12_comid, COMID)) %>%
- select(-huc12_comid)
-
- # Prepare TE for POI Creation
- TE_COMIDs <- TE_COMIDs_full %>%
- inner_join(select(st_drop_geometry(nhd), COMID, VPUID), by = "COMID") %>%
- filter(grepl(paste0("^", substr(vpu_codes, 1, 2), ".*"), .data$VPUID), COMID > 0) %>%
- switchDiv(., nhd) %>%
- group_by(COMID) %>%
- summarize(eia_id = paste0(unique(Plant.Code), collapse = " "), count = n()) %>%
- ungroup()
-
- if(nrow(TE_COMIDs) > 0){
- # Derive TE POIs
- tmp_POIs <- POI_creation(st_drop_geometry(TE_COMIDs), filter(nhd, poi == 1), "TE") %>%
- addType(., tmp_POIs, "TE", nexus = TRUE)
-
- # As a fail-safe, write out list of TE plants not assigned a POI
- if(nrow(filter(TE_COMIDs, !COMID %in% tmp_POIs$COMID)) > 0) {
- write_sf(filter(TE_COMIDs, !COMID %in% tmp_POIs$COMID),
- temp_gpkg, "unassigned_TE")
-
- # Write out geopackage layer representing POIs for given theme
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
- }
- } else {
- print ("no TEs in VPU")
- }
-
-} else {
- # Load TE POIs if they already exist
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-}
-
-#mapview(filter(tmp_POIs, Type_TE != ""), layer.name = "TE Plant POIs", col.regions = "blue")
-```
-
-```{r WB mods}
-
-# Waterbodies sourced from NHD waterbody layer for entire VPU
-WBs_VPU_all <- filter(readRDS(data_paths$waterbodies_path),
- COMID %in% nhd$WBAREACOMI) %>%
- filter(FTYPE != "SwampMarsh") %>%
- mutate(FTYPE = as.character(FTYPE),
- source = "NHDv2WB",
- wb_id = ifelse(GNIS_ID == " ", COMID, GNIS_ID)) %>%
- st_transform(geo_crs) %>%
- st_make_valid()
-
-write_sf(WBs_VPU_all, temp_gpkg, "wbs_draft")
-
-if(needs_layer(temp_gpkg, "WBs_layer_orig")){
- sf_use_s2(FALSE)
- ref_WB <- WBs_VPU_all %>%
- group_by(wb_id, GNIS_NAME) %>%
- sfheaders::sf_remove_holes() %>%
- summarize(do_union = TRUE, member_comid = paste(COMID, collapse = ",")) %>%
- st_make_valid() %>%
- ungroup() %>%
- mutate(area_sqkm = as.numeric(st_area(.)) / 1000000) %>%
- mutate(source = "NHDv2WB") %>%
- st_cast("MULTIPOLYGON")
- #sf_use_s2(TRUE)
-
- write_sf(ref_WB, temp_gpkg, "WBs_layer_orig")
-} else {
- ref_WB <- read_sf(temp_gpkg, "WBs_layer_orig")
-}
-```
-
-```{r resops}
-# # 1: Many based on original GNIS_ID
-if(!"Type_resops" %in% names(tmp_POIs)){
-
- # Unnest wb_poi_lst
- wb_table <- st_drop_geometry(ref_WB) %>%
- dplyr::mutate(member_comid = strsplit(member_comid, ",")) %>%
- tidyr::unnest(cols = member_comid) %>%
- mutate(member_comid = as.character(member_comid)) %>%
- distinct()
-
- # TODO: get this file from data.json
- # ResOpsUS locations with attributes from the VPU NHDv2 wb set
- resops_wb_df <- read.csv("data/reservoir_data/istarf_xwalk_final_48_gfv11.csv") %>%
- # Subset to VPU, only one DAMID per waterbody
- filter(flowlcomid %in% nhd$COMID |
- wbareacomi %in% wb_table$member_comid) %>%
- dplyr::select(grand_id, nid_source_featureid, source = comi_srclyr,
- flowlcomid, wbareacomi, hr_permid, onoffnet) %>%
- mutate(wbareacomi = as.character(wbareacomi)) %>%
- # Link to waterbody table
- left_join(distinct(wb_table, member_comid, wb_id),
- by = c("wbareacomi" = "member_comid")) %>%
- left_join(select(st_drop_geometry(WBs_VPU_all), wbareacomi = COMID,
- wb_id2 = wb_id) %>%
- mutate(wbareacomi = as.character(wbareacomi)),
- by = "wbareacomi") %>%
- mutate(wb_id = ifelse(is.na(wb_id), wb_id2, wb_id),
- wbareacomi = ifelse(wbareacomi == -2,
- hr_permid, wbareacomi)) %>%
- select(-c(wb_id2))
-
- # POIs with waterbody IDs in reference waterbodies
- resops_wb_pois <- filter(ref_WB, wb_id %in% resops_wb_df$wb_id) %>%
- inner_join(select(st_drop_geometry(resops_wb_df), grand_id,
- NID_ID = nid_source_featureid,
- resops_flowlcomid = flowlcomid, wb_id),
- by = "wb_id") %>%
- #mutate(source = "ref_WB") %>%
- group_by(wb_id) %>%
- filter(n() == 1) %>%
- st_as_sf()
-
- # Add ResOPsUS locations to waterbody list with attributed NID and resops data
- if(nrow(resops_wb_pois) > 0){
- wb_poi_lst_filtered <- filter(ref_WB, !wb_id %in% resops_wb_pois$wb_id,
- !is.na(wb_id))
-
- wb_poi_lst <- data.table::rbindlist(list(wb_poi_lst_filtered, resops_wb_pois),
- fill = TRUE) %>%
- mutate(accounted = 0) %>%
- st_as_sf()
-
- }
-
- # Reach resopsus
- reach_resops <- filter(resops_wb_df, !wb_id %in% wb_poi_lst$wb_id,
- !source %in% c("NHDAREA", "HR ID AVAILABLE") |
- onoffnet == 0) %>%
- mutate(Type_WBOut = NA) %>%
- select(COMID = flowlcomid, resops = grand_id, Type_WBOut)
-
- # Existing reservoir-waterbody outlet POIs
- exist_POIs_WB <- filter(resops_wb_df, flowlcomid %in% tmp_POIs$COMID) %>%
- mutate(wb_id = ifelse(source == "NHDAREA", wbareacomi, wb_id)) %>%
- filter(!is.na(wb_id)) %>%
- select(COMID = flowlcomid, resops = grand_id, Type_WBOut = wb_id)
-
- # Resops POIs
- resops_pois <- rbind(reach_resops, exist_POIs_WB)
-
- # Resops defined by reach
- if(nrow(resops_pois) > 0){
- # Resops POIs with no reservoirs, defined by reach
- tmp_POIs <- POI_creation(resops_pois, filter(nhd, poi == 1), "resops") %>%
- addType(., tmp_POIs, "resops", nexus = TRUE)
-
- # Add waterbody attribute
- tmp_POIs <- tmp_POIs %>%
- left_join(select(resops_pois, COMID, Type_WBOut),
- by = "COMID") %>%
- mutate(Type_WBOut = ifelse(!nexus,
- Type_WBOut, NA))
-
- # Resops with reservoirs
- resops_wb_df <- resops_wb_df %>%
- mutate(accounted = ifelse(grand_id %in% tmp_POIs$Type_resops, 1, 0),
- source = ifelse(grand_id %in% reach_resops$resops, "REACH", source))
- }
-
- write_sf(wb_poi_lst, temp_gpkg, WBs_layer)
- write_sf(resops_wb_df, temp_gpkg, "wb_resops")
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-} else {
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
- wb_poi_lst <- read_sf(temp_gpkg, WBs_layer)
- resops_wb_df <- read_sf(temp_gpkg, "wb_resops")
-}
-```
-
-```{r Hilari}
-if(!"Type_hilarri" %in% names(tmp_POIs)){
-
- # 1: Many based on original GNIS_ID
- wb_table <- st_drop_geometry(wb_poi_lst) %>%
- dplyr::mutate(member_comid = strsplit(member_comid, ",")) %>%
- tidyr::unnest(cols = member_comid) %>%
- mutate(member_comid = as.integer(member_comid)) %>%
- distinct()
-
- # Hilarri POI Points
- hilarri_points <- read.csv("data/HILARRI/HILARRI_v2.csv",
- colClasses = "character") %>%
- mutate(nhdwbcomid = as.integer(nhdwbcomid)) %>%
- filter(!dataset %in% c('Power plant only; no reservoir or inventoried dam'),
- nhdv2comid %in% nhd$COMID) %>%
- left_join(select(wb_table, member_comid, wb_id),
- by = c("nhdwbcomid" = "member_comid"))
-
- # Waterbodies linked to hilarri information
- hil_wb_pois_rf <- wb_poi_lst %>%
- inner_join(select(st_drop_geometry(hilarri_points), hilarriid, nhdwbcomid,
- nidid_hill = nidid, wb_id),
- by = "wb_id") %>%
- group_by(wb_id) %>%
- mutate(hilarriid = paste0(hilarriid, collapse = ","),
- nidid_hill = paste0(nidid_hill, collapse = ",")) %>%
- ungroup() %>%
- distinct() %>%
- st_as_sf() %>%
- st_compatibalize(wb_poi_lst)
-
- # Add ResOPsUS locations to waterbody list
- if(nrow(hil_wb_pois_rf) > 0){
- wb_poi_lst <- wb_poi_lst %>%
- #filter(!is.na(resops_FlowLcomid)) %>%
- left_join(select(st_drop_geometry(hil_wb_pois_rf), wb_id, hilarriid,
- nidid_hill), by = "wb_id") %>%
- mutate(NID_ID = ifelse(is.na(NID_ID), nidid_hill, NID_ID)) %>%
- select(-nidid_hill) %>%
- distinct()
- }
-
- # Reach POIs
- reach_pois <- filter(hilarri_points, !hilarriid %in% wb_poi_lst$hilarriid) %>%
- select(COMID = nhdv2comid, hilarriid) %>%
- mutate(COMID = as.integer(COMID)) %>%
- group_by(COMID) %>%
- filter(row_number() == 1) %>%
- st_drop_geometry()
-
- # Make POIs for reach POIs
- tmp_POIs <- POI_creation(reach_pois, filter(nhd, poi == 1), "hilarri") %>%
- addType(., tmp_POIs, "hilarri", nexus = TRUE)
-
- # Write out geopackage layer representing POIs for given theme
- write_sf(wb_poi_lst, temp_gpkg, WBs_layer)
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-} else {
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
- wb_poi_lst <- read_sf(temp_gpkg, WBs_layer)
-}
-```
-
-```{r waterbody outlet POIs}
-# Derive or load Waterbody POIs ----------------------
-if("Type_WBOut" %in% names(tmp_POIs)) {
-
- # Hillari IDs
- hilarri_pois <- filter(tmp_POIs, !is.na(Type_hilarri))
-
- # Get NHDARea polygons already define
- nhd_area_vpu <- read_sf(data_paths$nhdplus_gdb, "NHDArea") %>%
- select(COMID, AREASQKM) %>%
- mutate(source = "NHDv2Area") %>%
- st_compatibalize(wb_poi_lst)
-
- # Non-nhdv2 waterbodies from resopsus
- other_wbs <- filter(resops_wb_df, !wb_id %in% wb_poi_lst$wb_id,
- source != "REACH") %>%
- select(grand_id, member_comid = wbareacomi, resops_flowlcomid = flowlcomid,
- source, accounted) %>%
- left_join(select(st_drop_geometry(hilarri_pois),
- COMID, hilarriid = Type_hilarri),
- by = c("resops_flowlcomid" = "COMID")) %>%
- group_by(resops_flowlcomid) %>%
- filter(row_number() == 1) %>%
- ungroup() %>%
- left_join(select(nhd_area_vpu, COMID) %>%
- mutate(COMID = as.character(COMID)),
- by = c("member_comid" = "COMID")) %>%
- st_as_sf() %>%
- mutate(wb_id = ifelse(source == "NHDAREA", member_comid, NA)) %>%
- st_drop_geometry()
-
- # copy empty geometries to wb_list for use in wb POI function
- if(nrow(other_wbs) > 0){
- wb_poi_list_fn <- data.table::rbindlist(list(wb_poi_lst,
- other_wbs),
- fill = TRUE) %>%
- distinct() %>%
- st_as_sf()
- } else {
- wb_poi_list_fn <- wb_poi_lst
- }
-
- # Final waterbody list with all requirements
- final_wb_list <- filter(wb_poi_list_fn,
- area_sqkm > wb_area_thresh |
- !is.na(resops_flowlcomid) |
- !is.na(hilarriid))
-
- new_path <- "?prefix=StagedProducts/Hydrography/NHDPlusHR/VPU/Current/GDB/"
- assign("nhdhr_file_list", new_path, envir = nhdplusTools:::nhdplusTools_env)
-
- # fix something here
- wb_layers <- wbout_POI_creaton(filter(nhd, poi == 1), final_wb_list, data_paths, tmp_POIs, proj_crs)
-
- if(!is.data.frame(wb_layers)){
- # write_sf(wb_layers$nhd_wb_table, temp_gpkg, "wb_flowpaths")
- #
- wbs <- wb_layers$WBs
- wb_pois <- wb_layers$POIs %>% filter(!is.na(Type_WBOut))
- missing_wb <- filter(wbs, !wb_id %in% wb_pois$Type_WBOut)
- print(nrow(missing_wb))
-
- double_wb <- wbs %>% group_by(wb_id) %>% filter(n() > 1)
- write_sf(double_wb, temp_gpkg, "double_wb")
-
- double_wb_pois <- wb_pois %>%
- group_by(Type_WBOut) %>%
- filter(n() > 1)
- write_sf(double_wb_pois, temp_gpkg, "double_wb_pois")
-
- wb_missing_poi <- filter(wbs,
- !wb_id %in% wb_pois$Type_WBOut)
- write_sf(wb_missing_poi, temp_gpkg, "wb_missing_poi")
-
- poi_missing_wb <- filter(wb_pois, !Type_WBOut %in% wb_layers$WBs$wb_id)
- write_sf(poi_missing_wb, temp_gpkg, "poi_missing_wb")
-
- #*************************************
- if(all(!is.na(wb_layers$events))){
- wb_events <- wb_layers$events %>%
- select(COMID, REACHCODE, REACH_meas, POI_identifier = wb_id,
- event_type, nexus) %>%
- st_compatibalize(., events)
-
- events <- rbind(events, wb_events)
- write_sf(events, temp_gpkg, split_layer)
- }
-
- if(all(!is.na(wb_layers$events))){
- # Write out events and outlets
- if(!needs_layer(temp_gpkg, split_layer)){
- events <- read_sf(temp_gpkg, split_layer) %>%
- rbind(st_compatibalize(wb_events,.)) %>%
- unique()
-
- write_sf(events, temp_gpkg, split_layer)
- } else {
- write_sf(wb_events, nav_gpkg, split_layer)
- }
- }
-
- tmp_POIs <- wb_layers$POIs
- # Write out waterbodies
- WB_VPU_pois <- wb_layers$WBs
- } else {
- tmp_POIs <- wb_layers
- WB_VPU_pois <- final_wb_list
- }
-
- #12047928
- wb_out_col <- wb_poi_collapse(tmp_POIs, nhd, events)
- tmp_POIs <- wb_out_col$POIs
-
- if(all(!is.na(wb_out_col$events_ret))){
- write_sf(wb_out_col$events_ret, temp_gpkg, split_layer)
- }
-
- nhd <- wb_layers$nhd_wb_table
- write_sf(nhd, nav_gpkg, nhd_flowline)
- write_sf(WB_VPU_pois, temp_gpkg, WBs_layer)
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
- } else {
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
- nhd <- read_sf(nav_gpkg, nhd_flowline)
- WB_VPU_pois <- read_sf(temp_gpkg, WBs_layer)
-}
-
-#mapview(filter(tmp_POIs, !is.na(Type_WBOut)), layer.name = "Waterbody outlet POIs", col.regions = "red")
-```
-
-```{r AR POIs}
-# Derive or load Waterbody POIs ----------------------
-if(!"Type_AR" %in% names(tmp_POIs)) {
- # Sparrow AR Events
- ar_events <- readxl::read_xlsx(file.path(data_paths$USGS_IT_path, "ARevents_final7.xlsx")) %>%
- filter(fromComid > 0 | ToComid > 0) %>%
- select(COMID = NHDComid, AREventID) %>%
- filter(COMID %in% nhd$COMID)
-
- if(nrow(ar_events) > 0){
- # AR POIs
- tmp_POIs <- POI_creation(ar_events, nhd, "AR") %>%
- addType(., tmp_POIs, "AR", nexus = TRUE)
-
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
- } else {
- print(paste0("No AR events in ", vpu_codes))
- }
-
-} else {
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-}
-```
-
-```{r Terminal POIs}
-# Derive or load Terminal POIs ----------------------
-if(!"Type_Term" %in% names(tmp_POIs)) {
-
- # Terminal POIs on levelpaths with upstream POIs
- term_paths <- filter(st_drop_geometry(nhd), Hydroseq %in% filter(nhd, COMID %in% tmp_POIs$COMID)$TerminalPa)
-
- # Non-POI levelpath terminal pois, but meet size criteria
- terminal_POIs <- st_drop_geometry(nhd) %>%
- filter(Hydroseq == TerminalPa | (toCOMID == 0 | is.na(toCOMID) | !toCOMID %in% COMID)) %>%
- filter(!COMID %in% term_paths$COMID, TotDASqKM >= min_da_km_terminal) %>%
- 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", nexus = TRUE)
-
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-} else {
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
- nhd <- read_sf(nav_gpkg, nhd_flowline)
-}
-
-mapview(filter(tmp_POIs, !is.na(Type_Term)), layer.name = "Terminal POIs", col.regions = "blue")
-```
-
-```{r Confluence POIs}
-# Derive POIs at confluences where they are absent ----------------------
-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), \(x, nhdDF) get_UM(nhdDF, x), nhdDF = nhd)))
- finalNet <- unique(NetworkConnection(up_net, st_drop_geometry(nhd)))
-
- # Subset NHDPlusV2 flowlines to navigation results and write to shapefile
- nhd <- mutate(nhd, minNet = ifelse(COMID %in% finalNet, 1, 0))
-
- # Create new confluence POIs
- conf_COMIDs <- st_drop_geometry(filter(nhd, minNet == 1)) %>%
- # Downstream hydrosequence of 0 indicates
- # the flowline is terminating or
- # leaving the domain, so they
- # are excluded from this process
- filter(DnHydroseq > 0) %>%
- group_by(DnHydroseq) %>%
- filter(n()> 1) %>%
- mutate(Type_Conf = LevelPathI) %>%
- ungroup() %>%
- select(COMID, Type_Conf)
-
- tmp_POIs <- POI_creation(conf_COMIDs, filter(nhd, minNet == 1), "Conf") %>%
- addType(., tmp_POIs, "Conf", nexus = TRUE)
-
- write_sf(nhd, nav_gpkg, nhd_flowline)
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-} else {
- nhd <- read_sf(nav_gpkg, nhd_flowline)
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-}
-
-mapview(filter(tmp_POIs, !is.na(Type_Conf)), layer.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(filter(nhd, minNet == 1), WB_VPU_pois,
- data_paths, proj_crs, split_layer, tmp_POIs)
- wb_in_col <- wb_inlet_collapse(wb_layers$POIs, nhd, events)
- tmp_POIs <- wb_in_col$POIs
-
- if(!all(is.na(wb_layers$events))) {
- wb_inlet_events <- wb_layers$events %>%
- select(COMID, REACHCODE, REACH_meas, POI_identifier,
- event_type, nexus)
-
- # 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)
- #tmp_POIs <- wb_layers$POIs
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-} else {
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
- nhd <- read_sf(nav_gpkg, nhd_flowline)
-}
-
-mapview(filter(tmp_POIs, !is.na(Type_WBIn)), layer.name = "Waterbody inlet POIs", col.regions = "blue") +
- mapview(filter(tmp_POIs, !is.na(Type_WBOut)), layer.name = "Waterbody outlet POIs", col.regions = "red")
-```
-
-```{r NID POIs}
-# Derive or load NID POIs ----------------------
-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(COMID) %>%
- summarize(Type_NID = paste0(unique(NIDID), collapse = " "))
-
- # Derive other NID POIs
- tmp_POIs <- POI_creation(NID_COMIDs, nhd, "NID") %>%
- addType(., tmp_POIs, "NID", nexus = FALSE, bind = FALSE)
-
- # Write out geopackage layer representing POIs for given theme
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-} else {
- # Load NID POIs if they already exist
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-}
-
-mapview(filter(tmp_POIs, Type_NID != ""), layer.name = "NID POIs", col.regions = "blue")
-```
-
-```{r HW DA POIs}
-if(!"Type_DA" %in% names(tmp_POIs)){
-
- # derive incremental segments from POIs
- inc_segs <- make_incremental_segments(nhd,
- filter(st_drop_geometry(nhd),
- COMID %in% tmp_POIs$COMID)) %>%
- # bring over VAA data
- inner_join(select(atts_rds, COMID,
- DnHydroseq, VA_MA, TOTMA, LENGTHKM, MAXELEVSMO,
- MINELEVSMO, WBAREACOMI, WBAreaType, FTYPE, StartFlag,
- AreaSqKM, TotDASqKM), by = "COMID")
-
- hw_segs <- inc_segs %>%
- group_by(POI_ID) %>%
- filter(any(StartFlag == 1)) %>%
- filter(any(TotDASqKM > max_da_km_hw)) %>%
- ungroup()
-
- att_group <- function(a, athres) {
- #cumsum <- 0
- group <- 1
- result <- numeric()
- for (i in 1:length(a)) {
- #cumsum <- cumsum + a[i]
- tot_DA <- a[i]
- if (tot_DA > athres) {
- group <- group + 1
- athres <- athres + athres
- }
- result = c(result, group)
- }
- return (result)
- }
-
- #TODO Add magic numbers to config file
- hw_DA_splits <- hw_segs %>%
- st_drop_geometry() %>%
- #filter(.data$ID %in% cat$ID) %>%
- group_by(POI_ID) %>%
- arrange(-Hydroseq) %>%
- mutate(ind = att_group(TotDASqKM, min_da_km_hw)) %>%
- ungroup() %>%
- group_by(POI_ID, ind) %>%
- mutate(total_length = cumsum(LENGTHKM)) %>%
- ungroup() %>%
- group_by(POI_ID, ind) %>%
- mutate(set = cur_group_id()) %>%
- filter(TotDASqKM == max(TotDASqKM) &
- total_length > 5) %>%
- mutate(da_class = paste(POI_ID, ind, sep = "_")) %>%
- ungroup() %>%
- select(-ind) %>%
- filter(!COMID %in% tmp_POIs$COMID)
-
- if(nrow(hw_DA_splits) > 0) {
- tmp_POIs <- POI_creation(select(hw_DA_splits, COMID, da_class),
- filter(nhd, poi == 1), "DA") %>%
- addType(., tmp_POIs, "DA", nexus = TRUE)
- }
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-} else {
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-}
-```
-
-This chunk breaks up potential aggregated segments where the elevation range of the are contributing to the segment exceeds 500-m
-```{r Elevation Break POIs}
-if(!"Type_elev" %in% names(tmp_POIs)){
-
- inc_segs <- make_incremental_segments(filter(nhd, minNet == 1),
- filter(st_drop_geometry(nhd),
- COMID %in% tmp_POIs$COMID)) %>%
- # bring over VAA data
- inner_join(select(atts_rds, COMID,
- DnHydroseq, VA_MA, TOTMA, LENGTHKM, MAXELEVSMO,
- MINELEVSMO, WBAREACOMI, WBAreaType, FTYPE, StartFlag,
- AreaSqKM, TotDASqKM), by = "COMID")
-
- elev_fp <- inc_segs %>%
- group_by(POI_ID) %>%
- arrange(Hydroseq) %>%
- # Get elevation info
- mutate(MAXELEVSMO = na_if(MAXELEVSMO, -9998), MINELEVSMO = na_if(MINELEVSMO, -9998),
- elev_diff_seg = max(MAXELEVSMO) - min(MINELEVSMO),
- total_length = cumsum(LENGTHKM)) %>%
- filter((max(MINELEVSMO) - min(MINELEVSMO)) > elev_diff) %>%
- mutate(inc_elev_diff = c(MINELEVSMO[1], (MINELEVSMO - lag(MINELEVSMO))[-1])) %>%
- mutate(inc_elev_diff = ifelse(inc_elev_diff == MINELEVSMO, 0, inc_elev_diff)) %>%
- ungroup()
-
- cs_group <- function(a, athres) {
- cumsum <- 0
- group <- 1
- result <- numeric()
- for (i in 1:length(a)) {
- cumsum <- cumsum + a[i]
- if (cumsum > athres) {
- group <- group + 1
- cumsum <- a[i]
- }
- result = c(result, group)
- }
- return (result)
- }
-
- if(nrow(elev_fp) > 0){
- #TODO Add magic numbers to config file
- hw_elev_splits <- elev_fp %>%
- st_drop_geometry() %>%
- #filter(.data$ID %in% cat$ID) %>%
- group_by(POI_ID) %>%
- arrange(-Hydroseq) %>%
- mutate(ind = cs_group(inc_elev_diff, elev_diff/2)) %>%
- ungroup() %>%
- group_by(POI_ID, ind) %>%
- filter(TotDASqKM == max(TotDASqKM) &
- total_length > 5) %>%
- mutate(elev_class = paste(POI_ID, ind, sep = "_")) %>%
- ungroup() %>%
- select(-ind) %>%
- filter(!COMID %in% tmp_POIs$COMID)
-
- if(nrow(hw_elev_splits) > 0) {
- tmp_POIs <- POI_creation(select(hw_elev_splits, COMID, elev_class),
- filter(nhd, poi == 1),
- "elev") %>%
- addType(., tmp_POIs, "elev", nexus = TRUE)
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
- }
- }
-
-
-} else {
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-}
-
-```
-
-```{r Time of travel POIs}
-if(all(is.na(tmp_POIs$Type_Travel))) {
-
- # derive incremental segments from POIs
- inc_segs <- make_incremental_segments(nhd,
- filter(st_drop_geometry(nhd),
- COMID %in% tmp_POIs$COMID,
- COMID %in% nhd$COMID)) %>%
- # bring over VAA data
- inner_join(select(atts_rds, COMID, DnHydroseq, VA_MA, TOTMA, LENGTHKM,
- MAXELEVSMO, MINELEVSMO, WBAREACOMI, WBAreaType, FTYPE,
- AreaSqKM, TotDASqKM), by = "COMID")
-
- # TT POIs
- tt_pois_split <- inc_segs %>%
- arrange(-Hydroseq) %>%
- # Should we substitute with a very small value
- mutate(VA_MA = ifelse(VA_MA < 0, NA, VA_MA)) %>%
- mutate(FL_tt_hrs = (LENGTHKM * ft_per_km)/ VA_MA / s_per_hr ) %>%
- group_by(POI_ID) %>%
- filter(sum(FL_tt_hrs) > travt_diff, max(TotDASqKM) > max_elev_TT_DA) %>%
- mutate(cum_tt = cumsum(FL_tt_hrs),
- hrs = sum(FL_tt_hrs))
-
- #TODO Add magic numbers to config file
- tt_splits <- tt_pois_split %>%
- st_drop_geometry() %>%
- group_by(POI_ID) %>%
- arrange(-Hydroseq) %>%
- mutate(ind = cs_group(FL_tt_hrs, travt_diff * 0.75)) %>%
- ungroup() %>%
- group_by(POI_ID, ind) %>%
- filter(TotDASqKM == max(TotDASqKM)) %>%
- mutate(tt_class = paste(POI_ID, ind, sep = "_")) %>%
- ungroup() %>%
- select(-ind) %>%
- filter(!COMID %in% tmp_POIs$COMID)
-
- if(nrow(tt_splits) > 0) {
- tmp_POIs <- POI_creation(select(tt_splits, COMID, tt_class),
- filter(nhd, poi == 1), "Travel") %>%
- addType(., tmp_POIs, "Travel", nexus = TRUE)
-
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
- }
-
-}else{
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-}
-
-```
-
-```{r Final POIs}
-# Derive final POI set ----------------------
-if(needs_layer(temp_gpkg, pois_all_layer)) {
-
- unCon_POIs <- filter(tmp_POIs, COMID %in% filter(nhd, AreaSqKM == 0)$COMID)
-
- # If any POIs happened to fall on flowlines w/o catchment
- if (nrow(unCon_POIs) >0){
- # For confluence POIs falling on Flowlines w/o catchments, derive upstream valid flowline,
- poi_fix <- DS_poiFix(tmp_POIs, filter(nhd, minNet == 1))
- new_POIs <- st_compatibalize(poi_fix$new_POIs, tmp_POIs)
- xWalk <- poi_fix$xWalk
-
- # POIs that didn't need to be moved
- tmp_POIs_fixed <- filter(tmp_POIs, nexus == TRUE |
- !COMID %in% c(poi_fix$xWalk$oldPOI, poi_fix$xWalk$COMID))
- # bind together
- final_POIs <- bind_rows(tmp_POIs_fixed, new_POIs) %>%
- mutate(Type_Term = ifelse(nexus == 1, NA, Type_Term)) %>%
- select(-dplyr::any_of(c("ID")))
-
- # Write out fixes
- write_sf(new_POIs, temp_gpkg, poi_moved_layer)
- write_sf(xWalk, temp_gpkg, poi_xwalk_layer)
- # write out final POIs
- write_sf(final_POIs, temp_gpkg, pois_all_layer)
- } else {
- # If no fixes designate as NA
- poi_fix <- NA
-
- tmp_POIs$nexus <- as.integer(tmp_POIs$nexus)
-
- # if a POI will be a nexus, it can not be terminal.
- tmp_POIs <- mutate(tmp_POIs, Type_Term = ifelse(nexus == 1, NA, Type_Term))
-
- # write out final POIs
- write_sf(tmp_POIs, temp_gpkg, pois_all_layer)
-
- }
-
-} else {
- # Need all three sets for attribution below
- final_POIs <- read_sf(temp_gpkg, pois_all_layer)
- new_POIs <- if(!needs_layer(temp_gpkg, poi_moved_layer)) read_sf(temp_gpkg, poi_moved_layer) else (NA)
- xWalk <- if(!needs_layer(temp_gpkg, poi_xwalk_layer)) read_sf(temp_gpkg, poi_xwalk_layer) else (NA)
- unCon_POIs <- filter(st_drop_geometry(filter(nhd, minNet == 1)),
- COMID %in% tmp_POIs$COMID, AreaSqKM == 0)
-}
-```
-
-```{r Draft nsegment generation}
-### NEEDED?
-
-if(needs_layer(temp_gpkg, nsegments_layer)) {
-
- if("POI_ID" %in% colnames(nhd)){
- nhd <- select(nhd, -POI_ID)
- }
- # Sort POIs by Levelpath and Hydrosequence in upstream to downstream order
- seg_POIs <- filter(st_drop_geometry(nhd),
- COMID %in% tmp_POIs$COMID, COMID %in% filter(nhd, minNet == 1)$COMID)
- # derive incremental segments from POIs
- inc_segs <- make_incremental_segments(filter(nhd, minNet == 1), seg_POIs)
-
- nhd_Final <- nhd %>%
- left_join(select(inc_segs, COMID, POI_ID), by = "COMID")
-
- # create and write out final dissolved segments
- nsegments_fin <- segment_creation(nhd_Final, xWalk)
-
- nhd_Final <- select(nhd_Final, -POI_ID) %>%
- left_join(distinct(st_drop_geometry(nsegments_fin$raw_segs), COMID, POI_ID), by = "COMID")
- nsegments <- nsegments_fin$diss_segs
-
- write_sf(nhd_Final, nav_gpkg, nhd_flowline)
- write_sf(nsegments, temp_gpkg, nsegments_layer)
-} else {
- # Read in NHDPlusV2 flowline simple features and filter by vector processing unit (VPU)
- final_POIs <- read_sf(temp_gpkg, pois_all_layer)
- nhd_Final <- read_sf(nav_gpkg, nhd_flowline)
- nsegments <- read_sf(temp_gpkg, nsegments_layer)
-}
-
-# Ensure that all the problem POIs have been taken care of
-sub <- nhd_Final %>% filter(COMID %in% final_POIs$COMID)
-print (paste0(nrow(sub[sub$AreaSqKM == 0,]), " POIs on flowlines with no local drainage contributions"))
-noDA_pois <- filter(final_POIs, COMID %in% filter(sub, AreaSqKM == 0)$COMID)
-write_sf(noDA_pois, temp_gpkg, "noDA_pois")
-```
-
-```{r POI Collapse}
-# number POIs
-final_POIs_prec <- mutate(final_POIs, id = row_number(), moved = NA) %>%
- write_sf(temp_gpkg, pois_all_layer) %>%
- write_sf(nav_gpkg, nav_poi_layer)
-collapse <- TRUE
-
-# TODO: document what this is doing!
-source("R/poi_move.R")
-
-check_dups <- final_POIs %>%
- group_by(COMID) %>%
- filter(n() > 1)
-
-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, !id %in% no_dups$id)
- write_sf(dups, temp_gpkg, dup_pois)
-} else {
- print("All double COMIDs nexus for gage or WB splitting")
-}
-
-mapview(final_POIs, layer.name = "All POIs", col.regions = "red") +
- mapview(nsegments, layer.name = "Nsegments", col.regions = "blue")
-```
-
-```{r lookup}
-# Final POI layer
-POIs <- final_POIs %>%
- 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,
- -c(nexus, nonrefactor, id, moved))),
- 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()
-
-if(exists("moved_pois_table")){
- pois_data_moved <- select(st_drop_geometry(moved_pois_table),
- hy_id = COMID, hl_link = new_val, hl_reference = moved_value) %>%
- inner_join(distinct(pois_data_orig, hy_id, geom_id, poi_id), by = "hy_id")
-
- pois_data <- data.table::rbindlist(list(pois_data_moved, pois_data_orig), use.names = TRUE) %>%
- filter(!hl_reference %in% c("id", "moved"))
-} else {
- pois_data <- pois_data_orig
-}
-
-# POI Geometry table
-poi_geometry <- select(final_POIs_table, hy_id = COMID, geom_id) %>%
- 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,
- realized_catchmentID = COMID,
- mainstem = LevelPathI) %>%
- dplyr::mutate(realized_catchmentID = ifelse(realized_catchmentID %in% full_cats$FEATUREID,
- realized_catchmentID, NA)) %>%
- left_join(select(st_drop_geometry(poi_geometry), 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/02_NHD_navigate_Sp.Rmd b/workspace/02_NHD_navigate_Sp.Rmd
deleted file mode 100644
index a18d0f1..0000000
--- a/workspace/02_NHD_navigate_Sp.Rmd
+++ /dev/null
@@ -1,822 +0,0 @@
----
-title: "NHD Navigate"
-output: html_document
-editor_options:
- chunk_output_type: console
----
-This notebook Generate Segments and POIS from POI data sources and builds
-a minimally-sufficient stream network connecting all POis, as well as generating first-cut of
-national segment network.
-
-```{r setup_rmd, echo=FALSE, cache=FALSE}
-knitr::opts_chunk$set(
- collapse = TRUE,
- comment = "#>",
- fig.width=6,
- fig.height=4,
- cache=FALSE)
-```
-
-```{r setup}
-# Load custom functions
-source("R/utils.R")
-source("R/NHD_navigate.R")
-
-# increase timeout for data downloads
-options(timeout=600)
-
-# Load Configuration of environment
-source("R/config.R")
-
-# Gages output from Gage_selection
-gages <- read_sf(gage_info_gpkg, "Gages")
-# need some extra attributes for a few POI analyses
-atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.rds"))
-
-sparrow_dir <- "data/sparrow"
-
-collapse <- TRUE
-```
-
-```{r huc12 POIs}
-# Derive or load HUC12 POIs
-if(needs_layer(temp_gpkg, nav_poi_layer)) {
-
- nhd <- read_sf(nav_gpkg, nhd_flowline)
- try(nhd <- select(nhd, -c(minNet, WB, struct_POI, struct_Net, POI_ID, dend, poi)), silent = TRUE)
-
- # Some NHDPlus VPUs include HUCs from other VPUs
- if(vpu == "02"){
- grep_exp <-"^02|^04"
- } else if (vpu == "08") {
- grep_exp <- "^03|^08"
- } else {
- grep_exp <- paste0("^", substr(vpu, start = 1, stop = 2))
- }
-
- # Join HUC12 outlets with NHD
- HUC12_COMIDs <- read_sf(data_paths$hu12_points_path, "hu_points") %>%
- filter(grepl(grep_exp, .data$HUC12)) %>%
- select(COMID, HUC12) %>%
- # Remove this when HUC12 outlets finished
- group_by(COMID) %>%
- filter(row_number() == 1) %>%
- ungroup()
-
- # Create POIs - some r05 HUC12 POIs not in R05 NHD
- huc12_POIs <- POI_creation(st_drop_geometry(HUC12_COMIDs), nhd, "HUC12")
- tmp_POIs <- huc12_POIs
-
- # Write out geopackage layer representing POIs for given theme
- write_sf(huc12_POIs, temp_gpkg, nav_poi_layer)
- tmp_POIs <- huc12_POIs
-} else {
- # Load HUC12 POIs as the tmpPOIs if they already exist
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
- nhd <- read_sf(nav_gpkg, nhd_flowline)
-}
-
-mapview(filter(tmp_POIs, Type_HUC12 != ""), layer.name = "HUC12 POIs", col.regions = "red")
-```
-
-```{r streamgage POIs}
-if(!"Type_Gages" %in% names(tmp_POIs)) {
-
- ref_gages <- read_sf("data/ref_gages.gpkg")
-
- # Previously identified streamgages within Gage_Selection.Rmd
- streamgages_VPU <- gages %>%
- rename(COMID = comid) %>%
- filter(COMID %in% nhd$COMID) %>%
- #st_drop_geometry() %>%
- switchDiv(., nhd)
-
- streamgages <- streamgages_VPU %>%
- group_by(COMID) %>%
- # If multiple gages per COMID, pick one with highest rank as rep POI_ID
- filter(gage_score == max(gage_score), !is.na(drain_area)) %>%
- ungroup() %>%
- select(COMID, site_no)
-
- # Derive TE POIs
- tmp_POIs <- POI_creation(st_drop_geometry(streamgages), nhd, "Gages") %>%
- addType(., tmp_POIs, "Gages", nexus = FALSE)
-
- # Write out geopackage layer representing POIs for given theme
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-} else {
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-}
-
-mapview(filter(tmp_POIs, !is.na(Type_Gages)), layer.name = "Streamgage POIs", col.regions = "blue")
-```
-
-# R01 - missing one
-```{r Sparrow Cal POIs}
-if(!"Type_Sparrow" %in% names(tmp_POIs)) {
-
- # Load Sparrow calibration POIs
- sparrow_cal_pois <- read_sf("data/sparrow/sparrow_source.gpkg", "calibration_poi") %>%
- mutate(COMID = as.integer(COMID)) %>%
- filter(COMID %in% nhd$COMID) %>%
- select(COMID, MonitoringLocationIdentifier, site_no, station_nm,
- usgs_hy_id, distance_m, longitude, latitude, nat_id) %>%
- st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>%
- select(COMID, nat_id)
-
- # Derive TE POIs
- tmp_POIs <- POI_creation(st_drop_geometry(sparrow_cal_pois), nhd, "Sparrow") %>%
- addType(., tmp_POIs, "Sparrow", nexus = FALSE)
-
- missing_sparrow <- filter(sparrow_cal_pois, !nat_id %in% tmp_POIs$Type_Sparrow)
- if(nrow(missing_sparrow) > 0){write_sf(missing_sparrow, temp_gpkg, "missing_sparrow")}
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-} else {
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-}
-
-mapview(filter(tmp_POIs, !is.na(Type_Gages)), layer.name = "Streamgage POIs", col.regions = "blue")
-```
-
-```{r TE POIs}
-if(!"Type_TE" %in% names(tmp_POIs)) {
-
- if(vpu == "08"){
- nhd$VPUID <- "08"
- } else {
- nhd$VPUID <- substr(nhd$RPUID, 1, 2)
- }
-
- # Read in Thermoelectric shapefile
- TE_COMIDs <- read_sf(data_paths$TE_points_path, "2015_TE_Model_Estimates_lat.long_COMIDs") %>%
- mutate(COMID = as.integer(COMID)) %>%
- inner_join(., select(st_drop_geometry(nhd), COMID, VPUID), by = "COMID") %>%
- filter(grepl(paste0("^", substr(vpu, 1, 2), ".*"), .data$VPUID), COMID > 0) %>%
- switchDiv(., nhd) %>%
- group_by(COMID) %>%
- summarize(EIA_PLANT = paste0(unique(EIA_PLANT_), collapse = " "), count = n()) %>%
- ungroup() %>%
- select(COMID, EIA_PLANT)
-
- # Derive TE POIs
- tmp_POIs <- POI_creation(st_drop_geometry(TE_COMIDs), nhd, "TE") %>%
- addType(., tmp_POIs, "TE", nexus = FALSE, bind = TRUE)
-
- # As a fail-safe, write out list of TE plants not assigned a POI
- if(nrow(filter(TE_COMIDs, !COMID %in% tmp_POIs$COMID)) > 0) {
- write_sf(filter(TE_COMIDs, !COMID %in% tmp_POIs$COMID),
- temp_gpkg, "unassigned_TE")
- }
-
- # Write out geopackage layer representing POIs for given theme
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-} else {
- # Load TE POIs if they already exist
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-}
-
-mapview(filter(tmp_POIs, Type_TE != ""), layer.name = "TE Plant POIs", col.regions = "blue")
-```
-
-```{r waterbody outlet POIs}
-# Derive or load Waterbody POIs ----------------------
-if(!"Type_WBOut" %in% names(tmp_POIs)) {
-
- sparrow_wb_outlet <- read_sf("data/sparrow/sparrow_source.gpkg", "wb_outlets") %>%
- filter(COMID %in% nhd$COMID)
-
- sparrow_wb <- read_sf("data/sparrow/sparrow_source.gpkg", "waterbodies") %>%
- filter(COMID %in% nhd$WBAREACOMI)
-
- # Waterbodies sourced from NHD waterbody layer for entire VPU
- WBs_VPU_all <- filter(readRDS("data/NHDPlusNationalData/nhdplus_waterbodies.rds"),
- COMID %in% nhd$WBAREACOMI) %>%
- mutate(FTYPE = as.character(FTYPE)) %>%
- mutate(source = "NHDv2WB") %>%
- st_transform(crs)
-
- # Waterbody list that strictly meet the size criteria
- wb_lst <- st_drop_geometry(WBs_VPU_all) %>%
- dplyr::filter(FTYPE %in% c("LakePond", "Reservoir") & AREASQKM >= (min_da_km/2)) %>%
- dplyr::select(COMID) %>%
- filter(!COMID == 166410600) %>%
- rbind(dplyr::select(st_drop_geometry(sparrow_wb), COMID))
-
- if (file.exists(data_paths$resops_NID_CW)){
- # ResOpsUS locations in the VPU waterbody set
- resops_wb_df <- read.csv(data_paths$resops_NID_CW) %>%
- dplyr::filter(WBAREACOMI %in% WBs_VPU_all$COMID) %>%
- dplyr::select(DAM_ID, DAM_NAME, COMID = WBAREACOMI, FlowLcomid) %>%
- distinct()
-
- # Add ResOPsUS locations to waterbody list
- if(nrow(resops_wb_df) > 0){
- wb_lst <- rbind(wb_lst, select(resops_wb_df, COMID)) %>%
- distinct()
- }
- }
-
- # Waterbodies that meet the size criteria and/or are in ResOpsUS datasets
- WBs_VPU <- WBs_VPU_all %>%
- dplyr::filter(COMID %in% wb_lst$COMID)
-
- WBs_VPU <- filter(WBs_VPU, COMID != 666170) #1R16
-
- # Attribute flowlines that are in waterbodies
- nhd <- mutate(nhd, WB = ifelse(WBAREACOMI > 0 & WBAREACOMI %in% WBs_VPU$COMID, 1, 0))
-
- wbout_COMIDs <- nhd %>% #filter(nhd, dend == 1 & WB == 1)
- filter(WBAREACOMI %in% WBs_VPU$COMID) %>%
- filter(!WBAREACOMI %in% sparrow_wb_outlet$WB_COMID) %>%
- group_by(WBAREACOMI) %>%
- slice(which.min(Hydroseq)) %>%
- switchDiv(., nhd) %>%
- select(COMID, WBAREACOMI, LevelPathI) %>%
- st_drop_geometry() %>%
- rbind(select(st_drop_geometry(sparrow_wb_outlet), COMID, WBAREACOMI = WB_COMID))
-
- tmp_POIs <- POI_creation(wbout_COMIDs, nhd, "WBOut") %>%
- addType(., tmp_POIs, "WBOut", nexus = FALSE, bind = TRUE)
-
- # Write out waterbodies
- write_sf(WBs_VPU, nav_gpkg, WBs_layer)
-
- # Write specific ResOpsUS data to a separate sheet
- if(nrow(wb_lst) > 0){
- resops_POIs_df <- st_drop_geometry(tmp_POIs) %>%
- dplyr::select(COMID, Type_WBOut) %>%
- dplyr::filter(!is.na(Type_WBOut)) %>%
- dplyr::inner_join(select(resops_wb_df, -FlowLcomid), by = c("Type_WBOut" = "COMID"))
-
- write.csv(resops_POIs_df, file.path("data/reservoir_Data", paste0("resops_POIs_",vpu,".csv")))
- }
-
-
- wb_out <- filter(tmp_POIs, !is.na(Type_WBOut))
- missing_wb_out <- filter(sparrow_wb_outlet, !WB_COMID %in% wb_out$Type_WBOut)
- missing_wb_out_com <- filter(sparrow_wb_outlet,
- !COMID %in% wb_out$COMID)
- if(nrow(missing_wb_out) > 0){write_sf(missing_wb_out, temp_gpkg, "missing_sparrow_WB")}
- if(nrow(missing_wb_out_com) > 0){write_sf(missing_wb_out_com, temp_gpkg, "missing_sparrow_WBout")}
-
- sparrow_wb_ids <- select(sparrow_wb_outlet, COMID, WB_COMID, nat_id)
-
- tmp_POIs <- tmp_POIs %>%
- left_join(select(st_drop_geometry(sparrow_wb_ids), -COMID), by = c("Type_WBOut" = "WB_COMID")) %>%
- mutate(Type_Sparrow = ifelse(is.na(Type_Sparrow), nat_id, Type_Sparrow)) %>%
- select(-nat_id)
-
- write_sf(nhd, nav_gpkg, nhd_flowline)
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-} else {
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
- nhd <- read_sf(nav_gpkg, nhd_flowline)
- WBs_VPU <- read_sf(nav_gpkg, WBs_layer)
- resops_POIs_df <- read.csv(file.path("data/reservoir_Data", paste0("resops_POIs_",vpu,".csv")))
-}
-
-mapview(filter(tmp_POIs, !is.na(Type_WBOut)), layer.name = "Waterbody outlet POIs", col.regions = "red")
-```
-
-```{r NID POIs}
-# Derive or load NID POIs ----------------------
-if(all(is.na(tmp_POIs$Type_NID))) {
-
- sparrow_nid <- read_sf("data/sparrow/sparrow_source.gpkg", "nid_pois")
-
- sparrow_wb_nid <- sparrow_nid %>%
- #rename(COMID = ComID) %>%
- filter(COMID %in% nhd$COMID) %>%
- select(COMID, NIDID, nat_id)
-
- if(nrow(sparrow_wb_nid) > 0){
- # Determine which NID facilitites have been co-located with ResOps
- res_ops_table <- read.csv(data_paths$resops_NID_CW) %>%
- mutate(new_WBARECOMI = ifelse(is.na(WBAREACOMI) & is.na(NHDAREA), HR_NHDPLUSID, WBAREACOMI)) %>%
- mutate(new_WBARECOMI = ifelse(is.na(new_WBARECOMI) & is.na(HR_NHDPLUSID), NHDAREA, new_WBARECOMI)) %>%
- filter(new_WBARECOMI %in% tmp_POIs$Type_WBOut, !is.na(new_WBARECOMI)) %>%
- select(NID_ID = NID_Source_FeatureID, Type_WBOut = new_WBARECOMI) %>%
- filter()
-
- # Attribute those NID facilities precisely at the outlet
- if(nrow(res_ops_table) > 0){
- tmp_POIs <- tmp_POIs %>%
- left_join(res_ops_table, by = "Type_WBOut") %>%
- mutate(Type_NID = ifelse(!is.na(NID_ID), NID_ID, NA)) %>%
- select(-NID_ID)
-
- sparrow_wb_nid <- filter(sparrow_wb_nid, !NIDID %in% tmp_POIs$Type_NID)
- }
-
- # Derive other NID POIs
- tmp_POIs <- POI_creation(st_drop_geometry(sparrow_wb_nid), nhd, "NID") %>%
- addType(., tmp_POIs, "NID", nexus = FALSE, bind = TRUE) %>%
- left_join(select(st_drop_geometry(sparrow_wb_nid), NIDID, nat_id), by = c("Type_NID" = "NIDID")) %>%
- mutate(Type_Sparrow = ifelse(is.na(Type_Sparrow), nat_id, Type_Sparrow)) %>%
- select(-nat_id)
-
- sparrow_wb_nid_missing <- filter(sparrow_wb_nid, !nat_id %in% tmp_POIs$Type_Sparrow) %>%
- mutate(WB_Comid = paste("Sparrow_nid_", row_number()))
-
- if(nrow(sparrow_wb_nid_missing) > 0){write_sf(sparrow_wb_nid_missing, temp_gpkg, "missing_wb_nid")}
- # Write out geopackage layer representing POIs for given theme
- write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
- mapview(filter(tmp_POIs, Type_NID != ""), layer.name = "NID POIs", col.regions = "blue")
- }
-
-} else {
- # Load NID POIs if they already exist
- tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-}
-```
-
-
-```{r Diversion POIs}
-# # Derive POIs at confluences where they are absent ----------------------
-# if(needs_layer(temp_gpkg, "diversions")) {
-#
-# div_pois_source <- read_sf("data/sparrow/sparrow_source.gpkg", "div_poi") %>%
-# filter(COMID %in% nhd$COMID) %>%
-# select(COMID, nat_id)
-#
-# div_POIs <- POI_creation(st_drop_geometry(div_pois_source), nhd, "Div") %>%
-# mutate(Type_Sparrow = NA, nexus = NA) %>%
-# st_compatibalize(., tmp_POIs)
-#
-# tmp_POIs <- data.table::rbindlist(list(tmp_POIs,
-# select(div_POIs, -c(Type_Div))), fill = TRUE) %>%
-# st_as_sf() %>%
-# left_join(select(st_drop_geometry(div_pois_source), COMID, nat_id), by = "COMID") %>%
-# mutate(Type_Sparrow = ifelse(is.na(Type_Sparrow), nat_id, Type_Sparrow)) %>%
-# select(-nat_id)
-#
-# missing_div <- filter(div_pois_source, !nat_id %in% tmp_POIs$Type_Sparrow)
-# if(nrow(missing_div) > 0){write_sf(missing_div, temp_gpkg, "missing_div")}
-#
-# write_sf(div_pois_source, temp_gpkg, "diversions")
-# # Write out geopackage layer representing POIs for given theme
-# write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-# } else {
-# # Load NID POIs if they already exist
-# tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-# }
-```
-
-```{r Confluence POIs}
-# # Derive POIs at confluences where they are absent ----------------------
-# if(all(is.na(tmp_POIs$Type_Conf)) | !"Type_Conf" %in% colnames(tmp_POIs)) {
-#
-# # Navigate upstream from each POI and determine minimally-sufficient network between current POI sets
-# up_net <- unique(unlist(lapply(unique(tmp_POIs$COMID), \(x, nhdDF) get_UM(nhdDF, x), nhdDF = nhd)))
-# finalNet <- unique(NetworkConnection(up_net, st_drop_geometry(nhd)))
-#
-# # Subset NHDPlusV2 flowlines to navigation results and write to shapefile
-# nhd <- mutate(nhd, minNet = ifelse(COMID %in% finalNet, 1, 0))
-#
-# # Create new confluence POIs
-# conf_COMIDs <- st_drop_geometry(filter(nhd, minNet == 1)) %>%
-# # Downstream hydrosequence of 0 indicates
-# # the flowline is terminating or
-# # leaving the domain, so they
-# # are excluded from this process
-# filter(DnHydroseq > 0) %>%
-# group_by(DnHydroseq) %>%
-# filter(n()> 1) %>%
-# mutate(Type_Conf = LevelPathI) %>%
-# ungroup() %>%
-# select(COMID, Type_Conf)
-#
-# tmp_POIs <- POI_creation2(conf_COMIDs, filter(nhd, minNet == 1), "Conf") %>%
-# addType2(., tmp_POIs, "Conf", nexus = FALSE, bind = TRUE)
-#
-# write_sf(nhd, nav_gpkg, nhd_flowline)
-# write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-# } else {
-# nhd <- read_sf(nav_gpkg, nhd_flowline)
-# tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-# }
-#
-# mapview(filter(tmp_POIs, !is.na(Type_Conf)), layer.nameble = "Confluence POIs", col.regions = "blue")
-```
-
-```{r Terminal POIs}
-# # Derive or load Terminal POIs ----------------------
-# if(all(is.na(tmp_POIs$Type_Term))) {
-#
-# # Terminal POIs on levelpaths with upstream POIs
-# term_paths <- filter(st_drop_geometry(nhd), Hydroseq %in% filter(nhd, COMID %in% tmp_POIs$COMID)$TerminalPa)
-#
-# # Non-POI levelpath terminal pois, but meet size criteria
-# terminal_POIs <- st_drop_geometry(nhd) %>%
-# filter(Hydroseq == TerminalPa | (toCOMID == 0 | is.na(toCOMID) | !toCOMID %in% COMID)) %>%
-# filter(!COMID %in% term_paths$COMID, TotDASqKM >= min_da_km) %>%
-# bind_rows(term_paths) %>%
-# # Use level path identifier as Type ID
-# select(COMID, LevelPathI)
-#
-# tmp_POIs <- POI_creation2(terminal_POIs, nhd, "Term") %>%
-# addType2(., tmp_POIs, "Term", nexus = FALSE, bind = TRUE)
-#
-# write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-# } else {
-# tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-# nhd <- read_sf(nav_gpkg, nhd_flowline)
-# }
-#
-# mapview(filter(tmp_POIs, !is.na(Type_Term)), layer.name = "Terminal POIs", col.regions = "blue")
-```
-
-```{r waterbody inlet POIs}
-# # Derive or load Waterbody POIs ----------------------
-# if(all(is.na(tmp_POIs$Type_WBIn))) {
-#
-# sparrow_inlets <- read_sf("data/sparrow/sparrow_source.gpkg", "wb_inlets") %>%
-# filter(COMID %in% nhd$COMID)
-#
-# WBs_VPU <- filter(WBs_VPU, COMID != 666170) #1R16
-#
-# wbout_COMIDs <- filter(tmp_POIs, !is.na(Type_WBOut))
-# nhd_out <- filter(nhd, COMID %in% wbout_COMIDs$COMID)
-#
-# # Create waterbody inlet POIs
-# wbin_COMIDs <- filter(nhd, WB == 0,
-# DnHydroseq %in% filter(nhd, WB == 1)$Hydroseq) %>%
-# select(-WBAREACOMI) %>%
-# switchDiv(., nhd) %>%
-# inner_join(st_drop_geometry(filter(nhd, minNet == 1)) %>%
-# select(Hydroseq, WBAREACOMI), by = c("DnHydroseq" = "Hydroseq")) %>%
-# select(COMID, WBAREACOMI, minNet) %>%
-# group_by(COMID) %>%
-# # Ensure the inlets are on the network defined by confluence POIs
-# filter(minNet == 1) %>%
-# select(COMID, WBAREACOMI)
-#
-# # Headwater Waterbodies that may need the network extended through the inlet
-# need_wbin <- st_drop_geometry(filter(WBs_VPU, source == "NHDv2WB")) %>%
-# dplyr::select(COMID)%>%
-# dplyr::filter(!COMID %in% wbin_COMIDs$WBAREACOMI)
-#
-# up_lp <- filter(nhd, LevelPathI %in% nhd_out$LevelPathI,
-# WBAREACOMI %in% need_wbin$COMID) %>%
-# group_by(WBAREACOMI) %>%
-# filter(n() > 1, StartFlag != 1) %>%
-# slice(which.min(Hydroseq))
-#
-# wbin_pois <- select(nhd, COMID, DnHydroseq, LevelPathI) %>%
-# inner_join(select(st_drop_geometry(up_lp), Hydroseq, WBAREACOMI, LevelPathI),
-# by = c("DnHydroseq" = "Hydroseq", "LevelPathI" = "LevelPathI")) %>%
-# select(COMID, WBAREACOMI) %>%
-# rbind(wbin_COMIDs)
-#
-# tmp_POIs <- POI_creation2(st_drop_geometry(wbin_pois), nhd, "WBIn") %>%
-# addType2(., tmp_POIs, "WBIn", nexus = FALSE, bind = TRUE)
-#
-# tmp_POIs <- tmp_POIs %>%
-# left_join(select(st_drop_geometry(sparrow_inlets), WB_COMID, nat_id),
-# by = c("Type_WBIn" = "WB_COMID"), multiple = "all") %>%
-# mutate(Type_Sparrow = ifelse(is.na(Type_Sparrow), nat_id, Type_Sparrow)) %>%
-# select(-nat_id)
-#
-# missing_wb_in <- filter(sparrow_inlets, !nat_id %in% tmp_POIs$Type_Sparrow |
-# !WB_COMID %in% tmp_POIs$Type_WBIn)
-#
-# if(nrow(missing_wb_in) > 0){write_sf(missing_wb_in, temp_gpkg, "missing_wbin")}
-# #tmp_POIs$nexus <- ifelse(is.na(tmp_POIs$nexus), FALSE, tmp_POIs$nexus)
-# write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-# } else {
-# tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-# nhd <- read_sf(nav_gpkg, nhd_flowline)
-# }
-#
-# mapview(filter(tmp_POIs, !is.na(Type_WBIn)), layer.name = "Waterbody inlet POIs", col.regions = "blue") +
-# mapview(filter(tmp_POIs, !is.na(Type_WBOut)), layer.name = "Waterbody outlet POIs", col.regions = "red")
-```
-
-
-This chunk breaks up potential aggregated segments where the elevation range of the are contributing to the segment exceeds 500-m
-```{r Elevation Break POIs}
-# # derive incremental segments from POIs
-# inc_segs <- make_incremental_segments(filter(nhd, minNet == 1),
-# filter(st_drop_geometry(nhd),
-# COMID %in% tmp_POIs$COMID, COMID %in% filter(nhd, minNet == 1)$COMID)) %>%
-# # bring over VAA data
-# inner_join(select(atts, COMID, DnHydroseq, VA_MA, TOTMA, LENGTHKM, MAXELEVSMO,
-# MINELEVSMO, WBAREACOMI, WBAreaType, FTYPE,
-# AreaSqKM, TotDASqKM), by = "COMID")
-#
-# # Build dataframe for creation of points based on breaks
-# elev_pois_split <- inc_segs %>%
-# group_by(POI_ID) %>%
-# # Get elevation info
-# mutate(MAXELEVSMO = na_if(MAXELEVSMO, -9998), MINELEVSMO = na_if(MINELEVSMO, -9998),
-# elev_diff_seg = max(MINELEVSMO) - min(MINELEVSMO)) %>%
-# # Filter out to those incremental segments that need splitting above the defined elevation threshold
-# filter((max(MINELEVSMO) - min(MINELEVSMO)) > (elev_diff * 100)) %>%
-# # Do not aggregate on fake lfowlines
-# filter(AreaSqKM > 0, TotDASqKM > max_elev_TT_DA) %>%
-# # Order from upstream to downstream
-# arrange(desc(Hydroseq)) %>%
-# # Get attributes used for splitting
-# # csum_length = cumulative length US to DS along segment, cumsum_elev = cumulative US to DS elev diff along segment
-# # inc_elev = elevation range of each segment
-# # iter = estimated number of times we'd need to split existing agg flowlines, or number of rows in set
-# mutate(csum_length = cumsum(LENGTHKM),
-# inc_elev = cumsum(MINELEVSMO - MINELEVSMO),
-# #nc_elev_diff = c(inc_elev[1], (inc_elev - lag(inc_elev))[-1]),
-# iter = min(round(max(inc_elev) / (elev_diff * 100)), n()),
-# elev_POI_ID = NA) %>%
-# # Remove segments we can't break
-# filter(iter > 0, n() > 1) %>%
-# # Track original iteration number
-# mutate(orig_iter = iter) %>%
-# ungroup()
-#
-# if(nrow(elev_pois_split) > 0) {
-#
-# # For reach iteration
-# elev_POIs <- do.call("rbind", lapply(unique(elev_pois_split$POI_ID), split_elev,
-# elev_pois_split, elev_diff*100, max_elev_TT_DA)) %>%
-# filter(!elev_POI_ID %in% tmp_POIs$COMID, COMID == elev_POI_ID) %>%
-# mutate(Type_Elev = 1) %>%
-# select(COMID, Type_Elev) %>%
-# unique()
-#
-# if(nrow(elev_POIs) > 0){
-# tmp_POIs <- POI_creation(select(elev_POIs, COMID, Type_Elev), nhd, "Elev") %>%
-# addType(., tmp_POIs, "Elev")
-# }
-#
-#
-# } else {
-#
-# tmp_POIs$Type_Elev <- rep(NA, nrow(tmp_POIs))
-#
-# write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-#
-# tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-# nhd <- read_sf(nav_gpkg, nhd_flowline)
-# }
-#
-# if(!all(is.na(tmp_POIs$Type_Elev)))
-# mapview(filter(tmp_POIs, Type_Elev == 1), layer.name = "Elevation break POIs", col.regions = "blue")
-```
-
-```{r Time of travel POIs}
-# if(all(is.na(tmp_POIs$Type_Travel))) {
-#
-# # derive incremental segments from POIs
-# inc_segs <- make_incremental_segments(filter(nhd, minNet == 1),
-# filter(st_drop_geometry(nhd),
-# COMID %in% tmp_POIs$COMID, COMID %in% filter(nhd, minNet == 1)$COMID)) %>%
-# # bring over VAA data
-# inner_join(select(atts, COMID, DnHydroseq, VA_MA, TOTMA, LENGTHKM, MAXELEVSMO, MINELEVSMO,
-# WBAREACOMI, WBAreaType, FTYPE,
-# AreaSqKM, TotDASqKM), by = "COMID")
-#
-# # TT POIs
-# tt_pois_split <- inc_segs %>%
-# # Should we substitute with a very small value
-# mutate(VA_MA = ifelse(VA_MA < 0, NA, VA_MA)) %>%
-# mutate(FL_tt_hrs = (LENGTHKM * ft_per_km)/ VA_MA / s_per_hr ) %>%
-# group_by(POI_ID) %>%
-# filter(sum(FL_tt_hrs) > travt_diff, TotDASqKM > max_elev_TT_DA) %>%
-# # Sort upstream to downsream
-# arrange(desc(Hydroseq)) %>%
-# # Get cumulative sums of median elevation
-# mutate(csum_length = cumsum(LENGTHKM),
-# cumsum_tt = cumsum(FL_tt_hrs),
-# iter = min(round(sum(FL_tt_hrs) / travt_diff), n()),
-# tt_POI_ID = NA) %>%
-# # Only look to split segments based on travel time longer than 10 km
-# #filter(sum(LENGTHKM) > (split_meters/1000)) %>%
-# filter(iter > 0, n() > 1) %>%
-# # Track original iteration number
-# mutate(orig_iter = iter) %>%
-# ungroup()
-#
-# # For reach iteration
-# tt_POIs <- do.call("rbind", lapply(unique(tt_pois_split$POI_ID), split_tt,
-# tt_pois_split, travt_diff, max_elev_TT_DA)) %>%
-# filter(!tt_POI_ID %in% tmp_POIs$COMID, COMID == tt_POI_ID) %>%
-# mutate(Type_Travel = 1) %>%
-# select(COMID, Type_Travel) %>%
-# unique()
-#
-# tmp_POIs <- POI_creation2(select(tt_POIs, COMID, Type_Travel), nhd, "Travel") %>%
-# addType2(., tmp_POIs, "Travel", nexus = FALSE, bind = TRUE)
-#
-# write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
-# }else{
-# tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
-# nhd <- read_sf(nav_gpkg, nhd_flowline)
-# }
-#
-# mapview(filter(tmp_POIs, Type_Travel == 1), layer.name = "Travel time break POIs", col.regions = "blue")
-```
-
-```{r Final POIs}
-# # Derive final POI set ----------------------
-# if(needs_layer(nav_gpkg, pois_all_layer)) {
-#
-# unCon_POIs <- filter(nhd, COMID %in% tmp_POIs$COMID, AreaSqKM == 0) %>%
-# distinct()
-#
-# # If any POIs happened to fall on flowlines w/o catchment
-# if (nrow(unCon_POIs) >0){
-# # For confluence POIs falling on Flowlines w/o catchments, derive upstream valid flowline,
-# poi_fix <- DS_poiFix2(tmp_POIs, nhd %>% distinct())
-# new_POIs <- st_compatibalize(poi_fix$new_POIs, tmp_POIs)
-# xWalk <- distinct(poi_fix$xWalk)
-# not_matched <- distinct(poi_fix$unmatched)
-#
-# # POIs that didn't need to be moved
-# tmp_POIs_fixed <- filter(tmp_POIs, nexus == TRUE |
-# !COMID %in% c(poi_fix$xWalk$oldPOI, poi_fix$xWalk$COMID))
-# # bind together
-# final_POIs <- bind_rows(tmp_POIs_fixed, select(poi_fix$new_POIs, -c(oldPOI, to_com)))
-#
-# # if a POI will be a nexus, it can not be terminal.
-# final_POIs <- mutate(final_POIs, Type_Term = ifelse(nexus == 1, NA, Type_Term))
-#
-# # Write out fixes
-# write_sf(new_POIs, temp_gpkg, poi_moved_layer)
-# write_sf(xWalk, temp_gpkg, poi_xwalk_layer)
-# # write out final POIs
-# write_sf(final_POIs, nav_gpkg, pois_all_layer)
-# write_sf(not_matched, temp_gpkg, "on_zeroDA")
-# } else {
-# # If no fixes designate as NA
-# xWalk <- NA
-#
-# tmp_POIs$nexus <- as.integer(tmp_POIs$nexus)
-#
-# # if a POI will be a nexus, it can not be terminal.
-# tmp_POIs <- mutate(tmp_POIs, Type_Term = ifelse(nexus == 1, NA, Type_Term))
-#
-# # write out final POIs
-# write_sf(tmp_POIs, temp_gpkg, pois_all_layer)
-# final_POIs <- tmp_POIs
-# }
-#
-# } else {
-# # Need all three sets for attribution below
-# final_POIs <- read_sf(nav_gpkg, pois_all_layer)
-# new_POIs <- if(!needs_layer(temp_gpkg, poi_moved_layer)) read_sf(temp_gpkg, poi_moved_layer) else (NA)
-# xWalk <- if(!needs_layer(temp_gpkg, poi_xwalk_layer)) read_sf(temp_gpkg, poi_xwalk_layer) else (NA)
-# unCon_POIs <- filter(st_drop_geometry(filter(nhd, minNet == 1)),
-# COMID %in% tmp_POIs$COMID, AreaSqKM == 0)
-# }
-```
-
-```{r Draft nsegment generation}
-# if(needs_layer(temp_gpkg, nsegments_layer)) {
-#
-# # Sort POIs by Levelpath and Hydrosequence in upstream to downstream order
-# seg_POIs <- filter(st_drop_geometry(nhd),
-# COMID %in% tmp_POIs$COMID, COMID %in% filter(nhd, minNet == 1)$COMID)
-# # derive incremental segments from POIs
-# inc_segs <- make_incremental_segments(filter(nhd, minNet == 1), seg_POIs)
-#
-# nhd_Final <- nhd %>%
-# left_join(select(inc_segs, COMID, POI_ID), by = "COMID")
-#
-# # create and write out final dissolved segments
-# nsegments_fin <- segment_creation(nhd_Final, xWalk)
-#
-# nhd_Final <- select(nhd_Final, -POI_ID) %>%
-# left_join(select(st_drop_geometry(nsegments_fin$raw_segs), COMID, POI_ID), by = "COMID")
-# nsegments <- nsegments_fin$diss_segs
-#
-# # Produce the minimal POIs needed to derive the network based on LP and terminal outlets
-# strucFeat <- structPOIsNet(nhd_Final, final_POIs)
-# nhd_Final <- nhd_Final %>%
-# mutate(struct_POI = ifelse(COMID %in% strucFeat$struc_POIs$COMID, 1, 0),
-# struct_Net = ifelse(COMID %in% strucFeat$structnet$COMID, 1, 0))
-#
-# write_sf(nhd_Final, nav_gpkg, nhd_flowline)
-# write_sf(nsegments, temp_gpkg, nsegments_layer)
-# } else {
-# # Read in NHDPlusV2 flowline simple features and filter by vector processing unit (VPU)
-# final_POIs <- read_sf(nav_gpkg, pois_all_layer)
-# nhd_Final <- read_sf(nav_gpkg, nhd_flowline)
-# nsegments <- read_sf(temp_gpkg, nsegments_layer)
-# }
-#
-# # Ensure that all the problem POIs have been taken care of
-# sub <- nhd_Final %>% filter(COMID %in% final_POIs$COMID)
-# print (paste0(nrow(sub[sub$AreaSqKM == 0,]), " POIs on flowlines with no local drainage contributions"))
-# noDA_pois <- filter(final_POIs, COMID %in% filter(sub, AreaSqKM == 0)$COMID)
-# write_sf(noDA_pois, temp_gpkg, "noDA_pois")
-```
-
-# Start from here next time
-```{r POI Collapse}
-# # number POIs
-# final_POIs <- mutate(final_POIs, id = row_number(), moved = NA)
-# write_sf(final_POIs, nav_gpkg, pois_all_layer)
-#
-# # Load data
-# if(collapse){
-#
-# # Move HUC12 to other POIs
-# huc_to_poi <- poi_move(final_POIs, "Type_HUC12", .10, 2.5)
-# move_table <- huc_to_poi$moved_points %>%
-# mutate(move_cat = "HUC12 to other")
-#
-# # Gages to confluences
-# gage_to_conf <- poi_move(huc_to_poi$final_pois, "Type_Gages", .05, 2.5, "Type_Conf") # 47
-# move_table <- gage_to_conf$moved_points %>%
-# mutate(move_cat = "Gage to conf") %>%
-# rbind(move_table)
-#
-# # Gages to waterbody inlets
-# gage_to_wbin <- poi_move(gage_to_conf$final_pois, "Type_Gages", .05, 2.5, "Type_WBIn") # 2
-# move_table <- gage_to_wbin$moved_points %>%
-# mutate(move_cat = "Gage to WBin") %>%
-# rbind(move_table)
-#
-# # Gages to waterbody outlets
-# gage_to_wbout <- poi_move(gage_to_wbin$final_pois, "Type_Gages", .05, 2.5, "Type_WBOut") # 6
-# move_table <- gage_to_wbout$moved_points %>%
-# mutate(move_cat = "Gage to WBout") %>%
-# rbind(move_table)
-#
-# # Streamgage to term outlet
-# gage_to_term <- poi_move(gage_to_wbout$final_pois, "Type_Gages", .05, 2.5, "Type_Term")
-# # NID to waterbody outlet
-# nid_to_wbout <- poi_move(gage_to_wbout$final_pois, "Type_NID", .025, 1, "Type_WBOut")
-#
-# # Sparrow to confluence
-# sparrow_to_conf <- poi_move(gage_to_wbout$final_pois, "Type_Sparrow", .05, 2.5, "Type_Conf")
-# move_table <- sparrow_to_conf$moved_points %>%
-# mutate(move_cat = "Sparrow to Conf") %>%
-# rbind(move_table)
-#
-# # Sparrow to WBoutlet
-# sparrow_to_wbout <- poi_move(sparrow_to_conf$final_pois, "Type_Sparrow", .05, 2.5, "Type_WBOut")
-# # Sparrow to WBInlet
-# sparrow_to_wbin <- poi_move(gage_to_wbout$final_pois, "Type_Sparrow", .05, 2.5, "Type_WBIn")
-#
-# POIs <- sparrow_to_conf$final_pois
-# write_sf(POIs, nav_gpkg, final_poi_layer)
-# write_sf(move_table, temp_gpkg, "pois_collapsed")
-# } else {
-# write_sf(final_POIs, nav_gpkg, final_poi_layer)
-# POIs <- final_POIs
-# }
-```
-
-# ```{r lookup}
-# # Final POI layer
-# POIs <- POIs %>%
-# arrange(COMID) %>%
-# mutate(identifier = row_number())
-#
-# # Unique POI geometry
-# final_POI_geom <- POIs %>%
-# select(identifier) %>%
-# cbind(st_coordinates(.)) %>%
-# group_by(X, Y) %>%
-# mutate(geom_id = cur_group_id()) %>%
-# ungroup()
-#
-# final_POIs_table <- POIs %>%
-# inner_join(select(st_drop_geometry(final_POI_geom), -X, -Y), by = "identifier") %>%
-# select(-identifier)
-#
-# # POI data theme table
-# pois_data_table <- reshape2::melt(st_drop_geometry(select(final_POIs_table,
-# -nexus)),
-# id.vars = c("COMID", "geom_id")) %>%
-# filter(!is.na(value)) %>%
-# group_by(COMID, geom_id) %>%
-# mutate(identifier = cur_group_id()) %>%
-# rename(hy_id = COMID, poi_id = identifier, hl_reference = variable, hl_link = value) %>%
-# distinct()
-#
-# # POI Geometry table
-# poi_geometry_table <- select(final_POIs_table, hy_id = COMID, geom_id) %>%
-# inner_join(distinct(pois_data_table, hy_id, geom_id, poi_id),
-# by = c("hy_id" = "hy_id", "geom_id" = "geom_id")) %>%
-# distinct()
-#
-# write_sf(pois_data_table, nav_gpkg, "poi_data")
-# write_sf(poi_geometry_table, nav_gpkg, "poi_geometry")
-#
-# # Create xwalk table
-# if(needs_layer(nav_gpkg, lookup_table_refactor)) {
-# full_cats <- readRDS(data_paths$fullcats_table)
-#
-# lookup <- dplyr::select(sf::st_drop_geometry(nhd_Final),
-# NHDPlusV2_COMID = COMID,
-# realized_catchmentID = COMID,
-# mainstem = LevelPathI) %>%
-# dplyr::mutate(realized_catchmentID = ifelse(realized_catchmentID %in% full_cats$FEATUREID,
-# realized_catchmentID, NA)) %>%
-# left_join(select(st_drop_geometry(poi_geometry_table), hy_id, poi_geom_id = geom_id),
-# by = c("NHDPlusV2_COMID" = "hy_id"))
-#
-# sf::write_sf(lookup, nav_gpkg, lookup_table_refactor)
-# }
-#
-# ```
diff --git a/workspace/03_hyRefactor_flines.Rmd b/workspace/03_hyRefactor_flines.Rmd
deleted file mode 100644
index cabb6c3..0000000
--- a/workspace/03_hyRefactor_flines.Rmd
+++ /dev/null
@@ -1,198 +0,0 @@
----
-title: "NHD Navigate"
-output: html_document
-editor_options:
- chunk_output_type: console
----
-
-Project: GFv2.0
-Script purpose: Run hyRefactor workflow on flowline network
-Date: 2/29/2019
-Author: [Dave Blodgett](dblodgett@usgs.gov)
-Notes:
-
-```{r setup_0, echo=FALSE, cache=FALSE}
-knitr::opts_chunk$set(
- collapse = TRUE,
- comment = "#>",
- fig.width=6,
- fig.height=4,
- cache=FALSE)
-```
-
-Load functions and configuration.
-See `hyRefactory_config.R` for all constants.
-```{r setup}
-source("R/utils.R")
-source("R/config.R")
-source("R/user_vars.R")
-```
-
-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_data <- read_sf(nav_gpkg, poi_data_table)
-
-POIs <- read_sf(nav_gpkg, poi_geometry_table)
-
-POI_id <- st_drop_geometry(POIs)
-
-events <- read_sf(nav_gpkg, event_geometry_table)
-```
-
-Read POIs and add nhd outlets. Save to a non-spatial table.
-```{r refactor}
-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_ref$DnHydroseq, AreaSqKM > 0)
-
-# build final outlets set
-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)
-
-# Need to avoid refactoring catchments that are long and thing and
-# reasonably large. They cause problems in a downstream step.
-avoid <- dplyr::filter(nhd, (sqrt(AreaSqKM) / LENGTHKM) > 3 & AreaSqKM > 1)
-
-# attribute flowline network used for refactor
-nhd <- mutate(nhd, refactor = ifelse(TerminalPa %in% TerminalPaths, 1, 0))
-write_sf(nhd, out_refac_gpkg, nhd_flowline)
-
-# Prep flowlines for refactor
-nhdplus_flines <- sf::st_zm(filter(nhd, refactor == 1)) %>%
- st_as_sf()
-
-if(!needs_layer(nav_gpkg, event_geometry_table)){
- events <- read_sf(nav_gpkg, event_geometry_table)
-
- events_ref <- filter(events, hy_id %in% nhdplus_flines$COMID) %>%
- rename(COMID = hy_id) %>%
- distinct(COMID, REACHCODE, REACH_meas,
- poi_id, geom) %>%
- arrange(poi_id) %>%
- mutate(event_identifier = row_number())
-
- outlets <- rename(outlets, COMID = hy_id) %>%
- filter(!poi_id %in% events_ref$poi_id)
-} else {
- events_ref <- NULL
-}
-
-if(needs_layer(out_refac_gpkg, refactored_layer)) {
-
- tf <- paste0("temp/refactored_", rpu_code,".gpkg")
- tr <- paste0("temp/reconciled_", rpu_code, ".gpkg")
-
- unlink(list(tf, tr))
- 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,
- collapse_flines_main_meters = combine_meters,
- out_refactored = tf,
- out_reconciled = tr,
- three_pass = TRUE,
- purge_non_dendritic = FALSE,
- events = events_ref,
- exclude_cats = c(outlets$COMID,
- avoid$COMID, POI_downstream$COMID),
- warn = TRUE)
-
- write_sf(st_transform(read_sf(tf), proj_crs), out_refac_gpkg, refactored_layer)
- write_sf(st_transform(read_sf(tr), proj_crs), out_refac_gpkg, reconciled_layer)
-
- unlink(list(tf, tr))
- rm(tf, tr)
-}
-```
-
-```{r lookup}
-# create lookup for ref flowlines to use in the non-dendritic steps
-refactor_lookup <- st_drop_geometry(read_sf(out_refac_gpkg, reconciled_layer)) %>%
- dplyr::select(ID, member_COMID) %>%
- dplyr::mutate(member_COMID = strsplit(member_COMID, ",")) %>%
- tidyr::unnest(cols = member_COMID) %>%
- dplyr::mutate(NHDPlusV2_COMID = as.integer(member_COMID)) %>% # note as.integer truncates
- dplyr::rename(reconciled_ID = ID)
-
-if(is.character(refactor_lookup$reconciled_ID)) refactor_lookup$reconciled_ID <- as.integer(refactor_lookup$reconciled_ID)
-
-lookup_table <- tibble::tibble(NHDPlusV2_COMID = unique(as.integer(refactor_lookup$member_COMID))) %>%
- dplyr::left_join(refactor_lookup, by = "NHDPlusV2_COMID")
-
-readr::write_csv(lookup_table, lookup_table_file)
-write_sf(lookup_table, out_refac_gpkg, lookup_table_refactor)
-
-if(needs_layer(out_refac_gpkg, outlets_layer)) {
-
- # Join refactored to original NHD
- refactored <- read_sf(out_refac_gpkg, refactored_layer) %>%
- select(member_COMID = COMID, Hydroseq, event_identifier, event_REACHCODE) %>%
- inner_join(select(st_drop_geometry(nhd), orig_COMID = COMID, Hydroseq), by = "Hydroseq")
-
- if(nrow(events_ref) > 0){
- # Subset for events
- refactored_events <- refactored %>%
- filter(!is.na(event_REACHCODE), !is.na(event_identifier)) %>%
- mutate(event_identifier = as.numeric(event_identifier))
-
- event_outlets <- events_ref %>%
- inner_join(st_drop_geometry(refactored_events), by = "event_identifier") %>%
- select(hy_id = COMID, event_id = event_identifier, poi_id, member_COMID)%>%
- inner_join(select(lookup_table, -NHDPlusV2_COMID), by = "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 %>%
- filter(!poi_id %in% event_outlets$poi_id) %>%
- inner_join(lookup_table,
- by = c("COMID" = "NHDPlusV2_COMID")) %>%
- group_by(COMID) %>%
- filter(member_COMID == max(member_COMID)) %>%
- select(hy_id = COMID, poi_id, member_COMID, reconciled_ID, 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 %>%
- inner_join(lookup_table,
- by = c("COMID" = "NHDPlusV2_COMID")) %>%
- group_by(COMID) %>%
- filter(member_COMID == max(member_COMID)) %>%
- select(hy_id = COMID, poi_id, member_COMID, reconciled_ID, type)
- }
-
- write_sf(outlets_ref_COMID, out_refac_gpkg, outlets_layer)
-} else {
- outlets_ref_COMID <- read_sf(out_refac_gpkg, outlets_layer)
-}
-
-check_dups_poi <- outlets_ref_COMID %>%
- group_by(reconciled_ID) %>%
- filter(n() > 1) %>%
- ungroup()
-
-if(nrow(check_dups_poi) > 1){
- print("Double-check for double POIs")
- write_sf(check_dups_poi, out_refac_gpkg, paste0(dup_pois, "_", rpu_code))
-} else {
- print("no double POIs detected")
-}
-```
-
diff --git a/workspace/04-1_hyRefactor_cats.Rmd b/workspace/04-1_hyRefactor_cats.Rmd
deleted file mode 100644
index e920d7a..0000000
--- a/workspace/04-1_hyRefactor_cats.Rmd
+++ /dev/null
@@ -1,76 +0,0 @@
----
-title: "Reconcile Catchments"
-output: html_document
----
-
-Project: GFv2.0
-Script purpose: Run hyRefactor workflow on catchments
-Date: 2/29/2019
-Author: [Dave Blodgett](dblodgett@usgs.gov)
-Notes: Assumes that hyRefactor_flines.Rmd has previously been run.
-
-```{r setup_0, echo=FALSE, cache=FALSE}
-knitr::opts_chunk$set(
- collapse = TRUE,
- comment = "#>",
- fig.width=10,
- fig.height=8)
-```
-
-Load functions and configuration.
-See `config.R` for all constants.
-```{r libs}
-source("R/utils.R")
-source("R/config.R")
-source("R/user_vars_ww.R")
-source("R/NHD_navigate.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.
-```{r reconcile_catchments}
-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)
- # refactored has the original COMIDs and other attributes
- refactored <- st_transform(read_sf(out_refac_gpkg, refactored_layer), raster_crs) %>%
- distinct(.keep_all = T)
-
- cats <- st_transform(read_sf(out_refac_gpkg, nhd_catchment), raster_crs)
-
- divides <- reconcile_catchment_divides(catchment = cats,
- fline_ref = refactored,
- fline_rec = reconciled,
- fdr = fdr,
- fac = fac,
- para = para_reconcile,
- cache = cache_split,
- keep = NULL)
-
- if(exists("divides")) {
-
- load(cache_split)
-
- 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, proj_crs), out_refac_gpkg, split_divide_layer)
-
- rm(split_cats)
- }
-
- unlink(cache_split)
- write_sf(sf::st_transform(divides, proj_crs), out_refac_gpkg, divide_layer)
-}
-
-sf::st_layers(out_refac_gpkg)
-```
\ No newline at end of file
diff --git a/workspace/04-2_aggregate_cats.Rmd b/workspace/04-2_aggregate_cats.Rmd
deleted file mode 100644
index bb49aa1..0000000
--- a/workspace/04-2_aggregate_cats.Rmd
+++ /dev/null
@@ -1,222 +0,0 @@
----
-title: "hyRefactor Catchments"
-output: html_document
----
-
-Project: GFv2.0
-Script purpose: Run hyRefactor workflow on catchments
-Date: 2/29/2019
-Author: [Dave Blodgett](dblodgett@usgs.gov)
-Notes: Assumes that hyRefactor_flines.Rmd has previously been run.
-
-```{r setup_0, cache=FALSE, echo = FALSE}
-knitr::opts_chunk$set(
- collapse = TRUE,
- comment = "#>",
- fig.width=10,
- fig.height=8)
-```
-
-Load functions and configuration.
-See `hyRefactory_config.R` for all constants.
-```{r setup}
-source("R/utils.R")
-source("R/config.R")
-source("R/user_vars.R")
-source("R/NHD_navigate.R")
-source("R/user_vars.R")
-```
-
-Derive set of outlets to aggregate catchments for. There are three sets of
-locations loaded below: POIs, outlets, and events. These are reconciled
-into a single set of HRU outlets and all identifiers mapped
-to the id scheme used by the refactored network. "mapped_outlets_all" is the
-result of this block.
-```{r POIs/outlets}
-# Load nhd and outlets
-nhd <- read_sf(out_refac_gpkg, nhd_flowline)
-
-# Load outlets
-outlets_POI <- read_sf(out_refac_gpkg, outlets_layer) %>%
- 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), proj_crs)
-refactored <- st_drop_geometry(read_sf(out_refac_gpkg, refactored_layer)) %>%
- select(COMID, Hydroseq) %>%
- mutate(COMID = as.integer(COMID))
-
-divides <- st_transform(read_sf(out_refac_gpkg, divide_layer), proj_crs)
-
-lookup <- read_sf(out_refac_gpkg, lookup_table_refactor) %>%
- inner_join(select(st_drop_geometry(nhd), COMID, TotDASqKM), by = c("NHDPlusV2_COMID" = "COMID"))
-
-highest_DA <- lookup %>%
- inner_join(refactored, by = c("NHDPlusV2_COMID" = "COMID")) %>%
- group_by(reconciled_ID) %>%
- filter(Hydroseq == min(Hydroseq))
-```
-
-```{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_cast(sf::st_geometry(divides[sf::st_is_empty(divides),]), "MULTIPOLYGON")
-
- # 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))
-
- # 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")
-
- }
-
- # Identify reconciled flowpaths that are terminal and below DA threshold
- reconciled_DA <- filter(reconciled, TotDASqKM <= aggregate_da_thresh_sqkm, is.na(toID))
-
- # Ones not in existing pois
- extra_outlets <- filter(reconciled_DA, !ID %in% outlets$ID) %>%
- pull(ID)
-
- 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)
-
- outlets <- outlets %>%
- left_join(select(st_drop_geometry(reconciled), ID, toID), by = "ID") %>%
- mutate(type = ifelse(is.na(toID), "terminal", "outlet")) %>%
- mutate(type = ifelse(type == "terminal" & !is.na(toID), "outlet", type)) %>%
- select(-toID)
-
- outlets <- filter(outlets, ID %in% reconciled$ID) %>%
- distinct()
-
- agg_cats <- aggregate_catchments(flowpath = reconciled,
- divide = divides,
- outlets = outlets,
- 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)
-
- # Get physical geometry of reconciled FL end node whose ID is an aggregated outlet
- rec_outlets <- filter(st_drop_geometry(reconciled), ID %in% agg_cats$fline_sets$ID) %>%
- cbind(get_node(filter(reconciled, ID %in% agg_cats$fline_sets$ID))) %>%
- st_as_sf()
-
- # Build final POI set
- POIs_att <- select(rec_outlets, ID, toID) %>%
- # Bring over POI information
- 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()
-
-
- if(nrow(POIs_missing) > 0){
- final_POIs <- data.table::rbindlist(list(POIs_att, select(POIs_missing,
- -c(Type_Term, Type_Con))), fill = TRUE) %>%
- st_as_sf()
- } else {
- final_POIs <- POIs_att
- }
-
-
- 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))
-
- final_POIs <- read_sf(out_agg_gpkg, mapped_outlets_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, ",")) %>%
- tidyr::unnest(cols = member_COMID) %>%
- dplyr::mutate(NHDPlusV2_COMID = as.integer(member_COMID)) %>% # note as.integer truncates
- dplyr::rename(reconciled_ID = ID)
-
-aggregate_lookup_fline <- dplyr::select(st_drop_geometry(agg_cats$fline_sets), ID, set) %>%
- dplyr::mutate(set = strsplit(set, ",")) %>%
- tidyr::unnest(cols = set) %>%
- dplyr::rename(aggregated_flowpath_ID = ID, reconciled_ID = set)
-
-aggregate_lookup_catchment <- dplyr::select(st_drop_geometry(agg_cats$cat_sets), ID, set) %>%
- dplyr::mutate(set = strsplit(set, ",")) %>%
- tidyr::unnest(cols = set) %>%
- dplyr::rename(aggregated_divide_ID = ID, reconciled_ID = set)
-
-if(is.character(aggregate_lookup_catchment$reconciled_ID)) aggregate_lookup_catchment$reconciled_ID <- as.integer(aggregate_lookup_catchment$reconciled_ID)
-if(is.character(aggregate_lookup_fline$reconciled_ID)) aggregate_lookup_fline$reconciled_ID <- as.integer(aggregate_lookup_fline$reconciled_ID)
-
-lookup_table <- tibble::tibble(NHDPlusV2_COMID = unique(as.integer(refactor_lookup$member_COMID))) %>%
- dplyr::left_join(refactor_lookup, by = "NHDPlusV2_COMID") %>%
- dplyr::left_join(aggregate_lookup_fline, by = "reconciled_ID") %>%
- dplyr::left_join(aggregate_lookup_catchment, by = "reconciled_ID")
-
-readr::write_csv(lookup_table, lookup_table_file)
-write_sf(lookup_table, out_agg_gpkg, lookup_table_refactor)
-
-sf::st_layers(out_agg_gpkg)
-
-mapview::mapview(list(agg_cats$cat_sets, agg_cats$fline_sets, final_POIs))
-```
diff --git a/workspace/04_merge.Rmd b/workspace/04_merge.Rmd
deleted file mode 100644
index 1ca2783..0000000
--- a/workspace/04_merge.Rmd
+++ /dev/null
@@ -1,15 +0,0 @@
-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-1_merge.Rmd b/workspace/05-1_merge.Rmd
deleted file mode 100644
index 5a82023..0000000
--- a/workspace/05-1_merge.Rmd
+++ /dev/null
@@ -1,94 +0,0 @@
----
-title: "07-1_merge.Rmd"
-output: html_document
-editor_options:
- chunk_output_type: console
----
-title: 07-1_merge.Rmd
-Project: GFv2.0
-Date: 4/2022
-Author: [David Blodgett](dblodgett@usgs.gov)
-Script purpose: Merge refactored outputs into VPU basis.
-
-
-```{r rmd, echo=F, message=F, warning=FALSE, fig.width=4}
-knitr::opts_chunk$set(
- collapse = TRUE,
- comment = "#>",
- fig.width=10,
- fig.height=8)
-```
-
-```{r setup}
-# Load custom functions
-source("R/utils.R")
-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 refactor}
-if(needs_layer(rfc_gpkg, divide_layer)) {
-
- # Thematic POIs
- POIs <- read_sf(nav_gpkg, poi_data_table) %>%
- select(-geom_id)
-
- 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,
- split_divide_layer,
- agg_fline_layer,
- agg_cats_layer,
- mapped_outlets_layer)
-
- sf::write_sf(select(merged_layers[[lookup_table_refactor]],
- -aggregated_flowpath_ID,
- -aggregated_divide_ID),
- rfc_gpkg, lookup_table_refactor)
-
- sf::write_sf(merged_layers[[reconciled_layer]],
- rfc_gpkg, reconciled_layer)
-
- sf::write_sf(merged_layers[[divide_layer]],
- rfc_gpkg, divide_layer)
-
- sf::write_sf(merged_layers[[split_divide_layer]],
- rfc_gpkg, split_divide_layer)
-
- sf::write_sf(merged_layers[[mapped_outlets_layer]], gf_gpkg, mapped_outlets_layer)
-
- merged_layers[[agg_cats_layer]] <- merged_layers[[agg_cats_layer]] %>%
- mutate(areasqkm = as.numeric(units::set_units(sf::st_area(.), "km^2")))
-
- sf::write_sf(merged_layers[[agg_cats_layer]], gf_gpkg, agg_cats_layer)
-
- merged_layers[[agg_fline_layer]] <- merged_layers[[agg_fline_layer]] %>%
- mutate(lengthkm = as.numeric(units::set_units(sf::st_length(.), "km"))) %>%
- left_join(select(sf::st_drop_geometry(merged_layers[[agg_cats_layer]]), ID, areasqkm), by = "ID")
-
- sf::write_sf(merged_layers[[agg_fline_layer]], gf_gpkg, agg_fline_layer)
-
- sf::write_sf(merged_layers[[lookup_table_refactor]],
- gf_gpkg, lookup_table_refactor)
-
- st_drop_geometry(merged_layers[[reconciled_layer]]) %>%
- select(id = ID, toID = toID, lengthkm = LENGTHKM, levelpathid = LevelPathID) %>%
- left_join(select(st_drop_geometry(merged_layers[[divide_layer]]),
- id = ID, areasqkm), by = "id") %>%
- relocate(id, toID, lengthkm, areasqkm, levelpathid) %>%
- write_sf(rfc_gpkg, catchment_network_table)
-
- st_drop_geometry(merged_layers[[agg_fline_layer]]) %>%
- select(id = ID, toid = toID, lengthkm, areasqkm) %>%
- left_join(distinct(select(merged_layers[[lookup_table_refactor]],
- id = aggregated_flowpath_ID,
- levelpathid = LevelPathID))) %>%
- sf::write_sf(gf_gpkg, catchment_network_table)
-}
-```
diff --git a/workspace/05-2_NonDend.Rmd b/workspace/05-2_NonDend.Rmd
deleted file mode 100644
index e700887..0000000
--- a/workspace/05-2_NonDend.Rmd
+++ /dev/null
@@ -1,412 +0,0 @@
----
-title: "07-2_NonDend.Rmd"
-output: html_document
-editor_options:
- chunk_output_type: console
----
-title: NonDend.Rmd
-Project: GFv2.0
-Date: 7/2020
-Author: [Marilyn Santiago, Andrew Bock](msant@usgs.gov, abock@usgs.govb)
-Script purpose:: To provide workflow for aggregating non-dendritic catchments
-Coastal zero-order and frontal catchments are accounted for first, followed by
-endorheic basins and sinks.
-
-
-```{r rmd, echo=F, message=F, warning=FALSE, fig.width=4}
-knitr::opts_chunk$set(
- collapse = TRUE,
- comment = "#>",
- fig.width=10,
- fig.height=8)
-```
-
-```{r setup}
-# set to true to output extra layers.
-debug <- FALSE
-
-# Load custom functions
-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"){
- full_nhd <- readRDS(file.path(data_paths$islands_dir,
- "NHDPlusNationalData/nhdplus_flowline.rds")) %>%
- st_transform(crs)
-} else {
- full_nhd <- readRDS(data_paths$VAA_rds)
-}
-
-# Elevation Data
-elev <- data_paths$elev_cm[grepl(paste0("Ned", substr(vpu_codes, 1, 2)),
- data_paths$elev_cm, ignore.case = TRUE)]
-
-if (vpu_codes == "08"){
- elev$rpu_03g <- data_paths$elev_cm$rpu_03g
-}
-
-# HUC extraction for specific NHD+ vpus
-if(vpu_codes == "02"){
- grep_exp <-"^02|^04"
-} else if (vpu_codes == "08") {
- grep_exp <- "^03|^08"
-} else {
- grep_exp <- paste0("^", substr(vpu_codes, start = 1, stop = 2))
- elev <- append(elev, list(rpu_03g = data_paths$elev_cm$rpu_03g))
-}
-
-cat_rpu_table <- readRDS(data_paths$fullcats_table)
-
-# get full flowline set including coastlines
-vpu_nhd <- full_nhd %>%
- filter(grepl(paste0("^", grep_exp, ".*"), .data$VPUID)) %>%
- align_nhdplus_names(.)
-
-vpu_WBD <- readRDS(file.path(data_paths$nhdplus_dir, "HUC12.rds")) %>%
- filter(grepl(paste0("^", grep_exp, ".*"), .data$HUC_12))
-
-nhd <- st_transform(read_sf(nav_gpkg, nhd_flowline), proj_crs)
-cats <- st_transform(read_sf(nav_gpkg, nhd_catchment), proj_crs)
-divides <- st_transform(read_sf(rfc_gpkg, divide_layer), proj_crs)
-lookup <- read_sf(gf_gpkg, lookup_table_refactor)
-
-```
-
-This next section intersects HUC12 for a given VPU with NHD+ catchments and catchment divides.
-
-See __get_HUC12_Xwalk__ in __R/utils.R__
-```{r Load HUC12 / xWalk information}
-if(needs_layer(ND_gpkg, xwalk_layer)){
-
- # Read in full NHD cats
- cats <- filter(cats, full_cats == 1)
- cats <- sf::st_make_valid(cats)
- divides <- sf::st_make_valid(divides)
-
- # Intersect NHD catchments and divides with HUC12
- nhd_wbd_int <- get_HUC12_Xwalk(vpu_codes, cats, divides,
- file.path(data_paths$nhdplus_dir,
- "CrosswalkTable_NHDplus_HU12.csv"),
- file.path(data_paths$nhdplus_dir, "HUC12.rds"))
-
- # Bring over divides/HUC12 intersection information into divides layer
- xwalk_divides_wbd <- st_drop_geometry(nhd_wbd_int$divides_HUC12) %>%
- select(-c(ACRES, HU_12_MOD))
-
- divides <- divides %>%
- left_join(distinct(xwalk_divides_wbd) %>%
- group_by(ID) %>%
- filter(intArea == max(intArea)) %>%
- ungroup() %>%
- select(ID, HUC_12_int, intArea, AreaHUC12),
- by = "ID")
-
- # Bring over divides/HUC12 intersection information into divides layer
- xwalk_nhd_wbd <- st_drop_geometry(nhd_wbd_int$cats_HUC12) %>%
- select(-c(ACRES, HU_12_MOD))
-
- rm(nhd_wbd_int)
-
- cats <- cats %>%
- left_join(distinct(xwalk_nhd_wbd) %>%
- group_by(FEATUREID) %>%
- slice(which.max(intArea)) %>%
- ungroup() %>%
- select(FEATUREID, HUC_12_int, intArea, AreaHUC12),
- by = "FEATUREID") %>%
- st_transform(proj_crs)
-
- # All intersecting NHD cat / divides / HUC12 geometries (1: many)
- write_sf(xwalk_nhd_wbd, ND_gpkg, xwalk_layer)
- write_sf(xwalk_divides_wbd, ND_gpkg, divides_xwalk)
-
- # Cats with the HUC_12 value (HUC_12_int) that is the max overlap
- write_sf(cats, ND_gpkg, nhd_catchment)
- write_sf(divides, ND_gpkg, divide_layer)
-} else {
- xwalk_nhd_wbd <- read_sf(ND_gpkg, xwalk_layer)
- xwalk_divides_wbd <- read_sf(ND_gpkg, divides_xwalk)
- cats <- read_sf(ND_gpkg, nhd_catchment)
- divides <- read_sf(ND_gpkg, divide_layer)
-}
-```
-
-This next step takes catchment divides and:
-
-1. Adds a field that indicates how a divide without an aggregated catchment ID assigned in hyRefactor::aggregate_catchments was aggregated (*__aggStep__*)
-2. Identifies internal sinks and non-dendritic areas entirely surrounded by a single aggregated catchment, and assigns those the catchment ID (see __dissolve_holes__ in __non_dend.R__ function)
-
-```{r First-cut}
-if(needs_layer(ND_gpkg, divides_nd)){
-
- # Bind divide with lookup, identify aggregation step
- divides_lu <- divides %>%
- left_join(distinct(select(lookup, reconciled_ID, aggregated_divide_ID)),
- by = c("ID" = "reconciled_ID")) %>%
- filter(!is.na(ID) & !is.na(aggregated_divide_ID)) %>%
- rename(POI_ID = aggregated_divide_ID) %>%
- # attribute that hyrefactor was the step used to aggregate these catchments
- mutate(aggStep = "hyRef")
-
- # Identify cats witin the full catchment set not refactored or aggregated, add
- # to divides data frame
- cats_miss <- cats %>%
- left_join(select(lookup, NHDPlusV2_COMID, member_COMID, reconciled_ID ,
- POI_ID = aggregated_divide_ID),
- by = c("FEATUREID" = "NHDPlusV2_COMID")) %>%
- filter(is.na(POI_ID), !member_COMID %in% divides_lu$member_COMID) %>%
- mutate(aggStep = NA) %>%
- distinct() %>%
- left_join(select(cat_rpu_table, FEATUREID, rpu = RPUID)) %>%
- #filter(!reconciled_ID %in% divides_lu$ID) %>%
- select(-member_COMID) %>%
- mutate(areasqkm = as.numeric(st_area(.)/1000000)) %>%
- select(ID = reconciled_ID, member_COMID = FEATUREID, rpu, HUC_12_int, intArea, AreaHUC12, POI_ID, aggStep)
-
- divides_lu <- divides_lu %>%
- select(-areasqkm) %>%
- rbind(cats_miss) %>%
- # resolve terminal non-POI outlets where we can
- nhdplusTools:::check_valid() %>%
- st_cast("POLYGON")
-
- # clean_geometry(divides_lu, ID = "POI_ID", keep = 1.0, crs = st_crs(divides_lu))
- divides_lu <- dissolve_holes(divides_lu)
- rm(cats_miss)
-
- # DEBUG:
- # HRU layer
- if(debug) {
- protoHRUs <- divides_lu %>%
- group_by("POI_ID") %>%
- summarize(do_union = TRUE) %>%
- sfheaders::sf_remove_holes(.) %>%
- st_make_valid()
- write_sf(protoHRUs, ND_gpkg, HRU_layer)
- rm(protoHRUs)
- }
-
- write_sf(divides_lu, ND_gpkg, divides_nd)
-} else {
- divides_lu <- read_sf(ND_gpkg, divides_nd)
-}
-```
-
-This next section takes divides not assigned an aggregated ID (*__POI_ID__*) and determines if they can be grouped on the basis of their HUC12/HUC10 ID. See __assign_HUC10/12__ in __non_dend.R__. This assigns a temporary aggregated ID that reflects the HUC10/12 the divides overlap with.
-
-```{r HUC12 cats}
-# Check for HUC12 aggregation assignments
-if(!"HUC_12" %in% unique(divides_lu$aggStep)){
- print(paste0("Currently there are ", sum(is.na(divides_lu$POI_ID)), " cats without a POI_ID"))
-
- # This will crate aggregated IDs for NHD catchments that match HUC12 footprint if aggregated together
- # The Coefficient of Areal correpsondence - is the last argument
- if("AREASQKM" %in% names(xwalk_nhd_wbd)) {
- xwalk_nhd_wbd <- rename(xwalk_nhd_wbd,
- AreaSqKM = AREASQKM)
- }
-
- divides_lu <- assign_HUC12(divides_lu, xwalk_nhd_wbd,
- filter(vpu_nhd, FTYPE != "Coastline"), CAC_thresh)
-
- # Determine if any unaggregated catchments match HUC10 footprint (mostly in the west)
- divides_lu <- assign_HUC10(divides_lu, xwalk_nhd_wbd,
- filter(vpu_nhd, FTYPE != "Coastline"), CAC_thresh) %>%
- select(-comid_char)
-
- if(debug) {
- # DEBUG:
- # HRU layer
- protoHRUs <- divides_lu %>%
- group_by(POI_ID) %>%
- summarize(do_union = TRUE) %>%
- sfheaders::sf_remove_holes(.) %>%
- st_make_valid()
- write_sf(protoHRUs, ND_gpkg, HRU_layer)
- rm(protoHRUs)
- }
-
- # Update missing_cats
- write_sf(divides_lu, ND_gpkg, divides_nd)
-} else {
- divides_lu <- read_sf(ND_gpkg, divides_nd)
-}
-```
-
-Deal with unaggregated catchment divides that are coastal if they exist within VPU
-
-Note these are primarily assigned the HUC12 as the temporary aggregated ID, but they differ from above in that they don't meet the Areal Correspondence threshold.
-
-Note that in this case it can involve ugly multi-part HRUs which we may want to split to single part.
-
-```{r Coastal Cats}
-# aggregate frontal hucs of same type
-if("Coastline" %in% unique(vpu_nhd$FTYPE)){
- print(paste0("Currently there are ", sum(is.na(divides_lu$POI_ID)), " cats without a POI_ID"))
-
- # Only run through if coastal catchments are present
- if(!"coastal" %in% unique(divides_lu$aggStep)){
-
- # Function to create coastal catchments by HUC12 and catch inflows.
- divides_lu <- coastal_cats(divides_lu, vpu_nhd, vpu_WBD)
-
- if(debug) {
- # DEBUG:
- # HRU layer
- protoHRUs <- divides_lu %>%
- group_by(POI_ID) %>%
- summarize(do_union = TRUE) %>%
- sfheaders::sf_remove_holes(.) %>%
- st_make_valid()
- # Write out updates
- write_sf(protoHRUs, ND_gpkg, HRU_layer)
- rm(protoHRUs)
- }
-
-
- # Update missing_cats
- write_sf(divides_lu, ND_gpkg, divides_nd)
- } else {
- divides_lu <- read_sf(ND_gpkg, divides_nd)
- }
-} else {
- print ("No Coastlines in VPU")
-}
-vpu_nhd <- filter(vpu_nhd, FTYPE != "Coastline")
-```
-
-
-This is where alot of the more complex steps happen. For most VPUs, there are not of catchments left to aggregate after this step.
-
-1, Determine where network terminals below the drainage area threshold used to derive the hyRefactor outlets contribute to (__assignable_cats__ in __non_dend.R__)
-2. Assign other cats by the basis of HUC12 and HUC12 of bordering catchments (__assignable_cats__)
-3. Uses the FAC grid to assign the rest ( __assign_holes__ )
-```{r FixingtheHoles}
-if(needs_layer(ND_gpkg, missing_cats)){
- print(paste0("Currently there are ", sum(is.na(divides_lu$POI_ID)), " cats without a POI_ID"))
-
- if(sum(is.na(divides_lu$POI_ID)) > 0){
-
- # Handle nhd sinks
- # arrange divides_lu by member_COMID and populate with RowID for easy iteration within next two steps
- divides_lu <- divides_lu %>%
- arrange(member_COMID, POI_ID) %>%
- mutate(row_ID = row_number()) %>%
- arrange(row_ID)
-
- HUC_sinks <- NHD_sinks(divides_lu, area_thresh = min_da_km_terminal/2,
- HUC12_table =file.path(data_paths$nhdplus_dir, "HUC12.rds"),
- NHD_sinks = read_sf(data_paths$nhdplus_gdb, "Sink"))
-
- if(length(HUC_sinks) == 2){
- divides_lu <- HUC_sinks$divides_poi
- sinks_table <- HUC_sinks$sink_cats_table
- }
-
- # Scan for terminals that may have been refactored
- missing_ds <- filter(divides_lu, is.na(POI_ID))
- term_refactored <- lapply(missing_ds$member_COMID, function(x){
- ds_ref <- get_DM(vpu_nhd,x, include = F)
- if(length(ds_ref) == 0){
- return(filter(missing_ds, member_COMID == x))
- }
-
- lookup_test <- filter(lookup, NHDPlusV2_COMID %in% ds_ref)
- divides_test <- filter(divides_lu, ID %in% lookup_test$reconciled_ID) %>%
- filter(!is.na(POI_ID))
-
- print(unique(divides_test$POI_ID))
-
- if(length(unique(divides_test$POI_ID)) == 1){
- return(filter(missing_ds, member_COMID == x) %>%
- mutate(POI_ID = unique(divides_test$POI_ID)))
- } else {
- return(filter(missing_ds, member_COMID == x))
- }
-
- })
-
- term_refactored <- data.table::rbindlist(term_refactored[lengths(term_refactored) > 1],
- fill = TRUE) %>%
- st_as_sf()
-
- divides_lu <- filter(divides_lu, !member_COMID %in% term_refactored$member_COMID) %>%
- rbind(term_refactored)
-
- if(sum(is.na(divides_lu$POI_ID)) > 0) {
- divides_dem <- miss_term_assign(divides_lu, vpu_nhd, elev)
-
- divides_lu <- divides_lu %>%
- left_join(select(divides_dem, member_COMID, agg_ID),
- by = "member_COMID") %>%
- mutate(POI_ID = ifelse(!is.na(agg_ID), agg_ID, POI_ID),
- aggStep = ifelse(!is.na(agg_ID), "boundary DEM", aggStep)) %>%
- select(-agg_ID)
-
- if(exists("sinks_table")){
- sinks_table_fin <- filter(sinks_table, !member_COMID %in% divides_dem$member_COMID)
- sinks_table_fin <- data.table::rbindlist(list(sinks_table_fin,
- divides_dem), fill = TRUE)
- } else {
- sinks_table_fin <- divides_dem
- }
-
- write_sf(sinks_table_fin, ND_gpkg, ND_table)
- }
-
- } else {
- print ("all unaggregated catchments assigned")
- }
-
- divides_lu <- dissolve_holes(divides_lu)
- if(sum(is.na(divides_lu$POI_ID)) > 0){
- write_sf(filter(divides_lu, is.na(POI_ID)), ND_gpkg, missing_cats)
- }
-
- # Prob HRU - filter(all_hrus, POI_ID == 140402000209)
- all_hrus <- filter(divides_lu, !is.na(POI_ID)) %>%
- group_by(POI_ID) %>%
- summarize(do_union = TRUE) %>%
- sfheaders::sf_remove_holes(.) %>%
- nhdplusTools:::check_valid(.)
-
- write_sf(divides_lu, ND_gpkg, divides_nd)
- write_sf(all_hrus, ND_gpkg, HRU_layer)
-} else {
- divides_lu <- read_sf(ND_gpkg, divides_nd)
- all_hrus <- read_sf(ND_gpkg, HRU_layer)
-}
-```
-
-```{r Manual checks}
-noagg_divides <- read_sf(ND_gpkg, missing_cats) %>%
- select(row_ID, POI_ID_noagg = POI_ID) %>%
- st_drop_geometry()
-
-if(all(!is.na(noagg_divides$POI_ID_noagg))){
- print ("all unaggregated divides accounted for")
-
- divides_lu <- divides_lu %>%
- left_join(noagg_divides, by = "row_ID") %>%
- mutate(POI_ID = ifelse(!is.na(POI_ID_noagg), POI_ID_noagg, POI_ID),
- aggStep = ifelse(!is.na(aggStep), "manual", aggStep)) %>%
- select(-POI_ID_noagg)
-
- # Prob HRU - filter(all_hrus, POI_ID == 140402000209)
- all_hrus <- filter(divides_lu, !is.na(POI_ID)) %>%
- group_by(POI_ID) %>%
- summarize(do_union = TRUE) %>%
- sfheaders::sf_remove_holes(.) %>%
- nhdplusTools:::check_valid(.)
-
- write_sf(divides_lu, ND_gpkg, divides_nd)
- write_sf(all_hrus, ND_gpkg, HRU_layer)
-} else {
- print ("account for unaggregated divides")
-}
-
-```
diff --git a/workspace/runners/README.md b/workspace/runners/README.md
deleted file mode 100644
index d79d9db..0000000
--- a/workspace/runners/README.md
+++ /dev/null
@@ -1,19 +0,0 @@
-gfv2 runners
-
-## 1. make_docker_calls.Rmd
-
-Generates `navigate_calls.txt` and `hyRefactor_calls.txt` and checks already generated outputs for completeness.
-
-Intended to be run from within a running Docker / Jupyter instance started with docker-compose in the root of this project.
-
-## 2. make_shifter_calls.Rmd
-
-Clone of `make_docker_calls.Rmd` but tailored for the shifter-based containerization environment on Denali.
-
-## 3. batch.sh
-
-Useful on a linux environment with Docker available. Can run `navigate_calls.txt` and `hyRefactor_calls.txt` in parallel.
-
-## 4. run_all_R.R
-
-Simple runner to execute calls in parallel on Windows with Docker.
diff --git a/workspace/runners/batch.sh b/workspace/runners/batch.sh
deleted file mode 100644
index 5385014..0000000
--- a/workspace/runners/batch.sh
+++ /dev/null
@@ -1,13 +0,0 @@
-#!/bin/bash
-#SBATCH -A impd
-#SBATCH -N 1i
-#SBATCH -c 40
-#SBATCH -n 1
-#SBATCH --job-name=gfv2
-#SBATCH --time=24:00:00
-
-parallel --jobs 20 < workspace/runners/shifter_navigate_calls.txt
-
-shifter --volume="/absolute/path/gfv2/gfv2.0/workspace:/jupyter" --image=dblodgett/hydrogeoenv-custom:latest R -e "rmarkdown::render('/jupyter/04_merge.Rmd', output_file='/jupyter/04_merge.html')"
-
-parallel --tmpdir workspace/temp --delay 30s --jobs 10 < workspace/runners/shifter_hyRefactor_calls.txt
diff --git a/workspace/runners/cleaner.R b/workspace/runners/cleaner.R
deleted file mode 100644
index 304dea6..0000000
--- a/workspace/runners/cleaner.R
+++ /dev/null
@@ -1,70 +0,0 @@
-possible_vpu <- c("01", "08", "10L", "15", "02", "04", "05", "06",
- "07", "09", "03S", "03W",
- "03N", "10U", "11", "12",
- "13", "14", "16", "17", "18")
-
-
-vpu_codes <- possible_vpu
-
-# remove all reference fabric content
-reference <- FALSE
-# remove al POI content
-POI <- TRUE
-# remove refactored content leaving rpu-based reference content behind
-refactor <- TRUE
-# remove all aggregate content
-aggregate <- FALSE
-
-for(VPU in vpu_codes) {
-
- rm(rpu_code)
-
- source("R/config.R")
- source("R/utils.R")
-
- if(reference) {
- unlink(nav_gpkg)
- }
-
- if(POI) {
- sf::st_delete(nav_gpkg, final_poi_layer)
- sf::st_delete(nav_gpkg, lookup_table_refactor)
- unlink(temp_gpkg)
- }
-
- if(refactor) {
- unlink(rfc_gpkg)
- }
-
- if(aggregate) {
- unlink(gf_gpkg)
- unlink(ND_gpkg)
- }
-
- if(refactor | aggregate) {
- unlink(vpu_lookup_table_file)
- }
-
- for(r in rpu_codes$rpuid) {
- f <- paste0("cache/", r, "_refactor.gpkg")
- if(refactor & !reference) {
- sf::st_delete(f, lookup_table_refactor)
- sf::st_delete(f, refactored_layer)
- sf::st_delete(f, reconciled_layer)
- sf::st_delete(f, outlets_layer)
- sf::st_delete(f, divide_layer)
- unlink(paste0("cache/", r, "_lookup.csv"))
- }
-
- if(refactor & reference) {
- unlink(f)
- unlink(paste0("cache/", r, "_lookup.csv"))
- }
-
- if(aggregate) {
- unlink(paste0("cache/", r, "_aggregate.gpkg"))
- unlink(paste0("cache/", r, "_lookup.csv"))
- }
- }
-
-}
diff --git a/workspace/runners/local_batch.sh b/workspace/runners/local_batch.sh
deleted file mode 100644
index d3a1761..0000000
--- a/workspace/runners/local_batch.sh
+++ /dev/null
@@ -1,9 +0,0 @@
-#!/bin/bash
-
-docker run --mount type=bind,source=$PWD/workspace,target=/jupyter dblodgett/hydrogeoenv-custom:latest R -e "rmarkdown::render('/jupyter/01_NHD_prep.Rmd', output_file='/jupyter/reports/01_NHD_prep.html')"
-
-parallel --jobs 5 < workspace/runners/navigate_calls.txt
-
-parallel --jobs 5 < workspace/runners/hyRefactor_calls.txt
-
-parallel --jobs 5 < workspace/runners/finalize_calls.txt
diff --git a/workspace/runners/make_docker_calls.Rmd b/workspace/runners/make_docker_calls.Rmd
deleted file mode 100644
index 3d35450..0000000
--- a/workspace/runners/make_docker_calls.Rmd
+++ /dev/null
@@ -1,191 +0,0 @@
----
-title: "Make Docker Calls"
-output: html_document
----
-
-First chunk creates calls for NHD_navigate and POI_Collapse.
-
-Includes one environment variable for the Hydro Region.
-
-The output .txt file can be executed with GNU parallel like:
-`parallel --jobs 10 < workspace/navigate_calls.txt`
-
-```{r}
-data_paths <- jsonlite::read_json(file.path("../cache", "data_paths.json"))
-
-out_script <- "navigate_calls.txt"
-
-call <- "docker run"
-mount <- "--mount type=bind,source=$PWD/workspace,target=/jupyter"
-docker_image <- "dblodgett/hydrogeoenv-custom:latest"
-
-e1 <- "--env HYDREG=z!z"
-
-command <- paste("R -e \"rmarkdown::render('/jupyter/02_NHD_navigate.Rmd', output_file='/jupyter/reports/02_NHD_navigate_z!z.html')\" \n")
-
-hus <- c(paste0("0", c(1:2, 4:9)),
- "03S", "03W", "03N",
- "10U", "10L",
- paste0("1", c(1:8)))
-
-calls <- sapply(hus, function(hu) gsub("z!z", hu, paste(call, mount, e1, docker_image, command)))
-
-cat("", file = out_script, append = FALSE)
-for(i in 1:length(calls)) {
- cat(calls[i], file = out_script, append = TRUE)
-}
-
- hus
-```
-
-The second chunk checks for completion of the NHD_navigate and POI_Collapse steps.
-
-Any completely missing geopackages and specific missing layers are printed if any.
-
-```{r}
-expect_files <- paste0("reference_", hus, ".gpkg")
-
-have_files <- list.files("/jupyter/cache")
-
-(missing_files <- expect_files[!expect_files %in% have_files])
-
-have_files <- expect_files[expect_files %in% have_files]
-
-expect_layers <- c("POIs", "reference_flowline")
-
-layers <- lapply(file.path("/jupyter/cache", have_files), function(x) sf::st_layers(x)$name)
-
-check <- lapply(layers, function(x, expect_layers) {
- sapply(paste0("^", expect_layers, ".*"), function(y, x) {
- any(grepl(y, unlist(x)))
- }, x = x)
-}, expect_layers = expect_layers)
-
-for(f in seq_along(have_files)) {
- if(any(!check[[f]])) {
- print(paste0(have_files[f], ": ", paste(names(check[[f]])[which(!check[[f]])], collapse = ", ")))
- }
-}
-```
-
-The following creates docker commands for hyRefactor flines and cats.
-
-Three environment variables are set for the output geopackage pack, the Raster Processing Unit, and a cache to be used for error output environments.
-
-The txt file can be used with GNU parallel like:
-`parallel --jobs 10 < workspace/hyRefactor_calls.txt`
-
-```{r}
-out_script <- "hyRefactor_calls.txt"
-
-e1 <- "--env OUT_GPKG=/jupyter/cache/z!z_aggregate.gpkg"
-e11 <- "--env OUT_REFACTOR_GPKG=/jupyter/cache/z!z_refactor.gpkg"
-e2 <- "--env RPU_CODE=z!z"
-e3 <- "--env SPLIT_CACHE=/jupyter/temp/z!z_split.rda"
-e4 <- "--env LOOKUP_TABLE=/jupyter/cache/z!z_lookup.csv"
-
-command <- paste("R -e \"rmarkdown::render('/jupyter/05_hyRefactor_flines.Rmd', output_file='/jupyter/reports/05_hyRefactor_flines_z!z.html')\"",
- "-e \"rmarkdown::render('/jupyter/06-1_hyRefactor_cats.Rmd', output_file='/jupyter/reports/06-1_hyRefactor_cats_z!z.html')\"",
- "-e \"rmarkdown::render('/jupyter/06-2_aggregate_cats.Rmd', output_file='/jupyter/reports/06-2_aggregate_cats_z!z.html')\" \n")
-
-all_hus <- gsub("rpu_", "", names(data_paths$fdr))
-
-all_hus <- all_hus[!grepl("^2", all_hus)]
-
-errors <- list.files("temp", "*.rda$", full.names = TRUE)
-errors <- gsub(".rda", "", gsub("temp/", "", errors))
-
-hus <- all_hus[!all_hus %in% errors]
-
-calls <- sapply(hus, function(hu) gsub("z!z", hu, paste(call, mount, e1, e11, e2, e3, e4, docker_image, command)))
-
-cat("", file = out_script, append = FALSE)
-for(i in 1:length(calls)) {
- cat(calls[i], file = out_script, append = TRUE)
-}
-```
-
-The next chunk checks for the expected output geopackages from the above calls.
-
-```{r}
-agg <- paste0(all_hus, "_aggregate.gpkg")
-ref <- paste0(all_hus, "_refactor.gpkg")
-expect_files <- c(ref, agg)
-
-have_files <- list.files("/jupyter/cache")
-
-(missing_files <- expect_files[!expect_files %in% have_files])
-```
-
-The next chunk checks which existing geopackages are missing the reconciled output.
-
-```{r}
-ref <- ref[ref %in% have_files]
-
-missing_reconciled <- sapply(file.path("/jupyter/cache", ref),
- function(x) !"refactored_flowpaths" %in% sf::st_layers(x)$name)
-
-missing_reconciled <- missing_reconciled[missing_reconciled]
-
-dplyr::arrange(tibble::tibble(file = names(missing_reconciled)))
-```
-
-The next chunk checks which divides have been competed and which haven't. This is the first step in catchment refactoring.
-
-```{r}
-missing_div <- sapply(file.path("/jupyter/cache", ref),
- function(x) !"refactored_divides" %in% sf::st_layers(x)$name)
-
-has_div <- missing_div[!missing_div]
-
-missing_div <- missing_div[missing_div]
-
-dplyr::arrange(tibble::tibble(file = names(missing_div)), file)
-
-# dplyr::arrange(tibble::tibble(file = names(has_div)), file)
-```
-
-The next chunk checks which final aggregation have been completed and which haven't.
-
-```{r}
-agg <- agg[agg %in% have_files]
-
-missing_agg <- sapply(file.path("/jupyter/cache", agg),
- function(x) !"aggregated_divides" %in% sf::st_layers(x)$name)
-
-has_agg <- missing_agg[!missing_agg]
-
-missing_agg <- missing_agg[missing_agg]
-
-dplyr::arrange(tibble::tibble(file = names(missing_agg)), file)
-
-dplyr::arrange(tibble::tibble(file = names(has_agg)), file)
-```
-
-```{r}
-out_script <- "finalize_calls.txt"
-
-e1 <- "--env HYDREG=z!z"
-
-command <- paste("R -e \"rmarkdown::render('/jupyter/07-1_merge.Rmd', output_file='/jupyter/reports/07-1_merge_z!z.html')\"",
- "-e \"rmarkdown::render('/jupyter/07-2_NonDend.Rmd', output_file='/jupyter/reports/07-2_NonDend_z!z.html')\" \n")
-
-command <- paste("R -e \"rmarkdown::render('/jupyter/07-1_merge.Rmd', output_file='/jupyter/reports/07-1_merge_z!z.html')\"",
- "-e \"rmarkdown::render('/jupyter/07-2_NonDend.Rmd', output_file='/jupyter/reports/07-2_NonDend_z!z.html')\" \n")
-
-hus <- c(paste0("0", c(1:2, 4:9)),
- "03S", "03W", "03N",
- "10U", "10L",
- paste0("1", c(1:8)))
-
-calls <- sapply(hus, function(hu) gsub("z!z", hu, paste(call, mount, e1, docker_image, command)))
-
-cat("", file = out_script, append = FALSE)
-for(i in 1:length(calls)) {
- cat(calls[i], file = out_script, append = TRUE)
-}
-```
-
-```{r}
-
-```
diff --git a/workspace/runners/make_shifter_calls.Rmd b/workspace/runners/make_shifter_calls.Rmd
deleted file mode 100644
index 045fdd2..0000000
--- a/workspace/runners/make_shifter_calls.Rmd
+++ /dev/null
@@ -1,63 +0,0 @@
----
- title: "Make Shifter Calls"
- output: html_document
----
-
- First chunk creates calls for NHD_navigate and POI_Collapse.
-
-Includes one environment variable for the Hydro Region.
-
-The output .txt file can be executed with GNU parallel like:
- `parallel --jobs 10 < workspace/shifter_navigate_calls.txt`
-
-```{r}
-data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
-
-out_script <- "shifter_navigate_calls.txt"
-
-call <- "shifter"
-mount <- '--volume="/absolute/path/gfv2/gfv2.0/workspace:/jupyter"'
-docker_image <- '--image=dblodgett/hydrogeoenv-custom:latest'
-
-command <- paste("R -e \"Sys.setenv(HYDREG='z!z')\"",
- "-e \"rmarkdown::render('/jupyter/02_NHD_navigate.Rmd', output_file='/jupyter/reports/02_NHD_navigate_z!z.html')\"",
- "-e \"rmarkdown::render('/jupyter/03_POI_Collapse.Rmd', output_file='/jupyter/reports/03_POI_Collapse_z!z.html')\" \n")
-
-hus <- c(paste0("0", c(1:9)), paste0("1", c(1:8)), "10U", "10L")
-
-calls <- sapply(hus, function(hu) gsub("z!z", hu, paste(call, mount, docker_image, command)))
-
-cat("", file = out_script, append = FALSE)
-for(i in 1:length(calls)) {
- cat(calls[i], file = out_script, append = TRUE)
-}
-```
-
-The following creates docker commands for hyRefactor flines and cats.
-
-Three environment variables are set for the output geopackage pack, the Raster Processing Unit, and a cache to be used for error output environments.
-
-The txt file can be used with GNU parallel like:
- `parallel --jobs 10 < workspace/shifter_hyRefactor_calls.txt`
-
-```{r}
-out_script <- "shifter_hyRefactor_calls.txt"
-
-command <- paste("R -e \"Sys.setenv(OUT_GPKG='/jupyter/cache/z!z.gpkg', RPU_CODE='z!z', SPLIT_CACHE='/jupyter/temp/z!z_split.rda', LOOKUP_TABLE='/jupyter/cache/z!z_lookup.csv')\"",
- "-e \"rmarkdown::render('/jupyter/05_hyRefactor_flines.Rmd', output_file='/jupyter/reports/05_hyRefactor_flines_z!z.html')\"",
- "-e \"rmarkdown::render('/jupyter/07_hyRefactor_cats.Rmd', output_file='/jupyter/reports/07_hyRefactor_cats_z!z.html')\" \n")
-
-all_hus <- gsub("rpu_", "", names(data_paths$fdr))
-
-errors <- list.files("temp", "*.rda$", full.names = TRUE)
-errors <- gsub(".rda", "", gsub("temp/", "", errors))
-
-hus <- all_hus[!all_hus %in% errors]
-
-calls <- sapply(hus, function(hu) gsub("z!z", hu, paste(call, mount, docker_image, command)))
-
-cat("", file = out_script, append = FALSE)
-for(i in 1:length(calls)) {
- cat(calls[i], file = out_script, append = TRUE)
-}
-```
diff --git a/workspace/runners/run_all_R.R b/workspace/runners/run_all_R.R
deleted file mode 100644
index 228c631..0000000
--- a/workspace/runners/run_all_R.R
+++ /dev/null
@@ -1,16 +0,0 @@
-rmarkdown::render('01_NHD_prep.Rmd',
- output_file='temp/01_NHD_prep_test.html')
-
-setwd("../")
-
-calls <- readLines("workspace/runners/navigate_calls.txt")
-
-cl <- parallel::makeCluster(5, outfile = "workspace/temp/error.log")
-
-pbapply::pblapply(calls, shell, cl = cl)
-
-calls <- readLines("workspace/runners/hyRefactor_calls.txt")
-
-parallel::parLapply(cl = cl, X = calls, fun = shell)
-
-parallel::stopCluster(cl)
diff --git a/workspace/runners/run_one_R.R b/workspace/runners/run_one_R.R
deleted file mode 100644
index 1bfba3d..0000000
--- a/workspace/runners/run_one_R.R
+++ /dev/null
@@ -1,74 +0,0 @@
-possible_vpu <- c("01", "08", "10L", "15", "02", "04", "05", "06",
- "07", "09", "03S", "03W",
- "03N", "10U", "11", "12",
- "13", "14", "16", "17", "18")
-
-
-vpu_codes <- VPU <- vpu <- "01"
-
-check_file <- file.path("temp/", paste0(VPU, ".lock"))
-
-if(!file.exists(check_file)) {
- cat("...", file = check_file)
-
- source("R/config.R")
- source("R/utils.R")
-
- unlink(rfc_gpkg)
- unlink(nav_gpkg)
- unlink(ND_gpkg)
- unlink(gf_gpkg)
- unlink(temp_gpkg)
- unlink(vpu_lookup_table_file)
-
- for(r in rpu_codes$rpuid) {
- unlink(paste0("cache/", r, "_refactor.gpkg"))
- }
-
- rmarkdown::render('01_NHD_prep.Rmd',
- output_file=paste0('temp/01_NHD_prep_', VPU, '.html'))
-
- rmarkdown::render('02_NHD_navigate.Rmd',
- output_file=paste0('temp/02_NHD_navigate_', VPU, '.html'))
-
- for(rpu_code in rpu_codes$rpuid) {
-
- source("R/config.R")
-
- unlink(out_agg_gpkg)
- unlink(lookup_table_file)
-
- rmarkdown::render("05_hyRefactor_flines.Rmd",
- output_file = paste0("temp/05_hyRefactor_flines_", rpu_code, ".html"))
-
- rmarkdown::render("06-1_hyRefactor_cats.Rmd",
- output_file = paste0("temp/06-1_hyRefactor_cats_", rpu_code, ".html"))
-
- rmarkdown::render("06-2_aggregate_cats.Rmd",
- output_file = paste0("temp/06-2_aggregate_cats_", rpu_code, ".html"))
- }
-
- rmarkdown::render("07-1_merge.Rmd",
- output_file = "temp/07-1_merge_test.html")
-
- rmarkdown::render("07-2_NonDend.Rmd",
- output_file = "temp/07-2_NonDend_test.html")
-
- if(is.na(Sys.getenv("sb_user", unset=NA))){
- message("upload skipped due to lack of login info")
- } else {
-
- library(sbtools)
-
- authenticate_sb(Sys.getenv("sb_user", unset=""), Sys.getenv("sb_pass", unset=""))
-
- sb_reference <- "61295190d34e40dd9c06bcd7"
- sb_refactor <- "61fbfdced34e622189cb1b0a"
- sb_gf <- "60be1502d34e86b9389102cc"
-
- sbtools::item_replace_files(sb_reference, nav_gpkg)
- sbtools::item_replace_files(sb_refactor, rfc_gpkg)
- sbtools::item_replace_files(sb_gf, gf_gpkg)
- }
-
-}
--
GitLab