Get 3DHP Data

Calls the 3DHP_all web service and returns sf data.frames for the selected layers. See https://3dhp.nationalmap.gov/arcgis/rest/services/usgs_3dhp_all/FeatureServer for source data documentation.

get_3dhp(
  AOI = NULL,
  ids = NULL,
  type = NULL,
  universalreferenceid = NULL,
  t_srs = NULL,
  buffer = 0.5,
  page_size = 2000
)

Arguments

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.

ids

character vector of id3dhp ids, mainstem uris, or workunitid prefixed ids (e.g. "workunitid:300585")

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

universalreferenceid

character vector of hydrolocation universal reference ids. The service only populates these on sink, spring and waterbody outlet hydrolocations.

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.

buffer

numeric. The amount (in meters) to buffer a POINT AOI by for an extended search. Default = 0.5

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

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 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_3dhp(AOI = AOI, type = "flowline")
#> GET https://3dhp.nationalmap.gov/arcgis/rest/services/usgs_3dhp_all/FeatureServer/?f=json
#> POST https://3dhp.nationalmap.gov/arcgis/rest/services/usgs_3dhp_all/FeatureServer/50/query
#> POST https://3dhp.nationalmap.gov/arcgis/rest/services/usgs_3dhp_all/FeatureServer/50/query
hydrolocation <- get_3dhp(AOI = AOI, type = "hydrolocation")
#> POST https://3dhp.nationalmap.gov/arcgis/rest/services/usgs_3dhp_all/FeatureServer/20/query
#> Warning: No hydrolocation - sink, spring, waterbody outlet features found in area of interest.
#> POST https://3dhp.nationalmap.gov/arcgis/rest/services/usgs_3dhp_all/FeatureServer/30/query
#> Warning: No hydrolocation - headwater, terminus, divergence, confluence, catchment outlet features found in area of interest.
#> POST https://3dhp.nationalmap.gov/arcgis/rest/services/usgs_3dhp_all/FeatureServer/40/query
#> Warning: No hydrolocation - external connection features found in area of interest.
waterbody <- get_3dhp(AOI = AOI, type = "waterbody")
#> POST https://3dhp.nationalmap.gov/arcgis/rest/services/usgs_3dhp_all/FeatureServer/60/query
#> POST https://3dhp.nationalmap.gov/arcgis/rest/services/usgs_3dhp_all/FeatureServer/60/query

if(!is.null(waterbody) & !is.null(flowlines) & !is.null(hydrolocation)) {
plot(sf::st_geometry(waterbody), col = "lightblue", border = "lightgrey")
plot(sf::st_geometry(flowlines), col = "blue", add = TRUE)
plot(sf::st_geometry(hydrolocation), col = "grey", pch = "+", add = TRUE) }

# given mainstem ids from any source, can query for them in ids.

SU <- get_3dhp(ids = "https://geoconnex.us/ref/mainstems/194408",
               type = "flowline")
#> POST https://3dhp.nationalmap.gov/arcgis/rest/services/usgs_3dhp_all/FeatureServer/50/query
#> POST https://3dhp.nationalmap.gov/arcgis/rest/services/usgs_3dhp_all/FeatureServer/50/query

if(!is.null(SU))
  plot(sf::st_geometry(SU), col = "blue")


# get all the waterbodies along the Susquehanna river
SU_wb <- get_3dhp(ids = unique(SU$waterbodyid3dhp), type = "waterbody")
#> POST https://3dhp.nationalmap.gov/arcgis/rest/services/usgs_3dhp_all/FeatureServer/60/query
#> POST https://3dhp.nationalmap.gov/arcgis/rest/services/usgs_3dhp_all/FeatureServer/60/query

if(!is.null(SU_wb)) {
plot(sf::st_geometry(SU_wb[grepl("Otsego", SU_wb$gnisidlabel),]),
     col = "blue", border = "NA") }


# given a workunitid, can query for all features in that work unit with
# ids = "workunitid:300585" -- not run here as a work unit is tens of
# thousands of features.

# given universalreferenceid, can query for them but only
# for hydrolocations. This is useful for looking up mainstem ids.

urid <- unique(stats::na.omit(hydrolocation$universalreferenceid))

if(length(urid) > 0) {
get_3dhp(universalreferenceid = head(urid, 5),
         type = "hydrolocation")
}
# }