"""How to use the trained Taiwan intensity models. Three worked examples on one hypothetical earthquake: 1. predict at a single site 2. predict at every real P-Alert station and rank the strongest shaking 3. render a ShakeMap PNG over the Taiwan grid Run: python use_model.py python use_model.py --M 7.0 --depth 15 --ev-lat 23.9 --ev-lon 121.6 python use_model.py --shakemap data/shakemap.png Both models are pure functions of (M, depth, ev_lat, ev_lon, t_lat, t_lon) and are loaded from ONNX, so this needs only numpy + onnxruntime + matplotlib. """ import argparse import json import os import sys import numpy as np from palert.features import INTENSITY_LEVELS from palert.paths import DATA_DIR, RAW_DIR from palert.predict import IntensityClassifier, IntensityGMM STATIONS_PATH = os.path.join(RAW_DIR, "stations.json") COUNTIES_PATH = os.path.join(DATA_DIR, "twcounty.geojson") # DejaVu (matplotlib default) has no CJK glyphs, so plots use ASCII labels ASCII_LABELS = ["0", "1", "2", "3", "4", "5-", "5+", "6-", "6+", "7"] def load_stations(): """Real P-Alert site coordinates (skips entries missing lon/lat).""" with open(STATIONS_PATH, encoding="utf-8") as fh: stas = json.load(fh) stas = [s for s in stas if s.get("lon") is not None and s.get("lat") is not None] return stas def example_single(gmm, clf, M, depth, ev_lat, ev_lon): """1. One target site, a fixed offset from the epicenter.""" t_lat, t_lon = ev_lat + 0.1, ev_lon + 0.1 pga, pgv, level, label = gmm.predict(M, depth, ev_lat, ev_lon, t_lat, t_lon) print("[1] single site (regression GMM)") print(f" site ({t_lat:.4f}, {t_lon:.4f})") print(f" PGA = {pga:8.3f} gal PGV = {pgv:6.3f} cm/s intensity = {label} ({level})") clevel, clabel, probs = clf.predict(M, depth, ev_lat, ev_lon, t_lat, t_lon) print(f" classifier: intensity = {clabel} ({clevel}), p = {probs[clevel]:.3f}") def example_stations(gmm, M, depth, ev_lat, ev_lon, top=10): """2. Batch prediction at every real station, ranked by intensity then PGA.""" stas = load_stations() t_lat = np.array([s["lat"] for s in stas], dtype=np.float64) t_lon = np.array([s["lon"] for s in stas], dtype=np.float64) pga, pgv, level = gmm.predict_batch(M, depth, ev_lat, ev_lon, t_lat, t_lon) order = np.lexsort((-pga, -level)) print(f"[2] strongest shaking across {len(stas)} P-Alert stations (top {top})") print(f" {'code':>6} {'area':<12} {'PGA(gal)':>9} {'PGV(cm/s)':>10} intensity") for i in order[:top]: name = INTENSITY_LEVELS[int(level[i])] area = stas[i].get("area", "") or "" print(f" {stas[i]['station']:>6} {area:<12} {pga[i]:9.3f} {pgv[i]:10.3f} {name}") print(f" max intensity = {INTENSITY_LEVELS[int(level.max())]}, " f"stations at or above 5弱 = {int((level >= 5).sum())}") def plot_counties(ax, path=COUNTIES_PATH, color="k", linewidth=0.6): """Draw Taiwan county boundaries as black outlines (skips if file missing).""" if not os.path.exists(path): print(f" (county borders skipped: {path} not found)") return from matplotlib.collections import LineCollection with open(path, encoding="utf-8") as fh: gj = json.load(fh) rings = [] for feat in gj.get("features", []): geom = feat.get("geometry") or {} coords = geom.get("coordinates") if not coords: continue polys = coords if geom.get("type") == "MultiPolygon" else [coords] for poly in polys: for ring in poly: arr = np.asarray(ring, dtype=np.float64) if arr.shape[0] >= 2: rings.append(arr[:, :2]) ax.add_collection(LineCollection(rings, colors=color, linewidths=linewidth, zorder=6)) def example_shakemap(gmm, M, depth, ev_lat, ev_lon, out_path): """3. Render a ShakeMap grid over Taiwan and save it as a PNG.""" import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt from matplotlib.colors import ListedColormap, BoundaryNorm lats = np.arange(21.7, 25.5, 0.02) lons = np.arange(119.8, 122.2, 0.02) lon_grid, lat_grid = np.meshgrid(lons, lats) _, _, level = gmm.predict_batch(M, depth, ev_lat, ev_lon, lat_grid, lon_grid) level = level.reshape(lat_grid.shape) colors = ["#ffffff", "#c8e6f5", "#a5d8f3", "#7fc9ee", "#57b6e5", "#3b9ede", "#f5d000", "#f5a000", "#f56000", "#d00000"] cmap = ListedColormap(colors) norm = BoundaryNorm(np.arange(-0.5, 10.5, 1), cmap.N) fig, ax = plt.subplots(figsize=(7, 9), dpi=130) mesh = ax.pcolormesh(lon_grid, lat_grid, level, cmap=cmap, norm=norm, shading="auto") ax.plot(ev_lon, ev_lat, marker="*", color="k", markersize=16, label="epicenter") plot_counties(ax) ax.set_xlim(lons[0], lons[-1]) ax.set_ylim(lats[0], lats[-1]) ax.set_xlabel("longitude") ax.set_ylabel("latitude") ax.set_title(f"Predicted CWA intensity | M{M} depth {depth} km @ ({ev_lat:.3f}, {ev_lon:.3f})") ax.set_aspect("equal") cb = fig.colorbar(mesh, ax=ax, ticks=range(10), shrink=0.7) cb.ax.set_yticklabels(ASCII_LABELS) cb.set_label("intensity") ax.legend(loc="lower left") fig.tight_layout() fig.savefig(out_path) print(f"[3] shakemap saved -> {out_path}") def main(): try: sys.stdout.reconfigure(encoding="utf-8") except Exception: pass ap = argparse.ArgumentParser() ap.add_argument("--M", type=float, default=6.5) ap.add_argument("--depth", type=float, default=10.0) ap.add_argument("--ev-lat", type=float, default=23.9, dest="ev_lat") ap.add_argument("--ev-lon", type=float, default=121.6, dest="ev_lon") ap.add_argument("--gmm", default=None, help="path to GMM ONNX (default: TDN-IntensityGMM-TW-524k-v1.onnx)") ap.add_argument("--clf", default=None, help="path to classifier ONNX") ap.add_argument("--shakemap", default=None, help="if set, write a ShakeMap PNG here") args = ap.parse_args() gmm = IntensityGMM(args.gmm) if args.gmm else IntensityGMM() clf = IntensityClassifier(args.clf) if args.clf else IntensityClassifier() print(f"event: M{args.M} depth {args.depth} km epicenter ({args.ev_lat}, {args.ev_lon})\n") example_single(gmm, clf, args.M, args.depth, args.ev_lat, args.ev_lon) print() example_stations(gmm, args.M, args.depth, args.ev_lat, args.ev_lon) if args.shakemap: print() example_shakemap(gmm, args.M, args.depth, args.ev_lat, args.ev_lon, args.shakemap) if __name__ == "__main__": main()