Using SWOT data

Show/hide code

%load_ext autoreload 
%autoreload 2

Show/hide code

import cartopy.crs as ccrs
import cmocean.cm as cmo
import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt
import numpy as np
import optax
import xarray as xr

from jaxparrow.cyclogeostrophy import cyclogeostrophic_imbalance, gradient_wind, minimization_based
from jaxparrow.utils import kinematics, operators


jax.config.update("jax_enable_x64", True)

This notebook demonstrates how regularization can be added to the minimization-based formulation of the cyclogeostrophic inversion problem.

SWOT observation

We use the filtered \(u\) and \(v\) component of the geostrophic currents, distributed in version 3.0 of the product from Aviso. And we focus for the illustration on the Balearic Sea.

Show/hide code

def make_da(data, ds):
    return xr.DataArray(
        data, coords=ds.ssha_filtered.coords, dims=ds.ssha_filtered.dims
    )


def get_med_sea(swot_ds):
    med_mask = (
        (swot_ds["latitude"] > 30.0681) & (swot_ds["latitude"] < 47.3764) & 
        (swot_ds["longitude"] > -6.0326) & (swot_ds["longitude"] < 42.355)
    )
    swot_ds = swot_ds.where(med_mask, drop=True)

    sea_mask = ~np.isnan(swot_ds["ssha_filtered"])
    min_lat, max_lat = swot_ds["latitude"].where(sea_mask).min(), swot_ds["latitude"].where(sea_mask).max()
    min_lon, max_lon = swot_ds["longitude"].where(sea_mask).min(), swot_ds["longitude"].where(sea_mask).max()

    swot_ds = swot_ds.where(
        (swot_ds["latitude"] >= min_lat) & (swot_ds["latitude"] <= max_lat) &
        (swot_ds["longitude"] >= min_lon) & (swot_ds["longitude"] <= max_lon),
        drop=True
    )

    return swot_ds


swot_003_ds = xr.open_dataset("data/SWOT/SWOT_L3_LR_SSH_Expert_506_003_20230429T184429_20230429T193535_v3.0.nc")
swot_003_ds = get_med_sea(swot_003_ds)

swot_016_ds = xr.open_dataset("data/SWOT/SWOT_L3_LR_SSH_Expert_506_016_20230430T054843_20230430T063947_v3.0.nc")
swot_016_ds = get_med_sea(swot_016_ds)

lat_003 = jnp.asarray(swot_003_ds.latitude.values)
lon_003 = jnp.asarray(swot_003_ds.longitude.values)
ug_003 = jnp.asarray(swot_003_ds.ugos_filtered.values)
vg_003 = jnp.asarray(swot_003_ds.vgos_filtered.values)

lat_016 = jnp.asarray(swot_016_ds.latitude.values)
lon_016 = jnp.asarray(swot_016_ds.longitude.values)
ug_016 = jnp.asarray(swot_016_ds.ugos_filtered.values)
vg_016 = jnp.asarray(swot_016_ds.vgos_filtered.values)

Show/hide code

uvg_003_da = (swot_003_ds.ugos_filtered**2 + swot_003_ds.vgos_filtered**2)**0.5

rvg_003 = kinematics.vorticity(ug_003, vg_003, lat_003, lon_003)
rvg_003_da = make_da(rvg_003, swot_003_ds)

cg_uimbg_003, cg_vimbg_003 = cyclogeostrophic_imbalance(ug_003, vg_003, ug_003, vg_003, lat_003, lon_003)
cg_imbg_003 = (cg_uimbg_003**2 + cg_vimbg_003**2)**0.5
cg_imbg_003_da = make_da(cg_imbg_003, swot_003_ds)

uvg_016_da = (swot_016_ds.ugos_filtered**2 + swot_016_ds.vgos_filtered**2)**0.5

rvg_016 = kinematics.vorticity(ug_016, vg_016, lat_016, lon_016)
rvg_016_da = make_da(rvg_016, swot_016_ds)

cg_uimbg_016, cg_vimbg_016 = cyclogeostrophic_imbalance(ug_016, vg_016, ug_016, vg_016, lat_016, lon_016)
cg_imbg_016 = (cg_uimbg_016**2 + cg_vimbg_016**2)**0.5
cg_imbg_016_da = make_da(cg_imbg_016, swot_016_ds)

We can see that in the version 3.0 geostrophic currents are relatively smooth, and do not look too contaminated by noise and/or unbalanced signal.

Show/hide code

uvg_max = np.nanquantile(jnp.concatenate([uvg_003_da.values.flatten(), uvg_016_da.values.flatten()]), 0.98)
rvg_max = np.nanquantile(jnp.concatenate([rvg_003_da.values.flatten(), rvg_016_da.values.flatten()]), 0.98)
imbg_max = np.nanquantile(jnp.concatenate([cg_imbg_003_da.values.flatten(), cg_imbg_016_da.values.flatten()]), 0.98)

fig, axes = plt.subplots(1, 3, figsize=(20, 6), subplot_kw={"projection": ccrs.PlateCarree()})

_ = uvg_003_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.speed, add_colorbar=True, vmin=0, vmax=uvg_max,
    cbar_kwargs={"label": "Speed (m s$^{-1}$)", "extend": "max"},
    transform=ccrs.PlateCarree(), ax=axes[0]
)
_ = uvg_016_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.speed, add_colorbar=False, vmin=0, vmax=uvg_max,
    transform=ccrs.PlateCarree(), ax=axes[0]
)

_ = rvg_003_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.curl, add_colorbar=True, vmin=-rvg_max, vmax=rvg_max, 
    cbar_kwargs={"label": "Vorticity", "extend": "both"},
    transform=ccrs.PlateCarree(), ax=axes[1]
)

_ = rvg_016_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.curl, add_colorbar=False, vmin=-rvg_max, vmax=rvg_max, 
    transform=ccrs.PlateCarree(), ax=axes[1]
)

_ = cg_imbg_003_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.amp, add_colorbar=True, vmin=0, vmax=imbg_max,
    cbar_kwargs={"label": "Cyclogeostrophic imbalance (m s$^{-1}$)", "extend": "max"},
    transform=ccrs.PlateCarree(), ax=axes[2]
)

_ = cg_imbg_016_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.amp, add_colorbar=False, vmin=0, vmax=imbg_max,
    transform=ccrs.PlateCarree(), ax=axes[2]
)

for ax in axes:
    ax.coastlines()

fig.suptitle("SWOT L3 Geostrophic Currents")

plt.show()

However, when looking at the cyclogeostrophic currents estimated using the gradient_wind method we can notice a few empty areas, where the azimuthal cyclogeostrophic equation does not have a physical solution.

Show/hide code

gw_003_results = gradient_wind(lat_t=lat_003, lon_t=lon_003, ug_t=ug_003, vg_t=vg_003)

ucg_gw_003 = gw_003_results.ucg
vcg_gw_003 = gw_003_results.vcg

uvcg_gw_003 = kinematics.magnitude(ucg_gw_003, vcg_gw_003)
rvcg_gw_003 = kinematics.vorticity(ucg_gw_003, vcg_gw_003, lat_003, lon_003)

uvcg_gw_003_da = make_da(uvcg_gw_003, swot_003_ds)
rvcg_gw_003_da = make_da(rvcg_gw_003, swot_003_ds)

cg_uimb_gw_003, cg_vimb_gw_003 = cyclogeostrophic_imbalance(ug_003, vg_003, ucg_gw_003, vcg_gw_003, lat_003, lon_003)
cg_imb_gw_003 = (cg_uimb_gw_003**2 + cg_vimb_gw_003**2)**0.5
cg_imb_gw_003_da = make_da(cg_imb_gw_003, swot_003_ds)

gw_016_results = gradient_wind(lat_t=lat_016, lon_t=lon_016, ug_t=ug_016, vg_t=vg_016)

ucg_gw_016 = gw_016_results.ucg
vcg_gw_016 = gw_016_results.vcg

uvcg_gw_016 = kinematics.magnitude(ucg_gw_016, vcg_gw_016)
rvcg_gw_016 = kinematics.vorticity(ucg_gw_016, vcg_gw_016, lat_016, lon_016)

uvcg_gw_016_da = make_da(uvcg_gw_016, swot_016_ds)
rvcg_gw_016_da = make_da(rvcg_gw_016, swot_016_ds)

cg_uimb_gw_016, cg_vimb_gw_016 = cyclogeostrophic_imbalance(ug_016, vg_016, ucg_gw_016, vcg_gw_016, lat_016, lon_016)
cg_imb_gw_016 = (cg_uimb_gw_016**2 + cg_vimb_gw_016**2)**0.5
cg_imb_gw_016_da = make_da(cg_imb_gw_016, swot_016_ds)

Show/hide code

fig, axes = plt.subplots(1, 3, figsize=(20, 6), subplot_kw={"projection": ccrs.PlateCarree()})

_ = uvcg_gw_003_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.speed, add_colorbar=True, vmin=0, vmax=uvg_max,
    cbar_kwargs={"label": "Speed (m s$^{-1}$)", "extend": "max"},
    transform=ccrs.PlateCarree(), ax=axes[0]
)

_ = uvcg_gw_016_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.speed, add_colorbar=False, vmin=0, vmax=uvg_max,
    transform=ccrs.PlateCarree(), ax=axes[0]
)

_ = rvcg_gw_003_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.curl, add_colorbar=True, vmin=-rvg_max, vmax=rvg_max,
    cbar_kwargs={"label": "Vorticity", "extend": "both"},
    transform=ccrs.PlateCarree(), ax=axes[1]
)

_ = rvcg_gw_016_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.curl, add_colorbar=False, vmin=-rvg_max, vmax=rvg_max,
    transform=ccrs.PlateCarree(), ax=axes[1]
)

_ = cg_imb_gw_003_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.amp, add_colorbar=True, vmin=0, vmax=imbg_max,
    cbar_kwargs={"label": "Cyclogeostrophic imbalance (m s$^{-1}$)", "extend": "max"},
    transform=ccrs.PlateCarree(), ax=axes[2]
)

_ = cg_imb_gw_016_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.amp, add_colorbar=False, vmin=0, vmax=imbg_max,
    transform=ccrs.PlateCarree(), ax=axes[2]
)

for ax in axes:
    ax.coastlines()

fig.suptitle("Gradient-Wind (GW) Cyclogeostrophic Currents")

plt.show()

In those same regions, the minimization_based fields presents some suspicious patterns, and the cyclogeostrophic imbalance is not closed even when optimizing the currents with respect to its gradient.

But with this approach, we can also add a regularization term to the loss function being minimized!

Show/hide code

mb_003_results = minimization_based(
    lat_t=lat_003, lon_t=lon_003, ug_t=ug_003, vg_t=vg_003,
    optim=optax.chain(optax.clip(1), optax.sgd(learning_rate=5e-3))
)

ucg_mb_003 = mb_003_results.ucg
vcg_mb_003   = mb_003_results.vcg

uvcg_mb_003 = kinematics.magnitude(ucg_mb_003, vcg_mb_003)
rvcg_mb_003 = kinematics.vorticity(ucg_mb_003, vcg_mb_003, lat_003, lon_003)

uvcg_mb_003_da = make_da(uvcg_mb_003, swot_003_ds)
rvcg_mb_003_da = make_da(rvcg_mb_003, swot_003_ds)

cg_uimb_mb_003, cg_vimb_mb_003 = cyclogeostrophic_imbalance(ug_003, vg_003, ucg_mb_003, vcg_mb_003, lat_003, lon_003)
cg_imb_mb_003 = (cg_uimb_mb_003**2 + cg_vimb_mb_003**2)**0.5
cg_imb_mb_003_da = make_da(cg_imb_mb_003, swot_003_ds)

mb_016_results = minimization_based(
    lat_t=lat_016, lon_t=lon_016, ug_t=ug_016, vg_t=vg_016,
    optim=optax.chain(optax.clip(1), optax.sgd(learning_rate=5e-3))
)

ucg_mb_016 = mb_016_results.ucg
vcg_mb_016   = mb_016_results.vcg

uvcg_mb_016 = kinematics.magnitude(ucg_mb_016, vcg_mb_016)
rvcg_mb_016 = kinematics.vorticity(ucg_mb_016, vcg_mb_016, lat_016, lon_016)

uvcg_mb_016_da = make_da(uvcg_mb_016, swot_016_ds)
rvcg_mb_016_da = make_da(rvcg_mb_016, swot_016_ds)

cg_uimb_mb_016, cg_vimb_mb_016 = cyclogeostrophic_imbalance(ug_016, vg_016, ucg_mb_016, vcg_mb_016, lat_016, lon_016)
cg_imb_mb_016 = (cg_uimb_mb_016**2 + cg_vimb_mb_016**2)**0.5
cg_imb_mb_016_da = make_da(cg_imb_mb_016, swot_016_ds)

Show/hide code

fig, axes = plt.subplots(1, 3, figsize=(20, 6), subplot_kw={"projection": ccrs.PlateCarree()})

_ = uvcg_mb_003_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.speed, add_colorbar=True, vmin=0, vmax=uvg_max,
    cbar_kwargs={"label": "Speed (m s$^{-1}$)", "extend": "max"},
    transform=ccrs.PlateCarree(), ax=axes[0]
)

_ = uvcg_mb_016_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.speed, add_colorbar=False, vmin=0, vmax=uvg_max,
    transform=ccrs.PlateCarree(), ax=axes[0]
)

_ = rvcg_mb_003_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.curl, add_colorbar=True, vmin=-rvg_max, vmax=rvg_max,
    cbar_kwargs={"label": "Vorticity", "extend": "both"},
    transform=ccrs.PlateCarree(), ax=axes[1]
)

_ = rvcg_mb_016_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.curl, add_colorbar=False, vmin=-rvg_max, vmax=rvg_max,
    transform=ccrs.PlateCarree(), ax=axes[1]
)

_ = cg_imb_mb_003_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.amp, add_colorbar=True, vmin=0, vmax=imbg_max,
    cbar_kwargs={"label": "Cyclogeostrophic imbalance (m s$^{-1}$)", "extend": "max"},
    transform=ccrs.PlateCarree(), ax=axes[2]
)

_ = cg_imb_mb_016_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.amp, add_colorbar=False, vmin=0, vmax=imbg_max,
    transform=ccrs.PlateCarree(), ax=axes[2]
)

for ax in axes:
    ax.coastlines()

fig.suptitle("Minimization-Based (MB) Cyclogeostrophic Currents")

plt.show()

Regularization

Adding a regularization term simply resolves to writing a JAX-friendly regularization function. The signature (i.e. the parameters) of the regularization function must contain ucg_t and vcg_t as arguments for the components of the cyclogeostrophic velocity field, and lat_t and lon_t as arguments for the grid coordinate arrays. We can also used some additional arguments automatically passed to the function if present in its signature, such as dx_t, dy_t, and land_mask to avoid recomputing them at every minimization step, but those are not mandatory. The order of the arguments does not matter, but their exact spelling does. We refer the reader to the API documentation of the minimization_based function for the comprehensive description of the possible signature of the regularization function.

Here we use a biharmonic-like regularization on the SSH by penalizing large vorticity horizontal gradients. For this we use the vorticity and horizontal_derivatives functions of the public API.

In order to focus on the regions where it is likely that cyclogeostrophy is not a valid assumption, we apply a weight to each grid cell used in the computation of the regularization term, and this weight is taken to be the cyclogeostrophic imbalance after estimating cyclogeostrophic currents without a regularization term. The weights are provided to the regularisation function as an extra parameter, arbitrarly named reg_weight.

We can provide any valid Pytree as additional argument to the function, but if the grid is not rectilinear, vector fields should very probably be represented in grid coordinates for consistency with ucg_t and vcg_t. This can be achieved using the function rotate_to_grid.

Show/hide code

lambda_reg = 1e1

def biharmonic_like_regularization(ucg_t, vcg_t, lat_t, lon_t, dx_t, dy_t, land_mask, reg_weight):
    vort = kinematics.vorticity(ucg_t, vcg_t, lat_t, lon_t)
    dvort_x, dvort_y = operators.horizontal_derivatives(vort, dx=dx_t, dy=dy_t, land_mask=land_mask)
    reg = (dvort_x * dx_t)**2 + (dvort_y * dy_t)**2
    reg *= reg_weight
    return lambda_reg * jnp.nansum(reg)

We then pass this function as an argument to minimization_based when performing the inversion, along with the minimization-based cyclogeostrophic imbalance computed previously.

Show/hide code

mb_reg_003_results = minimization_based(
    lat_t=lat_003, lon_t=lon_003, ug_t=ug_003, vg_t=vg_003, 
    optim=optax.chain(optax.clip(1), optax.sgd(learning_rate=5e-3)), 
    regularization=biharmonic_like_regularization, reg_kwargs={"reg_weight": cg_imb_mb_003}
)

mb_reg_016_results = minimization_based(
    lat_t=lat_016, lon_t=lon_016, ug_t=ug_016, vg_t=vg_016, 
    optim=optax.chain(optax.clip(1), optax.sgd(learning_rate=5e-3)), 
    regularization=biharmonic_like_regularization, reg_kwargs={"reg_weight": cg_imb_mb_016}
)

Show/hide code

ucg_mb_reg_003 = mb_reg_003_results.ucg
vcg_mb_reg_003 = mb_reg_003_results.vcg

uvcg_mb_reg_003 = kinematics.magnitude(ucg_mb_reg_003, vcg_mb_reg_003)
rvcg_mb_reg_003 = kinematics.vorticity(ucg_mb_reg_003, vcg_mb_reg_003, lat_003, lon_003)

uvcg_mb_reg_003_da = make_da(uvcg_mb_reg_003, swot_003_ds)
rvcg_mb_reg_003_da = make_da(rvcg_mb_reg_003, swot_003_ds)

cg_uimb_mb_reg_003, cg_vimb_mb_reg_003 = cyclogeostrophic_imbalance(
    ug_003, vg_003, ucg_mb_reg_003, vcg_mb_reg_003, lat_003, lon_003
)
cg_imb_mb_reg_003 = (cg_uimb_mb_reg_003**2 + cg_vimb_mb_reg_003**2)**0.5
cg_imb_mb_reg_003_da = make_da(cg_imb_mb_reg_003, swot_003_ds)

ucg_mb_reg_016 = mb_reg_016_results.ucg
vcg_mb_reg_016 = mb_reg_016_results.vcg

uvcg_mb_reg_016 = kinematics.magnitude(ucg_mb_reg_016, vcg_mb_reg_016)
rvcg_mb_reg_016 = kinematics.vorticity(ucg_mb_reg_016, vcg_mb_reg_016, lat_016, lon_016)

uvcg_mb_reg_016_da = make_da(uvcg_mb_reg_016, swot_016_ds)
rvcg_mb_reg_016_da = make_da(rvcg_mb_reg_016, swot_016_ds)

cg_uimb_mb_reg_016, cg_vimb_mb_reg_016 = cyclogeostrophic_imbalance(
    ug_016, vg_016, ucg_mb_reg_016, vcg_mb_reg_016, lat_016, lon_016
)
cg_imb_mb_reg_016 = (cg_uimb_mb_reg_016**2 + cg_vimb_mb_reg_016**2)**0.5
cg_imb_mb_reg_016_da = make_da(cg_imb_mb_reg_016, swot_016_ds)

The result is visually very similar to the solution obtain without the regularization, but when looking closely we can notice some differences, in particular in regions where the balance was not closed without regularization.

Show/hide code

fig, axes = plt.subplots(1, 3, figsize=(20, 6), subplot_kw={"projection": ccrs.PlateCarree()})

_ = uvcg_mb_reg_003_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.speed, add_colorbar=True, vmin=0, vmax=uvg_max,
    cbar_kwargs={"label": "Speed (m s$^{-1}$)", "extend": "max"},
    transform=ccrs.PlateCarree(), ax=axes[0]
)

_ = uvcg_mb_reg_016_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.speed, add_colorbar=False, vmin=0, vmax=uvg_max,
    transform=ccrs.PlateCarree(), ax=axes[0]
)

_ = rvcg_mb_reg_003_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.curl, add_colorbar=True, vmin=-rvg_max, vmax=rvg_max,
    cbar_kwargs={"label": "Vorticity", "extend": "both"},
    transform=ccrs.PlateCarree(), ax=axes[1]
)

_ = rvcg_mb_reg_016_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.curl, add_colorbar=False, vmin=-rvg_max, vmax=rvg_max,
    transform=ccrs.PlateCarree(), ax=axes[1]
)

_ = cg_imb_mb_reg_003_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.amp, add_colorbar=True, vmin=0, vmax=imbg_max,
    cbar_kwargs={"label": "Cyclogeostrophic imbalance (m s$^{-1}$)", "extend": "max"},
    transform=ccrs.PlateCarree(), ax=axes[2]
)

_ = cg_imb_mb_reg_016_da.plot.pcolormesh(
    x="longitude", y="latitude", 
    cmap=cmo.amp, add_colorbar=False, vmin=0, vmax=imbg_max,
    transform=ccrs.PlateCarree(), ax=axes[2]
)

for ax in axes:
    ax.coastlines()

fig.suptitle("Regularized Minimization-Based (MB - Reg.) Cyclogeostrophic Currents")

plt.show()

Comparison

To better observe the qualitative discrepancy between the inversion results, we can visualize the difference of the fieds in addition to the fields themselfs.

Show/hide code

methods = ["Geos.", "GW", "MB", "MB - Reg."]

row_labels = methods[:-1]  # subtracted from
col_labels = methods[1:]   # subtracted


def build_diff_da(fields, ds):
    diff_data = np.full((len(row_labels), len(col_labels)) + fields[0].shape, np.nan)
    for i in range(len(row_labels)):
        for j in range(len(col_labels)):
            if i <= j:  # upper triangle
                diff_data[i, j] = np.array(fields[j + 1] - fields[i])
    return xr.DataArray(
        diff_data,
        dims=("row", "col", "num_lines", "num_pixels"),
        coords={
            "row": row_labels, "col": col_labels,
            "latitude": (("num_lines", "num_pixels"), ds["latitude"].values),
            "longitude": (("num_lines", "num_pixels"), ds["longitude"].values),
        },
    )


def plot_diff_triangle(da_003, da_016, cmap, cbar_label):
    da = xr.concat([da_003, da_016], dim="num_lines")

    fg = da.plot.pcolormesh(
        x="longitude", y="latitude", row="row", col="col",
        cmap=cmap, robust=True, add_colorbar=True,
        cbar_kwargs={"label": cbar_label, "aspect": 50},
        figsize=(13, 12),
        subplot_kws={"projection": ccrs.PlateCarree()},
    )
    fg.set_titles("{value}")
    for i, ax in enumerate(fg.axs.flat):
        row_idx, col_idx = divmod(i, len(col_labels))
        if row_idx > col_idx:
            ax.set_visible(False)
        else:
            ax.coastlines()
    plt.show()


uvcg_diff_003_da = build_diff_da([uvg_003_da, uvcg_gw_003_da, uvcg_mb_003_da, uvcg_mb_reg_003_da], swot_003_ds)
uvcg_diff_016_da = build_diff_da([uvg_016_da, uvcg_gw_016_da, uvcg_mb_016_da, uvcg_mb_reg_016_da], swot_016_ds)
rv_diff_003_da = build_diff_da([rvg_003_da, rvcg_gw_003_da, rvcg_mb_003_da, rvcg_mb_reg_003_da], swot_003_ds)
rv_diff_016_da = build_diff_da([rvg_016_da, rvcg_gw_016_da, rvcg_mb_016_da, rvcg_mb_reg_016_da], swot_016_ds)

In the speed difference fields of the regularization variant we clearly see the impact of the regularization in the regions where the gradient-wind solution does not exist and where the cyclogeostrophic imbalance was large.

Show/hide code

plot_diff_triangle(uvcg_diff_003_da, uvcg_diff_016_da, cmo.balance, "Speed difference (m/s)")

Those corrections are more pronounced in the vorticity fields, again in the areas where the gradient-wind reconstruction and the cyclogeostrophic imbalance suggest that cyclogeostrophy is not valid.

Show/hide code

plot_diff_triangle(rv_diff_003_da, rv_diff_016_da, cmo.balance, "Relative vorticity difference")