multi_altimeters — multi-altimeter OI & point-vortex SSH

Multi-altimeter mapping and inversion for the ODSL pipeline. This subpackage turns local SWOT L2 LR SSH Expert swaths and conventional RADS nadir altimetry into gridded sea-surface-height-anomaly (SSHA) products, and provides an experimental point-vortex inversion of SSH from a SWOT segment.

Overview

Four workflows live here:

• Optimal-interpolation (OI) mapping — map_multi_altimeters.py searches and loads SWOT Expert files under /spray/swot/versionD (via swot/utils.py) and RADS files under /spray/rads/data (via rads/utils.py), filters by time and bounding box, reduces both to a common observation table (time, lat, lon, ssha, uncertainty, plus source, mission, cycle, pass, path metadata), and grids SSHA with local ordinary kriging (covariance-based optimal interpolation). Each grid point is solved from its space-time nearest neighbors with a selectable spatial × temporal decorrelation model (Gaussian by default; exponential and second-order-autoregressive available) and uncertainty^2 on the observation-covariance diagonal — per-mission observation errors supported. Output is a NetCDF map with the analysis variable ssh_anomaly_oi.

• Regional space-time cube — map_region_ssha.py uses the same loaders and OI/tide numerics, but writes a lon/lat time cube: a 2nd-order spatial background and harmonic point-source internal tide are fit once, then the residual is mapped at every analysis time (default 6-hourly) onto a regular grid whose posting is dx_km (default 2 km) at the bbox-center latitude.

• Expanded swath-grid mapping — grid_swath_ssha.py maps one SWOT Expert segment plus nearby nadir altimetry on a grid that follows the SWOT swath but is widened by --expansion-km (default 100 km) per side at the same 2 km posting; nadir data are searched inside the expanded footprint within ±--delta-t-days (default 4 days) of the segment time. Same OI core, swath targets.

• Point-vortex SSH inversion — point_vortex_ssh_iteration0.py fits a field of point vortices on a regular lon/lat grid so their geostrophic streamfunction (eta = f0/g * psi) reproduces observed SSHA. It is a first-step, spatial-only weighted-least-squares/coordinate-descent solve; time is used only to select the input window. Reusable numerics live in point_vortex_core.py.

Shared orchestration, longitude handling, SWOT thinning, and the OI solver live in utils.py; background/internal-tide fitting in internal_tides.py. The mathematics of the swath-grid pipeline (grid geometry, thinning operators, background polynomial, point-source tide model, kriging equations) is written out in methodology.md.

Installation & dependencies

The subpackage ships with the odsl distribution. Core deps (numpy, pandas, xarray, netCDF4, scipy, pyproj, …) are always installed. Feature deps are extras:

python -m pip install -e ".[plot,kriging]"   # or ".[all]"

• plot extra → cartopy (map backgrounds; all plotting scripts degrade gracefully to plain lon/lat axes when cartopy is missing).
• kriging extra → pykrige (declared for kriging workflows; note the OI solver in utils.py is a self-contained numpy ordinary-kriging system and does not require pykrige at runtime).

Scripts & entry points

Script	Entry point	Purpose
OI mapping
map_multi_altimeters.py	odsl-multi-altimeter-oi	Main OI driver: load SWOT+RADS, grid SSHA with ordinary kriging, write NetCDF.
map_region_ssha.py	odsl-multi-altimeter-region	Regional lon/lat time cube: fit background + harmonic internal tide once, OI the residual at every dt (default 6 h) onto a dx_km (default 2 km) grid.
grid_swath_ssha.py	odsl-multi-altimeter-swath-grid	Grid one SWOT Expert segment + nearby nadir altimetry on an expanded swath-following grid (--expansion-km per side, or --target-width-km total width) and OI-map it.
grid_swath_ssha_pass.py	odsl-swot-swath-grid-300km	Pass-by-pass production driver: 300 km total width, 5×5 SWOT thinning with sample-std uncertainty, 0.02 m nadir error, 4-day time half-width, ±5.5-day search, other SWOT passes included.
plot_oi_result.py	odsl-multi-altimeter-plot	Plot/save an OI NetCDF result (default variable ssh_anomaly_oi).
Regional subsetting
(moved)	odsl-super-subsetter	The regional native-structure subsetter (eleven /spray archives → one HDF5) now lives in its own package — see ../super_subsetter.
Plotting / demos
demo_swath_grid_ssha.py	(run directly)	Run the expanded swath-grid pipeline for a California-coast SWOT segment and plot three figures: swath-coordinate panels (SWOT / nadir / OI), a lon/lat map with nadir ground tracks, and a cross-track profile.
demo_plot_thinned_swot.py	(run directly)	Inspect SWOT N×N median thinning: two-panel thinned SSHA + per-block uncertainty; uses common.plotting.add_map_context.
example_oi_gulf_america_20250110.py	(run directly)	Configured Gulf-of-America OI example (13×13 SWOT thinning, 70 km / 10 day half-widths) with an obs-overlay figure.
Point-vortex inversion
point_vortex_ssh_iteration0.py	(run directly)	Main inversion driver. Spatial-only point-vortex fit; writes reconstructed SSH/streamfunction/velocity/vorticity NetCDFs, obs-fit CSV, and a PNG.
point_vortex_ssh_from_swot_segment.py	(run directly)	Derive domain + target time from one SWOT Expert segment, then run the iteration-0 solver from nearby RADS; write a comparison CSV at the SWOT pixels.
experimental/point_vortex_ssh_iteration1.py	(run directly)	Quarantined / unfinished moving-center solve. The old path point_vortex_ssh_iteration1.py is a stub that exits with a pointer here.
point_vortex_core.py	(library)	Reusable VortexGrid, build_vortex_grid, load_observations, solve_vorticity*, kernel/evaluation, dataset validation.

Usage

Grid SSHA for a Gulf of America date, then plot it:

# OI map: SWOT + RADS, 3x3 SWOT median thinning, 0.1-degree grid
odsl-multi-altimeter-oi \
    --start 2025-01-01 --end 2025-01-10 \
    --bbox -98 18 -80 31 \
    --sources swot rads \
    --swot-median-window 3 \
    --grid-spacing 0.1 \
    --spatial-half-width-km 30 \
    --time-half-width-days 3 \
    --target-time 2025-01-05 \
    --out tmp/gom_multi_altimeter_oi.nc \
    --observations-out tmp/gom_multi_altimeter_obs.csv

# Plot the resulting field
odsl-multi-altimeter-plot tmp/gom_multi_altimeter_oi.nc \
    --var ssh_anomaly_oi \
    --out tmp/gom_multi_altimeter_oi.png

Both scripts also run directly, e.g. python multi_altimeters/map_multi_altimeters.py ….

Quick smoke test (cap files and points):

python multi_altimeters/map_multi_altimeters.py \
    --start 2025-01-01 --end 2025-01-01T12:00:00 \
    --bbox -98 18 -80 31 \
    --max-files 10 --max-points 20000 \
    --out tmp/test_oi.nc

--bbox also accepts 0..360 longitudes (--bbox 262 18 280 31). Instead of --start/--end, target a single SWOT granule with --swot-file <name> --delta-t-days 2.

Useful OI flags (defaults in parentheses): --grid-spacing (0.1), --spatial-half-width-km (30), --time-half-width-days (3), --covariance-model {gaussian,exponential,markov} (gaussian) and --time-covariance-model (defaults to the spatial choice) — the decorrelation shape; the half-width flags set its size, defined so correlation = 0.5 at the half-width for every shape (markov = second-order autoregressive, (1 + r/L) e^{-r/L}), --swot-obs-error (0.002 m), --nadir-obs-error (0.02 m), --nadir-obs-error-by-mission JASON-3=0.03 SNTNL-3A=0.015 (per-mission overrides, case-insensitive match on the RADS mission field), --rads-thin-window (10; thins nadir data into blocks of N consecutive along-track points per pass — median ssha at the block mean position/median time, uncertainty max(--rads-thin-error-floor, std of the block), floor default 0.03 m; 0 disables), --max-neighbors (80), --min-neighbors (3), --max-points (200000), --search-radius-km (off; space-time nearest-neighbor selection by default), --fit-rads-pass-bias, --fit-swot-cross-swath-poly, --allow-land, --rads-missions j2 j3 3a 3b 6a, --rads-index-source {tables,auto,netcdf}. Runs can also be driven from JSON via --config-json (CLI flags override JSON).

Regional space-time cube (map_region_ssha.py)

Map SWOT Expert + RADS onto a regular lon/lat time cube in a bbox. The observation loaders and the OI / point-source tide numerics are the same as the single-time mapper and the swath-grid mapper; this driver does not OI a full SWOT pass and does not add a new kriging solver.

odsl-multi-altimeter-region \
    --bbox -126 34.5 -122 36.5 \
    --start 2024-02-01 --end 2024-02-10 \
    --dx-km 2 --dt-hours 6 \
    --out tmp/region_ssha_cube.nc

The spatial grid covers --bbox with posting --dx-km (default 2 km) at the bbox-center latitude (dlat = dx_km / 110.57, dlon = dx_km / (111.32 cos lat0)). OI distances still use the local azimuthal-equidistant transformer. Analysis times are UTC instants in [start, end] snapped to --dt-hours (default 6 h: 00/06/12/18). The load window is padded by max(--time-half-width-days, --tide-window-days) on each side. A 2nd-order spatial background and sequential harmonic constituents (default M2) are fit once; each frame stores ssh_anomaly_total = background + tide(t_k) + OI residual. Frames with more than ~500k lon/lat cells are refused unless --allow-large-grid is set.

Defaults that differ from odsl-multi-altimeter-oi: SWOT 5×5 median thinning with --swot-block-error std, nadir error 0.02 m and no RADS thinning, OI time half-width 4 days, SWOT quality <= 0, and internal_tide_hret restored on KaRIn so it matches RADS ssha. Output dims are time, lat, lon with ssh_anomaly_oi, oi_error_variance, observation_count, ssh_background, ssh_internal_tide, and ssh_anomaly_total.

Observation times

Every SWOT and RADS sample keeps its native measurement time. The loaders do not stamp observations with a granule-center, pass-center, or analysis-frame time. This is the same observation table used by odsl-multi-altimeter-oi, odsl-multi-altimeter-region, and odsl-multi-altimeter-swath-grid.

• RADS. Each along-track sample keeps the NetCDF time value. The regional cube default is no RADS thinning, so that timestamp is unchanged. --rads-thin-window N (used by the single-time mapper) tags the block with the median time of those N samples.
• SWOT Expert. L2 time is along-track (num_lines): one timestamp per line, broadcast across that line’s cross-track pixels. A KaRIn line is acquired nearly simultaneously; Expert has no per-pixel clock.
• SWOT N×N median (cube default 5×5). Each block is tagged with the center line’s time, or with the median of valid times in the block if the center pixel is invalid. A 5-line block is only a few seconds of flight.

Those per-point times t_i are what the later stages use:

• Internal tide — harmonic in time: A cos(kr − ω t_i) + B sin(...). The 2nd-order background polynomial is spatial-only and discards time on purpose.
• OI — space-time covariance at analysis time t0 (or each cube frame t_k) uses Δt = t_i − t0. Neighbors several days away are down-weighted; they are not treated as if they occurred at the map time.

On the regional cube the 00/06/12/18 axis is only the map time. The data keep their own clocks. See methodology.md §3.

Regional subset — moved to ../super_subsetter

The regional native-structure subsetter (odsl-super-subsetter) that harvests eleven /spray archives into one HDF5 file now lives in its own package. See ../super_subsetter/README.md.

Swath-grid SSHA (grid_swath_ssha.py)

Map one SWOT Expert segment plus nearby nadir altimetry on a grid that follows the SWOT swath but is widened in the cross-track direction:

odsl-multi-altimeter-swath-grid \
    --swot-file SWOT_L2_LR_SSH_Expert_525_013_20230519T001722_20230519T010829_PGD0_01.nc \
    --line-range 7068 7182 \
    --expansion-km 100 --delta-t-days 4 \
    --out tmp/swath_ssha_oi.nc

For a latitude-delimited regional run, replace --line-range with, for example, --lat-bounds 25 35. If neither option is supplied, the full pass is automatically processed as 320 km along-track segments with 20 km overlaps. Segmented products use deterministic names such as swath_ssha_oi.lines_00000_00160.nc.

• Reusable coordinates. The full expanded longitude/latitude grid is saved once per pass, expansion, and grid spacing under LOOKUP_DIR/expanded_swath_coordinates (or --swath-coordinates-dir). Later cycles of the same pass reuse that geometry while retaining their own observation times and SSHA.

• Grid. One row per along-track num_lines line; the cross-track axis is a uniform --grid-spacing-km (2 km) comb centered on the nadir track. By default it extends to ±(swath edge + --expansion-km) — ≈ ±168 km, 169 columns. --target-width-km 300 instead locks the total width to 300 km (±150 km, 151 columns at 2 km posting). Node positions are placed geodesically (WGS84) perpendicular to the local along-track direction, so descending passes and the dateline need no special handling. The signed cross_track_distance coordinate is in meters, negative = left of the flight direction (SWOT convention).

• Nadir search window. RADS observations are selected inside the expanded swath footprint and within ±--delta-t-days (4 days) of the SWOT segment mid time (or --target-time). The production 300 km driver uses ±5.5 days.
• Other SWOT passes. The local archive is searched for other Expert granules crossing the same footprint/time window and their points join the OI (disable with --no-other-swot-passes). For cal/val-phase granules pass --no-swot-nadir-prefilter: the RADS SWOT-nadir tracks use a different cycle numbering, so the prefilter would drop every candidate. The 300 km pass driver turns this on automatically for cycle >= 474.
• Output (dims n_lines × n_pixels): ssha_swot (granule data placed on its native columns; the drifting native comb means colliding pixels are averaged, see swot_pixel_count), ssha_nadir_mean/nadir_count (binned nadir observations), and the OI map ssh_anomaly_oi, oi_error_variance, observation_count. SWOT data cover only the native swath (nadir gap excluded by quality flags); the expansion wings are constrained by the nadir altimeters through the OI.
• Internal-tide decomposition (--fit-internal-tides, opt-in). Before the OI, the observations are decomposed: (1) a 2nd-order 2D polynomial background (weighted, time discarded); (2) a coherent fixed-frequency internal tide (--tide-constituent, M2 default) modeled as --tide-n-sources (6) point sources around the region center at equally spaced azimuths — each with its own fitted distance-to-center, wavelength, amplitude, and phase (a near source produces curved wavefronts, a far one nearly plane waves; amplitudes/phases solve linearly inside a bounded nonlinear search over distances/wavelengths with a greedy multi-start). Both fits use a wide ±--time-half-width-days-internal-tides (10 days) window; the OI then maps the residual within ±--delta-t-days. The output gains ssh_background, ssh_internal_tide (at the analysis time), and ssh_anomaly_total = background + tide + OI residual — the field comparable to ssha_swot — plus tide_* attrs with per-source parameters and variance explained. Bounds: --tide-wavelength-bounds (80–250 km), --tide-distance-bounds (default 1.5× region radius to 3000 km, keeping sources outside the data region). Fitting logic lives in multi_altimeters/internal_tides.py (pure numpy/scipy, reusable).
• Shares the OI flags (including --covariance-model / --time-covariance-model decorrelation shapes and --nadir-obs-error-by-mission per-mission errors), bias corrections (--fit-rads-pass-bias, --fit-swot-cross-swath-poly), the thinning options (--swot-median-window N×N SWOT median blocks, --rads-thin-window N-point nadir blocks — both with scatter-based errors), and --config-json handling with map_multi_altimeters.py. The OI solves one local system per grid point; full passes therefore run as bounded, overlapping segments by default.

Demonstration with figures (runs the pipeline once, reuses the NetCDF on rerun; writes tmp/demo_swath_grid_ssha_{swath,map,profile}.png):

python multi_altimeters/demo_swath_grid_ssha.py

Production 300 km pass mapping (grid_swath_ssha_pass.py)

Map one SWOT Expert pass at a time onto a 300 km total cross-track width (±150 km from the swath centerline — not 300 km per side, and not centered on the KaRIn nadir gap). The driver pins the production knobs and calls the same segmented swath-grid OI as odsl-multi-altimeter-swath-grid:

• SWOT Expert thinned 5×5; uncertainty is the sample std of valid pixels in the block (--swot-block-error std). Blocks with fewer than two valid pixels are skipped; uncertainty is max(sample std, 1e-4 m) so a numerically zero std cannot blow up the OI.
• Conventional RADS nadir uncertainty 0.02 m, no along-track thinning.
• OI time half-width 4 days; observation search window ±5.5 days.
• Other SWOT Expert passes that fall in the footprint/time window are included.
• Full passes are still split into 320 km / 20 km-overlap segments, then stitched (overlap cells weighted by distance from the segment ends).

odsl-swot-swath-grid-300km --cycle 10 --pass 317
odsl-swot-swath-grid-300km --start-cycle 10 --end-cycle 10 --passes 1 2 3

Default output is /spray/swot/regrid_swath_300km/cycle_CCC/pass_PPP_ssha_300km.nc (--out-root overrides). Existing pass files are skipped unless --overwrite. Do not drive production from grid_swath_ssha.config — that file disagrees with these knobs; the pass driver sets them on the engine CLI explicitly. Unknown flags are forwarded to the engine (for example --max-neighbors 40).

A one-segment diagnostic (California coast, cycle 10 pass 317, lines 6900–6960) writes tmp/demo_swath_grid_ssha_300km_{swath,map,profile}.png:

python multi_altimeters/demo_swath_grid_ssha.py \
    --target-width-km 300 --swot-median-window 5 --swot-block-error std \
    --rads-thin-window 0 --nadir-obs-error 0.02 \
    --time-half-width-days 4 --delta-t-days 5.5 \
    --swot-nadir-prefilter \
    --grid-file tmp/demo_swath_grid_ssha_300km.nc \
    --out-prefix tmp/demo_swath_grid_ssha_300km

Configured example and SWOT thinning demo:

python multi_altimeters/example_oi_gulf_america_20250110.py

python multi_altimeters/demo_plot_thinned_swot.py \
    --start 2025-01-01 --end 2025-01-10 --bbox -98 18 -80 31 \
    --median-window 3 \
    --out tmp/swot_thinned_n3.png --table-out tmp/swot_thinned_n3.csv

Point-vortex inversion (reads point_vortex_ssh_inversion_config.json by default):

# Baseline iteration-0 fit
python multi_altimeters/point_vortex_ssh_iteration0.py \
    --config multi_altimeters/point_vortex_ssh_inversion_config.json

# Quick coarse override
python multi_altimeters/point_vortex_ssh_iteration0.py \
    --sources rads --grid-spacing 3 --search-radius-deg 15 \
    --max-observations 100 --max-iterations 5 \
    --out tmp/test_point_vortex.nc --out-png tmp/test_point_vortex.png

# Fit point vortices from RADS near a chosen SWOT segment
python multi_altimeters/point_vortex_ssh_from_swot_segment.py \
    --swot-file SWOT_L2_LR_SSH_Expert_..._nc \
    --line-range 0 2000 --delta-t-days 5

Conventions & notes

• Bounding boxes are west south east north in degrees; both [-180,180) and [0,360) longitudes are accepted and normalized.
• The OI solves a dense local kriging system per grid point. For each cell it gathers up to --max-neighbors space-time nearest observations, builds a spatial × temporal covariance from the selected decorrelation model (--covariance-model, default Gaussian; the half-width flags fix the scale so correlation = 0.5 at the half-width for every shape), adds uncertainty^2 to the diagonal, augments with the ordinary-kriging Lagrange row/column, and solves (falling back to least squares if singular). Cells with fewer than --min-neighbors are left unfilled.
• Observation errors default to 0.002 m (SWOT) and 0.02 m (nadir RADS), with per-mission overrides via --nadir-obs-error-by-mission. When thinning is enabled, uncertainties become scatter-based: each block carries the standard deviation of its member values, floored at --swot-obs-error (SWOT) or --rads-thin-error-floor (nadir, 0.03 m). SWOT defaults to ssha_karin_2 + height_cor_xover + internal_tide_hret, keeping ssha_karin_2_qual == 0 and crossover quality <= 1; RADS defaults to ssha. The internal tide is added back to SWOT by default (--no-swot-internal-tide disables): the L2 product subtracts HRET from ssha_karin_2, whereas the RADS sla formula never removes tide_internal (see /spray/rads/data/conf/rads.xml), so restoring it keeps both sources carrying the same internal-tide signal. Ocean-only filtering is on unless --allow-land.
• SWOT thinning: --swot-median-window N (odd N) reduces the KaRIn swath by non-overlapping N×N median blocks before mapping — each block contributes its median ssha. --swot-block-error floor (default) sets the uncertainty to max(--swot-obs-error, block std); std uses the sample std alone (skipping blocks with fewer than two valid pixels). SWOT file search uses RADS SWOT-nadir tracks as a prefilter (0.5° bbox padding).
• Nadir thinning: --rads-thin-window N (default 10) groups N consecutive along-track points per pass into one observation — median ssha at the block's mean position and median time, error max(--rads-thin-error-floor, block std).
• Point-vortex: iteration-0 ignores time except for input selection; the streamfunction kernel is softened at small radius; the fit is an exploratory basis-function inversion, not a full dynamical assimilation. iteration1.py is under development.
• Lookup tables in lookup_tables/ (RADS/SWOT file indices) are created on first run; rebuild with --refresh-index. --rads-missions are encoded into the RADS cache filename (e.g. rads_file_index_j2_j3.csv).
• Outputs go to --out / tmp/, never into the source tree (data files are gitignored).

Related

• ../common — shared swath→grid regridding and common.plotting.add_map_context.
• ../swot — SWOT L2 LR SSH Expert loading (swot/utils.py).
• ../rads — RADS nadir altimetry loading and indexing (rads/utils.py).
• ../README.md — repo-wide conventions (CLI structure, bbox/longitude, data roots, extras).
• point_vortex_ssh_readme.md — point-vortex iteration status notes.