Added waterbody inlet/outlet POIs

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)