Core Concepts

This page explains the key concepts underlying Nereus to help you use it effectively.

Unstructured vs. Regular Grids

Atmospheric and ocean models increasingly use unstructured grids (also called triangular meshes or irregular grids) instead of regular latitude-longitude grids. Examples include:

  • FESOM2: Finite Element Sea ice-Ocean Model

  • ICON: Icosahedral Nonhydrostatic model

  • MPAS: Model for Prediction Across Scales

Unstructured grids offer advantages:

  • Variable resolution (high resolution where needed)

  • Better representation of coastlines and topography

  • No convergence of grid lines at poles

However, they require special handling for visualization and analysis since most plotting and analysis tools expect regular grids.

Note

Unstructured grids use triangular elements that can vary in size, allowing higher resolution in areas of interest (like coastlines) while using larger elements in the open ocean to reduce computational cost.

The Regridding Process

Nereus converts unstructured data to regular grids using spatial interpolation. Two methods are available:

Nearest neighbor (default, method="nearest")

For each target grid point, find the nearest source point using a KD-tree built on 3D Cartesian coordinates. Fast, preserves original values.

Linear (method="linear")

Build a Delaunay triangulation of source points and interpolate using barycentric coordinates. Slower, but produces smooth results without the blocky patterns of nearest neighbor.

The process involves:

  1. Coordinate transformation: Convert lon/lat to 3D Cartesian coordinates on a unit sphere

  2. KD-tree construction: Build a spatial index of source points

  3. Interpolation: Nearest-neighbor lookup or Delaunay-based linear interpolation

  4. Distance validation: Mask points that are too far from any source data

Source longitude conventions (0-360 or -180-180) are handled automatically.

The RegridInterpolator

The RegridInterpolator is the core class for regridding:

interpolator = nr.RegridInterpolator(
    source_lon, source_lat,
    resolution=1.0,           # Target resolution in degrees
    method="nearest",         # "nearest" or "linear"
    influence_radius=80000.0, # Max distance in meters
)

# Apply to any data on the same source grid
regridded = interpolator(data)

Key parameters:

resolution

Can be a single number (e.g., 1.0 for 1-degree) or a tuple (nlon, nlat) for explicit grid dimensions.

method

"nearest" for fast nearest-neighbor lookup, "linear" for smooth Delaunay-based interpolation.

influence_radius

Points on the target grid farther than this distance (in meters) from any source point are masked. This prevents extrapolation into data-void regions.

lon_bounds / lat_bounds

Customize the target grid extent. Defaults to global coverage (-180 to 180, -90 to 90).

Caching and Performance

Building the interpolation weights (KD-tree and nearest neighbor lookup) is the expensive part. Nereus provides automatic caching:

# First call: builds interpolator (~1-2 seconds for large meshes)
fig, ax, interp = nr.plot(data1, lon, lat)

# Second call: reuses cached interpolator (instant)
fig, ax, interp = nr.plot(data2, lon, lat)

The cache is keyed by:

  • Source coordinates (lon, lat arrays)

  • Resolution

  • Influence radius

  • Bounds

You can control caching behavior:

# Configure cache size
nr.set_cache_options(max_memory_items=20)

# Disable caching for a specific plot
nr.plot(data, lon, lat, use_cache=False)

# Manually reuse an interpolator
nr.plot(data, lon, lat, interpolator=interp)

Flexible Input Handling

Nereus is designed to work with data in various formats without requiring manual preprocessing.

Automatic Array Reshaping

Functions like plot(), regrid(), and transect() accept data and coordinates in multiple formats:

  • 1D arrays: Traditional unstructured mesh format

  • 2D arrays: Regular grid format (automatically raveled)

  • Mixed formats: 2D data with 1D coordinates (meshgrid created automatically)

# All of these work:
nr.plot(data_1d, lon_1d, lat_1d)           # Unstructured mesh
nr.plot(data_2d, lon_2d, lat_2d)           # Curvilinear grid
nr.plot(data_2d, lon_1d, lat_1d)           # Regular grid

Warnings are issued when transformations are applied, so you know exactly what’s happening.

Automatic Coordinate Extraction

When using xarray DataArrays, coordinates can be extracted automatically:

import xarray as xr

ds = xr.open_dataset("model_output.nc")
temp = ds.temperature.isel(time=0)

# Coordinates extracted from xarray metadata
nr.plot(temp)
nr.regrid(temp, resolution=0.5)

Nereus recognizes common coordinate names used by various models:

  • FESOM/ICON: lon, lat, nod2d_lon, nod2d_lat

  • NEMO: nav_lon, nav_lat

  • MOM/GFDL: glon, glat, xt_ocean, yt_ocean

  • Generic: longitude, latitude, x, y

Coordinate Systems

Nereus uses standard geographic conventions:

Longitude

-180 to 180 degrees (or 0 to 360, automatically handled)

Latitude

-90 to 90 degrees

Depth

Positive downward in meters (following oceanographic convention)

Area

In square meters (m²)

Internally, Nereus converts lon/lat to Cartesian coordinates on a unit sphere for spatial operations:

\[\begin{split}x &= \cos(\text{lat}) \cdot \cos(\text{lon}) \\ y &= \cos(\text{lat}) \cdot \sin(\text{lon}) \\ z &= \sin(\text{lat})\end{split}\]

This avoids issues with the pole singularity and dateline crossings.

Area-Weighted Operations

Many diagnostics require proper area weighting. For unstructured meshes, each grid point has an associated area (often called “cluster area” or “dual cell area”).

# Total sea ice area
ice_area = sum(concentration * cell_area)

# Area-weighted mean
mean_temp = sum(temp * area) / sum(area)

# Volume-weighted mean (3D)
mean_temp = sum(temp * area * thickness) / sum(area * thickness)

Nereus handles these calculations correctly:

# Sea ice area
total_ice = nr.ice_area(concentration, area)

# Volume mean for upper ocean
upper_mean = nr.volume_mean(
    temp, area, thickness, depth,
    depth_max=500  # Upper 500m only
)

Model-Specific Support

Each model has its own mesh format and conventions. Nereus provides model-specific submodules:

# FESOM2
mesh = nr.fesom.load_mesh("/path/to/mesh/")

# Future: ICON-Ocean
mesh = nr.icono.load_mesh("icon_mesh.nc")

# Future: ICON-Atmosphere
mesh = nr.icona.load_mesh("icon_atmo_mesh.nc")

Each mesh object provides a consistent interface:

  • mesh.lon, mesh.lat: Coordinates

  • mesh.area: Cell areas

  • mesh.find_nearest(): Spatial queries

  • mesh.subset_by_bbox(): Geographic subsetting

The MeshProtocol

All mesh classes follow the MeshBase protocol:

class MeshProtocol:
    @property
    def lon(self) -> NDArray: ...

    @property
    def lat(self) -> NDArray: ...

    @property
    def area(self) -> NDArray: ...

    def find_nearest(self, lon, lat, k=1): ...
    def subset_by_bbox(self, lon_min, lon_max, lat_min, lat_max): ...

This allows writing generic code that works with any supported model.

Dask Integration

Nereus is designed to work with Dask arrays for out-of-core and parallel computation:

import dask.array as da

# Load data lazily with Dask
ds = xr.open_mfdataset("output_*.nc", chunks={"time": 10})

# Diagnostics work with Dask arrays
ice_area_timeseries = nr.ice_area(ds.a_ice, mesh.area)
# Returns a Dask array - computation is deferred

# Trigger computation
result = ice_area_timeseries.compute()

When working with Dask:

  • Regridding and plotting trigger immediate computation (needed for visualization)

  • Diagnostics preserve lazy evaluation where possible

  • Use .compute() or .values to trigger computation

Physical Constants

Nereus uses standard physical constants for ocean calculations:

Constant

Value

Description

Earth radius

6,371,000 m

WGS84 mean radius

Seawater density

1,025 kg/m³

Reference density

Specific heat

3,985 J/(kg·K)

Seawater heat capacity

These are used in calculations like ocean heat content:

\[OHC = \rho \cdot c_p \cdot \sum (T - T_{ref}) \cdot \Delta z \cdot A\]

Where: - \(\rho\) = seawater density - \(c_p\) = specific heat capacity - \(T\) = temperature - \(T_{ref}\) = reference temperature - \(\Delta z\) = layer thickness - \(A\) = cell area