oceanmesh: Automatic coastal ocean mesh generation

🌊 🌀

Coastal ocean mesh generation from vector and raster GIS data.

1. Quick Start

Get a mesh up and running in minutes. For a full walkthrough, see 5. Basic Usage.

import numpy as np
import oceanmesh as om

# 1) Define a region (WGS84 example)
region = om.Region(extent=(-75.00, -70.00, 40.00, 42.00), crs=4326)

# Alternatively, define an arbitrary polygon extent (lon, lat vertices)
poly_vertices = np.array(
  [
    [-74.2, 40.4],
    [-73.9, 40.4],
    [-73.8, 40.7],
    [-74.1, 40.8],
    [-74.2, 40.4],  # close polygon
  ]
)
poly_region = om.Region(extent=poly_vertices, crs=4326)

# 2) Build shoreline and signed distance function from a coastline shapefile
shore = om.Shoreline("path/to/coastline.shp", poly_region, min_edge_length=0.01)
sdf = om.signed_distance_function(shore)

# 3) Create a sizing function and generate the mesh
edge = om.distance_sizing_function(shore, rate=0.15)
points, cells = om.generate_mesh(sdf, edge)

# 4) Clean up common boundary issues
points, cells = om.make_mesh_boundaries_traversable(points, cells)
points, cells = om.delete_boundary_faces(points, cells, min_qual=0.15)

Back to top

Table of contents

• 1. Quick Start
• 2. Features
• 3. Installation
  • 3.1 Linux/Mac
  • 3.2 Windows
  • 3.3 Development installation
• 4. Support & Community
• 5. Basic Usage
  • 5.1 Setting the Region
  • 5.2 Reading Geophysical Data
  • 5.3 Defining the Domain
  • 5.4 Building Mesh Sizing Functions
  • 5.5 Cleaning up the Mesh
  • 5.6 Mesh Generation
• 6. Advanced Topics
  • 6.1 Multiscale Mesh Generation
  • 6.2 Global and Multiscale Meshing
• 7. Performance Optimization
• 8. Third-Party Code
• 9. Testing
• 10. Citation
• 11. License

2. Features

• A Python package for the development of unstructured triangular meshes used in coastal ocean circulation modeling. The software integrates mesh generation directly with geophysical datasets such as topo-bathymetric rasters/digital elevation models and shapefiles representing coastal features. It provides pre- and post-processing tools to enable successful numerical simulation with the developed model.
  • Automatically handles arbitrarily complex shoreline vector datasets and incorporates them into mesh generation.
  • A variety of commonly used mesh size functions with simple, scriptable controls.
  • Mesh checking and clean-up methods to avoid simulation problems.

Back to top

3. Installation

⚠️ OceanMesh 1.0 provides a stable public API, but the project is still under active development. Check the release notes for details of any breaking changes between minor versions.

The notes below refer to installation on platforms other than MS Windows. For Windows, see 3.2.

3.1 Linux/Mac

oceanmesh needs CGAL:

sudo apt install libcgal-dev

CGAL can also be installed with conda:

conda install -c conda-forge cgal

After that, install or update OceanMesh with pip (recommended for most users):

pip install -U oceanmesh

Prebuilt wheels are available for Apple Silicon (M1+) via pip.

On some clusters/HPC in order to install CGAL, you may need to load/install gmp and mpfr. For example:

sudo apt install libmpfr-dev libgmp3-dev

3.2 Windows

Python on Windows can encounter DLL conflicts due to version incompatibilities among required packages. We provide install_cgal.bat to build a CGAL development distribution separately as a prerequisite.

Prerequisites to build CGAL using the provided batch file:

• Windows 10 or later
• Visual Studio with C++
• CMake
• Git

After successful installation of a CGAL development package, proceed via one of the two options below to generate a Python environment with OceanMesh installed.
performance-critical operations:
If you are using a conda-based Python distribution, then install_oceanmesh.bat should take care of everything, provided no package conflicts arise.

• Point-in-polygon queries are implemented by
  :mod:oceanmesh.geometry.point_in_polygon using a pure-Python
  ray-casting backend with optional fast paths via Shapely and
  Matplotlib when available. When built, an optional Cython extension
  (:mod:oceanmesh.geometry.point_in_polygon_) accelerates the core
  ray-casting kernel and is enabled by default.

  Backend selection follows two environment variables:

  • OCEANMESH_INPOLY_METHOD chooses among the Python backends
    ("raycasting", "shapely", "matplotlib"). When this is
    set, the corresponding backend is used and the Cython kernel is
    not invoked.
  • OCEANMESH_INPOLY_ACCEL controls whether the compiled kernel is
    considered when OCEANMESH_INPOLY_METHOD is not set to a
    recognised method name. By default (unset or non-falsey), OceanMesh
    will attempt to import and use the Cython kernel; if the extension
    is missing or fails at runtime, a warning is logged and the
    code gracefully falls back to the pure-Python implementation.

  If you intend to run without acceleration (for example on a
  platform without a compiler toolchain), set
  OCEANMESH_INPOLY_ACCEL=0 to disable the compiled kernel and
  silence the warning. Otherwise, a warning is emitted to highlight
  that performance may be significantly degraded compared to the
  accelerated path.

Note: CMake is required by vcpkg to build CGAL dependencies, but is not used to build oceanmesh itself (which uses setuptools with pybind11).

3.3 Development installation

To install from source for development and testing:

pip install -e .

Back to top

4. Support & Community

• Questions or problems? Post issues on GitHub or ask in Slack: https://join.slack.com/t/oceanmesh2d/shared_invite/zt-su1q3lh3-C_j6AIOQPrewqZnanhzN7g
• Contact: Dr. Keith Roberts (keithrbt0@gmail.com)
• Version information: oceanmesh uses versioneer.

python -c "import oceanmesh; print(oceanmesh.__version__)"
python setup.py version

Logging during script execution (default is suppressed):

import logging, sys
logging.basicConfig(stream=sys.stdout, level=logging.WARNING)
# logging.basicConfig(stream=sys.stdout, level=logging.INFO)
# logging.basicConfig(stream=sys.stdout, level=logging.DEBUG)

Back to top

5. Basic Usage

5.1 Setting the Region

import oceanmesh as om

EPSG = 32619  # CRS (UTM19N here)
bbox = (-70.29637, -43.56508, -69.65537, 43.88338)
extent = om.Region(extent=bbox, crs=4326)  # bbox given in WGS84
extent = extent.transform_to(EPSG)         # transform to target CRS (UTM19N)
print(extent.bbox)                         # extents now in desired CRS

5.2 Reading Geophysical Data

Shoreline vector datasets (e.g., shapefiles) and digital elevation models (DEMs) are used to construct mesh size and signed distance functions. The dataset download and heavy plotting examples are skipped in CI.

import zipfile
import requests
import oceanmesh as om

url = "http://www.soest.hawaii.edu/pwessel/gshhg/gshhg-shp-2.3.7.zip"
with open("gshhg-shp-2.3.7.zip", "wb") as f:
    f.write(requests.get(url).content)
zipfile.ZipFile("gshhg-shp-2.3.7.zip").extractall("gshhg-shp-2.3.7")

fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
EPSG = 4326
extent = om.Region(extent=(-75.000, -70.001, 40.0001, 41.9000), crs=EPSG)
shoreline = om.Shoreline(fname, extent, 0.01)  # Preferred: pass Region

Working with Projected Coordinate Systems

When working in projected CRSs (e.g., UTM), prefer passing a Region object so both bbox and CRS travel together.

import oceanmesh as om

EPSG = 32610  # UTM Zone 10N
extent = om.Region(extent=(xmin, xmax, ymin, ymax), crs=EPSG)

shore = om.Shoreline(fname, extent, min_edge_length=15)               # carries CRS
shore = om.Shoreline(fname, extent.bbox, 15, crs=EPSG)                 # explicit CRS
# Wrong: bbox in UTM but default CRS=WGS84 (do NOT do this)
# shore = om.Shoreline(fname, extent.bbox, 15)

Best practice: Pass a Region object to Shoreline instead of just a bbox to ensure CRS matches automatically.

# DEM usage (example dataset from datasets/EastCoast.nc)
import oceanmesh as om
dem = om.DEM("datasets/EastCoast.nc", crs=4326)
dem.plot(title="SRTM 30m", vmin=-10, vmax=10)

5.3 Defining the Domain

import oceanmesh as om
fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
extent = om.Region(extent=(-75.00, -70.001, 40.0001, 41.9000), crs=4326)
shoreline = om.Shoreline(fname, extent, 0.01)
sdf = om.signed_distance_function(shoreline)

To flip the inside/outside definition:

sdf = om.signed_distance_function(shoreline, invert=True)

5.4 Building Mesh Sizing Functions

All mesh size functions are defined on regular Cartesian grids. See the Grid class for details.

import oceanmesh as om
fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
extent = om.Region(extent=(-75.00, -70.001, 40.0001, 41.9000), crs=4326)
shoreline = om.Shoreline(fname, extent, 0.01)
edge_length = om.distance_sizing_function(shoreline, rate=0.15)

import oceanmesh as om
fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
extent = om.Region(extent=(-75.00, -70.001, 40.0001, 41.9000), crs=4326)
shoreline = om.Shoreline(fname, extent, 0.01)
sdf = om.signed_distance_function(shoreline)
edge_length = om.feature_sizing_function(shoreline, sdf, max_edge_length=0.05)
edge_length = om.enforce_mesh_gradation(edge_length, gradation=0.15)

import oceanmesh as om
fdem = "datasets/EastCoast.nc"
fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
extent = om.Region(extent=(-74.3, -73.8, 40.3, 40.8), crs=4326)
dem = om.DEM(fdem, bbox=extent, crs=4326)
shoreline = om.Shoreline(fname, dem.bbox, 0.01)
sdf = om.signed_distance_function(shoreline)
edge1 = om.feature_sizing_function(shoreline, sdf, max_edge_length=0.05)
edge2 = om.wavelength_sizing_function(dem, wl=100, period=12.42 * 3600)
edge = om.enforce_mesh_gradation(om.compute_minimum([edge1, edge2]), gradation=0.15)

import oceanmesh as om
fdem = "datasets/EastCoast.nc"
fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
extent = om.Region(extent=(-74.4, -73.4, 40.2, 41.2), crs=4326)
dem = om.DEM(fdem, crs=4326)
shoreline = om.Shoreline(fname, extent, 0.0025)
sdf = om.signed_distance_function(shoreline)
edge_feat = om.feature_sizing_function(shoreline, sdf, max_edge_length=0.10, crs=4326)
edge_grad = om.bathymetric_gradient_sizing_function(
    dem, slope_parameter=5.0, filter_quotient=50, min_edge_length=0.0025, max_edge_length=0.10, crs=4326
)
edge = om.enforce_mesh_gradation(om.compute_minimum([edge_feat, edge_grad]), gradation=0.15)

5.5 Cleaning up the Mesh

points, cells = fix_mesh(points, cells)
points, cells = make_mesh_boundaries_traversable(points, cells)
points, cells = delete_faces_connected_to_one_face(points, cells)
points, cells = delete_boundary_faces(points, cells, min_qual=0.15)
points, cells = laplacian2(points, cells)

5.6 Mesh Generation

Mesh generation uses the DistMesh algorithm and requires only a signed distance and a sizing function.

import meshio
import oceanmesh as om
fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
extent = om.Region(extent=(-75.00, -70.001, 40.0001, 41.9000), crs=4326)
shore = om.Shoreline(fname, extent.bbox, 0.01)
edge = om.distance_sizing_function(shore, max_edge_length=0.05)
domain = om.signed_distance_function(shore)
points, cells = om.generate_mesh(domain, edge)
points, cells = om.make_mesh_boundaries_traversable(points, cells)
points, cells = om.delete_faces_connected_to_one_face(points, cells)
points, cells = om.delete_boundary_faces(points, cells, min_qual=0.15)
points, cells = om.laplacian2(points, cells)
meshio.write_points_cells("new_york.vtk", points, [("triangle", cells)], file_format="vtk")

Back to top

6. Advanced Topics

6.1 Multiscale Mesh Generation

Areas of finer refinement can be incorporated seamlessly by using generate_multiscale_mesh with lists of signed distance and edge length functions. Transitions are blended automatically.

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np

import oceanmesh as om

fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
EPSG = 4326  # EPSG:4326 or WGS84
extent1 = om.Region(extent=(-75.00, -70.001, 40.0001, 41.9000), crs=EPSG)
min_edge_length1 = 0.01  # minimum mesh size in domain in projection
bbox2 = np.array(
  [
    [-73.9481, 40.6028],
    [-74.0186, 40.5688],
    [-73.9366, 40.5362],
    [-73.7269, 40.5626],
    [-73.7231, 40.6459],
    [-73.8242, 40.6758],
    [-73.9481, 40.6028],
  ],
  dtype=float,
)
extent2 = om.Region(extent=bbox2, crs=EPSG)
min_edge_length2 = 4.6e-4  # minimum mesh size in domain in projection
s1 = om.Shoreline(fname, extent1.bbox, min_edge_length1)
sdf1 = om.signed_distance_function(s1)
el1 = om.distance_sizing_function(s1, max_edge_length=0.05)
s2 = om.Shoreline(fname, extent2.bbox, min_edge_length2)
sdf2 = om.signed_distance_function(s2)
el2 = om.distance_sizing_function(s2)
# Control the element size transition
# from coarse to fine with the kwargs prefixed with `blend`
points, cells = om.generate_multiscale_mesh(
  [sdf1, sdf2],
  [el1, el2],
)
# Remove degenerate mesh faces and other common problems in the mesh
points, cells = om.make_mesh_boundaries_traversable(points, cells)
# Remove singly connected elements (elements connected to only one other element)
points, cells = om.delete_faces_connected_to_one_face(points, cells)
# Remove poor boundary elements with quality < 15%
points, cells = om.delete_boundary_faces(points, cells, min_qual=0.15)
# Apply a Laplacian smoother that preservers the mesh size distribution
points, cells = om.laplacian2(points, cells)

# Plot it showing the different levels of resolution
triang = tri.Triangulation(points[:, 0], points[:, 1], cells)
gs = gridspec.GridSpec(2, 2)
gs.update(wspace=0.5)
plt.figure()

bbox3 = np.array(
  [
    [-73.78, 40.60],
    [-73.75, 40.60],
    [-73.75, 40.64],
    [-73.78, 40.64],
    [-73.78, 40.60],
  ],
  dtype=float,
)

ax = plt.subplot(gs[0, 0])
ax.set_aspect("equal")
ax.triplot(triang, "-", lw=1)
ax.plot(bbox2[:, 0], bbox2[:, 1], "r--")
ax.plot(bbox3[:, 0], bbox3[:, 1], "m--")

ax = plt.subplot(gs[0, 1])
ax.set_aspect("equal")
ax.triplot(triang, "-", lw=1)
ax.plot(bbox2[:, 0], bbox2[:, 1], "r--")
ax.set_xlim(np.amin(bbox2[:, 0]), np.amax(bbox2[:, 0]))
ax.set_ylim(np.amin(bbox2[:, 1]), np.amax(bbox2[:, 1]))
ax.plot(bbox3[:, 0], bbox3[:, 1], "m--")

ax = plt.subplot(gs[1, :])
ax.set_aspect("equal")
ax.triplot(triang, "-", lw=1)
ax.set_xlim(-73.78, -73.75)
ax.set_ylim(40.60, 40.64)
plt.show()

6.2 Global and Multiscale Meshing

Global meshes are defined in EPSG:4326 but meshed in a stereographic projection. Regional refinement can be added as additional domains.

Global mesh generation (two-step: EPSG:4326 sizing → stereographic meshing)

Global mesh generation is done in two steps:

1. Define the shoreline and sizing functions in EPSG:4326.
2. Generate the mesh in a stereographic projection using a stereographic coastline.

The repository includes example global shoreline shapefiles under tests/global/:

• tests/global/global_latlon.shp: shoreline in EPSG:4326 (lon/lat)
• tests/global/global_stereo.shp: shoreline already transformed for stereographic meshing

Note: global_stereo.shp can be produced using global_tag() in pyPoseidon:
https://github.com/ec-jrc/pyPoseidon/blob/9cfd3bbf5598c810004def83b1f43dc5149addd0/pyposeidon/boundary.py#L452

import numpy as np
import oceanmesh as om
from oceanmesh.region import to_lat_lon
import matplotlib.pyplot as plt


def crosses_dateline(lon1, lon2):
  return abs(lon1 - lon2) > 180


def filter_triangles(points_lonlat, cells):
  """Drop triangles that cross the dateline to avoid plot artifacts."""
  filtered = []
  for cell in cells:
    p1, p2, p3 = (
      points_lonlat[cell[0]],
      points_lonlat[cell[1]],
      points_lonlat[cell[2]],
    )
    if not (
      crosses_dateline(p1[0], p2[0])
      or crosses_dateline(p2[0], p3[0])
      or crosses_dateline(p3[0], p1[0])
    ):
      filtered.append(cell)
  return filtered


# WGS84 shoreline and stereographic shoreline
fname_wgs84 = "tests/global/global_latlon.shp"
fname_stereo = "tests/global/global_stereo.shp"

extent = om.Region(extent=(-180.0, 180.0, -89.0, 90.0), crs=4326)

# 1) Define sizing functions in WGS84
min_edge_length = 0.5
shoreline = om.Shoreline(fname_wgs84, extent.bbox, min_edge_length)
edge_length = om.distance_sizing_function(shoreline, rate=0.11)

# 2) Mesh in stereographic projection using a stereographic shoreline
shoreline_stereo = om.Shoreline(fname_stereo, extent.bbox, min_edge_length, stereo=True)
domain = om.signed_distance_function(shoreline_stereo)

points, cells = om.generate_mesh(domain, edge_length, stereo=True, max_iter=100)

# Clean up and smooth
points, cells = om.make_mesh_boundaries_traversable(points, cells)
points, cells = om.delete_faces_connected_to_one_face(points, cells)
points, cells = om.laplacian2(points, cells, max_iter=100)

# Convert back to lon/lat for plotting
lon, lat = to_lat_lon(points[:, 0], points[:, 1])
tri_cells = filter_triangles(np.array([lon, lat]).T, cells)

fig, ax, pc = edge_length.plot(
  holding=True,
  plot_colorbar=True,
  cbarlabel="Resolution in °",
  cmap="magma",
)
ax.triplot(lon, lat, tri_cells, color="w", linewidth=0.25)
plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import matplotlib.gridspec as gridspec
import oceanmesh as om
from oceanmesh.region import to_lat_lon

fname_global_latlon = "tests/global/global_latlon.shp"
fname_global_stereo = "tests/global/global_stereo.shp"

global_region = om.Region(extent=(-180.0, 180.0, -89.0, 90.0), crs=4326)
shoreline_global_latlon = om.Shoreline(fname_global_latlon, global_region, 1.0)
sdf_global_latlon = om.signed_distance_function(shoreline_global_latlon)
edge_global = om.enforce_mesh_gradation(
    om.compute_minimum([
        om.distance_sizing_function(shoreline_global_latlon, rate=0.11),
        om.feature_sizing_function(shoreline_global_latlon, sdf_global_latlon, max_edge_length=3.0),
    ]),
    gradation=0.15,
    stereo=True,
)

aus_region = om.Region(extent=(110.0, 160.0, -45.0, -10.0), crs=4326)
shoreline_regional = om.Shoreline(fname_global_latlon, aus_region, 0.25)
sdf_regional = om.signed_distance_function(shoreline_regional)
edge_regional = om.enforce_mesh_gradation(
    om.compute_minimum([
        om.distance_sizing_function(shoreline_regional, rate=0.13),
        om.feature_sizing_function(shoreline_regional, sdf_regional, max_edge_length=1.5),
    ]),
    gradation=0.12,
)

shoreline_global_stereo = om.Shoreline(fname_global_stereo, global_region, 1.0, stereo=True)
sdf_global_stereo = om.signed_distance_function(shoreline_global_stereo)

points, cells = om.generate_multiscale_mesh(
    [sdf_global_stereo, sdf_regional],
    [edge_global, edge_regional],
    blend_width=1.0e6,
    blend_max_iter=50,
    max_iter=75,
)

*The image shows the global mesh with a refined Australia region.

See the tests in the tests/ folder for more inspiration; work is ongoing on this package.

Back to top

7. Performance Optimization

OceanMesh uses efficient, GPL-compatible geometry backends for
performance-critical operations:

• Point-in-polygon queries use the new
  oceanmesh.geometry.inpoly2 implementation, which can automatically
  take advantage of Shapely prepared geometries or Matplotlib path
  operations when those libraries are available, falling back to a
  portable pure-Python ray-casting algorithm otherwise.
• Delaunay triangulation is provided by a pure-Python
  Bowyer–Watson implementation with an optional Cython-accelerated
  kernel, enabled automatically when built.

No special extras or build flags are required to enable these
optimizations; the fastest available backend is selected at runtime
based on the installed dependencies.

Back to top

8. Third-Party Code

OceanMesh relies on a number of well-established open-source
dependencies (see setup.cfg for the full list), but does not
currently vendor any third-party geometry libraries. Earlier releases
included a vendored copy of inpoly-python; this has been fully
replaced by the native GPL-compatible implementation in
oceanmesh.geometry.inpoly2.

Back to top

9. Testing

To run the oceanmesh unit tests (and turn off plots), check out this repository and run tox.

Back to top

10. Citation

[1] - Roberts, K. J., Pringle, W. J., and Westerink, J. J., 2019.
      OceanMesh2D 1.0: MATLAB-based software for two-dimensional unstructured mesh generation in coastal ocean modeling,
      Geoscientific Model Development, 12, 1847-1868. https://doi.org/10.5194/gmd-12-1847-2019.

Back to top

11. License

This software is published under the GPLv3 license