Merge pull request #116 from dblodgett-usgs/master

More mainstems functionality.

.Rbuildignore

@@ -13,3 +13,5 @@ code.json
 ^doc$
 ^Meta$
 .Renviron
+tests/testthat/data/03
+rosm.cache

.gitignore

@@ -10,3 +10,4 @@ inst/doc
 doc
 Meta
 .drake*
+tests/testthat/data/03

NAMESPACE

@@ -15,12 +15,14 @@ export(get_DM)
 export(get_UM)
 export(get_UT)
 export(get_flowline_index)
+export(get_hr_data)
 export(get_levelpaths)
 export(get_nhdplushr)
 export(get_nldi_basin)
 export(get_nldi_feature)
 export(get_pfaf)
 export(get_streamorder)
+export(make_standalone)
 export(navigate_nldi)
 export(nhdplus_path)
 export(plot_nhdplus)
@@ -49,11 +51,18 @@ importFrom(igraph,topo_sort)
 importFrom(jsonlite,fromJSON)
 importFrom(magrittr,"%>%")
 importFrom(methods,is)
+importFrom(sf,"st_geometry<-")
 importFrom(sf,read_sf)
+importFrom(sf,st_cast)
+importFrom(sf,st_crs)
 importFrom(sf,st_drop_geometry)
+importFrom(sf,st_geometry)
 importFrom(sf,st_layers)
 importFrom(sf,st_multilinestring)
 importFrom(sf,st_sf)
+importFrom(sf,st_simplify)
+importFrom(sf,st_transform)
+importFrom(sf,st_zm)
 importFrom(sf,write_sf)
 importFrom(tidyr,pivot_wider)
 importFrom(utils,download.file)

R/calc_network.R

@@ -286,14 +286,15 @@ get_streamorder <- function(fl) {
 #' @param max_level integer number of pfaf levels to attempt to calculate.
 #' If the network doesn't have resolution to support the desired level,
 #' unexpected behavior may occur.
+#' @param status boolean print status or not
 #' @return data.frame with ID and pfaf columns.
 #' @export
 #' @importFrom sf st_drop_geometry
 #' @importFrom methods is
 #' @examples
 #' library(dplyr)
 #' source(system.file("extdata/nhdplushr_data.R", package = "nhdplusTools"))
-#' hr_flowline <- nhdplusTools:::rename_nhdplus(hr_flowline)
+#' hr_flowline <- nhdplusTools:::rename_nhdplus(hr_data$NHDFlowline)
 #'
 #' fl <- prepare_nhdplus(hr_flowline, 0, 0, purge_non_dendritic = FALSE, warn = FALSE) %>%
 #'   left_join(select(hr_flowline, COMID, AreaSqKM), by = "COMID") %>%
@@ -311,9 +312,9 @@ get_streamorder <- function(fl) {
 #'
 #' plot(fl["pf_level_3"], lwd = 2)
 #'
-#' pfaf <- get_pfaf(fl, max_level = 7)
+#' pfaf <- get_pfaf(fl, max_level = 4)
 #'
-#' hr_catchment <- left_join(hr_catchment, pfaf, by = c("NHDPlusID" = "ID"))
+#' hr_catchment <- left_join(hr_data$NHDPlusCatchment, pfaf, by = c("FEATUREID" = "ID"))
 #'
 #' colors <- data.frame(pf_level_4 = unique(hr_catchment$pf_level_4),
 #'                      color = sample(terrain.colors(length(unique(hr_catchment$pf_level_4)))),
@@ -340,26 +341,29 @@ get_streamorder <- function(fl) {
 #'
 #' plot(fl["pf_level_2"], lwd = 2)
 #'
-get_pfaf <- function(fl, max_level = 2) {
+get_pfaf <- function(fl, max_level = 2, status = FALSE) {
   if(is(fl, "sf")) fl <- st_drop_geometry(fl)
   check_names(fl, "get_pfaf")
 
   mainstem_levelpath <- unique(fl$levelpath[fl$topo_sort == min(fl$topo_sort)])
 
   mainstem <- fl[fl$levelpath == mainstem_levelpath, ]
 
-  pfaf <- do.call(rbind, get_pfaf_9(fl, mainstem, max_level))
+  pfaf <- do.call(rbind, get_pfaf_9(fl, mainstem, max_level, status = status))
 
   return(cleanup_pfaf(pfaf))
 }
 
 #' @noRd
 #' @importFrom dplyr arrange left_join
 #' @importFrom methods is
-get_pfaf_9 <- function(fl, mainstem, max_level, pre_pfaf = 0, assigned = NA) {
+get_pfaf_9 <- function(fl, mainstem, max_level, pre_pfaf = 0, assigned = NA, status = FALSE) {
 
   if((pre_pfaf / 10^(max_level-1)) > 1) return()
 
+  if(status && ((pre_pfaf - 1111) %% 1000) == 0) {
+    message(paste("On level:", pre_pfaf - 1111))
+  }
   # Get all tributary outlets that go to the passed mainstem.
   trib_outlets <- fl[fl$toID %in% mainstem$ID &
                    fl$levelpath != mainstem$levelpath[1], ]
@@ -413,18 +417,18 @@ get_pfaf_9 <- function(fl, mainstem, max_level, pre_pfaf = 0, assigned = NA) {
   }
 
   c(out, unlist(lapply(c(1:9), apply_fun,
-                       p9 = out[[1]], fl = fl, max_level = max_level),
+                       p9 = out[[1]], fl = fl, max_level = max_level, status = status),
                 recursive = FALSE))
 }
 
-apply_fun <- function(p, p9, fl, max_level) {
+apply_fun <- function(p, p9, fl, max_level, status) {
   p_sub <- p9[p9$p_id == p, ]
   ms_ids <- p_sub$members
   pre_pfaf <- unique(p_sub$pfaf)
   mainstem <- fl[fl$ID %in% ms_ids, ]
 
   if(length(pre_pfaf) > 0) {
-    get_pfaf_9(fl, mainstem, max_level, pre_pfaf = pre_pfaf, assigned = p9)
+    get_pfaf_9(fl, mainstem, max_level, pre_pfaf = pre_pfaf, assigned = p9, status = status)
   } else {
     NULL
   }
@@ -443,8 +447,8 @@ cleanup_pfaf <- function(pfaf) {
   # Deduplicate problem tributaries
   remove <- do.call(c, lapply(1:length(unique(pfaf$level)), function(l, pfaf) {
     check <- pfaf[pfaf$level == l, ]
-    check <- dplyr::group_by(check, ID)
-    check <- dplyr::filter(check, n() > 1 & pfaf < max(pfaf))$uid
+    check <- dplyr::group_by(check, .data$ID)
+    check <- dplyr::filter(check, n() > 1 & .data$pfaf < max(.data$pfaf))$uid
   }, pfaf = pfaf))
 
   pfaf <- pivot_wider(select(pfaf[!pfaf$uid %in% remove, ], -.data$uid),

R/get_nhdplus.R

@@ -156,159 +156,6 @@ get_nhdplus_bybox <- function(box, layer) {
 
 }
 
-#' Download NHDPlus HiRes
-#' @param nhd_dir character directory to save output into
-#' @param hu_list character vector of hydrologic region(s) to download
-#' @param download_files boolean if FALSE, only URLs to files will be returned
-#' can be hu02s and/or hu04s
-#'
-#' @return Paths to geodatabases created.
-#' @importFrom xml2 read_xml xml_ns_strip xml_find_all xml_text
-#' @importFrom utils download.file unzip
-#' @export
-#' @examples
-#' \donttest{
-#' download_nhdplushr(tempdir(), c("01", "0203"), download_files = FALSE)
-#' }
-download_nhdplushr <- function(nhd_dir, hu_list, download_files = TRUE) {
-
-  nhdhr_bucket <- get("nhdhr_bucket", envir = nhdplusTools_env)
-  nhdhr_file_list <- get("nhdhr_file_list", envir = nhdplusTools_env)
-
-  hu02_list <- unique(substr(hu_list, 1, 2))
-  hu04_list <- hu_list[which(nchar(hu_list) == 4)]
-  subset_hu02 <- sapply(hu02_list, function(x)
-    sapply(x, function(y) any(grepl(y, hu04_list))))
-
-  out <- c()
-
-  for(h in 1:length(hu02_list)) {
-    hu02 <- hu02_list[h]
-
-    if(download_files) {
-      out <- c(out, file.path(nhd_dir, hu02))
-    }
-
-    if(download_files) {
-      dir.create(out[length(out)], recursive = TRUE, showWarnings = FALSE)
-    }
-
-    file_list <- read_xml(paste0(nhdhr_bucket, nhdhr_file_list,
-                                 "NHDPLUS_H_", hu02)) %>%
-      xml_ns_strip() %>%
-      xml_find_all(xpath = "//Key") %>%
-      xml_text()
-
-    file_list <- file_list[grepl("_GDB.zip", file_list)]
-
-    if(subset_hu02[h]) {
-      file_list <- file_list[sapply(file_list, function(f)
-        any(sapply(hu04_list, grepl, x = f)))]
-    }
-
-    for(key in file_list) {
-      out_file <- paste0(out[length(out)], "/", tail(strsplit(key, "/")[[1]], 1))
-      url <- paste0(nhdhr_bucket, key)
-
-      if(download_files & !dir.exists(gsub(".zip", ".gdb", out_file))) {
-        download.file(url, out_file)
-        unzip(out_file, exdir = out[length(out)])
-        unlink(out_file)
-      } else if(!download_files) {
-        out <- c(out, url)
-      }
-    }
-  }
-  return(out)
-}
-
-#' Get NHDPlus HiRes as single geopackage
-#' @param hr_dir character directory with geodatabases (gdb search is recursive)
-#' @param out_gpkg character path to write output geopackage
-#' @param layers character vector with desired layers to return.
-#' c("NHDFlowline", "NHDPlusCatchment") is default.
-#' Choose from:
-#' c("NHDFlowline", "NHDPlusCatchment", "NHDWaterbody", "NHDArea", "NHDLine",
-#' "NHDPlusSink", "NHDPlusWall", "NHDPoint", "NHDPlusBurnWaterbody",
-#' "NHDPlusBurnLineEvent", "HYDRO_NET_Junctions",
-#' "WBDHU2", "WBDHU4","WBDHU6", "WBDHU8" "WBDHU10", "WBDHU12", "WBDLine")
-#' @param pattern character optional regex to select certain files in hr_dir
-#'
-#' @details
-#' NHDFlowline is joined to value added attributes prior to being
-#' returned.
-#' Names are not modified from the NHDPlusHR geodatabase.
-#' Set layers to "NULL" to get all layers.
-#'
-#' @importFrom sf st_layers read_sf st_sf write_sf
-#' @export
-#' @examples
-#' \donttest{
-#' # Note this will download a lot of data to a temp directory.
-#' # Change 'tempdir()' to your directory of choice.
-#' download_dir <- download_nhdplushr(tempdir(), c("0302", "0303"))
-#' get_nhdplushr(download_dir, file.path(download_dir, "nhdplus_0302-03.gpkg"))
-#' }
-get_nhdplushr <- function(hr_dir, out_gpkg = NULL,
-                          layers = c("NHDFlowline", "NHDPlusCatchment"),
-                          pattern = ".*GDB.gdb$") {
-  gdb_files <- list.files(hr_dir, pattern = pattern,
-                          full.names = TRUE, recursive = TRUE, include.dirs = TRUE)
-
-  if(length(gdb_files) == 0) {
-    # For testing.
-    gdb_files <- list.files(hr_dir, pattern = "sub.gpkg", full.names = TRUE)
-  }
-  if(is.null(layers)) {
-    layers <- st_layers(gdb_files[1])
-
-    layers <- layers$name[!is.na(layers$geomtype) & layers$features > 0]
-  }
-
-  out_list <- list()
-
-  for(layer in layers) {
-    layer_set <- lapply(gdb_files, get_hr_data, layer = layer)
-
-    out <- do.call(rbind, layer_set)
-
-    try(out <- st_sf(out))
-
-    if(!is.null(out_gpkg)) {
-      write_sf(out, layer = layer, dsn = out_gpkg)
-    } else {
-      out_list[layer] <- list(out)
-    }
-  }
-  if(!is.null(out_gpkg)) {
-    return(out_gpkg)
-  } else {
-    return(out_list)
-  }
-}
-
-#' @noRd
-#' @importFrom sf st_multilinestring
-get_hr_data <- function(gdb, layer = NULL) {
-  if(layer == "NHDFlowline") {
-    hr_data <- suppressWarnings(read_sf(gdb, "NHDPlusFlowlineVAA"))
-    hr_data <- select(hr_data, -ReachCode, -VPUID)
-    hr_data <- left_join( sf::st_zm(read_sf(gdb, layer)), hr_data, by = "NHDPlusID")
-
-    fix <- which(!sapply(sf::st_geometry(hr_data), function(x) class(x)[2]) %in% c("LINESTRING", "MULTILINESTRING"))
-
-    for(f in fix) {
-      sf::st_geometry(hr_data)[[f]] <- st_multilinestring(lapply(sf::st_geometry(hr_data)[[f]][[1]],
-                                                                 sf::st_cast, to = "LINESTRING"), dim = "XY")
-    }
-
-    return(hr_data)
-
-  } else {
-    read_sf(gdb, layer)
-  }
-}
-
 make_web_sf <- function(content) {
   if(content$status_code == 200) {
     tryCatch(sf::read_sf(rawToChar(content$content)),