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 
##         connecting confluences; Run over Hydrologic Region 06
##         in this example
## NEEDED MODS
#  0 - a) create setup script for downloading and formatting data 
#      b) 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"){
    # Upstream network navigation
    Seg<-nhdplusTools::get_UM(nhdDF,inCom,include=T)
  }else{
    # Downstream network navigation (distance currently hard-coded)
    Seg<nhdplusTools::get_DM(nhdDF,inCom,distance=1050,include=T)
  }
  return(Seg)
}

# POI Creation
POI_creation<-function(inSegs){
  # break segs it to points
  Segpts <- st_geometry(inSegs) %>% 
    lapply(., function(x) {sf::st_sfc(x) %>% sf::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))]}) %>% sf::st_sfc() %>% sf::st_sf('geom' = .)
  return(Seg_ends)
}

library(nhdplusTools) #https://cran.r-project.org/web/packages/nhdplusTools/nhdplusTools.pdf
library(tidyverse) # a collection of libraries that work with eacher other (dplyr)
library(sf) # simple features - a library for handling spatial data
library(rgdal) # rgdal - an open source GIS library that's more akin to ESRI and handles shapefile

# Set working directory
setwd("D:/abock/NHM_FY20/gfv2.0")

#######################
# Read Region 06 HUC12 outlets, get this data at https://www.sciencebase.gov/catalog/item/5dbc53d4e4b06957974eddae
HUC12<-rgdal::readOGR("data/hu_points.gpkg","hu_points")
# Add new region field
HUC12@data$Region<-substr(HUC12@data$HUC12,1,2)
# Get the COMIDS of HUC12 outlets within Region 06
comIDs<-HUC12@data %>% dplyr::filter(Region == "06") %>% dplyr::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 seamless data to this location
  nhdplusTools::download_nhdplusv2("data")
  # Set the path
  nhdplusTools::nhdplus_path("data/NHDPlusNationalData/NHDPlusV21_National_Seamless_Flattened_Lower48.gdb")
  # Stage the data as RDS for easy access
  nhdplusTools::stage_national_data(include=c("attribute","flowline","catchment"),output_path="data") # Get catchments later?
}else{
  print ("NHDPlusv2 data already downloaded and staged")
}

# Read flowline attribute data frames for entire seamless database
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 %>% dplyr::filter(COMID %in% comIDs[,1]) # 15 missing?

#########################
# Create POIs from HUC12 segments
huc12POIs<-POI_creation(huc12POIsegs)
# Write out shapefile representing POIs for given theme
sf::st_write(huc12POIs,dsn="workspace/huc12POIs.shp",delete_dsn = T)
write_rds(huc12POIs,path="workspace/huc12POIS.rds")

##########################
# Navigate upstream along the mainstems from HUC12 POIS
upNet<-unique(unlist(lapply(comIDs,NetworkNav,nhdDF,"up")))
# Subset by dangling segments that need downstream navigation to be connected
#        reduces the number of features we need to connect
upNet_DF<-nhd %>% dplyr::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)
sf::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 %>% dplyr::filter(COMID %in% huc12Segs$COMID) # check on redundancy here

confluences<-nhdConf %>% dplyr::group_by(DnHydroseq) %>% dplyr::filter(n()>1) %>% 
  dplyr::filter(!COMID %in% huc12POIsegs$COMID)
  
# Create POIs from HUC12 segments
huc12POIs<-POI_creation(confluences)

# Write out shapefile representing POIs for given theme
sf::st_write(conf_ends,dsn="workspace/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')
# }