#!/usr/bin/env python3 """Render the verification metrics as one self-contained HTML report. Layout follows section 4.2.4 of HW_HPC_1: one row per metric, the six forecast weeks W1..W6 side by side and tight, all sharing a single colour scale so the lead-time degradation is readable straight off the row. Maps are drawn on the 541 land points of Mask_VN_S2S_EC_0.25.nc with the seven climatic sub-region boundaries of Vietnam7.shp on top. Figures are embedded as base64 PNG so the HTML file stands alone. """ from __future__ import annotations import argparse import base64 import io import os import sys from html import escape import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np import shapefile from matplotlib.colors import BoundaryNorm from netCDF4 import Dataset sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import ehd_common as C # metric -> (title, colormap, symmetric-about-zero, fixed range or None) SPEC = { "tx_me_daily": ("Tx daily ME (°C)", "RdBu_r", True, None), "tx_mae_daily": ("Tx daily MAE (°C)", "YlOrRd", False, None), "tx_rmse_daily": ("Tx daily RMSE (°C)", "YlOrRd", False, None), "tx_me_weekly_mean": ("Weekly-mean Tx ME (°C)", "RdBu_r", True, None), "tx_mae_weekly_mean": ("Weekly-mean Tx MAE (°C)", "YlOrRd", False, None), "tx_rmse_weekly_mean":("Weekly-mean Tx RMSE (°C)", "YlOrRd", False, None), "tx_me_weekly_max": ("Weekly-max Tx ME (°C)", "RdBu_r", True, None), "tx_mae_weekly_max": ("Weekly-max Tx MAE (°C)", "YlOrRd", False, None), "tx_rmse_weekly_max": ("Weekly-max Tx RMSE (°C)", "YlOrRd", False, None), "ehd_count_me": ("EHD day-count ME (days/week)", "RdBu_r", True, None), "ehd_count_mae": ("EHD day-count MAE (days/week)", "YlOrRd", False, None), "ehd_count_rmse": ("EHD day-count RMSE (days/week)", "YlOrRd", False, None), "base_rate": ("Observed EHD base rate", "YlOrRd", False, None), "pod": ("POD (target > 0.6)", "YlGnBu", False, (0, 1)), "pofd": ("POFD (target < 0.2)", "YlOrRd", False, (0, 0.2)), "far": ("FAR (target < 0.2)", "YlOrRd", False, (0, 1)), "fbi": ("FBI (target ~ 1)", "RdBu_r", False, (0, 2)), "csi": ("CSI (target > 0.4)", "YlGnBu", False, (0, 1)), "ets": ("ETS", "YlGnBu", False, None), "auc": ("ROC AUC (sweep on Sev_mod)", "YlGnBu", False, (0.5, 1)), } ORDER = list(SPEC) def load_regions(path=C.SHAPE_7REG): """Polygon parts of the seven sub-regions, as lists of (n,2) arrays.""" reader = shapefile.Reader(path) out = [] for shp in reader.shapes(): stops = list(shp.parts) + [len(shp.points)] for i in range(len(stops) - 1): out.append(np.asarray(shp.points[stops[i]:stops[i + 1]])) return out def _levels(data, symmetric, fixed): """Eleven contour levels covering the robust range of the whole row.""" if fixed is not None: lo, hi = fixed else: finite = data[np.isfinite(data)] if finite.size == 0: lo, hi = 0.0, 1.0 elif symmetric: hi = np.percentile(np.abs(finite), 98) hi = max(hi, 1e-6) lo = -hi else: lo = np.percentile(finite, 2) hi = np.percentile(finite, 98) if hi <= lo: hi = lo + 1e-6 return np.linspace(lo, hi, 11) def metric_row(field, grid, borders, title, cmap, symmetric, fixed): """One PNG: six maps in a row sharing a colour scale. field is (6, nlat, nlon).""" lev = _levels(field, symmetric, fixed) norm = BoundaryNorm(lev, ncolors=256, clip=False) fig, axes = plt.subplots(1, len(C.WEEKS), figsize=(13.2, 3.9), gridspec_kw=dict(wspace=0.02, left=0.012, right=0.90, top=0.90, bottom=0.04)) for k, ax in enumerate(axes): ax.set_facecolor("#f2f2f2") mesh = ax.pcolormesh(grid.lon, grid.lat, np.ma.masked_invalid(field[k]), cmap=cmap, norm=norm, shading="nearest") for poly in borders: ax.plot(poly[:, 0], poly[:, 1], color="0.25", lw=0.45) ax.set_xlim(grid.lon.min() - 0.125, grid.lon.max() + 0.125) ax.set_ylim(grid.lat.min() - 0.125, grid.lat.max() + 0.125) ax.set_aspect("equal") ax.set_xticks([]); ax.set_yticks([]) for s in ax.spines.values(): s.set_linewidth(0.5); s.set_color("0.6") ax.set_title(C.WEEK_NAMES[k], fontsize=11, pad=3) cax = fig.add_axes([0.915, 0.06, 0.014, 0.80]) cb = fig.colorbar(mesh, cax=cax, ticks=lev[::2]) cb.ax.tick_params(labelsize=8) cb.outline.set_linewidth(0.5) fig.suptitle(title, fontsize=12.5, y=0.985, x=0.012, ha="left") buf = io.BytesIO() fig.savefig(buf, format="png", dpi=125, facecolor="white") plt.close(fig) return base64.b64encode(buf.getvalue()).decode("ascii") PAGE = """ %(title)s

%(title)s

%(subtitle)s

%(body)s """ def main(): ap = argparse.ArgumentParser() ap.add_argument("--nc", default=os.path.join(C.OUT, "ehd_verification.nc")) ap.add_argument("--out", default=os.path.join(C.ROOT, "ehd_verification.html")) ap.add_argument("--title", default="EC S2S — EHD forecast verification (2003–2022)") args = ap.parse_args() grid = C.Grid() borders = load_regions() with Dataset(args.nc) as ds: avail = [m for m in ORDER if m in ds.variables] data = {m: np.ma.filled(ds[m][:], np.nan).astype("f4") for m in avail} years = getattr(ds, "years", "") body = [] # summary table: land-mean of every metric by week rows = ["" + "".join("" % w for w in C.WEEK_NAMES) + ""] for m in avail: v = [np.nanmean(data[m][k][grid.land2d]) for k in range(len(C.WEEKS))] rows.append("%s" % (escape(SPEC[m][0]), "".join("" % x for x in v))) rows.append("
Metric%s
%s
%.3f
") body.append("

Land-average of each metric by forecast week

") body.append("".join(rows)) for m in avail: title, cmap, sym, fixed = SPEC[m] png = metric_row(data[m], grid, borders, title, cmap, sym, fixed) body.append('

%s

%s
' % (escape(title), escape(title), png)) print(" row: %s" % m, flush=True) sub = ("Land points of Mask_VN_S2S_EC_0.25.nc (541 grid points), " "sub-region boundaries from Vietnam7.shp. " "EHD threshold Tx_thres = min(Tx90p, 35 °C), " "Tx90p from a 31-day window, computed separately for the model (per lead time) " "and for the observations. Weeks: W1 = lt 0–10, W2 = 11–17, " "W3 = 18–24, W4 = 25–31, W5 = 32–38, W6 = 39–45. " "Hindcast years: %s." % escape(str(years))) with open(args.out, "w") as fh: fh.write(PAGE % dict(title=escape(args.title), subtitle=sub, body="\n".join(body))) print("[done] -> %s" % args.out, flush=True) if __name__ == "__main__": main()