Diagnostics Module

The nereus.diag module provides diagnostic functions for analyzing ocean and sea ice model output.

Sea Ice Diagnostics

Core functions (ice_area, ice_volume, ice_extent) compute total metrics with optional masking. Hemisphere convenience functions (*_nh, *_sh) provide shortcuts for Northern and Southern Hemisphere calculations.

Sea ice diagnostics for nereus.

This module provides functions for computing sea ice metrics: - ice_area: Total sea ice area - ice_volume: Total sea ice volume - ice_extent: Sea ice extent (area with concentration above threshold)

Hemisphere convenience functions are also provided: - ice_area_nh, ice_area_sh: Northern/Southern Hemisphere ice area - ice_volume_nh, ice_volume_sh: Northern/Southern Hemisphere ice volume - ice_extent_nh, ice_extent_sh: Northern/Southern Hemisphere ice extent

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

nereus.diag.ice.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.diag.ice.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.diag.ice.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.diag.ice.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.diag.ice.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.diag.ice.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.diag.ice.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.diag.ice.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.diag.ice.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)

Vertical/Ocean Diagnostics

Functions for computing ocean diagnostics:

  • surface_mean: Area-weighted mean for 2D fields (SST, SSS, single depth levels)

  • volume_mean: Volume-weighted mean for 3D fields

  • heat_content: Ocean heat content (total or map)

  • find_closest_depth: Find index and value of closest depth level to target

  • interpolate_to_depth: Interpolate 3D data to target depth levels

Vertical/ocean diagnostics for nereus.

This module provides functions for computing ocean diagnostics: - surface_mean: Area-weighted mean for 2D fields (SST, SSS, etc.) - volume_mean: Volume-weighted mean in a depth range - heat_content: Ocean heat content - find_closest_depth: Find index and value of closest depth to target - interpolate_to_depth: Interpolate 3D data to target depths

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

nereus.diag.vertical.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.diag.vertical.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.diag.vertical.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
nereus.diag.vertical.find_closest_depth(depths, target)[source]

Find the index and value of the depth closest to a target depth.

This is useful when comparing multiple models with different depth levels to find corresponding levels, and to assess how far model depths are from target depths.

Parameters:

  • depths (array_like) – 1D array of depth values (typically positive downward in meters).

  • target (float) – Target depth value to find the closest match for.

Returns:

A tuple of (index, value) where index is the position of the closest depth in the input array, and value is the actual depth at that index.

Return type:

tuple[int, float]

Examples

>>> depths = [0, 10, 25, 50, 100, 200, 500, 1000]
>>> idx, val = nr.find_closest_depth(depths, 100)
>>> print(f"Index: {idx}, Depth: {val}m")
Index: 4, Depth: 100.0m

>>> # Check how far model depth is from target
>>> idx, val = nr.find_closest_depth(depths, 75)
>>> print(f"Closest depth: {val}m, difference: {abs(val - 75)}m")
Closest depth: 50.0m, difference: 25.0m
nereus.diag.vertical.interpolate_to_depth(data, lon, lat, model_depths, target_depths)[source]

Interpolate 3D data to target depth levels using linear interpolation.

Performs column-wise linear interpolation from model depth levels to specified target depths. Values outside the model depth range are extrapolated (with a warning for significant extrapolation).

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).

  • lon (array_like or None) – Longitude coordinates, shape (npoints,). If provided along with lat, these are returned with the result for convenience. Pass None if not needed.

  • lat (array_like or None) – Latitude coordinates, shape (npoints,). If provided along with lon, these are returned with the result for convenience. Pass None if not needed.

  • model_depths (array_like) – Depth levels of the input data in meters (positive downward), shape (nlevels,).

  • target_depths (array_like or float) – Target depth(s) to interpolate to in meters. Can be a single value or an array of depths.

Returns:

If lon and lat are None:

Interpolated data with shape (ntargets, npoints) or (…, ntargets, npoints) for higher-dimensional input. If target_depths is a scalar, ntargets=1.

If lon and lat are provided:

Tuple of (interpolated_data, lon, lat).

Return type:

ndarray or tuple

Examples

>>> # Interpolate temperature to 100m depth (without coordinates)
>>> temp_100m = nr.interpolate_to_depth(temp, None, None, mesh.depth, 100)

>>> # Interpolate to multiple standard depths
>>> standard_depths = [10, 50, 100, 200, 500, 1000]
>>> temp_interp = nr.interpolate_to_depth(temp, None, None, mesh.depth, standard_depths)

>>> # With coordinates for plotting
>>> temp_100m, lon, lat = nr.interpolate_to_depth(
...     temp, mesh.lon, mesh.lat, mesh.depth, 100
... )
>>> nr.plot(temp_100m.squeeze(), lon, lat)

>>> # Compare models at the same depth
>>> temp_model1 = nr.interpolate_to_depth(temp1, None, None, depths1, 100)
>>> temp_model2 = nr.interpolate_to_depth(temp2, None, None, depths2, 100)

Constants

  • RHO_SEAWATER: Reference seawater density: 1025.0 kg/m³

  • CP_SEAWATER: Seawater specific heat capacity: 3990.0 J/(kg·K) (consistent with FESOM2)

Hovmoller Diagrams

Hovmoller diagram generation for nereus.

This module provides functions for computing and plotting Hovmoller diagrams (time-depth or time-latitude plots).

The hovmoller function is dask-friendly: if inputs are dask arrays, the result will be lazy dask arrays for both mode=”depth” and mode=”latitude”.

nereus.diag.hovmoller.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.diag.hovmoller.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})

Region Masks

Region mask utilities for geographic data analysis.

This module provides functions for loading GeoJSON region definitions and creating boolean masks for points within named regions.

nereus.diag.regions.load_geojson(name)[source]

Load a GeoJSON file from the package data directory.

Parameters:

name (str) – Name of the GeoJSON file (without .geojson extension). Available files: MOCBasins, NinoRegions, oceanBasins

Returns:

Parsed GeoJSON data as a dictionary.

Return type:

dict

Raises:

FileNotFoundError – If the GeoJSON file does not exist.

nereus.diag.regions.list_available_regions(geojson_name='MOCBasins')[source]

List available region names in a GeoJSON file.

Parameters:

geojson_name (str, optional) – Name of the GeoJSON file (without extension). Default is “MOCBasins”.

Returns:

List of region names defined in the GeoJSON file.

Return type:

list[str]

nereus.diag.regions.get_region_mask(lon, lat, region, geojson_name='MOCBasins')[source]

Create a boolean mask for points within a named region.

Parameters:

  • lon (np.ndarray) – Array of longitude values.

  • lat (np.ndarray) – Array of latitude values.

  • region (str) – Name of the region to create a mask for.

  • geojson_name (str, optional) – Name of the GeoJSON file (without extension). Default is “MOCBasins”.

Returns:

Boolean array where True indicates points inside the region.

Return type:

np.ndarray

Raises:

  • ImportError – If shapely is not installed.

  • ValueError – If the region is not found in the GeoJSON file.

Built-in GeoJSON Files

Nereus includes the following GeoJSON files with predefined regions:

MOCBasins - Meridional Overturning Circulation regions:

  • Atlantic_MOC

  • IndoPacific_MOC

  • Pacific_MOC

  • Indian_MOC

oceanBasins - Major ocean basins:

  • Atlantic_Basin

  • Pacific_Basin

  • Indian_Basin

  • Arctic_Basin

  • Southern_Ocean_Basin

  • Mediterranean_Basin

  • Global Ocean

  • Global Ocean 65N to 65S

  • Global Ocean 15S to 15N

NinoRegions - El Nino monitoring regions:

  • Nino 3.4

  • Nino 3

  • Nino 4

Formulas

This section documents the mathematical formulas used in the diagnostic calculations.

Sea Ice Area

\[A_{ice} = \sum_i c_i \cdot A_i\]

where:

  • \(c_i\) = ice concentration at cell \(i\) (fraction, 0-1)

  • \(A_i\) = area of cell \(i\) (m²)

Sea Ice Extent

\[E_{ice} = \sum_i A_i \cdot \mathbf{1}_{c_i \geq \tau}\]

where:

  • \(\tau\) = concentration threshold (default 0.15)

  • \(\mathbf{1}_{(\cdot)}\) = indicator function

Sea Ice Volume

\[V_{ice} = \sum_i h_i \cdot c_i \cdot A_i\]

where:

  • \(h_i\) = ice thickness at cell \(i\) (m)

  • \(c_i\) = ice concentration (optional weighting)

  • \(A_i\) = cell area (m²)

Surface Mean

\[\bar{X} = \frac{\sum_i X_i \cdot A_i}{\sum_i A_i}\]

where:

  • \(X_i\) = quantity at horizontal cell \(i\)

  • \(A_i\) = horizontal cell area (m²)

Volume Mean

\[\bar{X} = \frac{\sum_{i,k} X_{i,k} \cdot A_i \cdot \Delta z_k}{\sum_{i,k} A_i \cdot \Delta z_k}\]

where:

  • \(X_{i,k}\) = quantity at horizontal cell \(i\), vertical level \(k\)

  • \(A_i\) = horizontal cell area (m²)

  • \(\Delta z_k\) = layer thickness at level \(k\) (m)

Ocean Heat Content

The heat_content function supports two output modes:

Total heat content (output="total", default):

\[OHC = \rho \cdot c_p \cdot \sum_{i,k} (T_{i,k} - T_{ref}) \cdot A_i \cdot \Delta z_k\]

Returns total ocean heat content in Joules.

Heat content map (output="map"):

\[OHC_i = \rho \cdot c_p \cdot \sum_{k} (T_{i,k} - T_{ref}) \cdot \Delta z_k\]

Returns heat content per unit area at each horizontal point in J/m² (consistent with FESOM2 output).

where:

  • \(\rho\) = seawater density (default 1025 kg/m³)

  • \(c_p\) = specific heat capacity (default 3990 J/(kg·K))

  • \(T_{i,k}\) = temperature at cell \(i\), level \(k\) (°C)

  • \(T_{ref}\) = reference temperature (default 0°C)

  • \(A_i\) = cell area (m²)

  • \(\Delta z_k\) = layer thickness (m)

Find Closest Depth

\[i^* = \arg\min_i |z_i - z_{target}|\]

Returns index \(i^*\) and value \(z_{i^*}\) of the closest depth level.

Linear Depth Interpolation

For a target depth \(z_{target}\) between model levels \(z_k\) and \(z_{k+1}\):

\[X(z_{target}) = X_k + \frac{z_{target} - z_k}{z_{k+1} - z_k} \cdot (X_{k+1} - X_k)\]

where:

  • \(X_k\) = value at depth level \(k\)

  • \(z_k\) = depth of level \(k\) (m)