USGS Stream Discharge Analysis

1. Packages

# Run once to install (remove the # to uncomment, run, then re-comment):
# install.packages(c("dataRetrieval", "tidyverse", "lubridate",
#                    "scales", "knitr", "kableExtra"))

library(dataRetrieval)  # official USGS R package — fetches NWIS data
library(tidyverse)      # ggplot2, dplyr, tidyr
library(lubridate)      # date arithmetic
library(scales)         # comma-formatted and log axis labels
library(knitr)          # kable() tables
library(kableExtra)     # styled HTML tables

📄 What is dataRetrieval? dataRetrieval is maintained by the USGS itself and provides direct R access to the National Water Information System (NWIS) — the same database behind waterdata.usgs.gov. Instead of navigating a web interface and downloading CSVs, you call one function and the data lands directly in a tidy R data frame. The package also has functions for water-quality data, rating curves, and site metadata.

2. Settings — The Only Lines You Need to Change

All user-defined parameters live in this one chunk. To analyse a different river, change only these values and re-knit. Everything else — titles, axis labels, table captions — updates automatically.

# ── What each parameter means ─────────────────────────────────
site_no   # 8-digit USGS site identifier
param_cd # 00060 = discharge (cfs) | 00065 = gage height (ft)
stat_cd   # 00003 = daily mean (the standard)
start_dt # YYYY-MM-DD — at least 2 full years recommended
end_dt   # YYYY-MM-DD

# ── Part 1: run as-is first ───────────────────────────────────────────────────
# The San Gabriel River at Azusa is a local reference site with a long,
# continuous record. Confirm you get a complete output before changing anything.
site_no  <- "11119750"   # Ventura River near Ventura, CA
param_cd <- "00060"      # discharge (cfs)
stat_cd  <- "00003"      # daily mean
start_dt <- "2020-01-01"
end_dt   <- "2023-12-31"

# ── Part 2 (assignment): replace with your chosen station ─────────────────────
# See Section 7 for instructions on finding a station in any state.

💡 Southern California default. The default site is the Ventura River near Ventura, CA — a flashy, Mediterranean-climate stream whose discharge is dominated by winter storm events and drops to near zero in summer. This makes it a useful contrast to the groundwater-sustained rivers of the eastern US. You should see a strongly seasonal signal in the monthly boxplot.

3. Data Download

3.1 Query the NWIS API

readNWISdv() sends an HTTP request to the USGS NWIS daily-values API and returns a data frame. renameNWISColumns() replaces the auto-generated column name (e.g. X_00060_00003) with the readable name Flow. readNWISsite() retrieves station metadata — official name, coordinates, drainage area — so plot titles are generated automatically.

raw  <- readNWISdv(
  siteNumbers = site_no,
  parameterCd = param_cd,
  startDate   = start_dt,
  endDate     = end_dt
)

raw  <- renameNWISColumns(raw)
meta <- readNWISsite(site_no)

cat("Rows downloaded:", nrow(raw), "\n")

## Rows downloaded: 1461

cat("Columns        :", paste(names(raw), collapse = ", "), "\n")

## Columns        : agency_cd, site_no, Date, Flow, Flow_cd

cat("Site name      :", meta$station_nm, "\n")

## Site name      : MISSION C NR MISSION ST NR SANTA BARBARA CA

cat("Drainage area  :", meta$drain_area_va, "sq mi\n")

## Drainage area  : 8.6 sq mi

⚠️ Network errors. readNWISdv() requires an internet connection. On a campus network with a VPN or firewall, the call occasionally times out even when the USGS server is healthy. If you get a connection error, try again in a moment. If the problem persists, use the CSV fallback in section 3.2 below.

3.2 Alternative: load a pre-downloaded CSV

If the API call fails, download the data manually from waterdata.usgs.gov and load it here. Change eval=FALSE to eval=TRUE on this chunk and eval=FALSE on the download chunk above.

# raw_csv <- read_csv("data/usgs_discharge.csv") %>%
#   mutate(time = ymd_hms(time)) %>%
#   filter(parameter_code == "00060") %>%
#   rename(Date = time, Flow = value) %>%
#   filter(approval_status == "A")
#
# raw <- raw_csv %>% select(Date, Flow)

4. Data Cleaning

This chunk does four things:

Selects only Date and Flow — all other NWIS metadata columns are dropped.
Removes rows where Flow is NA (equipment malfunctions and provisional data gaps appear as NA).
Derives three helper columns: year, month (Jan–Dec ordered factor), and doy (day of year, 1–366).
Stores the site name and unit label for automatic plot titles.

df <- raw %>%
  select(Date, Flow) %>%
  filter(!is.na(Flow)) %>%
  mutate(
    year  = year(Date),
    month = month(Date, label = TRUE, abbr = TRUE),
    doy   = yday(Date)
  )

site_name  <- meta$station_nm
units_lab  <- "Discharge (cfs)"

cat("Records after cleaning:", nrow(df), "\n")

## Records after cleaning: 1461

cat("Date range            :", as.character(min(df$Date)),
    "to", as.character(max(df$Date)), "\n")

## Date range            : 2020-01-01 to 2023-12-31

cat("Years covered         :", paste(sort(unique(df$year)), collapse = ", "), "\n")

## Years covered         : 2020, 2021, 2022, 2023

cat("Missing days removed  :", nrow(raw) - nrow(df), "\n")

## Missing days removed  : 0

📌 Why are there NA values in USGS data? USGS gauges measure water-surface elevation (stage) with a pressure transducer or float, then convert to discharge via a rating curve. NAs appear when the sensor malfunctions, when the gauge is underwater during an extreme flood (outside the rating curve range), or when the record is provisional and awaiting quality review. NWIS marks approved data with code A and provisional data with P — renameNWISColumns() retains this as a cd column you can filter on if needed.

5. Summary Statistics

stats <- df %>%
  group_by(year) %>%
  summarise(
    N      = n(),
    Mean   = round(mean(Flow),   1),
    Median = round(median(Flow), 1),
    SD     = round(sd(Flow),     1),
    Min    = round(min(Flow),    1),
    Max    = round(max(Flow),    1),
    .groups = "drop"
  )

stats %>%
  kbl(
    caption     = paste("Annual discharge summary \u2014", site_name),
    col.names   = c("Year", "N days", "Mean", "Median",
                    "SD", "Min", "Max"),
    format.args = list(big.mark = ",")
  ) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  add_header_above(c(" " = 2, "Discharge (cfs)" = 5))

Annual discharge summary — MISSION C NR MISSION ST NR SANTA BARBARA CA
		Discharge (cfs)
Year	N days	Mean	Median	SD	Min	Max
2,020	366	1.0	0	7.3	0	97.0
2,021	365	0.7	0	5.7	0	88.2
2,022	365	0.4	0	5.0	0	75.9
2,023	365	11.0	0	62.9	0	1,050.0

6. Standard Hydrological Plots

6.1 Daily time series — linear scale

A linear y-axis gives the most intuitive first view of the record. Large flood peaks stand out clearly. For a Mediterranean-climate stream like the Ventura River, you will see near-zero discharge through summer and sharp, narrow peaks in winter.

ggplot(df, aes(x = Date, y = Flow)) +
  geom_line(color = "#2e7d94", linewidth = 0.5) +
  scale_x_date(date_labels = "%b %Y", date_breaks = "6 months") +
  scale_y_continuous(labels = comma) +
  labs(
    title    = paste("Daily Mean Discharge \u2014", site_name),
    subtitle = paste("USGS site", site_no, "\u00b7 linear y-axis"),
    x        = NULL,
    y        = units_lab
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(family = "serif", face = "bold",
                                 color = "#1e3a2f", size = 13),
    plot.subtitle = element_text(color = "#7a7a72", size = 10),
    axis.text.x   = element_text(angle = 35, hjust = 1),
    panel.grid.minor = element_blank()
  )

6.2 Daily time series — log₁₀ scale

Stream discharge is highly skewed: flood peaks can be 100–1,000 times larger than base flows. A log₁₀ y-axis compresses the peaks and stretches the low-flow period so both extremes are readable simultaneously. annotation_logticks() adds the characteristic minor tick marks seen on logarithmic axes.

ggplot(df, aes(x = Date, y = Flow)) +
  geom_line(color = "#2e7d94", linewidth = 0.5) +
  scale_x_date(date_labels = "%b %Y", date_breaks = "6 months") +
  scale_y_log10(labels = comma) +
  annotation_logticks(sides = "l") +
  labs(
    title    = paste("Daily Mean Discharge (log scale) \u2014", site_name),
    subtitle = paste("USGS site", site_no, "\u00b7 log\u2081\u2080 y-axis"),
    x        = NULL,
    y        = "log\u2081\u2080 Discharge (cfs)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(family = "serif", face = "bold",
                                 color = "#1e3a2f", size = 13),
    plot.subtitle = element_text(color = "#7a7a72", size = 10),
    axis.text.x   = element_text(angle = 35, hjust = 1),
    panel.grid.minor = element_blank()
  )

📌 When to use log scale. Use linear scale when you want to show the magnitude of peaks (flood hazard, infrastructure design). Use log scale when you want to see low-flow variability (drought, ecological base flow, groundwater connectivity). For most stream analysis, showing both is standard.

6.3 Monthly discharge boxplots

Grouping by calendar month and pooling all years reveals the seasonal pattern of discharge, independent of year-to-year variability. Each box shows the interquartile range, the line is the median, and points beyond the whiskers are outliers.

ggplot(df, aes(x = month, y = Flow, fill = month)) +
  geom_boxplot(
    outlier.size  = 0.9,
    outlier.alpha = 0.4,
    show.legend   = FALSE
  ) +
  scale_y_log10(labels = comma) +
  annotation_logticks(sides = "l") +
  scale_fill_manual(
    values = colorRampPalette(c("#d4e6da", "#1e3a2f"))(12)
  ) +
  labs(
    title    = paste("Monthly Discharge Distribution \u2014", site_name),
    subtitle = "All years pooled  \u00b7  log\u2081\u2080 y-axis",
    x        = NULL,
    y        = "log\u2081\u2080 Discharge (cfs)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(family = "serif", face = "bold",
                                 color = "#1e3a2f", size = 13),
    plot.subtitle = element_text(color = "#7a7a72", size = 10),
    panel.grid.major.x = element_blank(),
    panel.grid.minor   = element_blank()
  )

7. Flow Duration Curve

The flow duration curve (FDC) is one of the most widely used tools in hydrology. It answers: “What discharge is equalled or exceeded X% of the time?” The curve is constructed by ranking all daily values from highest to lowest and computing each value’s exceedance probability.

fdc <- df %>%
  arrange(desc(Flow)) %>%
  mutate(exceed_pct = 100 * row_number() / n())

ggplot(fdc, aes(x = exceed_pct, y = Flow)) +
  geom_line(color = "#2e7d94", linewidth = 0.9) +
  scale_y_log10(labels = comma) +
  annotation_logticks(sides = "l") +
  scale_x_continuous(
    breaks = seq(0, 100, 10),
    labels = paste0(seq(0, 100, 10), "%")
  ) +
  labs(
    title    = paste("Flow Duration Curve \u2014", site_name),
    subtitle = paste("USGS site", site_no),
    x        = "Percent of time flow equalled or exceeded",
    y        = "log\u2081\u2080 Discharge (cfs)"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title    = element_text(family = "serif", face = "bold",
                                 color = "#1e3a2f", size = 13),
    plot.subtitle = element_text(color = "#7a7a72", size = 10),
    panel.grid.minor = element_blank()
  )

📌 Reading the flow duration curve.

Steep FDC — the curve drops sharply from high to low discharge over a narrow range of exceedance percentages. This indicates a flashy, rainfall-dominated stream: very high peak flows but very low base flows. Most Southern California streams look like this.
Flat FDC — discharge changes slowly across the full range of exceedance probabilities. This indicates a groundwater-sustained stream with stable year-round flow, typical of the Pacific Northwest or the Appalachians.
The Q90 value — the discharge exceeded 90% of the time — is a common index of ecological base flow and drought severity. The Q10 value is often used as a flood reference flow.

# Extract key exceedance quantiles
q_vals <- quantile(df$Flow, probs = c(0.90, 0.50, 0.10),
                   na.rm = TRUE)
cat(sprintf("Q90 (exceeded 90%% of days): %7.1f cfs\n", q_vals["10%"]))

## Q90 (exceeded 90% of days):     0.0 cfs

cat(sprintf("Q50 / median discharge:     %7.1f cfs\n", q_vals["50%"]))

## Q50 / median discharge:         0.0 cfs

cat(sprintf("Q10 (exceeded 10%% of days): %7.1f cfs\n", q_vals["90%"]))

## Q10 (exceeded 10% of days):     1.9 cfs

cat(sprintf("Flashiness ratio (Q10/Q90): %7.1f\n",
            q_vals["90%"] / q_vals["10%"]))

## Flashiness ratio (Q10/Q90):     Inf

⚙️ The flashiness ratio. Q10 / Q90 is a simple index of stream flashiness — the ratio of a high-flow reference to a low-flow reference. A ratio of 2 indicates modest variability (stable, groundwater-fed stream). A ratio of 100 or more indicates extreme flashiness (ephemeral or intermittent stream). Mediterranean-climate streams in Southern California commonly exceed 1,000.

8. Assignment: Analyse a Station of Your Choice

Step 1 — Confirm the default site works

Knit this document as provided and confirm you get a complete HTML output with all four plots and the summary table. The Ventura River example has a continuous record and is a reliable test case. Fix any errors before proceeding.

Step 2 — Find a USGS station

Option A — USGS WaterWatch map: Go to waterwatch.usgs.gov, zoom to your state, click a gauge dot, and note the 8-digit site number from the popup.

Option B — browse stations in R:

# Run this in the Console to browse sites in any state.
# Replace "CA" with your two-letter state abbreviation.
my_sites <- whatNWISsites(
  stateCd     = "CA",
  parameterCd = "00060",
  service     = "dv"
)
View(my_sites)   # opens a spreadsheet-style viewer in RStudio

Choose a site with:

Parameter 00060 (discharge) — not just stage height
At least 2 complete calendar years of approved (A) daily records
No major gaps — check the period of record shown in the popup

Step 3 — Update the Settings chunk

Return to Section 2 and replace the three values in the settings chunk:

site_no  <- "YOUR_8_DIGIT_SITE_NUMBER"
start_dt <- "YYYY-MM-DD"
end_dt   <- "YYYY-MM-DD"

Re-knit. All titles, labels, and the metadata line update automatically.

Step 4 — Write up your results

Add a new section below this one in your .Rmd and address these questions:

What river, state, and USGS site did you choose? What drew you to it — proximity, environmental significance, or something else?
What months have the highest and lowest median discharge? What physical mechanism drives the seasonality at your site — winter storms, snowmelt, hurricane season, monsoon?
Describe the shape of your flow duration curve compared to the Ventura River default. Is your stream more or less flashy? Compute the Q10/Q90 flashiness ratio for both sites and compare.
How many NA values did you remove during cleaning? What might cause data gaps at a USGS gauge?
Optional challenge: Use whatNWISdata() to find whether your site also has water-quality records (e.g., parameter 00010 = water temperature). If so, download a year of temperature data and add a fifth plot showing temperature alongside discharge.

Summary

💡 Key takeaways

readNWISdv() gives you programmatic access to the entire USGS daily discharge archive. The eight-digit site number is all you need to pull any station’s record.
Always inspect NAs after download. Equipment failures, rating-curve extrapolations, and provisional records all produce NAs in NWIS data. Document how many you removed and why.
Log scale is standard for discharge. The orders-of-magnitude range between base flow and flood peaks makes linear scale compress the low-flow period almost to invisibility.
The flow duration curve is the most information-dense single plot in hydrology. One curve encodes the full statistical distribution of daily discharge — base flow, seasonal variability, and flood frequency — in a single panel.
Mediterranean-climate streams are inherently flashy. The combination of seasonal rainfall concentrated in 3–4 winter months, steep terrain, and low-permeability soils produces Q10/Q90 ratios that are orders of magnitude larger than humid-continental or Pacific Northwest streams.

💬 Discussion questions

The Ventura River nearly ceases to flow in late summer (discharge approaches zero). What are the ecological consequences of intermittent flow for aquatic species — specifically steelhead trout, which are present in this watershed?
Compare the annual summary table across years. If mean annual discharge increased from one year to the next, does that necessarily mean it rained more? What other watershed changes could produce the same signal?
The FDC is constructed from the entire period of record pooled together. If you built separate FDCs for each year, what would a rightward shift in the high-exceedance tail (the dry end) tell you about long-term drought trends?
USGS converts stage height (feet) to discharge (cfs) using a site-specific rating curve. What events could invalidate a rating curve and require USGS to recalibrate it? How would you detect this in the data?
Southern California streams are heavily managed: Ventura River has a dam (Matilija Dam, currently being removed), diversions for agriculture, and groundwater pumping in its lower reach. How would you design a statistical test to detect the hydrological signature of dam removal in the discharge record?

USGS Stream Discharge Analysis

Marc Los Huertos

EA030 · Compiled 2026-04-24

USGS Stream Discharge: A Reproducible Hydrology Workflow

🎯 Learning Outcomes

The Analytical Pipeline

1. Packages

2. Settings — The Only Lines You Need to Change

3. Data Download

3.1 Query the NWIS API

3.2 Alternative: load a pre-downloaded CSV

4. Data Cleaning

5. Summary Statistics

6. Standard Hydrological Plots

6.1 Daily time series — linear scale

6.2 Daily time series — log₁₀ scale

6.3 Monthly discharge boxplots

7. Flow Duration Curve

8. Assignment: Analyse a Station of Your Choice

Step 1 — Confirm the default site works

Step 2 — Find a USGS station

Step 3 — Update the Settings chunk

Step 4 — Write up your results

Summary

💡 Key takeaways

💬 Discussion questions

References