First commit

NHD_navigate.R

+##################################################
+## Project: GFv2.0
+## Script purpose: Generate Segments and POIS from HUC12 outlets
+## Date: 12/4/2019
+## Author: abock@usgs.gov
+## Notes:  First cut, includes both HUC12 outlets as POIS and c
+##         connecting confluences; Run over Hydrologic Region 06
+## NEEDED MODS
+#  0 - a) create setup script for downloading and formatting data (commented out sections below)
+#      b) downloading data from SB into the proper formats
+#      c) ensuring the proper libraries are installed and loaded.
+#  1 - Functionalize script with input arguments to cover workspace, results, etc.
+#  2 - Add Gages to workflow
+#  3 - Add NID/Waterbodies COMIDs to workflow
+#  3 - Collapse Points close together
+#  4 - Spit out segment length, travel time
+#  5 - identify incremental catchments associated with each downstream POI
+##################################################
+# Network navigation for upstream/downstream from a COMID of interest
+NetworkNav<-function(inCom,nhdDF,navType){
+  if (navType=="up"){
+    Seg<-get_UM(nhdDF,inCom,include=T)
+  }else{
+    Seg<-get_DM(nhdDF,inCom,distance=1050,include=T)
+  }
+  return(Seg)
+}
+
+library(nhdplusTools)
+library(tidyverse)
+library(sf)
+library(rgdal)
+
+#https://cran.r-project.org/web/packages/nhdplusTools/nhdplusTools.pdf
+setwd("D:/abock/NHM_FY20/gfv2.0")
+
+#######################
+# Read Region 06 HUC12s
+HUC12<-readOGR("data/hu_points.gpkg","hu_points")
+HUC12@data$Region<-substr(HUC12@data$HUC12,1,2)
+comIDs<-HUC12@data %>% filter(Region == "06") %>% select (COMID)
+
+# Set path to national NHDPlus v2 seamless data; 
+Print ("This will take awhile, go take a walk")
+if (!file.exists("data/nhdplus_flowline_attributes.rds")){ 
+  download_nhdplusv2("data")
+  nhdplus_path("data/NHDPlusNationalData/NHDPlusV21_National_Seamless_Flattened_Lower48.gdb")
+  stage_national_data(include=c("attribute","flowline"),output_path="data") # Get catchments later?
+}else{
+  print ("NHDPlusv2 data downloaded and staged")
+}
+
+# Read flowline attribute data frames
+nhdDF<-readRDS("data/nhdplus_flowline_attributes.rds") # %>% select(COMID, Pathlength, LENGTHKM, LevelPathI, UpHydroseq, Hydroseq, DnHydroseq)
+# Read flowline simple features
+nhd<-readRDS("data/nhdplus_flowline.rds") 
+# Return flowlinest that are POIs
+huc12POIsegs<-nhd %>% filter(COMID %in% comIDs[,1]) # 15 missing?
+
+#########################
+# break segs it to points
+Segpts <- st_geometry(huc12POIsegs) %>% 
+  lapply(., function(x) {st_sfc(x) %>% st_cast(., 'POINT')})
+# subset the last point from each geometry, make a POINT sf object
+Seg_ends <- sapply(Segpts, function(p) {
+  p[c(length(p))]}) %>% st_sfc() %>%st_sf('geom' = .)
+# Write out shapefile representing POIs for given theme
+st_write(Seg_ends,dsn="workspace/huc12POIs.shp",delete_dsn = T)
+write_rds(Seg_ends,path="workspace/huc12POIS.rds")
+
+##########################
+# Navigate upstream along the mainstems
+upNet<-unique(unlist(lapply(comIDs,NetworkNav,nhdDF,"up")))
+# Subset by dangling segments that need downstream navigation to be connected
+upNet_DF<-nhd %>% filter(COMID %in% upNet) %>% 
+  filter(!DnHydroseq %in% Hydroseq)
+# Nagivate downstream along the mainstems
+downNet<-unique(unlist(lapply(upNet_DF$COMID,NetworkNav,nhdDF,"down")))
+# Combine upstream/downstream navigation results
+huc12Segs<-unique(c(upNet,downNet))
+
+# Subset nhd segs to navigation results and write to shapefile
+huc12Segs<-nhd %>% filter(COMID %in% huc12Segs)
+st_write(huc12Segs,dsn="workspace/huc12segs.shp",delete_dsn = T)
+write_rds(Seg_ends,path="workspace/huc12segs.rds")
+
+# Derive confluence POIs where they are absent
+nhdConf<-nhd %>% filter(COMID %in% huc12Segs$COMID)
+confluences<-nhdConf %>% group_by(DnHydroseq) %>% filter(n()>1) %>% filter(!COMID %in% huc12POIsegs$COMID)
+  
+confpts <- st_geometry(confluences) %>% 
+  lapply(., function(x) {st_sfc(x) %>% st_cast(., 'POINT')})
+# subset the last point from each geometry, make a POINT sf object
+conf_ends <- sapply(confpts, function(p) {
+  p[c(length(p))]}) %>% st_sfc() %>%st_sf('geom' = .)
+# Write out shapefile representing POIs for given theme
+st_write(conf_ends,dsn="results/confPOIs.shp",delete_dsn = T)
+write_rds(conf_ends,path="workspace/confPOIs.rds")
+
+##########################################
+##############DEPRECATED SECTION##################
+
+# # retrieve POI from seg (this is deprecated)
+# retrievePOI<-function(nhdPOI){
+#   #print (iter)
+#   #linepnts<-st_sfc(nhdPOI)
+#   linepnts<-st_cast(nhdPOI,"POINT")
+#   print(linepnts)
+#   endPoint<-linepnts[length(linepnts)]
+#   st_cast(endPoint,'POINT')
+# }