# --- # jupyter: # jupytext: # text_representation: # extension: .py # format_name: percent # format_version: '1.3' # --- # %% [markdown] # # Simple attribution: one HIST member against the NAT ensemble # # After Mira's SNR method (email 2026-10), on the attribution dataset # `round4/ensemble/output/lam8_jjas_dyn` (see the README there). # # For a quantity X (volume, area, length) and a time t: # # dX(t) = X(t) - X(baseline) baseline = mean 1820-1850 (Mira's PI) # # where X(t) is the mean over the WINDOW years ending at t (Mira's "modern" # period = the last 15 years). Then, against the NAT members that reach t: # # signal = dX_hist - median(dX_nat) # noise = std(dX_nat) (spread across NAT members) # SNR = -signal / noise (> 0: HIST lost MORE than NAT) # p_nat = fraction of NAT members with dX_nat <= dX_hist # (how often natural forcing alone gives a loss this large) # # Differences from Mira's snippet: # - the SNR is signed (losses count positive), not abs(); # - the noise is the spread across the NAT ensemble of the same change # dX, not the year-to-year std of one NAT run over years 1-1600 (which # includes the model spin-up and is a different quantity); # - NAT and HIST use the same calendar window. NAT members end in 2006-2021, # so after 2006 only the members that reach t are used (127 at 2021), # instead of "the last 15 years of each run". # # Part 1: regional sums (the CSVs of preprocess_mira.py). # Part 2: per glacier (Mira's histogram figure), at several times. # # Run from round4/attribution/ with the oggm_env python. Env vars: EXP, # MEMBER, WINDOW, REGIONAL_ONLY (=1: stop after part 1, no per-glacier cube). # %% import os import numpy as np import pandas as pd import xarray as xr import matplotlib.pyplot as plt from scipy.stats import norm EXP = os.environ.get('EXP', 'lam8_jjas_dyn') MEMBER = os.environ.get('MEMBER', 'ens183') WINDOW = int(os.environ.get('WINDOW', 15)) BASELINE = (1820, 1850) T_END = 2021 # last year of the longest NAT runs TIMES = [1875, 1900, 1925, 1950, 1975, 2000, 2021] # snapshot times SNR_THRESH = [(1, 'Likely (SNR>1)'), (1.6, 'Very likely (SNR>1.6)'), (2, 'Extremely likely (SNR>2)')] BASE = os.path.dirname(os.path.abspath(__file__)) ENS = os.path.join(BASE, '..', 'ensemble') FILTERED = os.path.join(ENS, 'filtered', EXP) OUTPUT = os.path.join(ENS, 'output', EXP) FIGS = os.path.join(BASE, 'figs') os.makedirs(FIGS, exist_ok=True) TAG = f'{EXP}_{MEMBER}_w{WINDOW}' HIST_C, NAT_C, NAT_DARK = '#c0392b', '#2e86c1', '#1b4f72' def delta(df): """X(t) as a trailing WINDOW-year mean, minus the baseline mean. df: (year, ...) DataFrame or (year, member, glacier) DataArray. A window with a missing year is NaN (no partial windows at run ends). """ if isinstance(df, pd.DataFrame): base = df.loc[BASELINE[0]:BASELINE[1]].mean() return df.rolling(WINDOW, min_periods=WINDOW).mean() - base base = df.sel(time=slice(*BASELINE)).mean('time') return df.rolling(time=WINDOW, min_periods=WINDOW).mean() - base def stats(d_hist, d_nat): """Attribution numbers for one value of dX_hist and an array of dX_nat.""" d_nat = d_nat[np.isfinite(d_nat)] med, sd = np.median(d_nat), np.std(d_nat, ddof=1) snr = -(d_hist - med) / sd p_nat = np.mean(d_nat <= d_hist) return dict(n_nat=len(d_nat), d_hist=d_hist, nat_median=med, nat_p05=np.percentile(d_nat, 5), nat_p95=np.percentile(d_nat, 95), noise=sd, signal=d_hist - med, snr=snr, p_nat=p_nat, p_gauss=norm.cdf(snr)) def fmt_p(p, n): return f'< 1/{n}' if p == 0 else f'{p:.1%}' # %% [markdown] # ## Part 1: regional volume and area # %% reg = {} for var, unit in [('volume', 'km$^3$'), ('area', 'km$^2$')]: short = 'km3' if var == 'volume' else 'km2' nat = pd.read_csv(f'{FILTERED}/{var}_{short}_nat_ens_reg11.csv', index_col=0) hist = pd.read_csv(f'{FILTERED}/{var}_{short}_hist_ens_reg11.csv', index_col=0)[MEMBER] nat.index = nat.index.astype(int) hist.index = hist.index.astype(int) d_nat, d_hist = delta(nat), delta(hist.to_frame())[MEMBER] years = np.arange(BASELINE[1] + 1, T_END + 1) ts = pd.DataFrame([stats(d_hist[y], d_nat.loc[y].values) for y in years], index=pd.Index(years, name='year')) ts['d_nat_same_member'] = d_nat[MEMBER].reindex(years).values reg[var] = dict(nat=nat, hist=hist, d_nat=d_nat, d_hist=d_hist, ts=ts, unit=unit) snap = ts.loc[TIMES].copy() snap.to_csv(f'{FIGS}/regional_{var}_{TAG}.csv') print(f'\n{var} ({MEMBER}, {WINDOW}-yr mean ending at t, minus {BASELINE[0]}-{BASELINE[1]}):') for y, r in snap.iterrows(): print(f' {y}: dHIST {r.d_hist:8.1f} NAT median {r.nat_median:8.1f} ' f'(P5-P95 {r.nat_p05:7.1f} to {r.nat_p95:7.1f}, N={r.n_nat:.0f}) ' f'SNR {r.snr:5.2f} p_nat {fmt_p(r.p_nat, int(r.n_nat))}') # %% # Figure 1: the change since the baseline over time, NAT spread vs HIST, # with the SNR below for var, R in reg.items(): ts, d_nat, d_hist = R['ts'], R['d_nat'], R['d_hist'] yrs = ts.index fig, (ax, ax2) = plt.subplots(2, 1, figsize=(10, 7.5), sharex=True, gridspec_kw=dict(height_ratios=[3, 2], hspace=0.08)) q = d_nat.loc[yrs].quantile([0.05, 0.25, 0.5, 0.75, 0.95], axis=1).T ax.fill_between(yrs, q[0.05], q[0.95], color=NAT_C, alpha=0.18, lw=0, label='NAT ensemble, 5–95%') ax.fill_between(yrs, q[0.25], q[0.75], color=NAT_C, alpha=0.30, lw=0, label='NAT ensemble, 25–75%') ax.plot(yrs, q[0.5], color=NAT_DARK, lw=2, label='NAT median') ax.plot(yrs, d_nat[MEMBER].loc[yrs], color=NAT_DARK, lw=1, ls='--', label=f'NAT {MEMBER} (same pre-1850 forcing)') ax.plot(yrs, d_hist.loc[yrs], color=HIST_C, lw=2.5, label=f'HIST {MEMBER}') ax.axhline(0, color='0.4', lw=0.8) ax.set_ylabel(f'Δ{var} vs {BASELINE[0]}–{BASELINE[1]} ({R["unit"]})') ax.set_title(f'Regional {var}, Alps ({EXP}): {WINDOW}-yr mean ending at t, ' f'minus {BASELINE[0]}–{BASELINE[1]}') ax.legend(loc='lower left', fontsize=9, frameon=False) ax.grid(alpha=0.3) ax2.plot(yrs, ts.snr, color='k', lw=2) for (th, lab), c in zip(SNR_THRESH, ['#f5b041', '#e67e22', '#c0392b']): ax2.axhline(th, color=c, ls='--', lw=1.5, label=lab) for y in TIMES: ax2.plot(y, ts.snr[y], 'o', color='k', ms=6) ax2.annotate(f'{ts.snr[y]:.1f}', (y, ts.snr[y]), textcoords='offset points', xytext=(0, 8), ha='center', fontsize=9) # Where the NAT ensemble shrinks for y in [2007, 2011, 2014, 2016]: ax2.axvline(y - 0.5, color='0.6', lw=0.6, ls=':') ax2.text(2004, ax2.get_ylim()[1] * 0.98, f'N NAT: 197 → {ts.n_nat.iloc[-1]}', ha='right', va='top', fontsize=8, color='0.4') ax2.set_ylim(top=ax2.get_ylim()[1] + 0.5) ax2.set_ylabel('SNR (HIST − NAT median) / NAT std') ax2.set_xlabel('Year t') ax2.legend(loc='upper left', fontsize=9, frameon=False) ax2.grid(alpha=0.3) ax2.set_xlim(yrs[0], yrs[-1]) fig.savefig(f'{FIGS}/regional_{var}_timeseries_{TAG}.png', dpi=150, bbox_inches='tight') plt.close(fig) # %% # Figure 2: the NAT distribution of the change at the snapshot times, with # the HIST value on top for var, R in reg.items(): ts, d_nat = R['ts'], R['d_nat'] fig, axes = plt.subplots(1, len(TIMES), figsize=(2.6 * len(TIMES), 3.6), sharex=True) allv = np.concatenate([d_nat.loc[TIMES].values.ravel(), ts.loc[TIMES, 'd_hist'].values]) lo, hi = np.nanmin(allv), np.nanmax(allv) pad = 0.05 * (hi - lo) bins = np.linspace(lo - pad, hi + pad, 50) for ax, y in zip(axes, TIMES): r = ts.loc[y] v = d_nat.loc[y].dropna().values ax.hist(v, bins=bins, color=NAT_C, alpha=0.7, edgecolor='white', linewidth=0.5) ax.axvline(r.nat_median, color=NAT_DARK, lw=1.5) ax.axvline(r.d_hist, color=HIST_C, lw=2.5) ax.set_title(f'{y}', fontweight='bold') ax.text(0.03, 0.97, f'SNR {r.snr:.1f}\np_nat {fmt_p(r.p_nat, int(r.n_nat))}' f'\nN={int(r.n_nat)}', transform=ax.transAxes, va='top', fontsize=9) ax.set_yticks([]) ax.set_xlabel(f'Δ{var} ({R["unit"]})', fontsize=9) ax.tick_params(labelsize=8) axes[0].set_ylabel('NAT members') fig.suptitle(f'Regional Δ{var} vs {BASELINE[0]}–{BASELINE[1]}: NAT ensemble ' f'(blue) vs HIST {MEMBER} (red)', y=1.02) fig.savefig(f'{FIGS}/regional_{var}_snapshots_{TAG}.png', dpi=150, bbox_inches='tight') plt.close(fig) # %% [markdown] # ## Part 2: per glacier # # Same numbers for every glacier of the common subset. The cube # (member, year, glacier) is cached in figs/, delete it to rebuild. # %% if os.environ.get('REGIONAL_ONLY') == '1': raise SystemExit(0) VARS = {'volume': 'volume', 'area': 'area_min_h', 'length': 'length'} def load_member(m, ids): with xr.open_dataset(f'{OUTPUT}/lmr_mira_nat_{m}_reg11.nc') as ds: return (ds[list(VARS.values())].sel(time=slice(1800, T_END), rgi_id=ids) .reindex(time=np.arange(1800, T_END + 1)).load()) ids = pd.read_csv(f'{FILTERED}/glacier_statistics_reg11.csv', index_col=0) rgi_area = ids['rgi_area_km2'] ids = ids.index cache = f'{FIGS}/cube_{EXP}.nc' if not os.path.exists(cache): members = sorted(f[len('lmr_mira_nat_'):-len('_reg11.nc')] for f in os.listdir(OUTPUT) if f.startswith('lmr_mira_nat_ens')) # Slow (~1 min per member, CPU-bound in netCDF decoding): in parallel from concurrent.futures import ProcessPoolExecutor with ProcessPoolExecutor(int(os.environ.get('NJOBS', 8))) as ex: arrs = list(ex.map(load_member, members, [ids] * len(members))) nat_cube = xr.concat(arrs, dim=pd.Index(members, name='member')) nat_cube.to_netcdf(cache) nat_cube = xr.open_dataset(cache).load() nat_cube['time'] = nat_cube.time.astype(int) with xr.open_dataset(f'{OUTPUT}/lmr_mira_hist_{MEMBER}_reg11.nc') as ds: hist_cube = ds[list(VARS.values())].sel(time=slice(1800, T_END), rgi_id=ids).load() hist_cube['time'] = hist_cube.time.astype(int) # %% glac = [] snr_at = {} for name, v in VARS.items(): dn = delta(nat_cube[v]) dh = delta(hist_cube[v]) for y in TIMES: dny = dn.sel(time=y) n = dny.notnull().sum('member') med = dny.median('member') sd = dny.std('member', ddof=1) dhy = dh.sel(time=y) snr = -(dhy - med) / sd p_nat = (dny <= dhy).sum('member') / n ok = np.isfinite(snr) & (sd > 0) snr_at[(name, y)] = snr.where(ok) s = snr.where(ok).to_series().dropna() w = rgi_area.reindex(s.index) row = dict(var=name, year=y, n_glaciers=len(s), median_snr=s.median()) for th, lab in SNR_THRESH: row[f'pct_glaciers_snr>{th}'] = 100 * (s > th).mean() row[f'pct_area_snr>{th}'] = 100 * w[s > th].sum() / w.sum() row['pct_glaciers_p_nat<5%'] = 100 * (p_nat.where(ok).to_series().dropna() < 0.05).mean() glac.append(row) glac = pd.DataFrame(glac).set_index(['var', 'year']) glac.to_csv(f'{FIGS}/glaciers_{TAG}.csv') with pd.option_context('display.width', 200, 'display.precision', 1): print('\nPer glacier (% of glaciers / % of RGI area above SNR thresholds):') print(glac) # %% # Figure 3: Mira's histogram, for the last time (modern period) y = T_END fig, axes = plt.subplots(1, 3, figsize=(18, 5.5)) for ax, name, c in zip(axes, VARS, ['steelblue', 'coral', 'seagreen']): s = snr_at[(name, y)].to_series().dropna() r = glac.loc[(name, y)] ax.hist(s.clip(-2, 10), bins=np.linspace(-2, 10, 97), color=c, alpha=0.75, edgecolor='white', linewidth=0.4) for (th, lab), lc in zip(SNR_THRESH, ['#f5b041', '#e67e22', '#c0392b']): ax.axvline(th, color=lc, ls='--', lw=2, label=lab) ax.axvline(s.median(), color='purple', lw=2, label=f'Median = {s.median():.2f}') ax.set_xlabel('Signal-to-noise ratio (clipped to [−2, 10])', fontsize=12) ax.set_ylabel('Number of glaciers', fontsize=12) ax.set_title(f'{name.upper()}', fontsize=14, fontweight='bold') ax.legend(fontsize=9, loc='upper right', frameon=False) ax.grid(alpha=0.3) txt = '\n'.join([f'N = {int(r.n_glaciers)}'] + [f'{lab.split(" (")[0]}: {r[f"pct_glaciers_snr>{th}"]:.0f}% of glaciers, ' f'{r[f"pct_area_snr>{th}"]:.0f}% of area' for th, lab in SNR_THRESH]) ax.text(0.97, 0.70, txt, transform=ax.transAxes, va='top', ha='right', fontsize=9.5, bbox=dict(boxstyle='round', facecolor='white', edgecolor='0.7', alpha=0.85)) fig.suptitle(f'Per-glacier SNR, HIST {MEMBER} vs NAT, {y - WINDOW + 1}–{y} mean minus ' f'{BASELINE[0]}–{BASELINE[1]} ({EXP})', fontsize=14) fig.savefig(f'{FIGS}/glaciers_snr_hist_{y}_{TAG}.png', dpi=150, bbox_inches='tight') plt.close(fig) # %% # Figure 4: share of glaciers above the thresholds, over the snapshot times fig, axes = plt.subplots(1, 3, figsize=(15, 4.2), sharey=True) for ax, name in zip(axes, VARS): g = glac.loc[name] for (th, lab), lc in zip(SNR_THRESH, ['#f5b041', '#e67e22', '#c0392b']): ax.plot(g.index, g[f'pct_glaciers_snr>{th}'], 'o-', color=lc, lw=2, label=lab) ax.plot(g.index, g[f'pct_area_snr>{th}'], 's:', color=lc, lw=1.5) ax.set_title(name.upper(), fontweight='bold') ax.set_xlabel('Year t') ax.grid(alpha=0.3) ax.set_ylim(0, 100) axes[0].set_ylabel('% of glaciers (solid) / of area (dotted)') axes[0].legend(fontsize=9, frameon=False, loc='upper left') fig.suptitle(f'Per-glacier attribution over time, HIST {MEMBER} ({EXP})', y=1.02) fig.savefig(f'{FIGS}/glaciers_snr_time_{TAG}.png', dpi=150, bbox_inches='tight') plt.close(fig) print('\nFigures in', FIGS)