From 63f8985dfac56659f3b807db9440b7f2520a7280 Mon Sep 17 00:00:00 2001
From: Bock <abock@usgs.gov>
Date: Mon, 8 May 2023 10:04:15 -0600
Subject: [PATCH] Some minor mods, and better POI collapse
---
workspace/R/NHD_navigate.R | 205 +++++++++++++++++++++++++++++++------
1 file changed, 171 insertions(+), 34 deletions(-)
diff --git a/workspace/R/NHD_navigate.R b/workspace/R/NHD_navigate.R
index dfeaba0..f2d45e0 100644
--- a/workspace/R/NHD_navigate.R
+++ b/workspace/R/NHD_navigate.R
@@ -98,7 +98,7 @@ segment_creation <- function(nhdDF, routing_fix){
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 (missing(routing_fix) == FALSE){
+ if (!is.na(routing_fix)){
routing_fix <- routing_fix %>%
rename(COMID = oldPOI, new_COMID = COMID)
@@ -207,19 +207,22 @@ DS_poiFix <- function(POIs_wgeom, nhd){
#'
#' @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 <- nhdplusTools::get_UM(nhdDF, poi_fix, sort=T)
- seg2fix <- filter(nhdDF, COMID == poi_fix)
+ 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(factor(COMID, levels = up_segs)) %>%
+ arrange(Hydroseq) %>%
filter(AreaSqKM > 0)
- down_segs <- nhdplusTools::get_DM(nhdDF, poi_fix, sort=T)
+ down_segs <- unique(nhdplusTools::get_DM(nhdDF, poi_fix, sort=T))
downstuff <- filter(nhdDF, COMID %in% unlist(down_segs)) %>%
- arrange(factor(COMID, levels = down_segs)) %>%
+ arrange(Hydroseq)%>%
filter(AreaSqKM > 0)
# combine into one dataframe, select up/downstream seg with least change in total drainage area
@@ -966,11 +969,16 @@ gage_POI_creation <- function(tmp_POIs, gages_info, nhd, combine_meters, reach_m
gage_POIs_nonevent <- filter(gage_POIs, !Type_Gages %in% events$Type_Gages) %>%
addType(., tmp_POIs, "Gages", nexus = FALSE, bind = TRUE)
- 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()
-
+ 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))
}
@@ -985,7 +993,7 @@ gage_POI_creation <- function(tmp_POIs, gages_info, nhd, combine_meters, reach_m
#'
wbout_POI_creaton <- function(nhd, WBs_VPU, data_paths, crs){
# Create waterbody outlet POIs for waterbodies that are in NHDv2 waterbody set
- wbout_COMIDs <- filter(nhd, dend == 1 & WB == 1) %>%
+ wbout_COMIDs <- nhd %>% #filter(nhd, dend == 1 & WB == 1)
group_by(WBAREACOMI) %>%
slice(which.min(Hydroseq)) %>%
switchDiv(., nhd) %>%
@@ -1074,7 +1082,7 @@ wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs,
filter(minNet == 1)
# Headwater Waterbodies that may need the network extended through the inlet
- need_wbin <- st_drop_geometry(filter(WBs_VPU, source == "NHDv2WB")) %>%
+ need_wbin <- st_drop_geometry(filter(WBs_VPU, source %in% c("ref_WB", "NHDv2WB"))) %>%
dplyr::select(COMID)%>%
dplyr::filter(!COMID %in% wbin_COMIDs$WBAREACOMI)
@@ -1151,7 +1159,7 @@ wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs,
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 = WBAREACOMI),
+ wbin_POIs <- bind_rows(wbin_POIs, POI_creation2(select(st_drop_geometry(wb_inlet_POIs), dsCOMID, Type_WBIn = WBAREACOMI),
nhd, "WBIn"))
wb_inlet_events <- filter(wb_inlet_events, !COMID %in% wb_inlet_POIs$dsCOMID)
@@ -1208,7 +1216,7 @@ wb_inlet_collapse <- function(tmp_POIs, nhd, events){
unique()
if(nrow(gage_reach) == 0){
- print("Wilton")
+ print("no gage reaches")
}
if(nrow(gage_event) == 0){
@@ -1333,11 +1341,9 @@ wb_inlet_collapse <- function(tmp_POIs, nhd, events){
#' @param events (sf data.frame) waterbody inlet events
#'
#' @return (sf data.frame) dataframe of wb inlet POIs collapse
-#'
-#' wb_poi_collapse <- function(tmp_POIs, nhd, events){
wb_poi_collapse <- function(tmp_POIs, nhd, events){
gage_dist_node <- function(x, wb_ds_ds, gage_add, events){
- print (x)
+ print (x) #6116850
wb_out_fl <- filter(wb_ds_ds, COMID == x)
gage_ds <- filter(wb_ds_ds, Hydroseq %in% wb_out_fl$Hydroseq |
Hydroseq %in% wb_out_fl$DnHydroseq)
@@ -1353,15 +1359,14 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
unique()
if(nrow(gage_reach) == 0){
- print("Wilton")
+ print("no gages")
}
if(nrow(gage_event) == 0){
- #print("akermayun")
- return("Akermayun")
+ return("no events")
} else if(gage_event$COMID != wb_out_fl$COMID) {
gage_reach <- gage_reach %>%
- filter(REACHCODE == gage_event$reachcode) %>%
+ filter(REACHCODE == unique(gage_event$reachcode)) %>%
mutate(gage_dist = ifelse(gage_event$nexus == TRUE,
total_length * (1 - (gage_event$reach_meas/100)),
total_length)) %>%
@@ -1369,10 +1374,10 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
wbout_comid = x)
} else if(gage_event$COMID == wb_out_fl$COMID){
if(nrow(wb_event) >0){
- wb_out_meas <- wb_event$REACH_meas
+ wb_out_meas <- min(wb_event$REACH_meas)
wb_RC <- wb_event$REACHCODE
} else {
- wb_out_meas <- wb_out_fl$FromMeas
+ wb_out_meas <- min(wb_out_fl$FromMeas)
wb_RC <- wb_out_fl$REACHCODE
}
@@ -1383,7 +1388,7 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
# wb info
wb_reach <- gage_reach %>%
- filter(REACHCODE == wb_RC) %>%
+ filter(REACHCODE == unique(wb_RC)) %>%
mutate(wb_dist = total_length * (1 - (wb_out_meas /100)))
gage_reach <- gage_reach %>%
@@ -1393,14 +1398,11 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
}
}
- #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
@@ -1421,10 +1423,16 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
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")
+ by = "site_no") %>%
+ distinct()
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) %>%
@@ -1438,12 +1446,11 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
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) %>%
+ 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_Term_ds = last(na.omit(Type_Term_ds)),
Type_TE_ds = last(na.omit(Type_TE_ds)),
nexus = last(na.omit(nexus)))
@@ -1452,9 +1459,8 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
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),
- 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))
+ 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))
@@ -1466,4 +1472,135 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
return(list(POIs = tmp_POIs_fin, events_ret = events))
+}
+
+#' Creates Waterbody POIs, calls a few other functions
+#' @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))
+
}
\ No newline at end of file
--
GitLab