How to fit return periods to annual maxima in Python

Problem statement

Fitting a flood frequency curve is a short piece of code wrapped in several decisions that change the answer: which distribution, which fitting method, whether to include historic peaks, and how to express the uncertainty. Get the code right and the curve still means nothing without those decisions written down beside it.

This guide fits annual peak flows for a long-record gauge, compares three distributions, adds a bootstrap confidence band, and produces a table that states its own assumptions. The gauge is the Delaware River at Trenton, 121 years of peaks, which is a long record by any standard — and the 1-in-100 estimate still spans 239,000 to 319,000 cfs between its 10th and 90th bootstrap percentile — 29% of the estimate.

Quick answer

Get annual maxima, fit, and report return levels with an interval:

import io, urllib.request
import numpy as np, pandas as pd
from scipy.stats import genextreme

def annual_peaks(site):
    url = (f"https://nwis.waterdata.usgs.gov/nwis/peak?site_no={site}"
           "&agency_cd=USGS&format=rdb")
    req = urllib.request.Request(url, headers={"User-Agent": "example/1.0"})
    text = urllib.request.urlopen(req, timeout=90).read().decode()
    rows = [ln for ln in text.splitlines() if not ln.startswith("#")]
    df = pd.read_csv(io.StringIO("\n".join(rows)), sep="\t").iloc[1:]
    df["q"] = pd.to_numeric(df.peak_va, errors="coerce")
    df["year"] = pd.to_datetime(df.peak_dt, errors="coerce", format="mixed").dt.year
    return df.dropna(subset=["q", "year"])[["year", "q", "peak_cd"]]

peaks = annual_peaks("01463500")                       # Delaware River at Trenton, NJ
q = peaks.q.to_numpy()

def fit_gev(x):
    return genextreme.fit(x, 0.1, loc=np.median(x), scale=x.std())

params = fit_gev(q)
def level(T, p=params):
    return genextreme.ppf(1 - 1 / T, *p)

print(f"{len(q)} annual peaks, {int(peaks.year.min())}–{int(peaks.year.max())}")
for T in (2, 10, 50, 100, 500):
    print(f"  1-in-{T:<4} {level(T):9,.0f} cfs")

121 annual peaks, 1897–2023
  1-in-2       92,855 cfs
  1-in-10     163,347 cfs
  1-in-50     240,980 cfs
  1-in-100    279,127 cfs
  1-in-500    381,646 cfs

The starting guess in fit_gev is not decoration: without it, scipy's optimiser can settle on a shape parameter that puts the 1-in-100 flood in the trillions.

Annual maxima, qualification codes, a chosen distribution

1. Get annual maxima, not all peaks

One value per water year. The USGS peak-flow service already returns annual maxima; other sources may need resampling from a continuous series.

2. Read the qualification codes

peak_cd flags peaks that are estimated, affected by dam regulation, or the maximum of a historic period rather than a gauged year. Codes 5, 6 and 7 mark regulated or urbanised conditions — a record that changes regime partway through is not one sample.

3. Choose a distribution deliberately

GEV, Gumbel and log-Pearson III are the common choices. Fit more than one and look at how far apart they are in the tail you care about.

4. Bootstrap the uncertainty

Resample the years with replacement, refit, and take percentiles of the return level. This captures sampling uncertainty, though not the uncertainty in the choice of distribution.

5. Check the fit against the record

Plot the fitted curve against plotting positions of the observed peaks. A systematic bend means the distribution is the wrong shape for this river.

6. Write the assumptions into the output

Site, record length, distribution, fitting method, and whether the record is regulated. A table of return levels without them is not reproducible.

Three distributions with bootstrap bands

Example 1 — the qualification codes

print("peaks by qualification code:")
print(peaks.peak_cd.fillna("(none)").value_counts().head(6).to_string())
regulated = peaks.peak_cd.fillna("").str.contains("5|6")
print(f"flagged as regulated or affected: {regulated.sum()} of {len(peaks)}")

peaks by qualification code:
peak_cd
5         69
(none)    48
2          3
2,5        1
flagged as regulated or affected: 70 of 121

Example 2 — three distributions and a bootstrap band

from scipy.stats import gumbel_r, pearson3

def levels(sample, T):
    g = genextreme.ppf(1 - 1 / T, *fit_gev(sample))
    gu = gumbel_r.ppf(1 - 1 / T, *gumbel_r.fit(sample))
    lp = 10 ** pearson3.ppf(1 - 1 / T, *pearson3.fit(np.log10(sample)))
    return g, gu, lp

g, gu, lp = levels(q, 100)
print(f"1-in-100: GEV {g:,.0f} | Gumbel {gu:,.0f} | log-Pearson III {lp:,.0f} cfs")

rng = np.random.default_rng(1)
boot = np.array([genextreme.ppf(0.99, *fit_gev(rng.choice(q, len(q), replace=True)))
                 for _ in range(500)])
lo, hi = np.percentile(boot, [10, 90])
print(f"GEV 1-in-100 bootstrap 10th–90th percentile: {lo:,.0f}–{hi:,.0f} cfs "
      f"({(hi - lo) / g:.0%} of the estimate)")

1-in-100: GEV 279,127 | Gumbel 236,587 | log-Pearson III 272,413 cfs
GEV 1-in-100 bootstrap 10th–90th percentile: 238,736–319,056 cfs (29% of the estimate)

Example 3 — does the fit follow the record?

ranked = np.sort(q)[::-1]
n = len(ranked)
weibull_T = (n + 1) / np.arange(1, n + 1)              # plotting positions
for i in (0, 1, 4, 9, 24, 59):
    T = weibull_T[i]
    print(f"observed rank {i + 1:3d}: {ranked[i]:8,.0f} cfs at T={T:6.1f} yr; "
          f"GEV says {level(T):8,.0f} cfs")

observed rank   1:  329,000 cfs at T= 122.0 yr; GEV says  290,690 cfs
observed rank   2:  295,000 cfs at T=  61.0 yr; GEV says  251,593 cfs
observed rank   5:  227,000 cfs at T=  24.4 yr; GEV says  204,699 cfs
observed rank  10:  161,000 cfs at T=  12.2 yr; GEV says  172,262 cfs
observed rank  25:  133,000 cfs at T=   4.9 yr; GEV says  132,152 cfs
observed rank  60:   90,700 cfs at T=   2.0 yr; GEV says   93,641 cfs

The middle of the record is followed closely — at T = 2 and T = 5 the curve is within a few per cent. The two largest peaks sit 13% and 17% above it, which is inside the bootstrap band but hints that the real tail is heavier than this fit. A curve that misses the top of the record in the same direction every time is the wrong distribution for the river.

The tail is where the choice matters

Why the distribution choice matters most in the tail

All three distributions are fitted to the same mean and spread, so they agree near the middle of the record. They differ in how fast the tail decays: Gumbel has a fixed light tail, GEV lets the data choose, and log-Pearson III works in log space with a skew parameter. Those differences are invisible at T = 10 and decide the answer at T = 500.

Why bootstrap intervals are the minimum

The fitted curve looks like a smooth, confident line. Resampling shows what a different 121 years from the same process could have produced: a band tens of thousands of cubic feet per second wide at T = 100. That band is sampling uncertainty only — the true uncertainty also includes the model choice and any change in the catchment.

Why regulation breaks the model

Extreme-value fitting assumes every year comes from the same distribution. A dam completed mid-record changes the distribution of peaks below it. The clean approach is to fit only the homogeneous part of the record, which trades sample size for validity; the honest minimum is to say which years were flagged.

Water years, historic peaks and regional skew

• Water years start in October in the US; calendar years split winter floods across two rows.
• Historic peaks before the gauge are censored data and need a fitting method that handles them.
• Regional skew is used in Bulletin 17C to stabilise log-Pearson III fits on short records.
• Partial-duration series give different low-T answers and the same high-T ones.
• Zero-flow years in arid catchments need a conditional model.
• Rating-curve extrapolation makes the largest peaks themselves uncertain; the flow is an estimate, not a measurement.

Internal links

• Return periods and annual exceedance probability explained — what the numbers mean
• A 100-year flood happened twice in ten years — the fix guide
• Flood hazard maps explained — where return levels become maps
• Hazard, exposure and vulnerability explained — the wider frame
• How to estimate b-value and completeness in Python — the seismic analogue
• Margins of error explained — reporting uncertainty

FAQ

How do I compute a 100-year flood in Python?

Fit an extreme-value distribution to annual maxima and take the quantile at 1 − 1/100.

Which distribution should I use?

GEV or Gumbel generally, log-Pearson III for US practice. Fit more than one and compare in the tail.

How do I get a confidence interval?

Bootstrap: resample years with replacement, refit, and take percentiles of the return level.

What do the USGS peak qualification codes mean?

They flag estimated, regulated and historic peaks; codes 5 and 6 mark regulation or urbanisation.

Do I need water years?

For US data, yes — calendar years can split a single flood season.

Why does my scipy fit produce absurd numbers?

The optimiser found a degenerate shape parameter. Give it a starting guess and sanity-check against the record.