drifterOI

Python code for 3D (space + time) optimal interpolation of sparse Lagrangian drifter velocities, using Gauss-Markov / Gaussian Process regression. Splits observed velocities into:

• a seasonal component: a smooth annual harmonic cycle per location, fit as a separable spatial-Matern x per-harmonic Gauss-Markov model (drifterOI.seasonal_gp, solved efficiently via datide.MultiscaleGP's block-diagonal decomposition rather than a dense joint solve), and
• an eddy/synoptic residual: fit as a Vecchia-approximated Gaussian Process over the seasonal-corrected residuals (via gpvecchia), to scale to the full drifter dataset without a dense covariance solve.

This is the companion code for the accompanying manuscript, submitted to Ocean Science: [preprint link]. See Citation below.

Installation

Local development install (no PyPI release yet):

git clone https://github.com/williamedge/drifterOI.git
cd drifterOI
pip install -e .

Requires Python >= 3.10 (set by the gpvecchia dependency's floor). Besides the usual scientific Python stack, this pulls in three sibling packages from the same research group: gpvecchia (Vecchia-approximated GPs), datide (Gauss-Markov / data assimilation routines), and gptide (GP kernels, MLE/MCMC fitting).

Quickstart

A minimal example of the separable seasonal Gauss-Markov model on synthetic data -- the same model seasonal_gp.init_gm_scales builds for the real pipeline, without needing the drifter dataset:

import numpy as np
from datide import MultiscaleGP

import drifterOI.seasonal_gp as sgp
from drifterOI.cov import matern32_covfunc

rng = np.random.default_rng(0)
n_obs = 200
X = np.column_stack([
    rng.uniform(-100, 100, n_obs),   # easting (km)
    rng.uniform(-100, 100, n_obs),   # northing (km)
    np.sort(rng.uniform(0, 200, n_obs)),  # time (days)
])
Y = rng.normal(size=n_obs)  # observed velocity component

# [sigma, lx, ly, lt]: spatial Matern prior for the seasonal coefficients
B_params = [1.0, 50.0, 50.0, 5.0]
R = np.eye(n_obs) * 0.05  # observation noise covariance

H_list, scales, design = sgp.init_gm_scales(X, matern32_covfunc, B_params, mean=True)
model = MultiscaleGP(H_list, scales, R)

param_stack = model.conditional_mean(Y)          # (n_params, n_obs) harmonic coefficients
seasonal_fit = sgp.apply_harmonic_design(design, param_stack)
eddy_residual = Y - seasonal_fit

Pipeline

The scripts/ directory holds the full analysis pipeline behind the manuscript, in stages (each stage reads the previous stage's saved output):

Example data (data/) is included in the repository so the L1-L2 scripts can be run directly (BODY_drifters_6_hours.nc drifter observations, SUNTANS_bathy_unstructured.nc bathymetry for plotting). data/, scripts/, and tests/ aren't part of the installable package (pip install drifterOI won't fetch them) -- clone the repo if you want the example data or to run the pipeline scripts yourself.

Citation

If you use this code, please cite it as described in CITATION.cff. Once published, the accompanying Ocean Science manuscript will be the preferred citation -- see [preprint link] in the meantime.

License

BSD-3-Clause -- see LICENSE.