| """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") |
|
|
| |
| 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() |
|
|
