Skip to content
Snippets Groups Projects
Commit 07a0090c authored by Blodgett, David L.'s avatar Blodgett, David L.
Browse files

conflict

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