NHD_navigate.Rmd

---
title: "NHD Navigate"
output: html_document
---

Project: GFv2.0  
Script purpose: Generate Segments and POIS from HUC12 outlets  
Date: 12/4/2019  
Author: [Andy Bock](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
1. Functionalize script with input arguments to cover workspace, results, etc.  
2. Add Gages to workflow  
3. Add NID/Waterbodies COMIDs to workflow  
4. Collapse Points close together  
5. Spit out segment length, travel time  
6. identify incremental catchments associated with each downstream POI  

```{r setup}
library(nhdplusTools)
library(dplyr)
library(sf)
library(hyRefactor)

source("R/utils.R")
source("R/NHD_navigate.R")
data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
nhdplus_path(data_paths$nhdplus_gdb)
staged_nhd <- stage_national_data(nhdplus_data = data_paths$nhdplus_gdb, 
                                  output_path = data_paths$nhdplus_dir)

# Output Geopackage
out_gpkg <- file.path("cache", "GF.gpkg")

# Output Layers
huc12_pois <- "huc12pois"
gagesIII_pois<-"gagesIIIpois"
huc12_segs <- "huc12segs"
hucgage_segs<-"hucgagesegs"
conf_pois <- "confpois"
huc12segs_all <- "huc12segs_all"
```

```{r huc12_pois}
if(needs_layer(out_gpkg, huc12_pois)) {
  HUC12 <- read_sf(data_paths$hu12_points_path, "hu_points") %>%
    filter(grepl("^06.*", .data$HUC12))
  
  # Read flowline simple features
  nhd<-readRDS(staged_nhd$flowline) 
  
  # Return flowlinest that are HUC12 POIs
  huc12POIsegs<-nhd %>% filter(COMID %in% HUC12$COMID) # 15 missing?
  
  # Sort data frame by highest-> lowest LP, then subsequently highest->lowest HS within each LP
  huc12POIsegs <- huc12POIsegs %>% arrange(desc(LevelPathI), desc(Hydroseq))
  
  #########################
  # Create POIs from HUC12 segments
  huc12POIs<-POI_creation(huc12POIsegs)
  # Write out shapefile representing POIs for given theme
  write_sf(huc12POIs, out_gpkg, huc12_pois)
} else {
  huc12POIs <- read_sf(out_gpkg, huc12_pois)   
}
```

```{r gagesIII_pois}
if(needs_layer(out_gpkg, gagesIII_pois)) {
  if(!exists("nhd")) nhd <- readRDS(staged_nhd$flowline) 
  
  gagesIII <- read_sf(data_paths$gagesiii_points_path, "GagesIII_gages") %>% 
    filter(NHDPlusReg == "06") %>%
    filter(!COMID %in% huc12POIs)
  
  # Return flowlinest that are HUC12 POIs
  gagesIIIsegs <- nhd %>% 
    filter(COMID %in% gagesIII$COMID) 
  
  # Sort data frame by highest -> lowest LP, then subsequently highest -> lowest HS within each LP
  gagesIIIsegs <- gagesIIIsegs %>% arrange(desc(LevelPathI), desc(Hydroseq))
  
  #########################
  # Create POIs from HUC12 segments
  gagesIIIPOIs <- POI_creation(gagesIIIsegs)
  
  # Write out shapefile representing POIs for given theme
  write_sf(gagesIIIPOIs, out_gpkg, gagesIII_pois)
} else {
  gagesIIIpois <- read_sf(out_gpkg, gagesIII_pois)   
}
```

```{r hucgages_segs}
# Navigate upstream along the mainstems from HUC12 POIS
if(needs_layer(out_gpkg, hucgage_segs)) {
  if(!exists("nhd")) nhd <- readRDS(staged_nhd$flowline)
  
  POIs <- rbind(huc12POIs, gagesIIIPOIs)
  
  # Navigate upstream from each POI
  upNet <- unique(unlist(lapply(POIs$COMID,
                                NetworkNav,
                                st_drop_geometry(nhd), "up")))
  
  # Connect POIs that don't have connecting levelpath looking downstream and 
  #         produce unique set of flowlines linking POIs
  finalNet <- unique(NetworkConnection(upNet, st_drop_geometry(nhd)))
  
  # Subset nhd segs to navigation results and write to shapefile
  hucgagessegs <- nhd %>% filter(COMID %in% finalNet)
  
  write_sf(hucgagessegs, out_gpkg, hucgage_segs)
} else {
  hucgagessegs <- read_sf(out_gpkg, hucgage_segs)
}
```

```{r }
if(needs_layer(out_gpkg, conf_pois)) {
  # Derive confluence POIs where they are absent
  nhdConf<-nhd %>% filter(COMID %in% hucgagessegs$COMID) # check on redundancy here
  
  # find confluences 
  confluences<-nhdConf %>% group_by(DnHydroseq) %>% filter(n()>1) %>% 
    filter(!COMID %in% huc12POIsegs$COMID)
  
  # Create POIs from HUC12 segments
  confPOIs<-POI_creation(confluences)
  
  # Write out shapefile representing POIs for given theme
  write_sf(confPOIs, out_gpkg, conf_pois)
} else {
  confPOIs <- read_sf(out_gpkg, conf_pois)
}
```

```{r eval = FALSE, echo = FALSE}
##############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 %>% 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)
# }
```

```{r}
# if(needs_layer(out_gpkg, huc12segs_all)) {
#   #################################################################
#   # Assign POI_ID to each flowline related to POI/Incremental catchment
#   #  Right now this is not working as I intended, need to 
#   #        debug on some areas that are not working (see POI_ID 14647681)
#   #  Maybe an issue with global assignment?
#   
#   # Column to hold the unique ID of flowlines locally-related to downstream POI
#   if(!exists("nhd")) {
#     # Doing this for global assignment in NetworkNav
#     nhd <- readRDS(staged_nhd$flowline)
#     nhdDF <- st_drop_geometry(nhd)
#   }
#   
#   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))
#   # Get flowlines (proto-segments) associated with each COMID
#   upNet<-unique(unlist(lapply(POIs_ordered$COMID,NetworkNav,nhdDF,"up")))
#   
#   # Subset nhd segs to navigation results and write to shapefile
#   allSegs<-nhdDF %>% filter(COMID %in% upNet)
#   
#   allSegs<-allSegs %>% 
#     inner_join(nhdDF, by=c("COMID","COMID")) %>%
#     select(COMID,LevelPathI.x,Hydroseq.x,POI_ID.x) %>%
#     left_join(select(nhd, COMID), by = "COMID")
#   
#   write_sf(allSegs, out_gpkg, huc12segs_all)
# } else {
#   allSegs <- read_sf(out_gpkg, huc12segs_all)
# }
```