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

Added structural POIs func., cleaned up others and added comments on input/output

parent 90c217db
No related branches found
No related tags found
1 merge request!31Bock
Hide whitespace changes
Inline Side-by-side
workspace/R/NHD_navigate.R
+ 236
− 80
View file @ 975155e8
......@@ -9,7 +9,7 @@ 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
gagesIII_pois <- paste0("GagesIII_", hydReg) # GAGESIII 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
......@@ -19,20 +19,38 @@ 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
# Network navigation for upstream/downstream from a COMID of interest
#
NetworkNav <- function(inCom, nhdDF, navType){
#print (inCom)
# Upstream network navigation is easier
##########################################
# Network navigation for upstream/downstream from a COMID of interest
#
# Arguments:
# inCOM : (list) list of input COMIDs that are 'dangles'
# nhdDF : (data frame) valid data frame of NHD flowlines
#
# Returns:
# mydata : (list) list of COMIDs to connect dangle to network
##########################################
Seg <- nhdplusTools::get_UM(nhdDF, inCom, include = TRUE)
return(Seg)
}
# Network navigation to connect dangles in the network
NetworkConnection <- function(inCom, nhdDF){
# Subset by dangling segments that need downstream navigation to be connected
# reduces the number of features we need to connect
upNet_DF <- nhdDF %>% dplyr::filter(COMID %in% inCom) %>%
NetworkConnection <- function(inCom, nhd){
##########################################
# Connects dangles in the network that are not
# terminal outlets
#
# Arguments:
# inCOM : (list) list of input COMIDs that are 'dangles'
# nhdDF : (data frame) valid data frame of NHD flowlines
#
# Returns:
# 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) %>%
filter(!DnHydroseq %in% Hydroseq)
# while the number of dangles is greater than 0
......@@ -44,13 +62,13 @@ NetworkConnection <- function(inCom, nhdDF){
DSLP <- upNet_DF$DnLevelPat[upNet_DF$COMID %in% inCom]
# Get the COMID of the hydroseq with level path value
# the lowest downstream flowline within the levelpath
inCom2 <- nhdDF$COMID[nhdDF$Hydroseq %in% DSLP]
inCom2 <- nhd$COMID[nhd$Hydroseq %in% DSLP]
# Run the upstream navigation code
upNet <- unique(unlist(lapply(inCom2, NetworkNav, nhdDF, "up")))
upNet <- unique(unlist(lapply(inCom2, NetworkNav, nhd, "up")))
# Append result to existing segment list
inCom <- append(inCom, upNet)
# Get the same variable as above
upNet_DF<-nhdDF %>% dplyr::filter(COMID %in% inCom) %>%
upNet_DF<-nhd %>% dplyr::filter(COMID %in% inCom) %>%
filter(!DnHydroseq %in% Hydroseq)
# Get the count
count2 <- dim(upNet_DF)[1]
......@@ -63,9 +81,19 @@ NetworkConnection <- function(inCom, nhdDF){
return(inCom)
}
# Swith divergence path for POIs from minor (2) to major (1)
# reduces the number of short, spurious segments and POIS as locations such as waterbody outlets
switchDiv <- function(inSegs, nhd){
##########################################
# Swith divergence path for POIs from minor (2) to major (1)
# reduces the number of short, spurious segments and POIS as locations such as waterbody outlets
#
# Arguments:
# inSegs : (list) list of input COMIDs that are 'dangles'
# nhd : (data frame) valid data frame of NHD flowlines
#
# Returns:
# mydata : (list) list of COMIDs for major diversions
##########################################
divSegs <- nhd %>% filter(COMID %in% inSegs$COMID)
if (2 %in% divSegs$Divergence){
print ("Switching divergence to other fork")
......@@ -86,54 +114,84 @@ switchDiv <- function(inSegs, nhd){
return(inSegs)
}
# POI Creation
POI_creation<-function(inSegs, nhd){
# break segs into points
Segpts <- st_geometry(inSegs) %>%
lapply(., function(x) {sf::st_sfc(x) %>% sf::st_cast(., 'POINT')})
POI_creation<-function(srcData, nhd, IDfield){
##########################################
# Creates POIs for a given data theme
#
# Arguments:
# srcData : (data frame) DF with two columns:
# 1 - COMID
# 2 - Unique ID value for POI (Streamgage ID, etc.)
# nhd : (data frame) valid data frame of NHD flowlines
# IDfield : (character) Name of 'Type' field to be modified in POI
#
# Returns:
# mydata : (simple features) POIs for the specific data theme
##########################################
# subset the last point from each geometry, make a POINT sf object
Seg_ends <- sapply(Segpts, function(p) {
p[c(length(p))]}) %>%
sf::st_sfc(crs = st_crs(inSegs)) %>%
sf::st_sf('geom' = .)
# Give generic ID to Identity field
colnames(srcData) <- c("COMID", "ID")
# Add the COMID, otherwise comes out as feature with no relating attributes
Seg_ends$COMID<-inSegs$COMID
if("geom" %in% colnames(Seg_ends)) Seg_ends <- Seg_ends %>% rename(geometry=geom)
return(Seg_ends)
subSegs <- nhd %>% filter(COMID %in% srcData$COMID)
#merge(st_drop_geometry(subSegs %>% select(COMID))
POIs <-subSegs %>% get_node(., position = "end") %>% mutate(COMID = subSegs$COMID) %>%
mutate(Type_HUC12 = NA, Type_WBIn = NA, Type_WBOut = NA, Type_Gages = NA, Type_TE = NA, Type_NID = NA, Type_Conf = NA) %>%
inner_join(srcData %>% select(COMID, ID), by = "COMID") %>% mutate(!!(paste0("Type_", IDfield)) := ID) %>%
select(COMID, Type_HUC12, Type_Gages, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf, geometry)
#if("geom" %in% colnames(Seg_ends)) Seg_ends <- Seg_ends %>% rename(geometry=geom)
return(POIs)
}
WB_area_width <- function(polygon){
addType<-function(POIs, newPOIs, IDfield){
##########################################
# Checks if existing POIs are co-located with new POI set
# Adds the type attribute for co-located POIs of multiple themes
#
# Arguments:
# POIs : (data frame) Existing POIs
# newPOIs: (data frame) newly-derived POIs of a given data theme
# IDfield : (character) Name of 'Type' field to be modified in POI
#
# Returns:
# mydata : (data frame) Existing POIs with Type field modified for
# duplicate POIs for given data theme
##########################################
POIs <- POIs %>% rename(ID = !!(paste0("Type_", IDfield)))
medWidth <- median(st_distance(st_cast(st_geometry(polygon),"POINT")))
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)) %>% rename(!!(paste0("Type_", IDfield)) := ID) %>% select(-c(ID2))
return (medWidth)
return(POIs_fin)
}
POI_move_down<-function(out_gpkg, pois, Seg, Seg2, exp){
POI_move_down<-function(out_gpkg, POIs, Seg, Seg2, exp){
##########################################
# Need to update & review this;
# seems really long and could be shortened
##########################################
# Segs to collapse downstream (fold POIs together downstream if valid catchment downstream)
segs_Down <- Seg %>% inner_join(select(st_drop_geometry(.),c(POI_ID, TotalDA, TotalLength)),
by=c("To_POI_ID" = "POI_ID")) %>%
inner_join(select(st_drop_geometry(poisAll),
c(COMID, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID)),by=c("POI_ID"="COMID")) %>%
segs_Down <- Seg %>% inner_join(select(st_drop_geometry(.), c(POI_ID, TotalDA, TotalLength)),
by = c("To_POI_ID" = "POI_ID")) %>%
inner_join(select(st_drop_geometry(POIs),
c(COMID, Type_HUC12, Type_Gages, Type_Conf, Type_TE, Type_NID)), by = c("POI_ID" = "COMID")) %>%
# Select POIs within the correct drainage area ratio and fit the size criteria
mutate(DAR = TotalDA.x/TotalDA.y, IncDA = TotalDA.y - TotalDA.x) %>%
# If the drainage area ratio is within acceptable range
filter(DAR < 1 & DAR >= 0.95 | TotalLength.y < 1, IncDA > 0)
# Filter by POI Type
Types <- c("Type_HUC12", "Type_GagesIII", "Type_TE", "Type_Conf", "Type_NID", "Type_TE")
Types <- c("Type_HUC12", "Type_Gages", "Type_TE", "Type_Conf", "Type_NID", "Type_TE")
Types <- Types[Types != exp]
# Divide to segments of interest based on POI Type
# Only move gages downstream if they are upstream of confluence
if (exp %in% c("Type_GagesIII", "Type_TE")){
if (exp %in% c("Type_Gages", "Type_TE")){
ConfPOIs <- pois %>% filter(Type_Conf == 1)
SegSub <- segs_Down %>% filter_at(Types, all_vars(is.na(.))) %>%
filter(To_POI_ID %in% ConfPOIs$COMID) %>%
......@@ -144,14 +202,14 @@ POI_move_down<-function(out_gpkg, pois, Seg, Seg2, exp){
}
#1 - POIs that need to be moved downstream
MoveDownPOI <- pois %>% filter(COMID %in% SegSub$POI_ID)
MoveDownPOI <- POIs %>% filter(COMID %in% SegSub$POI_ID)
# POIs that are stationary and their infomration will be merged with upstream POIs
stationaryPOI <- pois %>% filter(!COMID %in% SegSub$POI_ID)
stationaryPOI <- POIs %>% filter(!COMID %in% SegSub$POI_ID)
# Final Set to be merged with that don't involve either sets above
FinalPOI <- pois %>% filter(!COMID %in% SegSub$POI_ID & !COMID %in% SegSub$To_POI_ID) %>%
FinalPOI <- POIs %>% filter(!COMID %in% SegSub$POI_ID & !COMID %in% SegSub$To_POI_ID) %>%
mutate(merged_COMID = NA)
#2 - Join SegSub assignment to pois to bring over POI attributes
#2 - Join SegSub assignment to POIs to bring over POI attributes
MoveDownPOI_withAtt <- MoveDownPOI %>%
inner_join(st_drop_geometry(SegSub), by = c("COMID" = "POI_ID"), suffix = c(".x", ".y")) %>%
select(-c(geom, TotalLength.y, DAR, IncDA))
......@@ -162,22 +220,24 @@ POI_move_down<-function(out_gpkg, pois, Seg, Seg2, exp){
by = c("COMID" = "To_POI_ID"), suffix = c(".x", ".y")) %>%
# Bring over relevant attributes
mutate(Type_HUC12 = ifelse(!is.na(Type_HUC12.y), 1, Type_HUC12.x)) %>%
mutate(Type_GagesIII = ifelse(is.na(Type_GagesIII.x) & !is.na(Type_GagesIII.y), Type_GagesIII.y, NA)) %>%
mutate(Type_Gages = ifelse(is.na(Type_Gages.x) & !is.na(Type_Gages.y), Type_Gages.y, NA)) %>%
# Gages_B is incase there are two gages being merged together, not writing out for now
mutate(Type_GagesIII_B = ifelse(!is.na(Type_GagesIII.y), Type_GagesIII.y, NA)) %>%
mutate(Type_Gages_B = ifelse(!is.na(Type_Gages.y), Type_Gages.y, NA)) %>%
mutate(Type_Conf = ifelse(!is.na(Type_Conf.y), Type_Conf.y, Type_Conf.x)) %>%
mutate(Type_TE = ifelse(!is.na(Type_TE.y), Type_TE.y, Type_TE.x)) %>%
mutate(Type_NID = ifelse(!is.na(Type_NID.y), Type_NID.y, Type_NID.x)) %>%
mutate(merged_COMID = COMID.y) %>%
select(COMID, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID, merged_COMID, geom)
select(COMID, Type_HUC12, Type_Gages, Type_Conf, Type_TE, Type_NID, merged_COMID, geom)
# Write to regional geopackage
write_sf(MergedPOIs, 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(MergedPOIs),c(COMID, merged_COMID)),
nseg_raw <- Seg2 %>% left_join(select(st_drop_geometry(MergedPOIs),c(COMID, merged_COMID)),
by = c("POI_ID" = "merged_COMID")) %>%
mutate(POI_ID = ifelse(!is.na(COMID.y), COMID.y, POI_ID)) %>% rename(COMID = COMID.x) %>% select (-c(COMID.y))
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))
......@@ -194,8 +254,8 @@ POI_move_down<-function(out_gpkg, pois, Seg, Seg2, exp){
# 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)
select(POI_ID, TotalLength, TotalDA, HW, TT, Hydroseq.x, To_POI_ID) %>%
rename(Hydroseq = Hydroseq.x)
# Write out dissolved segments
write_sf(nsegments, out_gpkg, paste0("nsegment_", hydReg))
......@@ -208,28 +268,42 @@ POI_move_down<-function(out_gpkg, pois, Seg, Seg2, exp){
# Adjust confluence POIs based on if they are associated with flowline with no
# catchment/IncDA == 0
DS_poiFix <- function(tmpPOIs, nhd){
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
##########################################
# Have subset nhd with no geometry
# Remove geom from NHD
nhdDF <- st_drop_geometry(nhd)
# Include downstream HS into data frame
ToComDF <- nhdDF %>% select(COMID, DnHydroseq) %>%
inner_join(st_drop_geometry(nhd) %>% select(COMID, Hydroseq), by = c("DnHydroseq" = "Hydroseq"))
# Find segments with POIs where there is no corresponding catchment
unConPOIs <- nhdDF %>% filter(COMID %in% tmpPOIs$COMID, AreaSqKM == 0)
poiFix <- tmpPOIs %>% filter(COMID %in% unConPOIs$COMID) %>%
mutate(New_COMID = unlist(lapply(.$COMID, movePOI_NA_DA, nhdDF)), old_COMID = COMID)
unConPOIs <- nhdDF %>% filter(COMID %in% POIs$COMID, AreaSqKM == 0)
# apply initial move function to the whole POI Set
poiFix <- POIs %>% filter(COMID %in% unConPOIs$COMID) %>%
mutate(New_COMID = unlist(lapply(.$COMID, movePOI_NA_DA, nhdDF)), old_COMID = COMID)
# Convert to POI data frame format with appropriate fields
repPOIs <- nhd %>% filter(COMID %in% poiFix$New_COMID) %>% POI_creation(., nhd) %>%
# Convert to POI data frame format with appropriate fields to match input POIs
repPOIs <- nhd %>% filter(COMID %in% poiFix$New_COMID) %>% get_node(., position = "end") %>%
mutate(COMID = (nhd %>% filter(COMID %in% poiFix$New_COMID) %>% pull(COMID))) %>%
left_join(poiFix %>% 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
repPOIs_unique <- tmpPOIs %>% filter(COMID %in% repPOIs$COMID) %>% rbind(repPOIs %>% select(-old_COMID)) %>%
repPOIs_unique <- POIs %>% filter(COMID %in% repPOIs$COMID) %>% rbind(repPOIs %>% select(-old_COMID)) %>%
group_by(COMID) %>% summarise_each(funs(max(., na.rm=T)))
# Replace -Inf values
# Replace -Inf values that is an output of the 'summarise-each' application
repPOIs_unique[mapply(is.infinite, repPOIs_unique)] <- NA
# Get final list of POIs being moved, and the downstream POI they are linking to (many : 1 join)
......@@ -269,23 +343,34 @@ DS_poiFix <- function(tmpPOIs, nhd){
return (repPOIs_unique)
}
# Move POIs that fall on flowlines with no catchment upstream/downstream
# to adjacent flowline with most similar total drainage area
movePOI_NA_DA <- function(poiFix, nhdDF){
# by most reandom to least random
movePOI_NA_DA <- function(poiFix, 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:
# poiFix : (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
upSegs <- nhdplusTools::get_UM(nhdDF, poiFix, sort=T)
seg2fix <-nhdDF %>% filter(COMID == poiFix)
upSegs <- nhdplusTools::get_UM(nhd, poiFix, sort=T)
seg2fix <-nhd %>% filter(COMID == poiFix)
# Sorted results and filter out all flowlines w/o catchments
upStuff <- nhdDF %>% filter(COMID %in% unlist(upSegs)) %>% arrange(factor(COMID, levels=upSegs)) %>%
upStuff <- nhd %>% filter(COMID %in% unlist(upSegs)) %>% arrange(factor(COMID, levels=upSegs)) %>%
filter(AreaSqKM > 0)
downSegs <- nhdplusTools::get_DM(nhdDF, poiFix, sort=T)
downStuff <- nhdDF %>% filter(COMID %in% unlist(downSegs)) %>% arrange(factor(COMID, levels=downSegs)) %>%
downSegs <- nhdplusTools::get_DM(nhd, poiFix, sort=T)
downStuff <- nhd %>% filter(COMID %in% unlist(downSegs)) %>% arrange(factor(COMID, levels=downSegs)) %>%
filter(AreaSqKM > 0)
# combine into one dataframe
# combine into one dataframe, select up/downstream seg with least change in total drainage area
nearFL <- rbind(upStuff %>% select(COMID, TotDASqKM, AreaSqKM) %>% slice(1), downStuff %>% filter(AreaSqKM > 0) %>%
select(COMID, TotDASqKM, AreaSqKM) %>% slice(1))
......@@ -293,6 +378,7 @@ movePOI_NA_DA <- function(poiFix, nhdDF){
if (sum(nearFL$AreaSqKM) > 0){
newPOI <- nearFL$COMID[(which.min(abs(seg2fix$TotDASqKM - nearFL$TotDASqKM)))]
} else {
# Remove POI if not catchment associated with flowlines upstream or downstream
print (poiFix)
print ("US and DS flowlines also have no catchment, removing POI")
newPOI <- NA
......@@ -300,26 +386,96 @@ movePOI_NA_DA <- function(poiFix, nhdDF){
return(newPOI)
}
# Write out final POI datasets with information
# used to generate segments
writePOIs <- function(inDF, out_gpkg, Type){
writePOIs <- function(POIs, out_gpkg, Type){
##########################################
# Write out final POI datasets with information
#
# Arguments:
# POIs : (data frame) POI data set
# out_gkpg : (geopackage) Geopackage where final POI layers written
# Type : (character) Type of POI being written; default is write features for all types
#
# Returns:
# finPOIs : (data frame) DF of final POIs
##########################################
print ("Writing out final POIs")
# If type is missing write out all flowlines
if (missing(Type)){
print ("Writing out all POI Types")
lyrs <- st_layers(out_gpkg)
# get subcategory of POIs
POIs <- lyrs$name[lyrs$geomtype == "Point"]
for (poi in POIs[-(grep("POIs_*", POIs))]){
POInames <- lyrs$name[!is.na(lyrs$geomtype) & lyrs$geomtype== "Point"]
print (POInames)
for (poi in POInames){
print (poi)
poiType <- unlist(strsplit(poi, "_"))[1]
subPOIs <- inDF %>%
subPOIs <- POIs %>%
filter(!is.na(!!as.name(paste0("Type_", poiType))))
write_sf(subPOIs, out_gpkg, poi)
}
} else {
subPOIs <- inDF %>%
subPOIs <- POIs %>%
filter(!is.na(!!as.name(paste0("Type_", Type))))
write_sf(subPOIs, out_gpkg, poi)
}
return(inDF)
}
structPOIsNet <- function(ncombined, nhd, finalPOIs, out_gpkg){
##########################################
# Produce the structural POIs
#
# Arguments:
# ncombined : (data frame) final Segments
# nhd : (data frame) valid data frame of NHD flowlines
# finalPOIs : (data frame) final POIs
# out_gkpg : (geopackage) Geopackage where final POI layers written
# Type : (character) Type of POI being written; default is write features for all types
#
# Returns:
# writes Structural POIs and segments to geopackage
##########################################
# Terminal outlets from the initial network
termOut <- nhd %>% filter(!Hydroseq %in% DnHydroseq & TerminalFl == 1 & !COMID %in% finalPOIs$COMID)
# POIs that are also terminal outlets
outPOIs <- nhd %>% filter(COMID %in% finalPOIs$COMID & TerminalFl == 1)
# Confluence POIs
confPOIs <- ncombined %>% filter(COMID %in% finalPOIs$COMID[finalPOIs$Type_Conf == 1])
# Waterbody outlet POIs
WBPOIs <- ncombined %>% filter(COMID %in% finalPOIs$COMID[!is.na(finalPOIs$Type_WBOut) | !is.na(finalPOIs$Type_WBIn)])
# Waterbody inlets POIs
strucFlines <- termOut %>% bind_rows(outPOIs, confPOIs, WBPOIs) %>% arrange(COMID)
strucPOIs <- get_node(strucFlines, position = "end") %>% mutate(COMID = strucFlines$COMID, LevelPathI = strucFlines$LevelPathI) %>%
st_drop_geometry()
# Add back in any levelpaths missing (shouldn't be much)
miss_LP <- ncombined %>% filter(COMID %in% finalPOIs$COMID) %>% filter(!LevelPathI %in% strucPOIs$LevelPathI)
# Only bind if there are rows present
if (nrow(miss_LP) > 0){
LP_pois <- c(miss_LP$LevelPathI, strucPOIs$LevelPathI)
} else {
LP_pois <- strucPOIs$LevelPathI
}
# final stuctural POIs, order by COMID to match with strucPOIs
finalStruct <- nhd %>% filter(Hydroseq %in% LP_pois) %>% arrange(COMID)
structPOIs <- get_node(finalStruct, position = "end") %>% mutate(COMID = finalStruct$COMID)
# produce unique set of flowlines linking POIs
finalNet <- unique(unlist(lapply(unique(finalStruct$COMID), NetworkNav, st_drop_geometry(nhd), "up")))
# Subset NHDPlusV2 flowlines to navigation results and write to shapefile
StructNet <- nhd %>% filter(COMID %in% finalNet & grepl(paste0("^", hydReg, ".*"), .data$VPUID))
# Write out minimally-sufficient network
write_sf(structPOIs, out_gpkg, "struct_POIs2")
write_sf(StructNet, out_gpkg, "struct_Net2")
}
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