Merge branch 'main' into 'main'

Network Update

See merge request !167

workspace/00_enhd_nwm_network.Rmd

----
-title: "NHDPlusV2 Data Prep"
-output: html_document
----
-
-Project: GFv2.0  
-Script purpose: Get enhdplus and NWM topology incorporated.
-Date: 12/24/2019  
-Author: [Dave Blodgett](dblodgett@usgs.gov)
-Notes: This script updates the catchment topology of NHDPlusV2 with two data sources then recalculates NHDPlus attributes.
-
-Most functionality is encapsulated in a `run_plus_attributes` function.
-
-```{r}
-# sbtools::item_file_download("5dcd5f96e4b069579760aedb", names = "enhd_nhdplusatts.csv", 
-#                             dest_dir = "cache")
-
-# To create this file, run the following.
-
-library(dplyr)
-library(sf)
-
-source("R/utils.R")
-
-source("R/config.R")
-
-enhd <- data.table::fread(data_paths$e2nhd_network, sep = ",", integer64 = "character")
-
-enhd <- as.data.frame(enhd)
-
-nwm_nc <- RNetCDF::open.nc(data_paths$nwm_network)
-
-nwm <- data.frame(link = RNetCDF::var.get.nc(nwm_nc, "link"), 
-                  to = RNetCDF::var.get.nc(nwm_nc, "to"),
-                  Length = RNetCDF::var.get.nc(nwm_nc, "Length"))
-
-RNetCDF::close.nc(nwm_nc)
-  
-nhd_net <- nhdplusTools::get_vaa()
-
-names(nhd_net) <- tolower(names(nhd_net))
-names(enhd) <- tolower(names(enhd))
-
-enhd <- filter(enhd, comid %in% nhd_net$comid)
-
-enhd$divergence[enhd$comid == 11976913] <- 0 # not a divergence
-enhd$divergence[enhd$comid == 17309884] <- 0 # not a divergence
-enhd$divergence[enhd$comid == 7692005] <- 1 # this should be main path
-
-nwm_sub_nhd <- filter(nwm, link %in% nhd_net$comid) %>%
-  select(comid = link, tocomid = to)
-
-nwm <- nwm %>%
-  mutate(lengthkm = Length / 1000) %>%
-  select(comid = link, tocomid = to, lengthkm)
-
-```
-
-Now join enhd to NHDPlusV2 comids and join tonode to fromnode to create "tocomid" column.
-
-```{r}
-
-net_new <- nhd_net %>%
-  select(comid, fcode, nameID = levelpathi, areasqkm) %>%
-  left_join(select(enhd, comid, lengthkm, reachcode, 
-                   fromnode = cfromnode, tonode = ctonode, divergence, 
-                   frommeas, tomeas), by = "comid") %>%
-  left_join(select(., tocomid = comid, fromnode), 
-            by = c("tonode" = "fromnode"))
-
-div <- filter(net_new, divergence == 2)
-
-# find diverted node groups that do not contain a main path.
-all_div_nodes <- filter(net_new, fromnode %in% div$fromnode) %>%
-  group_by(fromnode) %>% 
-  filter(!any(divergence == 1))
-
-# for nodes found where there is no main path, the divergence flag
-# needs to be set to 0. This seems to be from changes applied in the 
-# enhd dataset where diversions leave a reservoir but have no flowline
-# representing the diversion.
-net_new$divergence[net_new$comid %in% all_div_nodes$comid] <- 0
-
-all_div_groups <- all_div_nodes %>%
-  filter(n() > 1)
-
-if(nrow(all_div_groups) > 0) stop()
-## If any divergence nodes don't have both a 1 and a 2, then the network is broken.
-
-if(!all(nhd_net$comid %in% net_new$comid)) stop()
-## If not all source comids in the new set, the network is broken.
-
-```
-
-Now work in tocomids from National Water Model topology. This avoids unwanted flowlines and replaces tocomids that the nwm topology has changed.
-
-```{r}
-
-net_new <- clean_net(nhd_net, net_new, nwm_sub_nhd)
-
-net_new <- select(net_new, comid, tocomid, fcode, nameID, lengthkm, reachcode, frommeas, tomeas, areasqkm)
-
-```
-
-Need to add a bunch of required attributes for NHDPlusTools.
-
-Attributes added are:  
-- Arbolate Sum
-- Level Path
-- Hydrosequence
-- Terminal Path
-- terminal flag
-- Path Length
-- Down Level Path
-- Down Hydrosequence
-
-```{r}
-net_new <- add_plus_network_attributes(net_new, 
-                                       override = 5, 
-                                       cores = 12)
-
-net_new <- dplyr::select(net_new, -nameID) %>%
-  dplyr::rename(arbolatesu = weight)
-
-net_new$dnlevelpat[is.na(net_new$dnlevelpat)] <- 0
-
-coastal_comids <- nhd_net$comid[nhd_net$terminalfl == 1 & nhd_net$streamleve == 1]
-
-net_new$streamleve <- select(net_new, levelpathi, dnlevelpat, comid) %>%
-  mutate(coastal = comid %in% coastal_comids) %>%
-  get_streamlevel()
-
-net_new$streamorde <- nhdplusTools::get_streamorder(select(net_new, ID = comid, toID = tocomid), status = TRUE)
-
-net_new <- left_join(net_new, 
-                     select(nhd_net, comid, vpuin, vpuout, wbareatype, 
-                            slope, slopelenkm, ftype, gnis_name, 
-                            gnis_id, wbareacomi, hwnodesqkm, rpuid, 
-                            vpuid, roughness), by = "comid")
-
-data.table::fwrite(net_new, 
-                   file = "cache/enhd_nhdplusatts.csv", 
-                   sep = ",")
-
-net_new <- as.data.frame(data.table::fread("cache/enhd_nhdplusatts.csv", integer64 = "character"))
-
-# Verify that levelpath and hydrosequence coded right.
-lp_out <- net_new %>%
-  group_by(levelpathi) %>%
-  filter(hydroseq == min(hydroseq)) %>%
-  ungroup()
-
-lp_out_2 <- filter(net_new, hydroseq == levelpathi)
-
-if(nrow(lp_out) != nrow(lp_out_2)) stop("Something wrong with hydrosequence and levelpath relationships.")
-
-fline <- sf::read_sf(data_paths$nhdplus_gdb, "NHDFlowline_Network")
-
-fline <- select(fline, COMID)
-
-fline <- right_join(fline, net_new, by = c("COMID" = "comid"))
-
-write_sf(fline, "cache/enhd_atts_fline.gpkg")
-
-fst::write_fst(net_new, "cache/enhd_nhdplusatts.fst")
-arrow::write_parquet(net_new, "cache/enhd_nhdplusatts.parquet")
-
-# use infrequently!! 
-# mainstems <- sf::read_sf("https://www.hydroshare.org/resource/4a22e88e689949afa1cf71ae009eaf1b/data/contents/mainstems.gpkg")
-# 
-# mainstems <- select(sf::st_drop_geometry(mainstems), 
-#                     id, uri)
-# 
-# mainstems <- mutate(mainstems, id = as.integer(id)) %>%
-#   rename(mainstem_id = id)
-# 
-# mainstem_lookup <- select(net_new, nhdpv2_comid = comid, mainstem_id = levelpathi)
-# mainstem_lookup <- left_join(mainstem_lookup, mainstems, by = "mainstem_id")
-# readr::write_csv(mainstem_lookup, "data/mainstem_lookup.csv")
-# R.utils::gzip("data/mainstem_lookup.csv")
-
-### ONLY RUN IF CHANGED ###
-# sbtools::authenticate_sb()
-# sbtools::item_replace_files("60c92503d34e86b9389df1c9", "cache/enhd_nhdplusatts.csv")
-# sbtools::item_replace_files("60c92503d34e86b9389df1c9", "cache/enhd_nhdplusatts.fst")
-# sbtools::item_replace_files("60c92503d34e86b9389df1c9", "cache/enhd_nhdplusatts.parquet")
-
-```

workspace/00_get_data.Rmd

@@ -49,7 +49,7 @@ if(!file.exists(hu12_points_path)) {
   if(is.null(sbtools::current_session()))
     sb <- authenticate_sb()
   
-  sbtools::item_file_download("60cb5edfd34e86b938a373f4", names = "102020wbd_outlets.gpkg", 
+  sbtools::item_file_download("63cb38b2d34e06fef14f40ad", names = "102020wbd_outlets.gpkg", 
                               destinations = hu12_points_path, session = sb)
 }
 
@@ -470,32 +470,6 @@ if(!file.exists(file.path(islands_dir, "hi.gpkg"))) {
 out_list <- c(out_list, out_hi)
 ```
 
-```{r e2nhd}
-
-zip_file <- "e2nhdplusv2_us_csv.zip"
-out_csv <- "e2nhdplusv2_us.csv"
-
-out_zip <- file.path(data_dir, out_csv)
-
-out <- list(e2nhd_network = file.path(data_dir, out_csv))
-
-if(!file.exists(out$e2nhd_network)) {
-  
-  if(is.null(sbtools::current_session()))
-    authenticate_sb()
-
-  sbtools::item_file_download("5d16509ee4b0941bde5d8ffe",
-                              names = zip_file,
-                              destinations = file.path(data_dir, zip_file))
-
-  zip::unzip(file.path(data_dir, zip_file), exdir = data_dir)
-
-}
-
-out_list <- c(out_list, out)
-
-```
-
 ```{r nwm_topology}
 nwm_targz_url <- "https://www.nohrsc.noaa.gov/pub/staff/keicher/NWM_live/NWM_parameters/NWM_parameter_files_v2.1.tar.gz"
 nwm_parm_url <- "https://www.nohrsc.noaa.gov/pub/staff/keicher/NWM_live/web/data_tools/NWM_v2.1_channel_hydrofabric.tar.gz"
@@ -534,13 +508,13 @@ out_list <- c(out_list, out)
 
 ```{r nhdplus_attributes}
 
-out <- list(new_nhdp_atts = file.path("cache", (sb_f <- "enhd_nhdplusatts.csv")))
+out <- list(new_nhdp_atts = file.path("data", (sb_f <- "enhd_nhdplusatts.csv")))
 
 if(!file.exists(out$new_nhdp_atts)) {
   if(is.null(sbtools::current_session()))
     authenticate_sb()
 
-  sbtools::item_file_download("60c92503d34e86b9389df1c9",
+  sbtools::item_file_download("63cb311ed34e06fef14f40a3",
                               names = sb_f,
                               destinations = out$new_nhdp_atts)
 }
@@ -548,22 +522,6 @@ if(!file.exists(out$new_nhdp_atts)) {
 out_list <- c(out_list, out)
 ```
 
-```{r nwm_nhdplus_attributes}
-
-out <- list(new_nwm_nhdp_atts = file.path("cache", (sb_f <- "nwm_enhd_nhdplusatts.csv")))
-
-if(!file.exists(out$new_nwm_nhdp_atts)) {
-  if(is.null(sbtools::current_session()))
-    authenticate_sb()
-
-  sbtools::item_file_download("60c92503d34e86b9389df1c9",
-                              names = sb_f,
-                              destinations = out$new_nwm_nhdp_atts)
-}
-
-out_list <- c(out_list, out)
-```
-
 ```{r GFv1.1}
 GFv11_dir <- file.path(data_dir, "GFv11")
 out <- list(GFv11_gages_lyr = file.path(data_dir, "GFv11/GFv11_gages.rds"),
@@ -640,6 +598,7 @@ if(!file.exists(out$ref_flowline)) {
     staged_nhd <- stage_national_data(nhdplus_data = out_list$nhdplus_gdb,
                                       output_path  = out_list$nhdplus_dir)
     
+    future::plan(future::multisession(workers = 16))
     # This will generate 
     fix_headwaters(staged_nhd$flowline, 
                    out$ref_flowline,

workspace/R/utils.R

@@ -161,73 +161,6 @@ gage_document <- function(data, source, drop_reason){
   }
 }
 
-#' clean network
-#' @param net data.frame nhdplus network
-#' @param net_new data.frame new nhdplus network as from e2nhd
-#' @param nwm data.frame with national water model new network
-#' @return data.frame fully reconciled attributes
-clean_net <- function(net, net_new, nwm) {
-  
-  # NOTE: need to disconnected diverted paths as indicated in the tonode / fromnode attributes.
-  
-  diverted_heads <- filter(net_new, divergence == 2)
-  
-  net_new <- filter(net_new, !tocomid %in% diverted_heads$comid)
-  
-  # NOTE: avoid fcode == 56600 (coastal)
-  avoid <- net$comid[net$fcode == 56600]
-  
-  # NOTE: avoid large groups where things have been redirected in an un-natural way.
-  groups <- group_by(nwm, tocomid)
-  
-  groups <- data.frame(tocomid = summarise(groups, .groups = "drop"), size = group_size(groups))
-  
-  groups <- filter(groups, size > 10 & tocomid != 0)
-  
-  # NOTE: also avoid all flowlines that go to coastal flowlines.
-  # This is due to unwanted great lakes connections in the NWM.
-  avoid <- c(avoid, 
-             net_new$comid[net_new$tocomid %in% avoid], 
-             net_new$comid[net_new$tocomid %in% groups$tocomid],
-             nwm$comid[nwm$tocomid %in% groups$tocomid])
-  
-  avoid <- unique(avoid)
-  
-  nwm <- filter(nwm, tocomid %in% c(net_new$comid, 0))
-  
-  net_new <- left_join(net_new, select(nwm, comid, 
-                                       nwm_tocomid = tocomid), 
-                       by = "comid") %>% 
-    mutate(tocomid = ifelse(is.na(tocomid), 0, tocomid))
-  ## NWM already uses 0 for outlets.
-  
-  # NOTE: just being explicit about this for clarity.
-  candidate <- select(net_new, comid, tocomid, nwm_tocomid) %>%
-    filter(!comid %in% avoid)
-  
-  mismatch <- candidate[candidate$tocomid != candidate$nwm_tocomid, ]
-  
-  net_new <- left_join(net_new, select(mismatch, comid, 
-                                       new_tocomid = nwm_tocomid), 
-                       by = "comid")
-  
-  to_change <- filter(net_new, !is.na(new_tocomid))
-  
-  net_new <- net_new %>%
-    mutate(divergence = ifelse(comid %in% to_change$new_tocomid & divergence == 2, 1,
-                               ifelse(comid %in% to_change$tocomid & divergence == 1, 2, 
-                                      divergence)),
-           tocomid = ifelse(comid %in% to_change$comid, new_tocomid, tocomid)) %>%
-    select(-nwm_tocomid, -new_tocomid)
-  
-  net_new <- filter(net_new, fcode != 56600)
-  
-  net_new <- mutate(net_new, tocomid = ifelse(!tocomid %in% comid, 0, tocomid))
-  
-  return(net_new)
-  
-}
-
 #' Merges geopackages together to create CONUs geopackage of features
 #' @param feat (character) Type of feature to merge (POIs, segments) character
 #' @param in_gpkg  (character) geopackage to merge (GF, ND_, etc.)

workspace/cache/data_paths.json

@@ -247,11 +247,9 @@
   "merit_fac": "data/merged_AK_MERIT_Hydro/ak_merit_fac.tif",
   "ak_source": "data/AK/ak.gpkg",
   "hi_source": "data/islands/hi.gpkg",
-  "e2nhd_network": "data/e2nhdplusv2_us.csv",
   "nwm_network": "data/NWM_parameters_v2.1/RouteLink_CONUS.nc",
   "nwm_parm": "data/NWM_v2.1_channel_hydrofabric_10262020/nwm_v2_1_hydrofabric.gdb",
-  "new_nhdp_atts": "cache/enhd_nhdplusatts.csv",
-  "new_nwm_nhdp_atts": "cache/nwm_enhd_nhdplusatts.csv",
+  "new_nhdp_atts": "data/enhd_nhdplusatts.csv",
   "GFv11_gages_lyr": "data/GFv11/GFv11_gages.rds",
   "GFv11_gdb": "data/GFv11/GFv1.1.gdb",
   "GFv11_tgf": "data/GFv11/TGF.gdb",