"""exercise_0b.py Exercise 0-B. The variogram and what its parameters do. Interpolates mean annual air temperature across the Arizona Meteorological Network and moves every parameter a commercial tool hides, scoring each change by leave one out cross validation. python exercise_0b.py """ from __future__ import annotations import sys from pathlib import Path import matplotlib import numpy as np matplotlib.use("Agg") import matplotlib.pyplot as plt # noqa: E402 sys.path.insert(0, str(Path(__file__).resolve().parent)) import common as C # noqa: E402 from kriging import (empirical_variogram, fit_variogram, grid, idw, # noqa loocv, model_function, ordinary_kriging, project_albers) VAR = "tmean_C" LAG_SWEEP = [40000, 80000, 150000] # meters SAMPLE_SIZES = [8, 12, 16, 20, 24] DIRECTIONS = [0, 45, 90, 135] def load(): d = C.load_azmet().dropna(subset=[VAR, "lat", "lon"]).reset_index(drop=True) x, y = project_albers(d.lon.to_numpy(), d.lat.to_numpy()) return d, x, y, d[VAR].to_numpy(float) def main(): C.ensure_dirs() d, x, y, z = load() C.banner("exercise 0-B: the variogram and its parameters", "AZMet daily observations, 2024, keyless", "azmet_2024_annual.csv", {"stations": len(d), "variable": VAR}) out = {"stations": int(len(d)), "elev_min_m": int(d.elev_m.min()), "elev_max_m": int(d.elev_m.max())} # ---- 2. the empirical variogram ------------------------------------ cen, val, cnt, lag, npairs = empirical_variogram(x, y, z, n_bins=12) p = fit_variogram(cen, val, cnt) print("\n[2] EMPIRICAL VARIOGRAM") print(f" pairs {npairs} bins {len(cen)} lag {lag/1000:.1f} km") print(f" fitted nugget {p['nugget']:.3f} sill {p['sill']:.3f} " f"range {p['rng']/1000:.1f} km") base = loocv(x, y, z, p) print(f" LOOCV rmse {base['rmse']:.3f} C mae {base['mae']:.3f} C") out["pairs"] = npairs out["base"] = {"lag_km": round(lag/1000, 1), "bins": int(len(cen)), "nugget": round(p["nugget"], 3), "sill": round(p["sill"], 3), "range_km": round(p["rng"]/1000, 1), "rmse": round(base["rmse"], 3), "mae": round(base["mae"], 3)} # ---- 3. lag and bin count ------------------------------------------ print("\n[3] LAG AND BIN COUNT") print(f" {'lag km':>8}{'bins':>7}{'nugget':>9}{'sill':>8}" f"{'range km':>10}{'rmse':>8}") out["lag_sweep"] = [] for L in LAG_SWEEP: c2, v2, n2, _, _ = empirical_variogram(x, y, z, lag=L) if len(c2) < 3: continue q = fit_variogram(c2, v2, n2) cv = loocv(x, y, z, q) print(f" {L/1000:>8.0f}{len(c2):>7d}{q['nugget']:>9.3f}" f"{q['sill']:>9.3f}{q['rng']/1000:>11.1f}{cv['rmse']:>9.3f}") out["lag_sweep"].append({"lag_km": L/1000, "bins": int(len(c2)), "range_km": round(q["rng"]/1000, 1), "rmse": round(cv["rmse"], 3)}) # ---- 4. nugget, sill, range overridden one at a time ---------------- print("\n[4] OVERRIDING ONE PARAMETER AT A TIME") overrides = [("nugget = 0", dict(nugget=0.0)), ("nugget = sill (pure nugget)", dict(nugget=p["sill"] * 0.999, sill=p["sill"])), ("range halved", dict(rng=p["rng"] * 0.5)), ("range doubled", dict(rng=p["rng"] * 2.0))] out["overrides"] = [] print(f" {'setting':<30}{'rmse':>8}{'delta':>9}") for name, kw in overrides: q = fit_variogram(cen, val, cnt, **kw) cv = loocv(x, y, z, q) print(f" {name:<30}{cv['rmse']:>8.3f}" f"{cv['rmse']-base['rmse']:>+9.3f}") out["overrides"].append({"setting": name, "rmse": round(cv["rmse"], 3), "delta": round(cv["rmse"]-base["rmse"], 3)}) # ---- 5. direction --------------------------------------------------- print("\n[5] DIRECTIONAL VARIOGRAMS") print(f" {'bearing':>8}{'bins':>7}{'range km':>10}{'sill':>9}") out["directional"] = [] for a in DIRECTIONS: c3, v3, n3, _, _ = empirical_variogram(x, y, z, n_bins=8, angle=a) if len(c3) < 3: print(f" {a:>8d}{len(c3):>7d} too few pairs to fit") continue q = fit_variogram(c3, v3, n3) print(f" {a:>8d}{len(c3):>7d}{q['rng']/1000:>10.1f}" f"{q['sill']:>9.3f}") out["directional"].append({"bearing": a, "range_km": round(q["rng"]/1000, 1), "sill": round(q["sill"], 3)}) # ---- 6. sample size ------------------------------------------------- print("\n[6] SAMPLE SIZE") rng_ = np.random.default_rng(C.SEED) print(f" {'n':>4}{'rmse':>9}{'range km':>11} (median of 15 draws)") out["sample_size"] = [] for n in SAMPLE_SIZES: if n >= len(z): continue skipped = 0 rr, gg = [], [] for _ in range(15): k = rng_.choice(len(z), n, replace=False) c4, v4, n4, _, _ = empirical_variogram(x[k], y[k], z[k], n_bins=6) if len(c4) < 3: skipped += 1 continue q = fit_variogram(c4, v4, n4) rr.append(loocv(x[k], y[k], z[k], q)["rmse"]) gg.append(q["rng"] / 1000) if not rr: print(f" {n:>4d} no draw produced enough bins to fit") out["sample_size"].append({"n": n, "rmse": None, "note": "unfittable at this n"}) else: note = f" ({skipped} of 15 draws unfittable)" if skipped else "" print(f" {n:>4d}{np.median(rr):>9.3f}{np.median(gg):>11.1f}{note}") out["sample_size"].append({"n": n, "rmse": round(float(np.median(rr)), 3), "range_km": round(float(np.median(gg)), 1)}) # ---- 7. kriging against inverse distance weighting ------------------ idw_cv = loocv(x, y, z, method="idw", power=2.0) print("\n[7] KRIGING AGAINST INVERSE DISTANCE WEIGHTING") print(f" ordinary kriging rmse {base['rmse']:.3f} mae {base['mae']:.3f}") print(f" IDW power 2 rmse {idw_cv['rmse']:.3f} " f"mae {idw_cv['mae']:.3f}") out["idw"] = {"rmse": round(idw_cv["rmse"], 3), "mae": round(idw_cv["mae"], 3)} # ---- figures --------------------------------------------------------- XX, YY, gx, gy = grid(x, y) kp, kv = ordinary_kriging(x, y, z, XX.ravel(), YY.ravel(), p) iw = idw(x, y, z, XX.ravel(), YY.ravel()) src = (f"Source: Arizona Meteorological Network, 2024 daily observations. " f"{len(d)} stations. EPSG:5070.") fig, ax = plt.subplots(figsize=(8, 5.5)) ax.plot(cen / 1000, val, "o", color="#1f3864", label="binned semivariance") hh = np.linspace(0, cen.max(), 250) f = model_function(p["model"]) ax.plot(hh / 1000, f(hh, p["nugget"], p["sill"], p["rng"]), "-", color="#b2182b", lw=2, label=f"fitted {p['model']}") ax.axhline(p["sill"], ls=":", color="#666") ax.axhline(p["nugget"], ls=":", color="#666") ax.axvline(p["rng"] / 1000, ls=":", color="#666") ax.annotate("sill", (0.02, p["sill"]), xycoords=("axes fraction", "data"), va="bottom", fontsize=9, color="#444") ax.annotate("nugget", (0.02, p["nugget"]), xycoords=("axes fraction", "data"), va="bottom", fontsize=9, color="#444") ax.annotate("range", (p["rng"] / 1000, 0.02), xycoords=("data", "axes fraction"), fontsize=9, color="#444") ax.set_xlabel("separation distance (km)") ax.set_ylabel("semivariance (degrees C squared)") ax.set_title("Anatomy of a variogram", loc="left", fontweight="bold") ax.legend(frameon=False, fontsize=9) ax.text(0, -0.14, src, transform=ax.transAxes, fontsize=7, color="#666") fig.savefig(C.FIG / "exercise_0b_variogram.png", dpi=150, bbox_inches="tight") plt.close(fig) fig, axes = plt.subplots(1, 3, figsize=(17, 5.2)) for ax, dat, t in ((axes[0], kp, "Ordinary kriging, prediction"), (axes[1], np.sqrt(kv), "Kriging standard error"), (axes[2], iw, "Inverse distance weighting")): im = ax.pcolormesh(gx / 1000, gy / 1000, dat.reshape(XX.shape), cmap="magma" if "error" in t else "RdYlBu_r", shading="auto") ax.scatter(x / 1000, y / 1000, s=12, c="k", zorder=3) ax.set_aspect("equal") ax.set_title(t, loc="left", fontweight="bold", fontsize=11) ax.set_xticks([]), ax.set_yticks([]) fig.colorbar(im, ax=ax, shrink=0.8) axes[0].text(0, -0.06, src, transform=axes[0].transAxes, fontsize=7, color="#666") fig.suptitle("The prediction, its uncertainty, and the method that " "reports none", fontsize=13, fontweight="bold") fig.tight_layout() fig.savefig(C.FIG / "exercise_0b_surfaces.png", dpi=140, bbox_inches="tight") plt.close(fig) # ---- handoff --------------------------------------------------------- import pandas as pd tbl = d[["station", "name", "elev_m", "lat", "lon", VAR]].copy() tbl["loocv_pred"] = base["pred"] tbl["loocv_error"] = base["error"] C.write_handoff("exercise_0b", tbl, { "model_form": p["model"], "lag_m": round(lag, 1), "bin_count": int(len(cen)), "nugget": round(p["nugget"], 4), "sill": round(p["sill"], 4), "range_m": round(p["rng"], 1), "anisotropy": "none assumed; the fit pooled all directions", "search_neighborhood": "all stations, no maximum", "grid_resolution": f"{XX.shape[1]} by {XX.shape[0]}", "color_ramp": "RdYlBu_r for prediction, magma for error, unjustified", "projection": "EPSG:5070", }, [ "Is a spherical model right for this phenomenon, or does the shape of " "the empirical cloud favor exponential?", "The fit pooled all directions. The directional results in Block 5 " "say whether that assumption survives.", "Should the surface be drawn beyond the convex hull of the stations " "at all, given the error surface there?", "Does the color ramp on the error surface imply a threshold that " "nobody set?", "Which stations are doing the work, and would losing any one of them " "change the published surface?", ]) C.save_checkvalues("exercise_0b", out) print(f"\nwrote {C.FIG / 'exercise_0b_variogram.png'}") print(f"wrote {C.FIG / 'exercise_0b_surfaces.png'}") if __name__ == "__main__": main()