"""Build the observed (forecast-independent) targets + climatology from VnGAD Tx.
Outputs (all derived; originals untouched):
data/processed/obs_daily_ehd_2003_2022.nc (ehd int8, exceedance f32)
data/processed/obs_climatology_2003_2022.nc (per-cell maps)
outputs/metrics/obs_episode_catalog_2003_2022.csv
outputs/metrics/obs_climatology_by_region.csv
outputs/figures/obs_{ehd_frequency,mean_episode_duration,mean_severity}.png
outputs/OBS_CLIMATOLOGY.html
"""
import os
import sys
import numpy as np
import xarray as xr
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
sys.path.insert(0, os.path.join(os.path.dirname(__file__), "..", "src"))
from targets.obs import build_daily, episode_catalog, cell_climatology # noqa: E402
ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
DATA = os.path.join(ROOT, "data")
PROC = os.path.join(DATA, "processed")
WIN = ("2003-01-01", "2022-12-31")
BUF = ("2002-12-20", "2023-01-12") # buffer for boundary episodes
N_YEARS = 20
def log(*a):
print(*a, flush=True)
log("[1/6] load Tx, threshold, region field")
tx = xr.open_dataset(os.path.join(DATA, "VnGAD_daily-0.25_1961-2023_tasmax.nc"))
thr = xr.open_dataset(os.path.join(PROC, "EHD_THRESHOLD_2003_2022.nc"))["EHD_THRESHOLD"].values
reg = xr.open_dataset(os.path.join(PROC, "region_codes_vngad.nc"))["region"].values
lat = tx["lat"].values; lon = tx["lon"].values
valid = reg >= 0 # national domain (0..7); -1 outside
log("[2/6] daily EHD/exceedance on buffered window")
sub = tx["tasmax"].sel(time=slice(BUF[0], BUF[1]))
dates_buf = sub["time"].values
ehd_buf, exc_buf = build_daily(sub.values, dates_buf, thr)
tx.close() # release the source HDF5 handle before writing new NetCDF files
log("[3/6] episode catalog (continuous, onsets counted in-window)")
catalog = episode_catalog(ehd_buf, exc_buf, dates_buf, valid, reg, WIN)
log(" episodes total:", len(catalog),
"| onsets in-window:", sum(r["onset_in_window"] for r in catalog))
log("[4/6] per-cell climatology")
# window slice for daily frequency
win_mask = (dates_buf.astype("datetime64[D]") >= np.datetime64(WIN[0])) & \
(dates_buf.astype("datetime64[D]") <= np.datetime64(WIN[1]))
ehd_win = ehd_buf[win_mask]
clim = cell_climatology(ehd_win, valid, catalog, valid.shape, N_YEARS)
log(" national mean daily EHD freq:", round(float(np.nanmean(clim['ehd_freq'])), 4))
log("[5/6] write NetCDF + CSV")
# daily fields (compressed)
dts_win = dates_buf[win_mask].astype("datetime64[ns]")
xr.Dataset(
{"ehd": (("time", "lat", "lon"), ehd_win.astype("int8")),
"exceedance": (("time", "lat", "lon"), exc_buf[win_mask].astype("float32"))},
coords={"time": dts_win, "lat": lat, "lon": lon},
attrs={"title": "Observed daily EHD and exceedance (VnGAD, 2003-2022)",
"note": "EHD=1 if Tx>=min(P90,35); exceedance=max(Tx-thr,0)"},
).to_netcdf(os.path.join(PROC, "obs_daily_ehd_2003_2022.nc"))
cmap_vars = {k: (("lat", "lon"), v.astype("float32")) for k, v in clim.items()}
xr.Dataset(cmap_vars, coords={"lat": lat, "lon": lon},
attrs={"title": "Observed extreme-heat climatology (VnGAD, 2003-2022)",
"reference_period": "2003-2022"}
).to_netcdf(os.path.join(PROC, "obs_climatology_2003_2022.nc"))
# episode catalog CSV
cols = ["lat_i", "lon_j", "region", "onset", "end", "duration", "i_mean",
"i_max", "severity", "left_censored", "right_censored",
"onset_in_window", "complete"]
with open(os.path.join(ROOT, "outputs/metrics/obs_episode_catalog_2003_2022.csv"), "w") as f:
f.write("lat,lon," + ",".join(cols[2:]) + "\n")
for r in catalog:
f.write(f"{lat[r['lat_i']]:.3f},{lon[r['lon_j']]:.3f}," +
",".join(str(r[c]) for c in cols[2:]) + "\n")
# per-region summary
reg_rows = ["region,n_cells,mean_ehd_freq,episodes_per_year_mean,mean_duration,mean_severity"]
for code in [0, 1, 2, 3, 4, 5, 6, 7]:
m = (reg == code) & valid
if not m.any():
continue
reg_rows.append(
f"{code},{int(m.sum())},{np.nanmean(clim['ehd_freq'][m]):.4f},"
f"{np.nanmean(clim['episodes_per_year'][m]):.3f},"
f"{np.nanmean(clim['mean_duration'][m]):.3f},"
f"{np.nanmean(clim['mean_severity'][m]):.3f}")
allm = valid
reg_rows.append(
f"ALL,{int(allm.sum())},{np.nanmean(clim['ehd_freq'][allm]):.4f},"
f"{np.nanmean(clim['episodes_per_year'][allm]):.3f},"
f"{np.nanmean(clim['mean_duration'][allm]):.3f},"
f"{np.nanmean(clim['mean_severity'][allm]):.3f}")
with open(os.path.join(ROOT, "outputs/metrics/obs_climatology_by_region.csv"), "w") as f:
f.write("\n".join(reg_rows) + "\n")
log("[6/6] figures + HTML")
def _map(field, title, fname, cmap="YlOrRd"):
fig, ax = plt.subplots(figsize=(5.2, 7))
im = ax.pcolormesh(lon, lat, field, cmap=cmap, shading="nearest")
ax.set_title(title); ax.set_xlabel("lon"); ax.set_ylabel("lat"); ax.set_aspect("equal")
fig.colorbar(im, ax=ax, shrink=0.8)
fig.tight_layout()
fig.savefig(os.path.join(ROOT, "outputs/figures", fname), dpi=120)
plt.close(fig)
_map(clim["ehd_freq"], "Observed daily EHD frequency (2003-2022)", "obs_ehd_frequency.png")
_map(clim["mean_duration"], "Observed mean HW episode duration (days)", "obs_mean_episode_duration.png")
_map(clim["mean_severity"], "Observed mean HW episode severity (degC day)", "obs_mean_severity.png")
reg_table = "\n".join(
"" + "".join(f"| {c} | " for c in row.split(",")) + "
"
for row in reg_rows)
html = f"""
S2S-FM-VN Observed Extreme-Heat Climatology (2003-2022)
Observed Extreme-Heat Climatology — VnGAD, 2003–2022
Forecast-independent observed truth from the locked P90/EHD/episode definitions.
National domain = all valid VnGAD cells (region 0–7); regional rows exclude the
unassigned region 0. Episodes detected on the continuous daily series with a boundary buffer.
By region
| region | n_cells | mean EHD freq |
episodes/yr | mean duration (d) | mean severity (degC·d) |
{reg_table}
Maps
Sources: data/processed/obs_climatology_2003_2022.nc,
outputs/metrics/obs_episode_catalog_2003_2022.csv, obs_climatology_by_region.csv.
"""
with open(os.path.join(ROOT, "outputs/OBS_CLIMATOLOGY.html"), "w") as f:
f.write(html)
log("DONE")