Grid operations, regridding and sigma coordinates

Grid operations, regridding and sigma coordinates#

Introduction#

The xoa library provides comprehensive tools for working with grids, performing regridding operations, and handling terrain-following parametric vertical coordinates (sigma coordinates). This guide covers three main modules:

  • xoa.grid: Grid utilities for 1D to nD grid operations on a single grid

  • xoa.regrid: Regridding utilities for operations between different grids

  • xoa.sigma: Terrain-following parametric vertical coordinates following CF conventions

Grid operations with xoa.grid#

The xoa.grid module provides utilities to get information or perform operations on a grid.

Note

For operations between different grids, use xoa.regrid instead.

Basic grid operations#

The xoa.grid module offers various functions to manipulate grids:

Getting grid centers and edges:

In [1]: import xarray as xr

In [2]: import numpy as np

In [3]: import xoa.grid as xgrid

# Create a simple 1D coordinate
In [4]: x_edges = xr.DataArray(np.arange(0, 5), dims='x', name='x')

# Get centers from edges
In [5]: x_centers = xgrid.get_centers(x_edges, "x")

In [6]: print(x_centers.values)
[0.5 1.5 2.5 3.5]

# Get edges from centers
In [7]: x_back = xgrid.get_edges(x_centers, "x", mode="linear_extrap")

In [8]: print(x_back.values)
[0. 1. 2. 3. 4.]

Applying operations along dimensions:

The apply_along_dim() function allows you to apply an xarray operator on data array or dataset dimensions, potentially changing the size of the array.

# Create sample dataset
In [9]: ds = xr.Dataset({
   ...:     'temp': (['x', 'y'], np.random.rand(4, 3)),
   ...: }, coords={
   ...:     'x': np.arange(4),
   ...:     'y': np.arange(3)
   ...: })
   ...: 

# Interpolate between grid points
In [10]: ds_mean = xgrid.apply_along_dim(ds, 'x', xgrid.get_centers)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[10], line 1
----> 1 ds_mean = xgrid.apply_along_dim(ds, 'x', xgrid.get_centers)

File ~/checkouts/readthedocs.org/user_builds/xoa/conda/v2026.2.0/lib/python3.14/site-packages/xoa/grid.py:115, in apply_along_dim(ds, dim, func, coord_func, data_kwargs, coord_kwargs, name_kwargs, **kwargs)
    113 if isinstance(dso, xr.Dataset):
    114     dso = dso.update(daos)
--> 115     dso.attrs = ds.attrs
    116     dso.encoding = ds.encoding
    117 else:

AttributeError: 'NoneType' object has no attribute 'attrs' and no __dict__ for setting new attributes

Staggered grids#

The xoa library supports staggered grids commonly used in ocean and atmospheric models. Grid location information can be encoded in variable names using the staggered grid locator (see Metadata, CF and naming conventions for more details on the SGLocator).

Regridding with xoa.regrid#

The xoa.regrid module provides utilities for regridding data from one grid to another. It supports various interpolation methods for both horizontal and vertical regridding.

1D regridding#

The regrid1d() function performs 1D interpolation/regridding with support for multiple methods:

Supported methods#

In [11]: from xoa.regrid import regrid1d_methods

In [12]: print(list(regrid1d_methods))
[<regrid1d_methods.linear: 1>, <regrid1d_methods.nearest: 0>, <regrid1d_methods.cubic: 2>, <regrid1d_methods.hermit: 3>, <regrid1d_methods.cellave: -1>, <regrid1d_methods.cellerr: -2>]

These include:

  • linear: Linear interpolation (default)

  • nearest: Nearest neighbor interpolation

  • cubic: Cubic interpolation

  • hermit/hermitian: Hermitian interpolation

  • cellave: Cell-averaging or conservative regridding

Example usage#

In [13]: from xoa import regrid

# Create source data
In [14]: x_old = xr.DataArray(np.arange(0, 10, 2.0), dims='x')

In [15]: data_old = xr.DataArray(np.sin(x_old), dims='x', coords={'x': x_old})

# Create new grid
In [16]: x_new = xr.DataArray(np.arange(0, 9, 0.5), dims='x')

# Regrid with linear interpolation
In [17]: data_new = regrid.regrid1d(data_old, x_new, method='linear')

Extrapolation modes#

You can control extrapolation behavior:

In [18]: from xoa.regrid import extrap_modes

In [19]: print(list(extrap_modes))
[<extrap_modes.no: 0>, <extrap_modes.top: 1>, <extrap_modes.bottom: -1>, <extrap_modes.both: 2>]

Available modes include:

  • no/none/false: No extrapolation (default)

  • top/above/after: Extrapolate toward the top

  • bottom/below: Extrapolate toward the bottom

  • both/all/yes/true: Extrapolate both directions

Horizontal regridding#

For 2D horizontal regridding, xoa provides the grid2loc() function to interpolate from a grid to specific locations.

Vertical regridding#

Vertical regridding is particularly important for ocean models. The isoslice() function can extract data on isobaric or isopycnal surfaces.

Sigma coordinates with xoa.sigma#

The xoa.sigma module handles terrain-following parametric vertical coordinates according to the CF conventions for Parametric Vertical Coordinates.

What are sigma coordinates?#

Sigma coordinates are terrain-following vertical coordinates commonly used in ocean and atmospheric models. Instead of using fixed depth or pressure levels, sigma coordinates follow the bathymetry (ocean models) or topography (atmospheric models), which allows better resolution near boundaries.

Supported coordinate types#

The module supports the following CF-compliant sigma coordinate types:

In [20]: from xoa.sigma import SIGMA_COORDINATE_FUNCTIONS

In [21]: for coord_type in SIGMA_COORDINATE_FUNCTIONS:
   ....:     print(f"- {coord_type}")
   ....: 
- atmosphere_sigma_coordinate
- ocean_sigma_coordinate
- ocean_s_coordinate
- ocean_s_coordinate_g1
- ocean_s_coordinate_g2

These include:

  • atmosphere_sigma_coordinate: For atmospheric models

  • ocean_sigma_coordinate: Basic ocean sigma coordinates

  • ocean_s_coordinate: Generalized ocean s-coordinates

  • ocean_s_coordinate_g1: Song & Haidvogel s-coordinates (generic form 1)

  • ocean_s_coordinate_g2: Shchepetkin s-coordinates (generic form 2)

Converting sigma to physical coordinates#

The main functionality is to decode sigma coordinates into physical depths (for ocean) or pressures (for atmosphere) using the appropriate formula terms.

Note

The core sigma conversion functions use Numba’s guvectorize decorator to create optimized universal functions (ufuncs). This provides:

  • Efficient computation through JIT compilation

  • Automatic broadcasting over horizontal dimensions

  • Compatibility with dask arrays for lazy evaluation

  • Seamless integration with xarray.apply_ufunc()

Ocean sigma coordinates#

For ocean models, sigma coordinates are converted to depths using formulas that depend on:

  • sigma (σ): The sigma level (typically from -1 to 0)

  • eta (η): Sea surface height

  • depth: Bathymetry (bottom depth)

  • Additional parameters depending on the coordinate type (e.g., hc, theta_s, theta_b for s-coordinates)

# Example with CROCO model output
In [22]: import xoa

# Open a sample with sigma coordinates
In [23]: ds = xoa.open_data_sample("MODELS/CROCO/SOUTH-AFRICA/croco.south-africa.meridional.nc")

# Check the sigma coordinate
In [24]: print(ds.s_rho.attrs.get('standard_name'))
ocean_s_coordinate_g2

# Decode sigma to depths using the sigma module
# The decode_cf_sigma function will automatically detect
# the coordinate type and apply the correct formula
In [25]: from xoa.sigma import decode_sigma

In [26]: dsz = decode_sigma(ds)

In [27]: print(dsz.depth)
<xarray.DataArray 'depth' (time: 1, s_rho: 32, eta_rho: 56, xi_rho: 1)> Size: 14kB
array([[[[-4.47877110e+03],
         [-4.45341971e+03],
         [-4.43332704e+03],
         ...,
         [-7.38942221e+01],
         [-7.34145067e+01],
         [-7.34091375e+01]],

        [[-4.10998820e+03],
         [-4.08677643e+03],
         [-4.06837911e+03],
         ...,
         [-7.04949884e+01],
         [-7.00423304e+01],
         [-7.00372638e+01]],

        [[-3.70399833e+03],
         [-3.68314549e+03],
         [-3.66661727e+03],
         ...,
...
         ...,
         [-4.29529272e+00],
         [-4.27484202e+00],
         [-4.27461272e+00]],

        [[-1.00794918e+01],
         [-1.00480588e+01],
         [-1.00177861e+01],
         ...,
         [-2.57388434e+00],
         [-2.56165088e+00],
         [-2.56151371e+00]],

        [[-3.17876268e+00],
         [-3.15406668e+00],
         [-3.12913607e+00],
         ...,
         [-8.56901999e-01],
         [-8.52836074e-01],
         [-8.52790485e-01]]]], shape=(1, 32, 56, 1))
Coordinates:
  * time     (time) float64 8B 2.592e+06
  * s_rho    (s_rho) float32 128B -0.9844 -0.9531 -0.9219 ... -0.04688 -0.01562
  * eta_rho  (eta_rho) float32 224B 6.0 7.0 8.0 9.0 10.0 ... 58.0 59.0 60.0 61.0
  * xi_rho   (xi_rho) float32 4B 131.0
    lat_rho  (eta_rho, xi_rho) float32 224B -37.67 -37.61 ... -34.03 -33.96
    lon_rho  (eta_rho, xi_rho) float32 224B 18.83 18.83 18.83 ... 18.83 18.83
    depth    (time, s_rho, eta_rho, xi_rho) float64 14kB -4.479e+03 ... -0.8528
Attributes:
    long_name:      Ocean s roms coordinate at rho point
    Vtransform:     2
    standard_name:  ocean_layer_depth
    units:          m

Atmosphere sigma coordinates#

For atmospheric models, sigma coordinates are converted to pressure levels:

# p(k) = ptop + sigma(k) * (ps - ptop)
# where:
# - p: pressure at level k
# - ptop: pressure at top of model
# - sigma: sigma coordinate value
# - ps: surface pressure

Working with vertical sections#

Once sigma coordinates are decoded to physical depths, you can work with vertical sections and perform operations like:

  • Extracting data at specific depth levels

  • Interpolating to regular depth levels using xoa.regrid

  • Computing vertical derivatives

  • Creating vertical transects

See the gallery example Interpolate a meridional section of CROCO outputs to regular depths for a practical demonstration of working with sigma coordinates and vertical sections.

Formula terms mapping#

The module provides automatic mapping from formula terms (as defined in the CF conventions) to the standard xoa CF names:

In [28]: from xoa.sigma import FORMULA_TERMS_TO_CF_NAMES

In [29]: from pprint import pprint

In [30]: pprint(FORMULA_TERMS_TO_CF_NAMES)
{'a': 'thetas',
 'b': 'thetab',
 'c': 'cs',
 'depth': 'bathy',
 'depth_c': 'hc',
 'eta': 'ssh',
 's': 'sig',
 'sigma': 'sig'}

This allows the module to find the required variables in your dataset regardless of their actual names, as long as they follow CF conventions.

Example: Working with CROCO model output#

Here’s a complete example showing how to work with grid operations, regridding, and sigma coordinates using CROCO model output:

In [31]: import xoa

In [32]: import xoa.meta

In [33]: import xoa.regrid

In [34]: import numpy as np

# To have 'xoa' and 'decode_sigma' accessors
In [35]: xoa.register_accessors()

# Load CROCO sample data
In [36]: ds = xoa.open_data_sample("MODELS/CROCO/SOUTH-AFRICA/croco.south-africa.meridional.nc")

# Set internal croco decoding rules and decode s-coordinate and dataset names
In [37]: xoa.meta.set_meta_specs('croco')
Out[37]: <xoa.meta.set_meta_specs at 0x759688e21e50>

In [38]: ds_decoded = ds.decode_sigma().xoa.decode()

# The sigma coordinate has been converted to depths
# You can now regrid to regular depth levels if needed
# Define regular depth levels
In [39]: regular_depths = xr.DataArray(np.linspace(-100, 0, 20), dims='depth', attrs={'units': 'm', 'positive': 'up'})

# Regrid temperature to regular depths
In [40]: temp_regular = xoa.regrid.regrid1d(ds_decoded.ptemp, regular_depths, extrap="top")

In [41]: print(temp_regular[0, :, 0, 0])
<xarray.DataArray 'ptemp' (depth: 20)> Size: 80B
array([15.634643, 15.841397, 16.087772, 16.342157, 16.596542, 16.90044 ,
       17.249428, 17.598415, 17.978926, 18.375126, 18.781078, 19.230177,
       19.676584, 19.866636, 20.036734, 20.02078 , 20.003649, 19.984797,
       19.963604, 19.955126], dtype=float32)
Coordinates:
  * depth    (depth) float64 160B -100.0 -94.74 -89.47 ... -10.53 -5.263 0.0
    y_rho    float32 4B 6.0
    lat_rho  float32 4B -37.67
    lon_rho  float32 4B 18.83
    time     float64 8B 2.592e+06
    x_rho    float32 4B 131.0
Attributes:
    long_name:      Potential temperature
    units:          Celsius
    field:          temperature, scalar, series
    standard_name:  sea_water_potential_temperature

See also#

Contents