Forecast lead-time example using real NOAA water level forecasts

[ecomodeller]

Draft — the data and example work end to end, but the caveats below want a decision before this is merge-ready.

Summary

A worked example of assessing forecast skill as a function of lead time, built on real forecasts. Using real data paid off immediately: it exposed two verification faults that synthetic data would have hidden, and both were caught by the reference forecast rather than the model.

The harmonic tide prediction cannot depend on lead time — it is not initialised — so a non-flat tide curve is proof the verification is wrong. It sagged, which led to:

1. The verification window slides with lead time. A 48 h forecast can only verify against times at least 48 h after the first model run, so every lead time is scored over a different period. The later days here happen to be easier, so long lead times looked flattered. Restricting to valid times shared by all lead times makes the tide exactly flat at 0.260 m, and the forecast degrades more steeply than first measured: 0.106 → 0.132 m.
2. Both series carry a constant datum offset — about -0.08 m for the forecast, -0.26 m for the tide. remove_bias() removes it, and this reverses the ranking: tide sits flat at 0.047 m against the forecast's 0.092–0.104 m. Lewes is tide-dominated and this week has no surge, which is exactly when harmonic prediction is hard to beat. An earlier version of this PR claimed the forecast beat the tide table at every horizon; that was the datum offset, not skill.

What survives is the lead-time signal, now clean: forecast error grows 0.092 → 0.104 m with residual bias drifting +0.020 → -0.005 m, while the reference stays flat.

Data — NOAA Delaware Bay OFS forecasts paired with observations and harmonic tide predictions from CO-OPS station 8557380 (Lewes, DE). 28 cycles over one week, lead times 0–48 h hourly, 1372 rows, 83 KB. Both sources are US federal works in the public domain.

The OFS station files are what make this cheap: 6.8 MB per cycle containing time series at 64 predefined stations rather than the full model grid, so no OPeNDAP or grid subsetting is involved. build_dataset.py regenerates the CSV, declares its dependencies inline (PEP 723), and reproduces the committed file byte for byte.

No new API required. Lead time rides along as an auxiliary variable and skill(by=["model", "lead_time"]) does the grouping. Worth noting for the lead-time roadmap item: overlapping forecast cycles mean 1161 of 1372 timestamps are duplicated, and from_matched() accepts that without complaint, as do .plot.timeseries() and .plot.scatter().

Suggested next step

This example shows the gap is ergonomics, not capability. Producing the central artefact of the whole feature — the skill-vs-lead-time curve — currently means creating the axes by hand and labelling them manually:

_, ax = plt.subplots()
sk.rmse.plot.line(ax=ax)
ax.set_xlabel("Lead time [hours]")
ax.set_ylabel("RMSE [m]")

The axis is also plotted as a category rather than a numeric axis, so all 49 tick labels render and overlap into an unreadable band until they are thinned by hand — visible in the first render of this page.

ModelSkill already knows the grouping variable is a lead time, that it is numeric and in hours, and what the quantity and its unit are, so the user shouldn't have to restate any of it. The suggested next step is a dedicated accessor owning this, following the pattern Comparer.vertical already uses for depth — data stays 1-D along time, a non-dimension coordinate carries the second axis, and the accessor provides the binned skill and the characteristic plot:

cmp.forecast.skill(bins=[0, 12, 24, 36, 48])   # binned skill by horizon
cmp.forecast.degradation(metric="rmse")        # labelled skill-vs-lead-time curve

Two design questions to settle before implementing:

• Alignment — when verifying at a long horizon, is the set of forecast reference times held fixed, or the set of valid times? The two give different skill curves. climpred makes this an explicit user choice and is worth reading first.
• Naming — CF uses forecast_reference_time and forecast_period, while lead_time is more widely understood and is what this example and the existing prototype notebooks already use.

A lead-time-aware persistence baseline depends on this, since a persistence forecast is a different series at every horizon.

Open questions for review

1. The model station is identified by position, not name. DBOFS station files carry only lon_rho/lat_rho. Index 5 was chosen as nearest to the Lewes gauge, ~1.4 km away. Someone should confirm this against a published DBOFS station list before we publish skill numbers attributed to Lewes.
2. Datum handling — the example now corrects the offsets with remove_bias() and explains why. The committed CSV keeps the raw values. Reviewers should sanity-check that the -0.26 m tide offset really is a datum/epoch mismatch and not a mistake in how the CO-OPS prediction product was requested.
3. No storm in the window. Ordinary tides only (-0.48 to 1.12 m) — and this now matters more than it did, since it is why the tide table wins. NCEI archives OFS output back to 2014, so a surge event is available at the cost of more work, and would let the forecast show actual value over harmonic prediction. I'd argue this is the strongest reason to extend the dataset before merging.
4. Should the synthetic tests/testdata/forecast_skill/ data and its two notebooks be retired once this lands? They cover the same ground with unattributed values, and nothing documents where those numbers came from.

Related: #692 corrects the roadmap entries this touches.

[Copilot]

Pull request overview

Adds a new end-to-end documentation example demonstrating how to assess forecast skill as a function of lead time using real NOAA DBOFS water level forecasts paired with Lewes, DE observations and a harmonic tide reference, and includes the accompanying reproducible dataset and build script.

Changes:

• Add the Lewes forecast/observation/tide dataset (CSV) plus provenance/limitations documentation and a standalone rebuild script.
• Introduce a new Quarto example page showing from_matched(..., aux_items=["lead_time"]) and skill(by=["model", "lead_time"]) workflows, including binning and horizon selection.
• Wire the new example into the docs navigation.

Reviewed changes

Copilot reviewed 5 out of 5 changed files in this pull request and generated 2 comments.

Show a summary per file

File	Description
tests/testdata/forecast_lewes/README.md	Documents dataset schema, sources, and known limitations for the Lewes lead-time example data.
tests/testdata/forecast_lewes/lewes_dbofs_forecast.csv	Adds the committed long-format forecast/obs/tide dataset used by the example.
tests/testdata/forecast_lewes/build_dataset.py	Standalone (PEP 723) script to regenerate the committed CSV from NOAA sources.
docs/examples/Forecast_lead_time.qmd	New worked example demonstrating lead-time skill grouping, binning, and filtering using auxiliary variables.
docs/_quarto.yml	Adds the new example to the documentation sidebar/navigation.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

Pull request overview

Copilot reviewed 5 out of 5 changed files in this pull request and generated no new comments.

Suppressed comments (1)

docs/examples/Forecast_lead_time.qmd:91

• pd.cut uses right-closed bins by default, so with edges [-1, 12, 24, 36, 48] the second bin is (12, 24] (i.e., 13–24), etc. The current labels (e.g., "12-24h") don’t match the actual bin membership at the boundaries (12, 24, 36). Updating the labels avoids misinterpretation of the binned skill table.

    horizon=pd.cut(
        df.lead_time, [-1, 12, 24, 36, 48], labels=["0-12h", "12-24h", "24-36h", "36-48h"]
    )

[Copilot]

Pull request overview

Copilot reviewed 5 out of 5 changed files in this pull request and generated no new comments.

Suppressed comments (2)

docs/examples/Forecast_lead_time.qmd:84

• groupby(...).apply(..., include_groups=False) is not available in the project’s supported pandas range (pyproject.toml allows pandas>=1.4). This example will fail on older supported pandas versions; compute the daily RMSE with a version-agnostic groupby instead.

unique_times = df.drop_duplicates("valid_time")
unique_times.assign(day=unique_times.valid_time.dt.date).groupby("day").apply(
    lambda d: ((d.tide - d["observed"]) ** 2).mean() ** 0.5, include_groups=False
).round(3)

tests/testdata/forecast_lewes/build_dataset.py:3

• The PEP 723 header declares requires-python = ">=3.10", but the package itself requires Python >=3.12 (pyproject.toml). Using a lower requirement here is misleading when running via uv in this repo; align it with the project requirement.

# /// script
# requires-python = ">=3.10"
# dependencies = ["xarray", "netCDF4", "pandas"]