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)

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"
huc12_segs <- "huc12segs"
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 huc12_segs}
# Navigate upstream along the mainstems from HUC12 POIS
if(needs_layer(out_gpkg, huc12_segs)) {
  if(!exists("nhd")) nhd<-readRDS(staged_nhd$flowline)
  
  # Navigate upstream from each POI
  upNet<-unique(unlist(lapply(HUC12$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
  huc12Segs<-nhd %>% filter(COMID %in% finalNet)
  write_sf(huc12Segs, out_gpkg, huc12_segs)
} else {
  huc12Segs <- read_sf(out_gpkg, huc12_segs)
}
```

```{r }
if(needs_layer(out_gpkg, conf_pois)) {
  # Derive confluence POIs where they are absent
  nhdConf<-nhd %>% filter(COMID %in% huc12Segs$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}
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)
}
```

```{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)
# }
```