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

Changes to functions associated with NHD_Navigate.R

parent 9845be59
No related branches found
No related tags found
1 merge request!43NHD_Navigate cleanup/simplicication
Hide whitespace changes
Inline Side-by-side
workspace/R/NHD_navigate.R
+ 334
− 275
View file @ b893d0db
# NHDPlus V2 Region
if(!exists("hydReg") && !nchar(hydReg <- Sys.getenv("HYDREG")) > 1)
hydReg <- "01"
if(!exists("VPU") && !nchar(VPU <- Sys.getenv("HYDREG")) > 1)
VPU <- "01"
# Source Data
crs <- 5070
data_paths <- data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
huc12_xWalk <- file.path(data_paths$nhdplus_dir, "CrosswalkTable_NHDplus_HU12.csv")
HUC12 <- file.path(data_paths$nhdplus_dir, "HUC12.rds")
# Defined during NHD_Navigate.Rmd
huc12_pois <- paste0("HUC12_", hydReg) # HUC12 POIs
dup_comids <- paste0("cache/dupCOMIDs_",hydReg, ".rds") # data frame of COMIDs with multiple HUC12 assignments
nhdflowline <- "NHDFlowline" # NHD flowlines pers VPU
WBs <- paste0("WB_", hydReg) # Waterbodies within VPU
WBout_POIs <- paste0("WBOut_", hydReg) # Waterbody Outlet POIs
WBin_POIs <- paste0("WBIn_", hydReg) # Waterbody Inlet POIs
gages_pois <- paste0("Gages_", hydReg) # GAGESIII POIs
gageLoc <- paste0("gageLoc_", hydReg) # GageLoc Files
TE_pois <- paste0("TE_", hydReg) # Thermoelectric POIs
NID_pois <- paste0("NID_", hydReg) # NID POIs
hucgage_segs <- paste0("hucgagesegs_", hydReg) # inimally-sufficient network
conf_pois <- paste0("Conf_", hydReg) # Confluence POIs
pois_all <- paste0("POIs_", hydReg) # All POIs binded together
nsegment_raw <- paste0("nsegment_raw_", hydReg) # Minimally-sufficient network attributed with POI_ID
n_segments <- paste0("Nsegment_", hydReg) # Minimally-sufficient network dissolved by POI_ID
struct_POIs <- paste0("structPOIs_", hydReg) # structual POIs
struc_Net <- paste0("structNet", hydReg) # structural network
coastFlowlines <- "coast_NCA" # Coastal dendritic segments (modify this)
# Defined during POI_Collapse.Rmd
nav_gpkg <- file.path("cache", paste0("GF_",hydReg,".gpkg"))
col_gpkg <- file.path("cache", paste0("GF_",hydReg,"_Collapse.gpkg"))
nav_gpkg <- file.path("cache", paste0("GF_", VPU,".gpkg"))
nhdflowline <- "NHDFlowline" # NHD flowlines per VPU
temp_POIs <- paste0("tmp_POIs_", VPU) # Running POI list, precedes "POIs_VPU"
WBs <- paste0("WB_", VPU) # Waterbodies within VPU
pois_all <- paste0("POIs_", VPU) # All POIs binded together
poi_moved <- paste0("POIs_mv_", VPU) # POIs moved from original COMID assignment
poi_xWalk <- paste0("poi_xWalk_", VPU) # POIs that changed COMIDS during the navigate part of the workflow
n_segments <- paste0("nsegment_", VPU) # Minimally-sufficient network dissolved by POI_ID
pois_merge <- paste0("merPOIs_", VPU) # All POIs binded together
# Making an "events" script to address "events" information
# i.e. gages, NID, TE Plants, waterbodies
gage_gpkg <- "cache/Gage_info.gpkg"
gage_cand <- "Potential_Gages"
crs <- 5070
# Output Layers
pois_merge <- paste0("merPOIs_", hydReg) # All POIs binded together
nsegment_raw_merge <- paste0("nsegment_rawMG_", hydReg)
nsegment_merge <- paste0("nsegment_MG_", hydReg)
hydReg_gage_Table <- file.path("cache", paste0("R", hydReg, "_Gages.csv"))
VPU_gage_Table <- file.path("cache", paste0("R", VPU, "_Gages.csv"))
CONUS_gage_Table <- "data/gages_MDA.rds"
# Defined during NonDend.Rmd
nd_gpkg <- file.path("cache", paste0("ND_",hydReg,".gpkg"))
coast_FL <- paste0("coastal_FL_", hydReg)
nhd_aggLine <- paste0("nhd_aggline_", hydReg)
agg_fline <- paste0("agg_fline_", hydReg)
reg_cats_lyr <- paste0("nhd_cats", hydReg)
agg_cats_lyr <- paste0("agg_cats", hydReg)
outlets_lyr <- paste0("outlets", hydReg)
lookup_table <- file.path("cache/", paste0(hydReg, "_lookup.csv"))
cat0_lyr <- paste0("cat0_", hydReg)
hru0_lyr <- paste0("hru0_", hydReg)
nd_table_out <- file.path("cache/", paste0(hydReg, "_ND_lookup.csv"))
nhdCoast_net <- paste0("nhdCoast_network_", hydReg)
catCoast_net <- paste0("catCoast_", hydReg)
nd_gpkg <- file.path("cache", paste0("ND_",VPU,".gpkg"))
coast_FL <- paste0("coastal_FL_", VPU)
nhd_aggLine <- paste0("nhd_aggline_", VPU)
agg_fline <- paste0("agg_fline_", VPU)
reg_cats_lyr <- paste0("nhd_cats", VPU)
agg_cats_lyr <- paste0("agg_cats", VPU)
outlets_lyr <- paste0("outlets", VPU)
lookup_table <- file.path("cache/", paste0(VPU, "_lookup.csv"))
cat0_lyr <- paste0("cat0_", VPU)
hru0_lyr <- paste0("hru0_", VPU)
nd_table_out <- file.path("cache/", paste0(VPU, "_ND_lookup.csv"))
nhdCoast_net <- paste0("nhdCoast_network_", VPU)
catCoast_net <- paste0("catCoast_", VPU)
cat0_lyr <- "catchment_coast" # modify this
hru0_lyr <- "hru_coast" # modify this
sink_lyr_out <- paste0("sink_", hydReg)
missing_cat_lyr <- paste0("missingCats_", hydReg)
sink_cats_out <- paste0("Terminal_Cats_", hydReg)
HUC12_lyr_out <- paste0("HUC12_", hydReg)
cat_final_out <- paste0("Final_Cats_", hydReg)
sink_net_out <- paste0("sink_Net_", hydReg)
rem_cat_out <- paste0("Rem_Cats_", hydReg)
hru_NCA_out <- paste0("HRU_NCA_", hydReg)
sink_lyr_out <- paste0("sink_", VPU)
missing_cat_lyr <- paste0("missingCats_", VPU)
sink_cats_out <- paste0("Terminal_Cats_", VPU)
HUC12_lyr_out <- paste0("HUC12_", VPU)
cat_final_out <- paste0("Final_Cats_", VPU)
sink_net_out <- paste0("sink_Net_", VPU)
rem_cat_out <- paste0("Rem_Cats_", VPU)
hru_NCA_out <- paste0("HRU_NCA_", VPU)
NetworkNav <- function(inCom, nhdDF, withTrib){
......@@ -101,8 +87,7 @@ NetworkConnection <- function(incom, nhd){
# mydata : (list) list of COMIDs to connect dangle to network
##########################################
# data frame of connections that need to be made
upnet_DF <- nhd %>%
dplyr::filter(COMID %in% incom) %>%
upnet_DF <- filter(nhd, COMID %in% incom) %>%
filter(!DnHydroseq %in% Hydroseq)
# while the number of dangles is greater than 0
......@@ -111,7 +96,7 @@ NetworkConnection <- function(incom, nhd){
count <- dim(upnet_DF)[1]
print (dim(upnet_DF))
# find out which level paths are downstream of dangling huc12 POIs
DSLP <- upnet_DF$DnLevelPat[upnet_DF$COMID %in% incom]
DSLP <- upnet_DF %>% pull(DnLevelPat)#[upnet_DF$COMID %in% incom]
# Get the COMID of the hydroseq with level path value
# the lowest downstream flowline within the levelpath
inCom2 <- nhd$COMID[nhd$Hydroseq %in% DSLP]
......@@ -152,11 +137,13 @@ switchDiv <- function(insegs, nhd){
upstream <- st_drop_geometry(nhd) %>%
inner_join(st_drop_geometry(div_segs) %>%
filter(Divergence == 2) %>%
select(COMID, Hydroseq), by = c("DnMinorHyd" = "Hydroseq"))
select(COMID, Hydroseq),
by = c("DnMinorHyd" = "Hydroseq"))
# From Upstream Segment switch POI to the downstream major/main path
insegs_maj <- st_drop_geometry(nhd) %>%
inner_join(upstream %>% select(COMID.y, DnHydroseq), by = c("Hydroseq" = "DnHydroseq")) %>%
inner_join(upstream %>% select(COMID.y, DnHydroseq),
by = c("Hydroseq" = "DnHydroseq")) %>%
select(COMID, COMID.y)
insegs <- insegs %>%
......@@ -202,7 +189,7 @@ POI_creation<-function(srcData, nhd, IDfield){
return(POIs)
}
addType<-function(POIs, newPOIs, IDfield){
addType<-function(new_POIs, POIs, IDfield, bind){
##########################################
# Checks if existing POIs are co-located with new POI set
# Adds the type attribute for co-located POIs of multiple themes
......@@ -216,20 +203,251 @@ addType<-function(POIs, newPOIs, IDfield){
# mydata : (data frame) Existing POIs with Type field modified for
# duplicate POIs for given data theme
##########################################
new_POIs <- st_compatibalize(new_POIs, POIs)
POIs <- POIs %>%
POIs_exist <- POIs %>%
rename(ID = !!(paste0("Type_", IDfield)))
POIs_fin <- POIs %>% left_join(filter(st_drop_geometry(newPOIs) %>%
select(COMID, !!paste0("Type_", IDfield)), COMID %in% POIs$COMID), by = "COMID", all.x = TRUE) %>%
rename(ID2 = !!(paste0("Type_", IDfield))) %>%
mutate(ID = ifelse(!is.na(ID2), ID2, NA)) %>%
POIs_newAtt <- filter(new_POIs, COMID %in% POIs$COMID) %>%
rename(ID2 = !!(paste0("Type_", IDfield)))
POIs_fin <- POIs_exist %>% left_join(st_drop_geometry(POIs_newAtt) %>%
select(COMID, ID2), by = "COMID", all.x = TRUE) %>%
mutate(ID = ifelse(!is.na(ID2), ID2, ID)) %>%
rename(!!(paste0("Type_", IDfield)) := ID) %>%
select(-c(ID2))
select(-ID2)
# Bind unless indicated not to
if(missing(bind)){
POIs_fin <- rbind(POIs_fin, filter(new_POIs, !COMID %in% POIs_newAtt$COMID))
}
return(POIs_fin)
}
segment_increment <- function(nhd, POIs){
##########################################
# Creates raw and dissolved segment layers with necessaary
# upstream/downstream routing attribution
#
# Arguments:
# nhd : (data frame) minimally-sufficient flowlines that connect POIs
# and meet requirements of NHDPlusTools
# POIs : (data frame) Raw POI data frame
# Field : (character) ID field of unique Segment identifier
#
# Returns:
# segList : (data frame) raw segments
#
##########################################
#nhd_orig <- filter(nhd, dend == 1)
#nhd <- nhd_orig
seg_POIs <- arrange(POIs, desc(LevelPathI), desc(Hydroseq))
nhd <- mutate(nhd, POI_ID = 0)
# Populate NHD with POI_ID
# Begin with upstream most Levelpath/Hydrosequenc
for (LP in unique(seg_POIs$LevelPathI)){
#print(LP)
# For each level path
seg_POIs_LP <- filter(seg_POIs, LevelPathI == LP)
# Go through sorted list of US -> DS comids
for (POI in seg_POIs_LP$COMID){
#print(POI)
# Re-iterate only on flowlines that have no POI_ID to increase speed
nhd <- filter(nhd, POI_ID == 0)
# Get upstream segments
Seg <- nhdplusTools::get_UM(nhd, POI, include = T)
# Assign POI_ID
nhd <- nhd %>%
mutate(POI_ID = ifelse(COMID %in% Seg & POI_ID == 0, POI, POI_ID))
# Derive incremental segment composed of flowlines
inc_seg <- filter(nhd, POI_ID == POI)
#print (dim(incSeg))
if(!exists("inc_segs")){
inc_segs <- inc_seg
}else{
inc_segs <- rbind(inc_segs, inc_seg)
}
}
}
print(proc.time()-ptm)
return(inc_segs)
}
segment_creation <- function(nhd, routing_fix){
##########################################
# Creates raw and dissolved segment layers with necessaary
# upstream/downstream routing attribution
#
# Arguments:
# inSegs : (data frame) segments with incremental IDs (POI_ID populated)
# nhd : (data frame) full nhd data frame
# routing_fix : (logical) routing fixes if available
#
# Returns:
# segList : (data frame) raw segments
#
##########################################
in_segs <- filter(nhd, !is.na(POI_ID))
routing_fix <- routing_fix %>%
rename(COMID = oldPOI, new_COMID = COMID)
# 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){
# Above we generated the network using the initial set of POIs; here we crosswalk over the old COMIDs to the new
nhd_fix <- nhd %>%
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), TT = sum(PathTimeMA),
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(nhd), !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, TT, To_POI_ID)
return(list(diss_segs = nsegments_fin, raw_segs = in_segs))
}
DS_poiFix <- function(POIs_wgeom, nhd){
##########################################
# Moves POI Upstream or downstream if it falls on COMID
# of flowline with no corresponding catchment
#
# Arguments:
# POIs : (data frame) POI data set
# nhd : (data frame) valid data frame of NHD flowlines
#
# Returns:
# repPOIs_unique : (data frame) DF of POIs with new COMID associated
##########################################
POIs <- st_drop_geometry(POIs_wgeom) %>%
arrange(COMID)
# DF of downstream segment
tocomDF <- select(nhd, COMID, DnHydroseq) %>%
inner_join(select(st_drop_geometry(nhd), COMID, Hydroseq), by = c("DnHydroseq" = "Hydroseq"))
# Find segments with POIs where there is no corresponding catchment
unCon_POIs <- filter(nhd, COMID %in% POIs$COMID, AreaSqKM == 0)
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"))
# 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_fix_unique <- filter(POIs, COMID %in% poi_fix$COMID) %>%
bind_rows(poi_fix) %>%
group_by(COMID) %>%
filter(n() > 1) %>%
summarise_each(funs(min(unique(.[!is.na(.)]), na.rm = T))) %>%
mutate_all(~na_if(., -Inf))
# Join new POI COMIDs and geometry with the old Type fields
new_POIs <- bind_rows(filter(poi_fix, !COMID %in% poi_fix_unique$COMID), poi_fix_unique) %>%
arrange(COMID) %>%
bind_cols(get_node(filter(nhd, COMID %in% .$COMID) %>% arrange(COMID), position = "end")) %>%
mutate(to_com = NA) %>%
st_sf() %>%
st_compatibalize(., POIs_wgeom)
# If the DS replacement flowlines with Inc Area > 0 equals the number of replacement comIDs
replacement_tocoms <- st_drop_geometry(new_POIs) %>%
inner_join(select(st_drop_geometry(nhd), COMID, Hydroseq), by = c("DnHydroseq" = "Hydroseq")) %>%
mutate(to_com = COMID.y) %>%
select(oldPOI, COMID.x, COMID.y, AreaSqKM)
# Add downstream COMIDS to the newCOMID
fin_POIs <- new_POIs %>%
left_join(replacement_tocoms, by = "oldPOI") %>%
mutate(to_com = ifelse(!is.na(COMID.y), COMID.y, NA)) %>%
select(-c(TotDASqKM, AreaSqKM.x, DnHydroseq, COMID.x, COMID.y, AreaSqKM.y))
return (list(xWalk = select(poi_fix, COMID, oldPOI), new_POIs = fin_POIs))
}
movePOI_NA_DA <- function(poi_fix, nhd){
##########################################
# 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
#
# Arguments:
# poi_fix : (data frame) POI data set of COMIDs to be changed
# nhd : (data frame) valid data frame of NHD flowlines
#
# Returns:
# newPOI : (data frame) DF of POIs with new COMID associated
##########################################
# Closest POI/US/DS
up_segs <- nhdplusTools::get_UM(nhd, poi_fix, sort=T)
seg2fix <- filter(nhd, COMID == poi_fix)
# Sorted results and filter out all flowlines w/o catchments
upstuff <- filter(nhd, COMID %in% unlist(up_segs)) %>%
arrange(factor(COMID, levels = up_segs)) %>%
filter(AreaSqKM > 0)
down_segs <- nhdplusTools::get_DM(nhd, poi_fix, sort=T)
downstuff <- filter(nhd, COMID %in% unlist(down_segs)) %>%
arrange(factor(COMID, levels = down_segs)) %>%
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)
}
#out_gages <- POI_move_down(nav_gpkg, final_POIs, nsegments_fin, filter(nhd_Final, !is.na(POI_ID)), "Type_Gages", .05)
POI_move_down<-function(out_gpkg, POIs, Seg, Seg2, exp, DA_diff){
##########################################
# Moves down/up POIs based on data theme
......@@ -252,40 +470,41 @@ POI_move_down<-function(out_gpkg, POIs, Seg, Seg2, exp, DA_diff){
# Don't want to move down to TE facility location
segs_down <- Seg %>%
inner_join(select(st_drop_geometry(.), c(POI_ID, TotalDA, TotalLength)),
by = c("To_POI_ID" = "POI_ID")) %>%
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(is.na(Type_WBOut), is.na(Type_WBIn), is.na(Type_Conf), is.na(Type_TE)) %>%
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) %>%
st_drop_geometry() %>%
# If the drainage area ratio is within acceptable range (1.05)
filter(DAR < (1 + DA_diff)) %>%
# 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, DAR, IncDA)
# Segs/POIS to collapose upstream if HUC12 or Gage
up_POIs <- POIs %>% filter(is.na(Type_Conf))
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(is.na(Type_WBOut), is.na(Type_WBIn), is.na(Type_Conf), is.na(Type_TE)) %>%
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) %>% st_drop_geometry() %>%
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)) %>%
# 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, DAR, IncDA)
# compiled list of POIs to move up or down
seg_choices <- POIs %>%
left_join(segs_down %>%
left_join(st_drop_geometry(segs_down) %>%
select(POI_ID, To_POI_ID, DAR, IncDA), by = c("COMID" = "POI_ID")) %>%
left_join(segs_up %>%
left_join(st_drop_geometry(segs_up) %>%
select(POI_ID, From_POI_ID, DAR, IncDA), by = c("COMID" = "POI_ID")) #%>%
#filter(IncDA.x < 10 | IncDA.y < 10)
#filter(IncDA.x < 10 | IncDA.y < 10)
# Filter by POI Type
Types <- c("Type_HUC12", "Type_Gages", "Type_TE", "Type_Conf", "Type_NID", "Type_WBIn", "Type_WBOut", "Type_TE")
Types <- Types[Types != exp]
......@@ -295,21 +514,24 @@ POI_move_down<-function(out_gpkg, POIs, Seg, Seg2, exp, DA_diff){
seg_sub <- seg_choices %>%
filter_at(Types, all_vars(is.na(.))) %>%
filter(!COMID %in% static_POIs$COMID) %>%
select (COMID, To_POI_ID, DAR.x, IncDA.x, From_POI_ID, DAR.y, DAR.y, IncDA.y, geom) %>%
st_as_sf()
select (COMID, To_POI_ID, DAR.x, IncDA.x, From_POI_ID, DAR.y, DAR.y, IncDA.y) %>%
st_sf()
st_write(seg_sub, out_gpkg, "Seg_choices", delete_layer = F, append = T)
# Commented for QAQCing
#st_write(seg_sub, out_gpkg, "Seg_choices", delete_layer = F, append = T)
segsub_dir <- seg_sub %>%
mutate(DirDA = ifelse(is.na(DAR.y), "Down",
# 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))
Move = ifelse(DirDA == "Down", To_POI_ID, From_POI_ID))
# POIs whose moves correspond to each others locations
Mult <- select(segsub_dir, COMID, Move) %>%
inner_join(st_drop_geometry(segsub_dir) %>%
select(COMID, Move), by = c("Move" = "COMID"), suffix = c(".A", ".B"))
inner_join(st_drop_geometry(.) %>%
select(COMID, Move), 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
......@@ -319,12 +541,12 @@ POI_move_down<-function(out_gpkg, POIs, Seg, Seg2, exp, DA_diff){
# 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
# 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")) %>%
......@@ -338,198 +560,32 @@ POI_move_down<-function(out_gpkg, POIs, Seg, Seg2, exp, DA_diff){
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(old_COMID = COMID.y) %>%
st_as_sf()
mutate(oldPOI = COMID.y, to_com = NA) %>%
st_sf()
# Bring new POI data over to new data
fin_POIs <- filter(POIs, !COMID %in% merged_POIs$old_COMID) %>%
left_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)
# Write to regional geopackage (x = new POI, y = old)
write_sf(segsub_dir, out_gpkg, paste0("POI_Change_", hydReg), append = T)
write_sf(merged_POIs, out_gpkg, paste0("MergedPOIs_", exp, "_", hydReg))
#***********************************************************************************************
# Maybe we don't need to this now?
# Raw nsegments (original resolution of NHDPlus flowlines)
nseg_raw <- Seg2 %>% left_join(select(st_drop_geometry(merged_POIs),c(COMID, old_COMID)),
by = c("POI_ID" = "old_COMID")) %>%
mutate(POI_ID = ifelse(!is.na(COMID.y), COMID.y, POI_ID)) %>%
rename(COMID = COMID.x) %>%
select(-c(COMID.y))
# Write to regional geopackage
#write_sf(nseg_raw, out_gpkg, paste0("nsegment_raw_", hydReg))
# Aggregate flowlines per POI_ID to a single segment, carrying over useful information
nsegments <- nseg_raw %>% group_by(POI_ID) %>%
mutate(PathTimeMA = na_if(PathTimeMA, -9999)) %>%
summarize(TotalLength = sum(LENGTHKM), TotalDA = max(TotDASqKM), HW = max(StartFlag),
TT = sum(PathTimeMA)) %>%
inner_join(st_drop_geometry(Seg2) %>%
select(COMID, Hydroseq, DnHydroseq), by=c("POI_ID"="COMID"))
changed_POIs <- POIs %>%
inner_join(select(st_drop_geometry(merged_POIs), COMID, oldPOI, to_com)) %>%
select(COMID, oldPOI, Type_HUC12, Type_Gages, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf, to_com) %>%
st_compatibalize(., read_sf(out_gpkg, poi_moved))
st_write(changed_POIs, out_gpkg, poi_moved, append = TRUE) # write_sf not appending?
# Format POI moves to table
xWalk <- select(st_drop_geometry(segsub_dir), Move, COMID) %>%
filter(!is.na(Move)) %>%
rename(COMID = Move, oldPOI = COMID)
st_write(xWalk, out_gpkg, poi_xWalk, append = TRUE) # write_sf not appending?
# produce a short data frame for populating TO_POI for downstream segment
to_from <- st_drop_geometry(Seg2) %>%
inner_join(st_drop_geometry(nseg_raw), by=c("DnHydroseq" = "Hydroseq")) %>%
select(COMID.x, Hydroseq, DnHydroseq, POI_ID.x, POI_ID.y) %>% rename(POI_ID = POI_ID.x, To_POI_ID = POI_ID.y)
# Add To_POI_ID to dissolved segments
nsegments <- nsegments %>%
left_join(to_from, by=c("POI_ID" = "COMID.x")) %>%
select(POI_ID, TotalLength, TotalDA, HW, TT, Hydroseq.x, To_POI_ID) %>%
rename(Hydroseq = Hydroseq.x)
newSegs <- segment_creation(mutate(Seg2, old_POI_ID = POI_ID), xWalk)
# Write out dissolved segments
#write_sf(nsegments, out_gpkg, paste0("nsegment_", hydReg))
# Write out new merged POIs
write_sf(fin_POIs, out_gpkg, paste0("allPOIs_", hydReg))
# Return POIs, segments, and raw segments
return (list(allPOIs = fin_POIs, segs = nsegments, segsRaw = nseg_raw))
}
DS_poiFix <- function(POIs, nhd){
##########################################
# Moves POI Upstream or downstream if it falls on COMID
# of flowline with no corresponding catchment
#
# Arguments:
# POIs : (data frame) POI data set
# nhd : (data frame) valid data frame of NHD flowlines
#
# Returns:
# repPOIs_unique : (data frame) DF of POIs with new COMID associated
##########################################
# Remove geom from NHD
nhdDF <- st_drop_geometry(nhd)
# Include downstream HS into data frame
tocomDF <- select(nhdDF, COMID, DnHydroseq) %>%
inner_join(select(nhdDF, COMID, Hydroseq), by = c("DnHydroseq" = "Hydroseq"))
# Find segments with POIs where there is no corresponding catchment
unCon_POIs <- filter(nhdDF, COMID %in% POIs$COMID, AreaSqKM == 0)
# apply initial move function to the whole POI Set
poi_fix <- filter(POIs, COMID %in% unCon_POIs$COMID) %>%
mutate(New_COMID = unlist(lapply(.$COMID, movePOI_NA_DA, nhdDF)), old_COMID = COMID)
# Convert to POI data frame format with appropriate fields to match input POIs
rep_POIs <- filter(nhd, COMID %in% poi_fix$New_COMID) %>%
get_node(., position = "end") %>%
mutate(COMID = (filter(nhd, COMID %in% poi_fix$New_COMID) %>%
pull(COMID))) %>%
left_join(poi_fix %>%
filter(!is.na(New_COMID)), by = c("COMID" = "New_COMID")) %>%
select(-COMID.y) %>%
st_drop_geometry()
# Fold in new POIs with existing POIs so all the "Type" attribution will carry over
rep_POIs_unique <- filter(POIs, COMID %in% rep_POIs$COMID) %>%
rbind(rep_POIs %>% select(-old_COMID)) %>%
group_by(COMID) %>%
summarise_each(funs(max(., na.rm=T)))
# Replace -Inf values that is an output of the 'summarise-each' application
rep_POIs_unique[mapply(is.infinite, rep_POIs_unique)] <- NA
# Get final list of POIs being moved, and the downstream POI they are linking to (many : 1 join)
rep_POIs_final <- rep_POIs_unique %>%
inner_join(poi_fix %>%
select(COMID, New_COMID), by = c("COMID" = "New_COMID")) %>%
inner_join(nhd %>%
select(COMID, DnHydroseq), by = c("COMID.y" = "COMID")) %>%
rename(new_COMID = COMID, old_COMID = COMID.y)
# Get the COMID of downstream flowlines of repPOIs_final still on flowline w/o catchment
more_replace <- filter(nhdDF, Hydroseq %in% rep_POIs_final$DnHydroseq, AreaSqKM == 0)
# Subset to produce replacement data frame of POIs that still need valid downstream COMID
move_down <- select(more_replace, COMID, DnHydroseq) %>%
inner_join(rep_POIs_final %>% select(new_COMID, old_COMID), by = c("COMID" = "old_COMID"))
# Column to hold the final toCOMID values
move_down$ToCom <- NA
# If the DS replacement flowlines with Inc Area > 0 equals the number of replacement comIDs
while(sum(is.na(move_down$ToCom)) > 0){
print(sum(is.na(move_down$ToCom)))
# DS replacement COMIDs
replacement_coms <- filter(st_drop_geometry(nhd), Hydroseq %in% move_down$DnHydroseq)#, AreaSqKM > 0)
# Add downstream COMIDS to the newCOMID
move_down <- move_down %>%
left_join(replacement_coms %>%
select(COMID, Hydroseq, DnHydroseq, AreaSqKM),
by = c("DnHydroseq" = "Hydroseq")) %>%
mutate(ToCom = ifelse(AreaSqKM > 0 & is.na(ToCom), COMID.y, ToCom), DnHydroseq = DnHydroseq.y) %>%
rename(COMID = COMID.x) %>%
select(COMID, DnHydroseq, new_COMID, ToCom, -c(COMID.y, AreaSqKM, DnHydroseq.y))
}
# Fold this information back in to the master data frame
rep_POIs_unique <- rep_POIs_unique %>%
left_join(move_down %>% select(new_COMID, ToCom), by = c("COMID" = "new_COMID")) %>%
right_join(poi_fix %>% select(New_COMID, old_COMID), by = c("COMID" = "New_COMID")) %>%
left_join(tocomDF, by = c("old_COMID" = "COMID.x")) %>%
mutate(ToCom = ifelse(!is.na(COMID) & is.na(ToCom), COMID.y, ToCom)) %>%
select(-c(DnHydroseq, COMID.y))
return (rep_POIs_unique)
}
movePOI_NA_DA <- function(poi_fix, nhd){
##########################################
# 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
#
# Arguments:
# poi_fix : (data frame) POI data set of COMIDs to be changed
# nhd : (data frame) valid data frame of NHD flowlines
#
# Returns:
# newPOI : (data frame) DF of POIs with new COMID associated
##########################################
# Closest POI/US/DS
up_segs <- nhdplusTools::get_UM(nhd, poi_fix, sort=T)
seg2fix <- filter(nhd, COMID == poi_fix)
# Sorted results and filter out all flowlines w/o catchments
upstuff <- filter(nhd, COMID %in% unlist(up_segs)) %>%
arrange(factor(COMID, levels = up_segs)) %>%
filter(AreaSqKM > 0)
down_segs <- nhdplusTools::get_DM(nhd, poi_fix, sort=T)
downstuff <- filter(nhd, COMID %in% unlist(down_segs)) %>%
arrange(factor(COMID, levels = down_segs)) %>%
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), downstuff %>%
filter(AreaSqKM > 0) %>%
select(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$COMID[(which.min(abs(seg2fix$TotDASqKM - near_FL$TotDASqKM)))]
} 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)
return (list(allPOIs = fin_POIs, segs = newSegs$diss_segs, FL = newSegs$raw_segs))
}
writePOIs <- function(POIs, out_gpkg, Type){
##########################################
# Write out final POI datasets with information
......@@ -567,7 +623,7 @@ writePOIs <- function(POIs, out_gpkg, Type){
}
structPOIsNet <- function(ncombined, nhd, final_POIs, out_gpkg){
structPOIsNet <- function(nhd, final_POIs){
##########################################
# Produce the structural POIs
#
......@@ -580,6 +636,10 @@ structPOIsNet <- function(ncombined, nhd, final_POIs, out_gpkg){
# Returns:
# writes Structural POIs and segments to geopackage
##########################################
# subset to flowlines used for nsegments
inc_segs <- filter(nhd, !is.na(POI_ID))
# Terminal outlets from the initial network
termout <- filter(nhd, !Hydroseq %in% DnHydroseq & TerminalFl == 1 & !COMID %in% final_POIs$COMID)
......@@ -587,21 +647,20 @@ structPOIsNet <- function(ncombined, nhd, final_POIs, out_gpkg){
out_POIs <- filter(nhd, COMID %in% final_POIs$COMID & TerminalFl == 1)
# Confluence POIs
conf_POIs <- filter(ncombined, COMID %in% final_POIs$COMID[final_POIs$Type_Conf == 1])
conf_POIs <- filter(inc_segs, COMID %in% final_POIs$COMID[final_POIs$Type_Conf == 1])
# Waterbody outlet POIs
wb_POIs <- filter(ncombined, COMID %in% final_POIs$COMID[!is.na(final_POIs$Type_WBOut) | !is.na(final_POIs$Type_WBIn)])
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_POIs <- get_node(struc_flines, position = "end") %>%
mutate(COMID = struc_flines$COMID, LevelPathI = struc_flines$LevelPathI) %>%
st_drop_geometry()
mutate(COMID = struc_flines$COMID, LevelPathI = struc_flines$LevelPathI)
# Add back in any levelpaths missing (shouldn't be much)
miss_LP <- filter(ncombined, COMID %in% final_POIs$COMID) %>%
miss_LP <- filter(inc_segs, COMID %in% final_POIs$COMID) %>%
filter(!LevelPathI %in% struc_POIs$LevelPathI)
# Only bind if there are rows present
......@@ -618,10 +677,10 @@ structPOIsNet <- function(ncombined, nhd, final_POIs, out_gpkg){
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(nhd), "up")))
final_net <- unique(unlist(lapply(unique(final_struct$COMID), NetworkNav, st_drop_geometry(nhd))))
# Subset NHDPlusV2 flowlines to navigation results and write to shapefile
structnet <- filter(nhd, COMID %in% final_net & grepl(paste0("^", hydReg, ".*"), .data$VPUID))
structnet <- filter(nhd, 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))
......
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