initial implementation of nhdplushr service function. fixes #415

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: nhdplusTools
 Type: Package
 Title: NHDPlus Tools
-Version: 1.3.0
+Version: 1.3.1
 Authors@R: c(person(given = "David", 
                     family = "Blodgett", 
                     role = c("aut", "cre"), 

NAMESPACE

@@ -32,6 +32,7 @@ export(get_huc)
 export(get_hydro_location)
 export(get_levelpaths)
 export(get_nhdarea)
+export(get_nhdphr)
 export(get_nhdplus)
 export(get_nhdplushr)
 export(get_nldi_basin)

R/arcrest_tools.R

@@ -1,15 +1,17 @@
-get_3dhp_service_info <- memoise::memoise(function() {
+get_arcrest_service_info <- memoise::memoise(function(service = "3DHP_all") {
 
-  # TODO: support more services?
-  server <- "3DHP_all"
+  stopifnot(service %in% c("3DHP_all", "NHDPlus_HR"))
 
   url_base <- paste0(get("arcrest_root", envir = nhdplusTools_env),
-                     server,
+                     service,
                      "/MapServer/")
 
   all_layers <- jsonlite::read_json(paste0(url_base, "?f=json"))
 
-  list(url_base = url_base, all_layers = all_layers)
+  id_name <- "id3dhp"
+  if(service == "NHDPlus_HR") id_name <- "nhdplusid"
+
+  list(url_base = url_base, all_layers = all_layers, id = id_name)
 
 })
 
@@ -31,9 +33,9 @@ get_3dhp_service_info <- memoise::memoise(function() {
 #' @param ids character. A set of identifier(s) from the data
 #' type requested, for 3dhp, this is id3dhp.
 #' @param type character. Type of feature to return
-#' ("hydrolocation", "flowline", "waterbody", "drainage area", "catchment").
-#' If NULL (default) a data.frame of available resources is returned
-#' @param where character. An where clause to pass to the server.
+#' If NULL (default) a data.frame of available types is returned
+#' @param service character chosen from "3DHP_all", "NHDPlus_HR"
+#' @param where character An where clause to pass to the server.
 #' @param t_srs  character (PROJ string or EPSG code) or numeric (EPSG code).
 #' A user specified - target -Spatial Reference System (SRS/CRS) for returned objects.
 #' Will default to the CRS of the input AOI if provided, and to 4326 for ID requests.
@@ -48,12 +50,13 @@ get_3dhp_service_info <- memoise::memoise(function() {
 #' @importFrom dplyr filter
 #' @importFrom methods as
 query_usgs_arcrest <- function(AOI = NULL,  ids = NULL,
-                               type = NULL, where = NULL,
+                               type = NULL, service = NULL,
+                               where = NULL,
                                t_srs = NULL,
                                buffer = 0.5,
-                               page_size = 2000){
+                               page_size = 2000) {
 
-  si <- get_3dhp_service_info()
+  si <- get_arcrest_service_info(service)
 
   source <- data.frame(user_call = sapply(si$all_layers$layers, \(x) tolower(x$name)),
                        layer = sapply(si$all_layers$layers, \(x) x$id))
@@ -73,7 +76,8 @@ query_usgs_arcrest <- function(AOI = NULL,  ids = NULL,
     return(NULL)
   }
 
-  if(grepl(paste(sapply(group_layers, \(x) x$name),
+  if(length(group_layers) > 0 &&
+     grepl(paste(sapply(group_layers, \(x) x$name),
                  collapse = "|"),
            type, ignore.case = TRUE)) {
     layer_id <- filter(source, .data$user_call == !!type)$layer
@@ -107,8 +111,14 @@ query_usgs_arcrest <- function(AOI = NULL,  ids = NULL,
 
       if(!is.null(where)) stop("can't specify ids and where")
 
-      where <- paste0("id3dhp IN ('",
-                      paste(ids, collapse = "', '"), "')")
+      if(si$id == "nhdplusid") {
+        where <- paste0(si$id, " IN (",
+                        paste(ids, collapse = ", "), ")")
+      } else {
+        where <- paste0(si$id, " IN ('",
+                        paste(ids, collapse = "', '"), "')")
+      }
+
     }
 
     post_body <- list(where = where,
@@ -186,8 +196,12 @@ query_usgs_arcrest <- function(AOI = NULL,  ids = NULL,
       if(inherits(out[[1]], "data.frame")) {
         out <- bind_rows(unify_types(out))
 
-        out <- check_valid(out[!duplicated(out[["id3dhp"]]), ],
-                           out_prj = t_srs)
+        if("id3dhp" %in% names(out)) {
+          out <- check_valid(out[!duplicated(out[["id3dhp"]]), ],
+                             out_prj = t_srs)
+        } else {
+          out <- check_valid(out, out_prj = t_srs)
+        }
 
       } else {
 

R/get_hydro.R

@@ -185,6 +185,9 @@ get_nwis <- function(AOI = NULL, t_srs = NULL, buffer = 20000){
 #' for source data documentation.
 #'
 #' @inherit query_usgs_arcrest details return params
+#' @param type character. Type of feature to return. e.g.
+#' ("hydrolocation", "flowline", "waterbody", "drainage area", "catchment").
+#' If NULL (default) a data.frame of available types is returned
 #' @param ids character vector of id3dhp ids or mainstem uris
 #' @param universalreferenceid character vector of hydrolocation universal
 #' reference ids such as reachcodes
@@ -248,7 +251,67 @@ get_3dhp <- function(AOI = NULL, ids = NULL, type = NULL,
     ids <- NULL
   }
 
-  query_usgs_arcrest(AOI, ids, type, where, t_srs, buffer, page_size)
+  query_usgs_arcrest(AOI, ids, type, "3DHP_all", where, t_srs, buffer, page_size)
 
 }
 
+#' Get NHDPlusHR Data
+#' @description
+#' Calls the NHDPlus_HR web service and returns sf data.frames for the selected
+#' layers. See https://hydro.nationalmap.gov/arcgis/rest/services/NHDPlus_HR/MapServer
+#' for source data documentation.
+#'
+#' @inherit query_usgs_arcrest details return params
+#'
+#' @param type character. Type of feature to return e.g.
+#' c("networknhdflowline", nonnetworknhdflowline", nhdwaterbody", "nhdpluscatchment").
+#' If NULL (default) a data.frame of available types is returned
+#'
+#' @param ids character vector of nhdplusid ids
+#'
+#' @param reachcode character vector of reachcodes
+#' NOTE: performance of this query is currently very poor,
+#' spatial queries are the primary use of this function.
+#'
+#' @export
+#' @examples
+#' \donttest{
+#' AOI <- sf::st_as_sfc(sf::st_bbox(c(xmin = -89.56684, ymin = 42.99816,
+#'                                    xmax = -89.24681, ymax = 43.17192),
+#'                                  crs = "+proj=longlat +datum=WGS84 +no_defs"))
+#'
+#' # get flowlines and hydrolocations
+#' flowlines <- get_nhdphr(AOI = AOI, type = "networknhdflowline")
+#' point <- get_nhdphr(AOI = AOI, type = "nhdpoint")
+#' waterbody <- get_nhdphr(AOI = AOI, type = "nhdwaterbody")
+#'
+#' if(!is.null(waterbody) & !is.null(flowlines) & !is.null(point)) {
+#' plot(sf::st_geometry(waterbody), col = "lightblue", border = "lightgrey")
+#' plot(sf::st_geometry(flowlines), col = "blue", add = TRUE)
+#' plot(sf::st_geometry(point), col = "grey", pch = "+", add = TRUE) }
+#'
+#' # given universalreferenceid (reachcodes), can query for them but only
+#' # for hydrolocations. This is useful for looking up mainstem ids.
+#'
+#' get_nhdphr(reachcode = "13020101021927", type = "networknhdflowline")
+#'}
+get_nhdphr <- function(AOI = NULL, ids = NULL, type = NULL,
+                       reachcode = NULL,
+                       t_srs = NULL, buffer = 0.5,
+                       page_size = 2000) {
+
+  if(!is.null(reachcode) && !isTRUE(grepl("nhdplusgage|nhdpoint|networknhdflowline|nonnetworknhdflowline|flowdirection|nhdwaterbody",
+                                          type))) {
+    stop("reachcode not defined for ", type)
+  }
+
+  where <- NULL
+  if(!is.null(reachcode)) {
+    where <- paste(paste0("reachcode IN ('",
+                          paste(reachcode, collapse = "', '"), "')"))
+    if(!is.null(ids)) stop("can not specify both reachcode and other ids")
+  }
+
+  query_usgs_arcrest(AOI, ids, type, "NHDPlus_HR", where, t_srs, buffer, page_size)
+
+}

nhdplusTools.Rproj

@@ -1,4 +1,5 @@
 Version: 1.0
+ProjectId: 903deedb-5779-4eeb-a8bc-71f5168cf6f9
 
 RestoreWorkspace: Default
 SaveWorkspace: Default

tests/testthat/test_arcrest.R

@@ -61,12 +61,14 @@ test_that("basic 3dhp service requests", {
                                      xmax = -89.4, ymax = 43.1),
                                    crs = "+proj=longlat +datum=WGS84 +no_defs"))
 
-  expect_message(expect_s3_class(nhdplusTools:::query_usgs_arcrest(AOI),
+  expect_message(expect_s3_class(nhdplusTools:::query_usgs_arcrest(AOI, service = "3DHP_all"),
                                  "data.frame"))
 
-  expect_warning(expect_warning(nhdplusTools:::query_usgs_arcrest(AOI, type = "hydrolocation")))
+  expect_warning(expect_warning(nhdplusTools:::query_usgs_arcrest(AOI, service = "3DHP_all",
+                                                                  type = "hydrolocation")))
 
-  test_data <- nhdplusTools:::query_usgs_arcrest(AOI, type = "reach code, external connection")
+  test_data <- nhdplusTools:::query_usgs_arcrest(AOI, , service = "3DHP_all",
+                                                 type = "reach code, external connection")
 
   expect_s3_class(test_data, "sf")
 })

tests/testthat/test_get_nhdphr.R

@@ -0,0 +1,23 @@
+test_that("get_nhdphr", {
+
+  skip_on_cran()
+
+  expect_error(get_nhdphr(reachcode = "01234", type = "test"),
+               "not defined")
+
+  expect_error(get_nhdphr(reachcode = "01234", type = "networknhdflowline",
+                        ids = "1"),
+               "can not specify both")
+
+  f <- get_nhdphr(ids = "50001000103671",
+                  type = "networknhdflowline")
+
+  expect_equal(f$nhdplusid, 50001000103671)
+
+  expect_s3_class(f, "sf")
+
+  skip("performance")
+  f2 <- get_nhdphr(reachcode = f$reachcode,
+                   type = "networknhdflowline")
+
+})