Formalized AK navigate

workspace/02_AK_navigate.Rmd

+---
+title: "NHD Navigate"
+output: html_document
+editor_options:
+  chunk_output_type: console
+---
+This notebook Generate Segments and POIS from POI data sources and builds
+a minimally-sufficient stream network connecting all POis, as well as generating first-cut of
+national segment network for Alaska (VPU 19)
+
+```{r}
+knitr::opts_chunk$set(error = TRUE)
+```
+
+```{r setup}
+#  Load data and libraries ----------------------
+library(nhdplusTools)
+library(dplyr)
+library(sf)
+library(hyRefactor)
+library(tidyr)
+library(hyfabric)
+
+VPU <- "19"
+
+# Load custom functions
+source("R/utils.R")
+source("R/NHD_navigate.R")
+source("R/ExCONUS_config.R")
+source("R/hyRefactor_config.R")
+
+# Load user-defined data paths
+data_paths <- jsonlite::read_json(file.path("cache", "19_data_paths.json"))
+```
+
+```{r huc12 POIs}
+#  Derive or load HUC12 POIs ----------------------
+# Still need to deal with some duplicate HUC12s
+if(needs_layer(nav_gpkg, temp_POIs)) {
+
+  if(needs_layer(flowline)){
+    # Subset NHD by VPU
+    nhd <- read_sf(data_paths$ak_hydro, flowline) %>% 
+    # Convert to albers
+    st_transform(crs) %>%
+    # Identify terminal flowlines
+    mutate(TerminalFl = ifelse(Hydroseq == TerminalPa, 1, 0), 
+           poi = ifelse(denTotalAreaSqKM >= min_da_km, 1, 0))
+    
+    write_sf(nhd, nav_gpkg, flowline)
+  } else {
+    nhd <- read_sf(nav_gpkg, flowline)
+  }
+
+  # Join Dave's HUC12 outlets with NHD
+  HUC12_COMIDs <- read_sf(data_paths$ak_hydro, "wbd_points") %>% 
+    # select just HUC12 guys
+    select(COMID, hu12) %>%
+    # Convert to albers
+    st_transform(crs) %>%
+    # Raname to generic for POI creation function
+    rename(ID = hu12) %>%
+    # There is some duplicity to the HUC12 COMIDs 
+    #  This slims it down to 1 record per COMID
+    group_by(COMID) %>% 
+    slice(1)
+
+  # Create POIs - some r05 HUC12 POIs not in R05 NHD
+  huc12_POIs <- POI_creation(st_drop_geometry(HUC12_COMIDs), filter(nhd, poi == 1), "HUC12")
+
+  # Write out geopackage layer representing POIs for given theme
+  write_sf(huc12_POIs, nav_gpkg, temp_POIs)
+  tmp_POIs <- huc12_POIs
+} else {
+  # Load HUC12 POIs as the tmpPOIs if they already exist
+  tmp_POIs <- read_sf(nav_gpkg, temp_POIs) 
+  nhd <- read_sf(nav_gpkg, flowline)
+}
+```
+
+```{r Gages POIs}
+if(all(is.na(tmp_POIs$Type_Gages))) {
+  
+  # the POI dataset
+  # Note that alot of HUC12 are missing from above
+  pois <- read_sf(data_paths$ak_poi, "pois") %>%
+    filter(!uniqueID %in% tmp_POIs$Type_HUC12,
+           type != "HUC12")
+  
+  ak_cats <- read_sf(data_paths$ak_hydro, "AK_cats") %>%
+    st_transform(crs)
+  
+  # Determine which merrit catchment non-HUC12 pois fall within
+  poi_types <- filter(pois) %>%
+    mutate(COMID = as.numeric(ak_cats[sapply(st_intersects(., ak_cats), `length<-`, 1), ]$COMID)) %>%
+    mutate(COMID = as.numeric(ifelse(is.na(.$COMID), 
+           ak_cats[st_nearest_feature(filter(., is.na(COMID)), ak_cats),]$COMID, COMID))) %>%
+    filter(!grepl("[a-zA-Z]", .data$uniqueID), nchar(uniqueID) ==8 ) %>%
+    st_drop_geometry() %>%
+    select(COMID, uniqueID)
+  
+  # Derive ML POIS
+  tmp_POIs <- POI_creation(poi_types, nhd %>% mutate(COMID = as.numeric(COMID)), "Gages") %>%
+    group_by(COMID) %>%
+    slice(1) %>%
+    ungroup() %>%
+    addType(., tmp_POIs, "Gages")
+    
+  write_sf(tmp_POIs, nav_gpkg, temp_POIs)
+} else {
+  # Load HUC12 POIs as the tmpPOIs if they already exist
+  tmp_POIs <- read_sf(nav_gpkg, temp_POIs) 
+  nhd <- read_sf(nav_gpkg, flowline)
+}
+```
+
+
+```{r WU/ML POIs}
+if(all(is.na(tmp_POIs$Type_WU))) {
+  
+  # the POI dataset
+  pois <- read_sf(data_paths$ak_poi, "pois") %>%
+    filter(!uniqueID %in% c(tmp_POIs$Type_HUC12, tmp_POIs$Type_Gages))
+  
+  ak_cats <- read_sf(data_paths$ak_hydro, "AK_cats") %>%
+    st_transform(crs)
+  
+  # Determine which merrit catchment non-HUC12 pois fall within
+  poi_types <- filter(pois) %>%
+    mutate(COMID = as.numeric(ak_cats[sapply(st_intersects(., ak_cats), `length<-`, 1), ]$COMID)) %>%
+    mutate(COMID = as.numeric(ifelse(is.na(.$COMID), 
+           ak_cats[st_nearest_feature(filter(., is.na(COMID)), ak_cats),]$COMID, COMID))) %>%
+    st_drop_geometry()
+  
+  # Derive ML POIS
+  tmp_POIs <- POI_creation(filter(poi_types, type == "ML") %>% select(COMID, uniqueID, -type), 
+                          nhd %>% mutate(COMID = as.numeric(COMID)), "ML") %>%
+    group_by(COMID) %>%
+    slice(1) %>%
+    ungroup() %>%
+    addType(., tmp_POIs, "ML")
+  
+  # Derive WU POIs
+  tmp_POIs <- POI_creation(filter(poi_types, type == "WU") %>% select(COMID, uniqueID, -type), 
+                          nhd %>% mutate(COMID = as.numeric(COMID)), "WU") %>%
+    group_by(COMID) %>%
+    slice(1) %>%
+    ungroup() %>%
+    mutate(Type_ML = NA) %>%
+    addType(., tmp_POIs, "WU")
+    
+  write_sf(tmp_POIs, nav_gpkg, temp_POIs)
+} else {
+  # Load HUC12 POIs as the tmpPOIs if they already exist
+  tmp_POIs <- read_sf(nav_gpkg, temp_POIs) 
+  nhd <- read_sf(nav_gpkg, flowline)
+}
+```
+
+```{r Confluence POIs}
+# Derive POIs at confluences where they are absent ----------------------
+if((needs_layer(nav_gpkg, pois_all))) {
+
+  # Gets the levelpahs of the HUC12 POIs
+  nhd_POIs <- filter(nhd, COMID %in% tmp_POIs$COMID)
+  
+  # NHD derived level path
+  nhd_LP <- nhd %>% 
+    filter(LevelPathI %in% nhd_POIs$LevelPathI)
+  
+  # NHD populated with uphydroseqeuence
+  nhd_Full <- nhd %>%
+    left_join(select(st_drop_geometry(nhd_LP), COMID, LevelPathI, Hydroseq, DnHydroseq), by = c("Hydroseq" = "DnHydroseq")) %>%
+    mutate(fromCOMID = ifelse(is.na(COMID.y), 0, COMID.y),
+           UpHydroseq = ifelse(is.na(Hydroseq.y), 0, Hydroseq.y)) %>%
+    rename(COMID = COMID.x, LevelPathI = LevelPathI.x) %>%
+    select (-c(COMID.y, Hydroseq.y)) %>%
+    group_by(COMID) %>%
+    filter(n() == 1 | (n() > 1 & LevelPathI == LevelPathI.y)) %>%
+    ungroup() %>%
+    select(-LevelPathI.y) 
+  
+  # Navigate upstream from each POI and determine minimally-sufficient network between current POI sets
+  up_net <- unique(unlist(lapply(unique(tmp_POIs$COMID), NetworkNav, nhd_Full))) 
+  finalNet <- unique(NetworkConnection(up_net, st_drop_geometry(nhd_Full))) 
+
+  # Subset NHDPlusV2 flowlines to navigation results and write to shapefile
+  nhd_Full <- mutate(nhd_Full, minNet = ifelse(COMID %in% finalNet, 1, 0))
+  
+  # Create new confluence POIs
+  conf_COMIDs <- st_drop_geometry(filter(nhd_Full, minNet == 1)) %>%
+    #st_drop_geometry(nhd_Full) %>%
+    # Downstream hydrosequence of 0 indicates
+    #   the flowline is terminating or
+    #   leaving the domain, so they
+    #   are excluded from this process
+    filter(DnHydroseq > 0) %>%
+    group_by(DnHydroseq) %>%
+    filter(n()> 1) %>%
+    mutate(Type_Conf = LevelPathI) %>%
+    ungroup() %>% 
+    select(COMID, Type_Conf) %>%
+    distinct()
+
+  tmp_POIs <- POI_creation(conf_COMIDs, filter(nhd_Full, minNet == 1), "Conf") %>%
+    addType(., tmp_POIs, "Conf")
+  
+  write_sf(nhd_Full, nav_gpkg, flowline)
+  write_sf(tmp_POIs, nav_gpkg, pois_all)
+} else {
+  nhd_Full <- read_sf(nav_gpkg, flowline)
+  pois_all <- read_sf(nav_gpkg, pois_all)
+}
+```
+
+```{r Final segments}
+# Derive first cut of segments ----------------------
+if(needs_layer(nav_gpkg, n_segments)) {
+
+  # Sort POIs by Levelpath and Hydrosequence in upstream to downstream order
+  seg_POIs <-  filter(st_drop_geometry(nhd_Full), COMID %in% tmp_POIs$COMID)
+  
+  # Attribute reaches that contribute locally to each POI
+  inc_segs <- segment_increment(nhd_Full, seg_POIs) 
+  nhd_Full <- nhd_Full %>%
+    left_join(select(inc_segs, COMID, POI_ID), by = "COMID") %>%
+    rename(TotDASqKM= denTotalAreaSqKM) %>%
+    mutate(StartFlag = ifelse(!Hydroseq %in% DnHydroseq, 1, 0),
+           minNet = 1, VPUID = "AK")
+  
+  # create and write out final dissolved segments
+  nsegments_fin <- segment_creation(nhd_Full)
+  nhd_Full <- select(nhd_Full, -POI_ID) %>% 
+    left_join(select(st_drop_geometry(nsegments_fin$raw_segs), COMID, POI_ID), by = "COMID")
+  nsegments <- nsegments_fin$diss_segs
+  
+  # Produce the minimal POIs needed to derive the network based on LP and terminal outlets
+  strucFeat <- structPOIsNet(nhd_Full, tmp_POIs)
+  nhd_Full <- nhd_Full %>% 
+    mutate(struct_POI = ifelse(COMID %in% strucFeat$struc_POIs$COMID, 1, 0),
+           struct_Net = ifelse(COMID %in% strucFeat$structnet$COMID, 1, 0))
+  
+  write_sf(nhd_Full, nav_gpkg, flowline)
+  write_sf(nsegments, nav_gpkg, n_segments)
+} else {
+  # Read in NHDPlusV2 flowline simple features and filter by vector processing unit (VPU)
+  tmp_POIs <- read_sf(nav_gpkg, pois_all)
+  nhd_Full <- read_sf(nav_gpkg, flowline)
+  nsegments <- read_sf(nav_gpkg, n_segments)
+}
+```
+
+```{r POI Collapse}
+#  Load data
+if(needs_layer(pois_merge)) {
+
+  #1 Move POIs downstream by category
+  out_gages <- POI_move_down(nav_gpkg, tmp_POIs, nsegments, filter(nhd_Final, !is.na(POI_ID)), "Type_Gages", .05)
+  out_HUC12 <- POI_move_down(nav_gpkg, tmp_POIs, nsegments, filter(nhd_Full, !is.na(POI_ID)), "Type_HUC12", .10)
+
+  nhd_Full <- select(nhd_Full, -POI_ID) %>%
+    left_join(st_drop_geometry(out_HUC12$FL) %>%
+                select(COMID, POI_ID), by = "COMID")
+
+  # Write out geopackage layer representing POIs for given theme
+  write_sf(out_HUC12$allPOIs, nav_gpkg, pois_merge)
+  write_sf(out_HUC12$FL, nav_gpkg, nhdflowline)
+  write_sf(out_HUC12$segs, nav_gpkg, n_segments)
+} else {
+  final_POIs <- read_sf(col_gpkg, pois_merge)
+  nhd_Final <- read_sf(nav_gpkg, nhdflowline)
+  nsegments <- read_sf(nav_gpkg, n_segments)
+}
+```