NHD_navigate.Rmd

---
title: "NHD Navigate"
output: html_document
---
This notebook Generate Segments and POIS from HUC12 outlets and GagesIII locations and builds
a minimally-sufficient stream network connecting all POis, as well as generating first-cut of
national segment network.

NEEDED MODS
1. HUC12 Mods (Missing HUC12, multiple oulets)
2. Decisions on NID
3. Waterbodies
```{r setup}
#  Load data and libraries ----------------------
 ptm<-proc.time()
library(nhdplusTools)
library(dplyr)
library(sf)
library(hyRefactor)

# Load custom functions
source("R/utils.R")
source("R/NHD_navigate.R")
# Load user-defined data paths
data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
# Load NHDPlus Data if precprocessed to RDS format
nhdplus_path(data_paths$nhdplus_gdb)
staged_nhd <- stage_national_data(nhdplus_data = data_paths$nhdplus_gdb, 
                                  output_path  = data_paths$nhdplus_dir)
# NHDPlus V2 Region
hydReg <-"01"

# Output Geopackage where we are writing output layers
out_gpkg <- file.path("cache", paste0("GF_",hydReg,".gpkg"))

# Output Layers
huc12_pois <- paste0("huc12pois_", hydReg) # HUC12 POIs
dup_comids <- paste0("cache/dupCOMIDs_",hydReg, ".rds") # data frame of COMIDs with multiple HUC12 assignments
nhdflowline <- "NHDFlowline"
gagesIII_pois <- paste0("gagesIIIpois_", hydReg) # GAGESIII POIs
gageLoc <- paste0("gageLoc_", hydReg) # GageLoc Files
TE_pois <- paste0("TEpois_", hydReg) # Thermoelectric POIs
NID_pois <- paste0("NID_", hydReg) # NID POIs
hucgage_segs <- paste0("hucgagesegs_", hydReg) # inimally-sufficient network
conf_pois <- paste0("confpois_", hydReg) # Confluence POIs
pois_all <- paste0("POIs_", hydReg) # All POIs binded together
nsegment_raw <- paste0("nsegment_raw_", hydReg) # Minimally-sufficient network attributed with POI_ID
n_segments <- paste0("Nsegment_", hydReg) # Minimally-sufficient network dissolved by POI_ID
```

```{r huc12_pois}
#  Derive or load HUC12 POIs ----------------------
if(needs_layer(out_gpkg, huc12_pois) | needs_layer(out_gpkg, nhdflowline)) {

  # Subset NHD by VPU
  nhd <- VPU_Subset(hydReg)
  write_sf(nhd, out_gpkg, nhdflowline)
  
  # Read in HUC12 outlets, get rid of duplicate COMID/HUC12 rows
  HUC12 <- read_sf(data_paths$hu12_points_path, "hu_points") %>% 
    filter(grepl(paste0("^", hydReg, ".*"), .data$HUC12)) %>% st_drop_geometry() %>%
    select(COMID, HUC12) %>% distinct()

  # Identify multiple COMIDs for HUC12
  dupCOMIDs<-HUC12 %>% group_by(COMID) %>% select(COMID, HUC12) %>% filter(n()>1)
  saveRDS(dupCOMIDs, dup_comids)
  
  # Return flowlinest that are HUC12 POIs
  huc12POIsegs <- nhd %>% filter(COMID %in% HUC12$COMID) 
  
  # Craete physical POIs from segments
  huc12POIs <- POI_creation(huc12POIsegs) %>%  
    #merge(st_drop_geometry(HUC12), by = "COMID") %>% #Need HUC12 mods before uncommenting
    mutate(Type_HUC12 = 1, Type_GagesIII = NA, Type_Conf = NA, Type_TE = NA, Type_NID = NA) %>% 
    select(COMID, geom, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
  
  # Write out geopackage layer representing POIs for given theme
  write_sf(huc12POIs, out_gpkg, huc12_pois)
} else {
  # Load HUC12 POIs if they already exist
  huc12POIs <- read_sf(out_gpkg, huc12_pois)   
  nhd <- read_sf(out_gpkg, nhdflowline)
}
```

```{r gagesIII_pois}
#  Derive or load GAGESIII POIs ----------------------
if(needs_layer(out_gpkg, gagesIII_pois)) {
  
  # Read in GAGESIII shapefile
  gagesIII <- read_sf(data_paths$gagesiii_points_path, "GagesIII_gages") %>% 
    filter(grepl(paste0("^", hydReg, ".*"), .data$NHDPlusReg)) %>% 
    mutate(COMID=as.integer(COMID))
  
  # Assign HUC12 POIs streamgage ID if they are co-located
  huc12POIs <- huc12POIs %>% merge(st_drop_geometry(gagesIII), all.x=TRUE) %>%
    mutate(Type_GagesIII = ifelse(!is.na(Gage_no), Gage_no, NA)) %>%
    select(COMID, geometry, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
  
  # Get NHD V2 flowlines that GAGESIII are assigned to, remove POIs already derived
  #     with HUC12 POIs
  gagesIIIsegs <- nhd %>% filter(COMID %in% gagesIII$COMID & !COMID %in% huc12POIs$COMID)
  
  # Derive GAGESIII POIs
  gagesIIIPOIs <- POI_creation(gagesIIIsegs) %>% 
    left_join(st_drop_geometry(gagesIII), by="COMID") %>% 
    mutate(Type_HUC12=NA, Type_GagesIII=Gage_no, Type_Conf=NA, Type_TE = NA, Type_NID = NA) %>%
    select(COMID, geom, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID) %>%
    rename(geometry=geom)
  
  # Write out geopackage layer representing POIs for given theme
  write_sf(gagesIIIPOIs, out_gpkg, gagesIII_pois)
  allPOIs <- do.call("rbind", list(huc12POIs, gagesIIIPOIs))
} else {
  # Load GAGESIII POIs if they already exist
  gagesIIIPOIs <- read_sf(out_gpkg, gagesIII_pois)   
  allPOIs <- do.call("rbind", list(huc12POIs, gagesIIIPOIs))
}
```

```{r TE_pois}
#  Derive or load Thermoelectric POIs ----------------------
if(needs_layer(out_gpkg, TE_pois)) {
  
  # Load POIs if they do not exist in workspace
  if(!exists("allPOIs")) allPOIs <- do.call("rbind", list(huc12POIs, gagesIIIPOIs))
  
  # Read in Thermoelectric shapefile
  TE_fac <- read_sf(data_paths$TE_points_path, "2015_TE_Model_Estimates_lat.long_COMIDs") %>% 
    mutate(COMID=as.integer(COMID)) %>% filter(COMID > 0) %>%
    inner_join(., select(st_drop_geometry(nhd), c(COMID, VPUID)), by = "COMID") %>%
    filter(grepl(paste0("^", hydReg, ".*"), .data$VPUID)) %>% select(COMID, VPUID) 
  
  # Get unique list of COMIDs
  TE_fac_unique <- TE_fac[!duplicated(TE_fac$COMID,fromLast=TRUE),]
  
  # Get NHD V2 flowlines that TE POIs are assigned to, remove POIs already derived
  #     with HUC12 and GAGESIII POIs
  TE_segs <- nhd %>% filter(COMID %in% TE_fac_unique$COMID & !COMID %in% allPOIs$COMID)
  
  # Assign existing POIs as a Type_TE if they are co-located
  allPOIs <- allPOIs %>% left_join(select(st_drop_geometry(TE_fac_unique), c(COMID, VPUID))) %>%
    mutate(Type_TE = ifelse(!is.na(VPUID), 1, NA))  %>%
    select(COMID, geometry, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
  
  # Derive TE POIs
  TE_POIs <- POI_creation(TE_segs) %>% rename(geometry = geom) %>%
    mutate(Type_HUC12=NA, Type_GagesIII=NA, Type_Conf=NA, Type_TE=1, Type_NID = NA) 
  
  # Write out geopackage layer representing POIs for given theme
  write_sf(TE_POIs, out_gpkg, TE_pois)
  allPOIs <- do.call("rbind", list(allPOIs, TE_POIs))
} else {
  # Load TE POIs if they already exist
  TE_POIs <- read_sf(out_gpkg, TE_pois)   
  allPOIs <- do.call("rbind", list(allPOIs, TE_POIs))  
}
```

```{r NID_pois}
#  Derive or load Thermoelectric POIs ----------------------
if(needs_layer(out_gpkg, NID_pois)) {
  
  # Load POIs if they do not exist in workspace
  if(!exists("allPOIs")) allPOIs <- do.call("rbind", list(huc12POIs, gagesIIIPOIs, TE_POIs))
  
  # Read in Thermoelectric shapefile
  NID_fac <- read.csv(paste0(data_paths$NID_points_path,"/NID_attributes_20170612.txt")) %>% 
    filter(ONNHDPLUS == 1, grepl(paste0("^",hydReg,".*"), .data$VPU))
  
  # Get unique list of COMIDs
  NID_fac <- NID_fac[!duplicated(NID_fac$FlowLcomid, fromLast=TRUE),]
  
  # Get NHD V2 flowlines that TE POIs are assigned to, remove POIs already derived
  #     with HUC12 and GAGESIII POIs
  NID_segs <- nhd %>% filter(COMID %in% NID_fac$FlowLcomid & !COMID %in% allPOIs$COMID)
  
  # Assign existing POIs as a Type_TE if they are co-located
  allPOIs <- allPOIs %>% left_join(select(NID_fac, c(FlowLcomid, VPU)),by = c("COMID" = "FlowLcomid")) %>%
    mutate(Type_NID = ifelse(!is.na(VPU), 1, NA))  %>% 
    select(COMID, geometry, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID) 
  
  # Derive TE POIs
  NID_POIs <- POI_creation(NID_segs) %>% rename(geometry = geom) %>%
    mutate(Type_HUC12=NA, Type_GagesIII=NA, Type_Conf=NA, Type_TE=NA, Type_NID=1) 
  
  # Write out geopackage layer representing POIs for given theme
  write_sf(NID_POIs, out_gpkg, NID_pois)
  allPOIs <- do.call("rbind", list(allPOIs, NID_POIs))
} else {
  # Load TE POIs if they already exist
  NID_POIs <- read_sf(out_gpkg, NID_pois)   
  allPOIs <- do.call("rbind", list(allPOIs, NID_POIs))  
}
```

```{r hucgages_segs}
# Build minimally-sufficient network to connect POIs ----------------------
if(needs_layer(out_gpkg, hucgage_segs)) {
  
  # Load POIs if they do not exist in workspace
  if(!exists("allPOIs")) allPOIs <- do.call("rbind", list(huc12POIs, gagesIIIPOIs, TE_POIs, NID_POIs))
  
  # Navigate upstream from each POI
  upNet <- unique(unlist(lapply(unique(allPOIs$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 NHDPlusV2 flowlines to navigation results and write to shapefile
  hucgagessegs <- nhd %>% filter(COMID %in% finalNet & grepl(paste0("^", hydReg, ".*"), .data$VPUID))
  
  # Write out minimally-sufficient network
  write_sf(hucgagessegs, out_gpkg, hucgage_segs) 
} else {
  # Load minimally-sufficient network
  hucgagessegs <- read_sf(out_gpkg, hucgage_segs)
}
```

```{r }
# Derive POIs at confluences where they are absent ----------------------
if(needs_layer(out_gpkg, conf_pois)) {
  
  # Load POIs if they do not exist in workspace
  if(!exists("allPOIs")) allPOIs <- do.call("rbind", list(huc12POIs, gagesIIIPOIs, TE_POIs, NID_POIs))
  
  # find confluences 
  confluences<-hucgagessegs %>% filter(DnHydroseq > 0) %>% group_by(DnHydroseq) %>% filter(n()>1) 
  
  # Mark existing POIs if they are a confluence
  allPOIs<-allPOIs %>% merge(st_drop_geometry(confluences), all.x=TRUE) %>%
    mutate(Type_Conf=ifelse(!is.na(RESOLUTION), 1, NA)) %>%
    select(COMID,geometry, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID)
    
  # Create new confluence POIs
  confPOIs<- POI_creation(confluences %>% filter(!COMID %in% allPOIs$COMID)) %>% 
    left_join(st_drop_geometry(confluences), by="COMID") %>% 
    mutate(Type_HUC12=NA, Type_GagesIII=NA, Type_Conf=1, Type_TE=NA, Type_NID=NA) %>%
    select(COMID, geom, Type_HUC12, Type_GagesIII, Type_Conf, Type_TE, Type_NID) %>%
    rename(geometry=geom) 
  
  # Write out shapefile representing POIs for given theme
  allPOIs<-rbind(allPOIs, confPOIs)
  write_sf(confPOIs, out_gpkg, conf_pois)
} else {
  confPOIs <- read_sf(out_gpkg, conf_pois)
  allPOIs <- do.call("rbind", list(huc12POIs, gagesIIIPOIs, TE_POIs, NID_POIs, confPOIs))
}
```

```{r }
# Derive first cut of segments ----------------------
if(needs_layer(out_gpkg, nsegment_raw)) {
  
  # Load POIs if they do not exist in workspace
  if(!exists("allPOIs")) allPOIs <- do.call("rbind", list(huc12POIs, gagesIIIPOIs, TE_POIs, NID_POIs, confPOIs))
  
  # Sort POIs by Levelpath and Hydrosequence in upstream to downstream order
  segPOIs <- st_drop_geometry(nhd) %>% filter(COMID %in% allPOIs$COMID) %>%
    arrange(desc(LevelPathI), desc(Hydroseq))
  hucgagessegs_ns <- st_drop_geometry(hucgagessegs)
  # Create column to hold  POI_ID values (most downstream COMID of each segment)
  hucgagessegs_ns$POI_ID<-0
  
  # Begin with upstream most Levelpath/Hydrosequence
  for (LP in unique(segPOIs$LevelPathI)){
    # For each level path
    segPOIs_LP <- segPOIs %>% filter(LevelPathI==LP)
    # Go through sorted list of US -> DS comids
    for (POI in segPOIs_LP$COMID){
      #print (POI)
      # Re-iterate only on flowlines that have no POI_ID to increase speed
      hucgagessegs_ns<-hucgagessegs_ns %>% filter(POI_ID==0)
      
      subPOI <- hucgagessegs_ns %>% filter (COMID == POI)
      
      # Get upstream segments
      Seg<-nhdplusTools::get_UM(hucgagessegs_ns, POI, include=T)
      # Assign POI_ID
      hucgagessegs_ns <- hucgagessegs_ns %>%  
        mutate(POI_ID = ifelse(COMID %in% Seg & POI_ID==0, POI, POI_ID))
      
      # Derive incremental segment composed of flowlines
      incSeg<-hucgagessegs_ns %>% filter(POI_ID==POI)
  
      #print (dim(incSeg))
      if(!exists("IncSegs")){
        IncSegs<-incSeg
      }else{
        IncSegs<-rbind(IncSegs, incSeg)
      }
    }
  }
  print(proc.time()-ptm)
  
  # Bring over NHDv2 attributes
  ncombined <- nhd %>% inner_join(select(IncSegs, c(COMID, POI_ID)), by = "COMID")
  
  # Write out shapefile representing POIs for given theme
  write_sf(allPOIs, out_gpkg, pois_all)
  # Write out nsegments composed of nhdflowlines
  write_sf(ncombined, out_gpkg, nsegment_raw)
  
  # Aggregate flowlines per POI_ID to a single segment, carrying over useful information
  nsegments<-ncombined %>% group_by(POI_ID) %>% mutate(PathTimeMA = na_if(PathTimeMA, -9999)) %>%
    summarize(TotalLength = sum(LENGTHKM), TotalDA = max(TotDASqKM), HW = max(StartFlag), 
              TT = sum(PathTimeMA)) %>% 
    inner_join(st_drop_geometry(hucgagessegs) %>% select(COMID, Hydroseq, DnHydroseq),   by=c("POI_ID"="COMID"))
  
  # produce a short data frame for populating TO_POI for downstream segment
  to_from<-segPOIs %>% inner_join(st_drop_geometry(ncombined), by=c("DnHydroseq" = "Hydroseq")) %>%
    select(COMID.x, Hydroseq, DnHydroseq, POI_ID) 
  
  # Add To_POI_ID to dissolved segments
  nsegments<-nsegments %>% left_join(to_from, by=c("POI_ID"="COMID.x")) %>% 
                                    mutate(To_POI_ID=POI_ID.y) %>%
    select(POI_ID,TotalLength,TotalDA,HW,TT,Hydroseq.x,To_POI_ID) %>% 
    rename(Hydroseq=Hydroseq.x)
  
  # Write out dissolved segments
  write_sf(nsegments, out_gpkg, n_segments)
} else {
  nsegmentraw <- read_sf(out_gpkg, nsegment_raw)
  nsegments <- read_sf(out_gpkg, n_segments)
  allPOIs <- read_sf(out_gpkg, pois_all)
}
```