diff --git a/DESCRIPTION b/DESCRIPTION
index 4c9ceeb..a748291 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -51,6 +51,7 @@ Imports:
sf,
stats,
stringr,
+ terra,
tibble,
tidycensus,
tidyr,
diff --git a/NAMESPACE b/NAMESPACE
index 79a9361..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)
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/_pkgdown.yml b/_pkgdown.yml
index 97a467e..1553fa3 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -9,6 +9,7 @@ reference:
- 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
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/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/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)
+```
+
+
+
+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)
+```
+
+
+
+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.