# Handoff: porting the dmdtda/mass-loss diagnostic into OGGM proper

Written for a fresh Claude session (or human dev) with no context on this conversation. Explains exactly
how `plot_mass_loss_by_region_*.png` and `plot_dmdtda_by_region_*.png` were built in the `glambie` project
(`plot_glambie_diagnostics.py`), so the idea can be reimplemented as a general OGGM diagnostic tool rather
than a GLAMBIE-submission-specific script.

## What the two plots are

Both compare an OGGM regional MB estimate against the Hugonnet et al. (2021) geodetic reference over the
standard 2000–2019 calibration period, per RGI region plus a global aggregate:

- **`plot_mass_loss_by_region`**: absolute mass loss rate, **Gt yr⁻¹** (log-scale y-axis, since regions span
  ~2 to ~2 orders of magnitude in size). Bars = OGGM; black error-barred dots = Hugonnet.
- **`plot_dmdtda_by_region`**: **specific** MB rate, **m w.e. yr⁻¹** (linear y-axis, values cluster in a
  narrow, comparable range regardless of region size — this is the more diagnostic of the two for spotting
  a model bias, since it isn't dominated visually by the biggest regions the way the Gt plot is). Bars =
  OGGM; black error-barred dots = Hugonnet regional (with area extrapolation to unmeasured glaciers); grey
  square markers = Hugonnet **per-glacier** data re-aggregated with Hugonnet's own area weights and **no**
  extrapolation — i.e. the same glacier population and aggregation method OGGM uses, for an apples-to-apples
  comparison independent of Hugonnet's own regional extrapolation choices.

## Inputs

In this project, the input is a GLAMBIE-format submission CSV: one row per (region, month), columns
`region_id`, `start_date` (DD/MM/YYYY), `glacier_change_observed` (specific MB, m w.e., **that one month**),
`glacier_area_reference_start` (km², the area backing that row — constant per region for a fixed-geometry
submission, but genuinely time-varying for a dynamic/spinup-geometry one).

**This CSV is GLAMBIE-specific glue, not something OGGM proper would have.** To generalize this as an OGGM
diagnostic, the natural input instead would be something already inside OGGM's own regional-compilation
pipeline — e.g. per-glacier `specific_mb` time series (from `compile_glacier_statistics` /
`compile_run_output`-style regional netcdfs, or directly per-gdir output) plus each glacier's/region's area.
The computation below only actually needs, per RGI region and per month or year: (a) monthly-or-annual
specific MB in m w.e., and (b) the area (km²) that specific MB is relative to, per timestep.

## Step-by-step computation

```python
# 1. Load and restrict to the standard reference period
mask = df["time"].dt.year.between(2000, 2019)
df_hp = df[mask]

# 2. Area for the Gt conversion: MEAN area over the comparison window itself,
#    not a single fixed value or the first row. This matters as soon as area
#    can vary in time (e.g. a dynamically-evolving-geometry run) - using a
#    stale/first-row area silently mis-scales the Gt conversion. This was a
#    real bug caught and fixed in this project: originally used
#    `.groupby("region_id")[...].first()` (the very first row of the whole
#    df, i.e. the 1975 area), which is invisible/harmless for constant-area
#    runs but wrong for a varying-geometry run evaluated at 2000-2019.
area = df_hp.groupby("region_id")["glacier_area_reference_start"].mean()  # km²

# 3. dmdtda: mean annual specific MB over the 20-year window (m w.e. yr⁻¹)
#    Sum 12 monthly values -> one annual number per (region, year), then
#    average those 20 annual numbers per region.
dmdtda_our = (
    df_hp.groupby(["region_id", df_hp["time"].dt.year])["glacier_change_observed"]
    .sum()
    .groupby(level=0)
    .mean()
)

# 4. dmdt: annual mass loss rate (Gt yr⁻¹)
#    1 m w.e. over 1 km² = 1e6 m³ water = 1e9 kg = 1e-3 Gt, hence the 1e-3.
dmdt_our = dmdtda_our * area * 1e-3

# 5. Global aggregate: area-weighted mean for the specific-rate plot,
#    plain sum for the mass-rate plot (Gt are already additive; m w.e. rates
#    are not - they must be weighted by the area they apply to).
total_area = area.sum()
dmdtda_global_our = (dmdtda_our * area).sum() / total_area
dmdt_global_our = dmdt_our.sum()
```

**Sign convention:** OGGM/GLAMBIE store MB as positive-accumulation, negative-loss (standard glaciology
convention). Hugonnet's `dmdt`/`dmdtda` columns use the same convention. Both are negated (`-dmdt_our`,
`-hug.loc[r, "dmdt"]`) only for **display**, so the bar/dot for a shrinking glacier points down — the sign
itself carries no other meaning and this negation is purely a plotting choice, not part of the physics.

## Fetching the Hugonnet et al. 2021 reference (already an OGGM utility — no porting needed)

This part is *already* general-purpose OGGM code (`oggm.utils.get_geodetic_mb_dataframe`), not something
built for this project:

```python
from oggm import utils

# Pre-aggregated regional values (with extrapolation to unmeasured glacier area, via Hugonnet's own is_cor
# corrected values). One row per (region, period); "period" is a string like "2000-01-01_2020-01-01".
# region 21 is the special "global" aggregate row Hugonnet provides.
hug_all = utils.get_geodetic_mb_dataframe(regional=True)
period = "2000-01-01_2020-01-01"
hug = hug_all[(hug_all["period"] == period) & (hug_all["reg"].isin(REGIONS))].copy()
hug.index = hug["reg"].astype(int)
hug_glob_row = hug_all[(hug_all["period"] == period) & (hug_all["reg"] == 21)].iloc[0]
# columns used: dmdt, err_dmdt, dmdtda, err_dmdtda

# Per-glacier values (NOT regionally extrapolated) - used to build the "grey square" series, an
# apples-to-apples comparison against OGGM's own glacier population/area weights.
gdf_all = utils.get_geodetic_mb_dataframe()  # regional=False (default) -> per-glacier parquet
gdf_p = gdf_all[gdf_all["period"] == period].copy()
gdf_p["reg"] = gdf_p["reg"].astype(int)
hug_pergla = gdf_p.groupby("reg").apply(
    lambda g: pd.Series({"dmdtda": (g["dmdtda"] * g["area"]).sum() / g["area"].sum()}),
    include_groups=False,
).reindex(REGIONS)
```

Both fetch from `https://cluster.klima.uni-bremen.de/~oggm/geodetic_ref_mb/` (one CSV for the regional
file, one parquet for the per-glacier file), cached via OGGM's standard `file_downloader` + an in-process
`cfg.DATA` cache keyed by file path — repeated calls in the same process are free. Nothing else to build
here; a generalized diagnostic can call this utility as-is.

## Plotting mechanics

Both are grouped bar charts, one bar/dot pair per RGI region plus a "Global" column at the end:

```python
labels = [f"RGI{r:02d}" for r in REGIONS] + ["Global"]
x = np.arange(len(labels))
bar_colors = [tab20_color_per_region] + ["#333333"]  # distinct color per region, dark grey for Global

fig, ax = plt.subplots(figsize=(13, 5))
ax.bar(x, our_values, color=bar_colors, alpha=0.8, zorder=2, label="OGGM ...")
ax.errorbar(x, hug_values, yerr=hug_err, fmt="o", color="k", ms=5, lw=1.2, capsize=3, zorder=3,
            label="Hugonnet et al. 2021 (...)")
# dmdtda plot additionally overlays the per-glacier grey squares and an axhline(0)
# mass-loss plot additionally uses ax.set_yscale("log")
ax.set_xticks(x); ax.set_xticklabels(labels, rotation=45, ha="right")
ax.grid(axis="y", alpha=0.3, zorder=1)
```

`zorder` matters: bars at 2, per-glacier markers at 3, the errorbar'd Hugonnet regional dots at 4 (topmost,
since that's the primary reference line readers compare against).

## What to change to make this a real OGGM diagnostic (not just a copy of this script)

1. **Input**: replace the GLAMBIE-CSV loader with whatever OGGM's own regional-compilation output looks like
   (per-glacier or per-region specific-MB time series + area per timestep). The downstream math (steps 1–5
   above) is unchanged — it only needs a `(region, year) -> specific_mb (m w.e.)` series and a
   `(region, year-or-comparison-window) -> area (km²)` series.
2. **Region list**: currently hardcoded to `range(1, 20)` (RGI 19 first-order regions) because that's what a
   GLAMBIE submission always contains; a general tool should derive `REGIONS` from whatever's actually in
   the input, and handle subsets (e.g. someone running just RGI11).
3. **Time window**: hardcoded to `2000-2019` here because that's Hugonnet's standard geodetic period; keep
   it as a parameter, since OGGM calibration workflows sometimes use other periods.
4. **The area-averaging bug** (step 2 above) generalizes as a rule worth keeping explicit in any OGGM
   version: **always average area over the exact comparison window, never take a value from outside it**,
   because OGGM increasingly supports non-constant-geometry runs (dynamic spinup, etc.) where a "first row"
   or "current" area shortcut silently produces a wrong Gt conversion with no error or warning.
5. Consider exposing both plots (and the underlying `dmdt_our`/`dmdtda_our`/`hug`/`hug_pergla` tables) from
   one function so a caller can get the numbers without being forced to also produce the plot, e.g.
   `oggm.utils.compare_to_geodetic_reference(mb_by_region, area_by_region, period=...) -> pd.DataFrame`,
   with `oggm.graphics.plot_geodetic_comparison(...)` as a thin plotting wrapper around it.

## Source

Full working code: `plot_glambie_diagnostics.py` in this repo (glambie project,
`/home/www/fmaussion/glambie/`). Run as `python plot_glambie_diagnostics.py <submission.csv> [suffix]`.