Parked. Not installed with odsl. See PARKED.md. The text below is the last active description of the module.

Pure-numpy quasi-geostrophic (QG) diagnostics computed from sea-surface height anomaly (SSHA). Given a 2-D SSHA field and a handful of explicit physical parameters, qg reconstructs the surface geostrophic flow, continues it into the interior with a Surface Quasi-Geostrophic (SQG) inversion, and solves the QG omega equation for vertical velocity, plus the usual strain / vorticity / Okubo-Weiss deformation diagnostics.

Overview

The physics, in one page:

  • Surface streamfunction. Geostrophy links SSHA η to a streamfunction ψ_s = (g/f₀) η; horizontal velocity follows from u = −∂ψ/∂y, v = ∂ψ/∂x.
  • SQG inversion. Assuming zero interior potential vorticity and constant stratification N, the streamfunction is downward-continued in spectral space, ψ̂(k, z) = ψ̂_s(k) · exp(−N|k|z / f₀), giving a full 3-D ψ(z, y, x) from surface data alone. The field is detrended and Tukey-tapered before the FFT.
  • Omega equation. With the 3-D ψ we form the Q-vector (Hoskins et al. 1978), take its divergence, and solve N²∇²w + f₀²∂²w/∂z² = 2∇·Q for the QG vertical velocity w, with rigid-lid w = 0 at the top and bottom.
  • Deformation diagnostics. Strain rate, relative vorticity, velocity-gradient magnitude, and the Okubo-Weiss parameter OW = S² − ζ² classify each point as eddy-dominated (OW < 0) or strain/front-dominated (OW > 0).

Design: every function takes plain numpy arrays and explicit physical parameters (f0, N/N², dx, dy, zlev). There is no I/O, no plotting, and no global state, so qg is completely data-source agnostic — feed it SWOT, multi-altimeter OI maps, model output, or an analytic field. Tests under qg/tests/ are analytic-solution checks (pytest); they are not part of the default odsl test run while this package is parked.

Public API

Exported from qg/__init__.py (__all__); implemented in qg/core.py. All 2-D arrays are (ny, nx), all 3-D arrays are (nz, ny, nx). Grid spacings dx, dy, dz are in metres; angular wavenumbers are used internally.

Function Signature (summary) Returns / units
coriolis coriolis(lat) Coriolis parameter f = 2Ω sin(lat) [s⁻¹]
ssha_to_streamfunction ssha_to_streamfunction(eta, f0) surface ψ_s = (g/f₀)η [m²/s]
geostrophic_velocity geostrophic_velocity(psi, dx, dy) (u, v) [m/s]; accepts 2-D or 3-D psi
relative_vorticity relative_vorticity(u, v, dx, dy) ζ = ∂v/∂x − ∂u/∂y [s⁻¹]
velocity_gradient_magnitude velocity_gradient_magnitude(u, v, dx, dy) |∇V| [s⁻¹]
strain_rate strain_rate(u, v, dx, dy) normal strain rate S [s⁻¹]
okubo_weiss okubo_weiss(u, v, dx, dy, f0=None) OW = S² − ζ² [s⁻²]; dimensionless (/f₀²) if f0 given
sqg_inversion sqg_inversion(eta, f0, N, dx, dy, zlev, taper_frac=0.25) (psi_3d, window) — psi_3d (nz, ny, nx) [m²/s], window (ny, nx)
qvector qvector(u_g, v_g, b, dx, dy) (Q1, Q2) [m/s³]
qvector_divergence qvector_divergence(Q1, Q2, dx, dy) ∇·Q [m/s⁴]
omega_equation omega_equation(psi_3d, zlev, f0, N2, dx, dy, taper_frac=0.25) dict (see below)
apply_taper apply_taper(field, frac=0.25) (tapered, window) — 2-D Tukey window
detrend_2d detrend_2d(field) field minus mean + linear plane (NaN-aware)

omega_equation returns a dict with keys w (vertical velocity, m/day), u_g, v_g [m/s], b (buoyancy b = f₀ ∂ψ/∂z, [m/s²]), zeta [s⁻¹], and div_Q [m/s⁴], all (nz, ny, nx). The spectral solve keeps the full complex Q-divergence forcing (dropping the imaginary part biases w for a general field).

Constants exported: OMEGA (Earth rotation rate), G (gravity), R_EARTH, DEG_KM.

Installation

qg needs only the core scientific stack that ships with odsl — numpy and scipy (scipy.ndimage.gaussian_filter). No extras are required.

python -m pip install -e .          # from the repo root

It imports two ways, import qg and import odsl_code.qg.

Usage

import numpy as np
from qg import (
    coriolis, ssha_to_streamfunction, geostrophic_velocity,
    relative_vorticity, strain_rate, okubo_weiss,
    sqg_inversion, omega_equation,
)

# --- Grid and physical parameters ---
lat = 34.0
f0 = coriolis(lat)          # ~8.1e-5 s^-1
N  = 5e-3                    # Brunt-Vaisala frequency [rad/s]
N2 = N**2
dx = dy = 4e3               # 4 km grid spacing [m]

nx = ny = 150
x = (np.arange(nx) - nx / 2) * dx
y = (np.arange(ny) - ny / 2) * dy
X, Y = np.meshgrid(x, y)

# --- Build an SSHA field (Gaussian eddy, 5 cm) ---
eta = 0.05 * np.exp(-X**2 / (60e3)**2 - Y**2 / (40e3)**2)   # [m], (ny, nx)

# --- Surface diagnostics ---
psi_s = ssha_to_streamfunction(eta, f0)          # (ny, nx) [m^2/s]
u, v  = geostrophic_velocity(psi_s, dx, dy)      # (ny, nx) [m/s]
zeta  = relative_vorticity(u, v, dx, dy)         # (ny, nx) [s^-1]
S     = strain_rate(u, v, dx, dy)                # (ny, nx) [s^-1]
ow    = okubo_weiss(u, v, dx, dy, f0=f0)         # (ny, nx) dimensionless (OW/f0^2)

# --- SQG inversion: 3-D streamfunction ---
zlev = np.linspace(0, 1000, 40)                  # depth levels [m], positive down
psi_3d, window = sqg_inversion(eta, f0, N, dx, dy, zlev, taper_frac=0.15)
print(psi_3d.shape)                              # (40, 150, 150)

# --- Omega equation: vertical velocity ---
res = omega_equation(psi_3d, zlev, f0, N2, dx, dy, taper_frac=0.15)
w = res["w"]                                     # (40, 150, 150) [m/day]
print(w.shape, w[0].max(), w[-1].max())          # ~0 at top/bottom (BC)

Expected shapes/units: surface fields are (ny, nx); sqg_inversion and omega_equation return (nz, ny, nx). OW/f₀² is O(1); w is in m/day and is identically ~0 at the first and last zlev (rigid-lid boundary condition).

Physical conventions

  • f0 — Coriolis parameter [s⁻¹] evaluated at the domain centre (f-plane). Use coriolis(lat); it is negative in the Southern Hemisphere.
  • Stratification — sqg_inversion takes N [rad/s]; omega_equation takes N2 = N² [s⁻²]. Constant N² (uniform stratification) is assumed.
  • Grid spacing — dx (x / cross-track / zonal) and dy (y / along-track / meridional) in metres; dz is taken from zlev[1] - zlev[0] (uniform levels expected). Derivatives use numpy.gradient (centred, one-sided at edges).
  • Depth zlev — positive downward, starting at 0 (surface). SQG decay uses exp(−N|k|z/f₀).
  • FFT / wavenumbers — real 2-D FFT (np.fft.rfft2); wavenumbers from rfftfreq/fftfreq are converted to angular wavenumbers (× 2π, rad/m). The k = 0 mode is zero after detrending.
  • Boundary conditions — the omega solve imposes w = 0 at the top (zlev[0]) and bottom (zlev[-1]) via a per-wavenumber tridiagonal solve.
  • Sign conventions — u = −∂ψ/∂y, v = ∂ψ/∂x; ζ = ∂v/∂x − ∂u/∂y; OW = S² − ζ² (negative → rotation-dominated, positive → strain-dominated); buoyancy b = f₀ ∂ψ/∂z.
  • Edge treatment — SSHA is detrended (detrend_2d) and Tukey-tapered (apply_taper, taper_frac) before spectral operations to suppress ringing. Trust the interior; gradient/taper artifacts contaminate the outer rows/cols.

Tests

pytest qg/tests

The suite validates the solvers against closed-form analytic solutions:

  • test_sqg_analytic.py — SQG inversion vs. the analytic ψ(x,y,d) = (B₀/NK) cos(kx)cos(ly) exp(−NKd/f₀); checks ψ, (u, v), the exponential decay slope −NK/f₀, and the e-folding depth H = f₀/(NK).
  • test_omega_analytic.py — omega solver vs. ω = −F₀ / [N²(k²+l²) + f₀²(mπ/H)²] · cos(kx)cos(ly)sin(mπz/H); checks amplitude, spectral purity, and the w = 0 rigid-lid BCs.
  • test_elliptical_eddy.py — an anisotropic Gaussian eddy exercising the full chain (velocity, vorticity, strain, Okubo-Weiss, velocity-gradient magnitude, SQG inversion, omega equation) against hand-derived closed forms.

Each test file is also runnable directly (python qg/tests/test_*.py) and writes a diagnostic PNG next to itself.

  • ../README.md — repository-wide layout and conventions (qg functions are intentionally data-source agnostic: plain numpy arrays with explicit f, N², and grid spacing).
  • multi_altimeters — a producer of gridded SSHA maps that pair naturally with these diagnostics.