"""Convert spinup-geometry regional netcdf files to the GLAMBIE submission CSV. Unlike the fixed-geometry contribution, glacier_area_reference_start/end vary row to row here (dynamic geometry), taken directly from each region's time-varying `covered_area_m2`. Checks run before writing: 1. All 19 RGI regions are present. 2. No NaN in glacier_change_observed (hard requirement per the submission guide). 3. Area check: total_region_area_m2 (full RGI catalogue, incl. glaciers that failed before ever reaching the hydro output) matches an independent glacier_statistics_RR.csv sum. 4. Self-consistency: the regional series is recomputed directly from the per-glacier file (same mass-sum formula) and compared to the stored region file. Not a scientific validation (same underlying computation), but catches file mismatches/regressions if the two are ever regenerated separately. 5. Optional, informational only: if a fixed-geometry v1.7a submission is already present, prints a rough regional-magnitude comparison. """ import glob import os import numpy as np import pandas as pd import xarray as xr SUMMARY_DIR = os.path.expanduser( "~/www_oggm/gdirs/oggm_v1.7a/tests/L3-L5_files/elev_bands/consensus/ERA5_Glambie/" "RGI62/b_160/L4/summary" ) NETCDF_DIR = os.path.join(os.path.dirname(__file__), "glambie_spinup_geometry") OUTPUT_CSV = os.path.join(os.path.dirname(__file__), "glambie_submission_spinup_geometry_v17a_TIModel.csv") FIXED_GEOM_CSV = os.path.join( os.path.dirname(__file__), "glambie_submission_fixed_geometry_v17a_TIModel.csv" ) COLUMNS = [ "region_id", "start_date", "end_date", "glacier_change_observed", "glacier_change_uncertainty", "glacier_change_observed_above_water_level", "glacier_change_observed_above_water_level_uncertainty", "glacier_change_observed_below_water_level", "glacier_change_observed_below_water_level_uncertainty", "glacier_frontal_ablation", "glacier_frontal_ablation_uncertainty", "unit", "glacier_area_reference_start", "glacier_area_reference_end", "observational_coverage_percentage", "remarks", ] def nc_to_rows(nc_path): with xr.open_dataset(nc_path) as ds: times = pd.DatetimeIndex(ds["time"].values) mb = ds["region_specific_mb_mwe"].values area_km2 = ds["covered_area_m2"].values / 1e6 # time-varying (dynamic geometry) coverage = float(ds["completion_rate_pct"].values) region_id = int(ds.attrs["region_id"]) rows = [] for i, t in enumerate(times): end = t + pd.offsets.MonthBegin(1) rows.append({ "region_id": region_id, "_sort_time": t, "start_date": t.strftime("%d/%m/%Y"), "end_date": end.strftime("%d/%m/%Y"), "glacier_change_observed": mb[i], "glacier_change_uncertainty": np.nan, "glacier_change_observed_above_water_level": np.nan, "glacier_change_observed_above_water_level_uncertainty": np.nan, "glacier_change_observed_below_water_level": np.nan, "glacier_change_observed_below_water_level_uncertainty": np.nan, "glacier_frontal_ablation": np.nan, "glacier_frontal_ablation_uncertainty": np.nan, "unit": "mwe", "glacier_area_reference_start": area_km2[i], "glacier_area_reference_end": area_km2[i], "observational_coverage_percentage": coverage, "remarks": np.nan, }) return rows # --------------------------------------------------------------------------- nc_files = sorted(glob.glob(os.path.join(NETCDF_DIR, "region_specific_mb_spinup_RGI*_v17a_TIModel.nc"))) if not nc_files: raise RuntimeError(f"No regional netcdf files found in {NETCDF_DIR}") # ---- Check 1: all 19 regions ----------------------------------------------- region_ids = [] for fp in nc_files: with xr.open_dataset(fp) as ds: region_ids.append(int(ds.attrs["region_id"])) missing = sorted(set(range(1, 20)) - set(region_ids)) if missing: raise RuntimeError(f"Missing RGI regions: {[f'{r:02d}' for r in missing]}") print("[OK] All 19 regions present\n") # ---- Build rows (needed for check 2, then reused for the final CSV) -------- all_rows = [] for fp in nc_files: all_rows.extend(nc_to_rows(fp)) df = pd.DataFrame(all_rows) df = df.sort_values(["region_id", "_sort_time"]).reset_index(drop=True) # ---- Check 2: no NaN in glacier_change_observed ----------------------------- nan_rows = df[df["glacier_change_observed"].isna()] if not nan_rows.empty: bad_regions = sorted(nan_rows["region_id"].unique()) raise RuntimeError( f"glacier_change_observed has {len(nan_rows)} NaN rows in regions {bad_regions} " "— every region-month must resolve to a finite value" ) print("[OK] No NaN in glacier_change_observed\n") # ---- Check 3: area vs glacier_statistics (full catalogue) ------------------ print("Area check (full RGI catalogue, incl. glaciers that never reached the hydro output):") area_warnings = [] for fp in nc_files: with xr.open_dataset(fp) as ds: region_id = f"{int(ds.attrs['region_id']):02d}" nc_area = float(ds["total_region_area_m2"].values) / 1e6 coverage = float(ds["completion_rate_pct"].values) stats_path = os.path.join(SUMMARY_DIR, f"glacier_statistics_{region_id}.csv") expected_area = pd.read_csv(stats_path, index_col=0, low_memory=False)["rgi_area_km2"].sum() diff = nc_area - expected_area flag = " *** MISMATCH" if abs(diff) > 1.0 else "" if flag: area_warnings.append(f"RGI{region_id}") print(f" RGI{region_id}: expected={expected_area:.3f} got={nc_area:.3f} diff={diff:+.3f} km²" f" coverage={coverage:.1f}%{flag}") if not area_warnings: print("[OK] All regional catalogue areas match glacier_statistics\n") else: print(f"\n[WARN] Area mismatch in: {', '.join(area_warnings)}\n") # ---- Check 4: self-consistency vs per-glacier file -------------------------- print("Self-consistency check (region file vs. recomputed from per-glacier file):") consistency_warnings = [] for fp in nc_files: with xr.open_dataset(fp) as ds: region_id = f"{int(ds.attrs['region_id']):02d}" region_mb = ds["region_specific_mb_mwe"].values region_area = ds["covered_area_m2"].values glac_fp = os.path.join(NETCDF_DIR, f"per_glacier_specific_mb_spinup_RGI{region_id}_v17a_TIModel.nc") with xr.open_dataset(glac_fp) as gds: succ = gds["task_success"].values.astype(bool) mb_mwe = gds["specific_mb_mwe"].values[succ] # (n_succ, n_months) area_m2 = gds["area_min_h_m2"].values[succ] recomputed_mb_kg = (mb_mwe * area_m2 * 1000.0).sum(axis=0) recomputed_area = area_m2.sum(axis=0) recomputed_mb = recomputed_mb_kg / recomputed_area / 1000.0 max_diff = np.nanmax(np.abs(recomputed_mb - region_mb)) area_diff = np.nanmax(np.abs(recomputed_area - region_area)) # Netcdf variables are stored as float32; tolerance is set well above that # noise floor (~1e-6 typical) but far below any real-bug-scale difference. flag = " *** MISMATCH" if max_diff > 1e-4 or area_diff > 1.0 else "" if flag: consistency_warnings.append(f"RGI{region_id}") print(f" RGI{region_id}: max mb diff={max_diff:.2e} m w.e. max area diff={area_diff:.2e} m²{flag}") if consistency_warnings: raise RuntimeError(f"Region/per-glacier mismatch in: {', '.join(consistency_warnings)}") print("[OK] Region files are consistent with per-glacier files\n") # ---- Check 5 (optional, informational): compare vs fixed-geometry ---------- if os.path.exists(FIXED_GEOM_CSV): print("Informational comparison vs fixed-geometry v17a_TIModel (cumulative 1975-2025, m w.e.):") fixed = pd.read_csv(FIXED_GEOM_CSV, parse_dates=["start_date"], dayfirst=True) for region_id in sorted(df["region_id"].unique()): spin_cum = df.loc[df["region_id"] == region_id, "glacier_change_observed"].sum() fix_cum = fixed.loc[fixed["region_id"] == region_id, "glacier_change_observed"].sum() print(f" RGI{region_id:02d}: spinup={spin_cum:+.2f} fixed={fix_cum:+.2f} diff={spin_cum - fix_cum:+.2f}") print() else: print(f"[SKIP] Fixed-geometry submission not found yet at {FIXED_GEOM_CSV} — skipping comparison\n") # ---- Write CSV --------------------------------------------------------------- df = df[COLUMNS] df["glacier_change_observed"] = df["glacier_change_observed"].round(4) df["glacier_change_observed_above_water_level"] = df["glacier_change_observed"] df["glacier_area_reference_start"] = df["glacier_area_reference_start"].round(3) df["glacier_area_reference_end"] = df["glacier_area_reference_end"].round(3) df.to_csv(OUTPUT_CSV, index=False, na_rep="") print(f"Wrote {len(df)} rows ({df['region_id'].nunique()} regions) → {OUTPUT_CSV}")