Drainage Area Estimation"
In nhdplusTools: NHDPlus Tools

```{css, echo = FALSE} / Override html_vignette narrow body to allow wider figures / body, div.main-container, div.body, div.contents { max-width: 1600px !important; width: 95% !important; } img { max-width: 100% !important; height: auto !important; } .vignette-note { background: #f8f8f8; border-left: 3px solid #ccc; padding: 8px 12px; margin: 8px 0 16px 0; font-size: 0.9em; }

```r
library(nhdplusTools)

#pkgdown build uses this
build_local <- (Sys.getenv("BUILD_VIGNETTES") == "TRUE")

#cran package uses this
cran_build <- (Sys.getenv("BUILD_VIGNETTES_CRAN") == "TRUE")

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 12,
  fig.height = 7,
  dpi = 150,
  out.width = "100%",
  eval = build_local && !cran_build,
  fig.path = "drainage_area_figures/"
)

huc12 <- NULL
wbd_gdb <- "../../Data/final_wbd/WBD_National_GDB/WBD_National_GDB.gdb"
if(file.exists(wbd_gdb) & !cran_build) {
  huc12 <- sf::read_sf(wbd_gdb, "WBDHU12")
  # Rename columns to the schema expected by
  # get_drainage_area_estimates() and the vignette plot helpers.
  rename_map <- c(
    huc12 = "huc_12",
    noncontributingareaacres = "ncontrb_a",
    hutype = "hu_12_type"
  )
  for(src in names(rename_map)) {
    if(src %in% names(huc12) && !(rename_map[[src]] %in% names(huc12)))
      names(huc12)[names(huc12) == src] <- rename_map[[src]]
  }
}

# Local NHDPlusV2 national snapshot. Gated behind USE_LOCAL_NHDP so the
# read only happens when explicitly requested. When enabled and the GDB
# is present, NHDWaterbody and CatchmentSP are passed through
# waterbody_data / catchment_data so the point-in-waterbody lookup
# (Malheur), the gap-zone catchment retrieval, and the full-network
# catchment retrieval are resolved locally instead of via the OGC API.
# NHDFlowline_Network is read for the vignette's flowline-plot fetch.
waterbody_data <- NULL
catchment_data <- NULL
flowlines_data <- NULL
if(Sys.getenv("USE_LOCAL_NHDP") == "TRUE") {
  nhdp_gdb <- paste0("../../Data/nhdp/NHDPlusNationalData/",
    "NHDPlusV21_National_Seamless_Flattened_Lower48.gdb")
  if(file.exists(nhdp_gdb)) {
    waterbody_data <- sf::read_sf(nhdp_gdb, "NHDWaterbody")
    catchment_data <- sf::read_sf(nhdp_gdb, "CatchmentSP")
    flowlines_data <- sf::read_sf(nhdp_gdb, "NHDFlowline_Network")
    flowlines_data <- rename_geometry(flowlines_data, "geom")
    names(flowlines_data) <- tolower(names(flowlines_data))
  }
}

has_davidson <- FALSE
oldoption <- options(scipen = 9999)

A fully rendered version of this vignette with figures is available at https://doi-usgs.github.io/nhdplusTools/articles/drainage_area_estimation.html.

get_drainage_area_estimates() combines Watershed Boundary Dataset (WBD) delineations and non-surface-contributing area estimates with National Hydrography Dataset Plus version 2 (NHDPlusV2) catchment areas for the gap between a basin outlet and the nearest upstream HU12. Non-contributing areas captured in HU12 boundaries are included.

This vignette shows output of the function on five basins that span a range of hydrologic settings -- from basins where topographic signal dominates to arid, glacial, and endorheic settings where it does not.

How `get_drainage_area_estimates()` Works

The function produces drainage area estimates by stitching together two kinds of spatial data: HU12 polygon areas for the bulk of the upstream basin and NHDPlusV2 catchment areas for the gap between the gage (or other outlet) and the nearest upstream HU12 boundaries. Because HU12 polygons carry a non-contributing area attribute (ncontrb_a) and non-surface-contributing HU12s can be retrieved from HU10 or HU8 groupings, the estimates separate total drainage area from surface network contributing drainage area.

Algorithm Steps

Resolve the start feature. A Network Linked Data Index (NLDI) feature list (e.g. featureSource = "nwissite", featureID = "USGS-05406500") is resolved to an NHDPlusV2 COMID via the NLDI. An sfc_POINT inside a waterbody triggers a waterbody lookup to find all outlet flowlines.
Negotiate the outlet catchment. For gage-based starts, the function computes where the gage sits along its outlet flowline as a 0--100 measure. This measure determines whether the outlet catchment needs to be split (see Catchment splitting below).
Fetch the upstream network and HU12 pour points. The full upstream flowline network is retrieved along with all HU12 pour points on that network. These pour points mark the downstream outlet of each HU12 watershed unit that drains to the gage.
Identify immediate HU12 outlets. A network navigation finds which HU12 outlets are directly upstream of the gage with no intervening HU12 outlet. These define the boundary between the "HU zone" (covered by HU12 polygon areas) and the "gap zone" (covered by individual catchment areas).
Fetch and assemble HU12 polygons. HU12 polygons are retrieved at three query levels: the specific NLDI-identified HU12 IDs, all HU12s within upstream HU10 boundaries, and all HU12s within upstream HU8 boundaries. The broader queries can capture HU12s that share a parent HU but lack an on-network pour point -- common in prairie-pothole and playa-dominated landscapes. A disconnect filter removes HU12s pulled incidentally by the broader queries that are not hydrologically connected. Each level is a superset of the narrower level, and each produces both a total and a contributing-only area estimate.
Compute gap area. Catchment areas between the gage and the immediate HU12 outlets are summed. Split-catchment logic is applied at both the HU12 outlets and (optionally) the gage point.
Assemble scalar estimates. Six drainage area estimates are computed: total and contributing at the HU12, HU10, and HU8 query levels. Each equals the HU12 polygon area sum plus the gap area. A seventh scalar, network_da_sqkm, is the NHDPlusV2 cumulative drainage area at the outlet -- a comparison baseline.
NHDPlus High Resolution (NHDPlusHR) estimate (optional). When nhdplushr = TRUE, the function fetches NHDPlusHR flowlines and catchments for the basin bounding box, indexes the start point to the HR network (disambiguating by drainage area), navigates upstream, and sums catchment areas.

Catchment Splitting

When a gage sits partway along a flowline rather than at its outlet, the downstream portion of the catchment does not contribute to the gage. The function calls the NLDI split-catchment service to divide the catchment at the gage point and counts only the upstream portion. The outlet_split_threshold_m parameter (default 100 m) sets the minimum gage-to-outlet distance before splitting is performed; if the gage is closer than this threshold, the full catchment is used.

A parallel split occurs at HU12 outlets: the HU boundary cuts through a catchment, and the portion upstream of the pour point overlaps with HU12 area already counted. The split-catchment service divides this catchment and the upstream overlap is excluded from the gap area to avoid double-counting.

The Outlet Catchment Splitting section later in this vignette illustrates both splits using Davidson Creek (USGS streamgage 08110075).

Running with Local Data

Three parameters reduce dependence on web services:

local_navigation = TRUE loads the NHDPlusV2 Value Added Attributes via get_vaa() for network navigation and flowline attribute lookups. Only HU12 pour points are still fetched from the NLDI.
huc12_data accepts a pre-loaded sf data.frame of HU12 polygons (e.g. from the WBD National Geodatabase). When provided, all HU12 polygon queries are resolved by subsetting this table locally instead of calling WBD web services. The table must include huc_12 and ncontrb_a columns (case-insensitive).
huc12_outlets accepts either a path to a GPKG or a pre-loaded sf data.frame of on-network HU12 pour points. A national CONUS file is available from Blodgett, D.L., 2022, Mainstem Rivers of the Conterminous United States (ver. 3.0, February 2026): U.S. Geological Survey data release, doi:10.5066/P13LNDDQ -- use the hu_points layer of finalwbd_outlets.gpkg. The function renames COMID and FinalWBD_HUC12 to comid and identifier and filters the table to the upstream network automatically. When supplied, no NLDI huc12pp queries are issued.

Using all three together eliminates all WBD and NLDI huc12pp calls and most flowline attribute calls. The offline call looks like:

result <- get_drainage_area_estimates(
  start, local_navigation = TRUE,
  huc12_data = wbd_sf,
  huc12_outlets = "path/to/finalwbd_outlets.gpkg"
)

Web Services and Performance

By default the function contacts several web services. Each adds latency and is subject to rate limits or outages:

NLDI (findNLDI): resolves the start feature, navigates the upstream network, and retrieves HU12 pour points. The huc12pp query is skipped when huc12_outlets is supplied.
NHDPlusV2 OGC API (get_nhdplus): fetches flowline attributes, catchment geometries in the gap zone, and (when catchments = TRUE) the full upstream catchment polygon set. Skipped for flowline attributes when local_navigation = TRUE.
WBD / HU service (get_huc, get_huc12_by_huc): fetches HU12 polygons by ID or by parent HU. Skipped when huc12_data is provided.
Split-catchment service (get_split_catchment): divides catchments at HU12 outlets and at the gage point. Always called.
NHDPlusHR service (get_nhdphr): fetches high-resolution flowlines and catchments for the basin bounding box. Skipped when nhdplushr = FALSE.

For large basins (e.g. the Brazos at Rosharon) the cumulative time for these calls can be substantial. Setting nhdplushr = FALSE eliminates the most expensive single call. Providing huc12_data and using local_navigation = TRUE together can reduce total run time to a fraction of the default, at the cost of requiring local data files. The vignette's fetch_or_load wrapper caches results as RDS files so that repeated runs avoid repeating the web service calls entirely.

Fetch and Cache Results

Each basin result is stored in a separate RDS file in the nhdplusTools data directory so that subsequent runs skip the web service calls.

library(sf)

data_dir <- nhdplusTools_data_dir()
dir.create(data_dir, recursive = TRUE, showWarnings = FALSE)

# On-network HU12 pour points from the Mainstem Rivers data release
# (Blodgett 2022, doi:10.5066/P13LNDDQ). Downloaded once and cached so
# subsequent builds reuse the local file instead of hitting the NLDI.
outlets_path <- file.path(data_dir, "finalwbd_outlets.gpkg")
outlets_url <- paste0(
  "https://prod-is-usgs-sb-prod-publish.s3.amazonaws.com/",
  "65cbc0b3d34ef4b119cb37e9/finalwbd_outlets.gpkg")
if(!file.exists(outlets_path)) {
  message("Downloading HU12 outlets GPKG to ", outlets_path)
  download.file(outlets_url, outlets_path, mode = "wb")
}

# Define the six study sites
sites <- list(
  black_earth = list(featureSource = "nwissite",
    featureID = "USGS-05406500"),
  french_broad = list(featureSource = "nwissite",
    featureID = "USGS-03451500"),
  brazos = list(featureSource = "nwissite",
    featureID = "USGS-08116650"),
  james = list(featureSource = "nwissite",
    featureID = "USGS-06468000"),
  malheur = st_sfc(st_point(c(-118.8, 43.3)), crs = 4326),
  purgatoire = list(featureSource = "nwissite",
    featureID = "USGS-07128500"),
  davidson = list(featureSource = "nwissite",
    featureID = "USGS-08110075")
)

# Skip NHDPlusHR for large basins to keep fetch times reasonable
hr_basins <- c("black_earth", "french_broad", "davidson", "purgatoire",
  "james")

fetch_or_load <- function(name, start) {
  rds_path <- file.path(data_dir,
    paste0("da_est_", name, ".rds"))
  if(file.exists(rds_path)) {
    message("Loading cached: ", name)
    return(readRDS(rds_path))
  }
  message("Fetching: ", name)
  result <- tryCatch(
    get_drainage_area_estimates(start, catchments = TRUE,
      huc12_data = huc12,
      huc12_outlets = outlets_path,
      waterbody_data = waterbody_data,
      catchment_data = catchment_data,
      local_navigation = TRUE,
      nhdplushr = name %in% hr_basins),
    error = function(e) {
      warning("Failed for ", name, ": ", conditionMessage(e),
        call. = FALSE)
      NULL
    })
  if(!is.null(result)) saveRDS(result, rds_path)
  result
}

da_results <- Map(fetch_or_load, names(sites), sites)
da_results <- Filter(Negate(is.null), da_results)

fetch_flowlines <- function(name, da_result) {
  rds_path <- file.path(data_dir, paste0("fl_", name, ".rds"))
  if(file.exists(rds_path)) {
    message("Loading cached flowlines: ", name)
    return(readRDS(rds_path))
  }
  message("Fetching flowlines: ", name)
  comids <- da_result$all_network$comid
  fl <- if(!is.null(flowlines_data)) {
    message("  Subsetting from local NHDFlowline_Network...")
    flowlines_data[flowlines_data$comid %in% comids, ]
  } else {
    tryCatch(
      get_nhdplus(comid = comids, realization = "flowline"),
      error = function(e) {
        warning("Flowline fetch failed for ", name, ": ",
          conditionMessage(e), call. = FALSE)
        NULL
      })
  }
  if(!is.null(fl) && nrow(fl) > 0) saveRDS(fl, rds_path)
  fl
}

flowlines <- Map(fetch_flowlines, names(da_results), da_results)
flowlines <- Filter(Negate(is.null), flowlines)

# Fetch NWIS drainage area for sites with an nwissite featureSource
nwis_ids <- vapply(sites, function(s) {
  if(is.list(s) && identical(s$featureSource, "nwissite"))
    s$featureID else NA_character_
}, character(1))
nwis_ids <- nwis_ids[!is.na(nwis_ids)]

nwis_da <- if(length(nwis_ids) > 0) {
  tryCatch({
    ml <- dataRetrieval::read_waterdata_monitoring_location(
      unname(nwis_ids))
    # Match returned rows back to site names by monitoring_location_id
    idx <- match(nwis_ids, ml$monitoring_location_id)
    # drainage_area and contributing_drainage_area are in sq miles
    data.frame(
      name = names(nwis_ids),
      nwis_da_sqmi = ml$drainage_area[idx],
      nwis_contrib_da_sqmi = ml$contributing_drainage_area[idx],
      nwis_da_sqkm = ml$drainage_area[idx] * 2.58999,
      nwis_contrib_da_sqkm = ml$contributing_drainage_area[idx] *
        2.58999,
      stringsAsFactors = FALSE
    )
  }, error = function(e) {
    warning("NWIS site fetch failed: ", conditionMessage(e),
      call. = FALSE)
    NULL
  })
} else NULL

Return Structure

Each result is a list with scalar drainage area estimates and spatial data frames. Here are the elements for Black Earth Creek:

names(da_results$black_earth)

The scalar estimates (square kilometers) across all basins:

basin_labels <- c(
  black_earth = "Black Earth Creek",
  french_broad = "French Broad",
  brazos = "Brazos at Rosharon",
  james = "James River",
  malheur = "Malheur Lake",
  davidson = "Davidson Creek",
  purgatoire = "Purgatoire River"
)

summary_df <- data.frame(
  basin = basin_labels[names(da_results)],
  network_da = vapply(da_results, \(x) x$network_da_sqkm, numeric(1)),
  da_huc12 = vapply(da_results, \(x) x$da_huc12_sqkm, numeric(1)),
  da_huc10 = vapply(da_results,
    \(x) ifelse(is.na(x$da_huc10_sqkm), NA_real_, x$da_huc10_sqkm),
    numeric(1)),
  da_huc08 = vapply(da_results,
    \(x) ifelse(is.na(x$da_huc08_sqkm), NA_real_, x$da_huc08_sqkm),
    numeric(1)),
  nhdplushr = vapply(da_results,
    \(x) ifelse(is.na(x$nhdplushr_network_dasqkm), NA_real_,
      x$nhdplushr_network_dasqkm),
    numeric(1))
)

# Add NWIS drainage areas where available
if(!is.null(nwis_da)) {
  summary_df$nwis_da <- ifelse(
    names(da_results) %in% nwis_da$name,
    nwis_da$nwis_da_sqkm[match(names(da_results), nwis_da$name)],
    NA_real_)
  summary_df$nwis_contrib_da <- ifelse(
    names(da_results) %in% nwis_da$name,
    nwis_da$nwis_contrib_da_sqkm[match(names(da_results), nwis_da$name)],
    NA_real_)
} else {
  summary_df$nwis_da <- NA_real_
  summary_df$nwis_contrib_da <- NA_real_
}

knitr::kable(summary_df, digits = 1,
  col.names = c("Basin", "Network DA", "HU12 DA",
    "HU10 DA", "HU8 DA", "NHDPlusHR DA",
    "NWIS DA", "NWIS Contributing DA"),
  caption = "Drainage area estimates (sq km) by basin and method")

library(ggplot2)

# USGS USTopo provider for maptiles. The "?f=jpg" suffix is a no-op
# query parameter that gives maptiles a recognized extension hint.
usgs_topo_provider <- list(
  src = "USGS_USTopo",
  q = paste0("https://basemap.nationalmap.gov/arcgis/rest/services/",
    "USGSTopo/MapServer/tile/{z}/{y}/{x}?f=jpg"),
  sub = NA,
  cit = "USGS The National Map: USGSTopo"
)

# Cache directory for tiles
tile_cache <- file.path(nhdplusTools_data_dir(), "da_v_tiles")
dir.create(tile_cache, recursive = TRUE, showWarnings = FALSE)

# Compute a target tile zoom level so the rendered tiles are sharp
# at the figure size we use. Web Mercator world width is 40075016 m
# and tile size is 256 px. The figure pixel width drives the
# resolution we want from the tile mosaic.
target_zoom <- function(sf_obj, fig_px = 1700) {
  bb <- st_bbox(sf_obj)
  span_m <- max(bb["xmax"] - bb["xmin"], bb["ymax"] - bb["ymin"])
  if(!is.finite(span_m) || span_m <= 0) return(10L)
  z <- log2(40075016 * fig_px / (256 * span_m))
  as.integer(min(max(round(z), 4L), 15L))
}

# Fetch USGS topo tiles for an sf bbox. Returns a list with a color
# matrix (suitable for annotation_raster) and the tile bounding box.
# Uses annotation_raster (not geom_raster) so the fill aesthetic
# stays free for choropleth use. When fig_ratio is provided, the
# bbox is expanded to match the figure aspect ratio via
# fit_bbox_to_fig so that tiles cover the full canvas.
fetch_topo_raster <- function(sf_obj, zoom = NULL,
  fig_ratio = NULL) {
  fit <- if(!is.null(fig_ratio))
    fit_bbox_to_fig(sf_obj, fig_ratio = fig_ratio) else NULL
  if(!is.null(fit)) {
    bb <- st_bbox(c(
      xmin = fit$xlim[1], xmax = fit$xlim[2],
      ymin = fit$ylim[1], ymax = fit$ylim[2]
    ), crs = st_crs(3857))
    bb_sf <- st_as_sfc(bb)
  } else {
    bb <- st_bbox(sf_obj)
    pad_x <- (bb["xmax"] - bb["xmin"]) * 0.05
    pad_y <- (bb["ymax"] - bb["ymin"]) * 0.05
    bb["xmin"] <- bb["xmin"] - pad_x
    bb["xmax"] <- bb["xmax"] + pad_x
    bb["ymin"] <- bb["ymin"] - pad_y
    bb["ymax"] <- bb["ymax"] + pad_y
    bb_sf <- st_as_sfc(bb, crs = st_crs(sf_obj))
  }
  if(is.null(zoom))
    zoom <- target_zoom(st_transform(bb_sf, 3857))
  args <- list(
    x = bb_sf, provider = usgs_topo_provider, crop = TRUE,
    cachedir = tile_cache, zoom = zoom
  )
  tiles <- tryCatch(
    do.call(maptiles::get_tiles, args),
    error = function(e) {
      warning("Tile fetch failed: ", conditionMessage(e),
        call. = FALSE)
      NULL
    })
  if(is.null(tiles)) return(NULL)
  rgb_arr <- terra::as.array(tiles)
  if(length(dim(rgb_arr)) != 3 || dim(rgb_arr)[3] < 3) return(NULL)
  col_mat <- matrix(
    rgb(rgb_arr[, , 1], rgb_arr[, , 2], rgb_arr[, , 3],
      maxColorValue = 255),
    nrow = dim(rgb_arr)[1], ncol = dim(rgb_arr)[2]
  )
  ext <- terra::ext(tiles)
  list(
    raster = col_mat,
    xmin = ext[1], xmax = ext[2],
    ymin = ext[3], ymax = ext[4]
  )
}

# Expand a bbox to match the figure aspect ratio so that coord_sf
# fills the canvas without white sidebars. The basin content is
# centered and the narrower axis is padded with additional basemap.
fit_bbox_to_fig <- function(sf_obj, fig_ratio = 12 / 7, pad = 0.08) {
  bb <- st_bbox(st_transform(sf_obj, 3857))
  dx <- unname(bb["xmax"] - bb["xmin"])
  dy <- unname(bb["ymax"] - bb["ymin"])
  if(!is.finite(dx) || !is.finite(dy) || dx <= 0 || dy <= 0)
    return(NULL)
  # Add uniform padding first
  cx <- unname((bb["xmin"] + bb["xmax"]) / 2)
  cy <- unname((bb["ymin"] + bb["ymax"]) / 2)
  dx <- dx * (1 + pad * 2)
  dy <- dy * (1 + pad * 2)
  # Expand the narrower axis to match the target aspect ratio
  current_ratio <- dx / dy
  if(current_ratio < fig_ratio) {
    dx <- dy * fig_ratio
  } else {
    dy <- dx / fig_ratio
  }
  list(
    xlim = c(cx - dx / 2, cx + dx / 2),
    ylim = c(cy - dy / 2, cy + dy / 2)
  )
}

# Add basemap as annotation_raster (does not consume fill scale)
add_topo <- function(p, sf_obj, fig_ratio = 12 / 7) {
  topo <- fetch_topo_raster(sf_obj, fig_ratio = fig_ratio)
  if(is.null(topo)) return(p)
  p + annotation_raster(
    topo$raster,
    xmin = topo$xmin, xmax = topo$xmax,
    ymin = topo$ymin, ymax = topo$ymax
  )
}

# Resolve gage point in target CRS
get_gage_pt <- function(da, target_crs) {
  pt <- st_geometry(da$start_feature)
  if(!inherits(pt, "sfc_POINT"))
    pt <- st_centroid(pt)
  st_transform(st_sf(geometry = pt), target_crs)
}

# Pick the spatially-largest of the available HUC12 sets returned by
# get_drainage_area_estimates: hu12_by_huc08 (if present), otherwise
# hu12_by_huc10, otherwise hu12_by_huc12. Largest is determined by
# the sum of HUC12 dasqkm in each set.
largest_hu12_set <- function(da) {
  candidates <- list(
    huc12 = da$hu12_by_huc12,
    huc10 = da$hu12_by_huc10,
    huc08 = da$hu12_by_huc08
  )
  candidates <- Filter(
    function(x) !is.null(x) && nrow(x) > 0, candidates)
  areas <- vapply(candidates,
    function(x) sum(x$dasqkm, na.rm = TRUE), numeric(1))
  candidates[[which.max(areas)]]
}

# Compute a single merged basin polygon for a given HUC level set
# (hu12_by_huc12, hu12_by_huc10, or hu12_by_huc08), unioned with
# extra_catchments and the local catchments around the gage.
# Uses hydroloom::dissolve_polygons() with gap_tolerance to absorb
# sliver gaps between polygons from different services (WBD HUC12s,
# NHDPlusV2 catchments, split-catchment service).
compute_basin_polygon <- function(da, hu_layer, target_crs,
  tol_m = 250) {
  common_crs <- st_crs(hu_layer)
  pieces <- list()
  pieces$hu <- st_geometry(hu_layer)

  if(!is.null(da$extra_catchments) && nrow(da$extra_catchments) > 0)
    pieces$extra <- st_geometry(
      st_transform(da$extra_catchments, common_crs))

  if(!is.null(da$split_catchment) && nrow(da$split_catchment) > 0) {
    sc <- da$split_catchment
    sc_full <- sc[sc$id == "catchment", ]
    if(nrow(sc_full) > 0)
      pieces$catch <- st_geometry(st_transform(sc_full, common_crs))
  }

  all_polys <- st_sf(
    geometry = do.call(c, pieces)
  ) |> st_set_crs(common_crs)

  hydroloom::dissolve_polygons(all_polys,
    gap_tolerance = tol_m,
    max_hole_area = Inf,
    single_polygon = TRUE) |>
    st_transform(target_crs)
}

# Common theme: no axis ticks, labels, or grid; legend sized for
# readability
map_theme <- function() {
  theme_void() +
    theme(
      legend.position = "bottom",
      legend.box = "horizontal",
      legend.text = element_text(size = 11),
      legend.title = element_text(size = 12),
      legend.key.size = unit(0.6, "cm"),
      plot.title = element_blank(),
      plot.margin = margin(2, 2, 2, 2)
    )
}

# Add the gray basin underlay to a plot. The polygon represents the
# full estimated drainage area envelope (largest available HUC level
# plus extra catchments and local catchments).
add_basin_underlay <- function(p, basin_poly) {
  p + geom_sf(data = basin_poly,
    fill = "gray30", color = NA, alpha = 0.15,
    inherit.aes = FALSE)
}

# Mode 1: Drainage Area Boundaries
plot_boundaries <- function(da, basin_label) {
  # Reproject working layers to Web Mercator (3857) so the basemap
  # tiles align cleanly.
  target_crs <- st_crs(3857)
  gage_pt <- get_gage_pt(da, target_crs)
  basin_poly <- compute_basin_polygon(da, largest_hu12_set(da), target_crs)

  # Dissolve boundaries at each HU level. Each HU level is combined
  # with extra_catchments and split_catchment via compute_basin_polygon
  # so every boundary reaches the gage instead of stopping at the
  # lowest complete HU outlet.
  boundaries <- list()
  boundaries[["HU12"]] <- st_geometry(
    compute_basin_polygon(da, da$hu12_by_huc12, target_crs))
  if(!is.null(da$hu12_by_huc10))
    boundaries[["HU10"]] <- st_geometry(
      compute_basin_polygon(da, da$hu12_by_huc10, target_crs))
  if(!is.null(da$hu12_by_huc08))
    boundaries[["HU8"]] <- st_geometry(
      compute_basin_polygon(da, da$hu12_by_huc08, target_crs))
  if(!is.null(da$all_catchments))
    boundaries[["Network catchments (NHDPlusV2)"]] <- st_union(
      st_transform(da$all_catchments, target_crs))
  if(!is.null(da$nhdplushr_boundary))
    boundaries[["NHDPlusHR"]] <- st_transform(
      da$nhdplushr_boundary, target_crs)

  # Linetype, color, and linewidth all keyed off the source factor so
  # one merged "DA boundary" legend distinguishes the layers.
  # Network catchments (NHDPlusV2) gets a saturated magenta dotted
  # line so it reads distinctly against the basemap and against HU10's
  # longdash. NHDPlusHR uses orange to separate from the HU layers.
  all_ltypes <- c(
    "HU12" = "solid",
    "HU10" = "longdash",
    "HU8" = "dotdash",
    "Network catchments (NHDPlusV2)" = "22",
    "NHDPlusHR" = "solid"
  )
  all_colors <- c(
    "HU12" = "gray10",
    "HU10" = "gray10",
    "HU8" = "gray10",
    "Network catchments (NHDPlusV2)" = "magenta",
    "NHDPlusHR" = "#D55E00"
  )
  all_lwds <- c(
    "HU12" = 0.8,
    "HU10" = 0.8,
    "HU8" = 0.8,
    "Network catchments (NHDPlusV2)" = 0.5,
    "NHDPlusHR" = 0.6
  )
  ltypes  <- all_ltypes[names(boundaries)]
  lcolors <- all_colors[names(boundaries)]
  lwds    <- all_lwds[names(boundaries)]

  bnd_sf <- do.call(rbind, lapply(seq_along(boundaries), function(i) {
    st_sf(source = factor(names(boundaries)[i],
        levels = names(boundaries)),
      geometry = boundaries[[i]])
  }))

  hu12_pp <- if(!is.null(da$hu12_outlet))
    st_transform(da$hu12_outlet, target_crs) else NULL

  # Combine gage + HU12 outlets into one points layer with a single
  # shape scale. Keeping shape only (no color aes) frees the color
  # scale for the boundary symbology.
  pts_list <- list(
    st_sf(label = "Stream gage / outlet",
      geometry = st_geometry(gage_pt))
  )
  if(!is.null(hu12_pp))
    pts_list <- c(pts_list, list(
      st_sf(label = "HU12 outlet",
        geometry = st_geometry(hu12_pp))))
  pts_sf <- do.call(rbind, pts_list)
  pts_sf$label <- factor(pts_sf$label,
    levels = c("Stream gage / outlet", "HU12 outlet"))

  fig_lims <- fit_bbox_to_fig(basin_poly)

  p <- ggplot() |>
    add_topo(basin_poly) |>
    add_basin_underlay(basin_poly)

  # HU12 outlets drawn first so basin boundaries overlay them.
  if(!is.null(hu12_pp))
    p <- p + geom_sf(
      data = pts_sf[pts_sf$label == "HU12 outlet", ],
      shape = 4, size = 1, color = "darkred", stroke = 0.6,
      alpha = 0.75, inherit.aes = FALSE)

  p <- p +
    geom_sf(data = bnd_sf,
      aes(linetype = source, color = source, linewidth = source),
      fill = NA, inherit.aes = FALSE) +
    scale_linetype_manual(values = ltypes,  name = "DA boundary") +
    scale_color_manual(   values = lcolors, name = "DA boundary") +
    scale_linewidth_manual(values = lwds,   name = "DA boundary") +
    geom_sf(
      data = pts_sf[pts_sf$label == "Stream gage / outlet", ],
      shape = 17, size = 5, color = "black", fill = "white",
      stroke = 1.4, inherit.aes = FALSE) +
    coord_sf(crs = target_crs, expand = FALSE,
      xlim = fig_lims$xlim, ylim = fig_lims$ylim) +
    labs(title = paste(basin_label, "-- Drainage Area Boundaries")) +
    map_theme()

  p
}

# Mode 2: Stream Network
plot_network <- function(da, fl, basin_label) {
  target_crs <- st_crs(3857)
  gage_pt <- get_gage_pt(da, target_crs)
  fl <- st_transform(fl, target_crs)
  basin_poly <- compute_basin_polygon(da, largest_hu12_set(da), target_crs)

  # fcode 46003 = intermittent, 46007 = ephemeral, 46006 = perennial
  fl$flow_class <- ifelse(
    fl$fcode %in% c(46003L, 46007L), "Intermittent/Ephemeral",
    "Perennial/Other"
  )
  fl$lwd <- pmin(pmax(fl$streamorde, 1) / 3, 1.5)

  pts <- st_sf(label = "Stream gage / outlet",
    geometry = st_geometry(gage_pt))

  fig_lims <- fit_bbox_to_fig(basin_poly)

  ggplot() |>
    add_topo(basin_poly) |>
    add_basin_underlay(basin_poly) +
    geom_sf(data = fl,
      aes(color = flow_class, linewidth = lwd),
      show.legend = "line", inherit.aes = FALSE) +
    scale_color_manual(
      values = c("Perennial/Other" = "steelblue",
        "Intermittent/Ephemeral" = "lightskyblue"),
      name = "Flow class") +
    scale_linewidth_identity() +
    ggnewscale_color(pts, target_crs) +
    coord_sf(crs = target_crs, expand = FALSE,
      xlim = fig_lims$xlim, ylim = fig_lims$ylim) +
    labs(title = paste(basin_label, "-- Stream Network")) +
    map_theme()
}

# Brazos-specific variant of plot_network that adds a Caprock Escarpment
# trace and label per LH-217625600 / LH-2005481448. The escarpment is
# the eastern margin of the Llano Estacado and the physiographic
# boundary between the High Plains and the Rolling Plains (Wermund,
# 1996). The polyline below approximates its trace through the Brazos
# basin's northwestern arm.
plot_network_brazos <- function(da, fl, basin_label) {
  target_crs <- st_crs(3857)
  caprock_coords <- matrix(c(
    -101.05, 34.50,
    -101.00, 33.70,
    -100.95, 32.85
  ), ncol = 2, byrow = TRUE)
  caprock_line <- sf::st_transform(
    sf::st_sfc(sf::st_linestring(caprock_coords), crs = 4326),
    target_crs)
  caprock_anchor <- sf::st_transform(
    sf::st_sfc(sf::st_point(c(-101.05, 34.50)), crs = 4326),
    target_crs)
  caprock_line_sf <- sf::st_sf(geometry = caprock_line)
  caprock_label_sf <- sf::st_sf(
    label = "Caprock Escarpment",
    geometry = caprock_anchor)

  plot_network(da, fl, basin_label) +
    geom_sf(data = caprock_line_sf,
      color = "black", linewidth = 0.8,
      inherit.aes = FALSE) +
    geom_sf_text(data = caprock_label_sf, aes(label = label),
      size = 3.8, fontface = "italic", color = "gray15",
      hjust = 0, nudge_x = 8000, nudge_y = 22000,
      inherit.aes = FALSE)
}

# Helper to add the gage point as a labeled layer using a shape scale
# that doesn't conflict with other color scales already in use.
ggnewscale_color <- function(pts, target_crs) {
  list(
    geom_sf(data = pts, aes(shape = label),
      size = 5, color = "black", fill = "white", stroke = 1.4,
      inherit.aes = FALSE),
    scale_shape_manual(
      values = c("Stream gage / outlet" = 17),
      name = NULL)
  )
}

# Mode 3: HUC12 Type (hu_12_type)
plot_type <- function(da, fl, basin_label,
  fig_ratio = 12 / 8) {
  target_crs <- st_crs(3857)
  largest <- largest_hu12_set(da)
  hu12 <- st_transform(largest, target_crs)

  type_labels <- c(
    "S" = "Standard",
    "C" = "Closed basin",
    "M" = "Multiple outlets",
    "F" = "Frontal",
    "W" = "Water",
    "D" = "Drainage",
    "I" = "Island"
  )
  hu12$type_label <- factor(
    type_labels[hu12$hu_12_type],
    levels = c("Standard", "Closed basin", "Multiple outlets",
      "Frontal", "Water", "Drainage", "Island")
  )

  gage_pt <- get_gage_pt(da, target_crs)
  fl <- st_transform(fl, target_crs)
  basin_poly <- compute_basin_polygon(da, largest, target_crs)

  pts <- st_sf(label = "Stream gage / outlet",
    geometry = st_geometry(gage_pt))

  fig_lims <- fit_bbox_to_fig(basin_poly, fig_ratio = fig_ratio)

  # 30% alpha baked into the fill so the legend swatches match the map.
  # Standard is light grey so the legend swatch reads against white.
  col30 <- function(col) adjustcolor(col, alpha.f = 0.3)
  type_colors <- c(
    "Standard" = adjustcolor("gray85", alpha.f = 0.4),
    "Closed basin" = col30("#8B4513"),
    "Multiple outlets" = col30("#555555"),
    "Frontal" = col30("#228B22"),
    "Water" = col30("#1E90FF"),
    "Drainage" = col30("#9467BD"),
    "Island" = col30("#FFC107")
  )

  p <- ggplot() |>
    add_topo(basin_poly, fig_ratio = fig_ratio) |>
    add_basin_underlay(basin_poly) +
    geom_sf(data = fl, color = "steelblue", linewidth = 0.25,
      alpha = 0.5, inherit.aes = FALSE) +
    geom_sf(data = hu12,
      aes(fill = type_label),
      color = "gray50", linewidth = 0.5,
      inherit.aes = FALSE) +
    scale_fill_manual(
      values = type_colors,
      drop = TRUE,
      name = "HU12 type")

  p +
    ggnewscale_color(pts, target_crs) +
    coord_sf(crs = target_crs, expand = FALSE,
      xlim = fig_lims$xlim, ylim = fig_lims$ylim) +
    labs(title = paste(basin_label, "-- HU12 Types")) +
    map_theme() +
    theme(legend.box = "vertical") +
    guides(
      fill = guide_legend(
        nrow = 1, order = 1,
        override.aes = list(color = "gray50", linewidth = 0.5)))
}

# Per-basin summary table
basin_summary_table <- function(da, basin_name = NULL) {
  vals <- data.frame(
    source = c("Network", "HU12",
      "HU10", "HU8", "NHDPlusHR"),
    total_sqkm = c(
      da$network_da_sqkm,
      da$da_huc12_sqkm,
      ifelse(is.na(da$da_huc10_sqkm), NA_real_, da$da_huc10_sqkm),
      ifelse(is.na(da$da_huc08_sqkm), NA_real_, da$da_huc08_sqkm),
      ifelse(is.na(da$nhdplushr_network_dasqkm), NA_real_,
        da$nhdplushr_network_dasqkm)
    )
  )
  # Add NWIS rows if this basin has NWIS data
  if(!is.null(nwis_da) && !is.null(basin_name) &&
    basin_name %in% nwis_da$name) {
    row <- nwis_da[nwis_da$name == basin_name, ]
    vals <- rbind(vals, data.frame(
      source = c("NWIS", "NWIS contributing"),
      total_sqkm = c(row$nwis_da_sqkm, row$nwis_contrib_da_sqkm)
    ))
  }
  knitr::kable(vals, digits = 1,
    col.names = c("Source", "Area (sq km)"),
    caption = "Drainage area estimates by source")
}

Basin Vignettes {.tabset}

French Broad River at Asheville

Well-determined basin with strong topographic signal (USGS streamgage 03451500). Contributing area essentially equals total drainage area across all sources.

Drainage Area Boundaries

plot_boundaries(da_results$french_broad, "French Broad")

basin_summary_table(da_results$french_broad, "french_broad")

::: {.vignette-note} - All boundary sources converge tightly around the same watershed outline. - The HU12-, HU10-, and HU8-derived boundaries are nearly identical, consistent with a basin whose divides are topographically unambiguous. - The gage (triangle) sits on the main stem French Broad River at Asheville. :::

Stream Network

plot_network(da_results$french_broad, flowlines$french_broad,
  "French Broad")

::: {.vignette-note} - The network is predominantly perennial with dense tributary coverage in the Appalachian headwaters. - Intermittent reaches are sparse, concentrated in low-order headwater channels. - HU12 outlets (x markers) align with major tributary junctions. :::

HU12 Types

plot_type(da_results$french_broad, flowlines$french_broad,
  "French Broad")

::: {.vignette-note} - Every HU12 in the French Broad basin is classified as Standard. - The uniform light-grey fill confirms that all units link to a downstream HU12 through normal surface drainage. :::

Brazos River, West Texas

Arid/ephemeral connectivity basin (USGS streamgage 08116650). Transmission losses and disconnected uplands create large differences between total and contributing drainage area.

Drainage Area Boundaries

plot_boundaries(da_results$brazos, "Brazos at Rosharon")

basin_summary_table(da_results$brazos, "brazos")

::: {.vignette-note} - Boundary sources diverge substantially. HU-derived boundaries extend further west into arid uplands than the network-derived boundary. - The gap between HU12 total area and contributing area reflects large noncontributing designations in western sub-basins. :::

Stream Network

plot_network_brazos(da_results$brazos, flowlines$brazos,
  "Brazos at Rosharon")

::: {.vignette-note} - Extensive intermittent and ephemeral reaches dominate the western (upstream) portion of the basin. - Perennial flow is concentrated in the lower main stem and major tributaries east of the Caprock Escarpment, the physiographic boundary between the High Plains (Llano Estacado) and the Rolling Plains traced on the figure. - The transition from ephemeral headwaters to perennial mainstem illustrates why contributing area is smaller than total area. :::

HU12 Types

plot_type(da_results$brazos, flowlines$brazos,
  "Brazos at Rosharon")

::: {.vignette-note} - Closed-basin HU12s (brown) concentrate in the arid western uplands of the Llano Estacado, where surface water never reaches the Brazos. - Multiple-outlet HU12s (grey) have more than one outlet to the same downstream system — for example, a HU12 with a braided channel where flow leaves through more than one point. - Frontal HU12s (green), where present, are coastal units that drain directly to the Gulf of Mexico rather than to a downstream HU12. - Standard HU12s (light grey) dominate the perennial eastern main stem. :::

James River / Cottonwood Lake

Glacial prairie basin (USGS streamgage 06468000). Prairie potholes create a large noncontributing fraction that varies by sub-basin.

Drainage Area Boundaries

plot_boundaries(da_results$james, "James River")

basin_summary_table(da_results$james, "james")

::: {.vignette-note} - Boundary estimates diverge in the upper basin where glacial topography creates ambiguous divides. - The HU12-derived total area exceeds the contributing area by a notable margin, reflecting prairie pothole storage. :::

Stream Network

plot_network(da_results$james, flowlines$james, "James River")

::: {.vignette-note} - The network includes a mix of perennial and intermittent reaches. - Intermittent channels are common in the upper basin where glacial drift creates closed depressions and episodic connectivity. - The lower main stem is perennial and well-defined. :::

HU12 Types

plot_type(da_results$james, flowlines$james, "James River")

::: {.vignette-note} - Two closed-basin HU12s (brown) appear in the upper prairie-pothole headwaters — these are units with no surface outlet to the James. - The remaining HU12s are Standard (light grey), including the lower perennial main stem. :::

Black Earth Creek

Small Driftless-Area basin in southern Wisconsin (USGS streamgage 05406500). Well-determined drainage area with minimal noncontributing fraction.

Drainage Area Boundaries

plot_boundaries(da_results$black_earth, "Black Earth Creek")

basin_summary_table(da_results$black_earth, "black_earth")

::: {.vignette-note} - All boundary sources converge on a compact watershed outline. - The basin is small enough to fall within a single HU10, so HU10- and HU8-level boundaries are not computed separately. :::

Stream Network

plot_network(da_results$black_earth, flowlines$black_earth,
  "Black Earth Creek")

::: {.vignette-note} - A short, predominantly perennial network drains the Driftless Area landscape. - Few intermittent headwater reaches are present. :::

HU12 Types

plot_type(da_results$black_earth, flowlines$black_earth,
  "Black Earth Creek")

::: {.vignette-note} - The single HU12 covering Black Earth Creek is Standard type (light-grey fill), consistent with a well-defined basin in Driftless-Area terrain. :::

Malheur Lake

Malheur Lake basin is an endorheic (closed) basin with no surface outlet. The terminal feature is a lake rather than a stream gage.

Because the start point falls in a closed-basin HU12 (Malheur Lake, 171200010710) with no flowline pour point on the NHD network, get_drainage_area_estimates() auto-promotes the parent HUC10 (1712000107) into HU_inclusion_override. The HU08 estimate then spans the full Harney-Malheur HU8 17120001 -- including the lake itself, Harney Lake, and the frontal HU12s along the lake margins -- rather than only the on-network HU12s reachable from the inflow tributaries.

Drainage Area Boundaries

plot_boundaries(da_results$malheur, "Malheur Lake")

basin_summary_table(da_results$malheur, "malheur")

::: {.vignette-note} - The basin is endorheic -- all boundaries terminate at Malheur Lake with no downstream outlet. - HU-derived and network-derived boundaries diverge in the surrounding high-desert uplands where surface connectivity is ambiguous. - The shaded grey area is the merged basin envelope used as a background underlay, not an HU8 boundary. - The start feature (triangle) marks the centroid of the Malheur Lake waterbody polygon rather than a stream gage. :::

Stream Network

plot_network(da_results$malheur, flowlines$malheur, "Malheur Lake")

::: {.vignette-note} - Intermittent and ephemeral reaches dominate the network, particularly in the southern and eastern tributaries. - Perennial flow is limited to the Silvies River and Donner und Blitzen River corridors draining into Malheur and Harney Lakes. - The network terminates at the lake with no surface outlet downstream. :::

HU12 Types

plot_type(da_results$malheur, flowlines$malheur, "Malheur Lake")

::: {.vignette-note} - Closed-basin HU12s (brown) in HU8 17120001 -- Malheur Lake itself, Harney Lake, Sunset Valley, and Lower Riddle Creek -- mark Harney-Malheur units where surface drainage terminates internally. These are now included in the HU08 estimate via the closed-basin auto-override. - The brown HU12s south and southwest of the lake (Hay Lake, Capehart Lake, Mule Springs Valley, Smoky Hollow, Lake On The Trail, Foster Lake) belong to HU8 17120004 Silver Creek-Silver Lake -- a separate endorheic basin pulled into the HU08 estimate because Silvies/Donner upstream HU12s share its parent HU8. - Water HU12s (blue) — for example, Mud Lake — designate large open water bodies coded as Water type in the current WBD. - Multiple-outlet HU12s (grey) appear on the low-gradient flats of the Harney Basin. - The surrounding Standard HU12s (light grey) supply the Silvies and Donner und Blitzen corridors that actually reach the lake. :::

Purgatoire River near Las Animas

Boundary placement uncertainty basin in southeastern Colorado (USGS streamgage 07128500). Flat terrain adjacent to the basin produces a drainage divide that different delineation methods place in different locations, resulting in divergent drainage area estimates for the gage and all downstream stations. Documented in Dupree and Crowfoot (2012, TM 11-C6).

Drainage Area Boundaries

plot_boundaries(da_results$purgatoire, "Purgatoire River")

basin_summary_table(da_results$purgatoire, "purgatoire")

::: {.vignette-note} - Boundary sources diverge in the flat terrain along the western and southern margins of the basin where the drainage divide is poorly defined. - Manual delineation from USGS topographic maps assigned a noncontributing area adjacent to the basin as part of the Purgatoire drainage; the WBD placed that area in the neighboring hydrologic unit to the west. - The resulting differences propagate to all downstream gages on the Purgatoire and Arkansas Rivers. :::

Stream Network

plot_network(da_results$purgatoire, flowlines$purgatoire,
  "Purgatoire River")

::: {.vignette-note} - The upper basin drains steep terrain along the Sangre de Cristo Range with predominantly perennial flow. - The lower basin crosses the high plains where intermittent and ephemeral channels are more common and terrain gradients weaken. - The transition from montane headwaters to plains illustrates why the divide becomes uncertain in the lower-gradient portions of the basin. :::

HU12 Types

plot_type(da_results$purgatoire, flowlines$purgatoire,
  "Purgatoire River")

::: {.vignette-note} - Every HU12 in the Purgatoire basin is Standard type (light-grey fill). - The boundary-placement-uncertainty story here is about where the WBD places the divide between Standard HU12s, not about closed-basin or multiple-outlet designations within the basin. :::

Outlet Catchment Splitting

When a gage sits partway along a flowline rather than at its outlet, the gage's outlet catchment is only partly upstream of the gage. The downstream portion does not contribute to the gage and should be excluded. A similar situation arises at HU12 outlets: the HU boundary cuts through a catchment and the portion upstream of the pour point overlaps with the HU area that is already counted.

get_drainage_area_estimates() handles both splits automatically. The outlet_split_threshold_m parameter (default 100 m) controls the minimum gage-to-outlet distance before the split is performed.

This example uses a small Texas gage (USGS streamgage 08110075, Davidson Creek) where the gage falls roughly two thirds of the way up its flowline.

da_dav <- da_results$davidson
fl_dav <- flowlines$davidson
has_davidson <- !is.null(da_dav) && !is.null(fl_dav) &&
  !is.null(da_dav$outlet_split_catchment)

if(has_davidson) {
  target_crs <- st_crs(3857)

  gage_pt <- get_gage_pt(da_dav, target_crs)
  gage_comid <- as.integer(da_dav$start_feature$comid)

  # all catchments in the gap between the gage and HUC12 outlets
  extra_cat <- st_transform(da_dav$extra_catchments, target_crs)
  all_cat <- st_transform(da_dav$all_catchments, target_crs)
  fl_sf <- st_transform(fl_dav, target_crs)
  hu12_pts <- st_transform(da_dav$hu12_outlet, target_crs)

  # HUC12 split catchment (split at the HUC12 pour point)
  hu12_sc <- st_transform(da_dav$split_catchment, target_crs)
  hu12_catch_full <- hu12_sc[hu12_sc$id == "catchment", ]
  hu12_catch_split <- hu12_sc[hu12_sc$id == "splitCatchment", ]
  hu12_comid <- as.integer(na.omit(hu12_sc$catchmentID))

  # gage outlet split catchment (split at the gage point)
  osc <- st_transform(da_dav$outlet_split_catchment, target_crs)
  osc_full <- osc[osc$id == "catchment", ]
  osc_split <- osc[osc$id == "splitCatchment", ]

  # flowlines in the gap area
  gap_comids <- c(gage_comid, extra_cat$featureid, hu12_comid)
  gap_fl <- fl_sf[fl_sf$comid %in% gap_comids, ]

  # focus bbox: union of gap catchments + split catchments + gage
  focus_geom <- st_union(c(
    st_geometry(extra_cat),
    st_geometry(hu12_catch_full),
    st_geometry(osc_full),
    st_geometry(gage_pt)
  ))
  focus_bb <- st_bbox(focus_geom)
  pad_x <- (focus_bb["xmax"] - focus_bb["xmin"]) * 0.05
  pad_y <- (focus_bb["ymax"] - focus_bb["ymin"]) * 0.05
  focus_xlim <- c(focus_bb["xmin"] - pad_x, focus_bb["xmax"] + pad_x)
  focus_ylim <- c(focus_bb["ymin"] - pad_y, focus_bb["ymax"] + pad_y)
}

Overview: Gage and HU Outlet Positions

The gage (triangle) and the nearest HU12 outlet (x) sit on different catchments with a gap between them. The gap catchments (light blue) and split catchment boundaries define the area between the gage and the HU12 boundary.

p_overview <- ggplot() |>
  add_topo(focus_geom) +
  # all catchments as light underlay
  geom_sf(data = all_cat, fill = "gray30", color = NA, alpha = 0.12,
    inherit.aes = FALSE) +
  # gap catchments
  geom_sf(data = extra_cat, fill = "lightblue", color = "steelblue",
    linewidth = 0.3, alpha = 0.5, inherit.aes = FALSE) +
  # HUC12 split catchment — full outline
  geom_sf(data = hu12_catch_full, fill = NA, color = "gray10",
    linewidth = 0.6, linetype = "dashed", inherit.aes = FALSE) +
  # gage outlet catchment — full outline
  geom_sf(data = osc_full, fill = NA, color = "gray10",
    linewidth = 0.6, inherit.aes = FALSE) +
  # flowlines in the gap
  geom_sf(data = gap_fl, color = "steelblue", linewidth = 0.5,
    inherit.aes = FALSE) +
  # HUC12 outlet
  geom_sf(data = hu12_pts[hu12_pts$comid == hu12_comid, ],
    shape = 4, color = "darkred", size = 1, stroke = 0.6, alpha = 0.5,
    inherit.aes = FALSE) +
  # gage point
  geom_sf(data = gage_pt, shape = 17, color = "black",
    fill = "white", size = 5, stroke = 1.4,
    inherit.aes = FALSE) +
  coord_sf(crs = target_crs, xlim = focus_xlim, ylim = focus_ylim,
    expand = FALSE) +
  labs(title = paste0("Davidson Creek (USGS streamgage 08110075)",
    " -- Gage and HU12 Outlet Positions")) +
  map_theme()
print(p_overview)

::: {.vignette-note} - The gage (triangle) and the HU12 pour point (x) sit on different catchments with several gap catchments (light blue) between them. - The dashed outline marks the catchment where the HU12 pour point falls; the solid outline marks the gage's catchment. - Both the gage and the HU12 outlet sit partway along their respective flowlines, so splitting is needed at both locations. :::

HU Outlet Detail

The HU12 split catchment is small enough that it is not visible in the overview. This view zooms to a 500 m square centered on the HU12 outlet.

hu12_outlet_pt <- hu12_pts[hu12_pts$comid == hu12_comid, ]
hu12_coords <- st_coordinates(hu12_outlet_pt)
half_side <- 250 # meters in EPSG:3857
zoom_xlim <- c(hu12_coords[1, "X"] - half_side,
  hu12_coords[1, "X"] + half_side)
zoom_ylim <- c(hu12_coords[1, "Y"] - half_side,
  hu12_coords[1, "Y"] + half_side)

p_huc_zoom <- ggplot() |>
  add_topo(hu12_outlet_pt) +
  # gap catchments
  geom_sf(data = extra_cat, fill = "lightblue", color = "steelblue",
    linewidth = 0.3, alpha = 0.5, inherit.aes = FALSE) +
  # HUC12 split catchment — full outline
  geom_sf(data = hu12_catch_full, fill = NA, color = "gray10",
    linewidth = 0.6, linetype = "dashed", inherit.aes = FALSE) +
  # HUC12 split portion
  geom_sf(data = hu12_catch_split, fill = "orange", color = "gray10",
    linewidth = 0.4, alpha = 0.5, inherit.aes = FALSE) +
  # flowlines in the gap
  geom_sf(data = gap_fl, color = "steelblue", linewidth = 0.5,
    inherit.aes = FALSE) +
  # HUC12 outlet
  geom_sf(data = hu12_outlet_pt,
    shape = 4, color = "darkred", size = 1, stroke = 0.6, alpha = 0.5,
    inherit.aes = FALSE) +
  coord_sf(crs = target_crs, xlim = zoom_xlim, ylim = zoom_ylim,
    expand = FALSE) +
  labs(title = paste0("Davidson Creek -- HU12 Outlet Detail",
    " (500 m view)")) +
  map_theme()
print(p_huc_zoom)

::: {.vignette-note} - The HU12 outlet (x) and its split catchment (orange fill, dashed outline) are visible at this scale. - The split catchment upstream of the HU outlet is excluded from the gap area because it overlaps with the HU12 drainage area. :::

Split Catchments: Upstream and Downstream Portions

Two catchments are split. At the gage, the downstream portion is removed because it does not contribute flow to the gage. At the HU12 outlet, the upstream portion is removed because it overlaps with the HU12 area already counted in the drainage area estimate.

# Build an sf with labeled polygons for a single legend
split_layers <- rbind(
  st_sf(
    role = "Upstream of gage (included)",
    geometry = st_geometry(osc_split)),
  st_sf(
    role = "Downstream of gage (excluded)",
    geometry = st_difference(
      st_geometry(osc_full), st_geometry(osc_split))),
  st_sf(
    role = "Upstream of HU outlet (excluded, overlaps HU)",
    geometry = st_geometry(hu12_catch_split)),
  st_sf(
    role = "Downstream of HU outlet (included as local area)",
    geometry = st_difference(
      st_geometry(hu12_catch_full), st_geometry(hu12_catch_split)))
)

split_colors <- c(
  "Upstream of gage (included)" = "#4DAF4A",
  "Downstream of gage (excluded)" = "#E41A1C",
  "Upstream of HU outlet (excluded, overlaps HU)" = "#FF7F00",
  "Downstream of HU outlet (included as local area)" = "#377EB8"
)

split_layers$role <- factor(split_layers$role,
  levels = names(split_colors))

p_split <- ggplot() |>
  add_topo(focus_geom) +
  # gap catchments as underlay
  geom_sf(data = extra_cat, fill = "gray80", color = "gray60",
    linewidth = 0.2, alpha = 0.3, inherit.aes = FALSE) +
  # split polygons with role-based fill
  geom_sf(data = split_layers,
    aes(fill = role), color = "gray20", linewidth = 0.4,
    alpha = 0.6, inherit.aes = FALSE) +
  scale_fill_manual(values = split_colors, name = NULL) +
  # flowlines
  geom_sf(data = gap_fl, color = "steelblue", linewidth = 0.5,
    inherit.aes = FALSE) +
  # HUC12 outlet
  geom_sf(data = hu12_pts[hu12_pts$comid == hu12_comid, ],
    shape = 4, color = "darkred", size = 1, stroke = 0.6, alpha = 0.5,
    inherit.aes = FALSE) +
  # gage
  geom_sf(data = gage_pt, shape = 17, color = "black",
    fill = "white", size = 5, stroke = 1.4,
    inherit.aes = FALSE) +
  coord_sf(crs = target_crs, xlim = focus_xlim, ylim = focus_ylim,
    expand = FALSE) +
  labs(title = paste0("Davidson Creek -- Split Catchment Roles"),
    subtitle = paste0(
      "Flowline measure at gage: ",
      round(da_dav$outlet_flowline_measure, 1),
      "; downstream removed: ",
      round(osc_full$dasqkm - osc_split$dasqkm, 2), " km\u00B2")) +
  map_theme() +
  guides(fill = guide_legend(ncol = 1))
print(p_split)

::: {.vignette-note} - Green: the portion of the gage's catchment upstream of the gage — this area is included in the drainage area estimate. - Red: the portion downstream of the gage — excluded because it does not contribute flow to the gage. - Orange: the portion of the HU12 outlet catchment upstream of the pour point — excluded because this area is already counted as part of the HU12 drainage area. - Blue: the portion downstream of the HU12 outlet — included as local area in the gap between the HU boundary and the gage. - The outlet_split_threshold_m parameter (default 100 m) controls whether the gage split is performed. If the gage is within the threshold distance of the catchment outlet, no split occurs. :::

options(oldoption)

Any scripts or data that you put into this service are public.

nhdplusTools documentation built on May 26, 2026, 5:06 p.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

nhdplusTools NHDPlus Tools

Drainage Area Estimation" In nhdplusTools: NHDPlus Tools

How get_drainage_area_estimates() Works

Algorithm Steps

Catchment Splitting

Running with Local Data

Web Services and Performance

Fetch and Cache Results

Return Structure

Basin Vignettes {.tabset}

French Broad River at Asheville

Drainage Area Boundaries

Stream Network

HU12 Types

Brazos River, West Texas

Drainage Area Boundaries

Stream Network

HU12 Types

James River / Cottonwood Lake

Drainage Area Boundaries

Stream Network

HU12 Types

Black Earth Creek

Drainage Area Boundaries

Stream Network

HU12 Types

Malheur Lake

Drainage Area Boundaries

Stream Network

HU12 Types

Purgatoire River near Las Animas

Drainage Area Boundaries

Stream Network

HU12 Types

Outlet Catchment Splitting

Overview: Gage and HU Outlet Positions

HU Outlet Detail

Split Catchments: Upstream and Downstream Portions

Try the nhdplusTools package in your browser

R Package Documentation

Browse R Packages

We want your feedback!

nhdplusTools
NHDPlus Tools

Drainage Area Estimation"
In nhdplusTools: NHDPlus Tools

How `get_drainage_area_estimates()` Works