pvtend 2.10.7

PV tendency decomposition diagnostics for blocking and weather extremes

pvtend

PV tendency decomposition for atmospheric blocking, propagating anticyclones, and all synoptic-scale cyclonic event lifecycle analysis.

pvtend diagnoses the growth, propagation, and decay of mid-latitude weather events by decomposing potential vorticity (PV) tendencies from ERA5 pressure-level data onto physically meaningful components using an orthogonal basis framework. This is the Part I work of Yan et al. (in prep.) about blocking lifecycle analyses on onset, peak, decay stages.

Gallery

Idealized validation — a Gaussian PV anomaly with prescribed propagation, intensification, and deformation is decomposed into six orthogonal bases and reconstructed with near-zero residual.	Real ERA5 blocking event (track 425) — animated lifecycle showing total PV on a cartopy map (left) and the six projected basis components (right) evolving from 13 h pre-onset to 12 h post-onset. The analysis is done on a weighted average surface across 300, 250, 200 hPa levels.
Geopotential-height (Z500) variant of the six-basis decomposition (track 425) — animated lifecycle showing Z anomaly from the 1990–2020 hourly climatology, with adaptive prenorm and blockid contour overlay. See notebook 03z_six_basis_projection_geopotential.

Event catalogues

Blocking and PRP-high events are identified as persistent anticyclonic anomalies in 500 hPa geopotential height. We are using TempestExtremes v2.1 to track contiguous Z500 anomaly features that exceed a fixed threshold for ≥5 days, producing CSV catalogues with columns for event ID, centre lat/lon, onset/peak/decay timestamps, and area. Following the threshold as in Drouard et al. (2021), we separate the tracked features into blocking and propagating (prp) high pressure systems.

Sample catalogue (ERA5, 1990–2020 blocking): ERA5_TempestExtremes_z500_anticyclone_blocking.csv

The CSVs are the inputs for pvtend-pipeline compute, which extracts event-centred patches and runs the full PV-tendency decomposition for each event in the blocking/prp catalogue.

Features

• Helmholtz-first architecture (v2.0): Helmholtz decomposition on the total wind field; climatological Helmholtz pre-computed as 24 monthly NetCDF files; anomaly Helmholtz = total − clim (no separate anomaly solve)
• 53 cross-term PV tendency budget: 20 primary + 16 alt-vertical + 16 div dry/moist horizontal + Q_LHR, all written per-timestep to NPZ
• PV tendency computation: RHS has zonal advection, baroclinic counter propagation, vertical advection, and approximated diabatic heating terms.
• QG omega solver: Hoskins Q-vector formulation with two methods: LOG20/SIP (default, Numba-accelerated 3-D elliptic, Li & O'Gorman 2020) and SP19 (Steinfeld & Pfahl 2019 empirical 1/3 scaling). Optional center_lat for dynamic f₀.
• Helmholtz decomposition: Spherical vorticity/divergence (with tan φ/a metric), conservative spherical Poisson solver (FFT in lon + tridiagonal in lat), spectral gradient for wind recovery — all on the full NH grid. Since v2.10 a spherical-harmonic backend (pvtend.sh_ops, via pyspharm) is the default (method="spectral"); legacy FFT path remains as method="fft". NH inputs are parity-mirrored to the global sphere, giving a near-zero harmonic residual.
• Four-way omega decomposition: ω_dry (QG A+B), ω_qg_moist (term C via ∂T/∂t), ω_emanuel_moist (Emanuel LHR), ω_moist (full residual), with corresponding divergent winds recovered by independent spherical Poisson inversion (verified linear to machine precision)
• Orthogonal basis decomposition: Projects PV tendency onto intensification (β), propagation (αx, αy), and deformation (γ) modes. Built-in temporal down-scaling (bi-linear interpolation, α = 0.75 by default) from hourly to 15-minute evaluation instants via _next keyword arguments. Single-blob selection: when the threshold mask produces multiple disconnected regions, only the connected component enclosing (or nearest to) the patch centre is retained
• RWB detection: Two classification methods — bay (path-order, recommended with circumpolar-cropped contours) and tilt (centerline slope ±0.15 dead zone). Circumpolar-first contour extraction for robust NH analysis.
• Composite lifecycle: Multi-stage ensemble averaging with onset/peak/decay staging
• NaN-safe throughout: All grid, derivative, solver, bootstrap, and plotting routines use nanmean/nanpercentile to handle partial-NaN edge events without corrupting composites or flipping projection signs
• ω blowup detection (v2.10.1): Post-hoc CSV scan flagging ERA5 timestamps whose |ω| at 300 hPa exceeds a fixed Pa s⁻¹ hard cutoff (default 10 Pa/s — synoptic 300 hPa ω rarely exceeds 2–3 Pa/s, so |ω|>10 Pa/s in the derived ω_dry/ω_moist/ω_LHR is essentially always solver blowup; configurable via --threshold). Replaces the silent in-solver clip used in pvtend ≤ 2.9, which was corrupting the affected events.
• CLI pipeline: End-to-end processing via pvtend-pipeline command (compute, clim-helmholtz, classify, composite, decompose, blowup-scan)

Installation

# From PyPI
pip install pvtend

# Or with uv (fast, Rust-based installer)
uv pip install pvtend

# From source (development)
git clone https://github.com/yanxingjianken/pvtend.git
cd pvtend
pip install -e ".[dev]"

# With micromamba
micromamba create -f environment.yml
micromamba activate pvtend_env
pip install -e ".[dev]"

Spherical-harmonic backend (recommended)

The Helmholtz / QG-ω velocity-potential / mixed-second-derivative operators default to a spherical-harmonic backend (pvtend.sh_ops, built on pyspharm / windspharm). These two packages compile Fortran code at install time, so they are kept as an optional extra:

# conda-forge ships pre-built wheels (recommended)
micromamba install -c conda-forge windspharm pyspharm

# Or via pip with numpy already installed and a Fortran compiler:
pip install --no-build-isolation pyspharm windspharm

With these installed, method="spectral" (the default) is used. Without them, callers fall back to the legacy FFT path or raise a clear error directing the user to install the SH extras.

Quick Start

import numpy as np
from pvtend import NHGrid, ddx, ddy, compute_orthogonal_basis, project_field

# Grid setup
grid = NHGrid(lat=np.linspace(90, 0, 61), lon=np.linspace(-180, 178.5, 240))
dx_arr = grid.dx_arr  # zonal spacing per latitude [m]
dy = grid.dy           # meridional spacing [m]

# Compute zonal derivative
dfdx = ddx(field, dx_arr, periodic=False)

# Orthogonal basis decomposition
basis = compute_orthogonal_basis(pv_anom, pv_dx, pv_dy, x_rel, y_rel)
result = project_field(tendency, basis)
print(f"β = {result['beta']:.3e}")  # intensification rate

CLI Pipeline

# Step 0a: Pre-compute Helmholtz climatology (once, ~5 min)
pvtend-pipeline clim-helmholtz \
    --clim-dir /data/climatology/ \
    --out-dir /data/climatology/

# Step 1: Compute PV tendencies → per-event NPZ files
pvtend-pipeline compute \
    --event-type blocking \
    --events-csv events.csv \
    --era5-dir /data/era5/ \
    --clim-path /data/climatology/era5_hourly_clim.nc \
    --clim-helmholtz-dir /data/climatology/ \
    --out-dir /data/composite_blocking_tempest/ \
    --dh-range=-49:25 --skip-existing

# Step 2: RWB classification → variant tracksets PKL
#   --levels accepts integer hPa values or 'wavg' (weighted-average Z)
pvtend-pipeline classify \
    --npz-dir /data/composite_blocking_tempest/ \
    --output /data/outputs/rwb_variant_tracksets.pkl \
    --stages onset peak decay \
    --levels 500 400 300 200 --threshold 3

# Step 3: Variant-aware composite accumulation → composite PKL
pvtend-pipeline composite \
    --npz-dir /data/composite_blocking_tempest/ \
    --rwb-pkl /data/outputs/rwb_variant_tracksets.pkl \
    --pkl-out /data/outputs/composite.pkl

# Step 4 (optional): Orthogonal-basis decomposition
pvtend-pipeline decompose \
    --pkl-in /data/outputs/composite.pkl \
    --out-dir /data/outputs/decomp/

Workflow

graph TD
    A[ERA5 Monthly NetCDF] --> B[pvtend.preprocessing]
    B --> C[Regridded NH Grid]
    C --> D[pvtend.climatology]
    D --> E[Monthly Climatology]
    E --> E2[pvtend.climatology — clim-helmholtz]
    E2 --> E3[24 Helmholtz Clim NetCDF]
    C & E & E3 --> F[pvtend.tendency.TendencyComputer]
    F --> G[PV Tendency Terms]
    F --> I[pvtend.helmholtz — Helmholtz on total wind]
    I --> I2[u_rot, u_div, v_rot, v_div]
    I2 & E3 --> I3[Anomaly Helmholtz = total − clim]
    G --> H[pvtend.omega — QG ω solver]
    H & I2 --> J[pvtend.moist_dry — ω splitting + 4 independent Poisson inversions]
    G & J & I3 --> K[Per-event NPZ patches — 53 cross-terms]
    K --> L1[pvtend.classify — RWB Pass 1]
    L1 --> L1a[rwb_variant_tracksets.pkl]
    K & L1a --> L2[pvtend.composite_builder — Pass 2]
    L2 --> L2a[composite.pkl]
    L2a --> W1[website/export_composite_web.py]
    L2a --> W2[website/build_and_export_clusters.py]
    W1 --> W3[blocking_export and prp_export JSON.gz]
    W2 --> W3
    W3 --> W4[blocking-plots viewer with contour and arrow overlays]
    L2a --> M[pvtend.decomposition — Orthogonal basis]
    M --> N[β, αx, αy, γ coefficients]
    N --> O[pvtend.plotting — Publication figures]

Website Export

The repository also contains a lightweight website export pipeline for the blocking/PRP composite viewer. The JSON exports under website/blocking_export/ and website/prp_export/ are built from the saved composite pickles and include grouped field metadata, basis payloads, and vector-overlay definitions so the viewer can switch between scalar contours and wind-arrow overlays.

graph TD
    A[composite_blocking.pkl or composite_prp.pkl] --> B[website/export_composite_web.py]
    C[composite_blocking_clusters.pkl or composite_prp_clusters.pkl] --> D[website/build_and_export_clusters.py]
    B --> E[Parent-category JSON.gz exports]
    D --> F[Cluster and subcluster JSON.gz exports]
    E --> G[Hugging Face dataset blocking-composites]
    F --> G
    G --> H[blocking-plots viewer]

Publication Figures

The paper/generals_paper_1/scripts/ directory contains the current paper Figure 5-8 generation workflow. Shared budget definitions live in _paper_budget.py, which centralizes the Eq. 11 term mapping, the pressure-weighted wavg level, the q'<0 projection basis, p99.9 track-exclusion lists, and the Figure 7 bootstrap cache location.

• fig5_beta_closure_qg_onset_wavg.py renders the blocking-onset raw PV-tendency closure and corresponding beta*Phi_1 projections in two four-column blocks. The stationary-flow and baroclinic advection terms are split into separate columns, with term subtitles, panel-specific colorbars, a shared range for the six process-projection colorbars, total-PV contours on the LHS panel, and wind arrows on the other raw-term panels.
• fig6_beta_stacked_qg_wavg.py renders the wavg-only onset lifecycle stacked-bar budget for beta, with the direct partial q / partial t projection overlaid as a black step curve and colors matched to examples/05_stacked_bar_beta.ipynb cell 14.
• fig7_ax_bootstrap_blocking_vs_prp_qg_wavg.py renders peak-stage blocking versus propagating-anticyclone alpha_x bootstrap bars from event NPZs after excluding p99.9 blowup tracks. Projections use the full event fields without a pre-projection significance mask, the pv_dt bars include simple mean-value labels, and repeated runs reuse paper/generals_paper_1/scripts/cache/ when the inventory and settings match.
• fig8_gamma_rwb_awb_cwb_onset_wavg.py renders AWB/CWB onset raw maps and -gamma_1*Phi_4 - gamma_2*Phi_5 deformation projections for four selected Eq. 11 terms. It uses QG-diabatic paper definitions, one colorbar per panel, and Figure-2-style stretch/compression arrows.

All four scripts write PNG and PDF files into paper/generals_paper_1/figures/, and paper/generals_paper_1/main.tex includes the generated PNGs.

graph TD
    A[outputs/blocking/composite_blocking.pkl] --> B[_paper_budget.py]
    B --> C[fig5_beta_closure_qg_onset_wavg.py]
    B --> D[fig6_beta_stacked_qg_wavg.py]
    B --> M[fig8_gamma_rwb_awb_cwb_onset_wavg.py]
    E[outputs/blocking and outputs/prp event NPZs] --> B
    F[p99.9 exclude CSVs] --> B
    B --> G[fig7_ax_bootstrap_blocking_vs_prp_qg_wavg.py]
    G --> H[scripts/cache/bootstrap NPZ]
    C --> I[figures/fig5_beta_closure_qg_onset_wavg PNG/PDF]
    D --> J[figures/fig6_beta_stacked_qg_wavg PNG/PDF]
    G --> K[figures/fig7_ax_bootstrap_blocking_vs_prp_qg_wavg PNG/PDF]
    M --> N[figures/fig8_gamma_rwb_awb_cwb_onset_wavg PNG/PDF]
    I & J & K & N --> L[paper/generals_paper_1/main.tex]

Package Structure

src/pvtend/
├── __init__.py          # Public API
├── _version.py          # Version
├── cli.py               # CLI entry point (clim-helmholtz, compute, classify, composite, decompose)
├── constants.py         # Physical constants
├── grid.py              # NH grid & event patches
├── preprocessing.py     # ERA5 loading & regridding
├── derivatives.py       # Finite difference operators
├── climatology.py       # Fourier-filtered climatology + Helmholtz climatology (compute/load)
├── omega.py             # QG omega solver (LOG20/SIP or SP19)
├── helmholtz.py         # Helmholtz decomposition (spherical Poisson + spectral gradient + laplacian_spherical_fft)
├── moist_dry.py         # Moist/dry omega split & independent Poisson wind recovery (solve_chi_from_omega, verify_div_additivity)
├── isentropic.py        # Isentropic PV-tendency diagnostics
├── tendency.py          # Main pipeline: Helmholtz-first, 53 cross-terms, data loading, derivatives, NPZ output
├── classify.py          # RWB classification Pass 1 (AWB/CWB/NEUTRAL → variant PKL)
├── composite_builder.py # Variant-aware composite accumulation Pass 2
├── rwb.py               # RWB detection (bay & tilt methods, circumpolar-first)
├── composites.py        # Legacy composite lifecycle
├── data/                # Bundled sample data
│   ├── __init__.py      # load_idealized_pv() loader
│   └── idealized_pv.npz # Synthetic Gaussian PV evolution
├── decomposition/       # Orthogonal six-basis framework
│   ├── __init__.py
│   ├── smoothing.py
│   ├── basis.py
│   ├── interpolation.py # Temporal bi-linear interpolation (lerp_fields)
│   └── projection.py
├── plotting/            # Visualization
│   ├── __init__.py
│   ├── basis_plots.py
│   ├── coefficient_plots.py
│   ├── field_plots.py
│   ├── composite_explorer.py  # plot_var: single-variable composite explorer with bootstrap
│   └── baroclinic.py          # plot_baroclinic_tilt: two-level v′ overlay
└── io/                  # File I/O
    ├── __init__.py
    ├── era5.py
    ├── npz.py
    └── pkl.py

website/
├── export_composite_web.py   # JSON.gz export for the blocking/PRP web viewer
├── build_and_export_clusters.py # Cluster-aware export and diagnostics generation
├── upload_to_hf.py           # Sync local website exports to Hugging Face
├── blocking_export/          # Exported blocking viewer payloads
├── prp_export/               # Exported PRP viewer payloads
└── cluster_diagnostics/      # Cluster composite diagnostics used for QA

Example Notebooks

Notebooks using real ERA5 blocking event data from the composite_blocking_tempest, composite_prp_tempest (single event npz), tempest_extreme_4_basis/outputs, and tempest_extreme_4_basis/outputs_prp (composite pkl) pipeline:

Notebook	Description
00_idealized_6basis	Idealized Gaussian PV anomaly: prescribed β/αx/αy/γ₁/γ₂/σ at two timesteps, 6-basis visualisation, Gram-Schmidt, projection & reconstruction
01_rwb_and_derivatives	Grid setup, ddx/ddy/ddp derivatives, RWB detection on a real event
02_helmholtz_and_qg_omega	3-D Helmholtz decomposition, QG omega (LOG20 vs SP19), moist/dry ω split
03_six_basis_projection	Orthogonal 6-basis (Φ₁–Φ₆), project dq'/dt → β/αx/αy/γ₁/γ₂/σ, lifecycle curves
03_six_basis_projection_composite	↳ Composite variant: 6-basis projection averaged across multiple events
03c_six_basis_cyclone	↳ Cyclone variant: 6-basis projection for a 300 hPa cyclone (PV > 0), lifecycle + budget closure
03prp_six_basis_anticyclone_timed_bases	↳ Anticyclone variant: 6-basis projection for a 300 hPa anticyclone (PV < 0, mask="< -2e-7"), current-time basis
03z_six_basis_projection_geopotential	↳ Supplement: same 6-basis projection using geopotential height Z instead of PV
04_single_var_composite	Single-variable composite explorer on pressure levels using pvtend.plotting.plot_var
04i_single_var_isentropic_composite	↳ Supplement: same as 04 but on isentropic (θ) surfaces
05_stacked_bar_beta	Stacked-bar β decomposition by PV-tendency term across lifecycle hours
05b_grouped_terms_bootstrap	↳ Supplement: grouped PV-tendency terms with bootstrap resampling & significance
06_baroclinic_structure	3-D composite PV anomaly, lon–p cross-sections, 2-PVU tropopause, v′ tilt via plot_baroclinic_tilt
07_facet_blocking_vs_prp	Facet comparison of blocking vs PRP: bar charts with bootstrap significance, shared-cbar spatial maps, baroclinic tilt

Testing

pytest tests/ -v

Documentation

Full documentation at pvtend.readthedocs.io.

Build locally:

cd docs && make html

Citation

If you use this package in your research, please cite:

@software{yan2025pvtend,
  author = {Yan, Xingjian and Tamarin-Brodsky, Talia},
  title = {pvtend: PV tendency decomposition for atmospheric blocking},
  year = {2025},
  url = {https://github.com/yanxingjianken/pvtend}
}

License

MIT — see LICENSE.