#' Network navigation for upstream/downstream from a COMID of interest
#' @param inCOM (list) list of input COMIDs
#' @param nhdDF (sf data.frame) (data frame) valid data frame of NHD flowlines
#' @param withTrib (logical) flag for if the upstream navigation should include tributaries
# or stick to mainstem level path
#
#' @return (list) list of COMIDs upstream of point
NetworkNav <- function(inCom, nhdDF, withTrib){
if (missing(withTrib)){
seg <- nhdplusTools::get_UM(nhdDF, inCom, include = TRUE)
} else {
seg <- nhdplusTools::get_UT(nhdDF, inCom)
}
return(seg)
}
#' Creates raw and dissolved segment layers with necessaary
# upstream/downstream routing attribution
#' @param nhdDF (sf data.frame) valid data frame of NHD flowlines
#' @param POIs (sf data.frame) Existing POIs
#'
#' @return (sf data.frame) data.frame of segments connecting POIs attributed
#' with POI_ID for each upstream flowpath
segment_increment <- function(nhdDF, POIs){
ptm<-proc.time()
seg_POIs <- arrange(POIs, desc(LevelPathI), desc(Hydroseq)) %>%
select(COMID, Hydroseq, LevelPathI) %>%
group_by(LevelPathI) %>%
# These next two levels arrange POIs hydrologically along the
# level path in order to meet the requirements of the code below
mutate(id = row_number(),
num = if(n() > 1) id[1L] + row_number()-1 else id) %>%
ungroup()
# Add an empty field for POI_Id population
nhdDF <- mutate(nhdDF, POI_ID = 0)
POI_ID_assign <- function(i, seg_POIs, nhd){
##########################################
# Populates POI_ID per segment
#
# Arguments:
# i : (integer) iterator
# seg_POIs : (data frame) POI data frame
# nhd : (dataframe) flowlines data frame (no geometry)
#
# Returns:
# nhd_sub : (data frame) raw segments with POI_ID populated
#
##########################################
library(dplyr)
# If POI is most upstream POI on levelpath
if (seg_POIs$num[i] == 1){
# Assign POI_ID
nhd_sub <- filter(nhd, Hydroseq >= seg_POIs$Hydroseq[i] &
LevelPathI == seg_POIs$LevelPathI[i]) %>%
mutate(POI_ID = seg_POIs$COMID[i])
# or as you travel downstream on set of POIs below level path
} else {
# Assign POI_ID
nhd_sub <- filter(nhd, LevelPathI == seg_POIs$LevelPathI[i] &
Hydroseq >= seg_POIs$Hydroseq[i] & Hydroseq < seg_POIs$Hydroseq[i-1]) %>%
mutate(POI_ID = seg_POIs$COMID[i])
}
# return assigned flowlines
return(select(nhd_sub, LevelPathI, Hydroseq, COMID, POI_ID) %>%
filter(POI_ID != 0))
}
library(parallel)
library(dplyr)
clust <- makeCluster(4)
POI_list <- parLapply(clust, c(1:nrow(seg_POIs)), POI_ID_assign, seg_POIs, st_drop_geometry(nhdDF))
#POI_list <- lapply(c(1:nrow(seg_POIs)), POI_ID_assign, seg_POIs, st_drop_geometry(nhdDF))
stopCluster(clust)
inc_segs <- data.table::rbindlist(POI_list)
print(proc.time()-ptm)
return(inc_segs)
}
#' Creates finalized segments and routing
#' @param nhdDF (sf data.frame) valid data frame of NHD flowlines
#' @param routing_fix (sf data.frame) any additional routing fixes
#'
#' @return (sf data.frame) data.frame of segments
segment_creation <- function(nhdDF, routing_fix = NULL){
if(!"StartFlag" %in% names(nhdDF)) {
nhdDF$StartFlag <- ifelse(nhdDF$Hydroseq %in% nhdDF$DnHydroseq, 0, 1)
}
in_segs <- filter(nhdDF, !is.na(POI_ID))
# If there are routing fixes to account for if a POI with a DA of 0 is moved upsream or downstream
if (is.data.frame(routing_fix)){
routing_fix <- routing_fix %>%
rename(COMID = oldPOI, new_COMID = COMID)
# Above we generated the network using the initial set of POIs; here we crosswalk over the old COMIDs to the new
nhd_fix <- nhdDF %>%
left_join(routing_fix %>%
select(COMID, new_COMID), by = c("POI_ID" = "COMID")) %>%
mutate(POI_ID = ifelse(is.na(new_COMID), POI_ID, new_COMID)) %>%
filter(!POI_ID %in% routing_fix$COMID) %>%
select(-new_COMID)
in_segs <- filter(nhd_fix, !is.na(POI_ID))
}
# Dissolve flowlines to aggregated segments
nsegments <- filter(in_segs, !is.na(POI_ID)) %>%
group_by(POI_ID) %>%
#arrange(desc(LevelPathI), desc(Hydroseq)) %>%
summarize(TotalLength = sum(LENGTHKM),TotalDA = max(TotDASqKM), HW = max(StartFlag),
do_union = FALSE) %>%
#st_cast("MULTILINESTRING") %>%
inner_join(st_drop_geometry(filter(in_segs, minNet == 1)) %>%
select(COMID, Hydroseq, DnHydroseq), by = c("POI_ID" = "COMID"))
# produce a short data frame for populating TO_POI for downstream segment
to_from <- filter(st_drop_geometry(in_segs)) %>%
left_join(filter(st_drop_geometry(nhdDF), !is.na(POI_ID)) %>%
select(COMID, Hydroseq, POI_ID), by = c("DnHydroseq" = "Hydroseq")) %>%
select(COMID.x, Hydroseq, DnHydroseq, POI_ID.y) %>%
rename(To_POI_ID = POI_ID.y)
# Add To_POI_ID to dissolved segments
nsegments_fin <- nsegments %>%
left_join(select(to_from, COMID.x, To_POI_ID), by = c("POI_ID" = "COMID.x")) %>%
select(POI_ID, TotalLength, TotalDA, HW, To_POI_ID)
return(list(diss_segs = nsegments_fin, raw_segs = in_segs))
}
#' Moves POI Upstream or downstream if it falls on COMID
# of flowline with no corresponding catchment
#' @param POIs_wgeom (sf data.frame) POIs
#' @param nhdDF (sf data.frame) valid data frame of NHD flowlines
#'
#' @return (sf data.frame) data.frame of POIs with new COMID associated
DS_poiFix <- function(POIs_wgeom, nhd){
nhd <- distinct(nhd)
POIs <- st_drop_geometry(POIs_wgeom) %>%
arrange(COMID) %>%
filter(nexus == FALSE)
# DF of downstream segment
tocomDF <- select(st_drop_geometry(nhd), COMID, Hydroseq, TotDASqKM,
DnHydroseq, WBAREACOMI) %>%
inner_join(select(st_drop_geometry(nhd), COMID_ds = COMID, Hydroseq,
WBAREACOMI_down = WBAREACOMI, totda_ds = TotDASqKM),
by = c("DnHydroseq" = "Hydroseq")) %>%
inner_join(select(st_drop_geometry(nhd), COMID_us = COMID, DnHydroseq,
WBAREACOMI_up = WBAREACOMI, totda_us = TotDASqKM),
by = c("Hydroseq" = "DnHydroseq"))
# Find segments with POIs where there is no corresponding catchment that are not terminal
unCon_fl <- filter(nhd, COMID %in% POIs$COMID, AreaSqKM == 0)# & Hydroseq != TerminalPa)
unCon_POIs <- filter(POIs, COMID %in% unCon_fl$COMID)
# Get specific fixes for waterbody inlets and outlets
wbout <- filter(unCon_POIs, !is.na(Type_WBOut)) %>%
inner_join(tocomDF, by = "COMID") %>%
mutate(nonrefactor = ifelse(WBAREACOMI %in% WBAREACOMI_up, COMID_ds, 0),
new_POI = COMID_us)
wb_pois <- filter(unCon_POIs, !is.na(Type_WBIn)) %>%
inner_join(tocomDF, by = "COMID") %>%
mutate(nonrefactor = ifelse(WBAREACOMI < 0, COMID_ds, 0),
new_POI = COMID_us) %>%
rbind(wbout) %>%
select(-c(nexus, Hydroseq, TotDASqKM, DnHydroseq, WBAREACOMI, COMID_ds,
WBAREACOMI_down, totda_ds, COMID_us, WBAREACOMI_up, totda_us)) %>%
rename(COMID = new_POI, oldPOI = COMID)
# The rest can be resolved with drainage are ratio
unCon_POIs <- filter(unCon_POIs, !COMID %in% wb_pois$oldPOI)
poi_fix <- as.data.frame(do.call("rbind", lapply(unCon_POIs$COMID, movePOI_NA_DA, st_drop_geometry(nhd)))) %>%
inner_join(POIs, by = c("oldPOI" = "COMID")) %>%
inner_join(select(st_drop_geometry(nhd), COMID), by = c("oldPOI" = "COMID")) %>%
select(-c(AreaSqKM, DnHydroseq, nexus, TotDASqKM)) %>%
distinct() %>%
bind_rows(wb_pois)
# Fold in new POIs with existing POIs so all the "Type" attribution will carry over
# using the minimum will ensure correct downstream hydrosequence gets carried over
poi_orig <- filter(POIs, COMID %in% poi_fix$COMID) %>%
bind_rows(poi_fix) %>%
select(-oldPOI)
list_df <- dplyr::group_by(poi_orig, COMID) |>
group_split()
compact <- function(l) {
if(nrow(l) == 1) return(as.data.frame(l))
lapply(names(l), \(x) {
out <- unique(l[[x]][!is.na(l[[x]])])
if(!length(out)) NA
else if(length(out) == 1) out
else {
cat(paste("duplicate id for", unique(l$COMID),
"column", x,
"values", paste(out, collapse = ", "),
"using", out[1]), file = "POI-issues.log")
out[1]
}
}) |>
setNames(names(l)) |>
as.data.frame()
}
poi_merged <- bind_rows(lapply(list_df, compact))
# # Combine POI information together for redundant pois
# poi_merged <- poi_orig %>%
# select(-c(nexus, AreaSqKM, oldPOI, DnHydroseq, TotDASqKM)) %>%
# group_by(COMID) %>%
# summarise_each(funs(toString(na.omit(.))))
# is.na(poi_merged) <- poi_merged == ""
# Join new POI COMIDs and geometry with the old Type fields
fin_POIs <- poi_merged %>%
arrange(COMID) %>%
bind_cols(get_node(filter(nhd, COMID %in% .$COMID) %>% arrange(COMID), position = "end")) %>%
st_sf() %>%
st_compatibalize(., POIs_wgeom)
return (list(xWalk = poi_fix, new_POIs = fin_POIs))
}
#' Move POIs that fall on flowlines with no catchment upstream/downstream
# to adjacent flowline with most similar total drainage area. Called from
# DS_poi_fix function above
#' @param poi_fix (data.frame) POI data set of COMIDs to be changed
#' @param nhdDF (sf data.frame) valid data frame of NHD flowlines
#'
#' @return (data frame) DF of POIs with new COMID associated
movePOI_NA_DA <- function(poi_fix, nhdDF){
#print(poi_fix)
nhdDF <- distinct(nhdDF)
# Closest POI/US/DS
up_segs <- unique(nhdplusTools::get_UM(nhdDF, poi_fix, sort=T))
seg2fix <- filter(nhdDF, COMID == poi_fix) %>%
distinct()
# Sorted results and filter out all flowlines w/o catchments
upstuff <- filter(nhdDF, COMID %in% unlist(up_segs)) %>%
arrange(Hydroseq) %>%
filter(AreaSqKM > 0)
down_segs <- unique(nhdplusTools::get_DM(nhdDF, poi_fix, sort=T))
downstuff <- filter(nhdDF, COMID %in% unlist(down_segs)) %>%
arrange(Hydroseq)%>%
filter(AreaSqKM > 0)
# combine into one dataframe, select up/downstream seg with least change in total drainage area
near_FL <- rbind(select(upstuff, COMID, TotDASqKM, AreaSqKM) %>%
slice(1),
select(downstuff, COMID, TotDASqKM, AreaSqKM) %>%
slice(1))
# If 1 or other adjacent flowlines are coupled with a catchment
if (sum(near_FL$AreaSqKM) > 0){
new_POI <- near_FL[(which.min(abs(seg2fix$TotDASqKM - near_FL$TotDASqKM))),] #near_FL$COMID
new_POI$oldPOI <- poi_fix
new_POI$DnHydroseq <-seg2fix$DnHydroseq
} else {
# Remove POI if not catchment associated with flowlines upstream or downstream
print (poi_fix)
print ("US and DS flowlines also have no catchment, removing POI")
new_POI <- NA
}
return(new_POI)
}
#' Collaposes/merges POIs together based on drainage area ratio and data theme
#' @param out_gpkg (geopackage) output geopackage to write intermediate results to
#' @param POIs (sf data.frame) Original POIs
#' @param Seg (sf data.frame) dissolved segments
#' @param Seg2 (sf data.frame) (data frame) raw segments/flowlines
#' @param exp (string) Type of thematic POI being moved round
#' @param DA_diff (float) threshold of drainage area difference to
# consider when merging two POIs
#'
#' @return (list of data.frames) 1 - New set of POIs
# 2/3 - correpsonding segments (both raw and dissolved)
POI_move_down<-function(out_gpkg, POIs, Seg, Seg2, exp, DA_diff){
nexus <- filter(POIs, nexus == TRUE)
POIs <- POIs %>% mutate_if(is.numeric, function(x) ifelse(is.infinite(x), NA, x)) %>%
filter(nexus != 1)
# downstream POIs and their drainage area ratios and incremental drainage areas
segs_down <- Seg %>%
inner_join(select(st_drop_geometry(.), c(POI_ID, TotalDA, TotalLength)),
by = c("To_POI_ID" = "POI_ID")) %>%
inner_join(st_drop_geometry(POIs), by = c("POI_ID" = "COMID")) %>%
# Keep WBOut and In to preserve NHD Waterbody geometry, Keep Confluences and TE where they are
filter_at(vars(Type_WBOut, Type_WBIn, Type_Conf, Type_TE), all_vars(is.na(.))) %>%
# Select POIs within the correct drainage area ratio and fit the size criteria
mutate(DAR = TotalDA.y/TotalDA.x, IncDA = TotalDA.y - TotalDA.x) %>%
# If the drainage area ratio is within acceptable range (1.05)
filter(DAR < (1 + DA_diff), POI_ID != To_POI_ID) %>%
# select and re-order
select(POI_ID, HW, To_POI_ID, Type_HUC12, Type_Gages, Type_TE, Type_NID,
Type_WBIn, Type_WBOut, Type_Conf, Type_Term, Type_Elev, Type_Travel, nexus, DAR, IncDA)
# upstream POIs and their drainage area ratios and incremental drainage areas
up_POIs <- filter(POIs, is.na(Type_Conf))
segs_up <- Seg %>% inner_join(select(filter(st_drop_geometry(.), POI_ID %in% up_POIs$COMID),
c(POI_ID, To_POI_ID, TotalDA, TotalLength)),
by = c("POI_ID" = "To_POI_ID")) %>%
inner_join(st_drop_geometry(POIs), by = c("POI_ID" = "COMID")) %>%
rename(From_POI_ID = POI_ID.y) %>%
# Don't want to move outlets upstream, don't move confluences upstream
filter_at(vars(Type_WBOut, Type_WBIn, Type_Conf, Type_TE), all_vars(is.na(.))) %>%
# Select POIs within the correct drainage area ratio and fit the size criteria
mutate(DAR = TotalDA.y/TotalDA.x, IncDA = TotalDA.x - TotalDA.y) %>%
# If the drainage area ratio is within acceptable range (0.95)
filter(DAR >= (1 - DA_diff), POI_ID != From_POI_ID) %>%
# select and re-order
select(POI_ID, HW, From_POI_ID, Type_HUC12, Type_Gages, Type_TE, Type_NID,
Type_WBIn, Type_WBOut, Type_Conf, Type_Term, Type_Elev, Type_Travel, nexus, DAR, IncDA)
# Filter by POI Type, this removes the POI theme
Types <- c("Type_Gages", "Type_HUC12", "Type_TE", "Type_Conf", "Type_NID", "Type_WBIn", "Type_WBOut",
"Type_Term", "Type_Elev", "Type_Travel")
Types <- Types[Types != exp]
# Gages coupled with confluences, waterbodies, do not move these
static_POIs <- POIs %>%
filter(!is.na(Type_Gages) & (!is.na(Type_Conf) | !is.na(Type_WBOut) | !is.na(Type_WBIn)))
# compiled list of POIs to move up or down
seg_choices <- filter(POIs, !COMID %in% static_POIs) %>%
left_join(st_drop_geometry(segs_down) %>%
select(POI_ID, To_POI_ID, DAR, IncDA), by = c("COMID" = "POI_ID")) %>%
left_join(st_drop_geometry(segs_up) %>%
select(POI_ID, From_POI_ID, DAR, IncDA), by = c("COMID" = "POI_ID")) %>%
filter(!is.na(To_POI_ID) | !is.na(From_POI_ID)) %>%
filter(IncDA.x < min_da_km | IncDA.y < min_da_km)
# We only want to move POIs that are not coupled within another theme
seg_sub <- seg_choices %>%
filter_at(Types, all_vars(is.na(.))) %>%
# Don't want to move the static POIS either
filter(!COMID %in% static_POIs$COMID) %>%
select (COMID, To_POI_ID, DAR.x, IncDA.x, From_POI_ID, DAR.y, DAR.y, IncDA.y) %>%
st_sf()
# Commented for QAQCing
#st_write(seg_sub, out_gpkg, "Seg_choices", delete_layer = F, append = T)
# Data frame to indicate whether POIs should be grouped upstream or downstream
segsub_dir <- mutate(seg_sub, DirDA = ifelse(is.na(DAR.y), "Down",
ifelse(is.na(DAR.x), "Up",
ifelse(IncDA.x < IncDA.y, "Down", "Up"))),
Move = ifelse(DirDA == "Down", To_POI_ID, From_POI_ID))
# POIs whose moves correspond to each others locations
Mult <- select(segsub_dir, -c(To_POI_ID, From_POI_ID)) %>% #select(segsub_dir, COMID, Move) %>%
inner_join(st_drop_geometry(.), by = c("Move" = "COMID"), suffix = c(".A", ".B"))
# Basicly just chose one of the pair in Mult.
segsub_dir <- filter(segsub_dir, !COMID %in% Mult$Move)
#1 - POIs that need to be moved downstream
move_POI <- filter(POIs, COMID %in% segsub_dir$COMID)
# POIs that are stationary and their infomration will be merged with upstream POIs
stationary_POI <- filter(POIs, !COMID %in% move_POI$COMID)
# Final Set to be merged with that don't involve either sets above
#FinalPOI <- POIs %>% filter(!COMID %in% SegSub_Dir$COMID & !COMID %in% SegSub_Dir$Move) %>%
# mutate(merged_COMID = NA)
#2 - Join SegSub assignment to POIs to bring over POI attributes
movedownPOI_withatt <- select(st_drop_geometry(segsub_dir), COMID, DirDA, Move) %>%
left_join(st_drop_geometry(move_POI), by = "COMID")
# Join SegSub mod to downstream POIs, COMID.x is the final COMID
merged_POIs <- stationary_POI %>%
inner_join(movedownPOI_withatt,
by = c("COMID" = "Move"), suffix = c(".x", ".y")) %>%
#merged_POIs[mapply(is.infinite, merged_POIs)] <- NA
#merged_POIs <- merged_POIs %>%
#mutate_all(., list(~na_if(.,""))) %>%
# Don't overwrite existing gages or HUC12s
filter(is.na(!!as.symbol(paste0(exp, ".x")))) %>%
# Bring over relevant attributes
mutate(Type_HUC12 = ifelse(is.na(Type_HUC12.y), Type_HUC12.x, Type_HUC12.y)) %>%
mutate(Type_Gages = ifelse(is.na(Type_Gages.x), Type_Gages.y, Type_Gages.x)) %>%
mutate(Type_WBIn = ifelse(is.na(Type_WBIn.y), Type_WBIn.x, Type_WBIn.y)) %>%
mutate(Type_WBOut = ifelse(is.na(Type_WBOut.y), Type_WBOut.x, Type_WBOut.y)) %>%
mutate(Type_TE = ifelse(is.na(Type_TE.y), Type_TE.x, Type_TE.y)) %>%
mutate(Type_NID = ifelse(is.na(Type_NID.y), Type_NID.x, Type_NID.y)) %>%
mutate(Type_Conf = ifelse(is.na(Type_Conf.y), Type_Conf.x, Type_Conf.y)) %>%
mutate(Type_Term = ifelse(is.na(Type_Term.y), Type_Term.x, Type_Term.y)) %>%
mutate(Type_Elev = ifelse(is.na(Type_Elev.y), Type_Elev.x, Type_Elev.y)) %>%
mutate(Type_Travel = ifelse(is.na(Type_Travel.y), Type_Travel.x, Type_Travel.y)) %>%
mutate(oldPOI = COMID.y, to_com = NA) %>%
select(COMID, Type_HUC12, Type_Gages, Type_WBIn, Type_WBOut,
Type_TE, Type_NID, Type_Conf, Type_Term,
Type_Elev, Type_Travel, DirDA, nexus = nexus.x, oldPOI = COMID.y) %>%
group_by(COMID) %>%
summarize_all(~paste(unique(na.omit(.)), collapse = ',')) %>%
st_sf()
print(colnames(POIs))
print(colnames(merged_POIs))
# Bring new POI data over to new data
fin_POIs <- filter(POIs, !COMID %in% c(merged_POIs$oldPOI, merged_POIs$COMID)) %>%
rbind(merged_POIs %>% select(-c(DirDA, oldPOI))) %>%
rbind(nexus)
#inner_join(st_drop_geometry(merged_POIs) %>% select(COMID, Type_Gages_A), by = "COMID") %>%
#mutate(Type_Gages = ifelse(!is.na(Type_Gages_A), Type_Gages_A, Type_Gages)) %>% select(-Type_Gages_A)
changed_POIs <- POIs %>%
inner_join(select(st_drop_geometry(merged_POIs), COMID, oldPOI)) %>%
select(COMID, oldPOI, Type_HUC12, Type_Gages, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf,
Type_Term, Type_Elev, Type_Travel) %>%
mutate(to_com = COMID) %>%
st_compatibalize(., read_sf(out_gpkg, poi_moved_layer))
st_write(changed_POIs, out_gpkg, poi_moved_layer, append = TRUE) # write_sf not appending?
# Format POI moves to table
xWalk <- select(st_drop_geometry(segsub_dir), Move, COMID) %>%
filter(!is.na(Move), Move %in% merged_POIs$COMID) %>%
rename(COMID = Move, oldPOI = COMID)
st_write(xWalk, out_gpkg, "collapse_xWalk", append = TRUE) # write_sf not appending?
newSegs <- segment_creation(mutate(Seg2, old_POI_ID = POI_ID), xWalk)
newSegs$raw_segs <- distinct(newSegs$raw_segs, .keep_all = T)
# Return POIs, segments, and raw segments
return (list(allPOIs = fin_POIs, segs = newSegs$diss_segs, FL = newSegs$raw_segs))
}
#' Identifies and attributes structural POIs
#' @param nhdDF (sf data.frame) valid data frame of NHD flowlines
#' @param final_POIs (sf data.frame) final POIs
#'
#' @return (data.frame columns) 1/0 Columns indicating structural POIs and the structural network
structPOIsNet <- function(nhdDF, final_POIs){
# subset to flowlines used for nsegments
inc_segs <- filter(nhdDF, !is.na(POI_ID))
vpu <- inc_segs$VPUID
# Terminal outlets from the initial network
termout <- filter(nhdDF, !Hydroseq %in% DnHydroseq & TerminalFl == 1 & !COMID %in% final_POIs$COMID)
# POIs that are also terminal outlets
out_POIs <- filter(nhdDF, COMID %in% final_POIs$COMID & TerminalFl == 1)
# Confluence POIs
conf_POIs <- filter(inc_segs, COMID %in% final_POIs$COMID[!is.na(final_POIs$Type_Conf)])
# Waterbody outlet POIs
wb_POIs <- filter(inc_segs, COMID %in% final_POIs$COMID[!is.na(final_POIs$Type_WBOut) | !is.na(final_POIs$Type_WBIn)])
# Waterbody inlets POIs
struc_flines <- termout %>%
bind_rows(out_POIs, conf_POIs, wb_POIs) %>%
arrange(COMID)
struc_flines <- struc_flines[!st_is_empty(struc_flines), , drop = F] %>%
st_cast("LINESTRING")
struc_POIs <- get_node(struc_flines, position = "end") %>%
mutate(COMID = struc_flines$COMID, LevelPathI = struc_flines$LevelPathI)
# Add back in any levelpaths missing (shouldn't be much)
miss_LP <- filter(inc_segs, COMID %in% final_POIs$COMID) %>%
filter(!LevelPathI %in% struc_POIs$LevelPathI)
# Only bind if there are rows present
if (nrow(miss_LP) > 0){
lp_POIs <- c(miss_LP$LevelPathI, struc_POIs$LevelPathI)
} else {
lp_POIs <- struc_POIs$LevelPathI
}
# final stuctural POIs, order by COMID to match with strucPOIs
final_struct <- filter(nhdDF, Hydroseq %in% lp_POIs) %>%
arrange(COMID)
struct_POIs <- get_node(final_struct, position = "end") %>%
mutate(COMID = final_struct$COMID)
# produce unique set of flowlines linking POIs
final_net <- unique(unlist(lapply(unique(final_struct$COMID), NetworkNav, st_drop_geometry(nhdDF))))
# Subset NHDPlusV2 flowlines to navigation results and write to shapefile
structnet <- filter(nhdDF, COMID %in% final_net & grepl(paste0("^", vpu, ".*"), .data$VPUID))
# Write out minimally-sufficient network and POIs as a list
return(list(struc_POIs = struc_POIs, structnet = structnet))
}
#' Splits a network based elevation change and size
#' @param in_POI_ID (integer) POI_ID of aggregated flowline that needs to be
#' split based on elevation
#' @param split_DF (sf data.frame) flowlines attributed with POI_ID
#' @param elev_diff (numeric) max elevation change threshold within a segment
#' @param max_DA (numeric) minimum drainage area for resulting split
#'
#' @return (sf data.frame) flowlines with new POI_IDs identified (elev_POI_ID)
split_elev <- function(in_POI_ID, split_DF, elev_diff, max_DA){
# subset to a given POI_ID
split_elev_DF <- filter(split_DF, POI_ID == in_POI_ID) %>%
filter(AreaSqKM > 0) %>%
mutate(inc_seg_area = c(TotDASqKM[1], (TotDASqKM - lag(TotDASqKM))[-1]))
first_iter <- unique(split_elev_DF$iter)
# Iterate through segs that need splitting
for (i in c(first_iter:(max(split_DF$iter, na.rm = T) + 10))){
#print(i)
# first split
elev_pois_init_iter <- split_elev_DF %>%
# filter by iteration, comes into play when multiple splits
# are required per single POI_ID
filter(iter == i) %>%
# Recalculate inc elev, length, and area based on the last split
mutate(csum_length = cumsum(LENGTHKM),
inc_elev = cumsum(MAXELEVSMO - MINELEVSMO),
sum_area = cumsum(inc_seg_area)) %>%
# Determine if split actually necessary on this iteration
filter((max(MAXELEVSMO) - min(MINELEVSMO)) > (elev_diff), sum_area > max_DA) %>%
# This is an important step, identify which row to split on based on incremental
# elevaton value. Some flowlines exceed the elevation difference threshold
# The 'split_index' values determiens which row to create a new POI on
mutate(split_index = ifelse(nrow(filter(., inc_elev > (elev_diff))) == nrow(.),
1, nrow(filter(., inc_elev < (elev_diff))))) %>%
# Take out if at last iteration to avoid duplicating existng POIs
filter(!COMID == POI_ID)
if(nrow(elev_pois_init_iter) == 0){
#print("done")
# Done iterating on possible splits
return(split_elev_DF)}
# Get the flowline POI
elev_pois_init_pre <- elev_pois_init_iter %>%
filter(row_number() == split_index) %>%
mutate(new_POI_ID = COMID) %>%
select(COMID, POI_ID, new_POI_ID, Hydroseq_cut = Hydroseq)
# Remaining rows downstream after split
elev_pois_init_post <- elev_pois_init_iter %>%
left_join(select(elev_pois_init_pre, -COMID), by = "POI_ID") %>%
filter(Hydroseq < Hydroseq_cut) %>%
ungroup()
# Test for if we need to split again
if(nrow(elev_pois_init_post) > 0){
#print("yope")
# Re-cacl the iter based on results above
split_elev_DF <- split_elev_DF %>%
left_join(select(elev_pois_init_pre, -COMID), by = "POI_ID") %>%
mutate(elev_POI_ID = ifelse((is.na(elev_POI_ID) & Hydroseq >= Hydroseq_cut), new_POI_ID, elev_POI_ID),
iter = ifelse((!is.na(new_POI_ID) & Hydroseq < Hydroseq_cut), iter + 1, orig_iter)) %>%
select(-c(new_POI_ID, Hydroseq_cut)) %>%
ungroup()
} else {
split_elev_DF <- split_elev_DF %>%
left_join(select(elev_pois_init_pre, -COMID), by = "POI_ID") %>%
mutate(elev_POI_ID = ifelse((is.na(elev_POI_ID) & Hydroseq >= Hydroseq_cut), new_POI_ID, elev_POI_ID),
iter = ifelse((!is.na(new_POI_ID) & Hydroseq < Hydroseq_cut), 0, orig_iter)) %>%
select(-c(new_POI_ID, Hydroseq_cut)) %>%
ungroup()
}
}
return(split_elev_DF)
}
#' Retrieves waterbodies from NHDArea layer or NHD Hi-Res for ResOpsUs
#' locations if absent from NHD waterbody
#' @param nhd (sf data.frame) Fowlines for given VPU
#' @param WBs (sf data.frame) waterbodiess for discretization within VPU
#' @param data_paths (list) data paths to data layers
#'
#' @return (list) wbs - sf data frame for NHDArea and HR waterbodies
#' wb_table - table of flowlines and outlet info for each
#' feature in wb
HR_Area_coms <- function(nhd, WBs, data_paths, crs){
# Pull out rows for VPU that are NHDArea
nhd_area_resops <- WBs %>%
filter(resops_flowlcomid %in% nhd$COMID, source == "NHDAREA")
# Pull out rows for VPU that are NHD HR
nhd_hr_wb_resops <- WBs %>%
filter(member_comid != 65000300139130) %>% # R09, in Canada
filter(resops_flowlcomid %in% nhd$COMID, source == "HR ID AVAILABLE") %>%
st_drop_geometry()
# Get reachcodes for waterbody outlets, so we have an idea of which
# NHD HR 4-digit geodatabase we may need to retrieve
RC <- filter(nhd, COMID %in% nhd_hr_wb_resops$resops_flowlcomid)$REACHCODE
# If no NHDArea or HR waterbodies needed return NULL
if (nrow(nhd_area_resops) == 0 & nrow(nhd_hr_wb_resops) == 0){
return(NA)
}
# If NHDArea feature needed retrieve from National GDB
if (nrow(nhd_area_resops) > 0){
nhd_area_vpu <- read_sf(data_paths$nhdplus_gdb, "NHDArea") %>%
filter(COMID %in% nhd_area_resops$member_comid) %>%
mutate(source = "NHDv2Area")
wb <- st_transform(nhd_area_vpu, crs)
}
# If NHDHR feature needed
if (nrow(nhd_hr_wb_resops) > 0){
# HUC04 we need to download
huc04 <- substr(RC, 1, 4)
# Download NHD HR HUC04 if we dont' have it, other wise load and
# Bind NHDHR waterbodies into one layer
hr_wb <- do.call("rbind", lapply(unique(huc04), function(x){
print(x)
vpu <- unique(substr(x, 1, 2))
if(!file.exists(file.path(data_paths$nhdplus_dir, vpu,
paste0("NHDPLUS_H_", x, "_HU4_GDB.gdb")))){
download_nhdplushr(data_paths$nhdplus_dir, unique(huc04))
}
gdb <- list.files(file.path(data_paths$nhdplus_dir, substr(vpu, 1, 2)),
pattern = paste0("_", x, "_.+gdb"), full.names = TRUE)
# Format to similar to NHDArea/Waterbody layers
read_sf(gdb, "NHDWaterbody")})) #%>%
# convert to lower case and designate new shape field
names(hr_wb) <- tolower(names(hr_wb))
st_geometry(hr_wb) <- "shape"
hr_wb <- filter(hr_wb, permanent_identifier %in% WBs$member_comid) %>%
rename(GNIS_NAME = gnis_name, GNIS_ID = gnis_id,
Permanent_Identifier = permanent_identifier, AREASQKM = areasqkm,
FTYPE = ftype, FCODE = fcode, FDATE = fdate,
REACHCODE = reachcode) %>%
#select(-c(permanent_identifier, visibilityfilter, vpuid)) %>%
st_zm(.) %>%
st_as_sf() %>%
mutate(source = "NHDHR",
wb_id = ifelse(!is.na(GNIS_ID), GNIS_ID, Permanent_Identifier)) %>%
st_transform(crs)
# Bind or create new object
if(exists("wb")){
hr_wb <- st_compatibalize(hr_wb, wb)
wb <- data.table::rbindlist(list(wb, hr_wb), fill = TRUE) %>%
st_as_sf()
} else {
wb <- hr_wb
}
}
# # get the outlt rows from the table
# resops_outlet <- read.csv("data/reservoir_data/ISTARF-CONUS.csv") %>%
# filter(WBAREACOMI %in% wb$COMID | HR_NHDPLUSID %in% wb$COMID)
# Clear out columns not needed
WBs_sub <- WBs %>%
select(-any_of(c("COMID", "GNIS_ID")))
WBs_fin <- st_drop_geometry(WBs_sub) %>%
mutate(member_comid = member_comid) %>%
inner_join(select(wb, Permanent_Identifier, GNIS_ID, GNIS_NAME2 = GNIS_NAME),
by = c("member_comid" = "Permanent_Identifier")) %>%
mutate(wb_id = ifelse(GNIS_ID %in% c(NA, " "), member_comid, GNIS_ID),
GNIS_NAME = ifelse(is.na(GNIS_NAME), GNIS_NAME2, GNIS_NAME)) %>%
select(-c(GNIS_ID, GNIS_NAME2)) %>%
st_as_sf()
# Create table of all flowlines that intersect the waterbody
nhd_wb <- st_intersects(st_transform(nhd, st_crs(wb)), wb)
comid <- nhd[lengths(nhd_wb) > 0,]$COMID
nhd_wb_all <- nhd_wb[lengths(nhd_wb) > 0] %>%
purrr::set_names(comid) %>%
stack() %>%
# Ind is the default name for the set_names
rename(comid = ind, nhd_wb = values) %>%
mutate(wb_comid = wb[nhd_wb,]$Permanent_Identifier,
outlet = ifelse(comid %in% WBs_fin$resops_flowlcomid, "outlet", NA),
comid = as.integer(as.character(comid))) %>%
#left_join(select(wb_table, DAM_ID, DAM_NAME, resops_FlowLcomid),
# by = c("comid" = "resops_FlowLcomid")) %>%
left_join(select(st_drop_geometry(wb), wb_id, Permanent_Identifier, GNIS_ID,
source),
by = c("wb_comid" = "Permanent_Identifier"))
return(list(nhd_wb_table = nhd_wb_all, wb = WBs_fin))
}
#' Creates wb inlet and outlet events for splitting in hyRefactor
#' for waterbodies derived from NHDArea and NHDHR waterbodies
#' @param WBs (sf data.frame) return from HR_Area_coms
#' @param nhd_wb (sf data.frame) flowlines that intersect waterbodies
#' @param type (character) whether to derive inlet or outlet points
#'
#' @return (sf data.frame) dataframe of WB inlet and outlet points to split
#'
WB_event <- function(WBs, nhd_wb, type){
# split into features and table
WBs_table <- WBs$nhd_wb_table
WBs_layer <- WBs$wb
if (type == "outlet"){
# get the outlet comid from the ResOps Table
outlet_fl <- filter(nhd_wb, COMID %in% filter(WBs_table, outlet == "outlet")$comid)
# Get the downstream flowline for continuity
ds_fl <- filter(nhd_wb, Hydroseq %in% outlet_fl$DnHydroseq,
LevelPathI %in% outlet_fl$LevelPathI) %>%
rbind(outlet_fl) %>%
arrange(desc(Hydroseq)) %>%
group_by(LevelPathI) %>%
# union together
summarize(do_union = TRUE) %>%
st_cast("LINESTRING")
WBs_HR <- filter(WBs_layer, source == "HR ID AVAILABLE")
# For NHD HR waterbody outlets its a bit more complicated
if (nrow(WBs_HR) > 0){
# Get the flowlines intersecting the HR waterbody and find one with the
# max DA
outlet_wb_int <- nhd_wb[lengths(st_intersects(st_transform(nhd_wb, st_crs(WBs_HR)),
WBs_HR)) > 0,] %>%
group_by(wb_id) %>%
filter(TotDASqKM == max(TotDASqKM)) %>%
ungroup()
# get the ds flo with the same levepath (JIC)
ds_fl <- filter(nhd_wb, DnHydroseq %in% outlet_wb_int$Hydroseq,
LevelPathI %in% outlet_wb_int$LevelPathI)
outlet_fl <- rbind(outlet_wb_int, ds_fl)
# Cast flowlines within NHDHR waterbody to point
WB_FL_pnts <- outlet_wb_int %>%
st_cast("POINT") %>%
group_by(wb_id) %>%
mutate(pnt_id = row_number()) %>%
ungroup()
outlet_pnts <- do.call("rbind", lapply(unique(WB_FL_pnts$wb_id), function(x){
fl <- filter(WB_FL_pnts, wb_id == x)
wb <- filter(WBs_HR, wb_id == x)
# Determine which points intersect waterbody
WB_outlet_pnts <- fl[lengths(st_intersects(fl, wb)) > 0,] %>%
st_drop_geometry() %>%
group_by(wb_id) %>%
mutate(within_wb_id = row_number()) %>%
filter(within_wb_id >= max(within_wb_id)) %>%
ungroup() %>%
select(wb_id, orig_pnt_id = pnt_id, within_wb_id)
# Deriv new linestring by concating points from most upstream point
# within waterbody to downstream so we can split at FL/waterbody
# nexus
outlet_FL <- fl %>%
inner_join(WB_outlet_pnts, by = "wb_id") %>%
select(wb_id, pnt_id, orig_pnt_id, within_wb_id) %>%
filter(pnt_id >= orig_pnt_id) %>%
group_by(wb_id) %>%
summarize(do_union = F) %>%
st_cast("LINESTRING") %>%
filter(wb_id %in% wb$wb_id)
outlet_pnt <- sf::st_intersection(outlet_FL, WBs_HR) %>%
st_cast("MULTIPOINT") %>%
st_cast("POINT") %>% # was point
group_by(wb_id) %>%
filter(row_number() == max(row_number(), na.rm = T)) %>%
ungroup() %>%
mutate(id = row_number()) %>%
filter(wb_id == wb_id.1)
}))
# Derive the events
if(exists("wb_events")){
hr_events <- get_flowline_index(st_transform(nhd_wb, st_crs(outlet_pnts)),
outlet_pnts) %>%
inner_join(select(st_drop_geometry(nhd_wb), COMID, wb_id), by = "COMID") %>%
filter(wb_id %in% WBs_HR$wb_id) %>%
mutate(event_type = type) %>%
cbind(select(outlet_pnts, geom)) %>%
st_as_sf() %>%
st_transform(st_crs(wb_events)) %>%
st_compatibalize(., wb_events)
wb_events <- rbind(wb_events, hr_events) %>%
distinct()
} else {
wb_events <- do.call("rbind", lapply(unique(outlet_pnts$wb_id), function(x){
outlet_pnt <- filter(outlet_pnts, wb_id == x)
get_flowline_index(nhd_wb, outlet_pnt)})) %>%
left_join(distinct(st_drop_geometry(outlet_pnts), resops_flowlcomid, wb_id),
by = c("COMID" = "resops_flowlcomid")) %>%
mutate(event_type = type) %>%
cbind(select(outlet_pnts, geom)) %>%
st_as_sf()
}
}
# For inlet points its alot easier for both NHDARea and NHDHR
} else {
start_pts <- get_node(nhd_wb, position = "start") %>%
cbind(st_drop_geometry(nhd_wb))
inlet_FL <- nhd_wb[lengths(st_intersects(
st_transform(start_pts, st_crs(WBs_layer)), WBs_layer)) == 0,] %>%
rbind(filter(nhd_wb, Hydroseq %in% .$DnHydroseq,
LevelPathI %in% .$LevelPathI)) %>%
arrange(desc(Hydroseq)) %>%
unique()
inlet_ls <- inlet_FL %>%
group_by(LevelPathI) %>%
st_cast("POINT") %>%
summarize(do_union = F) %>%
st_cast("LINESTRING") %>%
ungroup()
inlet_pnts_int <- sf::st_intersection(
st_transform(inlet_ls, st_crs(WBs_layer)), WBs_layer) #%>%
#st_cast("LINESTRING")
int_point <- inlet_pnts_int[st_geometry_type(inlet_pnts_int) %in% c("POINT"),]
if(nrow(int_point) > 0){
inlet_fl <- inlet_pnts_int[st_geometry_type(inlet_pnts_int) %in% c("LINESTRING", "MULTILINESTRING"),] %>%
st_cast("LINESTRING")
inlet_pnts_int <- rbind(int_point, inlet_fl)
}
inlet_pnts <- inlet_pnts_int %>%
group_by(LevelPathI) %>%
st_cast("POINT") %>%
mutate(id = row_number()) %>%
filter(id == min(id), !is.na(wb_id)) %>%
ungroup()
wb_events <- get_flowline_index(st_transform(nhd_wb, st_crs(inlet_pnts)),
inlet_pnts) %>%
inner_join(select(st_drop_geometry(nhd_wb), COMID, wb_id, LevelPathI), by = "COMID") %>%
mutate(event_type = type) %>%
inner_join(select(inlet_pnts, LevelPathI), by = "LevelPathI") %>%
#group_by(COMID, LEVELPATHI) %>%
mutate(nexus = T) %>%
rename(POI_identifier = wb_id) %>%
select(-c(id, offset)) %>%
st_as_sf()
}
return(wb_events)
}
#' Creates POIs for gages using refactor criteria
#' @param temp_POIs (sf data.frame) current data frame of POIs
#' @param gages_info (sf data.frame) VPU specific streamgages from 01_gage_selection
#' @param nhd (sf data.frame) flowline data.frame
#' @param combine_meters (integer) refactor threshold (m) for if if two adjacent fl should be combined
#' @reach_meas_thresh (integer) threshold added or substracted from reach_meas to determine if a gage will split fl
#'
#' @return (sf data.frame) dataframe of gage POIs
#'
gage_POI_creation <- function(tmp_POIs, gages_info, nhd, combine_meters, reach_meas_thresh){
# Create streamgage POIs
gage_POIs <- POI_creation(select(st_drop_geometry(gages_info), COMID, site_no), nhd, "Gages") %>%
st_compatibalize(., tmp_POIs)
# Avoid these for refactoring
avoid <- dplyr::filter(nhd, (sqrt(AreaSqKM) / LENGTHKM) > 3 & AreaSqKM > 1)
# Create events for streamgages for flowline splitting
events <- gages_info %>%
# bring over NHD information
inner_join(select(st_drop_geometry(gage_POIs), Type_Gages), by = c("site_no" = "Type_Gages")) %>%
# Get over the NHD attributes
#inner_join(select(st_drop_geometry(nhd), AreaSqKM, REACHCODE, LENGTHKM, COMID, FromMeas, ToMeas),
# by = "COMID") %>%
# Apply reach measure thresholds
filter(reach_meas - FromMeas > reach_meas_thresh & AreaSqKM > 2 &
ToMeas - reach_meas > reach_meas_thresh & LENGTHKM > (combine_meters / m_per_km)) %>%
select(COMID, REACHCODE = reachcode, REACH_meas = reach_meas, Type_Gages = site_no) %>%
filter(!COMID %in% avoid$COMID) %>%
mutate(event_type = "streamgage") %>%
st_compatibalize(., tmp_POIs) %>%
mutate(nexus = TRUE)
# Reset gage flag if even is created
gage_POIs_nonevent <- filter(gage_POIs, !Type_Gages %in% events$Type_Gages) %>%
addType(., tmp_POIs, "Gages", nexus = FALSE, bind = TRUE)
if(nrow(events) > 0){
tmp_POIs <- data.table::rbindlist(list(gage_POIs_nonevent,
select(events, COMID, Type_Gages, nexus)), fill = TRUE) %>%
mutate(nexus = ifelse(is.na(nexus), FALSE, nexus)) %>%
st_as_sf()
} else {
events <- NA
tmp_POIs <- mutate(tmp_POIs, nexus = FALSE)
}
return(list(events = events, tmp_POIs = tmp_POIs))
}
#' Creates Waterbody POIs, calls a few other functions
#' @param nhd (sf data.frame) flowline data.frame
#' @param WBs_VPU (sf data.frame)
#' @param data_paths (list) data paths to data layers
#' @param crs (integer) CRS to use (epsg code)
#'
#' @return (sf data.frame) dataframe of waterbody outlet POIs
#'
wbout_POI_creaton <- function(nhd, wb_table, data_paths, tmp_POIs, crs){
wb_components_table <- st_drop_geometry(wb_table) %>%
dplyr::mutate(member_comid = strsplit(member_comid, ",")) %>%
tidyr::unnest(cols = member_comid) %>%
mutate(member_comid = as.integer(member_comid)) %>%
distinct()
# R17, some wierd stuff going on in the waterbodies
if("947070203" %in% wb_components_table$member_comid){
wb_components_table <- wb_components_table %>%
filter(!wb_id %in% c(1503847, 1513298))
}
# Bring wb_id to NHD data frame
if(("wb_id") %in% colnames(nhd)){
nhd <- select(nhd, -wb_id)
}
nhd <- nhd %>%
left_join(distinct(st_drop_geometry(wb_components_table), wb_id, member_comid) %>%
mutate(member_comid = as.numeric(member_comid)),
by = c("WBAREACOMI" = "member_comid"))
# Resops that are either NHDArea or MHD HR source and are missing in the
# nhd attributiong table
resops_wb_other_sources <- filter(wb_table, source ==
"HR ID AVAILABLE" | (source == "NHDAREA" &
!wb_id %in% nhd$wb_id) &
accounted == 0)
# Resops that are not accounted for and NHDv2 waterbody or NHDArea sourced
resops_main <- filter(wb_table, !is.na(resops_flowlcomid) &
!wb_id %in% resops_wb_other_sources$wb_id,
accounted == 0)
# Remaining WB that don't have outlet COMID defined and aren't in ResOPsUS
WB_outlet <- filter(wb_table, is.na(resops_flowlcomid),
!wb_id %in% tmp_POIs$Type_WBOut)
# Size-based waterbody outlet POIs ot in NHDArea, HR, or
# manually defined with ResOPS
wbout_COMIDs <- nhd %>%
filter(wb_id %in% WB_outlet$wb_id) %>%
group_by(wb_id) %>%
slice(which.min(Hydroseq)) %>%
ungroup() %>%
switchDiv(., nhd) %>%
select(COMID, wb_id) %>%
mutate(wbtype = "size") %>%
st_drop_geometry() %>%
distinct()
# ResOpsUS defined outlets
resops_main_outlets <- select(st_drop_geometry(resops_main), COMID = resops_flowlcomid,
wb_id) %>%
mutate(wbtype = "resops")
# combine together
wbout_coms <- distinct(rbind(wbout_COMIDs, resops_main_outlets))
# Get NHDArea and HI-Res waterbodies for ResOpsUS locaitons
WBs_VPU_areaHR <- HR_Area_coms(nhd, resops_wb_other_sources, data_paths, crs)
# if there are no HR/Area waterbodies generated
if(all(is.na(WBs_VPU_areaHR))){
tmp_POIs <- POI_creation(wbout_coms,
nhd, "wb_id") %>%
addType(., tmp_POIs, "wb_id", nexus = TRUE) %>%
mutate(Type_WBOut = ifelse(is.na(Type_WBOut), Type_wb_id, Type_WBOut)) %>%
select(-Type_wb_id) %>%
left_join(select(st_drop_geometry(resops_main), GRAND_ID, resops_FlowLcomid),
by = c("COMID" = "resops_FlowLcomid")) %>%
mutate(Type_resops = ifelse(!is.na(GRAND_ID), GRAND_ID, Type_resops)) %>%
select(-GRAND_ID)
} else {
# WBarea out comids
wbareaout_coms <- select(st_drop_geometry(nhd), COMID, Hydroseq) %>%
inner_join(filter(WBs_VPU_areaHR$nhd_wb_table, !is.na(outlet)),
by = c("COMID" = "comid")) %>%
group_by(wb_id) %>%
filter(Hydroseq == min(Hydroseq)) %>% #, row_number() == 1
ungroup() %>%
inner_join(select(WBs_VPU_areaHR$wb, -source), by = "wb_id") %>%
filter(source != "NHDHR") %>%
select(COMID, wb_id) %>%
mutate(wbtype = "resops")
if(nrow(wbareaout_coms) > 0){
wbout_coms <- rbind(wbout_coms, wbareaout_coms)
}
tmp_POIs <- POI_creation(select(wbout_coms, COMID, wb_id),
nhd, "wb_id") %>%
addType(., tmp_POIs, "wb_id", nexus = TRUE) %>%
mutate(Type_WBOut = ifelse(is.na(Type_WBOut) & (!nexus | is.na(nexus)),
Type_wb_id, Type_WBOut)) %>%
select(-Type_wb_id) %>%
left_join(select(st_drop_geometry(resops_main), grand_id, resops_flowlcomid),
by = c("COMID" = "resops_flowlcomid")) %>%
mutate(Type_resops = ifelse(!is.na(grand_id) & (!nexus | is.na(nexus)),
grand_id, Type_resops)) %>%
select(-grand_id)
}
# If the return is not NA
if(inherits(WBs_VPU_areaHR, "list")){
# Attribute flowlines that intersect NHDARea and/or NHD HR waterbodies with
# waterbody COMIDs
# For NHDHR, two waterbodies might be back to back and share an intersecting
# flowline, this may cause problems below, so deal with this case
WBs_VPU_areaHR$nhd_wb_table <- WBs_VPU_areaHR$nhd_wb_table %>%
left_join(select(st_drop_geometry(resops_wb_other_sources),
resops_flowlcomid, member_comid, grand_id),
by = c("comid" = "resops_flowlcomid",
"wb_comid" = "member_comid")) %>%
group_by(comid) %>%
filter(case_when(n() == 2 ~ !is.na(grand_id),
TRUE ~ n() == 1)) %>%
ungroup()
nhd <- nhd %>%
left_join(distinct(WBs_VPU_areaHR$nhd_wb_table, comid, wb_id_hr = wb_id),
by = c("COMID" = "comid")) %>%
mutate(wb_id = ifelse(!is.na(wb_id_hr), wb_id_hr, wb_id)) %>%
select(-wb_id_hr)
# NHD flowlines within the HR waterbodies
nhd_WBs <- filter(nhd, wb_id %in%
filter(WBs_VPU_areaHR$nhd_wb_table, source == "NHDHR")$wb_id)
if(nrow(nhd_WBs) > 0){
# Create outlet events for splitting HR and NHDArea waterbodies
wb_outlet_events <- WB_event(WBs_VPU_areaHR, nhd_WBs, "outlet") %>%
st_compatibalize(., tmp_POIs) %>%
mutate(nexus = TRUE, Type_resops = 2)
tmp_POIs <- data.table::rbindlist(list(tmp_POIs,
select(wb_outlet_events, COMID, Type_WBOut = wb_id, Type_resops, nexus)), fill = TRUE) %>%
st_as_sf()
} else {
wb_outlet_events <- NA
}
WBs_VPU_areaHR$wb <- WBs_VPU_areaHR$wb %>%
sfheaders::sf_remove_holes(.) %>%
st_compatibalize(wb_poi_lst)
WBs_VPU_areaHR$wb <- nhdplusTools::st_compatibalize(
WBs_VPU_areaHR$wb, wb_table)
wb_table_filtered <- wb_table %>%
filter(!st_is_empty(wb_table))
# Append NHDArea and NHDHR waterbodies to waterbody layer
WBs_VPU <- data.table::rbindlist(list(wb_table_filtered, WBs_VPU_areaHR$wb),
fill = TRUE) %>%
st_as_sf()
st_geometry(WBs_VPU) <- st_make_valid(st_geometry(WBs_VPU))
return(list(POIs = tmp_POIs, events = wb_outlet_events, WBs = WBs_VPU,
nhd_wb_table = nhd))
} else {
return(list(POIs = tmp_POIs, events = NA, WBs = wb_table,
nhd_wb_table = nhd))
}
}
#' Creates Waterbody POIs, calls a few other functions
#' @param nhd (sf data.frame) flowline data.frame
#' @param WBs_VPU (sf data.frame) waterbodies to create inlet pois from
#' @param data_paths (list) list of data paths
#' @param crs (integer) CRS to use (epsg code)
#' @param split_layer (sf data.frame) events to split flowlines with
#' @param tmp_POIs (sf data.frame) rolling POI data frame
#'
#' @return (sf data.frame) dataframe of WB inlet POIs
#'
wbin_POIcreation <- function(nhd, wb_lyr, data_paths, crs,
split_layer, tmp_POIs){
wbout_COMIDs <- filter(tmp_POIs, !is.na(Type_WBOut))
WBs_in_POIs <- filter(wb_lyr, wb_id %in% wbout_COMIDs$Type_WBOut)
# Unnest to get list of waterbody COMIDS (member_comid)
wb_table <- st_drop_geometry(wb_lyr) %>%
dplyr::mutate(member_comid = strsplit(member_comid, ",")) %>%
tidyr::unnest(cols = member_comid) %>%
mutate(member_comid = as.integer(member_comid)) %>%
distinct()
# # Attribute flowline if waterbody exists for POIs
# nhd_wb <- nhd %>%
# left_join(distinct(wb_table, wb_id, member_comid),
# by = c("WBAREACOMI" = "member_comid"))
# Create waterbody inlet POIs
wbin_COMIDs <- filter(nhd, is.na(wb_id) &
DnHydroseq %in% filter(nhd, !is.na(wb_id))$Hydroseq) %>%
select(-wb_id) %>%
switchDiv(., nhd) %>%
inner_join(st_drop_geometry(filter(nhd, minNet == 1)) %>%
select(wb_id, Hydroseq), by = c("DnHydroseq" = "Hydroseq")) %>%
select(COMID, wb_id, minNet) %>%
group_by(COMID) %>%
# Ensure the inlets are on the network defined by confluence POIs
filter(minNet == 1)
# Headwater Waterbodies that may need the network extended through the inlet
need_wbin <- st_drop_geometry(wb_lyr) %>%
#dplyr::select(COMID)%>%
dplyr::filter(!wb_id %in% wbin_COMIDs$wb_id)
if(nrow(need_wbin)>0){
nhd_inlet <- mutate(nhd, WB = ifelse(wb_id > 0 & wb_id %in% need_wbin$wb_id, 1, 0))
missing_wbin_COMIDs <- filter(nhd_inlet, WB == 0,
DnHydroseq %in% filter(nhd_inlet, WB == 1)$Hydroseq) %>%
select(-wb_id) %>%
switchDiv(., nhd_inlet) %>%
inner_join(st_drop_geometry(filter(nhd_inlet, minNet == 1)) %>%
select(Hydroseq, wb_id),
by = c("DnHydroseq" = "Hydroseq")) %>%
select(COMID, wb_id) %>%
group_by(COMID) %>%
filter(row_number() == 1) %>%
ungroup()
missing_wbin_COMIDs <- missing_wbin_COMIDs %>%
left_join(select(st_drop_geometry(nhd_inlet), COMID, TotDASqKM)
, by = c("COMID")) %>%
group_by(COMID) %>%
slice(which.max(TotDASqKM))%>%
ungroup() %>%
select(-TotDASqKM)
if(nrow(missing_wbin_COMIDs) > 0){
wbin_COMIDs <- data.table::rbindlist(list(wbin_COMIDs,
missing_wbin_COMIDs),
fill = TRUE) %>%
select(COMID, wb_id)
}
}
wbin_POIs <- POI_creation(filter(st_drop_geometry(wbin_COMIDs), !COMID %in% wbout_COMIDs$COMID),
nhd, "WBIn")
# Get NHDArea and HR waterbodies
WBs_VPU_areaHR <- HR_Area_coms(nhd,
filter(wb_lyr, source == "HR ID AVAILABLE"),
data_paths, crs)
if(is.list(WBs_VPU_areaHR)) {
nhd_WBs <- filter(nhd, minNet == 1,
COMID %in% WBs_VPU_areaHR$nhd_wb_table$comid,
AreaSqKM > 0)
} else {
nhd_WBs <- NA
}
if(is.data.frame(nhd_WBs)){
# Get the outlet events again
wb_outlet_events <- read_sf(temp_gpkg, split_layer) %>%
filter(event_type == "outlet")
# Drive the inlet events
nhd_inlet <- filter(nhd_WBs, !COMID %in% wb_outlet_events$COMID)
if("12" %in% unique(nhd$VPUID)){
nhd_inlet <- rbind(nhd_inlet, filter(nhd_WBs, COMID == 1304709))
}
if (nrow(nhd_inlet) > 0){
# Get inlet events and bind with outlets
wb_inlet_events <- WB_event(WBs_VPU_areaHR, nhd_inlet, "inlet") %>%
mutate(nexus = TRUE) %>%
inner_join(select(st_drop_geometry(nhd), COMID, Hydroseq, ToMeas), by = "COMID")
# Determine which events are close enough to the end of an upstream flowline to just use that geometry for the POI
wbin_fl_upstream <- filter(nhd, DnHydroseq %in% filter(nhd, COMID %in% wb_inlet_events$COMID)$Hydroseq) %>%
group_by(DnHydroseq) %>%
filter(n() == 1) %>%
ungroup()
if(nrow(wbin_fl_upstream) > 0){
wb_inlet_POIs <- filter(wb_inlet_events, REACH_meas < 5 | as.integer(REACH_meas) == as.integer(ToMeas),
COMID %in% wbin_fl_upstream$toCOMID) %>%
inner_join(select(st_drop_geometry(wbin_fl_upstream), usCOMID = COMID, toCOMID), by = c("COMID" = "toCOMID")) %>%
mutate(dsCOMID = COMID, COMID = usCOMID)
if(nrow(wb_inlet_POIs) > 0) {
wbin_POIs <- bind_rows(wbin_POIs, POI_creation(select(st_drop_geometry(wb_inlet_POIs), dsCOMID,
Type_WBIn = POI_identifier),
nhd, "WBIn"))
wb_inlet_events <- filter(wb_inlet_events, !COMID %in% wb_inlet_POIs$dsCOMID)
}
}
if(nrow(wb_inlet_events) > 0){
wbin_POIs <- addType(wbin_POIs, tmp_POIs, "WBIn", nexus = TRUE)
wb_inlet_events <- st_compatibalize(wb_inlet_events, wbin_POIs)
wbin_POIs_fin <- data.table::rbindlist(list(wbin_POIs,
select(wb_inlet_events, COMID,
Type_WBIn = POI_identifier, nexus)),
fill = TRUE) %>%
st_as_sf()
}
return(list(POIs = wbin_POIs_fin, events = wb_inlet_events))
} else {
print("no waterbody inlets events")
wbin_POIs <- addType(wbin_POIs, tmp_POIs, "WBIn")
return(list(POIs = wbin_POIs, events = NA))
}
} else {
wbin_POIs <- addType(wbin_POIs, tmp_POIs, "WBIn")
return(list(POIs = wbin_POIs, events = NA))
}
}
#' collapses close upstream gages with wb inlet events.
#' @param temp_POIs (sf data.frame) rolling POI data frame
#' @param nhd (sf data.frame) nhd flowlines
#' @param events (list) waterbody inlet events
#'
#' @return (sf data.frame) dataframe of gage POIs
#'
wb_inlet_collapse <- function(tmp_POIs, nhd, events){
gage_dist_node <- function(x, wb_ds_ds, gage_add, events){
print (x)
wb_in_fl <- filter(wb_ds_ds, COMID == x)
gage_ds <- filter(wb_ds_ds, Hydroseq %in% wb_in_fl$Hydroseq |
Hydroseq %in% wb_in_fl$DnHydroseq) %>%
filter(TotDASqKM == max(TotDASqKM))
gage_reach <- gage_ds %>%
group_by(REACHCODE) %>%
summarize(do_union = TRUE,
total_length = sum(LENGTHKM))
#print(nrow(gage_reach))
gage_event <- filter(gage_add, COMID %in% gage_ds$COMID)
wb_event <- filter(events, COMID == x, event_type == "inlet") %>%
unique()
if(nrow(gage_reach) == 0){
print("no gage reaches")
}
if(nrow(gage_event) == 0){
return("No gage events")
} else if(gage_event$COMID != wb_in_fl$COMID) {
gage_reach <- gage_reach %>%
filter(REACHCODE == gage_event$reachcode) %>%
mutate(gage_dist = ifelse(gage_event$nexus == TRUE,
total_length * (1 - (gage_event$reach_meas/100)),
total_length)) %>%
mutate(gage_comid = gage_event$COMID,
wbin_comid = x)
} else if(gage_event$COMID == wb_in_fl$COMID){
if(nrow(wb_event) >0){
wb_in_meas <- wb_event$REACH_meas
wb_RC <- wb_event$REACHCODE
} else {
wb_out_meas <- wb_in_fl$FromMeas
wb_RC <- wb_in_fl$REACHCODE
}
# gage reach
gage_reach <- gage_reach %>%
filter(REACHCODE == gage_event$reachcode) %>%
mutate(gage_dist = total_length * (1 - (gage_event$reach_meas/100)))
# wb info
wb_reach <- gage_reach %>%
filter(REACHCODE == wb_RC) %>%
mutate(wb_dist = total_length * (1 - (wb_out_meas /100)))
gage_reach <- gage_reach %>%
mutate(gage_dist = abs(wb_reach$wb_dist - gage_dist),
gage_comid = gage_event$COMID,
wbin_comid = x)
}
}
#events <- read_sf(temp_gpkg, split_layer) %>%
# rbind(st_compatibalize(wb_,.))
# Previously identified streamgages within Gage_Selection.Rmd
streamgages_VPU <- gages %>%
rename(COMID = comid) %>%
filter(COMID %in% nhd$COMID) %>%
#st_drop_geometry() %>%
switchDiv(., nhd)
# get waterbody outlets
wb_in <- filter(tmp_POIs, !is.na(Type_WBIn), is.na(Type_Gages))
wb_in_nexus <- filter(wb_in, nexus == TRUE)
wb_in_nhd <- filter(nhd, COMID %in% wb_in$COMID)
wb_ds_ds <- wb_in_nhd %>%
rbind(filter(nhd, LevelPathI %in% .$LevelPathI & DnHydroseq %in% .$Hydroseq)) %>%
distinct()
# get gages
gage_wb <- filter(tmp_POIs, !is.na(Type_Gages) & is.na(Type_WBIn)) #%>%
# filter(COMID %in% wb_ds_ds$COMID)
#gage_nhd <- filter(nhd, COMID %in% gage_wb$COMID)
#gage_node <- filter(gage_wb, nexus == FALSE)
#gage_nexus <- filter(gage_wb, nexus == TRUE)
gage_add <- filter(streamgages_VPU, site_no %in% gage_wb$Type_Gages) %>%
select(COMID, reachcode, reach_meas, site_no) %>%
inner_join(select(st_drop_geometry(gage_wb), site_no = Type_Gages, nexus),
by = "site_no") %>%
filter(!nexus)
output <- lapply(wb_in$COMID,
gage_dist_node, wb_ds_ds, gage_add, events)
output_full <- do.call("rbind", output[lengths(output) > 1])
if(!is.null(output_full)){
output_full <- output_full %>%
filter(gage_dist < 1) %>%
left_join(select(st_drop_geometry(nhd), COMID, TotDASqKM_WB = TotDASqKM),
by = c("wbin_comid" = "COMID")) %>%
left_join(select(st_drop_geometry(nhd), COMID, TotDASqKM_gage = TotDASqKM),
by = c("gage_comid" = "COMID")) %>%
mutate(DAR = TotDASqKM_gage / TotDASqKM_WB) %>%
filter(gage_dist < 1, DAR > .975) %>%
st_drop_geometry()
gage_POI <- filter(tmp_POIs, COMID %in% output_full$gage_comid) %>%
select(COMID, Type_HUC12_ds = Type_HUC12, Type_Gages_ds = Type_Gages,
Type_TE_ds = Type_TE, Type_Term_ds = Type_Term, nexus) %>%
st_drop_geometry() %>%
group_by(COMID) %>%
summarise(Type_HUC12_ds = last(na.omit(Type_HUC12_ds)),
Type_Gages_ds = last(na.omit(Type_Gages_ds)),
Type_Term_ds = last(na.omit(Type_Term_ds)),
Type_TE_ds = last(na.omit(Type_TE_ds)),
nexus = last(na.omit(nexus)))
WB_POI <- filter(tmp_POIs, COMID %in% output_full$wbin_comid, !is.na(Type_WBIn)) %>%
left_join(select(output_full, gage_comid, wbin_comid), by = c("COMID" = "wbin_comid")) %>%
left_join(select(gage_POI, -nexus), by = c("gage_comid" = "COMID")) %>%
mutate(Type_HUC12 = ifelse(!is.na(Type_HUC12_ds), Type_HUC12_ds, Type_HUC12),
Type_Gages = ifelse(!is.na(Type_Gages_ds), Type_Gages_ds, Type_Gages),
Type_TE = ifelse(!is.na(Type_TE_ds), Type_TE_ds, Type_TE),
Type_Term = ifelse(!is.na(Type_Term_ds), Type_Term_ds, Type_Term)) %>%
select(-c(Type_HUC12_ds, Type_Gages_ds, Type_TE_ds, Type_Term_ds))
tmp_POIs_fin <- filter(tmp_POIs, !COMID %in% c(WB_POI$COMID, WB_POI$gage_comid)) %>%
rbind(select(WB_POI, -gage_comid))
if(any(gage_POI$nexus == TRUE)){
gage_events <- filter(gage_POI, nexus == TRUE)
events <- filter(events, !POI_identifier %in% gage_events$Type_Gages_ds)
}
return(list(POIs = tmp_POIs_fin, events_ret = events))
} else {
print ("no points collapse")
return(list(POIs = tmp_POIs, events_ret = NA))
}
}
#' Collapses POIs us/ds of waterbody POIs
#' @param (sf data.frame) rolling POI data frame
#' @param nhd (sf data.frame) nhd flowlines
#' @param events (sf data.frame) waterbody events
#'
#' @return (sf data.frame) dataframe of wb inlet POIs collapsed
wb_poi_collapse <- function(tmp_POIs, nhd, events){
gage_dist_node <- function(x, wb_ds_ds, gage_add, events){
print (x)
wb_out_fl <- distinct(filter(wb_ds_ds, COMID == x))
gage_ds <- filter(wb_ds_ds, Hydroseq %in% wb_out_fl$Hydroseq |
Hydroseq %in% wb_out_fl$DnHydroseq)
gage_reach <- gage_ds %>%
group_by(REACHCODE) %>%
summarize(do_union = TRUE,
total_length = sum(LENGTHKM))
#print(nrow(gage_reach))
gage_event <- filter(gage_add, COMID %in% gage_ds$COMID) %>%
filter(reach_meas == min(reach_meas))
wb_event <- filter(events, COMID == x, event_type == "outlet") %>%
unique()
if(nrow(gage_reach) == 0){
print("no gages")
}
if(nrow(gage_event) == 0){
return("no events")
} else if(gage_event$COMID != wb_out_fl$COMID) {
gage_reach <- gage_reach %>%
filter(REACHCODE == unique(gage_event$reachcode)) %>%
mutate(gage_dist = ifelse(gage_event$nexus == TRUE,
total_length * (1 - (gage_event$reach_meas/100)),
total_length)) %>%
mutate(gage_comid = gage_event$COMID,
wbout_comid = x)
} else if(gage_event$COMID == wb_out_fl$COMID){
if(nrow(wb_event) >0){
wb_out_meas <- min(wb_event$REACH_meas)
wb_RC <- wb_event$REACHCODE
} else {
wb_out_meas <- min(wb_out_fl$FromMeas)
wb_RC <- wb_out_fl$REACHCODE
}
# gage reach
gage_reach <- gage_reach %>%
filter(REACHCODE == gage_event$reachcode) %>%
mutate(gage_dist = total_length * (1 - (gage_event$reach_meas/100)))
# wb info
wb_reach <- gage_reach %>%
filter(REACHCODE == unique(wb_RC)) %>%
mutate(wb_dist = total_length * (1 - (wb_out_meas /100)))
gage_reach <- gage_reach %>%
mutate(gage_dist = abs(wb_reach$wb_dist - gage_dist),
gage_comid = gage_event$COMID,
wbout_comid = x)
}
}
# Previously identified streamgages within Gage_Selection.Rmd
streamgages_VPU <- gages %>%
rename(COMID = comid) %>%
filter(COMID %in% nhd$COMID) %>%
switchDiv(., nhd)
# get waterbody outlets
wb_out <- filter(tmp_POIs, !is.na(Type_WBOut), is.na(Type_Gages))
wb_out_nexus <- filter(wb_out, nexus)
wb_out_nhd <- filter(nhd, COMID %in% wb_out$COMID)
wb_ds_ds <- wb_out_nhd %>%
rbind(filter(nhd, LevelPathI %in% .$LevelPathI & Hydroseq %in% .$DnHydroseq))
# get gages
gage_wb <- filter(tmp_POIs, !is.na(Type_Gages) & is.na(Type_WBOut)) #%>%
# filter(COMID %in% wb_ds_ds$COMID)
#gage_nhd <- filter(nhd, COMID %in% gage_wb$COMID)
#gage_node <- filter(gage_wb, nexus == FALSE)
#gage_nexus <- filter(gage_wb, nexus == TRUE)
gage_add <- filter(streamgages_VPU, site_no %in% gage_wb$Type_Gages) %>%
select(COMID, reachcode, reach_meas, site_no) %>%
inner_join(select(st_drop_geometry(gage_wb), site_no = Type_Gages, nexus),
by = "site_no") %>%
distinct()
# 8693865
output <- lapply(wb_out$COMID,
gage_dist_node, wb_ds_ds, gage_add, events)
output_length <- output[lengths(output) > 1]
if(length(output_length) == 0){
return(list(POIs = tmp_POIs, events_ret = NA))
}
output_full <- do.call("rbind", output[lengths(output) > 1]) %>%
filter(gage_dist < 1) %>%
left_join(select(st_drop_geometry(nhd), COMID, TotDASqKM_WB = TotDASqKM),
by = c("wbout_comid" = "COMID")) %>%
left_join(select(st_drop_geometry(nhd), COMID, TotDASqKM_gage = TotDASqKM),
by = c("gage_comid" = "COMID")) %>%
mutate(DAR = TotDASqKM_WB / TotDASqKM_gage) %>%
filter(gage_dist < 1, DAR > .975) %>%
st_drop_geometry()
gage_POI <- filter(tmp_POIs, COMID %in% output_full$gage_comid) %>%
select(COMID, Type_HUC12_ds = Type_HUC12, Type_Gages_ds = Type_Gages,
Type_TE_ds = Type_TE, nexus) %>%
st_drop_geometry() %>%
group_by(COMID) %>%
summarise(Type_HUC12_ds = last(na.omit(Type_HUC12_ds)),
Type_Gages_ds = last(na.omit(Type_Gages_ds)),
Type_TE_ds = last(na.omit(Type_TE_ds)),
nexus = last(na.omit(nexus)))
WB_POI <- filter(tmp_POIs, COMID %in% output_full$wbout_comid, !is.na(Type_WBOut)) %>%
inner_join(select(output_full, gage_comid, wbout_comid), by = c("COMID" = "wbout_comid")) %>%
inner_join(select(gage_POI, -nexus), by = c("gage_comid" = "COMID")) %>%
mutate(Type_HUC12 = ifelse(!is.na(Type_HUC12_ds), Type_HUC12_ds, Type_HUC12),
Type_Gages = ifelse(!is.na(Type_Gages_ds), Type_Gages_ds, Type_Gages),
Type_TE = ifelse(!is.na(Type_TE_ds), Type_TE_ds, Type_TE)) %>%
select(-c(Type_HUC12_ds, Type_Gages_ds, Type_TE_ds))
tmp_POIs_fin <- filter(tmp_POIs, !COMID %in% c(WB_POI$COMID, WB_POI$gage_comid)) %>%
rbind(select(WB_POI, -gage_comid))
if(any(gage_POI$nexus == TRUE)){
gage_events <- filter(gage_POI, nexus == TRUE)
events <- filter(events, !POI_identifier %in% gage_events$Type_Gages_ds)
}
return(list(POIs = tmp_POIs_fin, events_ret = events))
}
#' Collapses POIs together based on criteria
#' @param (pois) sf data frame of POIs
#' @param move_category (character) POI data theme to move
#' @param DAR (numeric) drainage area threshold to move within
#' @param dist (numeric) maximum river distance between two points to move within
#' @param keep_category (character) POI data themes to keep static
#'
#' @return (sf data.frame, table) dataframe of pois, table of points that have moved
poi_move <- function(pois, move_category, DAR, dist, keep_category) {
# filter out features with identical geometry
# Add row_number
pois_edit <- pois %>%
mutate(nexus = ifelse(is.na(nexus), 0, nexus))
# Don't consider points already moved
if("moved" %in% colnames(pois_edit)){
pois_tomove <- filter(pois_edit, is.na(moved)) # change from poi_edit
pois_moved_pre <- filter(pois_edit, !is.na(moved))}
# If 'keep' category included
if(!missing(keep_category)){
poi2move <- filter(pois_tomove, !is.na(.data[[move_category]]) & nexus == FALSE) %>%
filter(if_all(!!as.symbol(keep_category), function(x) is.na(x))) %>%
# Never move these
filter_at(vars(Type_WBOut, Type_WBIn, Type_Conf, Type_Term), all_vars(is.na(.)))
pois2keep <- filter(pois_tomove, !id %in% poi2move$id)
#is.na(.data[[move_category]]) & nexus == FALSE) #%>%
#filter(if_all(!!as.symbol(keep_category), function(x) is.na(x)))
} else {
# POIs to move
poi2move <- pois_tomove %>%
filter_at(vars(Type_WBOut, Type_WBIn, Type_Conf, Type_Term), all_vars(is.na(.))) %>%
filter(nexus == 0) %>%
filter(!is.na(.data[[move_category]]))
pois2keep <- filter(pois_tomove, !id %in% poi2move$id)
}
# Get relevant NHD data
nhd_poi1 <- filter(nhd, COMID %in% pois2keep$COMID)
nhd_poi2 <- filter(nhd, COMID %in% poi2move$COMID)
# Ensure they are on same level path
nhd_poi2 <- filter(nhd_poi2, LevelPathI %in% nhd_poi1$LevelPathI)
# Join NHD data
pois2keep_nhd <- pois2keep %>%
inner_join(select(st_drop_geometry(nhd_poi1), COMID, LevelPathI, Hydroseq,
DA_keep = TotDASqKM, Pathlength_keep = Pathlength), by = "COMID") %>%
rename(COMID_keep = COMID)
# Join NHD data
pois2move_nhd <- select(poi2move, COMID, !!as.symbol(move_category), id_move = id) %>%
inner_join(select(st_drop_geometry(nhd_poi2), COMID, LevelPathI, Hydroseq, TotDASqKM, Pathlength),
by = "COMID")
# Candidates to move
pois2move_cand <-pois2move_nhd %>%
inner_join(select(st_drop_geometry(pois2keep_nhd), COMID_keep, DA_keep, LevelPathI,
Pathlength_keep, id_keep = id, nexus),
by = "LevelPathI") %>%
mutate(river_dist = abs(Pathlength - Pathlength_keep), DAR_poi = abs(DA_keep/TotDASqKM),
move_dir = ifelse(Pathlength < Pathlength_keep, "Up", "Down")) %>%
group_by(id_move, move_dir) %>%
ungroup() %>%
filter((river_dist < dist) & (DAR_poi > (1 - DAR)) & (DAR_poi < (1 + DAR))) %>%
select(!!as.symbol(move_category), id_move, COMID, id_keep, COMID_keep, river_dist, DAR_poi, move_dir, nexus) %>%
st_drop_geometry()
move_distinct <- pois2move_cand %>%
group_by(id_keep) %>%
filter(row_number() == 1) %>%
ungroup() %>%
distinct(id_move, COMID_move = COMID, id_keep, COMID_keep, river_dist, DAR_poi, move_dir, nexus) %>%
group_by(id_move) %>%
slice(which.min(abs(1 - DAR_poi)))
if(nrow(move_distinct) == 0){
print("no POIs to move")
return(pois)
}
pois2_move <- filter(st_drop_geometry(pois_tomove), id %in% move_distinct$id_move) %>%
select_if(~sum(!is.na(.)) > 0) %>%
select(-c(COMID, nexus)) %>%
inner_join(select(move_distinct, id_move, id_keep), by = c("id" = "id_move"))
move_fields <- colnames(select(pois2_move, -c(id, id_keep)))
if(length(move_fields) == 1){
pois2_keep <- filter(pois_tomove, id %in% pois2_move$id_keep, !id %in% pois2_move$id) %>%
inner_join(select(pois2_move, id_move = id, id_keep,
new_val = !!as.symbol(move_category)), by = c("id" = "id_keep")) %>%
mutate(moved := ifelse(is.na(!!as.symbol(move_category)),
id_move, moved),
!!as.symbol(move_category) := ifelse(is.na(!!as.symbol(move_category)),
new_val, !!as.symbol(move_category)))
moved_points <- filter(pois2_keep, !is.na(new_val), !is.na(moved)) %>%
mutate(moved_value = move_category)
} else {
for (field in move_fields){
pois2_keep <- filter(pois_tomove, id %in% pois2_move$id_keep, !id %in% pois2_move$id) %>%
inner_join(select(pois2_move, id_move = id, id_keep, new_val = !!as.symbol(field)),
by = c("id" = "id_keep")) %>%
mutate(moved := ifelse(is.na(!!as.symbol(field)),
id_move, moved),
!!as.symbol(field) := ifelse(is.na(!!as.symbol(field)),
new_val, !!as.symbol(field)))
pois_moved <- filter(pois2_keep, !is.na(new_val), !is.na(moved)) %>%
mutate(moved_value = field)
if(!exists("moved_points")){
moved_points <- pois_moved
} else {
moved_points <- rbind(moved_points, pois_moved)
}
}
}
pois_final <- data.table::rbindlist(list(filter(pois_edit, !id %in% c(moved_points$id_move, pois2_keep$id)),
select(pois2_keep, -c(new_val, id_move, new_val))), fill = TRUE) %>%
st_as_sf()
return(list(final_pois = pois_final, moved_points = moved_points))
}
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)
}
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)
}