Working with Meshes

FESOMP provides a comprehensive Mesh class for working with FESOM2 unstructured meshes.

Loading Meshes

From NetCDF

The recommended format is NetCDF, typically the fesom.mesh.diag.nc file:

import fesomp

mesh = fesomp.load_mesh("path/to/fesom.mesh.diag.nc")

From ASCII

You can also load meshes from ASCII files (the traditional FESOM format):

# Point to the directory containing nod2d.out, elem2d.out, etc.
mesh = fesomp.load_mesh("path/to/mesh/directory/")

The loader will automatically detect the format based on whether you provide a file or directory.

Mesh Properties

Basic Properties

mesh.n2d       # Number of 2D nodes
mesh.nelem     # Number of triangular elements
mesh.nlev      # Number of vertical levels

mesh.lon       # Node longitudes (n2d,)
mesh.lat       # Node latitudes (n2d,)

mesh.depth_levels  # Depth at level interfaces (nlev,)
mesh.depth_layers  # Depth at layer centers (nlev-1,)

Element Coordinates (Lazy)

Element center coordinates are computed on first access:

mesh.lon_elem  # Element center longitudes (nelem,)
mesh.lat_elem  # Element center latitudes (nelem,)

Triangulation

mesh.elem      # Triangle connectivity (nelem, 3)

Topology (Lazy)

Mesh topology is computed on first access for efficiency:

topology = mesh.topology

topology.edges           # Unique edges (nedges, 2)
topology.edge_indices    # Edge index for each triangle edge
topology.face_neighbors  # Adjacent triangles (-1 for boundary)
topology.boundary_edges  # Indices of boundary edges
topology.n_edges         # Total number of edges

Geometry (Lazy)

Geometric properties are computed on demand:

geometry = mesh.geometry

geometry.elem_area   # Triangle areas in m² (nelem,)
geometry.node_area   # Control volumes around nodes (3, n2d)
                     # Three computation methods
geometry.gradients   # Gradient operators for elements

Spatial Queries

FESOMP provides efficient spatial indexing for queries.

Find Nearest Nodes

# Find the 5 nearest nodes to a point
indices = mesh.find_nearest(lon=10.5, lat=54.3, k=5)

# Find nearest element centers
indices = mesh.find_nearest(lon=10.5, lat=54.3, k=5, on_elements=True)

Radius Search

# Find all nodes within 100 km of a point
indices = mesh.find_in_radius(lon=0, lat=0, radius_km=100)

Bounding Box

# Find nodes in a region
indices = mesh.subset_by_bbox(
    lon_min=-10, lon_max=10,
    lat_min=40, lat_max=60
)

Automatic Data Location Detection

FESOMP automatically detects whether data is on nodes or elements:

# Data on nodes (n2d points)
temp_nodes = np.random.rand(mesh.n2d)
# Functions will use mesh.lon, mesh.lat

# Data on elements (nelem points)
velocity = np.random.rand(mesh.nelem)
# Functions will use mesh.lon_elem, mesh.lat_elem

# 3D data on levels (nlev, n2d)
w_velocity = np.random.rand(mesh.nlev, mesh.n2d)
# Uses mesh.depth_levels

# 3D data on layers (nlev-1, n2d)
temperature = np.random.rand(mesh.nlev - 1, mesh.n2d)
# Uses mesh.depth_layers

This detection happens automatically in functions like transect() and plot().

Coordinate Transformations

FESOM2 meshes often use rotated coordinate systems to avoid singularities at the poles. FESOMP provides functions to convert between rotated and geographical coordinates.

Scalar Coordinate Conversion

Convert coordinates between rotated and geographical systems:

from fesomp.mesh import scalar_r2g, scalar_g2r
import numpy as np

# Euler angles defining the rotation (typical FESOM values)
alpha, beta, gamma = 50, 15, -90

# Convert from rotated to geographical coordinates
rlon = np.array([0.0, 10.0, 20.0])
rlat = np.array([45.0, -30.0, 60.0])
lon, lat = scalar_r2g(alpha, beta, gamma, rlon, rlat)

# Convert from geographical to rotated coordinates
rlon_back, rlat_back = scalar_g2r(alpha, beta, gamma, lon, lat)

Vector Rotation

Rotate velocity or other vector fields between coordinate systems:

from fesomp.mesh import vec_rotate_r2g, vec_rotate_g2r

# Euler angles
alpha, beta, gamma = 50, 15, -90

# Coordinates where vectors are defined
lon = np.array([10.0, 20.0, -30.0])
lat = np.array([45.0, -30.0, 60.0])

# Vector components in rotated coordinates
u_rot = np.array([1.0, 2.0, -1.0])  # eastward component
v_rot = np.array([0.5, -0.5, 1.5])  # northward component

# Rotate vectors to geographical coordinates
# flag=1 means lon/lat are in geographical coordinates
# flag=0 means lon/lat are in rotated coordinates
u_geo, v_geo = vec_rotate_r2g(alpha, beta, gamma, lon, lat, u_rot, v_rot, flag=1)

# Rotate vectors from geographical to rotated coordinates
u_rot_back, v_rot_back = vec_rotate_g2r(alpha, beta, gamma, lon, lat, u_geo, v_geo, flag=1)

The flag parameter specifies the coordinate system of the input lon/lat:

flag=1: lon/lat are in geographical coordinates
flag=0: lon/lat are in rotated coordinates

The functions will automatically compute the coordinates in the other system as needed.