I'm trying to download watershed boundary dataset geometries from the USGS National Map. Right now, I'm just trying to download 2 digit HU layers. You can see what I'm looking for here and here is a low res image of what geometries I'm after (white outlines).

Here's the code I have so far:


# Base URL for the WBD service
base_url <- "https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer/1/query?"

# Define query parameters
query_params <- list(
  where = "1=1",      # Query to select all features
  outFields = "*",    # Return all fields
  f = "json"          # Output format (e.g., geojson, json, etc.)

# Construct the full URL with query parameters
full_url <- modify_url(base_url, query = query_params)

# Make the GET request
response <- GET(full_url)

# Check for a successful response and parse the data
if (http_type(response) == "application/json") {
  # Parse the response
  data <- content(response, "parsed")
  # Process the data as needed
} else {
  cat("Unexpected response type:", http_type(response))

The query appears to work as I get a 200 response when I run response$status. But when I go to parse the data using the content function, I get a 500 response.

Anything I can do to make this work?

  • The 200 means the server returned ok, the 500 is the service saying you've got an error, in this case the JSON is {"error":{"code":500,"message":"Error performing query operation","details":[]}} - maybe you are asking for too much information? Can you do a more restrictive query rather than asking for all features and all fields? Does this work through a web browser?
    – Spacedman
    Nov 5, 2023 at 19:24
  • Thank you @Spacedman! I spent some more time working through this and I think you are correct. It appears to be a size related issue because I can call for one HU at a time and it works. I'll update this with what worked for me shortly. If only 500 gave a little more info for the error...
    – setty
    Nov 5, 2023 at 20:52

1 Answer 1


I spent a few more hours on this and I think it ended up being a size issue (thank you @Spacedman). Additionally, parsing the geoJSON was done best with the geojsonsf package as opposed to using the httr content function.

Here's the code that ended up working for me to download the 2 digit HU layers:


# load the ids for all objects
id_url <- "https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer/1/query?where=1%3D1&text=&objectIds=&time=&geometry=&geometryType=esriGeometryPolygon&inSR=&spatialRel=esriSpatialRelIntersects&distance=&units=esriSRUnit_Foot&relationParam=&outFields=&returnGeometry=true&returnTrueCurves=false&maxAllowableOffset=&geometryPrecision=&outSR=&havingClause=&returnIdsOnly=true&returnCountOnly=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&returnZ=false&returnM=false&gdbVersion=&historicMoment=&returnDistinctValues=false&resultOffset=&resultRecordCount=&returnExtentOnly=false&datumTransformation=&parameterValues=&rangeValues=&quantizationParameters=&featureEncoding=esriDefault&f=json"

response <- GET(id_url)

# Check for a successful response and parse the data
if (http_type(response) == "application/json") {
    # Parse the response to text then use fromJSON to convert read json
    ids <- content(response, "text") %>%
        fromJSON() %>%
    # Process the data as needed
} else {
    cat("Unexpected response type:", http_type(response))

# define the base url for the layer queries
base_url <- "https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer/1/query?"

# Define query parameters using map to iterate over the different object id's
# obtained via the query above
query_params <- purrr::map(ids$objectIds, ~ list(
    # objectIds filters to one HU at a time from the map call
    objectIds = .,      
    # return polygons
    geometryType = "esriGeometryPolygon",
    # return geometries
    returnGeometry = "true",
    # Output format (e.g., geojson, json, etc.)
    f = "geojson" 

# create list of full urls to map over
full_urls <- purrr::map(query_params, ~ modify_url(base_url, query = .))

# get a subset of the full urls to load for smaller sample
# full dataset appears to be 55.6 MB in size
full_urls_small <- full_urls[1:5]

# map over the full urls with geojson_sf
wbd_data <- purrr::map(full_urls_small, st_read)

# reduce to a dataframe where each row is an sf object
wbd_data_df <- wbd_data %>%
    purrr::reduce(bind_rows) %>%
    # texas has some weird lines going on inside the shape, but this does not
    # remove them like I thought it would

# plot
leaflet() %>%
    addTiles() %>%
    addPolygons(data=wbd_data_df, weight=2, color="blue")

I think my initial issue was not fully understanding how to build the query. Using the query tool found at the bottom of this page was very helpful for debugging. This article here pointed me in the right direction to make sure I was actually generating a correct query.

  • You need to include purrr or something that provides the map function to make this work, but I think that's making heavy work of sticking the id numbers into a URL string (maybe later work will be changing the template URL in more complex ways...idk) and I'd just use paste0. Also you can use sf::st_read on a URL and it will get the geoJSON - its possible what geojson_sf is doing anyway, but takes another potential dependency (failure point) out of your code.
    – Spacedman
    Nov 6, 2023 at 12:43
  • Thanks for catching the missing purrr. I had it loaded in my environment so I didn't notice that. As for the heavy work, I will at some point be making more complicated calls, but I agree, if I can rework this it would be better. It seems overdone to map and apply a new id to each url. And you are correct on geojson_sf. I switched it to just be st_read and removed that dependency. Thank you!
    – setty
    Nov 6, 2023 at 20:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.