Convert X,Y coordinates to hydrologic unit boundaries at the specified level

point_to_huc.R

#' Convert X,Y coordinates to hydrologic unit boundaries at the specified level.
#'
#' Hydrologic unit boundaries are retrieved from the USGS NationalMap ArcGIS MapServer web services: `"https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer/"`.
#'
#' @param x A `SpatVector` object or `numeric` vector.
#' @param y `NULL` (when x is a spatial object); otherwise a `numeric` vector equal in length to `x`
#' @param layerid numeric. `1:8` where `2*layerid` is equal to the number of digits in the Hydrologic Unit Code (HUC). Default: `5` for "HUC-10" watershed. See details.
#' @param sr_in integer. Spatial Reference System of input (`x`, `y`) as specified with numeric code. Default: `4326` for `"EPSG:4326"`.
#' @param sr_out integer. Spatial Reference System of result as specified with numeric code. Default: `sr_in`; equivalent to input.
#' @param base_url Default: `"https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer/"`
#' @details
#' The levels of `layerid` correspond to the numeric codes for the following hydrologic units (HU):
#' - WBDLine (0)
#' - 2-digit HU (Region) (1)
#' - 4-digit HU (Subregion) (2)
#' - 6-digit HU (Basin) (3)
#' - 8-digit HU (Subbasin) (4)
#' - 10-digit HU (Watershed) (5)
#' - 12-digit HU (Subwatershed) (6)
#' - 14-digit HU (7)
#' - 16-digit HU (8)
#' @return A `SpatVector` object derived from GeoJSON.
#' @export
#' @examplesIf !inherits(requireNamespace("terra", quietly = TRUE), 'try-error')
#' @examples
#' \dontrun{
#' point_to_huc(-120, 36:39)
#' }
point_to_huc <- function(x, y = NULL, layerid = 5, sr_in = 4326, sr_out = sr_in,
                         base_url = "https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer/") {
  
  stopifnot(requireNamespace("terra"))
  
  # convert to SpatVector points; polygons -> centroids  
  if (!is.numeric(x)) {
    if (!inherits(x, 'SpatVector')) {
      x <- terra::vect(x)
    }
    if (!terra::is.points(x)) {
      x <- terra::centroids(x)
    }
    xx <- terra::crds(x)
    x <- xx[,1]
    y <- xx[,2]
  }
  
  stopifnot(all(layerid %in% 0:8))
  
  if (length(y) < length(x)) {
    y <- rep(y, length(x))
  } else if (length(x) < length(y)) {
    x <- rep(x, length(y))
  }
  
  if (length(x) != length(y)) {
    stop("`x` should have length 1 or length equal to `y`", call. = FALSE)
  }
  
  if (length(layerid) < length(x)) {
    layerid <- rep(layerid, length(x))
  }
  
  if (length(layerid) != length(x)) {
    stop("`layerid` should have length 1 or length equal to `x`", call. = FALSE)
  }
  
  res <- lapply(seq_along(x), function(i) {
    in_geom <- utils::URLencode(sprintf("%%28%s,%s%%29", x[i], y[i]))
    queryurl <- paste0(base_url, layerid[i], "/query")
    urx <- paste0(queryurl, "?geometry=", in_geom, 
                  "&geometryType=esriGeometryPoint&inSR=", sr_in,
                  "&spatialRel=esriSpatialRelIntersects&outSR=", sr_out,
                  "&f=geojson")
    terra::vect(urx)
  })
  
  do.call('rbind', res)
}