diff --git a/workspace/NHD_navigate.Rmd b/workspace/NHD_navigate.Rmd
index 156e20d2af866eb17e9e9ba3b0a1bf541c5772e6..d2e8c54115c8d2de028a694a1599e2b024b2f4bd 100644
--- a/workspace/NHD_navigate.Rmd
+++ b/workspace/NHD_navigate.Rmd
@@ -6,10 +6,6 @@ This notebook Generate Segments and POIS from HUC12 outlets and GagesIII locatio
a minimally-sufficient stream network connecting all POis, as well as generating first-cut of
national segment network.
-NEEDED MODS
-1. HUC12 Mods (Missing HUC12, multiple oulets)
-2. Decisions on NID
-3. Waterbodies
```{r setup}
# Load data and libraries ----------------------
ptm<-proc.time()
@@ -28,7 +24,7 @@ nhdplus_path(data_paths$nhdplus_gdb)
staged_nhd <- stage_national_data(nhdplus_data = data_paths$nhdplus_gdb,
output_path = data_paths$nhdplus_dir)
# NHDPlus V2 Region
-hydReg <-"01"
+hydReg <-"03"
# Output Geopackage where we are writing output layers
out_gpkg <- file.path("cache", paste0("GF_",hydReg,".gpkg"))
@@ -37,6 +33,9 @@ out_gpkg <- file.path("cache", paste0("GF_",hydReg,".gpkg"))
huc12_pois <- paste0("huc12pois_", hydReg) # HUC12 POIs
dup_comids <- paste0("cache/dupCOMIDs_",hydReg, ".rds") # data frame of COMIDs with multiple HUC12 assignments
nhdflowline <- "NHDFlowline"
+WBs <- paste0("WB_", hydReg)
+WBout_POIs <- paste0("WBout_", hydReg)
+WBin_POIs <- paste0("WBin_", hydReg)
gagesIII_pois <- paste0("gagesIIIpois_", hydReg) # GAGESIII POIs
gageLoc <- paste0("gageLoc_", hydReg) # GageLoc Files
TE_pois <- paste0("TEpois_", hydReg) # Thermoelectric POIs
@@ -57,7 +56,7 @@ if(needs_layer(out_gpkg, huc12_pois) | needs_layer(out_gpkg, nhdflowline)) {
write_sf(nhd, out_gpkg, nhdflowline)
# Read in HUC12 outlets, get rid of duplicate COMID/HUC12 rows
- HUC12 <- read_sf(data_paths$hu12_points_path, "hu_points") %>%
+ HUC12 <- read_sf(data_paths$hu12_points_path, "hu_outlets") %>%
filter(grepl(paste0("^", hydReg, ".*"), .data$HUC12)) %>% st_drop_geometry() %>%
select(COMID, HUC12) %>% distinct()
@@ -71,8 +70,8 @@ if(needs_layer(out_gpkg, huc12_pois) | needs_layer(out_gpkg, nhdflowline)) {
# Craete physical POIs from segments
huc12POIs <- POI_creation(huc12POIsegs) %>%
#merge(st_drop_geometry(HUC12), by = "COMID") %>% #Need HUC12 mods before uncommenting
- mutate(Type_HUC12 = 1, Type_GagesIII = NA, Type_Conf = NA, Type_TE = NA, Type_NID = NA) %>%
- select(COMID, geom, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
+ mutate(Type_HUC12 = 1, Type_WBIn = NA, Type_WBOut = NA, Type_GagesIII = NA, Type_Conf = NA, Type_TE = NA, Type_NID = NA) %>%
+ select(COMID, geometry, Type_HUC12, Type_WBIn, Type_WBOut, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
# Write out geopackage layer representing POIs for given theme
write_sf(huc12POIs, out_gpkg, huc12_pois)
@@ -83,6 +82,51 @@ if(needs_layer(out_gpkg, huc12_pois) | needs_layer(out_gpkg, nhdflowline)) {
}
```
+```{r waterbody_pois}
+# Derive or load GAGESIII POIs ----------------------
+if(needs_layer(out_gpkg, WBin_POIs)) {
+
+ # Read in waterobdies
+ wb <- readRDS("data/NHDPlusNationalData/nhdplus_waterbodies.rds")
+ # Get Waterbodies that intersect with Region flowlines
+ WBs_VPU <- wb %>% filter(COMID %in% nhd$WBAREACOMI)
+ # Write out waterbodies
+ write_sf(WBs_VPU, out_gpkg, WBs)
+
+ # Segments that are in waterbodies
+ wbSegs <- nhd %>% filter(WBAREACOMI > 0)
+ wbOut <- nhd %>% filter(WBAREACOMI %in% WBs_VPU$COMID) %>% group_by(WBAREACOMI) %>% slice(which.min(Hydroseq))
+ wbIn <- nhd %>% filter(!WBAREACOMI %in% WBs_VPU$COMID) %>% filter(DnHydroseq %in% wbSegs$Hydroseq)
+
+ # Assign HUC12 POIs WB in/out if they are co-located
+ huc12POIs <- huc12POIs %>% merge(st_drop_geometry(wbIn), all.x=TRUE) %>%
+ mutate(Type_WBIn = ifelse(!is.na(GNIS_ID), 1, NA)) %>%
+ select(COMID, geometry, Type_HUC12, Type_WBIn, Type_WBOut, Type_GagesIII, Type_Conf, Type_TE, Type_NID) %>%
+ merge(st_drop_geometry(wbOut), all.x=TRUE) %>%
+ mutate(Type_WBOut = ifelse(!is.na(GNIS_ID), 1, NA)) %>%
+ select(COMID, geometry, Type_HUC12, Type_WBIn, Type_WBOut, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
+
+ # Craete physical POIs from segments
+ WBoutPOIs <- POI_creation(wbOut %>% filter(!COMID %in% huc12POIs$COMID)) %>%
+ mutate(Type_HUC12 = NA, Type_WBIn = NA, Type_WBOut = 1, Type_GagesIII = NA, Type_Conf = NA, Type_TE = NA, Type_NID = NA) %>%
+ select(COMID, geometry, Type_HUC12, Type_WBIn, Type_WBOut, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
+
+ WBinPOIs <- POI_creation(wbIn %>% filter(!COMID %in% huc12POIs$COMID)) %>%
+ mutate(Type_HUC12 = NA, Type_WBIn = 1, Type_WBOut = NA, Type_GagesIII = NA, Type_Conf = NA, Type_TE = NA, Type_NID = NA) %>%
+ select(COMID, geometry, Type_HUC12, Type_WBIn, Type_WBOut, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
+
+ allPOIs <- do.call("rbind",list(huc12POIs, WBinPOIs, WBoutPOIs))
+ write_sf(WBoutPOIs, out_gpkg, WBout_POIs)
+ write_sf(WBinPOIs, out_gpkg, WBin_POIs)
+} else {
+ # Load HUC12 POIs if they already exist
+ WBinPOIs <- read_sf(out_gpkg, WBin_POIs)
+ WBoutPOIs <- read_sf(out_gpkg, WBout_POIs)
+ allPOIs <- do.call("rbind",list(huc12POIs, WBinPOIs, WBoutPOIs))
+ nhd <- read_sf(out_gpkg, nhdflowline)
+}
+```
+
```{r gagesIII_pois}
# Derive or load GAGESIII POIs ----------------------
if(needs_layer(out_gpkg, gagesIII_pois)) {
@@ -93,9 +137,9 @@ if(needs_layer(out_gpkg, gagesIII_pois)) {
mutate(COMID=as.integer(COMID))
# Assign HUC12 POIs streamgage ID if they are co-located
- huc12POIs <- huc12POIs %>% merge(st_drop_geometry(gagesIII), all.x=TRUE) %>%
+ allPOIs <- allPOIs %>% merge(st_drop_geometry(gagesIII), all.x=TRUE) %>%
mutate(Type_GagesIII = ifelse(!is.na(Gage_no), Gage_no, NA)) %>%
- select(COMID, geometry, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
+ select(COMID, geometry, Type_HUC12, Type_WBIn, Type_WBOut, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
# Get NHD V2 flowlines that GAGESIII are assigned to, remove POIs already derived
# with HUC12 POIs
@@ -104,17 +148,16 @@ if(needs_layer(out_gpkg, gagesIII_pois)) {
# Derive GAGESIII POIs
gagesIIIPOIs <- POI_creation(gagesIIIsegs) %>%
left_join(st_drop_geometry(gagesIII), by="COMID") %>%
- mutate(Type_HUC12=NA, Type_GagesIII=Gage_no, Type_Conf=NA, Type_TE = NA, Type_NID = NA) %>%
- select(COMID, geom, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID) %>%
- rename(geometry=geom)
+ mutate(Type_HUC12=NA, Type_WBIn = NA, Type_WBOut = NA, Type_GagesIII=Gage_no, Type_Conf=NA, Type_TE = NA, Type_NID = NA) %>%
+ select(COMID, geometry, Type_HUC12, Type_WBIn, Type_WBOut, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
# Write out geopackage layer representing POIs for given theme
write_sf(gagesIIIPOIs, out_gpkg, gagesIII_pois)
- allPOIs <- do.call("rbind", list(huc12POIs, gagesIIIPOIs))
+ allPOIs <- do.call("rbind", list(allPOIs, gagesIIIPOIs))
} else {
# Load GAGESIII POIs if they already exist
gagesIIIPOIs <- read_sf(out_gpkg, gagesIII_pois)
- allPOIs <- do.call("rbind", list(huc12POIs, gagesIIIPOIs))
+ allPOIs <- do.call("rbind", list(allPOIs, gagesIIIPOIs))
}
```
@@ -141,11 +184,11 @@ if(needs_layer(out_gpkg, TE_pois)) {
# Assign existing POIs as a Type_TE if they are co-located
allPOIs <- allPOIs %>% left_join(select(st_drop_geometry(TE_fac_unique), c(COMID, VPUID))) %>%
mutate(Type_TE = ifelse(!is.na(VPUID), 1, NA)) %>%
- select(COMID, geometry, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
+ select(COMID, geometry, Type_HUC12, Type_WBIn, Type_WBOut, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
# Derive TE POIs
- TE_POIs <- POI_creation(TE_segs) %>% rename(geometry = geom) %>%
- mutate(Type_HUC12=NA, Type_GagesIII=NA, Type_Conf=NA, Type_TE=1, Type_NID = NA)
+ TE_POIs <- POI_creation(TE_segs) %>%
+ mutate(Type_HUC12=NA, Type_WBIn=NA, Type_WBOut=NA, Type_GagesIII=NA, Type_Conf=NA, Type_TE=1, Type_NID = NA)
# Write out geopackage layer representing POIs for given theme
write_sf(TE_POIs, out_gpkg, TE_pois)
@@ -178,11 +221,11 @@ if(needs_layer(out_gpkg, NID_pois)) {
# Assign existing POIs as a Type_TE if they are co-located
allPOIs <- allPOIs %>% left_join(select(NID_fac, c(FlowLcomid, VPU)),by = c("COMID" = "FlowLcomid")) %>%
mutate(Type_NID = ifelse(!is.na(VPU), 1, NA)) %>%
- select(COMID, geometry, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
+ select(COMID, geometry, Type_HUC12, Type_WBIn, Type_WBOut, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
# Derive TE POIs
- NID_POIs <- POI_creation(NID_segs) %>% rename(geometry = geom) %>%
- mutate(Type_HUC12=NA, Type_GagesIII=NA, Type_Conf=NA, Type_TE=NA, Type_NID=1)
+ NID_POIs <- POI_creation(NID_segs) %>%
+ mutate(Type_HUC12=NA, Type_WBIn=NA, Type_WBOut=NA, Type_GagesIII=NA, Type_Conf=NA, Type_TE=NA, Type_NID=1)
# Write out geopackage layer representing POIs for given theme
write_sf(NID_POIs, out_gpkg, NID_pois)
@@ -233,14 +276,13 @@ if(needs_layer(out_gpkg, conf_pois)) {
# Mark existing POIs if they are a confluence
allPOIs<-allPOIs %>% merge(st_drop_geometry(confluences), all.x=TRUE) %>%
mutate(Type_Conf=ifelse(!is.na(RESOLUTION), 1, NA)) %>%
- select(COMID,geometry, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
+ select(COMID,geometry, Type_HUC12, Type_WBIn, Type_WBOut, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
# Create new confluence POIs
confPOIs<- POI_creation(confluences %>% filter(!COMID %in% allPOIs$COMID)) %>%
left_join(st_drop_geometry(confluences), by="COMID") %>%
- mutate(Type_HUC12=NA, Type_GagesIII=NA, Type_Conf=1, Type_TE=NA, Type_NID=NA) %>%
- select(COMID, geom, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID) %>%
- rename(geometry=geom)
+ mutate(Type_HUC12=NA, Type_WBIn=NA, Type_WBOut=NA, Type_GagesIII=NA, Type_Conf=1, Type_TE=NA, Type_NID=NA) %>%
+ select(COMID, geometry, Type_HUC12, Type_WBIn, Type_WBOut, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
# Write out shapefile representing POIs for given theme
allPOIs<-rbind(allPOIs, confPOIs)