Skip to content
Snippets Groups Projects
Commit 2a830f5d authored by Bock, Andy's avatar Bock, Andy
Browse files

Major mods see MR for details

parent 131b63a1
No related branches found
No related tags found
1 merge request!169Updates through 07_merge
Hide whitespace changes
Inline Side-by-side
workspace/02_NHD_navigate.Rmd
+ 265
− 146
View file @ 2a830f5d
......@@ -22,8 +22,6 @@ knitr::opts_chunk$set(
source("R/utils.R")
source("R/NHD_navigate.R")
source("R/hyrefactor_funs.R")
#source("R/wb_poi_collapse.R")
#source("R/wb_inlet_collapse.R")
# increase timeout for data downloads
options(timeout=600)
......@@ -35,15 +33,15 @@ source("R/config.R")
gages <- read_sf(gage_info_gpkg, "Gages")
# NWM network
nc <- RNetCDF::open.nc(data_paths$nwm_network)
#nc <- RNetCDF::open.nc(data_paths$nwm_network)
# NWM gages
nwm_gages <- data.frame(
comid = RNetCDF::var.get.nc(nc, "link"),
gage_id = RNetCDF::var.get.nc(nc, "gages")) %>%
mutate(gage_id = gsub(" ", "", gage_id)) %>%
mutate(gage_id = ifelse(gage_id == "", NA, gage_id))
#nwm_gages <- data.frame(
# comid = RNetCDF::var.get.nc(nc, "link"),
# gage_id = RNetCDF::var.get.nc(nc, "gages")) %>%
# mutate(gage_id = gsub(" ", "", gage_id)) %>%
# mutate(gage_id = ifelse(gage_id == "", NA, gage_id))
RNetCDF::close.nc(nc)
#RNetCDF::close.nc(nc)
# need some extra attributes for a few POI analyses
atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.rds"))
......@@ -54,7 +52,7 @@ atts <- readRDS(file.path(data_paths$nhdplus_dir, "nhdplus_flowline_attributes.r
if(needs_layer(temp_gpkg, nav_poi_layer)) {
nhd <- read_sf(nav_gpkg, nhd_flowline)
try(nhd <- select(nhd, -c(minNet, WB, struct_POI, struct_Net, POI_ID, dend, poi)), silent = TRUE)
#try(nhd <- select(nhd, -c(minNet, WB, struct_POI, struct_Net, POI_ID, dend)), silent = TRUE)
# Some NHDPlus VPUs include HUCs from other VPUs
if(vpu == "02"){
......@@ -90,13 +88,12 @@ mapview(filter(tmp_POIs, Type_HUC12 != ""), layer.name = "HUC12 POIs", col.regio
```
```{r streamgage POIs}
if(all(is.na(tmp_POIs$Type_Gages))) {
if(!"Type_Gages" %in% names(tmp_POIs)) {
# Previously identified streamgages within Gage_Selection.Rmd
streamgages_VPU <- gages %>%
rename(COMID = comid) %>%
filter(COMID %in% nhd$COMID) %>%
#st_drop_geometry() %>%
switchDiv(., nhd)
streamgages <- streamgages_VPU %>%
......@@ -106,8 +103,9 @@ if(all(is.na(tmp_POIs$Type_Gages))) {
ungroup()
# Derive GAGE POIs; use NHD as we've already filtered by NWIS DA in the Gage selection step
gages_POIs <- gage_POI_creation(tmp_POIs, streamgages, nhd, combine_meters, reach_meas_thresh)
gages_POIs <- gage_POI_creation(tmp_POIs, streamgages, filter(nhd, poi == 1), combine_meters, reach_meas_thresh)
tmp_POIs <- gages_POIs$tmp_POIs
events <- rename(gages_POIs$events, POI_identifier = Type_Gages)
# As a fail-safe, write out list of gages not assigned a POI
if(nrow(filter(streamgages_VPU, !site_no %in% tmp_POIs$Type_Gages)) > 0) {
......@@ -122,7 +120,7 @@ if(all(is.na(tmp_POIs$Type_Gages))) {
}
# Write out events and outelts
write_sf(rename(gages_POIs$events, POI_identifier = Type_Gages), temp_gpkg, split_layer)
write_sf(events, temp_gpkg, split_layer)
# Write out geopackage layer representing POIs for given theme
write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
} else {
......@@ -134,7 +132,7 @@ mapview(filter(tmp_POIs, !is.na(Type_Gages)), layer.name = "Streamgage POIs", co
```
```{r TE POIs}
if(all(is.na(tmp_POIs$Type_TE))) {
if(!"Type_TE" %in% names(tmp_POIs)) {
if(vpu == "08"){
nhd$VPUID <- "08"
......@@ -153,7 +151,7 @@ if(all(is.na(tmp_POIs$Type_TE))) {
ungroup()
# Derive TE POIs
tmp_POIs <- POI_creation(st_drop_geometry(TE_COMIDs), filter(nhd, dend == 1), "TE") %>%
tmp_POIs <- POI_creation(st_drop_geometry(TE_COMIDs), filter(nhd, poi == 1), "TE") %>%
addType(., tmp_POIs, "TE")
# As a fail-safe, write out list of TE plants not assigned a POI
......@@ -172,38 +170,9 @@ if(all(is.na(tmp_POIs$Type_TE))) {
mapview(filter(tmp_POIs, Type_TE != ""), layer.name = "TE Plant POIs", col.regions = "blue")
```
```{r Terminal POIs}
# Derive or load Terminal POIs ----------------------
if(all(is.na(tmp_POIs$Type_Term))) {
# Terminal POIs on levelpaths with upstream POIs
term_paths <- filter(st_drop_geometry(nhd), Hydroseq %in% filter(nhd, COMID %in% tmp_POIs$COMID)$TerminalPa)
# Non-POI levelpath terminal pois, but meet size criteria
terminal_POIs <- st_drop_geometry(nhd) %>%
filter(Hydroseq == TerminalPa | (toCOMID == 0 | is.na(toCOMID) | !toCOMID %in% COMID)) %>%
filter(!COMID %in% term_paths$COMID, TotDASqKM >= min_da_km) %>%
bind_rows(term_paths) %>%
# Use level path identifier as Type ID
select(COMID, LevelPathI)
tmp_POIs <- POI_creation(terminal_POIs, nhd, "Term") %>%
addType(., tmp_POIs, "Term")
write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
} else {
tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
nhd <- read_sf(nav_gpkg, nhd_flowline)
}
mapview(filter(tmp_POIs, !is.na(Type_Term)), layer.name = "Terminal POIs", col.regions = "blue")
```
```{r waterbody outlet POIs}
# Derive or load Waterbody POIs ----------------------
if(all(is.na(tmp_POIs$Type_WBOut))) {
events <- read_sf(temp_gpkg, split_layer)
if(!"Type_WBOut" %in% names(tmp_POIs)) {
# Waterbodies sourced from NHD waterbody layer for entire VPU
WBs_VPU_all <- filter(readRDS("data/NHDPlusNationalData/nhdplus_waterbodies.rds"),
......@@ -211,6 +180,10 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
mutate(FTYPE = as.character(FTYPE)) %>%
mutate(source = "NHDv2WB") %>%
st_transform(crs)
ref_WB <- read_sf(nav_gpkg, nhd_waterbody) %>%
filter(COMID %in% nhd$WBAREACOMI) %>%
mutate(source = "Ref_WB")
# Waterbody list that strictly meet the size criteria
wb_lst <- st_drop_geometry(WBs_VPU_all) %>%
......@@ -239,14 +212,15 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
# Attribute flowlines that are in waterbodies
nhd <- mutate(nhd, WB = ifelse(WBAREACOMI > 0 & WBAREACOMI %in% WBs_VPU$COMID, 1, 0))
wb_layers <- wbout_POI_creaton(nhd, WBs_VPU, data_paths, crs)
wb_layers <- wbout_POI_creaton(filter(nhd, WB == 1), WBs_VPU, data_paths, crs)
if(!is.na(wb_layers$events) > 0) {events <- rbind(events, wb_layers$events)}
WBs_VPU <- wb_layers$WBs
#tmp_POIs <- wb_layers$POIs
wb_out_col <- wb_poi_collapse(wb_layers$POIs, nhd, events)
ho <- filter(wb_layers$POIs, !COMID %in% wb_out_col$POIs$COMID)
write_sf(wb_out_col$events_ret, temp_gpkg, split_layer)
wb_out_col <- wb_poi_collapse(wb_layers$POIs, filter(nhd, WB == 1), events)
tmp_POIs <- wb_out_col$POIs
if(!is.na(wb_out_col$events_ret)){
write_sf(wb_out_col$events_ret, temp_gpkg, split_layer)
}
if(!all(is.na(wb_layers$events))){
wb_events <- select(wb_layers$events, -c(id, offset)) %>%
......@@ -263,8 +237,22 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
}
}
ref_WB <- filter(ref_WB, COMID %in% WBs_VPU$COMID)
WBs_VPU <- filter(WBs_VPU, !COMID %in% ref_WB$COMID) %>%
group_by(GNIS_ID, GNIS_NAME, COMID, FTYPE, source) %>%
summarize(do_union = T) %>%
ungroup() %>%
st_make_valid() %>%
sfheaders::sf_remove_holes(.) %>%
mutate(member_comid = NA,
area_sqkm = as.numeric(st_area(.))/1000000) %>%
st_compatibalize(., ref_WB)
ref_WB <- rbind(ref_WB, WBs_VPU)
# Write out waterbodies
write_sf(WBs_VPU, nav_gpkg, WBs_layer)
write_sf(ref_WB, temp_gpkg, WBs_layer)
# Write specific ResOpsUS data to a separate sheet
if(nrow(wb_lst) > 0){
......@@ -281,16 +269,43 @@ if(all(is.na(tmp_POIs$Type_WBOut))) {
} else {
tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
nhd <- read_sf(nav_gpkg, nhd_flowline)
WBs_VPU <- read_sf(nav_gpkg, WBs_layer)
ref_WB <- read_sf(temp_gpkg, WBs_layer)
resops_POIs_df <- read.csv(file.path("data/reservoir_Data", paste0("resops_POIs_",vpu,".csv")))
}
mapview(filter(tmp_POIs, !is.na(Type_WBOut)), layer.name = "Waterbody outlet POIs", col.regions = "red")
```
```{r Terminal POIs}
# Derive or load Terminal POIs ----------------------
if(!"Type_Term" %in% names(tmp_POIs)) {
# Terminal POIs on levelpaths with upstream POIs
term_paths <- filter(st_drop_geometry(nhd), Hydroseq %in% filter(nhd, COMID %in% tmp_POIs$COMID)$TerminalPa)
# Non-POI levelpath terminal pois, but meet size criteria
terminal_POIs <- st_drop_geometry(nhd) %>%
filter(Hydroseq == TerminalPa | (toCOMID == 0 | is.na(toCOMID) | !toCOMID %in% COMID)) %>%
filter(!COMID %in% term_paths$COMID, TotDASqKM >= min_da_km) %>%
bind_rows(term_paths) %>%
# Use level path identifier as Type ID
select(COMID, LevelPathI)
tmp_POIs <- POI_creation(terminal_POIs, filter(nhd, poi == 1), "Term") %>%
addType(., tmp_POIs, "Term")
write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
} else {
tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
nhd <- read_sf(nav_gpkg, nhd_flowline)
}
mapview(filter(tmp_POIs, !is.na(Type_Term)), layer.name = "Terminal POIs", col.regions = "blue")
```
```{r Confluence POIs}
# Derive POIs at confluences where they are absent ----------------------
if(all(is.na(tmp_POIs$Type_Conf))) {
if(!"Type_Conf" %in% names(tmp_POIs)) {
# Navigate upstream from each POI and determine minimally-sufficient network between current POI sets
up_net <- unique(unlist(lapply(unique(tmp_POIs$COMID), NetworkNav, nhd)))
......@@ -325,16 +340,53 @@ if(all(is.na(tmp_POIs$Type_Conf))) {
mapview(filter(tmp_POIs, !is.na(Type_Conf)), layer.name = "Confluence POIs", col.regions = "blue")
```
```{r waterbody inlet POIs}
# Derive or load Waterbody POIs ----------------------
if(!"Type_WBIn" %in% names(tmp_POIs)) {
wb_layers <- wbin_POIcreation(nhd, ref_WB, data_paths, crs, split_layer, tmp_POIs)
wb_in_col <- wb_inlet_collapse(wb_layers$POIs, nhd, events)
#ho <- filter(wb_layers$POIs, !COMID %in% wb_in_col$POIs$COMID)
#write_sf(wb_in_col$events_ret, temp_gpkg, split_layer)
tmp_POIs <- wb_in_col$POIs
if(!all(is.na(wb_layers$events))) {
wb_inlet_events <- select(wb_layers$events, -c(id, offset, Hydroseq, ToMeas)) %>%
rename(POI_identifier = WBAREACOMI)
# Write out events and outlets
if(!needs_layer(temp_gpkg, split_layer)){
events <- read_sf(temp_gpkg, split_layer) %>%
rbind(st_compatibalize(wb_inlet_events,.))
write_sf(events, temp_gpkg, split_layer)
} else {
write_sf(wb_inlet_events, nav_gpkg, split_layer)
}
}
tmp_POIs$nexus <- ifelse(is.na(tmp_POIs$nexus), FALSE, tmp_POIs$nexus)
write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
} else {
tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
nhd <- read_sf(nav_gpkg, nhd_flowline)
}
mapview(filter(tmp_POIs, !is.na(Type_WBIn)), layer.name = "Waterbody inlet POIs", col.regions = "blue") +
mapview(filter(tmp_POIs, !is.na(Type_WBOut)), layer.name = "Waterbody outlet POIs", col.regions = "red")
```
```{r NID POIs}
# Derive or load NID POIs ----------------------
if(all(is.na(tmp_POIs$Type_NID))) {
if(!"Type_NID" %in% names(tmp_POIs)) {
# Read in NID shapefile
NID_COMIDs <- read_sf(data_paths$NID_points_path, "Final_NID_2018") %>%
st_drop_geometry() %>%
filter(EROM != 0, FlowLcomid %in% filter(nhd, dend ==1)$COMID) %>%
rename(COMID = FlowLcomid) %>%
switchDiv(., nhd) %>%
group_by(FlowLcomid) %>%
group_by(COMID) %>%
summarize(Type_NID = paste0(unique(NIDID), collapse = " "))
# Determine which NID facilitites have been co-located with ResOps
......@@ -349,14 +401,14 @@ if(all(is.na(tmp_POIs$Type_NID))) {
if(nrow(res_ops_table) > 0){
tmp_POIs <- tmp_POIs %>%
left_join(res_ops_table, by = "Type_WBOut") %>%
mutate(Type_NID = ifelse(!is.na(NID_ID), NID_ID, Type_NID)) %>%
mutate(Type_NID = ifelse(!is.na(NID_ID), NID_ID, NA)) %>%
select(-NID_ID)
NID_COMIDs <- filter(NID_COMIDs, !Type_NID %in% res_ops_table$NID_ID)
}
# Derive other NID POIs
tmp_POIs <- POI_creation(NID_COMIDs, filter(nhd, dend ==1), "NID") %>%
tmp_POIs <- POI_creation(NID_COMIDs, filter(nhd, poi == 1), "NID") %>%
addType(., tmp_POIs, "NID", nexus = TRUE, bind = FALSE)
# Write out geopackage layer representing POIs for given theme
......@@ -369,43 +421,6 @@ if(all(is.na(tmp_POIs$Type_NID))) {
mapview(filter(tmp_POIs, Type_NID != ""), layer.name = "NID POIs", col.regions = "blue")
```
```{r waterbody inlet POIs}
# Derive or load Waterbody POIs ----------------------
if(all(is.na(tmp_POIs$Type_WBIn))) {
wb_layers <- wbin_POIcreation(nhd, WBs_VPU, data_paths, crs, split_layer, tmp_POIs)
#tmp_POIs <- wb_layers$POIs
wb_in_col <- wb_inlet_collapse(wb_layers$POIs, nhd, events)
#ho <- filter(wb_layers$POIs, !COMID %in% wb_in_col$POIs$COMID)
#write_sf(wb_in_col$events_ret, temp_gpkg, split_layer)
tmp_POIs <- wb_in_col$POIs
if(!all(is.na(wb_layers$events))) {
wb_inlet_events <- select(wb_layers$events, -c(id, offset, Hydroseq, ToMeas)) %>%
rename(POI_identifier = WBAREACOMI)
# Write out events and outlets
if(!needs_layer(temp_gpkg, split_layer)){
events <- read_sf(temp_gpkg, split_layer) %>%
rbind(st_compatibalize(wb_inlet_events,.))
write_sf(events, temp_gpkg, split_layer)
} else {
write_sf(wb_inlet_events, nav_gpkg, split_layer)
}
}
tmp_POIs$nexus <- ifelse(is.na(tmp_POIs$nexus), FALSE, tmp_POIs$nexus)
write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
} else {
tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
nhd <- read_sf(nav_gpkg, nhd_flowline)
}
mapview(filter(tmp_POIs, !is.na(Type_WBIn)), layer.name = "Waterbody inlet POIs", col.regions = "blue") +
mapview(filter(tmp_POIs, !is.na(Type_WBOut)), layer.name = "Waterbody outlet POIs", col.regions = "red")
```
This chunk breaks up potential aggregated segments where the elevation range of the are contributing to the segment exceeds 500-m
```{r Elevation Break POIs}
# derive incremental segments from POIs
......@@ -424,7 +439,7 @@ elev_pois_split <- inc_segs %>%
mutate(MAXELEVSMO = na_if(MAXELEVSMO, -9998), MINELEVSMO = na_if(MINELEVSMO, -9998),
elev_diff_seg = max(MAXELEVSMO) - min(MINELEVSMO)) %>%
# Filter out to those incremental segments that need splitting above the defined elevation threshold
filter((max(MAXELEVSMO) - min(MINELEVSMO)) > (elev_diff * 100)) %>%
filter((max(MINELEVSMO) - min(MINELEVSMO)) > (elev_diff * 100)) %>%
# Do not aggregate on fake lfowlines
filter(AreaSqKM > 0, TotDASqKM > max_elev_TT_DA) %>%
# Order from upstream to downstream
......@@ -506,18 +521,18 @@ if(all(is.na(tmp_POIs$Type_Travel))) {
mutate(orig_iter = iter) %>%
ungroup()
# For reach iteration
tt_POIs <- do.call("rbind", lapply(unique(tt_pois_split$POI_ID), split_tt,
tt_pois_split, travt_diff, max_elev_TT_DA)) %>%
filter(!tt_POI_ID %in% tmp_POIs$COMID, COMID == tt_POI_ID) %>%
mutate(Type_Travel = 1) %>%
select(COMID, Type_Travel) %>%
unique()
tmp_POIs <- POI_creation(select(tt_POIs, COMID, Type_Travel), nhd, "Travel") %>%
addType(., tmp_POIs, "Travel")
# For reach iteration
tt_POIs <- do.call("rbind", lapply(unique(tt_pois_split$POI_ID), split_tt,
tt_pois_split, travt_diff, max_elev_TT_DA)) %>%
filter(!tt_POI_ID %in% tmp_POIs$COMID, COMID == tt_POI_ID) %>%
mutate(Type_Travel = 1) %>%
select(COMID, Type_Travel) %>%
unique()
write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
tmp_POIs <- POI_creation(select(tt_POIs, COMID, Type_Travel), nhd, "Travel") %>%
addType(., tmp_POIs, "Travel")
write_sf(tmp_POIs, temp_gpkg, nav_poi_layer)
}else{
tmp_POIs <- read_sf(temp_gpkg, nav_poi_layer)
nhd <- read_sf(nav_gpkg, nhd_flowline)
......@@ -611,48 +626,89 @@ if(needs_layer(temp_gpkg, nsegments_layer)) {
# Ensure that all the problem POIs have been taken care of
sub <- nhd_Final %>% filter(COMID %in% final_POIs$COMID)
print (paste0(nrow(sub[sub$AreaSqKM == 0,]), " POIs on flowlines with no local drainage contributions"))
noDA_pois <- filter(final_POIs, COMID %in% filter(sub, AreaSqKM == 0)$COMID)
write_sf(noDA_pois, temp_gpkg, "noDA_pois")
```
```{r POI Collapse}
# Load data
if(needs_layer(nav_gpkg, final_poi_layer)) {
# number POIs
final_POIs <- mutate(final_POIs, id = row_number(), moved = NA)
collapse <- TRUE
if(!"Type_Elev" %in% names(final_POIs)) final_POIs$Type_Elev <- rep(NA_real_, nrow(final_POIs))
# Convert nexus to 0,1
final_POIs <- mutate(final_POIs, nexus = ifelse(is.na(nexus), 0, nexus))
# Load data
if(collapse){
# Move HUC12 to other POIs
moved_pois <- poi_move(final_POIs, "Type_HUC12", .05, 2.5, c("Type_Gages", "Type_TE", "Type_NID"))
if(!is.data.frame(moved_pois)){
final_POIs <- moved_pois$final_pois
moved_pois_table <- moved_pois$moved_points %>%
mutate(move_type = "huc12 to other")
} else {
final_POIs <- moved_POIs
}
# Gages to confluences
moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_Conf") # 47
if(!is.data.frame(moved_pois)){
final_POIs <- moved_pois$final_pois
moved_pois_table <- moved_pois_table %>%
rbind(moved_pois$moved_points %>%
mutate(move_type = "gages to conf"))
} else {
final_POIs <- moved_POIs
}
#1 Move POIs downstream by category, nexus POIs get removed
out_HUC12 <- POI_move_down(temp_gpkg, final_POIs, nsegments, filter(nhd_Final, !is.na(POI_ID)), "Type_HUC12",
poi_dar_move *2)
out_HUC12$allPOIs$nexus <- as.numeric(out_HUC12$allPOIs$nexus)
out_gages <- POI_move_down(temp_gpkg, out_HUC12$allPOIs, out_HUC12$segs, out_HUC12$FL,
"Type_Gages", poi_dar_move)
# Convert empty strings to NA for handling within R
out_gages$allPOIs <- mutate_all(out_gages$allPOIs, list(~na_if(.,"")))
# Gages to waterbody inlets
moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_WBIn") # 2
if(!is.data.frame(moved_pois)){
final_POIs <- moved_pois$final_pois
moved_pois_table <- moved_pois_table %>%
rbind(moved_pois$moved_points %>%
mutate(move_type = "gages to wbin"))
} else {
final_POIs <- moved_pois
}
out_gages$allPOIs$snapped <- out_gages$allPOIs$COMID %in% new_POIs$COMID
out_gages$allPOIs$identifier <- seq(1, nrow(out_gages$allPOIs))
# Gages to waterbody outlets
moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_WBOut") # 6
if(!is.data.frame(moved_pois)){
final_POIs <- moved_pois$final_pois
moved_pois_table <- moved_pois_table %>%
rbind(moved_pois$moved_points %>%
mutate(move_type = "gages to wbout"))
} else {
final_POIs <- moved_pois
}
out_gages$allPOIs <- select(out_gages$allPOIs, identifier, everything())
nhd_Final <- select(nhd_Final, -POI_ID) %>%
left_join(st_drop_geometry(out_HUC12$FL) %>%
select(COMID, POI_ID), by = "COMID")
# Write out geopackage layer representing POIs for given theme
write_sf(out_gages$allPOIs, nav_gpkg, final_poi_layer)
write_sf(out_gages$segs, temp_gpkg, nsegments_layer)
write_sf(nhd_Final, temp_gpkg, nhd_flowline)
# Streamgage to term outlet
moved_pois <- poi_move(final_POIs, "Type_Gages", .05, 2.5, "Type_Term")
if(!is.data.frame(moved_pois)){
final_POIs <- moved_pois$final_pois
moved_pois_table <- moved_pois_table %>%
rbind(moved_pois$moved_points %>%
mutate(move_type = "gages to term"))
} else {
final_POIs <- moved_pois
}
nsegments <- out_gages$segs
final_POIs <- out_gages$allPOIs
} else {
final_POIs <- read_sf(nav_gpkg, final_poi_layer)
nhd_Final <- read_sf(nav_gpkg, nhd_flowline)
nsegments <- read_sf(temp_gpkg, nsegments_layer)
}
# NID to waterbody outlet
moved_pois <- poi_move(final_POIs, "Type_NID", .025, 1, "Type_WBOut")
if(!is.data.frame(moved_pois)){
final_POIs <- moved_pois$final_pois
moved_pois_table <- moved_pois_table %>%
rbind(moved_pois$moved_points %>%
mutate(move_type = "nid to wnpit"))
} else {
final_POIs <- moved_pois
}
write_sf(moved_pois_table, temp_gpkg, "pois_collapsed")
}
write_sf(final_POIs, nav_gpkg, pois_all_layer)
check_dups <- final_POIs %>%
group_by(COMID) %>%
......@@ -661,7 +717,7 @@ check_dups <- final_POIs %>%
if(nrow(filter(check_dups, all(c(0,1) %in% nexus))) != nrow(check_dups)){
print("Multiple POI ids at same geometric location")
no_dups <- filter(check_dups, all(c(0,1) %in% nexus))
dups <- filter(check_dups, !identifier %in% no_dups$identifier)
dups <- filter(check_dups, !id %in% no_dups$id)
write_sf(dups, temp_gpkg, dup_pois)
} else {
print("All double COMIDs nexus for gage or WB splitting")
......@@ -672,19 +728,82 @@ mapview(final_POIs, layer.name = "All POIs", col.regions = "red") +
```
```{r lookup}
# Final POI layer
POIs <- final_POIs %>%
arrange(COMID) %>%
mutate(identifier = row_number())
# Unique POI geometry
final_POI_geom <- POIs %>%
select(identifier) %>%
cbind(st_coordinates(.)) %>%
group_by(X, Y) %>%
mutate(geom_id = cur_group_id()) %>%
ungroup()
final_POIs_table <- POIs %>%
inner_join(select(st_drop_geometry(final_POI_geom), -X, -Y), by = "identifier") %>%
select(-identifier)
# POI data theme table
pois_data <- 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 <- select(final_POIs_table, hy_id = COMID, geom_id) %>%
inner_join(distinct(pois_data, hy_id, geom_id, poi_id),
by = c("hy_id" = "hy_id", "geom_id" = "geom_id")) %>%
distinct() %>%
st_as_sf()
write_sf(pois_data, nav_gpkg, poi_data_table)
write_sf(poi_geometry, nav_gpkg, poi_geometry_table)
poi_geom_xy <- cbind(poi_geometry, st_coordinates(poi_geometry)) %>%
st_drop_geometry()
events_data <- events %>%
arrange(COMID) %>%
cbind(st_coordinates(.)) %>%
st_drop_geometry() %>%
group_by(COMID, REACHCODE, REACH_meas) %>%
mutate(event_id = cur_group_id()) %>%
rename(hy_id = COMID) %>%
ungroup()
nexi <- filter(final_POIs_table, nexus == 1) %>%
cbind(st_coordinates(.)) %>%
select(hy_id = COMID, X, Y) %>%
inner_join(poi_geom_xy, by = c("hy_id" = "hy_id", "X" = "X", "Y" = "Y")) %>%
inner_join(events_data, by = c("hy_id" = "hy_id", "X" = "X", "Y" = "Y"), multiple = "all") %>%
select(hy_id, REACHCODE, REACH_meas, event_id, poi_id) %>%
group_by(hy_id, REACHCODE) %>%
filter(REACH_meas == min(REACH_meas)) %>%
ungroup()
#distinct(hy_id, REACHCODE, REACH_meas, event_id, poi_id)
write_sf(select(events_data, -c(nexus, X, Y)), nav_gpkg, event_table)
write_sf(nexi, nav_gpkg, event_geometry_table)
# Load data
if(needs_layer(nav_gpkg, lookup_table_refactor)) {
full_cats <- readRDS(data_paths$fullcats_table)
lookup <- dplyr::select(sf::st_drop_geometry(nhd_Final),
NHDPlusV2_COMID = COMID,
lookup <- dplyr::select(sf::st_drop_geometry(nhd_Final),
NHDPlusV2_COMID = COMID,
realized_catchmentID = COMID,
mainstem = LevelPathI) %>%
dplyr::mutate(realized_catchmentID = ifelse(realized_catchmentID %in% full_cats$FEATUREID,
realized_catchmentID, NA)) %>%
left_join(nexus_from_poi(dplyr::rename(final_POIs, ID = COMID)), by = c("NHDPlusV2_COMID" = "ID"))
left_join(select(st_drop_geometry(poi_geometry), hy_id, poi_geom_id = geom_id),
by = c("NHDPlusV2_COMID" = "hy_id"))
sf::write_sf(lookup, nav_gpkg, lookup_table_refactor)
}
```
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