Forecast lead-time example using real NOAA water level forecasts

docs/_quarto.yml

@@ -66,6 +66,8 @@ website:
           text: Gridded NetCDF ModelResult
         - href: examples/Prematched_with_auxiliary.qmd
           text: Prematched with auxiliary
+        - href: examples/Forecast_lead_time.qmd
+          text: Forecast skill by lead time
         - href: examples/Skill_vs_dummy.qmd
           text: Compare with dummy results
         - href: examples/Metrics_custom_metric.qmd

docs/examples/Forecast_lead_time.qmd

@@ -0,0 +1,147 @@
+---
+title: Forecast skill by lead time
+jupyter: python3
+---
+
+A forecast issued 6 hours ahead is usually better than one issued 48 hours ahead. This example
+quantifies that degradation using real water level forecasts from the NOAA Delaware Bay
+Operational Forecast System (DBOFS) at Lewes, Delaware.
+
+A forecast dataset needs **two** time dimensions: the time the forecast is *for* (the valid time)
+and how far ahead it was issued (the lead time). ModelSkill carries the lead time as an auxiliary
+variable, which can then be used for grouping.
+
+```{python}
+import matplotlib.pyplot as plt
+import pandas as pd
+import modelskill as ms
+```
+
+## The data
+
+Each row is one forecast: issued at `reference_time`, valid at `valid_time`, at a horizon of
+`lead_time` hours. The model runs four times a day out to 48 hours, so successive runs overlap —
+the same valid time is forecast repeatedly, at a shorter lead time each run.
+
+```{python}
+df = pd.read_csv(
+    "../data/forecast_lewes/lewes_dbofs_forecast.csv",
+    parse_dates=["reference_time", "valid_time"],
+)
+df.head()
+```
+
+The `tide` column is the harmonic tide prediction for the same station. It is not initialised from
+recent conditions, so its error *must* be independent of lead time. That makes it a control: if
+the tide curve is not flat, the verification is at fault, not the tide.
+
+## Matching
+
+```{python}
+cmp = ms.from_matched(
+    df.set_index("valid_time").drop(columns="reference_time"),
+    obs_item="observed",
+    mod_items=["dbofs", "tide"],
+    aux_items=["lead_time"],
+    quantity=ms.Quantity("Water Level", "m"),
+)
+cmp
+```
+
+## A first look
+
+```{python}
+sk = cmp.skill(by=["model", "lead_time"], metrics=["bias", "rmse"])
+
+_, ax = plt.subplots()
+sk.rmse.plot.line(ax=ax)
+ax.set_xticks(range(0, 49, 6))  # lead time is plotted as a category, so thin the ticks
+ax.set_xlabel("Lead time [hours]")
+ax.set_ylabel("RMSE [m]")
+ax.set_title("Skill by lead time — before aligning verification times");
+```
+
+The forecast error grows with the horizon, as expected. But the tide reference *falls* from 0.26 m
+to 0.24 m, which is impossible — so the verification needs fixing first.
+
+## The verification window slides
+
+A forecast at lead time 48 h can only verify against times at least 48 h after the first model
+run, so each lead time is scored over a different period:
+
+```{python}
+for lead in [0, 24, 48]:
+    d = df[df.lead_time == lead]
+    print(f"lead {lead:2d} h: {d.valid_time.min().date()} to {d.valid_time.max().date()}")
+```
+
+And the days are not equally easy — the tide alone fits better late in this record:
+
+```{python}
+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)
+```
+
+Longer lead times slide into the easier days, so skill differences between horizons are confounded
+with differences between periods.
+
+## Holding the verification times fixed
+
+Keep lead times on the model's 6-hourly cycle spacing so all of them sample the same hours of the
+day, then restrict to the valid times they share.
+
+```{python}
+on_cycle = df[df.lead_time % 6 == 0]
+shared = set.intersection(*[set(g.valid_time) for _, g in on_cycle.groupby("lead_time")])
+aligned = on_cycle[on_cycle.valid_time.isin(shared)]
+
+cmp_aligned = ms.from_matched(
+    aligned.set_index("valid_time").drop(columns="reference_time"),
+    obs_item="observed",
+    mod_items=["dbofs", "tide"],
+    aux_items=["lead_time"],
+    quantity=ms.Quantity("Water Level", "m"),
+)
+cmp_aligned.skill(by=["model", "lead_time"], metrics=["bias", "rmse"]).round(3)
+```
+
+The tide is now flat at 0.260 m, and the forecast rises 0.106 → 0.132 m — a steeper degradation
+than the first plot showed. The cost is sample size: 20 verification times per lead instead of 28.
+
+## Removing the datum offset
+
+Both series carry a large constant bias, about -0.08 m for the forecast and -0.26 m for the tide.
+Errors that size at lead time 0 indicate reference levels, not forecast quality. `remove_bias()`
+subtracts each model's mean residual, leaving the time-varying error untouched.
+
+```{python}
+unbiased = cmp_aligned.remove_bias()
+sk_unbiased = unbiased.skill(by=["model", "lead_time"], metrics=["bias", "rmse"])
+
+_, ax = plt.subplots()
+sk_unbiased.rmse.plot.line(ax=ax)
+ax.set_xlabel("Lead time [hours]")
+ax.set_ylabel("RMSE [m]")
+ax.set_title("Skill by lead time — datum offsets removed");
+```
+
+This reverses the ranking. The tide sits flat at 0.047 m, well below the forecast's 0.092–0.104 m;
+the earlier impression that the forecast beat the tide table was the tide's datum offset, not
+skill. Lewes is strongly tide-dominated and this week has no storm surge, which is exactly when
+harmonic prediction is hard to beat.
+
+The lead-time signal survives and is cleaner: error grows 0.092 → 0.104 m, and residual bias
+drifts from +0.020 m to -0.005 m, while the reference stays flat.
+
+## Selecting a single horizon
+
+Lead time is an ordinary auxiliary variable, so `where()` isolates one horizon and everything else
+behaves normally.
+
+```{python}
+unbiased.where(unbiased.data["lead_time"] == 48).plot.scatter();
+```
+
+See `data/forecast_lewes/README.md` for provenance and known limitations.

tests/testdata/forecast_lewes/README.md

@@ -0,0 +1,63 @@
+# Forecast lead-time example data: Lewes, Delaware
+
+Water level forecasts at multiple lead times, paired with observations, for use in the
+forecast lead-time example. All data originate from NOAA and are in the public domain.
+
+## `lewes_dbofs_forecast.csv`
+
+Long format, one row per (forecast cycle, valid time):
+
+| column | description |
+| --- | --- |
+| `reference_time` | when the forecast was issued (the model cycle, 4 per day) |
+| `valid_time` | the time the forecast is *for* |
+| `lead_time` | forecast horizon in hours, `valid_time - reference_time`, 0–48 |
+| `dbofs` | forecast water level from NOAA Delaware Bay OFS \[m] |
+| `observed` | observed water level, CO-OPS station 8557380 \[m, MSL] |
+| `tide` | harmonic tide prediction for the same station \[m, MSL] — a reference forecast with no initialisation, so its error does not grow with lead time |
+
+28 cycles covering 2026-07-26 to 2026-08-01, hourly, 1372 rows.
+
+## Sources
+
+- **Model**: NOAA Delaware Bay Operational Forecast System (DBOFS) station files, from the
+  [NOAA OFS data on AWS](https://registry.opendata.aws/noaa-ofs-pds/)
+  (`s3://noaa-ofs-pds/dbofs.<date>/dbofs.t<cycle>z.<date>.stations.forecast.nc`).
+  The station files contain time series at 64 predefined stations rather than the full model
+  grid, which is why they are small enough to subset quickly.
+  The AWS bucket holds only a trailing 30-day window; longer archives (back to 2014) are at
+  [NCEI](https://www.ncei.noaa.gov/products/weather-climate-models/co-ops-operational-forecast).
+- **Observations and tide predictions**: NOAA CO-OPS station 8557380 (Lewes, DE) via the
+  [CO-OPS data API](https://api.tidesandcurrents.noaa.gov/api/prod/).
+
+Both are works of the US federal government and are in the public domain.
+
+## Known limitations
+
+These are real caveats, not artefacts of the export — an example built on this data should
+acknowledge them rather than paper over them.
+
+- **The model station is identified by position, not by name.** DBOFS station files carry only
+  `lon_rho`/`lat_rho`, with no station identifiers. Station index 5 was selected as the nearest
+  to the Lewes gauge, about 1.4 km away. This has not been verified against a published DBOFS
+  station list.
+- **Datums are not reconciled in the data.** Model `zeta` is free-surface relative to the model's
+  own reference level; observations are relative to the station MSL datum. Both model columns
+  carry a large constant bias as a result — about -0.08 m for `dbofs` and -0.26 m for `tide`. The
+  example corrects this with `remove_bias()`, which matters: with the offsets left in, the
+  tide-only prediction looks worse than the forecast, and with them removed it is clearly better.
+- **No storm event.** Observed water level spans -0.48 to 1.12 m over this week — ordinary
+  tidal conditions. A surge event would show the value of the forecast over the tide-only
+  reference far more clearly.
+- **The 48 h lead time has only 28 points** (one per cycle), so the end of any skill-vs-lead-time
+  curve is noisier than the rest.
+
+## Regenerating
+
+`build_dataset.py` re-downloads and rebuilds the CSV. It declares its own dependencies inline
+(PEP 723), so it runs standalone. It requires network access, and will only work for dates still
+inside the AWS 30-day window — the committed CSV is the durable artefact.
+
+```bash
+uv run tests/testdata/forecast_lewes/build_dataset.py
+```

tests/testdata/forecast_lewes/build_dataset.py

@@ -0,0 +1,108 @@
+# /// script
+# requires-python = ">=3.10"
+# dependencies = ["xarray", "netCDF4", "pandas"]
+# ///
+"""Rebuild lewes_dbofs_forecast.csv from NOAA sources.
+
+Pairs NOAA Delaware Bay OFS (DBOFS) water level forecasts at multiple lead times with
+observations from CO-OPS station 8557380 (Lewes, DE), plus harmonic tide predictions as a
+reference forecast.
+
+Requires network access. The AWS bucket keeps only a trailing 30-day window, so DATES must
+be recent; see README.md for the NCEI archive if you need an older period.
+
+    uv run build_dataset.py
+"""
+
+from __future__ import annotations
+
+import tempfile
+from pathlib import Path
+from urllib.request import urlretrieve
+
+import pandas as pd
+import xarray as xr
+
+DATES = [f"202607{d:02d}" for d in range(26, 32)] + ["20260801"]
+CYCLES = ["00", "06", "12", "18"]
+
+STATION_ID = "8557380"  # CO-OPS Lewes, DE
+STATION_IX = 5  # nearest DBOFS station; see "Known limitations" in README.md
+
+S3 = "https://noaa-ofs-pds.s3.amazonaws.com"
+API = "https://api.tidesandcurrents.noaa.gov/api/prod/datagetter"
+
+OUT = Path(__file__).parent / "lewes_dbofs_forecast.csv"
+
+
+def read_cycle(date: str, cycle: str) -> pd.DataFrame:
+    """Read one forecast cycle and return its hourly water level by lead time."""
+    url = f"{S3}/dbofs.{date}/dbofs.t{cycle}z.{date}.stations.forecast.nc"
+    with tempfile.NamedTemporaryFile(suffix=".nc") as tmp:
+        urlretrieve(url, tmp.name)  # noqa: S310 - fixed https host
+        with xr.open_dataset(tmp.name, decode_timedelta=False) as ds:
+            time = pd.DatetimeIndex(ds.ocean_time.values)
+            zeta = pd.Series(ds.zeta.isel(station=STATION_IX).values, index=time)
+
+    zeta = zeta[zeta.index.minute == 0]  # 6-min output -> hourly
+    reference_time = time[0]
+
+    return pd.DataFrame(
+        {
+            "reference_time": reference_time,
+            "valid_time": zeta.index,
+            "lead_time": ((zeta.index - reference_time).total_seconds() / 3600).astype(
+                int
+            ),
+            "dbofs": zeta.values,
+        }
+    )
+
+
+def read_coops(product: str, column: str, **extra: str) -> pd.DataFrame:
+    """Read a CO-OPS product for the observation station."""
+    params = {
+        "product": product,
+        "application": "modelskill",
+        "begin_date": DATES[0],
+        "end_date": "20260804",
+        "datum": "MSL",
+        "station": STATION_ID,
+        "time_zone": "gmt",
+        "units": "metric",
+        "format": "csv",
+        **extra,
+    }
+    url = API + "?" + "&".join(f"{k}={v}" for k, v in params.items())
+    df = pd.read_csv(url, parse_dates=["Date Time"])
+    df = df.rename(columns={c: c.strip() for c in df.columns})
+    return df.rename(columns={"Date Time": "valid_time"})[["valid_time", column]]
+
+
+def main() -> None:
+    forecasts = pd.concat(
+        [read_cycle(date, cycle) for date in DATES for cycle in CYCLES],
+        ignore_index=True,
+    )
+
+    observed = read_coops("water_level", "Water Level").rename(
+        columns={"Water Level": "observed"}
+    )
+    observed = observed[observed.valid_time.dt.minute == 0]
+    tide = read_coops("predictions", "Prediction", interval="h").rename(
+        columns={"Prediction": "tide"}
+    )
+
+    df = forecasts.merge(observed, on="valid_time").merge(tide, on="valid_time")
+    df = df.sort_values(["reference_time", "lead_time"]).reset_index(drop=True)
+    df[["dbofs", "observed", "tide"]] = df[["dbofs", "observed", "tide"]].round(3)
+
+    df.to_csv(OUT, index=False)
+    print(
+        f"wrote {OUT} - {len(df)} rows, {df.reference_time.nunique()} cycles, "
+        f"lead times {df.lead_time.min()}-{df.lead_time.max()} h"
+    )
+
+
+if __name__ == "__main__":
+    main()