diff --git a/DESCRIPTION b/DESCRIPTION
index a0c93f19..14956775 100644
--- a/DESCRIPTION
+++ b/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"), 
diff --git a/NAMESPACE b/NAMESPACE
index 44924ca8..537594de 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/arcrest_tools.R b/R/arcrest_tools.R
index 47a6f6bd..407a031a 100644
--- a/R/arcrest_tools.R
+++ b/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 {
 
diff --git a/R/get_hydro.R b/R/get_hydro.R
index 5e63a71f..4d9a65d4 100644
--- a/R/get_hydro.R
+++ b/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)
+
+}
diff --git a/man/get_3dhp.Rd b/man/get_3dhp.Rd
index 2653e0fd..95035b63 100644
--- a/man/get_3dhp.Rd
+++ b/man/get_3dhp.Rd
@@ -21,9 +21,9 @@ in any Spatial Reference System.}
 
 \item{ids}{character vector of id3dhp ids or mainstem uris}
 
-\item{type}{character. Type of feature to return
+\item{type}{character. Type of feature to return. e.g.
 ("hydrolocation", "flowline", "waterbody", "drainage area", "catchment").
-If NULL (default) a data.frame of available resources is returned}
+If NULL (default) a data.frame of available types is returned}
 
 \item{universalreferenceid}{character vector of hydrolocation universal
 reference ids such as reachcodes}
diff --git a/man/get_nhdphr.Rd b/man/get_nhdphr.Rd
new file mode 100644
index 00000000..7b450f93
--- /dev/null
+++ b/man/get_nhdphr.Rd
@@ -0,0 +1,80 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/get_hydro.R
+\name{get_nhdphr}
+\alias{get_nhdphr}
+\title{Get NHDPlusHR Data}
+\usage{
+get_nhdphr(
+  AOI = NULL,
+  ids = NULL,
+  type = NULL,
+  reachcode = NULL,
+  t_srs = NULL,
+  buffer = 0.5,
+  page_size = 2000
+)
+}
+\arguments{
+\item{AOI}{sf (MULTI)POINT or (MULTI)POLYGON. An 'area of interest' can
+be provided as either a location (sf POINT) or area (sf POLYGON)
+in any Spatial Reference System.}
+
+\item{ids}{character vector of nhdplusid ids}
+
+\item{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}
+
+\item{reachcode}{character vector of reachcodes
+NOTE: performance of this query is currently very poor,
+spatial queries are the primary use of this function.}
+
+\item{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.}
+
+\item{buffer}{numeric. The amount (in meters) to buffer a POINT AOI by for an
+extended search. Default = 0.5}
+
+\item{page_size}{numeric default number of features to request at a time. Reducing
+may help if 500 errors are experienced.}
+}
+\value{
+a simple features (sf) object or valid types if no type supplied
+}
+\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.
+}
+\details{
+The returned object(s) will have the same
+Spatial Reference System (SRS) as the input AOI. If a individual or set of
+IDs are used to query, then the default CRS of EPSG:4269 is
+preserved. In all cases, a user-defined SRS can be passed to \code{t_srs}
+which will override all previous SRS (either input or default).
+All buffer and distance operations are handled internally using in
+EPSG:5070 Albers Equal Area projection
+}
+\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")
+}
+}
diff --git a/man/query_usgs_arcrest.Rd b/man/query_usgs_arcrest.Rd
index 98d23004..a12046a7 100644
--- a/man/query_usgs_arcrest.Rd
+++ b/man/query_usgs_arcrest.Rd
@@ -8,6 +8,7 @@ query_usgs_arcrest(
   AOI = NULL,
   ids = NULL,
   type = NULL,
+  service = NULL,
   where = NULL,
   t_srs = NULL,
   buffer = 0.5,
@@ -23,10 +24,11 @@ in any Spatial Reference System.}
 type requested, for 3dhp, this is id3dhp.}
 
 \item{type}{character. Type of feature to return
-("hydrolocation", "flowline", "waterbody", "drainage area", "catchment").
-If NULL (default) a data.frame of available resources is returned}
+If NULL (default) a data.frame of available types is returned}
 
-\item{where}{character. An where clause to pass to the server.}
+\item{service}{character chosen from "3DHP_all", "NHDPlus_HR"}
+
+\item{where}{character An where clause to pass to the server.}
 
 \item{t_srs}{character (PROJ string or EPSG code) or numeric (EPSG code).
 A user specified - target -Spatial Reference System (SRS/CRS) for returned objects.
diff --git a/nhdplusTools.Rproj b/nhdplusTools.Rproj
index 29736fb3..35418cc2 100644
--- a/nhdplusTools.Rproj
+++ b/nhdplusTools.Rproj
@@ -1,4 +1,5 @@
 Version: 1.0
+ProjectId: 903deedb-5779-4eeb-a8bc-71f5168cf6f9
 
 RestoreWorkspace: Default
 SaveWorkspace: Default
diff --git a/tests/testthat/test_arcrest.R b/tests/testthat/test_arcrest.R
index a078950f..0fcb6b63 100644
--- a/tests/testthat/test_arcrest.R
+++ b/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")
 })
diff --git a/tests/testthat/test_get_nhdphr.R b/tests/testthat/test_get_nhdphr.R
new file mode 100644
index 00000000..33e31950
--- /dev/null
+++ b/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")
+
+})