diff --git a/DESCRIPTION b/DESCRIPTION
index fd6e0d4..a748291 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -24,6 +24,7 @@ Imports:
     censusapi,
     crosswalk,
     curl,
+    dataRetrieval (>= 2.7.25),
     sfarrow,
     dplyr,
     esri2sf (>= 0.1.1),
@@ -50,9 +51,11 @@ Imports:
     sf,
     stats,
     stringr,
+    terra,
     tibble,
     tidycensus,
     tidyr,
+    tidylog,
     tidytable,
     tigris,
     utils,
diff --git a/NAMESPACE b/NAMESPACE
index ad7996c..130b59f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -15,6 +15,7 @@ export(get_fema_disaster_declarations)
 export(get_fema_floodplain)
 export(get_government_finances)
 export(get_hazard_mitigation_assistance)
+export(get_hrrr_smoke)
 export(get_hud_api_key)
 export(get_ihp_registrations)
 export(get_lodes)
@@ -33,6 +34,7 @@ export(get_sheldus)
 export(get_spatial_extent_census)
 export(get_structures)
 export(get_system_username)
+export(get_usgs_gage)
 export(get_wildfire_burn_zones)
 export(inflation_adjust)
 export(list_openfema_endpoints)
diff --git a/R/get_hrrr_smoke.R b/R/get_hrrr_smoke.R
new file mode 100644
index 0000000..8c7b386
--- /dev/null
+++ b/R/get_hrrr_smoke.R
@@ -0,0 +1,208 @@
+#' Get hourly wildfire smoke concentrations from the HRRR-Smoke model
+#'
+#' @description
+#' Retrieves hourly near-surface wildfire smoke concentrations (micrograms per
+#' cubic meter) from NOAA's High-Resolution Rapid Refresh (HRRR) model, cropped
+#' to an area of interest, and returns them as a single multi-layer raster
+#' (one layer per hour).
+#'
+#' @details
+#' HRRR is NOAA's 3-kilometer, hourly-updating weather model for the
+#' conterminous United States. Since late 2020 it has carried smoke as a
+#' modeled quantity, driven by satellite detections of active fires. This
+#' function returns the "analysis" field for each requested hour -- the
+#' model's real-time estimate for the hour it was issued. Data
+#' are downloaded on demand from NOAA's free public archive. 
+#' A two-week window at hourly resolution takes roughly a few minutes.
+#'
+#' Two smoke quantities are available via `variable`:
+#' \describe{
+#'   \item{`"surface"`}{Smoke mass density 8 meters above ground, in
+#'     micrograms per cubic meter (ug/m^3). This approximates what people at
+#'     ground level are breathing and is directly comparable to PM2.5 air
+#'     quality readings, which use the same unit. For reference, the EPA's
+#'     24-hour PM2.5 standard is 35 ug/m^3.}
+#'   \item{`"column"`}{Vertically integrated smoke -- all smoke in the
+#'     atmospheric column above each cell -- in milligrams per square meter
+#'     (mg/m^2). This corresponds to what satellites see and includes
+#'     high-altitude smoke that may never reach the ground.}
+#' }
+#'
+#' Because HRRR covers only the conterminous United States, Alaska, Hawaii,
+#' and the territories are unsupported. Note also that these are model
+#' estimates, not directly-measure smoke concentration observations.
+#'
+#' @param geometries An `sf`-formatted dataframe (or an `sfc` geometry column)
+#'   defining the area of interest, in any defined coordinate reference
+#'   system. The returned raster is cropped to this area's bounding box.
+#' @param start_date The first day to retrieve, as a `Date` or a
+#'   "YYYY-MM-DD" string.
+#' @param end_date The last day to retrieve (inclusive), as a `Date` or a
+#'   "YYYY-MM-DD" string. Defaults to `start_date`. HRRR-Smoke is archived
+#'   from 2021 onward; the most recent hours may not yet be posted.
+#' @param variable Which smoke quantity to retrieve: `"surface"` (default;
+#'   near-surface concentration) or `"column"` (vertically integrated smoke).
+#'   See Details.
+#' @param hours Which hours of each day (UTC, 0-23) to retrieve. Defaults to
+#'   all 24; for lighter-temporal-weight coverage, pass e.g. `seq(0, 21, by = 3)`.
+#'
+#' @return A `terra::SpatRaster` with one layer per successfully retrieved
+#'   hour, cropped to the bounding box of `geometries` (buffered by one
+#'   3-kilometer cell). Hours missing from the archive are dropped with a
+#'   single summary warning. The raster's components:
+#' \describe{
+#'   \item{cell values}{Numeric. The smoke quantity selected by `variable`:
+#'     near-surface smoke concentration in micrograms per cubic meter
+#'     (ug/m^3) when `variable = "surface"`, or vertically integrated
+#'     column smoke in milligrams per square meter (mg/m^2) when
+#'     `variable = "column"`.}
+#'   \item{layers}{One layer per hour, in chronological order. Convert to a
+#'     one-row-per-cell-per-hour tibble with
+#'     `terra::as.data.frame(x, xy = TRUE, wide = FALSE)`.}
+#'   \item{layer names}{Character. The layer's timestamp in UTC, formatted
+#'     "YYYY-MM-DD HH:00" (e.g. "2025-08-01 12:00").}
+#'   \item{time}{POSIXct. The same UTC timestamps, retrievable with
+#'     `terra::time()`; used directly by `tidyterra` and
+#'     `terra::animate()`.}
+#'   \item{coordinate reference system}{The HRRR model's native projection
+#'     (Lambert conformal conic), with 3-kilometer cells. Reproject with
+#'     `terra::project()`, or transform vector layers to it with
+#'     `sf::st_transform(x, sf::st_crs(raster))` before mapping.}
+#' }
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' county = tigris::counties(state = "CA", cb = TRUE) %>%
+#'   dplyr::filter(NAME == "Butte")
+#'
+#' smoke = get_hrrr_smoke(
+#'   geometries = county,
+#'   start_date = "2025-07-20",
+#'   end_date = "2025-08-03")
+#'
+#' # quick look at one hour, and a simple animation across all hours
+#' terra::plot(smoke[[1]])
+#' terra::animate(smoke, pause = 0.1)
+#' }
+get_hrrr_smoke = function(
+    geometries,
+    start_date,
+    end_date = start_date,
+    variable = c("surface", "column"),
+    hours = 0:23) {
+
+  variable = match.arg(variable)
+
+  start_date = as.Date(start_date)
+  end_date = as.Date(end_date)
+  if (is.na(start_date) || is.na(end_date)) {
+    stop("`start_date` and `end_date` must be Dates or 'YYYY-MM-DD' strings.") }
+  if (end_date < start_date) {
+    stop("`end_date` must not be earlier than `start_date`.") }
+  ## HRRR added smoke fields with the model's version 4 upgrade in December 2020;
+  ## the AWS archive holds them reliably from 2021 onward
+  if (start_date < as.Date("2021-01-01")) {
+    stop("HRRR-Smoke fields are available in the archive from 2021-01-01 onward.") }
+  if (!all(hours %in% 0:23)) {
+    stop("`hours` must contain only integers between 0 and 23.") }
+
+  if (inherits(geometries, "sfc")) { geometries = sf::st_as_sf(geometries) }
+  if (!inherits(geometries, "sf")) {
+    stop("`geometries` must be a simple features (sf) object.") }
+  if (is.na(sf::st_crs(geometries))) {
+    stop("`geometries` must have a defined coordinate reference system (CRS).") }
+
+  ## the string that identifies the smoke field on a line of the .idx sidecar
+  ## file, e.g. "MASSDEN:8 m above ground" (GRIB shorthand for smoke mass density)
+  idx_pattern = switch(
+    variable,
+    surface = "MASSDEN:8 m above ground",
+    column = "COLMD:entire atmosphere")
+
+  base_url = "https://noaa-hrrr-bdp-pds.s3.amazonaws.com"
+
+  ## one row per requested hour: the archive path of that hour's analysis file
+  requests = tidyr::expand_grid(
+      date = seq(start_date, end_date, by = "day"),
+      hour = sort(unique(as.integer(hours)))) %>%
+    dplyr::mutate(
+      timestamp = as.POSIXct(
+        stringr::str_c(date, " ", hour, ":00"), tz = "UTC"),
+      grib_url = stringr::str_c(
+        base_url, "/hrrr.", format(date, "%Y%m%d"), "/conus/hrrr.t",
+        sprintf("%02d", hour), "z.wrfsfcf00.grib2"))
+
+  ## fetch one hour's smoke field: read the .idx to find the field's byte range,
+  ## download only those bytes, read as a raster. Returns NULL if the hour is
+  ## not (yet) in the archive.
+  fetch_hour = function(grib_url, timestamp) {
+    idx_lines = tryCatch(
+      readLines(stringr::str_c(grib_url, ".idx"), warn = FALSE),
+      error = function(e) NULL)
+    if (is.null(idx_lines)) { return(NULL) }
+
+    ## .idx lines look like "37:24296434:d=2025072012:MASSDEN:8 m above ground:anl:"
+    ## -- field 2 is the field's starting byte; the next line's start is its end
+    line_number = stringr::str_which(idx_lines, stringr::fixed(idx_pattern))
+    if (length(line_number) != 1) { return(NULL) }
+
+    byte_starts = as.numeric(stringr::str_split_i(idx_lines, ":", 2))
+    range_start = byte_starts[line_number]
+    ## for the last field in the file there is no next line; curl accepts an
+    ## open-ended range ("start-"), which reads through the end of the file
+    range_end = dplyr::if_else(
+      line_number < length(idx_lines),
+      as.character(byte_starts[line_number + 1] - 1),
+      "")
+
+    grib_file = tempfile(fileext = ".grib2")
+    fetch = tryCatch({
+        curl::curl_download(
+          grib_url,
+          grib_file,
+          handle = curl::new_handle(
+            range = stringr::str_c(range_start, "-", range_end)))
+        terra::rast(grib_file) },
+      error = function(e) NULL)
+    if (is.null(fetch)) { return(NULL) }
+
+    names(fetch) = format(timestamp, "%Y-%m-%d %H:00")
+    terra::time(fetch) = timestamp
+    fetch
+  }
+
+  hourly_rasters = purrr::map2(
+    requests$grib_url,
+    requests$timestamp,
+    fetch_hour) %>%
+    purrr::compact()
+
+  if (length(hourly_rasters) == 0) {
+    stop(
+      "No HRRR-Smoke fields could be retrieved for the requested window. ",
+      "Check that the dates are not in the future and that you are online.") }
+
+  missing_count = nrow(requests) - length(hourly_rasters)
+  if (missing_count > 0) {
+    warning(
+      missing_count, " of ", nrow(requests),
+      " requested hours were not available in the HRRR archive and were dropped.") }
+
+  ## GRIB files store smoke in kilograms (per cubic meter for "surface", per
+  ## square meter for "column"); convert to the micrograms / milligrams
+  ## documented above, which match how air quality figures are usually reported
+  unit_factor = switch(variable, surface = 1e9, column = 1e6)
+  smoke_stack = terra::rast(hourly_rasters) * unit_factor
+
+  ## crop to the area of interest in the model's native projection. The bounding
+  ## box is buffered by one cell (3 km) so the area's edges are fully covered.
+  area_of_interest = geometries %>%
+    sf::st_transform(sf::st_crs(smoke_stack)) %>%
+    sf::st_bbox() %>%
+    sf::st_as_sfc() %>%
+    sf::st_buffer(3000) %>%
+    terra::vect()
+
+  terra::crop(smoke_stack, area_of_interest)
+}
diff --git a/R/get_preliminary_damage_assessments.R b/R/get_preliminary_damage_assessments.R
index 1e76e42..a12b256 100644
--- a/R/get_preliminary_damage_assessments.R
+++ b/R/get_preliminary_damage_assessments.R
@@ -1286,7 +1286,10 @@ correct_duplicate_disaster_numbers = function(pda_df) {
   corrected = pda_df %>%
     ## coerce so the if_else() below is type-stable regardless of how a cached CSV
     ## parsed the column (readr may guess double; the extracted value is character)
-    dplyr::mutate(disaster_number = as.character(disaster_number)) %>%
+    dplyr::mutate(
+      dplyr::across(
+        dplyr::any_of(c("disaster_number", "disaster_number_filename")),
+        as.character)) %>%
     dplyr::add_count(disaster_number, name = "disaster_number_count") %>%
     dplyr::mutate(
       disaster_number_from_text = stringr::str_extract(text, "FEMA-[0-9]{4}") %>%
diff --git a/R/get_usgs_gage.R b/R/get_usgs_gage.R
new file mode 100644
index 0000000..1430a75
--- /dev/null
+++ b/R/get_usgs_gage.R
@@ -0,0 +1,318 @@
+#' @importFrom magrittr %>%
+
+#' @title Acquire daily stream-gage readings from USGS gages
+#'
+#' @description Pulls daily stream-gage readings, over each gage's full period of
+#'   record by default, for USGS gages in one or more counties via the
+#'   dataRetrieval package (USGS Water Data OGC APIs; these replace the
+#'   now-decommissioning NWIS web services). Two statistics are supported:
+#'   \itemize{
+#'     \item "daily_mean": the published daily-mean series. Available for a
+#'       century-plus at many gages.
+#'     \item "daily_max": the maximum reading each day, computed here from the
+#'       continuous (15-minute) record, because USGS publishes no daily-maximum
+#'       gage-height series. Continuous records only begin in the mid-1990s at
+#'       the earliest. These pulls are slow (a minute or more per long-record
+#'       gage), so each gage's aggregated result is cached to its own parquet
+#'       file and the pull resumes wherever it left off.
+#'   }
+#'
+#' @param counties Character vector of five-digit county FIPS codes (e.g.,
+#'   "54097" for Upshur County, WV). Every stream gage in these counties with
+#'   data for the requested measure and statistic is pulled.
+#' @param measure One of "height" (gage height in feet, the default; USGS
+#'   parameter code 00065) or "discharge" (streamflow in cubic feet per second;
+#'   USGS parameter code 00060).
+#' @param statistic One of "daily_max" (the default; computed from continuous
+#'   readings) or "daily_mean" (the published daily-value series). See the
+#'   description for the record-length and runtime trade-offs.
+#' @param start_date,end_date Character "YYYY-MM-DD" bounds on the readings. The
+#'   defaults ("" for both) request each gage's full period of record; either
+#'   bound may be supplied alone.
+#' @param refresh_cache When TRUE, ignore cached parquet files (including the
+#'   per-site continuous-record caches) and pull fresh data. Defaults to FALSE.
+#'   Note that without a refresh, previously cached gages are frozen at the time
+#'   they were pulled.
+#' @param cache_dir Directory for cached parquet files. Defaults to a
+#'   session-specific temporary directory; supply a persistent directory (e.g.,
+#'   `tools::R_user_dir("climateapi", which = "cache")`) to keep the cache
+#'   across sessions, which is strongly recommended for
+#'   `statistic = "daily_max"` pulls so they can resume.
+#'
+#' @details Site metadata (name, county, coordinates, drainage area) is attached
+#'   to every reading. Data are from the USGS Water Data APIs; see
+#'   \url{https://api.waterdata.usgs.gov/}. USGS asks that heavy users register
+#'   a free API key (see `dataRetrieval::setAccess()` documentation); unkeyed
+#'   access is rate-limited but sufficient for modest pulls.
+#'
+#' @returns A tibble with one row per gage-day. Columns include:
+#'   \describe{
+#'     \item{site_number}{USGS site number.}
+#'     \item{gage_name}{USGS station name.}
+#'     \item{county_geoid}{Five-digit county FIPS code.}
+#'     \item{county_name}{County name.}
+#'     \item{state_abbreviation}{Two-letter state abbreviation.}
+#'     \item{latitude, longitude}{Gage coordinates (decimal degrees).}
+#'     \item{drainage_area_sqmi}{Upstream drainage area, in square miles.}
+#'     \item{date}{Calendar day. For `statistic = "daily_max"`, the day is
+#'       defined in the gage's local (Eastern) time zone.}
+#'     \item{value}{The reading, in feet ("height") or cubic feet per second
+#'       ("discharge").}
+#'     \item{approval_status}{"approved" or "provisional". For
+#'       `statistic = "daily_max"`, a day is "provisional" if any reading that
+#'       day is provisional.}
+#'   }
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' ## daily maximum gage heights in Upshur and Lewis Counties, WV
+#' get_usgs_gage(
+#'   counties = c("54097", "54041"),
+#'   cache_dir = tools::R_user_dir("climateapi", which = "cache"))
+#'
+#' get_usgs_gage(
+#'   counties = "54063",
+#'   measure = "discharge",
+#'   statistic = "daily_mean")
+#' }
+get_usgs_gage = function(
+    counties,
+    measure = c("height", "discharge"),
+    statistic = c("daily_max", "daily_mean"),
+    start_date = "",
+    end_date = "",
+    refresh_cache = FALSE,
+    cache_dir = file.path(tempdir(), "usgs-gage-history")) {
+
+  measure = match.arg(measure)
+  statistic = match.arg(statistic)
+  parameter_code = dplyr::if_else(measure == "height", "00065", "00060")
+
+  if (missing(counties) || !is.character(counties) || length(counties) == 0 ||
+      any(!stringr::str_detect(counties, "^[0-9]{5}$"))) {
+    stop("`counties` must be a character vector of five-digit county FIPS codes (e.g., '54097').") }
+
+  purrr::walk2(
+    list(start_date, end_date), c("start_date", "end_date"),
+    function(date_value, date_name) {
+      if (!is.character(date_value) || length(date_value) != 1 ||
+          (date_value != "" && !stringr::str_detect(date_value, "^\\d{4}-\\d{2}-\\d{2}$"))) {
+        stop("`", date_name, "` must be '' or a 'YYYY-MM-DD' string.") } })
+
+  dir.create(cache_dir, showWarnings = FALSE, recursive = TRUE)
+  cache_path = file.path(
+    cache_dir,
+    stringr::str_c(
+      "gage_history_",
+      stringr::str_c(sort(unique(counties)), collapse = "-"), "_",
+      measure, "_", statistic,
+      dplyr::if_else(start_date == "", "", stringr::str_c("_from", start_date)),
+      dplyr::if_else(end_date == "", "", stringr::str_c("_to", end_date)),
+      ".parquet"))
+
+  if (file.exists(cache_path) && !refresh_cache) {
+    message("Reading cached file: ", basename(cache_path))
+    return(arrow::read_parquet(cache_path)) }
+
+  ## the monitoring-locations endpoint filters by two-digit state and
+  ## three-digit county FIPS codes, so the five-digit codes are split and the
+  ## endpoint queried once per state. Filtering to site_type_code "ST"
+  ## restricts results to stream gages (the endpoint also serves wells, lakes,
+  ## etc.). County FIPS codes, coordinates, drainage area, and station names
+  ## all come from this endpoint; drainage area lets readers normalize
+  ## discharge across basins of different sizes
+  county_geoids = sort(unique(counties))
+  site_metadata = county_geoids %>%
+    split(stringr::str_sub(., 1, 2)) %>%
+    purrr::imap(
+      ~ dataRetrieval::read_waterdata_monitoring_location(
+        state_code = .y,
+        county_code = stringr::str_sub(.x, 3, 5),
+        site_type_code = "ST")) %>%
+    purrr::list_rbind() %>%
+    dplyr::mutate(
+      latitude = sf::st_coordinates(geometry)[, "Y"],
+      longitude = sf::st_coordinates(geometry)[, "X"]) %>%
+    sf::st_drop_geometry() %>%
+    tibble::as_tibble() %>%
+    dplyr::transmute(
+      monitoring_location_id,
+      site_number = monitoring_location_number,
+      gage_name = monitoring_location_name,
+      county_geoid = stringr::str_c(
+        stringr::str_pad(state_code, 2, pad = "0"),
+        stringr::str_pad(county_code, 3, pad = "0")),
+      drainage_area_sqmi = drainage_area,
+      latitude,
+      longitude) %>%
+    ## a state-plus-county query can only match the requested counties, but the
+    ## endpoint treats each filter as an independent OR-list, so a multi-state
+    ## request (e.g. "54097" and "39041") could otherwise return unrequested
+    ## state-county combinations
+    dplyr::filter(county_geoid %in% county_geoids)
+
+  ## of the counties' stream sites, keep those with a series of the requested
+  ## measure and statistic. daily_mean uses the published daily "Mean" series;
+  ## daily_max must be computed from the continuous ("Instantaneous") record.
+  ## Each series' period-of-record start is retained: the continuous endpoint
+  ## needs explicit time bounds (see below)
+  site_inventory = dataRetrieval::read_waterdata_ts_meta(
+      monitoring_location_id = site_metadata$monitoring_location_id,
+      parameter_code = parameter_code,
+      computation_identifier = dplyr::if_else(
+        statistic == "daily_mean", "Mean", "Instantaneous"),
+      computation_period_identifier = dplyr::if_else(
+        statistic == "daily_mean", "Daily", "Points"),
+      skipGeometry = TRUE) %>%
+    tibble::as_tibble() %>%
+    ## a site can have several series (e.g. sublocations); keep the earliest
+    dplyr::summarize(
+      record_begin_date = as.Date(min(begin, na.rm = TRUE)),
+      .by = monitoring_location_id)
+
+  site_metadata = site_metadata %>%
+    dplyr::filter(monitoring_location_id %in% site_inventory$monitoring_location_id)
+
+  if (nrow(site_metadata) == 0) {
+    stop(
+      "No USGS stream gages with ", statistic, " ", measure,
+      " data found in counties: ", stringr::str_c(county_geoids, collapse = ", ")) }
+
+  county_names = tigris::fips_codes %>%
+    tibble::as_tibble() %>%
+    dplyr::transmute(
+      county_geoid = stringr::str_c(state_code, county_code),
+      county_name = stringr::str_remove(county, " County$"),
+      state_abbreviation = state)
+
+  ## the OGC APIs express time bounds as a start/end interval; ".." leaves an
+  ## end open, and NA requests each gage's full period of record
+  time_bounds = if (start_date == "" && end_date == "") {
+    NA_character_
+  } else {
+    c(
+      dplyr::if_else(start_date == "", "..", start_date),
+      dplyr::if_else(end_date == "", "..", end_date)) }
+
+  daily_readings = if (statistic == "daily_mean") {
+    ## published daily-mean values; read_waterdata_daily() batches large site
+    ## lists into multiple requests internally
+    message(
+      "Pulling daily means for ", nrow(site_metadata),
+      " gages; this can take several minutes.")
+
+    dataRetrieval::read_waterdata_daily(
+      monitoring_location_id = site_metadata$monitoring_location_id,
+      parameter_code = parameter_code,
+      statistic_id = "00003", ## daily mean
+      time = time_bounds,
+      skipGeometry = TRUE) %>%
+      tibble::as_tibble() %>%
+      dplyr::transmute(
+        monitoring_location_id,
+        date = time,
+        value,
+        approval_status = dplyr::if_else(
+          approval_status == "Approved", "approved", "provisional"))
+  } else {
+    ## daily maxima computed from the continuous (15-minute) record, one site
+    ## at a time. Each site's full continuous record is a multi-minute
+    ## download, so the aggregated daily maxima are cached per site and the pull
+    ## resumes wherever it left off (which also makes it safe to run several
+    ## worker processes over disjoint site lists against the same cache)
+    site_cache_dir = file.path(
+      cache_dir, stringr::str_c("by-site-", measure, "-daily-max"))
+    dir.create(site_cache_dir, showWarnings = FALSE, recursive = TRUE)
+    site_cache_paths = file.path(
+      site_cache_dir, stringr::str_c(site_metadata$site_number, ".parquet"))
+    uncached_ids = site_metadata$monitoring_location_id[
+      refresh_cache | !file.exists(site_cache_paths)]
+
+    if (length(uncached_ids) > 0) {
+      message(
+        "Pulling full continuous records for ", length(uncached_ids),
+        " gages (", nrow(site_metadata) - length(uncached_ids),
+        " already cached). Expect a minute or more per long-record gage.") }
+
+    ## unlike the daily endpoint, the continuous endpoint silently returns only
+    ## the most recent year when no time interval is given, and rejects
+    ## intervals much beyond three years with an HTTP 400 -- so each site's
+    ## record is requested in two-year chunks anchored at its period-of-record
+    ## start
+    site_record_begins = site_inventory %>%
+      dplyr::filter(monitoring_location_id %in% uncached_ids)
+
+    for (site_id in uncached_ids) {
+      site_start = dplyr::if_else(
+        start_date == "",
+        site_record_begins$record_begin_date[
+          site_record_begins$monitoring_location_id == site_id],
+        as.Date(start_date))
+      site_end = dplyr::if_else(
+        end_date == "", Sys.Date() + 1, as.Date(end_date))
+
+      chunk_starts = seq(site_start, site_end, by = "2 years")
+      chunk_ends = c(utils::tail(chunk_starts, -1), site_end)
+
+      site_daily_maxima = tryCatch(
+        purrr::map2(
+          chunk_starts, chunk_ends,
+          ~ dataRetrieval::read_waterdata_continuous(
+            monitoring_location_id = site_id,
+            parameter_code = parameter_code,
+            time = c(as.character(.x), as.character(.y)))) %>%
+          purrr::list_rbind() %>%
+          tibble::as_tibble() %>%
+          dplyr::filter(!is.na(value)) %>%
+          ## readings arrive in UTC; the calendar day is defined in the gage's
+          ## local (Eastern) time zone. A day is provisional if any reading
+          ## that day is provisional
+          dplyr::mutate(
+            date = as.Date(lubridate::with_tz(time, "America/New_York"))) %>%
+          dplyr::summarize(
+            value = max(value),
+            approval_status = dplyr::if_else(
+              all(approval_status == "Approved"), "approved", "provisional"),
+            .by = c(monitoring_location_id, date)),
+        error = function(e) {
+          warning(
+            "Continuous-record request failed for site ",
+            stringr::str_remove(site_id, "^USGS-"), ": ", conditionMessage(e))
+          NULL })
+      if (!is.null(site_daily_maxima)) {
+        arrow::write_parquet(
+          site_daily_maxima,
+          file.path(
+            site_cache_dir,
+            stringr::str_c(stringr::str_remove(site_id, "^USGS-"), ".parquet"))) } }
+
+    available_cache_paths = file.path(
+      site_cache_dir, stringr::str_c(site_metadata$site_number, ".parquet"))
+    available_cache_paths[file.exists(available_cache_paths)] %>%
+      purrr::map(arrow::read_parquet) %>%
+      purrr::list_rbind()
+  }
+
+  gage_history = daily_readings %>%
+    dplyr::filter(!is.na(value)) %>%
+    tidylog::left_join(
+      site_metadata, by = "monitoring_location_id", relationship = "many-to-one") %>%
+    tidylog::left_join(county_names, by = "county_geoid", relationship = "many-to-one") %>%
+    dplyr::select(
+      site_number, gage_name, county_geoid, county_name, state_abbreviation,
+      latitude, longitude, drainage_area_sqmi, date, value, approval_status)
+
+  arrow::write_parquet(gage_history, cache_path)
+
+  return(gage_history)
+}
+
+utils::globalVariables(c(
+  "monitoring_location_id", "monitoring_location_number",
+  "monitoring_location_name", "site_type_code", "state_code", "county_code",
+  "drainage_area", "county", "state", "geometry", "time",
+  "begin", "record_begin_date", ".",
+  "site_number", "gage_name", "county_geoid", "county_name",
+  "state_abbreviation", "latitude", "longitude", "drainage_area_sqmi",
+  "value", "approval_status", "date"))
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 2e2d1c9..1553fa3 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -5,8 +5,11 @@ reference:
 - title: Disaster events
   contents: 
   - get_fema_disaster_declarations
+  - get_fema_floodplain
+  - get_usgs_gage
   - get_current_fire_perimeters
   - get_wildfire_burn_zones
+  - get_hrrr_smoke
 - title: Disaster-related damages and funding
   contents: 
   - get_national_risk_index
@@ -53,4 +56,5 @@ reference:
   - get_naics_codes
   - get_hud_api_key
   - register_hud_api_key
+  - scrape_pda_pdfs
 
diff --git a/man/get_hrrr_smoke.Rd b/man/get_hrrr_smoke.Rd
new file mode 100644
index 0000000..a68d567
--- /dev/null
+++ b/man/get_hrrr_smoke.Rd
@@ -0,0 +1,105 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/get_hrrr_smoke.R
+\name{get_hrrr_smoke}
+\alias{get_hrrr_smoke}
+\title{Get hourly wildfire smoke concentrations from the HRRR-Smoke model}
+\usage{
+get_hrrr_smoke(
+  geometries,
+  start_date,
+  end_date = start_date,
+  variable = c("surface", "column"),
+  hours = 0:23
+)
+}
+\arguments{
+\item{geometries}{An \code{sf}-formatted dataframe (or an \code{sfc} geometry column)
+defining the area of interest, in any defined coordinate reference
+system. The returned raster is cropped to this area's bounding box.}
+
+\item{start_date}{The first day to retrieve, as a \code{Date} or a
+"YYYY-MM-DD" string.}
+
+\item{end_date}{The last day to retrieve (inclusive), as a \code{Date} or a
+"YYYY-MM-DD" string. Defaults to \code{start_date}. HRRR-Smoke is archived
+from 2021 onward; the most recent hours may not yet be posted.}
+
+\item{variable}{Which smoke quantity to retrieve: \code{"surface"} (default;
+near-surface concentration) or \code{"column"} (vertically integrated smoke).
+See Details.}
+
+\item{hours}{Which hours of each day (UTC, 0-23) to retrieve. Defaults to
+all 24; for lighter-temporal-weight coverage, pass e.g. \code{seq(0, 21, by = 3)}.}
+}
+\value{
+A \code{terra::SpatRaster} with one layer per successfully retrieved
+hour, cropped to the bounding box of \code{geometries} (buffered by one
+3-kilometer cell). Hours missing from the archive are dropped with a
+single summary warning. The raster's components:
+\describe{
+\item{cell values}{Numeric. The smoke quantity selected by \code{variable}:
+near-surface smoke concentration in micrograms per cubic meter
+(ug/m^3) when \code{variable = "surface"}, or vertically integrated
+column smoke in milligrams per square meter (mg/m^2) when
+\code{variable = "column"}.}
+\item{layers}{One layer per hour, in chronological order. Convert to a
+one-row-per-cell-per-hour tibble with
+\code{terra::as.data.frame(x, xy = TRUE, wide = FALSE)}.}
+\item{layer names}{Character. The layer's timestamp in UTC, formatted
+"YYYY-MM-DD HH:00" (e.g. "2025-08-01 12:00").}
+\item{time}{POSIXct. The same UTC timestamps, retrievable with
+\code{terra::time()}; used directly by \code{tidyterra} and
+\code{terra::animate()}.}
+\item{coordinate reference system}{The HRRR model's native projection
+(Lambert conformal conic), with 3-kilometer cells. Reproject with
+\code{terra::project()}, or transform vector layers to it with
+\code{sf::st_transform(x, sf::st_crs(raster))} before mapping.}
+}
+}
+\description{
+Retrieves hourly near-surface wildfire smoke concentrations (micrograms per
+cubic meter) from NOAA's High-Resolution Rapid Refresh (HRRR) model, cropped
+to an area of interest, and returns them as a single multi-layer raster
+(one layer per hour).
+}
+\details{
+HRRR is NOAA's 3-kilometer, hourly-updating weather model for the
+conterminous United States. Since late 2020 it has carried smoke as a
+modeled quantity, driven by satellite detections of active fires. This
+function returns the "analysis" field for each requested hour -- the
+model's real-time estimate for the hour it was issued. Data
+are downloaded on demand from NOAA's free public archive.
+A two-week window at hourly resolution takes roughly a few minutes.
+
+Two smoke quantities are available via \code{variable}:
+\describe{
+\item{\code{"surface"}}{Smoke mass density 8 meters above ground, in
+micrograms per cubic meter (ug/m^3). This approximates what people at
+ground level are breathing and is directly comparable to PM2.5 air
+quality readings, which use the same unit. For reference, the EPA's
+24-hour PM2.5 standard is 35 ug/m^3.}
+\item{\code{"column"}}{Vertically integrated smoke -- all smoke in the
+atmospheric column above each cell -- in milligrams per square meter
+(mg/m^2). This corresponds to what satellites see and includes
+high-altitude smoke that may never reach the ground.}
+}
+
+Because HRRR covers only the conterminous United States, Alaska, Hawaii,
+and the territories are unsupported. Note also that these are model
+estimates, not directly-measure smoke concentration observations.
+}
+\examples{
+\dontrun{
+county = tigris::counties(state = "CA", cb = TRUE) \%>\%
+  dplyr::filter(NAME == "Butte")
+
+smoke = get_hrrr_smoke(
+  geometries = county,
+  start_date = "2025-07-20",
+  end_date = "2025-08-03")
+
+# quick look at one hour, and a simple animation across all hours
+terra::plot(smoke[[1]])
+terra::animate(smoke, pause = 0.1)
+}
+}
diff --git a/man/get_preliminary_damage_assessments.Rd b/man/get_preliminary_damage_assessments.Rd
index ae648e9..b6cccf0 100644
--- a/man/get_preliminary_damage_assessments.Rd
+++ b/man/get_preliminary_damage_assessments.Rd
@@ -28,7 +28,13 @@ A dataframe of preliminary damage assessment reports. Columns include:
 \describe{
 \item{path}{The local file path to the source PDA PDF.}
 \item{disaster_number}{FEMA disaster number.}
-\item{event_type}{Type of decision: "approved", "denial", "appeal_approved", or "appeal_denial".}
+\item{event_type}{Type of decision: "approved", "denial", "appeal_approved", or
+"appeal_denial". The denial classes are read from FEMA's filename convention and the
+report title. "appeal_approved" is read from the report body instead, because an
+approved appeal is titled and named exactly like a first-instance approval and carries
+a disaster number; what identifies it is a narrative of a denied request that was
+subsequently appealed. Both halves of that narrative are required, so ordinary
+approvals that merely describe the appeals process are not misclassified.}
 \item{event_title}{Title/description of the disaster event.}
 \item{event_date_determined}{Date the PDA determination was made.}
 \item{event_native_flag}{1 if tribal request, 0 otherwise.}
@@ -37,12 +43,16 @@ A dataframe of preliminary damage assessment reports. Columns include:
 \item{pa_primary_impact}{The primary type of impact described for Public Assistance purposes.}
 \item{pa_cost_estimate_total}{Estimated total Public Assistance cost.}
 \item{pa_per_capita_impact_statewide}{Statewide (or territory/commonwealth) per capita impact amount.}
-\item{pa_per_capita_impact_indicator_statewide}{Numeric ratio of the statewide per capita impact
-to the applicable threshold -- a decimal ratio (e.g. 1.5, 1.89), not a "Met"/"Not Met"
-categorical indicator despite the field's FEMA-assigned name.}
+\item{pa_per_capita_impact_indicator_statewide}{FEMA's statutory statewide per capita
+\emph{threshold} in dollars for the relevant year (observed range 1.24--1.94), not a ratio
+and not a "Met"/"Not Met" categorical despite the field's FEMA-assigned name. Compare
+it against \code{pa_per_capita_impact_statewide}, which is the estimated per capita impact
+in the same units; the ratio of the two is what indicates whether the threshold was met.}
 \item{pa_per_capita_impact_countywide}{Raw text of countywide per capita impact ratios (may list
 multiple values across affected counties for a multi-county event).}
-\item{pa_per_capita_impact_indicator_countywide}{Truncated text of the countywide per capita impact indicator.}
+\item{pa_per_capita_impact_indicator_countywide}{FEMA's statutory countywide per capita
+threshold in dollars (observed range 3.11--4.60), on the same basis as the statewide
+indicator above.}
 \item{pa_per_capita_impact_countywide_max}{Maximum countywide per capita impact ratio parsed
 from \code{pa_per_capita_impact_countywide}.}
 \item{pa_per_capita_impact_countywide_min}{Minimum countywide per capita impact ratio parsed
diff --git a/man/get_usgs_gage.Rd b/man/get_usgs_gage.Rd
new file mode 100644
index 0000000..3b46f12
--- /dev/null
+++ b/man/get_usgs_gage.Rd
@@ -0,0 +1,99 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/get_usgs_gage.R
+\name{get_usgs_gage}
+\alias{get_usgs_gage}
+\title{Acquire daily stream-gage readings from USGS gages}
+\usage{
+get_usgs_gage(
+  counties,
+  measure = c("height", "discharge"),
+  statistic = c("daily_max", "daily_mean"),
+  start_date = "",
+  end_date = "",
+  refresh_cache = FALSE,
+  cache_dir = file.path(tempdir(), "usgs-gage-history")
+)
+}
+\arguments{
+\item{counties}{Character vector of five-digit county FIPS codes (e.g.,
+"54097" for Upshur County, WV). Every stream gage in these counties with
+data for the requested measure and statistic is pulled.}
+
+\item{measure}{One of "height" (gage height in feet, the default; USGS
+parameter code 00065) or "discharge" (streamflow in cubic feet per second;
+USGS parameter code 00060).}
+
+\item{statistic}{One of "daily_max" (the default; computed from continuous
+readings) or "daily_mean" (the published daily-value series). See the
+description for the record-length and runtime trade-offs.}
+
+\item{start_date, end_date}{Character "YYYY-MM-DD" bounds on the readings. The
+defaults ("" for both) request each gage's full period of record; either
+bound may be supplied alone.}
+
+\item{refresh_cache}{When TRUE, ignore cached parquet files (including the
+per-site continuous-record caches) and pull fresh data. Defaults to FALSE.
+Note that without a refresh, previously cached gages are frozen at the time
+they were pulled.}
+
+\item{cache_dir}{Directory for cached parquet files. Defaults to a
+session-specific temporary directory; supply a persistent directory (e.g.,
+\code{tools::R_user_dir("climateapi", which = "cache")}) to keep the cache
+across sessions, which is strongly recommended for
+\code{statistic = "daily_max"} pulls so they can resume.}
+}
+\value{
+A tibble with one row per gage-day. Columns include:
+\describe{
+\item{site_number}{USGS site number.}
+\item{gage_name}{USGS station name.}
+\item{county_geoid}{Five-digit county FIPS code.}
+\item{county_name}{County name.}
+\item{state_abbreviation}{Two-letter state abbreviation.}
+\item{latitude, longitude}{Gage coordinates (decimal degrees).}
+\item{drainage_area_sqmi}{Upstream drainage area, in square miles.}
+\item{date}{Calendar day. For \code{statistic = "daily_max"}, the day is
+defined in the gage's local (Eastern) time zone.}
+\item{value}{The reading, in feet ("height") or cubic feet per second
+("discharge").}
+\item{approval_status}{"approved" or "provisional". For
+\code{statistic = "daily_max"}, a day is "provisional" if any reading that
+day is provisional.}
+}
+}
+\description{
+Pulls daily stream-gage readings, over each gage's full period of
+record by default, for USGS gages in one or more counties via the
+dataRetrieval package (USGS Water Data OGC APIs; these replace the
+now-decommissioning NWIS web services). Two statistics are supported:
+\itemize{
+\item "daily_mean": the published daily-mean series. Available for a
+century-plus at many gages.
+\item "daily_max": the maximum reading each day, computed here from the
+continuous (15-minute) record, because USGS publishes no daily-maximum
+gage-height series. Continuous records only begin in the mid-1990s at
+the earliest. These pulls are slow (a minute or more per long-record
+gage), so each gage's aggregated result is cached to its own parquet
+file and the pull resumes wherever it left off.
+}
+}
+\details{
+Site metadata (name, county, coordinates, drainage area) is attached
+to every reading. Data are from the USGS Water Data APIs; see
+\url{https://api.waterdata.usgs.gov/}. USGS asks that heavy users register
+a free API key (see \code{dataRetrieval::setAccess()} documentation); unkeyed
+access is rate-limited but sufficient for modest pulls.
+}
+\examples{
+\dontrun{
+## daily maximum gage heights in Upshur and Lewis Counties, WV
+get_usgs_gage(
+  counties = c("54097", "54041"),
+  cache_dir = tools::R_user_dir("climateapi", which = "cache"))
+
+get_usgs_gage(
+  counties = "54063",
+  measure = "discharge",
+  statistic = "daily_mean")
+}
+}
diff --git a/man/scrape_pda_pdfs.Rd b/man/scrape_pda_pdfs.Rd
index 21e3fc6..89bb5f0 100644
--- a/man/scrape_pda_pdfs.Rd
+++ b/man/scrape_pda_pdfs.Rd
@@ -7,18 +7,40 @@
 scrape_pda_pdfs(
   cache_directory = file.path(climateapi::get_box_path(), "hazards", "urban",
     "preliminary-damage-assessments", "pdfs"),
+  pages = NULL,
   max_pages = 200,
-  attempts_per_page = 3
+  attempts_per_page = 5,
+  delay_seconds = 2,
+  quiet = FALSE
 )
 }
 \arguments{
 \item{cache_directory}{The folder where scraped PDFs are written.}
 
+\item{pages}{Which listing pages to read, as a numeric vector. The default
+\code{NULL} walks the whole listing until a page returns no links, which is the
+only setting that guarantees complete coverage. FEMA lists newest first, so
+\code{pages = 0:2} is enough to pick up recently published reports and is much
+faster. Two caveats when this is set: cached files that are no longer listed
+cannot be detected, because the full listing was never read; and filename
+collisions are resolved only against the pages requested, so a name could be
+assigned that a page outside the range also claims.}
+
 \item{max_pages}{A guard against an unbounded walk if the listing ever stops
 returning empty pages. Raises an error if reached.}
 
 \item{attempts_per_page}{How many times to try a listing page before treating
 it as a failure.}
+
+\item{delay_seconds}{Seconds to pause between listing pages, to avoid
+hammering FEMA's site. This is what dominates the running time of a full
+walk -- roughly 70 pages -- not the downloads, which are skipped for reports
+already in \code{cache_directory}.}
+
+\item{quiet}{Suppress progress messages? The walk covers roughly 70 pages with
+a pause between each and is otherwise silent for several minutes, so progress
+is reported by default. Warnings about failed downloads and about cached
+files no longer listed are always raised, regardless of this setting.}
 }
 \value{
 Invisibly, a tibble with one row per report found on the site,
diff --git a/tests/testthat/test-get_usgs_gage.R b/tests/testthat/test-get_usgs_gage.R
new file mode 100644
index 0000000..7b296a5
--- /dev/null
+++ b/tests/testthat/test-get_usgs_gage.R
@@ -0,0 +1,128 @@
+# Tests for get_usgs_gage.R
+
+test_that("get_usgs_gage validates counties parameter", {
+  expect_error(get_usgs_gage(counties = 54097), "five-digit county FIPS")
+  expect_error(get_usgs_gage(counties = character(0)), "five-digit county FIPS")
+  expect_error(get_usgs_gage(counties = "WV"), "five-digit county FIPS")
+  expect_error(get_usgs_gage(counties = c("54097", "3904")), "five-digit county FIPS")
+})
+
+test_that("get_usgs_gage validates measure and statistic parameters", {
+  expect_error(get_usgs_gage(counties = "54097", measure = "00065"))
+  expect_error(get_usgs_gage(counties = "54097", measure = "gage_height"))
+  expect_error(get_usgs_gage(counties = "54097", statistic = "daily_median"))
+})
+
+test_that("get_usgs_gage validates date parameters", {
+  expect_error(
+    get_usgs_gage(counties = "54097", start_date = "01/01/2020"),
+    "`start_date` must be")
+  expect_error(
+    get_usgs_gage(counties = "54097", end_date = "2020-1-1"),
+    "`end_date` must be")
+  expect_error(
+    get_usgs_gage(counties = "54097", start_date = NULL),
+    "`start_date` must be")
+})
+
+test_that("get_usgs_gage function signature is correct", {
+  expect_true(is.function(get_usgs_gage))
+
+  params <- names(formals(get_usgs_gage))
+  expect_equal(
+    params,
+    c("counties", "measure", "statistic", "start_date", "end_date",
+      "refresh_cache", "cache_dir"))
+
+  f <- formals(get_usgs_gage)
+  expect_equal(eval(f$measure), c("height", "discharge"))
+  expect_equal(eval(f$statistic), c("daily_max", "daily_mean"))
+  expect_equal(f$start_date, "")
+  expect_equal(f$end_date, "")
+  expect_false(f$refresh_cache)
+})
+
+test_that("get_usgs_gage returns the cached file without hitting the API", {
+  cache_dir <- withr::local_tempdir()
+  cached_history <- tibble::tibble(
+    site_number = "03183500",
+    gage_name = "Greenbrier River at Alderson, WV",
+    county_geoid = "54063",
+    county_name = "Monroe",
+    state_abbreviation = "WV",
+    latitude = 37.7,
+    longitude = -80.6,
+    drainage_area_sqmi = 1364,
+    date = as.Date("2016-06-23"),
+    value = 15.5,
+    approval_status = "approved")
+
+  ## the cache file name encodes the counties, measure, and statistic; a
+  ## matching file must be returned as-is, with no network calls
+  arrow::write_parquet(
+    cached_history,
+    file.path(cache_dir, "gage_history_54063_height_daily_mean.parquet"))
+
+  result <- get_usgs_gage(
+    counties = "54063", measure = "height", statistic = "daily_mean",
+    cache_dir = cache_dir)
+  expect_equal(result, cached_history)
+})
+
+test_that("get_usgs_gage returns daily means with expected structure", {
+  skip_if_offline()
+
+  ## one month of daily-mean discharge across one county's gages keeps the
+  ## live pull small; Monroe County, WV (54063) contains the long-record
+  ## Greenbrier River at Alderson gage (03183500)
+  result <- tryCatch(
+    suppressWarnings(suppressMessages(get_usgs_gage(
+      counties = "54063",
+      measure = "discharge",
+      statistic = "daily_mean",
+      start_date = "2020-01-01",
+      end_date = "2020-01-31",
+      cache_dir = withr::local_tempdir()))),
+    error = function(e) NULL)
+  skip_if(is.null(result) || nrow(result) == 0, "Live USGS API data not available")
+
+  expect_equal(
+    names(result),
+    c("site_number", "gage_name", "county_geoid", "county_name",
+      "state_abbreviation", "latitude", "longitude", "drainage_area_sqmi",
+      "date", "value", "approval_status"))
+  expect_true("03183500" %in% result$site_number)
+  expect_true(all(result$county_geoid == "54063"))
+  expect_true(all(result$state_abbreviation == "WV"))
+  expect_true(all(result$date >= as.Date("2020-01-01") & result$date <= as.Date("2020-01-31")))
+  expect_true(all(result$approval_status %in% c("approved", "provisional")))
+  expect_true(all(!is.na(result$value)))
+})
+
+test_that("get_usgs_gage computes daily maxima from the continuous record", {
+  skip_if_offline()
+
+  ## a few days of the county's 15-minute records keeps the live pull small
+  cache_dir <- withr::local_tempdir()
+  result <- tryCatch(
+    suppressWarnings(suppressMessages(get_usgs_gage(
+      counties = "54063",
+      measure = "height",
+      statistic = "daily_max",
+      start_date = "2020-01-01",
+      end_date = "2020-01-05",
+      cache_dir = cache_dir))),
+    error = function(e) NULL)
+  skip_if(is.null(result) || nrow(result) == 0, "Live USGS API data not available")
+
+  expect_true(all(result$county_geoid == "54063"))
+  ## one row per gage-day
+  expect_equal(
+    nrow(result),
+    nrow(dplyr::distinct(result, site_number, date)))
+  expect_true(all(result$approval_status %in% c("approved", "provisional")))
+
+  ## the per-site parquet caches enable resumable pulls
+  expect_gt(
+    length(list.files(file.path(cache_dir, "by-site-height-daily-max"))), 0)
+})
diff --git a/vignettes/figure/get_hrrr_smoke-animation.gif b/vignettes/figure/get_hrrr_smoke-animation.gif
new file mode 100644
index 0000000..e119178
Binary files /dev/null and b/vignettes/figure/get_hrrr_smoke-animation.gif differ
diff --git a/vignettes/figure/get_usgs_gage-plot-daily-maxima-1.png b/vignettes/figure/get_usgs_gage-plot-daily-maxima-1.png
new file mode 100644
index 0000000..ffc9973
Binary files /dev/null and b/vignettes/figure/get_usgs_gage-plot-daily-maxima-1.png differ
diff --git a/vignettes/get_hrrr_smoke.Rmd b/vignettes/get_hrrr_smoke.Rmd
new file mode 100644
index 0000000..b8390f7
--- /dev/null
+++ b/vignettes/get_hrrr_smoke.Rmd
@@ -0,0 +1,169 @@
+---
+title: "Animating Wildfire Smoke"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Animating Wildfire Smoke}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+
+
+`get_hrrr_smoke()` retrieves hourly near-surface wildfire smoke concentrations
+from NOAA's High-Resolution Rapid Refresh (HRRR) model for any area in the
+conterminous United States, on a 3-kilometer grid. Here we build a two-week
+look-behind animation of smoke over Washington State: where smoke traveled,
+when it arrived, and how severe it was at ground level.
+
+
+``` r
+library(climateapi)
+library(dplyr)
+library(stringr)
+library(sf)
+library(ggplot2)
+library(urbnthemes)
+library(gganimate)
+```
+
+## Pulling two weeks of hourly smoke
+
+We ask for the "surface" variable -- smoke mass density 8 meters above ground,
+in micrograms per cubic meter, the quantity most comparable to PM2.5 air
+quality readings. The result is a `terra::SpatRaster` with one layer per
+retrieved hour, timestamped in UTC. Three-hourly steps (eight frames per day)
+keep the animation responsive over a two-week window; pass `hours = 0:23` for
+the full hourly record.
+
+
+``` r
+projection = 5070
+
+area_of_interest = tigris::states(cb = TRUE, year = 2023, progress_bar = FALSE) %>%
+  filter(str_detect(NAME, "Washington")) %>%
+  st_transform(projection)
+
+smoke_data = get_hrrr_smoke(
+  geometries = area_of_interest,
+  start_date = Sys.Date() - 14,
+  end_date = Sys.Date(),
+  variable = "surface",
+  hours = seq(0, 21, by = 3))
+#> 
|---------|---------|---------|---------|
=========================================
                                          

+
+smoke_data
+#> class       : SpatRaster 
+#> size        : 150, 200, 117  (nrow, ncol, nlyr)
+#> resolution  : 3000, 3000  (x, y)
+#> extent      : -2036020, -1436020, 991193.8, 1441194  (xmin, xmax, ymin, ymax)
+#> coord. ref. : +proj=lcc +lat_0=38.5 +lon_0=-97.5 +lat_1=38.5 +lat_2=38.5 +x_0=0 +y_0=0 +R=6371229 +units=m +no_defs 
+#> source(s)   : memory
+#> names       :  2026-~00:00,  2026-~03:00,  2026-~06:00,  2026-~09:00,  2026-~12:00,  2026-~15:00, ... 
+#> min values  : 2.988507e-15, 1.077663e-16, 1.053054e-16, 9.260972e-17, 3.603632e-14, 2.684881e-08, ... 
+#> max values  : 4.555040e+03, 4.114960e+03, 1.005960e+03, 6.965000e+02, 1.654240e+03, 1.312560e+03, ... 
+#> time        : 2026-07-23 to 2026-08-06 12:00:00 UTC (117 steps)
+```
+
+## Preparing the data for mapping
+
+The raster arrives in the HRRR model's native projection; we reproject it to
+match the vector layers, then melt it into a one-row-per-cell-per-hour tibble,
+which is the shape `gganimate` needs.
+
+Rather than mapping concentrations to a continuous color ramp, we bin them at
+the EPA's 24-hour PM2.5 Air Quality Index breakpoints (rounded for display):
+9, 35, 55, 125, and 225 micrograms per cubic meter mark the transitions from
+"good" air through "moderate", "unhealthy for sensitive groups", "unhealthy",
+"very unhealthy", and "hazardous". Cells below 1 microgram per cubic meter are
+dropped entirely so that clean air stays transparent and the basemap shows
+through.
+
+
+``` r
+smoke_data_projected = terra::project(smoke_data, paste0("EPSG:", projection))
+
+counties_context = tigris::counties(
+    cb = TRUE, year = 2023, state = "WA", progress_bar = FALSE) %>%
+  st_transform(projection) %>%
+  select(NAME)
+
+smoke_breaks = c(1, 9, 35, 55, 125, 225, Inf)
+smoke_labels = c("1–9", "9–35", "35–55", "55–125", "125–225", "225+")
+smoke_colors = palette_urbn_cyan[2:7]
+names(smoke_colors) = smoke_labels
+
+smoke_df = smoke_data_projected %>%
+  as.data.frame(xy = TRUE, wide = FALSE) %>%
+  as_tibble() %>%
+  left_join(
+    tibble(
+      layer = names(smoke_data_projected),
+      timestamp = terra::time(smoke_data_projected)),
+    by = "layer",
+    relationship = "many-to-one") %>%
+  rename(smoke_ug_m3 = values) %>%
+  mutate(
+    smoke_level = cut(
+      smoke_ug_m3,
+      breaks = smoke_breaks,
+      labels = smoke_labels,
+      right = FALSE)) %>%
+  filter(!is.na(smoke_level))
+```
+
+## Animating
+
+The plot is an ordinary ggplot -- light-filled counties first, the smoke
+raster on top, an outline-only state border last -- plus
+`transition_time(timestamp)`, which turns the hourly layers into frames.
+`animate()` renders one frame per model hour, and `anim_save()` writes the
+result next to the vignette's other figures so it can be embedded below.
+
+
+``` r
+smoke_animation = ggplot() +
+  geom_sf(data = counties_context, fill = "#f5f5f5", color = "#d2d2d2", linewidth = 0.35) +
+  geom_raster(data = smoke_df, aes(x = x, y = y, fill = smoke_level), alpha = 0.8) +
+  geom_sf(data = area_of_interest, fill = NA, color = "#696969", linewidth = 0.35) +
+  scale_fill_manual(
+    values = smoke_colors,
+    ## keep every severity bin in the legend even in frames where no cell
+    ## reaches it, so the legend does not change size between frames
+    drop = FALSE,
+    name = expression("Near-surface smoke (" * mu * "g/m"^3 * ")")) +
+  labs(
+    title = "Smoke over Washington State",
+    subtitle = "Time: {format(frame_time, '%B %d, %H:%M UTC')}") +
+  urbnthemes::theme_urbn_map() +
+  theme(
+    legend.position = "bottom",
+    legend.justification = "left",
+    legend.direction = "horizontal",
+    legend.title.position = "top",
+    legend.key.spacing.x = unit(0.25, "line")) +
+  transition_time(timestamp) +
+  ease_aes("linear")
+
+smoke_gif = animate(
+  smoke_animation,
+  renderer = gifski_renderer(),
+  device = "ragg_png",
+  nframes = n_distinct(smoke_df$timestamp),
+  fps = 8,
+  width = 900,
+  height = 700,
+  units = "px",
+  res = 120,
+  end_pause = 8)
+
+anim_save("figure/get_hrrr_smoke-animation.gif", smoke_gif)
+```
+
+<img src="figure/get_hrrr_smoke-animation.gif" alt="An animated map of Washington State showing modeled near-surface wildfire smoke concentrations hour by hour over a two-week period. Smoke plumes, shaded from light to dark blue by severity, drift and pool across the state as the timestamp advances." />
+
+The model's spatial patterns -- where plumes travel and pool -- are generally
+reliable, but the concentrations are estimates driven by satellite fire
+detections, so they can be biased where fire emission estimates are wrong.
+For observed ground truth at specific locations, compare against PM2.5
+monitor readings (for example, via the AirNow API), which share the same
+unit.
diff --git a/vignettes/get_hrrr_smoke.Rmd.orig b/vignettes/get_hrrr_smoke.Rmd.orig
new file mode 100644
index 0000000..2dacee7
--- /dev/null
+++ b/vignettes/get_hrrr_smoke.Rmd.orig
@@ -0,0 +1,164 @@
+---
+title: "Animating Wildfire Smoke"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Animating Wildfire Smoke}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  warning = FALSE,
+  message = FALSE,
+  fig.width = 8,
+  fig.height = 6,
+  dpi = 150
+)
+```
+
+`get_hrrr_smoke()` retrieves hourly near-surface wildfire smoke concentrations
+from NOAA's High-Resolution Rapid Refresh (HRRR) model for any area in the
+conterminous United States, on a 3-kilometer grid. Here we build a two-week
+look-behind animation of smoke over Washington State: where smoke traveled,
+when it arrived, and how severe it was at ground level.
+
+```{r setup}
+library(climateapi)
+library(dplyr)
+library(stringr)
+library(sf)
+library(ggplot2)
+library(urbnthemes)
+library(gganimate)
+```
+
+## Pulling two weeks of hourly smoke
+
+We ask for the "surface" variable -- smoke mass density 8 meters above ground,
+in micrograms per cubic meter, the quantity most comparable to PM2.5 air
+quality readings. The result is a `terra::SpatRaster` with one layer per
+retrieved hour, timestamped in UTC. Three-hourly steps (eight frames per day)
+keep the animation responsive over a two-week window; pass `hours = 0:23` for
+the full hourly record.
+
+```{r pull-smoke}
+projection = 5070
+
+area_of_interest = tigris::states(cb = TRUE, year = 2023, progress_bar = FALSE) %>%
+  filter(str_detect(NAME, "Washington")) %>%
+  st_transform(projection)
+
+smoke_data = get_hrrr_smoke(
+  geometries = area_of_interest,
+  start_date = Sys.Date() - 14,
+  end_date = Sys.Date(),
+  variable = "surface",
+  hours = seq(0, 21, by = 3))
+
+smoke_data
+```
+
+## Preparing the data for mapping
+
+The raster arrives in the HRRR model's native projection; we reproject it to
+match the vector layers, then melt it into a one-row-per-cell-per-hour tibble,
+which is the shape `gganimate` needs.
+
+Rather than mapping concentrations to a continuous color ramp, we bin them at
+the EPA's 24-hour PM2.5 Air Quality Index breakpoints (rounded for display):
+9, 35, 55, 125, and 225 micrograms per cubic meter mark the transitions from
+"good" air through "moderate", "unhealthy for sensitive groups", "unhealthy",
+"very unhealthy", and "hazardous". Cells below 1 microgram per cubic meter are
+dropped entirely so that clean air stays transparent and the basemap shows
+through.
+
+```{r prepare-smoke}
+smoke_data_projected = terra::project(smoke_data, paste0("EPSG:", projection))
+
+counties_context = tigris::counties(
+    cb = TRUE, year = 2023, state = "WA", progress_bar = FALSE) %>%
+  st_transform(projection) %>%
+  select(NAME)
+
+smoke_breaks = c(1, 9, 35, 55, 125, 225, Inf)
+smoke_labels = c("1–9", "9–35", "35–55", "55–125", "125–225", "225+")
+smoke_colors = palette_urbn_cyan[2:7]
+names(smoke_colors) = smoke_labels
+
+smoke_df = smoke_data_projected %>%
+  as.data.frame(xy = TRUE, wide = FALSE) %>%
+  as_tibble() %>%
+  left_join(
+    tibble(
+      layer = names(smoke_data_projected),
+      timestamp = terra::time(smoke_data_projected)),
+    by = "layer",
+    relationship = "many-to-one") %>%
+  rename(smoke_ug_m3 = values) %>%
+  mutate(
+    smoke_level = cut(
+      smoke_ug_m3,
+      breaks = smoke_breaks,
+      labels = smoke_labels,
+      right = FALSE)) %>%
+  filter(!is.na(smoke_level))
+```
+
+## Animating
+
+The plot is an ordinary ggplot -- light-filled counties first, the smoke
+raster on top, an outline-only state border last -- plus
+`transition_time(timestamp)`, which turns the hourly layers into frames.
+`animate()` renders one frame per model hour, and `anim_save()` writes the
+result next to the vignette's other figures so it can be embedded below.
+
+```{r animate-smoke, results = "hide"}
+smoke_animation = ggplot() +
+  geom_sf(data = counties_context, fill = "#f5f5f5", color = "#d2d2d2", linewidth = 0.35) +
+  geom_raster(data = smoke_df, aes(x = x, y = y, fill = smoke_level), alpha = 0.8) +
+  geom_sf(data = area_of_interest, fill = NA, color = "#696969", linewidth = 0.35) +
+  scale_fill_manual(
+    values = smoke_colors,
+    ## keep every severity bin in the legend even in frames where no cell
+    ## reaches it, so the legend does not change size between frames
+    drop = FALSE,
+    name = expression("Near-surface smoke (" * mu * "g/m"^3 * ")")) +
+  labs(
+    title = "Smoke over Washington State",
+    subtitle = "Time: {format(frame_time, '%B %d, %H:%M UTC')}") +
+  urbnthemes::theme_urbn_map() +
+  theme(
+    legend.position = "bottom",
+    legend.justification = "left",
+    legend.direction = "horizontal",
+    legend.title.position = "top",
+    legend.key.spacing.x = unit(0.25, "line")) +
+  transition_time(timestamp) +
+  ease_aes("linear")
+
+smoke_gif = animate(
+  smoke_animation,
+  renderer = gifski_renderer(),
+  device = "ragg_png",
+  nframes = n_distinct(smoke_df$timestamp),
+  fps = 8,
+  width = 900,
+  height = 700,
+  units = "px",
+  res = 120,
+  end_pause = 8)
+
+anim_save("figure/get_hrrr_smoke-animation.gif", smoke_gif)
+```
+
+<img src="figure/get_hrrr_smoke-animation.gif" alt="An animated map of Washington State showing modeled near-surface wildfire smoke concentrations hour by hour over a two-week period. Smoke plumes, shaded from light to dark blue by severity, drift and pool across the state as the timestamp advances." />
+
+The model's spatial patterns -- where plumes travel and pool -- are generally
+reliable, but the concentrations are estimates driven by satellite fire
+detections, so they can be biased where fire emission estimates are wrong.
+For observed ground truth at specific locations, compare against PM2.5
+monitor readings (for example, via the AirNow API), which share the same
+unit.
diff --git a/vignettes/get_usgs_gage.Rmd b/vignettes/get_usgs_gage.Rmd
new file mode 100644
index 0000000..527ed32
--- /dev/null
+++ b/vignettes/get_usgs_gage.Rmd
@@ -0,0 +1,112 @@
+---
+title: "USGS Stream-Gage Histories"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{USGS Stream-Gage Histories}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+
+
+`get_usgs_gage()` pulls daily stream-gage readings for every USGS stream gage
+in one or more counties, attaching each gage's name, county, coordinates, and
+drainage area to every reading. Under the hood it uses the USGS Water Data
+APIs via the `dataRetrieval` package. Here we ask a concrete question: at the
+stream gages in Upshur and Lewis Counties, West Virginia, how do daily maximum
+water levels in 2026 compare against each gage's own history?
+
+
+``` r
+library(climateapi)
+library(dplyr)
+library(ggplot2)
+library(lubridate)
+library(stringr)
+library(urbnthemes)
+
+set_urbn_defaults(style = "print")
+```
+
+## Pulling daily maximum gage heights
+
+Counties are identified by their five-digit FIPS codes: "54097" is Upshur
+County and "54041" is Lewis County. The daily maximum height is computed from
+each gage's continuous (15-minute) record -- a download of a minute or more
+per gage -- so each gage's aggregated daily maxima are cached to their own
+parquet file: an interrupted pull resumes where it left off, and a re-run
+reads entirely from disk. We cache under the standard per-user cache
+directory; any persistent directory works.
+
+
+``` r
+gage_history = get_usgs_gage(
+  counties = c("54097", "54041"),
+  measure = "height",
+  statistic = "daily_max",
+  cache_dir = tools::R_user_dir("climateapi", which = "cache"))
+
+gage_history %>% dplyr::glimpse()
+#> Rows: 41,572
+#> Columns: 11
+#> $ site_number        <chr> "03052120", "03052120", "03052120", "03052120", "03…
+#> $ gage_name          <chr> "BUCKHANNON RIVER AT ALTON, WV", "BUCKHANNON RIVER …
+#> $ county_geoid       <chr> "54097", "54097", "54097", "54097", "54097", "54097…
+#> $ county_name        <chr> "Upshur", "Upshur", "Upshur", "Upshur", "Upshur", "…
+#> $ state_abbreviation <chr> "WV", "WV", "WV", "WV", "WV", "WV", "WV", "WV", "WV…
+#> $ latitude           <dbl> 38.81964, 38.81964, 38.81964, 38.81964, 38.81964, 3…
+#> $ longitude          <dbl> -80.21389, -80.21389, -80.21389, -80.21389, -80.213…
+#> $ drainage_area_sqmi <dbl> 94.7, 94.7, 94.7, 94.7, 94.7, 94.7, 94.7, 94.7, 94.…
+#> $ date               <date> 2011-10-12, 2011-10-13, 2011-10-14, 2011-10-15, 20…
+#> $ value              <dbl> 5.64, 5.66, 7.89, 7.52, 6.65, 6.22, 5.95, 5.77, 5.8…
+#> $ approval_status    <chr> "approved", "approved", "approved", "approved", "ap…
+```
+
+## How does 2026 compare to each gage's history?
+
+We overlay one line per year at each gage, aligning every year on a shared
+January-through-December axis. The 2026 line (partial, through early August)
+is drawn in the Urban Institute's main blue over each gage's earlier years in
+light gray.
+
+
+``` r
+plot_data = gage_history %>%
+  dplyr::mutate(
+    year = lubridate::year(date),
+    ## a common x-axis date so that years plot atop one another; the (arbitrary)
+    ## year 2026 is used only for axis labeling
+    common_date = as.Date(lubridate::yday(date) - 1, origin = "2026-01-01"),
+    ## shorten station names so facet labels fit their strips
+    gage_label = gage_name %>%
+      stringr::str_to_title() %>%
+      stringr::str_remove(",?\\s*(Wv|WV)$") %>%
+      stringr::str_trunc(30) %>%
+      stringr::str_c(" (", site_number, ")"))
+
+plot_data %>%
+  ggplot(aes(x = common_date, y = value, group = year)) +
+  geom_line(
+    data = ~ dplyr::filter(.x, year != 2026),
+    color = palette_urbn_gray[5],
+    linewidth = 0.3) +
+  geom_line(
+    data = ~ dplyr::filter(.x, year == 2026),
+    color = palette_urbn_main[["cyan"]],
+    linewidth = 0.7) +
+  facet_wrap(~ gage_label, ncol = 2, scales = "free_y") +
+  scale_x_date(date_breaks = "3 months", date_labels = "%b") +
+  labs(
+    x = NULL,
+    y = "Maximum daily gage height (feet)",
+    title = "Daily maximum water levels in Upshur and Lewis Counties, WV",
+    subtitle = "One line per year; 2026 (through early August) in blue")
+```
+
+<img src="figure/get_usgs_gage-plot-daily-maxima-1.png" alt="Faceted line charts, one per stream gage in Upshur and Lewis Counties, West Virginia, each showing daily maximum gage height across the calendar year with one line per year since the late 2000s in light gray and 2026 highlighted in blue." />
+
+Because each series reflects the maximum reading on each day, short-lived
+flood crests are visible even when they lasted only hours. Note that the
+vertical axis varies by gage: gage height is measured relative to a
+gage-specific datum, so heights are comparable across years at the same gage
+but not across gages.
diff --git a/vignettes/get_usgs_gage.Rmd.orig b/vignettes/get_usgs_gage.Rmd.orig
new file mode 100644
index 0000000..08d0658
--- /dev/null
+++ b/vignettes/get_usgs_gage.Rmd.orig
@@ -0,0 +1,104 @@
+---
+title: "USGS Stream-Gage Histories"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{USGS Stream-Gage Histories}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  warning = FALSE,
+  message = FALSE,
+  fig.width = 8,
+  fig.height = 9,
+  dpi = 150
+)
+```
+
+`get_usgs_gage()` pulls daily stream-gage readings for every USGS stream gage
+in one or more counties, attaching each gage's name, county, coordinates, and
+drainage area to every reading. Under the hood it uses the USGS Water Data
+APIs via the `dataRetrieval` package. Here we ask a concrete question: at the
+stream gages in Upshur and Lewis Counties, West Virginia, how do daily maximum
+water levels in 2026 compare against each gage's own history?
+
+```{r setup}
+library(climateapi)
+library(dplyr)
+library(ggplot2)
+library(lubridate)
+library(stringr)
+library(urbnthemes)
+
+set_urbn_defaults(style = "print")
+```
+
+## Pulling daily maximum gage heights
+
+Counties are identified by their five-digit FIPS codes: "54097" is Upshur
+County and "54041" is Lewis County. The daily maximum height is computed from
+each gage's continuous (15-minute) record -- a download of a minute or more
+per gage -- so each gage's aggregated daily maxima are cached to their own
+parquet file: an interrupted pull resumes where it left off, and a re-run
+reads entirely from disk. We cache under the standard per-user cache
+directory; any persistent directory works.
+
+```{r pull-gage-history}
+gage_history = get_usgs_gage(
+  counties = c("54097", "54041"),
+  measure = "height",
+  statistic = "daily_max",
+  cache_dir = tools::R_user_dir("climateapi", which = "cache"))
+
+gage_history %>% dplyr::glimpse()
+```
+
+## How does 2026 compare to each gage's history?
+
+We overlay one line per year at each gage, aligning every year on a shared
+January-through-December axis. The 2026 line (partial, through early August)
+is drawn in the Urban Institute's main blue over each gage's earlier years in
+light gray.
+
+```{r plot-daily-maxima, fig.alt = "Faceted line charts, one per stream gage in Upshur and Lewis Counties, West Virginia, each showing daily maximum gage height across the calendar year with one line per year since the late 2000s in light gray and 2026 highlighted in blue."}
+plot_data = gage_history %>%
+  dplyr::mutate(
+    year = lubridate::year(date),
+    ## a common x-axis date so that years plot atop one another; the (arbitrary)
+    ## year 2026 is used only for axis labeling
+    common_date = as.Date(lubridate::yday(date) - 1, origin = "2026-01-01"),
+    ## shorten station names so facet labels fit their strips
+    gage_label = gage_name %>%
+      stringr::str_to_title() %>%
+      stringr::str_remove(",?\\s*(Wv|WV)$") %>%
+      stringr::str_trunc(30) %>%
+      stringr::str_c(" (", site_number, ")"))
+
+plot_data %>%
+  ggplot(aes(x = common_date, y = value, group = year)) +
+  geom_line(
+    data = ~ dplyr::filter(.x, year != 2026),
+    color = palette_urbn_gray[5],
+    linewidth = 0.3) +
+  geom_line(
+    data = ~ dplyr::filter(.x, year == 2026),
+    color = palette_urbn_main[["cyan"]],
+    linewidth = 0.7) +
+  facet_wrap(~ gage_label, ncol = 2, scales = "free_y") +
+  scale_x_date(date_breaks = "3 months", date_labels = "%b") +
+  labs(
+    x = NULL,
+    y = "Maximum daily gage height (feet)",
+    title = "Daily maximum water levels in Upshur and Lewis Counties, WV",
+    subtitle = "One line per year; 2026 (through early August) in blue")
+```
+
+Because each series reflects the maximum reading on each day, short-lived
+flood crests are visible even when they lasted only hours. Note that the
+vertical axis varies by gage: gage height is measured relative to a
+gage-specific datum, so heights are comparable across years at the same gage
+but not across gages.