Merge pull request #419 from DOI-USGS/nhdplushr-client

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)

NEWS.md

@@ -1,3 +1,8 @@
+nhdplusTools 1.3.1
+==========
+
+- added a web service client for the NHDPlusHR argis server `get_nhdphr()`. #415
+
 nhdplusTools 1.3.0
 ==========
 

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)
+
+}

_pkgdown.yml

@@ -6,6 +6,7 @@ reference:
       Functions that help get access data. (all U.S. context)
   - contents:
       - '`get_3dhp`'
+      - '`get_nhdphr`'
       - '`discover_geoconnex_reference`'
       - '`get_geoconnex_reference`'
       - '`download_nhdplushr`'