07 — Facet Comparison: Blocking vs Propagating Anticyclones (PRP)
Compares PV-budget grouped-term decomposition between blocking and propagating anticyclone (PRP) events using bootstrap significance testing and orthogonal basis projection.
Config options:
LEVEL:"wavg"for pre-computed wavg fields, or an integer hPa (e.g.250) for a single isobaric level from the_3darrays.When
LEVEL="wavg": LHS tendency =pv_anom_dt + pv_bar_dt, RHS terms use 2-D wavg fields.When
LEVEL=<int>: LHS tendency =pv_dt_3d[k], RHS terms use_3d[k]arrays directly.
Data:
Blocking:
/net/flood/data2/users/x_yan/composite_blocking_tempestPRP:
/net/flood/data2/users/x_yan/composite_prp_tempest
[1]:
import numpy as np
import matplotlib.pyplot as plt
import os, glob
from concurrent.futures import ThreadPoolExecutor
from zipfile import BadZipFile
from pvtend import compute_orthogonal_basis, project_field
from pvtend.plotting import plot_var, plot_baroclinic_tilt
from pvtend.decomposition.smoothing import gaussian_smooth_nan
1 Configuration
[2]:
# ═══════════════════════════════════════════════════════════
# Configuration
# ═══════════════════════════════════════════════════════════
DATA_ROOTS = {
"blocking": "/net/flood/data2/users/x_yan/pvtend/outputs/blocking",
"prp": "/net/flood/data2/users/x_yan/pvtend/outputs/prp",
}
STAGE = "peak"
DH = 0
DH_BASIS = max(DH - 1, -13)
SMOOTH_DEG = 3.0 # 3° Gaussian smoothing
GRID_SP = 1.5
N_BOOT = 200
ALPHA = 0.05
SEED = 42
# Level: "wavg" for pre-computed weighted-average, or int hPa (e.g. 250)
LEVEL = 250 #"wavg
# ═══════════════════════════════════════════════════════════
# Projection masking config (toggle on/off)
# ═══════════════════════════════════════════════════════════
USE_SIG_MASK = True # zero out non-significant pixels before projection
MASK_SPEC = "< 0" # restrict basis to q' < 0 region (mask)
# ═══════════════════════════════════════════════════════════
# I/O helpers
# ═══════════════════════════════════════════════════════════
def _load_npz(path):
try:
return dict(np.load(path))
except (BadZipFile, EOFError, OSError):
return None
def load_events(etype, dh, stage=STAGE):
"""Load all events for an event type at a given dh."""
sign = "+" if dh >= 0 else ""
d = f"{DATA_ROOTS[etype]}/{stage}/dh={sign}{dh}"
files = sorted(glob.glob(os.path.join(d, "track_*.npz")))
if not files:
return []
with ThreadPoolExecutor(max_workers=8) as pool:
results = list(pool.map(_load_npz, files))
good = [r for r in results if r is not None]
n_bad = len(results) - len(good)
if n_bad:
print(f" ⚠ {etype} dh={dh}: skipped {n_bad} corrupt NPZ files")
return good
def get_field(event, name, level=LEVEL):
"""Extract a 2-D field, respecting LEVEL config.
If LEVEL="wavg", returns the pre-computed 2-D wavg field.
If LEVEL=<int hPa>, slices the _3d array at that level.
"""
if isinstance(level, str) and level.lower() == "wavg":
return event[name]
else:
key_3d = name if name.endswith("_3d") else name + "_3d"
arr_3d = event[key_3d]
lvls = event["levels"]
idx = int(np.argmin(np.abs(lvls - float(level))))
return arr_3d[idx]
def get_lhs_tendency(event, level=LEVEL):
"""Return full-field dq/dt (LHS tendency).
wavg: pv_dt
3D: pv_dt_3d[k] (full field exists)
"""
if isinstance(level, str) and level.lower() == "wavg":
return event["pv_dt"]
else:
lvls = event["levels"]
idx = int(np.argmin(np.abs(lvls - float(level))))
return event["pv_dt_3d"][idx]
smooth = lambda f: gaussian_smooth_nan(f, smoothing_deg=SMOOTH_DEG, grid_spacing=GRID_SP)
print(f"Stage={STAGE} DH={DH} DH_BASIS={DH_BASIS} "
f"LEVEL={LEVEL} Smooth={SMOOTH_DEG}° N_BOOT={N_BOOT}")
print(f" USE_SIG_MASK={USE_SIG_MASK} USE_NEG_PV_MASK={MASK_SPEC}")
Stage=peak DH=0 DH_BASIS=-1 LEVEL=250 Smooth=3.0° N_BOOT=200
USE_SIG_MASK=True USE_NEG_PV_MASK=< 0
2 Load events & build composite basis
[3]:
# ── Load events for both types ──
events = {}
events_basis = {}
for etype in DATA_ROOTS:
events[etype] = load_events(etype, DH)
events_basis[etype] = load_events(etype, DH_BASIS)
print(f"{etype}: N={len(events[etype])} (dh={DH}), "
f"N_basis={len(events_basis[etype])} (dh={DH_BASIS})")
# Grid coordinates (same for both types)
X_rel = events["blocking"][0]["X_rel"]
Y_rel = events["blocking"][0]["Y_rel"]
x_rel = X_rel[0, :]
y_rel = Y_rel[:, 0]
# ── Build composite-mean basis from dh-1 for each type ──
basis = {}
for etype in DATA_ROOTS:
# evs_b = events_basis[etype]
# pv_b = np.nanmean([e["pv_anom"] for e in evs_b], axis=0)
# dx_b = np.nanmean([e["pv_dx"] for e in evs_b], axis=0)
# dy_b = np.nanmean([e["pv_dy"] for e in evs_b], axis=0)
# dh composite means (current-time positional fields)
evs = events[etype]
pv_n = np.nanmean([e["pv_anom"] for e in evs], axis=0)
dx_n = np.nanmean([e["pv_dx"] for e in evs], axis=0)
dy_n = np.nanmean([e["pv_dy"] for e in evs], axis=0)
dx_dx_n = np.nanmean([e["pv_dx_dx"] for e in evs], axis=0)
dy_dy_n = np.nanmean([e["pv_dy_dy"] for e in evs], axis=0)
dx_dy_n = np.nanmean([e["pv_dx_dy"] for e in evs], axis=0)
basis[etype] = compute_orthogonal_basis(
pv_n, dx_n, dy_n, x_rel, y_rel,
mask=MASK_SPEC,
apply_smoothing=True, smoothing_deg=SMOOTH_DEG, grid_spacing=GRID_SP,
# pv_anom_prev=pv_b, pv_dx_prev=dx_b, pv_dy_prev=dy_b,
pv_dx_dx_prev=dx_dx_n, pv_dy_dy_prev=dy_dy_n, pv_dx_dy_prev=dx_dy_n
)
print(f"{etype} basis norms: "
+ " ".join(f"{k}={v:.3e}" for k, v in basis[etype].norms.items()))
blocking: N=1260 (dh=0), N_basis=1260 (dh=-1)
prp: N=2097 (dh=0), N_basis=2096 (dh=-1)
blocking basis norms: beta=2.699e+02 ax=2.959e+02 ay=1.653e+02 gamma1=1.120e+02 gamma2=6.998e+01 sigma=4.141e+01
prp basis norms: beta=1.975e+02 ax=1.749e+02 ay=1.061e+02 gamma1=5.176e+01 gamma2=5.184e+01 sigma=1.812e+01
/tmp/ipykernel_1373389/2842660874.py:33: UserWarning: compute_orthogonal_basis: grid_spacing=1.5°, center_lat=60.0°N → dx(center)=83.4 km, dy=166.8 km
basis[etype] = compute_orthogonal_basis(
3 Bootstrap helpers
[4]:
def bootstrap_field(events_list, field_func, n_boot=N_BOOT, alpha=ALPHA, seed=SEED):
"""Bootstrap composite mean + significance mask for a field.
Returns (mean_2d, sig_mask) where sig_mask=True for significant pixels.
"""
stack = np.array([field_func(e) for e in events_list]) # (N, Y, X)
N = stack.shape[0]
rng = np.random.default_rng(seed)
boot = np.empty((n_boot, *stack.shape[1:]))
for b in range(n_boot):
idx = rng.integers(0, N, size=N)
boot[b] = np.nanmean(stack[idx], axis=0)
lo = np.nanpercentile(boot, 100 * alpha / 2, axis=0)
hi = np.nanpercentile(boot, 100 * (1 - alpha / 2), axis=0)
mean = np.nanmean(stack, axis=0)
sig_mask = ~((lo <= 0) & (hi >= 0))
return mean, sig_mask
def bootstrap_project(events_list, field_func, basis_obj,
n_boot=N_BOOT, seed=SEED):
"""Bootstrap the projected coefficients (sig-only field).
Per-event: zero out non-sig pixels, smooth, project.
Bootstrap the mean coefficients.
Returns dict with keys beta, ax, ay, gamma1, gamma2, sigma — each (n_boot,).
"""
# First get composite sig mask
_, sig_mask = bootstrap_field(events_list, field_func)
# Per-event projections using sig-only composite field
per_event = []
for e in events_list:
raw = field_func(e)
if USE_SIG_MASK:
sig_field = np.where(sig_mask, raw, 0.0)
else:
sig_field = raw.copy()
sig_s = smooth(sig_field)
p = project_field(sig_s, basis_obj)
per_event.append({k: p[k] for k in ["beta", "ax", "ay", "gamma1", "gamma2", "sigma"]})
per_arr = {k: np.array([pe[k] for pe in per_event])
for k in ["beta", "ax", "ay", "gamma1", "gamma2", "sigma"]}
N = len(events_list)
rng = np.random.default_rng(seed)
boot = {k: np.empty(n_boot) for k in per_arr}
for b in range(n_boot):
idx = rng.integers(0, N, size=N)
for k in boot:
boot[k][b] = per_arr[k][idx].mean()
return boot
print("Bootstrap helpers defined.")
Bootstrap helpers defined.
4 Define grouped terms for cell 6 (4-group) and cell 7 (7-term)
[5]:
# ═══════════════════════════════════════════════════════════
# 4-GROUP decomposition (Cell 6 bar plot)
# ═══════════════════════════════════════════════════════════
GROUPS_4 = {
r"$-$Lead. adv": lambda e: -(
get_field(e, "u_rot_bar_pv_anom_dx") + get_field(e, "v_rot_anom_pv_bar_dy")
),
r"$-$Eddy rot": lambda e: -(
get_field(e, "u_rot_anom_pv_anom_dx") + get_field(e, "v_rot_anom_pv_anom_dy")
),
r"$-u_{\chi, \text{adiabatic}} \nabla q'$": lambda e: -(
get_field(e, "u_div_adiabatic_pv_anom_dx") + get_field(e, "v_div_adiabatic_pv_anom_dy")
),
r"$-u_{\chi, \text{diabatic}} \nabla q'$": lambda e: -(
get_field(e, "u_div_diabatic_pv_anom_dx") + get_field(e, "v_div_diabatic_pv_anom_dy")
),
}
# ═══════════════════════════════════════════════════════════
# 7-TERM decomposition (Cell 7 bar plot)
# ═══════════════════════════════════════════════════════════
GROUPS_7 = {
r"$-\bar{u}_\Psi\,q'_x$": lambda e: -get_field(e, "u_rot_bar_pv_anom_dx"),
r"$-v'_\Psi\,\bar{q}_y$": lambda e: -get_field(e, "v_rot_anom_pv_bar_dy"),
r"$-$Eddy rot": lambda e: -(
get_field(e, "u_rot_anom_pv_anom_dx") + get_field(e, "v_rot_anom_pv_anom_dy")
),
r"$-$Eddy adiabatic div": lambda e: -(
get_field(e, "u_div_adiabatic_pv_anom_dx") + get_field(e, "v_div_adiabatic_pv_anom_dy")
),
r"$-$Eddy diabatic div": lambda e: -(
get_field(e, "u_div_diabatic_pv_anom_dx") + get_field(e, "v_div_diabatic_pv_anom_dy")
),
r"$-\omega\,\bar{q}_p$": lambda e: -(
get_field(e, "w_adiabatic_pv_bar_dp") + get_field(e, "w_diabatic_pv_bar_dp")
),
r"$Q$": lambda e: (
get_field(e, "Q")
),
}
print(f"4-group terms: {list(GROUPS_4.keys())}")
print(f"7-term terms: {list(GROUPS_7.keys())}")
4-group terms: ['$-$Lead. adv', '$-$Eddy rot', "$-u_{\\chi, \\text{adiabatic}} \\nabla q'$", "$-u_{\\chi, \\text{diabatic}} \\nabla q'$"]
7-term terms: ["$-\\bar{u}_\\Psi\\,q'_x$", "$-v'_\\Psi\\,\\bar{q}_y$", '$-$Eddy rot', '$-$Eddy adiabatic div', '$-$Eddy diabatic div', '$-\\omega\\,\\bar{q}_p$', '$Q$']
5 Run bootstrap for all terms (both event types)
[6]:
# Also bootstrap the LHS tendency (pv_dt)
ALL_TERMS = {"dq/dt": get_lhs_tendency}
ALL_TERMS.update(GROUPS_4)
print("═══ Bootstrap for 4-group decomposition ═══")
boot_4 = {} # boot_4[etype][term_name] = {beta: (N_BOOT,), ...}
for etype in DATA_ROOTS:
boot_4[etype] = {}
evs = events[etype]
b = basis[etype]
for tname, tfunc in ALL_TERMS.items():
boot_4[etype][tname] = bootstrap_project(evs, tfunc, b)
lo, hi = np.nanpercentile(boot_4[etype][tname]["ax"], [2.5, 97.5])
sig = "***" if lo * hi > 0 else "n.s."
print(f" {etype:10s} {tname:35s} αx 95%CI: [{lo:.3e}, {hi:.3e}] {sig}")
print("\n═══ Bootstrap for 7-term decomposition ═══")
boot_7 = {}
ALL_TERMS_7 = {"dq/dt": get_lhs_tendency}
ALL_TERMS_7.update(GROUPS_7)
for etype in DATA_ROOTS:
boot_7[etype] = {}
evs = events[etype]
b = basis[etype]
for tname, tfunc in ALL_TERMS_7.items():
boot_7[etype][tname] = bootstrap_project(evs, tfunc, b)
lo, hi = np.nanpercentile(boot_7[etype][tname]["ax"], [2.5, 97.5])
sig = "***" if lo * hi > 0 else "n.s."
print(f" {etype:10s} {tname:35s} αx 95%CI: [{lo:.3e}, {hi:.3e}] {sig}")
print("\nDone.")
═══ Bootstrap for 4-group decomposition ═══
blocking dq/dt αx 95%CI: [1.787e-01, 3.937e-01] ***
blocking $-$Lead. adv αx 95%CI: [6.164e+00, 7.690e+00] ***
blocking $-$Eddy rot αx 95%CI: [-3.305e-01, -5.820e-02] ***
blocking $-u_{\chi, \text{adiabatic}} \nabla q'$ αx 95%CI: [-7.527e-01, -5.935e-01] ***
blocking $-u_{\chi, \text{diabatic}} \nabla q'$ αx 95%CI: [-8.040e-01, -5.852e-01] ***
prp dq/dt αx 95%CI: [1.309e+00, 1.830e+00] ***
prp $-$Lead. adv αx 95%CI: [8.970e+00, 1.154e+01] ***
prp $-$Eddy rot αx 95%CI: [4.993e-01, 9.183e-01] ***
prp $-u_{\chi, \text{adiabatic}} \nabla q'$ αx 95%CI: [-8.538e-01, -6.429e-01] ***
prp $-u_{\chi, \text{diabatic}} \nabla q'$ αx 95%CI: [-1.121e+00, -7.810e-01] ***
═══ Bootstrap for 7-term decomposition ═══
blocking dq/dt αx 95%CI: [1.787e-01, 3.937e-01] ***
blocking $-\bar{u}_\Psi\,q'_x$ αx 95%CI: [1.783e+01, 2.035e+01] ***
blocking $-v'_\Psi\,\bar{q}_y$ αx 95%CI: [-1.181e+01, -1.052e+01] ***
blocking $-$Eddy rot αx 95%CI: [-3.305e-01, -5.820e-02] ***
blocking $-$Eddy adiabatic div αx 95%CI: [-7.527e-01, -5.935e-01] ***
blocking $-$Eddy diabatic div αx 95%CI: [-8.040e-01, -5.852e-01] ***
blocking $-\omega\,\bar{q}_p$ αx 95%CI: [-4.151e+00, -3.458e+00] ***
blocking $Q$ αx 95%CI: [-6.404e-02, -4.759e-02] ***
prp dq/dt αx 95%CI: [1.309e+00, 1.830e+00] ***
prp $-\bar{u}_\Psi\,q'_x$ αx 95%CI: [2.206e+01, 2.637e+01] ***
prp $-v'_\Psi\,\bar{q}_y$ αx 95%CI: [-1.438e+01, -1.226e+01] ***
prp $-$Eddy rot αx 95%CI: [4.993e-01, 9.183e-01] ***
prp $-$Eddy adiabatic div αx 95%CI: [-8.538e-01, -6.429e-01] ***
prp $-$Eddy diabatic div αx 95%CI: [-1.121e+00, -7.810e-01] ***
prp $-\omega\,\bar{q}_p$ αx 95%CI: [-6.217e+00, -5.053e+00] ***
prp $Q$ αx 95%CI: [-1.202e-01, -8.168e-02] ***
Done.
6 Bar plot — 4-group comparison (Blocking vs PRP)
[7]:
fig, axes = plt.subplots(2, 3, figsize=(16, 10))
coef_names = ["beta", "ax", "ay", "gamma1", "gamma2", "sigma"]
coef_labels = [r"$\beta$ (intensification) [s$^{-1}$]",
r"$\alpha_x$ (zonal propagation) [m/s]",
r"$\alpha_y$ (meridional propagation) [m/s]",
r"$\gamma_1$ (shear def.) [km$^2$/s]",
r"$\gamma_2$ (strain def.) [km$^2$/s]",
r"$\sigma$ (Laplacian) [km$^2$/s]"]
term_names_4 = list(ALL_TERMS.keys())
n_terms = len(term_names_4)
x_pos = np.arange(n_terms)
width = 0.35
colors_etype = {"blocking": "steelblue", "prp": "coral"}
for ax, cname, clabel in zip(axes.flat, coef_names, coef_labels):
bar_info = {}
# ── Draw bars ──
for i_et, etype in enumerate(["blocking", "prp"]):
means, los, his, sigs = [], [], [], []
for tname in term_names_4:
vals = boot_4[etype][tname][cname]
m = vals.mean()
lo, hi = np.nanpercentile(vals, [2.5, 97.5])
means.append(m)
los.append(m - lo)
his.append(hi - m)
sigs.append(lo * hi > 0)
bar_info[etype] = (means, los, his, sigs)
offset = -width / 2 if i_et == 0 else width / 2
ax.bar(x_pos + offset, means, width,
yerr=[los, his], capsize=4,
label=etype.capitalize(),
color=colors_etype[etype], alpha=0.85)
# ── Set y-limits from actual bar extents (not forced 0-centered) ──
ymin_all, ymax_all = 0, 0
for etype in ["blocking", "prp"]:
ms, ls, hs, _ = bar_info[etype]
for m, l, h in zip(ms, ls, hs):
ymax_all = max(ymax_all, m + h)
ymin_all = min(ymin_all, m - l)
pad = (ymax_all - ymin_all) * 0.15
ax.set_ylim(ymin_all - pad, ymax_all + pad)
# ── Significance stars ──
ylim = ax.get_ylim()
for i_et, etype in enumerate(["blocking", "prp"]):
means, los, his, sigs = bar_info[etype]
offset = -width / 2 if i_et == 0 else width / 2
for j, (s, m) in enumerate(zip(sigs, means)):
if s:
y_off = his[j] if m >= 0 else -los[j]
y_star = m + y_off * 1.1 + np.sign(m) * abs(m) * 0.05
y_star = np.clip(y_star, ylim[0] * 0.95, ylim[1] * 0.95)
ax.text(x_pos[j] + offset, y_star,
"***", ha="center", fontsize=7, color="red",
fontweight="bold", clip_on=True)
# ── Between-group significance: blocking − PRP ──
yr = ylim[1] - ylim[0]
for j, tname in enumerate(term_names_4):
diff = boot_4["blocking"][tname][cname] - boot_4["prp"][tname][cname]
lo_d, hi_d = np.nanpercentile(diff, [2.5, 97.5])
if lo_d * hi_d > 0:
top_j = max(
bar_info["blocking"][0][j] + bar_info["blocking"][2][j],
bar_info["prp"][0][j] + bar_info["prp"][2][j])
bot_j = min(
bar_info["blocking"][0][j] - bar_info["blocking"][1][j],
bar_info["prp"][0][j] - bar_info["prp"][1][j])
if abs(top_j) >= abs(bot_j):
ym = min(top_j + yr * 0.05, ylim[1] - yr * 0.01)
else:
ym = max(bot_j - yr * 0.05, ylim[0] + yr * 0.01)
ax.text(x_pos[j], ym, "†",
ha="center", fontsize=11, color="black", fontweight="bold")
ax.set_xticks(x_pos)
ax.set_xticklabels([t.replace("\n", " ") for t in term_names_4], fontsize=8, rotation=15, ha="right")
ax.axhline(0, color="k", lw=0.5, ls="--")
ax.set_ylabel(clabel, fontsize=9)
if ax is axes.flat[0]:
ax.legend(fontsize=9)
level_str = f"wavg" if isinstance(LEVEL, str) else f"{LEVEL} hPa"
fig.text(0.5, -0.02,
r"Red *** = significantly $\neq$ 0 (p < 0.05) "
r"Black † = blocking vs PRP difference significant (p < 0.05)",
ha="center", fontsize=9, style="italic")
fig.suptitle(
f"4-Group Decomposition — Blocking vs PRP\n"
f"{STAGE} dh={DH} Level={level_str} Smooth={SMOOTH_DEG}° "
f"N_block={len(events['blocking'])} N_prp={len(events['prp'])} B={N_BOOT}",
fontsize=12, fontweight="bold", y=1.02)
fig.tight_layout()
plt.show()
7 Bar plot — 7-term comparison (Blocking vs PRP)
[8]:
fig, axes = plt.subplots(2, 2, figsize=(18, 18))
term_names_7 = list(ALL_TERMS_7.keys())
n_terms_7 = len(term_names_7)
x_pos_7 = np.arange(n_terms_7)
for ax, cname, clabel in zip(axes.flat, coef_names, coef_labels):
bar_info = {} # etype -> (means, lo_err, hi_err, sig_neq0)
# ── Draw bars ──
for i_et, etype in enumerate(["blocking", "prp"]):
means, los, his, sigs = [], [], [], []
for tname in term_names_7:
vals = boot_7[etype][tname][cname]
m = vals.mean()
lo, hi = np.nanpercentile(vals, [2.5, 97.5])
means.append(m)
los.append(m - lo)
his.append(hi - m)
sigs.append(lo * hi > 0)
bar_info[etype] = (means, los, his, sigs)
offset = -width / 2 if i_et == 0 else width / 2
ax.bar(x_pos_7 + offset, means, width,
yerr=[los, his], capsize=3,
label=etype.capitalize(),
color=colors_etype[etype], alpha=0.85)
# Red *** = CI excludes 0 (significantly != 0)
for j, (s, m) in enumerate(zip(sigs, means)):
if s:
y_off = his[j] if m >= 0 else -los[j]
ax.text(x_pos_7[j] + offset,
m + y_off * 1.1 + np.sign(m) * abs(m) * 0.05,
"***", ha="center", fontsize=7,
color="red", fontweight="bold")
# ── Between-group significance: blocking − PRP ──
ylim = ax.get_ylim()
yr = ylim[1] - ylim[0]
for j, tname in enumerate(term_names_7):
diff = boot_7["blocking"][tname][cname] - boot_7["prp"][tname][cname]
lo_d, hi_d = np.nanpercentile(diff, [2.5, 97.5])
if lo_d * hi_d > 0: # significant difference between groups
top_j = max(
bar_info["blocking"][0][j] + bar_info["blocking"][2][j],
bar_info["prp"][0][j] + bar_info["prp"][2][j])
bot_j = min(
bar_info["blocking"][0][j] - bar_info["blocking"][1][j],
bar_info["prp"][0][j] - bar_info["prp"][1][j])
if abs(top_j) >= abs(bot_j):
ym = min(top_j + yr * 0.05, ylim[1] - yr * 0.01)
else:
ym = max(bot_j - yr * 0.05, ylim[0] + yr * 0.01)
ax.text(x_pos_7[j], ym, "†",
ha="center", fontsize=11, color="black", fontweight="bold")
ax.set_xticks(x_pos_7)
ax.set_xticklabels(term_names_7, fontsize=7, rotation=20, ha="right")
ax.axhline(0, color="k", lw=0.5, ls="--")
ax.set_ylabel(clabel, fontsize=9)
ax.set_yscale("symlog")
if ax is axes.flat[0]:
ax.legend(fontsize=9)
fig.text(0.5, -0.02,
r"Red *** = significantly $\neq$ 0 (p < 0.05) "
r"Black † = blocking vs PRP difference significant (p < 0.05)",
ha="center", fontsize=9, style="italic")
fig.suptitle(
f"7-Term Decomposition — Blocking vs PRP\n"
f"{STAGE} dh={DH} Level={level_str} Smooth={SMOOTH_DEG}° "
f"N_block={len(events['blocking'])} N_prp={len(events['prp'])} B={N_BOOT}",
fontsize=12, fontweight="bold", y=1.02)
fig.tight_layout()
plt.show()
[10]:
# ── αx-only bar plot (7-term) — CI bars only, no significance markers ──
fig, ax = plt.subplots(figsize=(10, 5))
cname = "ax"
clabel = r"$\alpha_x$ (zonal propagation) [m/s]"
for i_et, etype in enumerate(["blocking", "prp"]):
means, los, his = [], [], []
for tname in term_names_7:
vals = boot_7[etype][tname][cname]
m = vals.mean()
lo, hi = np.nanpercentile(vals, [2.5, 97.5])
means.append(m)
los.append(m - lo)
his.append(hi - m)
offset = -width / 2 if i_et == 0 else width / 2
ax.bar(x_pos_7 + offset, means, width,
yerr=[los, his], capsize=4,
label=etype.capitalize(),
color=colors_etype[etype], alpha=0.85)
ax.set_xticks(x_pos_7)
ax.set_xticklabels(term_names_7, fontsize=9, rotation=20, ha="right")
ax.axhline(0, color="k", lw=0.5, ls="--")
ax.set_ylabel(clabel, fontsize=11)
ax.set_yscale("symlog")
ax.legend(fontsize=10)
ax.set_title(
f"$\\alpha_x$ (zonal propagation) — 7-Term Blocking vs PRP\n"
f"{STAGE} dh={DH} Level={level_str} Smooth={SMOOTH_DEG}° "
f"N_block={len(events['blocking'])} N_prp={len(events['prp'])} B={N_BOOT}",
fontsize=11, fontweight="bold")
fig.tight_layout()
plt.show()
[11]:
# ── αx from composite-mean pkl files (no bootstrap, no CI) ──
import pickle
PKL_PATHS = {
"blocking": "/net/flood/data2/users/x_yan/pvtend/outputs/blocking/composite_blocking.pkl",
"prp": "/net/flood/data2/users/x_yan/pvtend/outputs/prp/composite_prp.pkl",
}
# Load composites
composites = {}
for etype, path in PKL_PATHS.items():
with open(path, "rb") as f:
composites[etype] = pickle.load(f)
print(f"{etype}: N={composites[etype]._pick('original', STAGE, DH)[2]}")
# Helper: extract 2-D composite-mean field from pkl
def comp_field(comp, name, stage=STAGE, dh=DH, level=LEVEL):
key = name if name.endswith("_3d") else name + "_3d"
return comp.reduce_2d(key, stage, dh, level_mode=level)
# Build basis from composite means (same recipe as cell 6)
basis_comp = {}
for etype in PKL_PATHS:
c = composites[etype]
pv_n = comp_field(c, "pv_anom")
dx_n = comp_field(c, "pv_dx")
dy_n = comp_field(c, "pv_dy")
dx_dx = comp_field(c, "pv_dx_dx")
dy_dy = comp_field(c, "pv_dy_dy")
dx_dy = comp_field(c, "pv_dx_dy")
xr = c.x_rel[0, :] if c.x_rel.ndim == 2 else c.x_rel
yr = c.y_rel[:, 0] if c.y_rel.ndim == 2 else c.y_rel
basis_comp[etype] = compute_orthogonal_basis(
pv_n, dx_n, dy_n, xr, yr,
mask=MASK_SPEC,
apply_smoothing=True, smoothing_deg=SMOOTH_DEG, grid_spacing=GRID_SP,
pv_dx_dx_prev=dx_dx, pv_dy_dy_prev=dy_dy, pv_dx_dy_prev=dx_dy,
)
print(f"{etype} composite basis norms: "
+ " ".join(f"{k}={v:.3e}" for k, v in basis_comp[etype].norms.items()))
# 7-term composite fields and projection
COMP_GROUPS_7 = {
r"$-\bar{u}_\Psi\,q'_x$": lambda c: -comp_field(c, "u_rot_bar_pv_anom_dx"),
r"$-v'_\Psi\,\bar{q}_y$": lambda c: -comp_field(c, "v_rot_anom_pv_bar_dy"),
r"$-$Eddy rot": lambda c: -(
comp_field(c, "u_rot_anom_pv_anom_dx") + comp_field(c, "v_rot_anom_pv_anom_dy")),
r"$-$Eddy adiabatic div": lambda c: -(
comp_field(c, "u_div_adiabatic_pv_anom_dx") + comp_field(c, "v_div_adiabatic_pv_anom_dy")),
r"$-$Eddy diabatic div": lambda c: -(
comp_field(c, "u_div_diabatic_pv_anom_dx") + comp_field(c, "v_div_diabatic_pv_anom_dy")),
r"$-\omega\,\bar{q}_p$": lambda c: -(
comp_field(c, "w_adiabatic_pv_bar_dp") + comp_field(c, "w_diabatic_pv_bar_dp")),
r"$Q$": lambda c: comp_field(c, "Q"),
}
COMP_ALL_7 = {"dq/dt": lambda c: comp_field(c, "pv_dt")}
COMP_ALL_7.update(COMP_GROUPS_7)
# Project each composite-mean term
comp_ax = {} # comp_ax[etype][tname] = scalar αx
for etype in PKL_PATHS:
comp_ax[etype] = {}
c = composites[etype]
b = basis_comp[etype]
for tname, tfunc in COMP_ALL_7.items():
fld = smooth(tfunc(c))
p = project_field(fld, b)
comp_ax[etype][tname] = p["ax"]
# ── Plot ──
fig, ax = plt.subplots(figsize=(10, 5))
comp_term_names = list(COMP_ALL_7.keys())
n_t = len(comp_term_names)
x_p = np.arange(n_t)
for i_et, etype in enumerate(["blocking", "prp"]):
vals = [comp_ax[etype][t] for t in comp_term_names]
offset = -width / 2 if i_et == 0 else width / 2
ax.bar(x_p + offset, vals, width,
label=etype.capitalize(),
color=colors_etype[etype], alpha=0.85)
ax.set_xticks(x_p)
ax.set_xticklabels(comp_term_names, fontsize=9, rotation=20, ha="right")
ax.axhline(0, color="k", lw=0.5, ls="--")
ax.set_ylabel(r"$\alpha_x$ (zonal propagation) [m/s]", fontsize=11)
ax.set_yscale("symlog")
ax.legend(fontsize=10)
N_blk = composites["blocking"]._pick("original", STAGE, DH)[2]
N_prp = composites["prp"]._pick("original", STAGE, DH)[2]
ax.set_title(
f"$\\alpha_x$ — Composite-mean projection (no bootstrap)\n"
f"{STAGE} dh={DH} Level={level_str} Smooth={SMOOTH_DEG}° "
f"N_block={N_blk} N_prp={N_prp}",
fontsize=11, fontweight="bold")
fig.tight_layout()
plt.show()
blocking: N=1260
prp: N=2097
blocking composite basis norms: beta=2.831e+02 ax=2.733e+02 ay=1.639e+02 gamma1=9.718e+01 gamma2=6.578e+01 sigma=3.624e+01
prp composite basis norms: beta=1.998e+02 ax=1.678e+02 ay=9.868e+01 gamma1=4.714e+01 gamma2=4.432e+01 sigma=1.418e+01
/tmp/ipykernel_1373389/659284345.py:34: UserWarning: compute_orthogonal_basis: grid_spacing=1.5°, center_lat=60.0°N → dx(center)=83.4 km, dy=166.8 km
basis_comp[etype] = compute_orthogonal_basis(
