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){
  # Upstream network navigation is easier
  Seg<-nhdplusTools::get_UM(nhdDF,inCom,include=T)
  # Assign POI_ID to flowlines associated with each POI/Incremental Segment
  #nhdDF$POI_ID[nhdDF$COMID %in% Seg & nhdDF$POI_ID==0]<<-inCom
  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) %>% 
    filter(!DnHydroseq %in% Hydroseq)
  
  # while the number of dangles is greater than 0
  while (length(upNet_DF$COMID)>0){
    # create item for number of dangles
    count=dim(upNet_DF)[1]
    print (dim(upNet_DF))
    # find out which level paths are downstream of dangling huc12 POIs
    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<-nhd$COMID[nhd$Hydroseq %in% DSLP]
    # Run the upstream navigation code
    upNet<-unique(unlist(lapply(inCom2,NetworkNav,nhdDF,"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) %>% 
      filter(!DnHydroseq %in% Hydroseq)
    # Get the count
    count2=dim(upNet_DF)[1]
    # if the count has remained the same we are done and return the flowline list
    if (count == count2){
      return (inCom)
    }
  }
  # Not sure this other return is needed
  return(inCom)
}

# POI Creation
POI_creation<-function(inSegs){
  # break segs into 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' = .)
  # Add the COMID, otherwise comes out as feature with no relating attributes
  Seg_ends$COMID<-inSegs$COMID
  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") 
}else{
  print ("NHDPlusv2 data already downloaded and staged, reading in flowlines information; will do catchment work in another step")
  # 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 HUC12 POIs
huc12POIsegs<-nhd %>% dplyr::filter(COMID %in% comIDs[,1]) # 15 missing?
# Sort data frame by highest-> lowest LP, then subsequently highest->lowest HS within each LP
huc12POIsegs<-huc12POIsegs %>% dplyr::arrange(desc(LevelPathI),desc(Hydroseq))

#########################
# 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

# Navigate upstream from each POI
upNet<-unique(unlist(lapply(comIDs[,1],NetworkNav,nhdDF,"up")))
# Connect POIs that don't have connecting levelpath looking downstream and 
#         produce unique set of flowlines linking POIs
finalNet<-unique(NetworkConnection(upNet,nhdDF))

# Subset nhd segs to navigation results and write to shapefile
huc12Segs<-nhd %>% filter(COMID %in% finalNet)
sf::st_write(huc12Segs,dsn="workspace/huc12segs_all.shp",delete_dsn = T)
write_rds(huc12Segs,path="workspace/huc12segs_US.rds")

##########################
# Derive confluence POIs where they are absent
nhdConf<-nhd %>% dplyr::filter(COMID %in% huc12Segs$COMID) # check on redundancy here

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

# Write out shapefile representing POIs for given theme
sf::st_write(confPOIs,dsn="workspace/confPOIs.shp",delete_dsn = T)
write_rds(confPOIs,path="workspace/confPOIs.rds")

#################################################################
# Assign POI_ID to each flowline related to POI/Incremental catchment
# Column to hold the unique ID
nhdDF$POI_ID<-0
# Combine POI Types
allPOIs<-rbind(huc12POIs,confPOIs)
# sort data frame by
POIs_ordered<-nhdDF %>% filter(COMID %in% allPOIs$COMID) %>% arrange(desc(LevelPathI,Hydroseq))

# from here go into funciton that gives ID to each segment that relates it to POI_ID

##############DEPRECATED/In Development SECTION##################

# segs<-nhdDF %>% filter(POI_ID > 0)
# nhdSegs<-nhd %>% inner_join(segs,by=c('COMID'='COMID')) %>% select(COMID,POI_ID,Shape)

# # retrieve POI from seg (this is deprecated)
# retrievePOI<-function(nhdPOI)
# # Network navigation for upstream/downstream from a COMID of interest
# NetworkNav<-function(inCom,nhdDF,navType){
#   if (navType=="up"){
#     # Upstream network navigation is easier
#     Seg<-nhdplusTools::get_UM(nhdDF,inCom,include=T)
#     # Just assign POI_ID to flowlines that are 0 and upstream
#     #nhdDF$POI_ID[nhdDF$COMID %in% Seg & nhdDF$POI_ID==0]<<-inCom
#     
#     # 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)
#     
#   }else{
#     # Downstream network navigation (distance currently hard-coded)
#     #Seg<-nhdplusTools::get_DM(nhdDF,inCom,distance=1050,include=T)
#     # subset all segs not equal to 0, assign the lowest HS value
#     #lowHS<-min(nhdDF$Hydroseq[nhdDF$COMID %in% Seg & nhdDF$POI_ID==0])
#     # Get the COMID of the lowest HS value
#     #comID<-nhdDF$COMID[nhdDF$Hydroseq==lowHS]
#     # Assign these segs the COMID
#     #nhdDF$POI_ID[nhdDF$COMID %in% Seg & nhdDF$POI_ID==0]<<-comID
#   }
#   return(Seg)
# }