Retrieve NHDPlus catchment FEATUREID (COMID) for a set of point locations by batching points to an ArcGIS Feature Service /query endpoint and mapping the returned catchment geometries back to the input points via a local spatial intersection. This avoids per-point HTTP requests and is faster for moderate-to-large point sets.
sc_get_comid(
points_sf,
service_layer_url = paste0("https://watersgeo.epa.gov/arcgis/rest/services/",
"NHDPlus_NP21/Catchments_NP21_Simplified/MapServer/0"),
layer_id = NULL,
chunk_size = 1000,
feature_id_field = "FEATUREID",
verbose = TRUE,
crs = NULL,
coords = NULL
)An `sf` POINT object containing the input locations. If you supply a plain `data.frame`, convert it to an `sf` object first.
URL to the ArcGIS service layer. Defaults to the EPA NHDPlus catchments endpoint.
Deprecated. This argument is kept only for backward compatibility. If supplied, it is appended to `service_layer_url` after trimming any trailing slash.
Number of points to send per batch request. Default is 1000.
Name of the feature attribute to return as the COMID. Default is `"FEATUREID"`.
Logical; print progress messages for each chunk. Default is `TRUE`.
Optional CRS for the input points when a `data.frame` is supplied or when an `sf` object has no attached CRS. Use an EPSG code such as 4326.
Optional character vector of length 2 giving x/y coordinate column names when `points_sf` is a data frame. Example: `c("lon", "lat")`.
An integer vector of FEATUREID (COMID) values in the same order as the input points. Missing values are returned as `NA`.
if (FALSE) { # \dontrun{
pts <- sf::st_as_sf(
data.frame(
id = 1:3,
x = c(-90, -90.1, -90.2),
y = c(40, 40.1, 40.2)
),
coords = c("x", "y"),
crs = 4326
)
comids <- sc_get_comid(points_sf = pts, chunk_size = 500)
} # }