Overview

The USGS National Hydrologic Geospatial Fabric (NHGF) exposes its vector geospatial layers through a pygeoapi server at https://api.water.usgs.gov/fabric/pygeoapi. The service provides the spatial backbone of the NHGF — stream networks, catchments, waterbodies, hydrologic unit boundaries, and reference gage locations — all queryable over standard HTTP.

It implements two complementary OGC API standards:

• OGC API – Features — query and retrieve vector geospatial features (points, lines, polygons) as GeoJSON. This is the primary interface for data access.
• OGC API – Tiles — access pre-rendered Mapbox Vector Tiles (MVT) for performant map display at multiple zoom levels.

If you are new to OGC API: it is a family of RESTful standards that replace older OGC services (WMS, WFS, WCS) with modern JSON/HTTP patterns. Every endpoint returns standard JSON by default and supports content negotiation via a ?f= parameter — more on that below .

This service replaces the legacy GeoServer WFS at https://api.water.usgs.gov/geoserver/wmadata/ows, which is being retired. See Migration from Legacy GeoServer for details.

Available Collections

The collections below are fetched directly from the live API at /collections .

library(httr2)
library(jsonlite)

resp <- request(paste0(FABRIC_BASE, "/collections?f=json")) |>
  req_perform()

collections_json <- resp_body_json(resp)

# Extract collection metadata into a data frame
collections <- do.call(rbind, lapply(collections_json$collections, function(col) {
  # Get keywords if available
  keywords <- if (!is.null(col$keywords)) paste(col$keywords, collapse = ", ") else ""

  data.frame(
    id = col$id %||% "",
    title = col$title %||% col$id,
    description = col$description %||% "",
    keywords = keywords,
    stringsAsFactors = FALSE
  )
}))

# Assign groups based on collection ID patterns
classify_group <- function(id) {
  if (grepl("^(nhd(flowline|area|waterbody)|catchmentsp)", id)) return("NHDPlusV2 Hydrography")
  if (grepl("^gagesii", id)) return("Gages-II Reference Gages")
  if (grepl("^nhdplusv2-huc", id)) return("NHDPlusV2 Hydrologic Units")
  if (grepl("^nhdplushr-huc", id)) return("NHDPlusHR Hydrologic Units")
  if (grepl("2020", id)) return("WBD 2020 Hydrologic Units")
  if (grepl("2025", id)) return("WBD 2025 Hydrologic Units")
  return("Other")
}

collections$group <- sapply(collections$id, classify_group)

# Define group order
group_order <- c(
  "NHDPlusV2 Hydrography",
  "Gages-II Reference Gages",
  "NHDPlusV2 Hydrologic Units",
  "NHDPlusHR Hydrologic Units",
  "WBD 2020 Hydrologic Units",
  "WBD 2025 Hydrologic Units",
  "Other"
)

collections$group <- factor(collections$group, levels = group_order)
collections <- collections[order(collections$group, collections$id), ]

# Build linked collection IDs
collections$Collection <- paste0(
  "[`", collections$id, "`](",
  FABRIC_BASE, "/collections/", collections$id, "?f=html)"
)

# Build all output as a single string to keep code-fold block intact
out <- character()

for (grp in levels(collections$group)) {
  subset_df <- collections[collections$group == grp, ]
  if (nrow(subset_df) == 0) next
  out <- c(out, paste0("\n### ", grp, "\n"))
  tbl <- knitr::kable(
    subset_df[, c("Collection", "title", "description")],
    col.names = c("Collection", "Title", "Description"),
    row.names = FALSE
  )
  out <- c(out, tbl, "")
}

out <- c(out, paste0("\n*Table generated: ", format(Sys.time(), "%Y-%m-%d"),
  " from `", FABRIC_BASE, "/collections`*\n"))

cat(paste(out, collapse = "\n"))

NHDPlusV2 Hydrography

Gages-II Reference Gages

NHDPlusV2 Hydrologic Units

NHDPlusHR Hydrologic Units

WBD 2020 Hydrologic Units

WBD 2025 Hydrologic Units

Table generated: 2026-04-03 from https://api.water.usgs.gov/fabric/pygeoapi/collections

Connecting to the NLDI: The comid attribute in the NHDPlusV2 hydrography collections is the same identifier used by the NLDI to index features to the hydrologic network. If the NLDI returns a comid, you can retrieve the corresponding geometry here — for example, /collections/nhdflowline_network/items?comid=13297246 or /collections/catchmentsp/items?comid=13297246. See Connecting Fabric Data to the NLDI for a walkthrough.

Querying Features

Response Headers and Pagination

Every API response includes HTTP headers that help you understand the result set and navigate pages. Here is a request for a single flowline by comid, showing the most useful headers:

resp <- request(paste0(
  FABRIC_BASE,
  "/collections/nhdflowline_network/items?comid=13297246&limit=1&f=json"
)) |>
  req_perform()

headers <- resp_headers(resp)

# Build all output as a single string to avoid fragmented code-fold blocks
out <- character()
out <- c(out, "Request URL:")
out <- c(out, paste0("  ", FABRIC_BASE,
    "/collections/nhdflowline_network/items?comid=13297246&limit=1&f=json"), "")

out <- c(out, "Selected response headers:")
useful <- c("content-type", "access-control-allow-origin",
            "content-encoding", "x-powered-by")
for (h in useful) {
  val <- headers[[h]]
  if (!is.null(val)) out <- c(out, sprintf("  %s: %s", h, val))
}

# Show pagination info from the JSON body
body <- resp_body_json(resp)
out <- c(out, "", "Pagination fields in JSON response:")
out <- c(out, sprintf("  numberMatched: %s", body$numberMatched %||% "not provided"))
out <- c(out, sprintf("  numberReturned: %s", body$numberReturned %||% "not provided"))

# Check for next link
links <- body$links
next_link <- Filter(function(l) identical(l$rel, "next"), links)
if (length(next_link) > 0) {
  out <- c(out, sprintf("  next: %s", next_link[[1]]$href))
} else {
  out <- c(out, "  next: (none — all results fit on this page)")
}

cat(paste(out, collapse = "\n"), "\n")

Request URL:
  https://api.water.usgs.gov/fabric/pygeoapi/collections/nhdflowline_network/items?comid=13297246&limit=1&f=json

Selected response headers:
  content-type: application/json
  content-encoding: gzip

Pagination fields in JSON response:
  numberMatched: 1
  numberReturned: 1
  next: https://api.water.usgs.gov:443/fabric/pygeoapi/collections/nhdflowline_network/items?offset=1&comid=13297246&limit=1 

The key pagination fields are in the JSON response body, not HTTP headers:

• numberMatched — total features matching your query
• numberReturned — features on this page
• links — follow the rel: "next" link to get the next page

Use the limit parameter to control page size (default varies by collection).

Listing and Retrieving Items

Request features from any collection as a paginated GeoJSON FeatureCollection:

GET /collections/{collectionId}/items

To retrieve a single feature by its server-assigned ID:

GET /collections/{collectionId}/items/{featureId}

Filtering with Queryables

Each collection exposes filterable properties at its queryables endpoint. For example, see the flowline network queryables for the full list of available filters.

Filter by appending query parameters:

/collections/nhdflowline_network/items?comid=13297246
/collections/nhdflowline_network/items?gnis_name=Connecticut%20River
/collections/nhdplusv2-huc08/items?huc8=01080205

Bounding Box Queries

Use the bbox parameter to retrieve features within a geographic extent (min longitude, min latitude, max longitude, max latitude):

# Retrieve flowlines in a bounding box around Hartford, CT
resp <- request(paste0(
  FABRIC_BASE,
  "/collections/nhdflowline_network/items",
  "?bbox=-72.75,41.7,-72.6,41.8&limit=5&f=json"
)) |>
  req_perform()

body <- resp_body_json(resp)
out <- character()
out <- c(out, sprintf("numberMatched: %s", body$numberMatched))
out <- c(out, sprintf("numberReturned: %s (limit=5)", body$numberReturned))

# Show the first few feature IDs and names
for (feat in body$features[1:min(3, length(body$features))]) {
  out <- c(out, sprintf("  comid: %s  gnis_name: %s",
    feat$properties$comid %||% "—",
    feat$properties$gnis_name %||% "(unnamed)"))
}
cat(paste(out, collapse = "\n"), "\n")

numberMatched: 56
numberReturned: 5 (limit=5)
  comid: 7700130  gnis_name: Burnham Brook
  comid: 7700132  gnis_name:  
  comid: 7700136  gnis_name: Burnham Brook 

Coordinate Reference Systems

The service supports CRS negotiation per OGC API – Features Part 2 . The default storage CRS is NAD83. Use the crs parameter to request output in a different coordinate reference system. Supported CRS options are listed in each collection’s metadata.

Vector Tiles

Collections also expose Mapbox Vector Tile (MVT) tilesets for use in web mapping libraries such as MapLibre GL JS, Mapbox GL JS, or OpenLayers.

Tile URL pattern:

/collections/{collectionId}/tiles/{tileMatrixSetId}/{z}/{y}/{x}

Tileset metadata (available matrix sets and zoom levels) is at:

/collections/{collectionId}/tiles

Use vector tiles when you need to display features on an interactive map without downloading full GeoJSON datasets. Note that tiles may not include all attributes available through the Features interface — use Features queries when you need the complete attribute set.

Client Libraries

Python — requests

For Python, query the Fabric OGC API directly with requests. This teaches the standard OGC API – Features pattern that works with any OGC-compliant server.

Query features as GeoJSON

import requests
import geopandas as gpd

base = "https://api.water.usgs.gov/fabric/pygeoapi"

# Get a flowline by comid
resp = requests.get(
    f"{base}/collections/nhdflowline_network/items",
    params={"comid": 13297246, "f": "json", "limit": 1}
)
data = resp.json()

# Convert to GeoDataFrame
gdf = gpd.GeoDataFrame.from_features(data["features"], crs="EPSG:4269")
print(f"Columns: {list(gdf.columns[:6])} ...")

Columns: ['geometry', 'objectid', 'comid', 'fdate', 'resolution', 'gnis_id'] ...

gdf[["comid", "gnis_name", "lengthkm"]].head()

      comid gnis_name  lengthkm
0  13297246               0.108

Bounding box query

# Get flowlines in a bounding box around Hartford, CT
resp = requests.get(
    f"{base}/collections/nhdflowline_network/items",
    params={
        "bbox": "-72.75,41.7,-72.6,41.8",
        "limit": 5,
        "f": "json",
    }
)
data = resp.json()
print(f"numberMatched: {data['numberMatched']}")

numberMatched: 56

print(f"numberReturned: {data['numberReturned']}")

numberReturned: 5

gdf = gpd.GeoDataFrame.from_features(data["features"], crs="EPSG:4269")
gdf[["comid", "gnis_name", "lengthkm"]].head()

     comid      gnis_name  lengthkm
0  7700130  Burnham Brook     2.688
1  7700132                    1.291
2  7700136  Burnham Brook     2.671
3  7700140                    2.575
4  7700144                    0.138

Paginating through results

# Follow the "next" link to paginate
resp = requests.get(
    f"{base}/collections/nhdflowline_network/items",
    params={"bbox": "-72.75,41.7,-72.6,41.8", "limit": 2, "f": "json"}
)
data = resp.json()

print(f"Page 1: {data['numberReturned']} of {data['numberMatched']} features")

Page 1: 2 of 56 features

# Find the "next" link
next_links = [l for l in data.get("links", []) if l.get("rel") == "next"]
if next_links:
    print(f"Next page URL: {next_links[0]['href'][:80]}...")
    # Follow it:
    resp2 = requests.get(next_links[0]["href"])
    data2 = resp2.json()
    print(f"Page 2: {data2['numberReturned']} features")
else:
    print("No next page — all results fit on one page")

Next page URL: https://api.water.usgs.gov:443/fabric/pygeoapi/collections/nhdflowline_network/i...
Page 2: 2 features

Discovering collections

# List all collections and their types
resp = requests.get(f"{base}/collections", params={"f": "json"})
collections = resp.json()["collections"]

print(f"Total collections: {len(collections)}")

Total collections: 31

for c in collections[:5]:
    print(f"  {c['id']}: {c.get('title', c['id'])}")

  catchmentsp: NHDPlusV2 Catchmentsp Features
  nhdarea: NHDPlusV2 NHD Area Features
  nhdflowline_network: NHDPlusV2 NHD Flowline Network
  nhdflowline_nonnetwork: NHDPlusV2 NHD Flowline Nonnetwork Features
  nhdwaterbody: NHDPlusV2 NHD Waterbody Features

print("  ...")

  ...

NLDI → Fabric (Python)

Use PyNHD for NLDI navigation, then query Fabric directly with requests:

from pynhd import NLDI

nldi = NLDI()

# Step 1: Navigate the NLDI — find upstream flowlines
flowlines = nldi.navigate_byid(
    fsource="nwissite",
    fid="USGS-05428500",
    navigation="upstreamTributaries",
    source="flowlines",
    distance=10,
)
comids = flowlines["nhdplus_comid"].tolist()
print(f"Found {len(comids)} upstream flowlines")

Found 10 upstream flowlines

print(f"First few comids: {comids[:5]}")

First few comids: ['13293750', '13294312', '13294308', '13294310', '13293634']

# Step 2: Query Fabric for a catchment by comid
comid = comids[0]
resp = requests.get(
    f"{base}/collections/catchmentsp/items",
    params={"comid": comid, "f": "json", "limit": 1}
)
data = resp.json()
print(f"\nCatchment for comid {comid}:")

Catchment for comid 13293750:

print(f"  numberMatched: {data.get('numberMatched')}")

  numberMatched: 2647455

if data.get("features"):
    props = data["features"][0]["properties"]
    print(f"  areasqkm: {props.get('areasqkm')}")

  areasqkm: 2.310929

R — nhdplusTools

nhdplusTools provides functions that query the Fabric pygeoapi service directly. Output is always an sf data frame.

Query by comid

library(nhdplusTools)

# Get a specific flowline by comid
flowline <- get_nhdplus(comid = 13297246)
cat("Columns:", paste(names(flowline)[1:8], collapse = ", "), "...\n\n")

Columns: objectid, comid, fdate, resolution, gnis_id, gnis_name, lengthkm, reachcode ...

print(sf::st_drop_geometry(flowline)[, 1:6])

# A tibble: 1 × 6
  objectid    comid fdate               resolution gnis_id gnis_name
     <int>    <int> <dttm>              <chr>      <chr>   <chr>    
1   968679 13297246 2009-04-30 23:00:00 Medium     " "     " "      

Query by bounding box

library(sf)

# Create a bounding box around Hartford, CT
bbox <- st_bbox(c(xmin = -72.75, ymin = 41.7,
                   xmax = -72.6,  ymax = 41.8),
                crs = 4326)
aoi <- st_as_sfc(bbox)

# Get flowlines in the AOI
flowlines <- get_nhdplus(AOI = aoi)
cat(sprintf("Retrieved %d flowlines in AOI\n\n", nrow(flowlines)))

Retrieved 56 flowlines in AOI

print(head(sf::st_drop_geometry(flowlines)[, c("comid", "gnis_name", "lengthkm")]))

# A tibble: 6 × 3
    comid gnis_name        lengthkm
    <int> <chr>               <dbl>
1 7700130 "Burnham Brook"     2.69 
2 7700132 " "                 1.29 
3 7700136 "Burnham Brook"     2.67 
4 7700140 " "                 2.58 
5 7700144 " "                 0.138
6 7700146 "Hockanum River"    0.222

Hydrologic units

# Get a HUC8 boundary by ID
huc8 <- get_huc(id = "01080205", type = "huc08")
cat(sprintf("HUC8 %s: %s\n  Area: %s sq km\n",
    huc8$huc8, huc8$name,
    round(as.numeric(sf::st_area(huc8)) / 1e6)))

HUC8 01080205: Outlet Connecticut River
  Area: 2805 sq km

Gages-II

# Get GAGES-II gages by ID
gages <- get_gagesII(id = "01011000")
sf::st_drop_geometry(gages)[, c("staid", "staname", "drain_sqkm", "huc02")]

# A tibble: 1 × 4
  staid    staname                             drain_sqkm huc02
  <chr>    <chr>                                    <dbl> <chr>
1 01011000 Allagash River near Allagash, Maine      3187. 01   

Connecting Fabric Data to the NLDI

The NLDI uses NHDPlusV2 comid values as its spatial backbone. Every feature in the NLDI is indexed to a catchment and its associated flowline via comid. This means you can move between the two services:

1. NLDI → Fabric: Use the NLDI to navigate the network (e.g., find all flowlines upstream of a gage), then query Fabric for the corresponding catchment or waterbody geometries using the comid values from the NLDI response. See the Python and R examples in the Client Libraries section above.

2. Fabric → NLDI: Start with a spatial query in Fabric (e.g., flowlines in a bounding box), then use the returned comid values to look up linked data sources in the NLDI (gages, water quality sites, etc.) via /linked-data/comid/{comid}.

# Get a flowline from Fabric by comid
flowline <- get_nhdplus(comid = 13297246)
cat("Flowline from Fabric:\n")

Flowline from Fabric:

print(sf::st_drop_geometry(flowline)[, c("comid", "gnis_name", "lengthkm", "totdasqkm")])

# A tibble: 1 × 4
     comid gnis_name lengthkm totdasqkm
     <int> <chr>        <dbl>     <dbl>
1 13297246 " "          0.108     1305.

# Use this comid to look up the NLDI feature
nldi_comid <- list(featureSource = "comid", featureID = "13297246")
feature <- get_nldi_feature(nldi_comid)
cat("\nNLDI feature for same comid:\n")

NLDI feature for same comid:

print(sf::st_drop_geometry(feature))

# A tibble: 1 × 3
  sourceName    identifier comid   
* <chr>         <chr>      <chr>   
1 NHDPlus comid 13297246   13297246

Using the Built-In Developer Tools

The HTML views that pygeoapi generates serve as browser-based API introspection tools for developers during integration and testing, not as a general-purpose user interface. They render the same data the JSON API returns, formatted for human reading.

Every endpoint supports a ?f= parameter for content negotiation:

• ?f=json — GeoJSON (default for programmatic access)
• ?f=jsonld — JSON-LD (linked data format)
• ?f=html — human-readable HTML view (auto-generated by pygeoapi)

Useful developer entry points:

• Swagger UI — /openapi?f=html — interactive endpoint testing; try queries and inspect request/response shapes
• HTML collection browser — append ?f=html to any endpoint (e.g., /collections/nhdflowline_network?f=html ) to browse metadata and available query options
• Queryables viewer — /collections/nhdflowline_network/queryables?f=html — see all filterable properties for a collection

Key Endpoints

All paths are relative to https://api.water.usgs.gov/fabric/pygeoapi.

Migration from Legacy GeoServer

The legacy GeoServer WFS at https://api.water.usgs.gov/geoserver/wmadata/ows served the same core data as roughly 15 WFS layers (NHDPlusV2 hydrography, gagesii, and WBD layers with date-stamped names like wbd02_20201006). That service is being retired.

Key differences in the Fabric pygeoapi replacement:

• RESTful JSON instead of XML-based WFS GetFeature requests
• Stable collection IDs without date stamps (e.g., nhdplusv2-huc08 instead of wbd08_20201006)
• Vector tile support for web mapping use cases
• Queryables-based filtering instead of CQL/OGC filter encoding
• NHDPlusHR hydrologic units added (not available in the legacy service)