Top-Level API

The top-level nereus module provides the main user-facing API for quick data exploration and analysis.

import nereus as nr

Mesh Functions

Universal Mesh Loading

nereus.load_mesh(path, *, mesh_type=None, use_dask=None, **kwargs)[source]

Universal mesh loader with auto-detection.

Parameters:

  • path (str, Path, or int) – Path to mesh directory/file, or number of points for HEALPix.

  • mesh_type (str, optional) – Mesh type: “fesom”, “healpix”, “nemo”, “ifs_tco”, or “auto”. If None or “auto”, attempts to auto-detect.

  • use_dask (bool, optional) – Whether to use dask arrays. Auto-detects if None.

  • **kwargs – Additional arguments passed to model-specific loader.

Returns:

Standardized mesh dataset.

Return type:

xr.Dataset

Examples

>>> # Auto-detect FESOM mesh
>>> mesh = nr.load_mesh("/path/to/mesh")

>>> # Explicit type
>>> mesh = nr.load_mesh("/path/to/mesh_mask.nc", mesh_type="nemo")

>>> # HEALPix from npoints
>>> mesh = nr.load_mesh(3145728, mesh_type="healpix")

Mesh Creation

nereus.create_lonlat_mesh(resolution, *, lon_bounds=(-180, 180), lat_bounds=(-90, 90), use_dask=None)[source]

Create regular lon-lat mesh.

Parameters:

  • resolution (float or tuple) – Grid resolution in degrees. If tuple, (dlon, dlat).

  • lon_bounds (tuple) – Longitude bounds (min, max).

  • lat_bounds (tuple) – Latitude bounds (min, max).

  • use_dask (bool, optional) – Whether to use dask arrays. Auto-detects if None.

Returns:

Mesh dataset with flattened lon, lat, area.

Return type:

xr.Dataset

nereus.mesh_from_arrays(lon, lat, *, area=None, use_dask=None)[source]

Create mesh from existing coordinate arrays.

Handles 1D, 2D, and regular grid side coordinates. If lon and lat are both 1D but have different sizes, they are treated as regular grid side coordinates and expanded via meshgrid before flattening. 2D arrays are flattened directly.

Parameters:

  • lon (array_like) – Longitude coordinates in degrees.

  • lat (array_like) – Latitude coordinates in degrees.

  • area (array_like, optional) – Cell areas in m^2. If None, estimates from grid spacing.

  • use_dask (bool, optional) – Whether to use dask arrays. Auto-detects if None.

Returns:

Mesh dataset with lon, lat, area.

Return type:

xr.Dataset

Spatial Queries

nereus.find_nearest(lon, lat, query_lon, query_lat, k=1, *, return_distance=False)[source]

Find nearest mesh points to query locations.

Uses a KDTree built on Cartesian coordinates for efficient spherical nearest-neighbor search.

Parameters:

  • lon (array_like) – Longitude coordinates of mesh points in degrees.

  • lat (array_like) – Latitude coordinates of mesh points in degrees.

  • query_lon (float or array_like) – Query longitude(s) in degrees.

  • query_lat (float or array_like) – Query latitude(s) in degrees.

  • k (int) – Number of nearest neighbors to find.

  • return_distance (bool) – If True, also return distances in meters.

Returns:

  • indices (ndarray) – Indices of nearest mesh points. Shape depends on inputs: - Scalar query, k=1: scalar int - Scalar query, k>1: (k,) array - Array query, k=1: (n_queries,) array - Array query, k>1: (n_queries, k) array

  • distances (ndarray, optional) – Distances in meters. Returned only if return_distance=True. Same shape as indices.

Return type:

ndarray[tuple[Any, …], dtype[int64]] | tuple[ndarray[tuple[Any, …], dtype[int64]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> mesh = nr.fesom.load_mesh(path)
>>> idx = nr.find_nearest(mesh["lon"].values, mesh["lat"].values, -30.5, 60.2)
>>> print(f"Nearest point: ({mesh['lon'].values[idx]}, {mesh['lat'].values[idx]})")

>>> # Find 3 nearest neighbors with distances
>>> indices, distances = nr.find_nearest(
...     mesh["lon"].values, mesh["lat"].values,
...     [-30, -31], [60, 61],
...     k=3,
...     return_distance=True,
... )
nereus.subset_by_bbox(lon, lat, lon_min, lon_max, lat_min, lat_max)[source]

Get mask for points within bounding box.

Parameters:

  • lon (array_like) – Longitude coordinates in degrees.

  • lat (array_like) – Latitude coordinates in degrees.

  • lon_min (float) – Longitude bounds in degrees.

  • lon_max (float) – Longitude bounds in degrees.

  • lat_min (float) – Latitude bounds in degrees.

  • lat_max (float) – Latitude bounds in degrees.

Returns:

mask – Boolean mask of points within bounds.

Return type:

ndarray

Examples

>>> mesh = nr.fesom.load_mesh(path)
>>> mask = nr.subset_by_bbox(
...     mesh["lon"].values, mesh["lat"].values,
...     lon_min=-10, lon_max=10,
...     lat_min=-5, lat_max=5,
... )
>>> equatorial_area = mesh["area"].values[mask].sum()
nereus.points_in_polygon(lon, lat, polygon_lon, polygon_lat)[source]

Get mask for points inside polygon.

Uses matplotlib.path for point-in-polygon testing.

Parameters:

  • lon (array_like) – Longitude coordinates of mesh points in degrees.

  • lat (array_like) – Latitude coordinates of mesh points in degrees.

  • polygon_lon (array_like) – Longitude coordinates of polygon vertices.

  • polygon_lat (array_like) – Latitude coordinates of polygon vertices.

Returns:

mask – Boolean mask of points inside polygon.

Return type:

ndarray

Notes

The polygon is assumed to be defined in Cartesian lon-lat space (not geodesic). For large-scale polygons crossing the dateline, consider normalizing longitudes first.

Examples

>>> # Select points in a triangular region
>>> poly_lon = [-10, 10, 0, -10]
>>> poly_lat = [0, 0, 10, 0]
>>> mask = nr.points_in_polygon(lon, lat, poly_lon, poly_lat)

Plotting Functions

nereus.plot(data, lon=None, lat=None, *, projection='pc', extent=None, resolution=1.0, interpolator=None, method='nearest', influence_radius=80000.0, cmap='viridis', vmin=None, vmax=None, coastlines=True, land=False, gridlines=False, colorbar=True, colorbar_label=None, title=None, figsize=None, ax=None, use_cache=True, **kwargs)[source]

Plot 2D map of unstructured data.

This function regrids unstructured data to a regular grid and plots it on a map with the specified projection.

The function accepts various input formats and automatically transforms them to 1D arrays for plotting:

  • All 1D arrays of same size: used directly (no warning)

  • 2D data with 2D lon/lat (same shape): all raveled to 1D

  • 1D data with 2D lon/lat: lon/lat raveled to match data

  • 2D data with 1D lon/lat: meshgrid created, then all raveled

A warning is issued whenever array transformations are applied.

If lon/lat are not provided and data is an xarray DataArray, the function will attempt to extract coordinates automatically by looking for common coordinate names (lon/lat, longitude/latitude, x/y, etc.).

Parameters:

  • data (array_like) – Data values at unstructured points. Can be 1D or 2D array. If xarray DataArray, coordinates may be extracted automatically.

  • lon (array_like, optional) – Longitude coordinates. Can be 1D or 2D array. If None, will attempt to extract from data (xarray only).

  • lat (array_like, optional) – Latitude coordinates. Can be 1D or 2D array. If None, will attempt to extract from data (xarray only).

  • projection (str or Projection) – Map projection. Options: “pc”, “rob”, “merc”, “npstere”, “spstere”, “moll”, “ortho”, “lcc”, or a Cartopy Projection.

  • extent (tuple of float, optional) – Map extent (lon_min, lon_max, lat_min, lat_max).

  • resolution (float or tuple of int) – Grid resolution for regridding.

  • interpolator (RegridInterpolator, optional) – Pre-computed interpolator. If None, one will be created.

  • method ({"nearest", "idw", "linear", "cubic"}) – Interpolation method. “nearest” uses nearest-neighbor lookup (fast). “idw” uses inverse distance weighting (fast, smooth). “linear” uses Delaunay triangulation with barycentric interpolation. “cubic” uses Clough-Tocher C1 interpolation (smoothest). Default is “nearest”.

  • influence_radius (float) – Maximum influence radius in meters for interpolation. Default is 80 km.

  • cmap (str) – Colormap name.

  • vmin (float, optional) – Color scale limits.

  • vmax (float, optional) – Color scale limits.

  • coastlines (bool) – Whether to draw coastlines.

  • land (bool) – Whether to fill land areas.

  • gridlines (bool) – Whether to draw gridlines.

  • colorbar (bool) – Whether to add a colorbar.

  • colorbar_label (str, optional) – Label for the colorbar.

  • title (str, optional) – Plot title.

  • figsize (tuple of float, optional) – Figure size (width, height) in inches.

  • ax (Axes, optional) – Existing axes to plot on. If None, creates new figure.

  • use_cache (bool) – Whether to use the interpolator cache.

  • **kwargs (Any) – Additional arguments passed to pcolormesh.

Returns:

  • fig (Figure) – The matplotlib Figure.

  • ax (Axes) – The matplotlib Axes (GeoAxes).

  • interpolator (RegridInterpolator) – The interpolator used (can be reused).

Return type:

tuple[‘Figure’, ‘Axes’, RegridInterpolator]

Examples

>>> fig, ax, interp = nr.plot(temp, mesh.lon, mesh.lat)
>>> fig, ax, _ = nr.plot(salinity, mesh.lon, mesh.lat, interpolator=interp)
nereus.transect(data, lon=None, lat=None, depth=None, start=None, end=None, *, n_points=100, cmap='viridis', vmin=None, vmax=None, depth_lim=None, invert_depth=True, colorbar=True, colorbar_label=None, title=None, figsize=None, ax=None, **kwargs)[source]

Plot vertical transect along a great circle path.

The function accepts various coordinate formats and automatically transforms them to 1D arrays:

  • Both 1D with same size: used directly (no warning)

  • Both 2D with same shape: raveled to 1D

  • Both 1D with different sizes: meshgrid created, then raveled

A warning is issued whenever coordinate transformations are applied.

If lon/lat are not provided and data is an xarray DataArray, the function will attempt to extract coordinates automatically by looking for common coordinate names (lon/lat, longitude/latitude, x/y, etc.).

Parameters:

  • data (array_like) – Data values with shape (nlevels, npoints) for 2D, or (nlevels, nlat, nlon) for 3D regular grids. 3D data is automatically reshaped to 2D. If xarray DataArray, coordinates may be extracted automatically.

  • lon (array_like, optional) – Longitude coordinates. Can be 1D or 2D array. If None, will attempt to extract from data (xarray only).

  • lat (array_like, optional) – Latitude coordinates. Can be 1D or 2D array. If None, will attempt to extract from data (xarray only).

  • depth (array_like) – 1D array of depth levels (positive downward).

  • start (tuple of float) – Start point (lon, lat).

  • end (tuple of float) – End point (lon, lat).

  • n_points (int) – Number of points along the transect.

  • cmap (str) – Colormap name.

  • vmin (float, optional) – Color scale limits.

  • vmax (float, optional) – Color scale limits.

  • depth_lim (tuple of float, optional) – Depth/height limits (min, max). If None, uses data range.

  • invert_depth (bool) – Whether to invert vertical axis. Default True for ocean (0 at top, depth increases downward). Set False for atmosphere (height increases upward).

  • colorbar (bool) – Whether to add a colorbar.

  • colorbar_label (str, optional) – Label for the colorbar.

  • title (str, optional) – Plot title.

  • figsize (tuple of float, optional) – Figure size.

  • ax (Axes, optional) – Existing axes to plot on.

  • **kwargs (Any) – Additional arguments passed to pcolormesh.

Returns:

  • fig (Figure) – The matplotlib Figure.

  • ax (Axes) – The matplotlib Axes.

Return type:

tuple[‘Figure’, ‘Axes’]

Examples

>>> fig, ax = nr.transect(
...     temp, mesh.lon, mesh.lat, depth,
...     start=(-30, 60), end=(30, 60)
... )

Regridding Functions

nereus.regrid(data, lon=None, lat=None, resolution=1.0, method='nearest', influence_radius=80000.0, fill_value=nan, lon_bounds=(-180.0, 180.0), lat_bounds=(-90.0, 90.0), as_xarray=False)[source]

Regrid unstructured data to regular grid.

This is a convenience function that creates a RegridInterpolator and applies it. For repeated regridding with the same source grid, create a RegridInterpolator once and reuse it.

Supports multi-dimensional data where the last axis contains the spatial points. For example:

  • 1D data (npoints,): single field

  • 2D data (nlevels, npoints): multi-level unstructured data (e.g., FESOM, ICON)

  • ND data (*dims, npoints): arbitrary leading dimensions

Coordinate arrays can be:

  • 1D arrays of same size: unstructured mesh coordinates (used directly)

  • 1D arrays of different sizes: regular grid side coordinates (meshgrid created)

  • 2D arrays of same shape: full coordinate arrays (raveled to 1D)

A warning is issued whenever coordinate transformations are applied.

If lon/lat are not provided and data is an xarray DataArray, the function will attempt to extract coordinates automatically by looking for common coordinate names (lon/lat, longitude/latitude, x/y, etc.).

Parameters:

  • data (array_like) – Data to interpolate. Last axis must be npoints (matching coordinates). Can be 1D (npoints,), 2D (nlevels, npoints), or ND (*dims, npoints). If xarray DataArray, coordinates may be extracted automatically.

  • lon (array_like, optional) – Source grid longitude coordinates. Can be 1D or 2D array. If None, will attempt to extract from data (xarray only).

  • lat (array_like, optional) – Source grid latitude coordinates. Can be 1D or 2D array. If None, will attempt to extract from data (xarray only).

  • resolution (float or tuple of int) – Target grid resolution.

  • method ({"nearest", "idw", "linear", "cubic"}) – Interpolation method. “nearest” uses nearest-neighbor lookup. “idw” uses inverse distance weighting (fast, smooth). “linear” uses Delaunay triangulation with barycentric interpolation. “cubic” uses Clough-Tocher C1 interpolation (smoothest).

  • influence_radius (float) – Maximum influence radius in meters.

  • fill_value (float) – Value for invalid points.

  • lon_bounds (tuple of float) – Target grid longitude bounds.

  • lat_bounds (tuple of float) – Target grid latitude bounds.

  • as_xarray (bool, default False) – If True, wrap the regridded array in an xr.DataArray with lat and lon as 1-D dimension coordinates. Leading dimensions (e.g. time, depth) and their coordinates are preserved when the input is an xr.DataArray. The return type of the tuple’s first element changes from NDArray to xr.DataArray.

Returns:

  • regridded (ndarray or xr.DataArray) – Regridded data. Returns xr.DataArray when as_xarray=True, otherwise ndarray.

  • interpolator (RegridInterpolator) – The interpolator used (can be reused for other variables).

Return type:

tuple[NDArray[np.floating], RegridInterpolator]

class nereus.RegridInterpolator(source_lon, source_lat, resolution=1.0, method='nearest', influence_radius=80000.0, lon_bounds=(-180.0, 180.0), lat_bounds=(-90.0, 90.0))[source]

Bases: object

Pre-computed interpolation for fast repeated regridding.

This class computes and stores interpolation weights for regridding unstructured data to a regular grid. The computation is done once during initialization, allowing fast repeated application.

Parameters:

  • source_lon (array_like) – Source grid longitude coordinates in degrees.

  • source_lat (array_like) – Source grid latitude coordinates in degrees.

  • resolution (float or tuple of int) – Target grid resolution. If float, specifies degrees per cell. If tuple (nlon, nlat), specifies number of grid points.

  • method ({"nearest", "idw", "linear", "cubic"}) – Interpolation method. “nearest” uses nearest-neighbor lookup via KDTree (fast). “idw” uses inverse distance weighting with 8 nearest neighbors (fast, smooth). “linear” uses Delaunay triangulation with barycentric interpolation (slower but smoother). “cubic” uses Clough-Tocher C1 interpolation for the smoothest results. Source longitudes are automatically normalized to match the target grid’s lon_bounds so that any input convention (0-360 or -180-180) works transparently.

  • influence_radius (float) – Maximum influence radius in meters. Points beyond this distance from any source point are masked. Default is 80 km.

  • lon_bounds (tuple of float) – Target grid longitude bounds. Default is (-180, 180).

  • lat_bounds (tuple of float) – Target grid latitude bounds. Default is (-90, 90).

target_lon

Target grid longitude coordinates (2D).

Type:

ndarray

target_lat

Target grid latitude coordinates (2D).

Type:

ndarray

indices

Source indices for each target point.

Type:

ndarray

distances

Distances from target to source points (in chord units).

Type:

ndarray

valid_mask

Boolean mask of valid target points within influence radius.

Type:

ndarray

Examples

>>> interpolator = RegridInterpolator(mesh_lon, mesh_lat, resolution=1.0)
>>> regridded = interpolator(data)
>>> regridded.shape
(180, 360)

Use linear interpolation for smoother results:

>>> interpolator = RegridInterpolator(
...     mesh_lon, mesh_lat, resolution=1.0, method="linear"
... )
source_lon: ndarray[tuple[Any, ...], dtype[floating]]
source_lat: ndarray[tuple[Any, ...], dtype[floating]]
resolution: float | tuple[int, int] = 1.0
method: Literal['nearest', 'idw', 'linear', 'cubic'] = 'nearest'
influence_radius: float = 80000.0
lon_bounds: tuple[float, float] = (-180.0, 180.0)
lat_bounds: tuple[float, float] = (-90.0, 90.0)
target_lon: ndarray[tuple[Any, ...], dtype[floating]]
target_lat: ndarray[tuple[Any, ...], dtype[floating]]
indices: ndarray[tuple[Any, ...], dtype[int64]]
distances: ndarray[tuple[Any, ...], dtype[floating]]
valid_mask: ndarray[tuple[Any, ...], dtype[bool]]
__post_init__()[source]

Initialize interpolation weights.

__call__(data, fill_value=nan)[source]

Apply interpolation to data.

Parameters:

  • data (array_like) – Data to interpolate. Can be: - 1D array of shape (npoints,) - 2D array of shape (nlevels, npoints) or (ntime, npoints) - ND array with last axis = npoints

  • fill_value (float) – Value for invalid points outside influence radius.

Returns:

Regridded data. Shape depends on input: - 1D input: (nlat, nlon) - 2D input: (extra_dim, nlat, nlon) - ND input: (*leading_dims, nlat, nlon)

Return type:

ndarray

property shape: tuple[int, int]

Shape of target grid (nlat, nlon).

__init__(source_lon, source_lat, resolution=1.0, method='nearest', influence_radius=80000.0, lon_bounds=(-180.0, 180.0), lat_bounds=(-90.0, 90.0))

Sea Ice Diagnostics

nereus.ice_area(concentration, area, *, mask=None, as_xarray=False)[source]

Compute total sea ice area.

Sea ice area is the sum of grid cell areas weighted by ice concentration.

This function is dask-friendly: if inputs are dask arrays, the result will be a lazy dask array that can be computed later with .compute().

Parameters:

  • concentration (array_like) – Sea ice concentration (fraction, 0-1). Can be 1D (npoints,) or ND with the last axis being npoints.

  • area (array_like) – Grid cell areas in m^2.

  • mask (array_like, optional) – Boolean mask (True = include). If None, all points are included.

  • as_xarray (bool) – If True, return the result as an xarray DataArray with dimension names and coordinates preserved from the input (default False).

Returns:

Total sea ice area in m^2. Returns float for 1D numpy input, ndarray for ND numpy input, dask array if inputs are dask, or xr.DataArray if as_xarray=True.

Return type:

float or ndarray or dask.array or xr.DataArray

Examples

>>> # 1D concentration array
>>> total_area = nr.ice_area(sic, mesh.area)

>>> # With time dimension (time, npoints)
>>> area_timeseries = nr.ice_area(sic, mesh.area)

>>> # With hemisphere mask
>>> nh_area = nr.ice_area(sic, mesh.area, mask=mesh.lat > 0)

>>> # With dask arrays (lazy computation)
>>> area = nr.ice_area(sic_dask, mesh.area)
>>> area.compute()  # triggers actual computation
nereus.ice_area_nh(concentration, area, lat, *, as_xarray=False)[source]

Compute Northern Hemisphere sea ice area.

Convenience function that calls ice_area with a Northern Hemisphere mask.

Parameters:

  • concentration (array_like) – Sea ice concentration (fraction, 0-1).

  • area (array_like) – Grid cell areas in m^2.

  • lat (array_like) – Latitude of grid points in degrees.

  • as_xarray (bool) – If True, return the result as an xarray DataArray (default False).

Returns:

Northern Hemisphere sea ice area in m^2.

Return type:

float or ndarray or dask.array or xr.DataArray

Examples

>>> nh_area = nr.ice_area_nh(sic, mesh.area, mesh.lat)
nereus.ice_area_sh(concentration, area, lat, *, as_xarray=False)[source]

Compute Southern Hemisphere sea ice area.

Convenience function that calls ice_area with a Southern Hemisphere mask.

Parameters:

  • concentration (array_like) – Sea ice concentration (fraction, 0-1).

  • area (array_like) – Grid cell areas in m^2.

  • lat (array_like) – Latitude of grid points in degrees.

  • as_xarray (bool) – If True, return the result as an xarray DataArray (default False).

Returns:

Southern Hemisphere sea ice area in m^2.

Return type:

float or ndarray or dask.array or xr.DataArray

Examples

>>> sh_area = nr.ice_area_sh(sic, mesh.area, mesh.lat)
nereus.ice_volume(thickness, area, concentration=None, *, mask=None, as_xarray=False)[source]

Compute total sea ice volume.

This function handles two types of thickness definitions:

Effective thickness (default, when concentration=None):

Thickness averaged over the entire grid cell, including open water. This is common in model output (e.g., CICE volume-per-area variables). Formula: V = sum(h_eff * cell_area)

Real thickness (when concentration is provided):

Thickness averaged only over ice-covered area (e.g., CMIP6 sithick). Formula: V = sum(h_ice * concentration * cell_area)

This function is dask-friendly: if inputs are dask arrays, the result will be a lazy dask array that can be computed later with .compute().

Parameters:

  • thickness (array_like) –

    Sea ice thickness in meters. Can be 1D (npoints,) or ND with the last axis being npoints.

    • If concentration is None: interpreted as effective thickness (grid-cell mean, already includes open water as zero).

    • If concentration is provided: interpreted as real thickness (ice-area mean, the physical thickness of the ice itself).

  • area (array_like) – Grid cell areas in m^2.

  • concentration (array_like, optional) – Sea ice concentration (fraction, 0-1). Provide this when thickness is “real thickness” (ice-area mean, like CMIP6 sithick). Leave as None when thickness is “effective thickness” (grid-cell mean).

  • mask (array_like, optional) – Boolean mask (True = include).

  • as_xarray (bool) – If True, return the result as an xarray DataArray with dimension names and coordinates preserved from the input (default False).

Returns:

Total sea ice volume in m^3. Returns a dask array if inputs are dask arrays (call .compute() to get the result), or xr.DataArray if as_xarray=True.

Return type:

float or ndarray or dask.array or xr.DataArray

Examples

>>> # Effective thickness (e.g., model output with volume-per-area)
>>> volume = nr.ice_volume(h_eff, mesh.area)

>>> # Real thickness (e.g., CMIP6 sithick) - must provide concentration
>>> volume = nr.ice_volume(sithick, mesh.area, concentration=siconc)

>>> # With dask arrays (lazy computation)
>>> volume = nr.ice_volume(h_eff_dask, mesh.area)
>>> volume.compute()  # triggers actual computation

Notes

Common pitfall: Using real thickness (sithick) without concentration will overestimate volume by a factor of 1/concentration, because it assumes ice covers the entire grid cell.

nereus.ice_volume_nh(thickness, area, lat, concentration=None, *, as_xarray=False)[source]

Compute Northern Hemisphere sea ice volume.

Convenience function that calls ice_volume with a Northern Hemisphere mask.

Parameters:

  • thickness (array_like) – Sea ice thickness in meters.

  • area (array_like) – Grid cell areas in m^2.

  • lat (array_like) – Latitude of grid points in degrees.

  • concentration (array_like, optional) – Sea ice concentration (fraction, 0-1). Required if thickness is “real thickness” (ice-area mean).

  • as_xarray (bool) – If True, return the result as an xarray DataArray (default False).

Returns:

Northern Hemisphere sea ice volume in m^3.

Return type:

float or ndarray or dask.array or xr.DataArray

Examples

>>> nh_volume = nr.ice_volume_nh(sit, mesh.area, mesh.lat)
nereus.ice_volume_sh(thickness, area, lat, concentration=None, *, as_xarray=False)[source]

Compute Southern Hemisphere sea ice volume.

Convenience function that calls ice_volume with a Southern Hemisphere mask.

Parameters:

  • thickness (array_like) – Sea ice thickness in meters.

  • area (array_like) – Grid cell areas in m^2.

  • lat (array_like) – Latitude of grid points in degrees.

  • concentration (array_like, optional) – Sea ice concentration (fraction, 0-1). Required if thickness is “real thickness” (ice-area mean).

  • as_xarray (bool) – If True, return the result as an xarray DataArray (default False).

Returns:

Southern Hemisphere sea ice volume in m^3.

Return type:

float or ndarray or dask.array or xr.DataArray

Examples

>>> sh_volume = nr.ice_volume_sh(sit, mesh.area, mesh.lat)
nereus.ice_extent(concentration, area, *, threshold=0.15, mask=None, as_xarray=False)[source]

Compute sea ice extent.

Sea ice extent is the total area of grid cells where ice concentration exceeds a threshold (typically 15%).

This function is dask-friendly: if inputs are dask arrays, the result will be a lazy dask array that can be computed later with .compute().

Parameters:

  • concentration (array_like) – Sea ice concentration (fraction, 0-1). Can be 1D (npoints,) or ND with the last axis being npoints.

  • area (array_like) – Grid cell areas in m^2.

  • threshold (float) – Concentration threshold (default 0.15 = 15%).

  • mask (array_like, optional) – Boolean mask (True = include).

  • as_xarray (bool) – If True, return the result as an xarray DataArray with dimension names and coordinates preserved from the input (default False).

Returns:

Total sea ice extent in m^2. Returns a dask array if inputs are dask arrays (call .compute() to get the result), or xr.DataArray if as_xarray=True.

Return type:

float or ndarray or dask.array or xr.DataArray

Examples

>>> extent = nr.ice_extent(sic, mesh.area)
>>> extent_nh = nr.ice_extent(sic, mesh.area, mask=mesh.lat > 0)

>>> # With dask arrays (lazy computation)
>>> extent = nr.ice_extent(sic_dask, mesh.area)
>>> extent.compute()  # triggers actual computation
nereus.ice_extent_nh(concentration, area, lat, *, threshold=0.15, as_xarray=False)[source]

Compute Northern Hemisphere sea ice extent.

Convenience function that calls ice_extent with a Northern Hemisphere mask.

Parameters:

  • concentration (array_like) – Sea ice concentration (fraction, 0-1).

  • area (array_like) – Grid cell areas in m^2.

  • lat (array_like) – Latitude of grid points in degrees.

  • threshold (float) – Concentration threshold (default 0.15 = 15%).

  • as_xarray (bool) – If True, return the result as an xarray DataArray (default False).

Returns:

Northern Hemisphere sea ice extent in m^2.

Return type:

float or ndarray or dask.array or xr.DataArray

Examples

>>> nh_extent = nr.ice_extent_nh(sic, mesh.area, mesh.lat)
nereus.ice_extent_sh(concentration, area, lat, *, threshold=0.15, as_xarray=False)[source]

Compute Southern Hemisphere sea ice extent.

Convenience function that calls ice_extent with a Southern Hemisphere mask.

Parameters:

  • concentration (array_like) – Sea ice concentration (fraction, 0-1).

  • area (array_like) – Grid cell areas in m^2.

  • lat (array_like) – Latitude of grid points in degrees.

  • threshold (float) – Concentration threshold (default 0.15 = 15%).

  • as_xarray (bool) – If True, return the result as an xarray DataArray (default False).

Returns:

Southern Hemisphere sea ice extent in m^2.

Return type:

float or ndarray or dask.array or xr.DataArray

Examples

>>> sh_extent = nr.ice_extent_sh(sic, mesh.area, mesh.lat)

Ocean Diagnostics

nereus.surface_mean(data, area, *, mask=None, as_xarray=False)[source]

Compute area-weighted mean of a 2D field (single level).

This is commonly used for surface fields like SST, SSS, or for analyzing a single depth level.

This function is dask-friendly: if inputs are dask arrays, the result will be a lazy dask array that can be computed later with .compute().

Parameters:

  • data (array_like) – 2D data with shape (npoints,) or higher-dimensional with the last axis being npoints. For time series, shape would be (ntime, npoints). Regular-grid data with 2D spatial dimensions (e.g. (nlat, nlon) or (ntime, nlat, nlon)) is automatically flattened when nlat*nlon matches the area size.

  • area (array_like) – Grid cell areas in m^2, shape (npoints,).

  • mask (array_like, optional) – Boolean mask for horizontal points, shape (npoints,). True = include.

  • as_xarray (bool) – If True, return the result as an xarray DataArray with dimension names and coordinates preserved from the input (default False).

Returns:

Area-weighted mean. Returns float for 1D numpy input (npoints,), ndarray for higher-dimensional numpy input, dask array if inputs are dask, or xr.DataArray if as_xarray=True.

Return type:

float or ndarray or dask.array or xr.DataArray

Examples

>>> # Mean SST
>>> mean_sst = nr.surface_mean(sst, mesh.area)

>>> # Mean SST in a region
>>> mean_sst = nr.surface_mean(sst, mesh.area, mask=region_mask)

>>> # With dask arrays (lazy computation)
>>> mean_sst = nr.surface_mean(sst_dask, mesh.area)
>>> mean_sst.compute()  # triggers actual computation
nereus.volume_mean(data, area, thickness, depth=None, *, depth_min=None, depth_max=None, mask=None, as_xarray=False)[source]

Compute volume-weighted mean of a quantity.

This function is dask-friendly: if inputs are dask arrays, the result will be a lazy dask array that can be computed later with .compute().

Parameters:

  • data (array_like) – 3D data with shape (nlevels, npoints) or higher-dimensional with the last two axes being (nlevels, npoints). For time series, shape would be (ntime, nlevels, npoints).

  • area (array_like) – Grid cell areas in m^2. Can be either: - 1D array of shape (npoints,) for surface area (uniform across depth) - 2D array of shape (nlevels, npoints) for depth-dependent area If 2D and has one extra level compared to data layers, the extra level is dropped with a warning (levels vs layers).

  • thickness (array_like) – Layer thicknesses in meters, shape (nlevels, npoints) or (nlevels,) if uniform across points.

  • depth (array_like, optional) – Depth of layer centers in meters (positive downward), shape (nlevels,). Required if depth_min or depth_max are specified.

  • depth_min (float, optional) – Minimum depth to include (meters, positive downward).

  • depth_max (float, optional) – Maximum depth to include (meters, positive downward).

  • mask (array_like, optional) – Boolean mask for horizontal points, shape (npoints,). True = include.

  • as_xarray (bool) – If True, return the result as an xarray DataArray with dimension names and coordinates preserved from the input (default False).

Returns:

Volume-weighted mean. Returns float for 2D numpy input (nlevels, npoints), ndarray for higher-dimensional numpy input, dask array if inputs are dask, or xr.DataArray if as_xarray=True.

Return type:

float or ndarray or dask.array or xr.DataArray

Examples

>>> # Mean temperature in upper 500m
>>> mean_temp = nr.volume_mean(
...     temp, mesh.area, mesh.layer_thickness, mesh.depth,
...     depth_max=500
... )

>>> # Mean salinity over full depth
>>> mean_sal = nr.volume_mean(sal, mesh.area, mesh.layer_thickness)

>>> # With dask arrays (lazy computation)
>>> mean_temp = nr.volume_mean(temp_dask, mesh.area, mesh.layer_thickness)
>>> mean_temp.compute()  # triggers actual computation
nereus.heat_content(temperature, area, thickness, depth=None, *, depth_min=None, depth_max=None, reference_temp=0.0, mask=None, rho=1025.0, cp=3990.0, output='total', as_xarray=False)[source]

Compute ocean heat content.

Heat content can be computed as either: - Total (default): OHC = rho * cp * sum(T * thickness * area) in Joules - Map: OHC = rho * cp * sum_z(T * thickness) in J/m² at each point

This function is dask-friendly: if inputs are dask arrays, the result will be a lazy dask array that can be computed later with .compute().

Parameters:

  • temperature (array_like) – Temperature in degrees Celsius, shape (nlevels, npoints) or higher dimensional with the last two axes being (nlevels, npoints).

  • area (array_like) – Grid cell areas in m^2. Can be either: - 1D array of shape (npoints,) for surface area (uniform across depth) - 2D array of shape (nlevels, npoints) for depth-dependent area If 2D and has one extra level compared to data layers, the extra level is dropped with a warning (levels vs layers). Note: area is not used when output=”map”.

  • thickness (array_like) – Layer thicknesses in meters, shape (nlevels, npoints) or (nlevels,).

  • depth (array_like, optional) – Depth of layer centers in meters (positive downward). Required if depth_min or depth_max are specified.

  • depth_min (float, optional) – Minimum depth to include (meters, positive downward).

  • depth_max (float, optional) – Maximum depth to include (meters, positive downward).

  • reference_temp (float) – Reference temperature for heat content calculation. Default 0°C.

  • mask (array_like, optional) – Boolean mask for horizontal points. True = include.

  • rho (float) – Seawater density in kg/m^3. Default 1025.

  • cp (float) – Specific heat capacity in J/(kg·K). Default 3990.

  • output (str) – Output type. One of: - “total”: Total heat content in Joules (scalar per timestep) - “map”: Heat content per unit area in J/m² (2D field at each point) Default is “total”.

  • as_xarray (bool) – If True, return the result as an xarray DataArray with dimension names and coordinates preserved from the input (default False).

Returns:

If output=”total”: Ocean heat content in Joules. If output=”map”: Heat content per unit area in J/m², shape (npoints,) or (…, npoints) for higher-dimensional input. Returns a dask array if inputs are dask, or xr.DataArray if as_xarray=True.

Return type:

float or ndarray or dask.array or xr.DataArray

Examples

>>> # Total ocean heat content
>>> ohc = nr.heat_content(temp, mesh.area, mesh.layer_thickness)

>>> # Heat content in upper 700m
>>> ohc_700 = nr.heat_content(
...     temp, mesh.area, mesh.layer_thickness, mesh.depth,
...     depth_max=700
... )

>>> # Heat content map (J/m² at each point, like FESOM2 output)
>>> ohc_map = nr.heat_content(
...     temp, mesh.area, mesh.layer_thickness,
...     output="map"
... )

>>> # With dask arrays (lazy computation)
>>> ohc = nr.heat_content(temp_dask, mesh.area, mesh.layer_thickness)
>>> ohc.compute()  # triggers actual computation

Hovmoller Diagrams

nereus.hovmoller(data, area, time=None, depth=None, lat=None, *, mode='depth', lat_bins=None, mask=None, as_xarray=False)[source]

Compute Hovmoller diagram data.

Computes area-weighted means at each time step, binned by either depth level or latitude.

This function is dask-friendly: if inputs are dask arrays, the result will be lazy dask arrays.

If lat is not provided and data is an xarray DataArray, the function will attempt to extract latitude coordinates automatically by looking for common coordinate names (lat, latitude, y, etc.).

Parameters:

  • data (array_like) – Data array. For depth mode: shape (ntime, nlevels, npoints). For latitude mode: shape (ntime, npoints) or (ntime, nlevels, npoints). 4D arrays with shape (ntime, nlevels, nlat, nlon) are automatically flattened to (ntime, nlevels, nlat*nlon) in both modes. If xarray DataArray, lat coordinates may be extracted automatically.

  • area (array_like) – Grid cell areas in m^2. Can be either: - 1D array of shape (npoints,) for surface area (uniform across depth) - 2D array of shape (nlevels, npoints) for depth-dependent area If 2D and has one extra level compared to data layers, the extra level is dropped with a warning (levels vs layers).

  • time (array_like, optional) – Time coordinates. If None, uses integer indices.

  • depth (array_like, optional) – Depth levels in meters. Required for mode=”depth”.

  • lat (array_like, optional) – Latitude coordinates in degrees. Required for mode=”latitude”. If None and mode=”latitude”, will attempt to extract from data (xarray only).

  • mode ({"depth", "latitude"}) – Type of Hovmoller diagram.

  • lat_bins (array_like, optional) – Latitude bin edges for mode=”latitude”. Default is 5-degree bins.

  • mask (array_like, optional) – Boolean mask for horizontal points, shape (npoints,). True = include.

  • as_xarray (bool) – If True, return the result as an xarray DataArray with time and depth/latitude dimension coordinates instead of a 3-tuple (default False).

Returns:

If as_xarray=False (default):

Tuple of (time_out, y_out, data_out) where data_out has shape (ntime, ny).

If as_xarray=True:

xr.DataArray with dims (“time”, “depth”) or (“time”, “latitude”) and corresponding coordinates.

Return type:

tuple or xr.DataArray

Examples

>>> # Time-depth Hovmoller
>>> t, z, hov = nr.hovmoller(temp, mesh.area, time, depth, mode="depth")

>>> # Time-latitude Hovmoller
>>> t, lat, hov = nr.hovmoller(sst, mesh.area, time, lat=mesh.lat, mode="latitude")

>>> # With dask arrays (depth mode)
>>> t, z, hov = nr.hovmoller(temp_dask, mesh.area, time, depth, mode="depth")
>>> hov.compute()  # triggers actual computation
nereus.plot_hovmoller(time, y, data, *, mode='depth', cmap='RdBu_r', vmin=None, vmax=None, colorbar=True, colorbar_label=None, title=None, figsize=None, ax=None, invert_y=None, anomaly=False, y_scale='linear', y_scale_kw=None, **kwargs)[source]

Plot a Hovmoller diagram.

Parameters:

  • time (array_like) – Time coordinates.

  • y (array_like) – Depth or latitude coordinates.

  • data (array_like) – Hovmoller data, shape (ntime, ny).

  • mode ({"depth", "latitude"}) – Type of diagram (affects axis labels and orientation).

  • cmap (str) – Colormap name.

  • vmin (float, optional) – Color scale limits.

  • vmax (float, optional) – Color scale limits.

  • colorbar (bool) – Whether to add a colorbar.

  • colorbar_label (str, optional) – Label for the colorbar.

  • title (str, optional) – Plot title.

  • figsize (tuple of float, optional) – Figure size.

  • ax (Axes, optional) – Existing axes to plot on.

  • invert_y (bool, optional) – Whether to invert y-axis. Default True for depth, False for latitude.

  • anomaly (bool) – If True and mode=”depth”, plot anomaly relative to first time step (data - data[0, :]). Default False.

  • y_scale ({"linear", "sqrt", "power", "symlog"}) – Vertical axis scaling for depth mode. Options: - “linear”: No transformation (default) - “sqrt”: Square root transform, gives more space to surface layers - “power”: Power transform with configurable exponent (see y_scale_kw) - “symlog”: Symmetric log scale, linear near zero then logarithmic

  • y_scale_kw (dict, optional) – Additional parameters for y_scale: - For “power”: {“exponent”: 0.4} (default 0.4, smaller = more surface detail) - For “symlog”: {“linthresh”: 10} (linear threshold in meters, default 10)

  • **kwargs (Any) – Additional arguments passed to pcolormesh.

Returns:

  • fig (Figure) – The matplotlib Figure.

  • ax (Axes) – The matplotlib Axes.

Return type:

tuple[‘Figure’, ‘Axes’]

Examples

>>> # Square root scaling for more surface detail
>>> fig, ax = plot_hovmoller(time, depth, data, y_scale="sqrt")

>>> # Power scaling with custom exponent (smaller = more surface detail)
>>> fig, ax = plot_hovmoller(time, depth, data, y_scale="power",
...                          y_scale_kw={"exponent": 0.3})

>>> # Symmetric log: linear in upper 20m, log below
>>> fig, ax = plot_hovmoller(time, depth, data, y_scale="symlog",
...                          y_scale_kw={"linthresh": 20})

Cache Configuration

nereus.set_cache_options(max_memory_items=10, disk_path=None)[source]

Configure the global interpolator cache.

Parameters:

  • max_memory_items (int) – Maximum number of interpolators to keep in memory.

  • disk_path (str or Path, optional) – Directory for disk cache. If None, disk caching is disabled.

Model Submodules

Access model-specific functionality through these namespaces:

  • nr.fesom - FESOM2 ocean model (Models Module)

  • nr.healpix - HEALPix grids (Models Module)

  • nr.nemo - NEMO ocean model (Models Module)

  • nr.ifs_tco - IFS TCO mesh (Models Module)

  • nr.icono - ICON-Ocean model (planned)

  • nr.icona - ICON-Atmosphere model (planned)

  • nr.ifs - IFS/ECMWF model (planned)

Version

nereus.__version__

The current version of Nereus.