#!/usr/bin/env python """Stage 26 — item 8: real-time case study, initialization 2026-08-29 (../HW_EC_10). Runs the deployable chain on a genuine real-time forecast and produces exactly three products: 1. P(yEHDE | g, Wk) maps over Vietnam, one panel per week W1..W6 (region boundaries overlaid). 2. 7x6 table of the regional event yEHDE(Rj,Wk) (yes/no). 3. 7x6 table of P(yEHDE | Rj,Wk). Each is delivered as a machine-readable file and a report-ready figure/CSV. REAL-TIME vs REFORECAST handling (item 8 / known_failure_modes), checked, not assumed: * The ../HW real-time file `txtm_cf_2026-08-29.nc` is CONTROL-ONLY (1 member, `_cf_`), matching the single-member reforecast the model was calibrated on — so the single-member MOS applies. Any ensemble would be reduced to the control first; there is none here. * It is on a COARSER grid (10x13, ~1.5 deg) than the QG_1 0.25 deg reforecast grid, so daily Tx is bilinearly regridded to the 0.25 deg grid. This upsampling is a documented limitation (the reforecast was pre-interpolated upstream; the real-time file is not). * The v3 EXPANSION members are NOT applied here: their persistence predictors need observations dated before the init, which the archive (ends 2023) does not contain for 2026. The real-time product therefore uses the deployable PRODUCTION model (all-years HGB + isotonic, forecast-only features; outputs/production_model). The v3 ensemble/expansion is evaluated on the reforecast CV. * VERIFICATION: the gridded obs archive ends 2023-12-31, so NO observation exists for this 2026 case. All six weeks are labelled NOT-YET-VERIFIABLE; nothing is presented as validated. Outputs (outputs/realtime_2026-08-29/ + results/): ehd_realtime_2026-08-29.nc prob_ehd + ehd_flag (block,lat,lon), CF-1.8 v8_regional_event_2026-08-29.csv 7x6 yes/no v8_regional_prob_2026-08-29.csv 7x6 P results/figures/v8_realtime_maps_2026-08-29.png """ from __future__ import annotations import os, sys, json, logging import numpy as np import pandas as pd import xarray as xr import joblib sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..", "src")) from ehd.config import load_config, lead_blocks from ehd import features as feat, area_weights as aw logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s") log = logging.getLogger("stage26") RT_FILE = "/home/tanpv/Progs/claude/HW_EC_10/data/nc_rt/txtm/2026/txtm_cf_2026-08-29.nc" INIT = "2026-08-29" STEPS_PER_DAY = 4 NDAYS = 46 def build_daily_tx_on_grid(rt_path, ref): """Coarse 6-hourly mx2t6 -> daily Tx (max of 4) -> bilinear regrid to the QG_1 grid. Returns a DataArray (init=1, lead_day=0..45, lat, lon) with a valid_date coord.""" ds = xr.open_dataset(rt_path) mx = ds["mx2t6"] # (time=184, lat, lon) nt = mx.sizes["time"]; nday = min(NDAYS, nt // STEPS_PER_DAY) daily = mx.isel(time=slice(0, nday * STEPS_PER_DAY)).coarsen( time=STEPS_PER_DAY, boundary="trim").max() # daily max of the four 6-h values daily = daily.rename({"latitude": "lat", "longitude": "lon"}).sortby("lat").sortby("lon") daily = daily.interp(lat=ref.lat, lon=ref.lon, method="linear") # regrid to 0.25 deg vdates = pd.to_datetime(INIT) + pd.to_timedelta(np.arange(nday), unit="D") out = xr.DataArray(daily.values[None].astype("float32"), dims=("init", "lead_day", "lat", "lon"), coords={"init": [pd.to_datetime(INIT)], "lead_day": np.arange(nday), "lat": ref.lat, "lon": ref.lon}) out = out.assign_coords(valid_date=(("init", "lead_day"), vdates.values[None])) return out - 273.15 if float(np.nanmax(out)) > 100 else out # K -> degC if needed def main(): cfg = load_config() blocks = lead_blocks(cfg) block_names = list(blocks.keys()) region_ids = cfg["regions"]["ids"] mdl_dir = os.path.join(cfg["paths"]["outputs"], "production_model") art = joblib.load(os.path.join(mdl_dir, "ehd_production_model.joblib")) thr = xr.open_dataarray(os.path.join(mdl_dir, "tx90p_threshold.nc")) iy, ix = art["land_iy"], art["land_ix"] lat, lon = art["lat"], art["lon"] region_of_cell = art["region_of_cell"] ref = xr.open_dataarray(os.path.join(cfg["paths"]["outputs"], "forecast_daily_tx", "fc_daily_tx_2003.nc")).isel(init=0) fc = build_daily_tx_on_grid(RT_FILE, ref) log.info("real-time daily Tx built & regridded: %s, Tx range %.1f..%.1f degC", dict(fc.sizes), float(np.nanmin(fc)), float(np.nanmax(fc))) thr_i = thr.reindex(lat=fc.lat, lon=fc.lon, method="nearest") feats = feat.forecast_features(fc, thr_i, {k: tuple(v) for k, v in blocks.items()}, iy, ix) # area weights for regional aggregation wcell, rule = aw.cell_area_weights(lat, lon, region_of_cell, iy, ix, os.path.join(cfg["paths"]["shapefile_dir"], "Vietnam7.shp")) log.info("area-weight rule: %s", rule) ncell = iy.size prob = np.full((len(block_names), lat.size, lon.size), np.nan, np.float32) flag = np.full((len(block_names), lat.size, lon.size), -1, np.int8) reg_prob = {}; reg_event = {} for k, b in enumerate(block_names): X = feats[b]["X"] p = art["per_block"][b]["iso"].transform(art["per_block"][b]["hgb"].predict_proba(X)) cthr = np.array([art["per_block"][b]["region_thr"].get(int(r), art["per_block"][b]["global_thr"]) for r in region_of_cell]) fl = (p >= cthr).astype(np.int8) prob[k, iy, ix] = p; flag[k, iy, ix] = fl # regional: area-weighted mean probability; event = area-weighted flag coverage >= 0.5 for rj in region_ids: cm = region_of_cell == rj; wj = wcell[cm]; sw = wj.sum() reg_prob[(rj, b)] = float((p[cm] * wj).sum() / sw) reg_event[(rj, b)] = int((fl[cm] * wj).sum() / sw >= 0.5) outdir = os.path.join(cfg["paths"]["outputs"], f"realtime_{INIT}") os.makedirs(outdir, exist_ok=True) ds = xr.Dataset( {"prob_ehd": (("block", "lat", "lon"), prob, {"long_name": "P(>=1 EHD in block), real-time production MOS", "units": "1"}), "ehd_flag": (("block", "lat", "lon"), flag, {"long_name": "operational EHD flag (prob >= region threshold); -1 = sea", "units": "1"})}, coords={"block": np.arange(1, len(block_names) + 1), "lat": lat.astype("float32"), "lon": lon.astype("float32"), "init": np.datetime64(INIT)}, attrs={"Conventions": "CF-1.8", "title": f"Real-time EHD forecast, init {INIT} (../HW_EC_10 control, regridded)", "verification_status": "NOT YET VERIFIABLE — obs archive ends 2023; no 2026 obs.", "real_time_notes": ("control-only single member; coarse forecast bilinearly regridded " "to 0.25deg; v3 expansion/persistence not applicable real-time; " "product = deployable production model."), "block_names": " ".join(block_names)}) ds.to_netcdf(os.path.join(outdir, f"ehd_realtime_{INIT}.nc")) # 7x6 tables ev = pd.DataFrame([[reg_event[(r, b)] for b in block_names] for r in region_ids], index=[f"R{r}" for r in region_ids], columns=block_names) pr = pd.DataFrame([[round(reg_prob[(r, b)], 3) for b in block_names] for r in region_ids], index=[f"R{r}" for r in region_ids], columns=block_names) tdir = os.path.join(cfg["paths"]["results"], "tables") ev.to_csv(os.path.join(tdir, f"v8_regional_event_{INIT}.csv")) pr.to_csv(os.path.join(tdir, f"v8_regional_prob_{INIT}.csv")) # maps import matplotlib; matplotlib.use("Agg"); import matplotlib.pyplot as plt reg = xr.open_dataset(cfg["paths"]["mask"])["reg"].reindex(lat=fc.lat, lon=fc.lon, method="nearest") fig, axes = plt.subplots(2, 3, figsize=(13, 9)) for k, (ax, b) in enumerate(zip(axes.ravel(), block_names)): im = ax.pcolormesh(lon, lat, prob[k], cmap="YlOrRd", vmin=0, vmax=1, shading="nearest") ax.contour(lon, lat, np.nan_to_num(reg.values), levels=np.arange(1.5, 7.5), colors="k", linewidths=0.4) ax.set_title(f"{b} (days {blocks[b][0]}-{blocks[b][1]})", fontsize=10) ax.set_xlabel("lon"); ax.set_ylabel("lat") fig.colorbar(im, ax=axes, shrink=0.6, label="P(>=1 EHD in block)") fig.suptitle(f"Real-time EHD probability, init {INIT} — NOT YET VERIFIABLE (no 2026 obs)") fig.savefig(os.path.join(cfg["paths"]["results"], "figures", f"v8_realtime_maps_{INIT}.png"), dpi=130, bbox_inches="tight") print("=== regional event yEHDE(Rj,Wk) ==="); print(ev.to_string()) print("\n=== regional P(yEHDE|Rj,Wk) ==="); print(pr.to_string()) print(f"\nverification: NOT YET VERIFIABLE (obs archive ends 2023). products -> {outdir}") if __name__ == "__main__": main()