From f138752f48d72f1e4001809c8e173614afb7b29a Mon Sep 17 00:00:00 2001 From: Keith Roberts Date: Mon, 1 Dec 2025 21:08:44 -0500 Subject: [PATCH 1/4] Add optional Cartopy/pyproj-backed StereoProjection for global meshing, refactor region transforms to use it, and fully propagate scale_factor (k0) from Shoreline/Domain into generate_mesh and stereographic distortion routines, fixing syntax issues in mesh_generator. --- oceanmesh/__init__.py | 9 +- oceanmesh/edgefx.py | 53 +++++++-- oceanmesh/geodata.py | 25 +++- oceanmesh/grid.py | 4 + oceanmesh/mesh_generator.py | 110 ++++++++++++++--- oceanmesh/projections.py | 165 ++++++++++++++++++++++++++ oceanmesh/region.py | 57 ++++++++- oceanmesh/signed_distance_function.py | 61 +++++++--- setup.cfg | 3 + tests/test_global_stereo.py | 68 +++++++++++ 10 files changed, 502 insertions(+), 53 deletions(-) create mode 100644 oceanmesh/projections.py diff --git a/oceanmesh/__init__.py b/oceanmesh/__init__.py index e42a6a5..3f5c19f 100644 --- a/oceanmesh/__init__.py +++ b/oceanmesh/__init__.py @@ -126,6 +126,11 @@ __version__ = _version.get_versions()["version"] -from . import _version +try: # Optional global-stereo helpers (import may fail generically) + from .projections import StereoProjection, CARTOPY_AVAILABLE -__version__ = _version.get_versions()["version"] + __all__.extend(["StereoProjection", "CARTOPY_AVAILABLE"]) +except ImportError: + # Projections module not importable; expose CARTOPY_AVAILABLE flag as False + CARTOPY_AVAILABLE = False + __all__.append("CARTOPY_AVAILABLE") diff --git a/oceanmesh/edgefx.py b/oceanmesh/edgefx.py index 0f0cddd..1bb6a91 100644 --- a/oceanmesh/edgefx.py +++ b/oceanmesh/edgefx.py @@ -88,7 +88,9 @@ def enforce_mesh_size_bounds_elevation(grid, dem, bounds): return grid -def enforce_mesh_gradation(grid, gradation=0.15, crs="EPSG:4326", stereo=False): +def enforce_mesh_gradation( + grid, gradation=0.15, crs="EPSG:4326", stereo=False, scale_factor=1.0 +): """Enforce a mesh size gradation bound `gradation` on a :class:`grid` Parameters @@ -124,17 +126,40 @@ def enforce_mesh_gradation(grid, gradation=0.15, crs="EPSG:4326", stereo=False): tmp = gradient_limit([*sz], elen, gradation, 10000, cell_size) tmp = np.reshape(tmp, (sz[0], sz[1]), "F") if stereo: - logger.info("Global mesh: fixing gradient on the north pole...") - # max distortion at the pole: 2 / 180 * PI / (1 - cos(lat))**2 - dx_stereo = grid.dx * 1 / 180 * np.pi / 2 - # in stereo projection, all north hemisphere is contained in the unit sphere - # we want to fix the gradient close to the north pole, - # so we extract all the coordinates between -1 and 1 in stereographic projection + from .projections import CARTOPY_AVAILABLE + from .region import _get_stereo_projection + + logger.info( + "Global mesh: fixing gradient on the north pole (using %s projection, k0=%s)...", + "cartopy" if CARTOPY_AVAILABLE else "hardcoded", + scale_factor, + ) + # Obtain a stereographic projection configured with the + # caller-provided scale_factor so that k0 meaningfully affects + # the polar correction. + proj = _get_stereo_projection(scale_factor=scale_factor) + + # Choose a reasonable physical radius (in metres) around the + # north pole within which to regularise the gradient. Here we + # use ~2000 km. + r_max = 2_000_000.0 + + # Resolution in stereographic metres, based on underlying lon/lat step. + dx_deg = grid.dx + dx_stereo = dx_deg * np.pi / 180.0 * 6_371_000.0 us, vs = np.meshgrid( - np.arange(-1, 1, dx_stereo), np.arange(-1, 1, dx_stereo), indexing="ij" + np.arange(-r_max, r_max, dx_stereo), + np.arange(-r_max, r_max, dx_stereo), + indexing="ij", ) - ulon, vlat = to_lat_lon(us.ravel(), vs.ravel()) + + if proj is not None: + ulon, vlat = proj.to_lat_lon(us.ravel(), vs.ravel()) + else: + # Fall back to module-level transform if projection could + # not be constructed. + ulon, vlat = to_lat_lon(us.ravel(), vs.ravel()) utmp = grid.eval((ulon, vlat)) utmp = np.reshape(utmp, us.shape) szs = utmp.shape @@ -143,9 +168,9 @@ def enforce_mesh_gradation(grid, gradation=0.15, crs="EPSG:4326", stereo=False): # this is merely to fix the north pole gradient vtmp = gradient_limit([*szs], dx_stereo, 10, 10000, utmp.flatten("F")) vtmp = np.reshape(vtmp, (szs[0], szs[1]), "F") - # construct stereo interpolating function + # construct stereo interpolating function in projection metres grid_stereo = Grid( - bbox=(-1, 1, -1, 1), + bbox=(-r_max, r_max, -r_max, r_max), dx=dx_stereo, values=vtmp, hmin=grid.hmin, @@ -155,7 +180,11 @@ def enforce_mesh_gradation(grid, gradation=0.15, crs="EPSG:4326", stereo=False): grid_stereo.build_interpolant() # reinject back into the original grid and redo the gradient computation xg, yg = grid.create_grid() - tmp[yg > 0] = grid_stereo.eval(to_stereo(xg[yg > 0], yg[yg > 0])) + if proj is not None: + xs, ys = proj.to_stereo(xg[yg > 0], yg[yg > 0]) + else: + xs, ys = to_stereo(xg[yg > 0], yg[yg > 0]) + tmp[yg > 0] = grid_stereo.eval((xs, ys)) logger.info( "Global mesh: reinject back stereographic gradient and recomputing gradient..." ) diff --git a/oceanmesh/geodata.py b/oceanmesh/geodata.py index c98cfc3..15b8ebc 100644 --- a/oceanmesh/geodata.py +++ b/oceanmesh/geodata.py @@ -577,6 +577,7 @@ def __init__( minimum_area_mult=4.0, smooth_shoreline=True, stereo=False, + scale_factor=1.0, ): if isinstance(shp, str): shp = Path(shp) @@ -634,14 +635,23 @@ def __init__( self.shp = shp self.h0 = h0 - # Retain stereo flag for downstream validation (e.g., multiscale mixing) + # Retain stereo flag for downstream validation (e.g., multiscale mixing). + # When True, stereographic transformations will use cartopy's + # NorthPolarStereo when available (see oceanmesh.projections). self.stereo = bool(stereo) + if scale_factor <= 0: + raise ValueError("scale_factor must be > 0") + self.__scale_factor = float(scale_factor) self.inner = [] self.outer = [] self.mainland = [] self.boubox = _boubox self.refinements = refinements self.minimum_area_mult = minimum_area_mult + if self.stereo: + from oceanmesh.projections import check_cartopy_available + + check_cartopy_available() polys = self._read() @@ -708,6 +718,19 @@ def h0(self, value): raise ValueError("h0 must be > 0") self.__h0 = value + @property + def scale_factor(self): + """Reference stereographic scale factor k0 for global meshes.""" + + return self.__scale_factor + + @scale_factor.setter + def scale_factor(self, value): + value = float(value) + if value <= 0: + raise ValueError("scale_factor must be > 0") + self.__scale_factor = value + @staticmethod def transform_to(gdf, dst_crs): """Transform geodataframe ``gdf`` representing diff --git a/oceanmesh/grid.py b/oceanmesh/grid.py index 4ff7c78..0af69de 100644 --- a/oceanmesh/grid.py +++ b/oceanmesh/grid.py @@ -455,6 +455,10 @@ def plot( """ _xg, _yg = self.create_grid() if stereo: + # Transform to stereographic coordinates for plotting. + # Uses cartopy's NorthPolarStereo when available + # (pip install oceanmesh[global]) and falls back to + # analytic formulas otherwise. _xg, _yg = to_stereo(_xg, _yg) if ax is None: fig, ax = plt.subplots() diff --git a/oceanmesh/mesh_generator.py b/oceanmesh/mesh_generator.py index d866b2f..4595d11 100644 --- a/oceanmesh/mesh_generator.py +++ b/oceanmesh/mesh_generator.py @@ -104,7 +104,6 @@ def write_to_fort14( np.column_stack((k + 1, points[k][0], points[k][1], topobathymetry[k])), delimiter=" ", fmt="%i %f %f %f", - newline="\n", ) # Write element connectivity @@ -726,9 +725,15 @@ def generate_mesh(domain, edge_length, **kwargs): opts.update(kwargs) _parse_kwargs(kwargs) - fd, bbox = _unpack_domain(domain, opts) + fd, bbox, scale_factor = _unpack_domain(domain, opts) fh, min_edge_length = _unpack_sizing(edge_length, opts) + # Propagate stereographic scale factor (k0) from the signed-distance + # domain into the meshing options so that global stereo workflows can + # correctly apply distortion corrections in both the initial point + # generation and force computation routines. + opts["scale_factor"] = scale_factor + _check_bbox(bbox) bbox = np.array(bbox).reshape(-1, 2) @@ -756,6 +761,7 @@ def generate_mesh(domain, edge_length, **kwargs): fd, pfix, opts["stereo"], + k0=scale_factor, ) else: p = opts["points"] @@ -863,17 +869,35 @@ def _unpack_sizing(edge_length, opts): def _unpack_domain(domain, opts): + """Unpack a domain into (fd, bbox, scale_factor). + + Parameters + ---------- + domain : Domain or callable + Signed distance function domain object or a raw callable. + opts : dict + Meshing options, expected to contain ``bbox`` when ``domain`` is a + callable. + + Returns + ------- + fd : callable + The signed distance function evaluator. + bbox : tuple + The bounding box associated with the domain. + scale_factor : float + The stereographic reference scale factor k0 (defaults to 1.0 when + not provided on the domain). + """ + if isinstance(domain, Domain): - bbox = domain.bbox - fd = domain.eval + return domain.eval, domain.bbox, getattr(domain, "scale_factor", 1.0) elif callable(domain): - bbox = opts["bbox"] - fd = domain + return domain, opts["bbox"], 1.0 else: raise ValueError( "`domain` must be a function or a :class:`signed_distance_function object" ) - return fd, bbox def _get_bars(t): @@ -891,13 +915,14 @@ def _compute_forces(p, t, fh, min_edge_length, L0mult, opts): L = np.sqrt((barvec**2).sum(1)) # L = Bar lengths L[L == 0] = np.finfo(float).eps if opts["stereo"]: + k0 = opts.get("scale_factor", 1.0) # For global+regional multiscale meshes, this branch handles the global stereo case. # Regional sizing functions have been wrapped or transformed earlier so fh(p2) # evaluates correctly on lat/lon even though points are maintained in stereo space. p1 = p[bars].sum(1) / 2 x, y = to_lat_lon(p1[:, 0], p1[:, 1]) p2 = np.asarray([x, y]).T - hbars = fh(p2) * _stereo_distortion_dist(y) + hbars = fh(p2) * _stereo_distortion_dist(y, k0=k0) else: hbars = fh(p[bars].sum(1) / 2) # Guard against non-finite or non-positive sizing values that can poison forces @@ -1004,26 +1029,71 @@ def _deps_vec(i): return p -def _stereo_distortion(lat): - # we use here Stereographic projection of the sphere - # from the north pole onto the plane - # https://en.wikipedia.org/wiki/Stereographic_projection +def _stereo_distortion(lat, k0=1.0): + """Return stereographic scale factor at latitude. + + When cartopy is available, this uses + :class:`oceanmesh.projections.StereoProjection` to compute the + scale factor via + + .. math:: k(\phi) = 2 k_0 / (1 + \sin \phi), + + with :math:`k_0 = 1`. Otherwise, it falls back to the historical + analytic expression based on a north-polar stereographic + projection of the unit sphere. + """ + + try: + from .region import _get_stereo_projection + + proj = _get_stereo_projection(scale_factor=k0) + if proj is not None: + return proj.get_scale_factor(lat) + except Exception: + # Fall back to legacy formula if projections cannot be used + pass + lat0 = 90 ll = lat + lat0 lrad = ll / 180 * np.pi - res = 2 / (1 + np.sin(lrad)) + res = 2 * k0 / (1 + np.sin(lrad)) return res -def _stereo_distortion_dist(lat): +def _stereo_distortion_dist(lat, k0=1.0): + """Return stereographic scale factor in distance units. + + This is analogous to :func:`_stereo_distortion` but scaled by + ``pi/180`` to convert degrees to radians. When cartopy is + available, it reuses :class:`StereoProjection` and otherwise + falls back to the legacy implementation. + """ + + try: + from .region import _get_stereo_projection + + proj = _get_stereo_projection(scale_factor=k0) + if proj is not None: + return proj.get_scale_factor(lat) * np.pi / 180.0 + except Exception: + pass + lrad = np.radians(lat) - # Calculate the scale factor for the stereographic projection - res = 2 / (1 + np.sin(lrad)) / 180 * np.pi + res = 2 * k0 / (1 + np.sin(lrad)) / 180 * np.pi return res -def _generate_initial_points(min_edge_length, geps, bbox, fh, fd, pfix, stereo=False): - """Create initial distribution in bounding box (equilateral triangles)""" +def _generate_initial_points( + min_edge_length, geps, bbox, fh, fd, pfix, stereo=False, k0=1.0 +): + """Create initial distribution in bounding box (equilateral triangles). + + When ``stereo`` is True, points are generated in lon/lat and + then projected to stereographic coordinates. The local spacing is + adjusted by the stereographic distortion factor using the + configurable reference scale factor ``k0``. + """ + if stereo: bbox = np.array([[-180, 180], [-89, 89]]) p = np.mgrid[ @@ -1039,7 +1109,9 @@ def _generate_initial_points(min_edge_length, geps, bbox, fh, fd, pfix, stereo=F p0 = p.reshape(2, -1).T x, y = to_stereo(p0[:, 0], p0[:, 1]) p = np.asarray([x, y]).T - r0 = fh(to_lat_lon(p[:, 0], p[:, 1])) * _stereo_distortion(p0[:, 1]) + r0 = fh(to_lat_lon(p[:, 0], p[:, 1])) * _stereo_distortion( + p0[:, 1], k0=k0 + ) else: p = p.reshape(2, -1).T r0 = fh(p) diff --git a/oceanmesh/projections.py b/oceanmesh/projections.py new file mode 100644 index 0000000..16c1a40 --- /dev/null +++ b/oceanmesh/projections.py @@ -0,0 +1,165 @@ +"""Projection utilities for global stereographic meshing. + +This module provides a small wrapper around cartopy's stereographic +projection for use in global meshes. The primary entry point is the +:class:`StereoProjection` class, which exposes methods for converting +between geographic (longitude/latitude) coordinates and a polar +stereographic Cartesian system. + +Cartopy is treated as an optional dependency. The :data:`CARTOPY_AVAILABLE` +flag can be used to check availability at runtime, and the +:func:`check_cartopy_available` helper raises a clear, informative error +when global stereographic meshing is requested without cartopy being +installed. + +The design of this module and the choice of defaults follow the +"global stereographic meshing" discussion from PR #87, with a +north-polar stereographic projection centred on 0° longitude and +true scale at the pole. +""" + +from __future__ import annotations + +from typing import Tuple + +import numpy as np +from pyproj import CRS, Transformer + +try: # pragma: no cover - exercised indirectly in integration tests + import cartopy.crs as ccrs + + CARTOPY_AVAILABLE: bool = True +except Exception: # pragma: no cover - when cartopy is not installed + ccrs = None # type: ignore[assignment] + CARTOPY_AVAILABLE = False + + +class StereoProjection: + """North-polar stereographic projection helper. + + This class wraps cartopy's :class:`~cartopy.crs.NorthPolarStereo` + projection and exposes helpers for forward and inverse coordinate + transforms, as well as a convenience method for computing the local + scale factor of the projection. + + Parameters + ---------- + scale_factor: + The reference scale factor :math:`k_0` used in the analytic + expression for the stereographic scale distortion. + For a standard north-polar stereographic projection this is + typically 1.0. + """ + + def __init__(self, scale_factor: float = 1.0) -> None: + if not CARTOPY_AVAILABLE: # Defensive check; public callers + check_cartopy_available() # will usually call this first. + + # Geographic CRS for lon/lat coordinates. + self._pc = CRS.from_epsg(4326) + # North-polar stereographic CRS used for global meshes. Use a + # definition aligned with cartopy's NorthPolarStereo defaults. + self._stereo = CRS.from_proj4( + "+proj=stere +lat_0=90 +lat_ts=90 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" + ) + # pyproj transformer for forward/inverse operations. + self._fwd = Transformer.from_crs(self._pc, self._stereo, always_xy=True) + self._inv = Transformer.from_crs(self._stereo, self._pc, always_xy=True) + self.scale_factor = float(scale_factor) + + @property + def crs_geographic(self): + """Return the underlying geographic CRS (lon/lat).""" + + return self._pc + + @property + def crs_stereo(self): + """Return the underlying stereographic CRS.""" + + return self._stereo + + def to_stereo(self, lon: np.ndarray, lat: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + """Project longitude/latitude to stereographic x/y. + + Parameters + ---------- + lon, lat: + Arrays of longitudes and latitudes in degrees. + + Returns + ------- + x, y: + Arrays of stereographic coordinates in the projection's + native units (typically metres). + """ + + x, y = self._fwd.transform(lon, lat) + return x, y + + def to_lat_lon(self, x: np.ndarray, y: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + """Inverse transform from stereographic x/y to lon/lat. + + Parameters + ---------- + x, y: + Arrays of stereographic coordinates in the projection's + native units (typically metres). + + Returns + ------- + lon, lat: + Arrays of longitudes and latitudes in degrees. + """ + + lon, lat = self._inv.transform(x, y) + return lon, lat + + def get_scale_factor(self, lat: np.ndarray) -> np.ndarray: + """Return the stereographic scale factor at the given latitude. + + The scale factor :math:`k(\phi)` for a polar stereographic + projection is given by + + .. math:: + + k(\phi) = \frac{2 k_0}{1 + \sin \phi}, + + where :math:`k_0` is a reference scale factor. Here ``lat`` is + given in degrees and internally converted to radians. + + Parameters + ---------- + lat: + Latitude or array of latitudes in degrees. + + Returns + ------- + k: + The corresponding scale factor(s). + """ + + lat_rad = np.deg2rad(lat) + return 2.0 * self.scale_factor / (1.0 + np.sin(lat_rad)) + + +def check_cartopy_available() -> None: + """Raise an informative error if cartopy is not installed. + + This helper is intended to be called from high-level entry points + such as :class:`oceanmesh.geodata.Shoreline` when ``stereo=True`` + is requested. It ensures that users receive a clear error message + explaining how to enable the optional global meshing dependencies + instead of encountering opaque import or runtime failures later in + the workflow. + """ + + if not CARTOPY_AVAILABLE: + msg = ( + "Global stereographic meshing (stereo=True) requires cartopy. " + "Install it with: pip install oceanmesh[global] or pip install cartopy" + ) + raise ImportError(msg) + + +__all__ = ["StereoProjection", "check_cartopy_available", "CARTOPY_AVAILABLE"] diff --git a/oceanmesh/region.py b/oceanmesh/region.py index 221d628..1b75eb0 100644 --- a/oceanmesh/region.py +++ b/oceanmesh/region.py @@ -1,7 +1,9 @@ import numpy as np from pyproj import CRS, Transformer -__all__ = ["Region", "warp_coordinates"] +from .projections import CARTOPY_AVAILABLE, StereoProjection + +__all__ = ["Region", "warp_coordinates", "to_stereo", "to_lat_lon", "stereo_to_3d", "to_3d"] # ---- Helpers for CRS and bbox validation in multiscale meshes ---- @@ -104,6 +106,27 @@ def warp_coordinates(points, src_crs, dst_crs): return np.asarray(points).T +_STEREO_PROJ_CACHE = {} + + +def _get_stereo_projection(scale_factor=1.0): + """Return a cached :class:`StereoProjection` for a given k0. + + Parameters + ---------- + scale_factor: + Reference stereographic scale factor :math:`k_0` used in the + distortion formula. Projections are cached per distinct k0 + value when cartopy/pyproj are available. + """ + + global _STEREO_PROJ_CACHE + k0_key = float(scale_factor) + if k0_key not in _STEREO_PROJ_CACHE and CARTOPY_AVAILABLE: + _STEREO_PROJ_CACHE[k0_key] = StereoProjection(scale_factor=k0_key) + return _STEREO_PROJ_CACHE.get(k0_key, None) + + def stereo_to_3d(u, v, R=1): # to 3D # c=4*R**2/(u**2+v**2+4*R**2) @@ -120,6 +143,20 @@ def stereo_to_3d(u, v, R=1): def to_lat_lon(x, y, z=None, R=1): + """Convert stereographic or 3D Cartesian coordinates to lon/lat. + + When cartopy is available and ``z is None`` with the default + radius ``R == 1``, this uses the cartopy-based + :class:`StereoProjection` inverse transform. Otherwise it falls + back to the legacy analytic formulas for backward compatibility. + """ + + if z is None and R == 1 and CARTOPY_AVAILABLE: + proj = _get_stereo_projection() + if proj is not None: + lon, lat = proj.to_lat_lon(np.asarray(x), np.asarray(y)) + return lon, lat + if z is None: x, y, z = stereo_to_3d(x, y, R=R) @@ -150,11 +187,23 @@ def to_3d(x, y, R=1): def to_stereo(x, y, R=1): + """Project lon/lat to stereographic coordinates. + + When cartopy is available and the default radius ``R == 1`` is + used, this delegates to the cartopy-based :class:`StereoProjection` + helper. Otherwise, it falls back to the legacy analytic formulas + based on a unit sphere, preserving historical behaviour. + """ + + if R == 1 and CARTOPY_AVAILABLE: + proj = _get_stereo_projection() + if proj is not None: + u, v = proj.to_stereo(np.asarray(x), np.asarray(y)) + return u, v + kx, ky, kz = to_3d(x, y, R) - # to 2D in stereo - # u = 2*R*kx/(R+kz) - # v = 2*R*ky/(R+kz) + # to 2D in stereo (legacy formulation) u = -kx / (R + kz) v = -ky / (R + kz) diff --git a/oceanmesh/signed_distance_function.py b/oceanmesh/signed_distance_function.py index 88447a1..15cc34e 100644 --- a/oceanmesh/signed_distance_function.py +++ b/oceanmesh/signed_distance_function.py @@ -116,13 +116,22 @@ def _plot(geo, filename=None, samples=100000): class Domain: - def __init__(self, bbox, func, covering=None, crs=None, stereo=False): + def __init__( + self, + bbox, + func, + covering=None, + crs=None, + stereo=False, + scale_factor=1.0, + ): self.bbox = bbox self.domain = func self.covering = covering # Optional projection metadata for validation/rich workflows self.crs = crs self.stereo = bool(stereo) + self.scale_factor = float(scale_factor) def eval(self, x): return self.domain(x) @@ -150,7 +159,11 @@ def __init__(self, domains): crs_vals = [d.crs for d in domains if getattr(d, "crs", None) is not None] crs = crs_vals[0] if len(crs_vals) > 0 else None stereo = any(getattr(d, "stereo", False) for d in domains) - super().__init__(bbox, domains, crs=crs, stereo=stereo) + scale_factor_vals = [getattr(d, "scale_factor", 1.0) for d in domains] + scale_factor = max(scale_factor_vals) if scale_factor_vals else 1.0 + super().__init__( + bbox, domains, crs=crs, stereo=stereo, scale_factor=scale_factor + ) def eval(self, x): d = [d.eval(x) for d in self.domain] @@ -163,7 +176,11 @@ def __init__(self, domains): crs_vals = [d.crs for d in domains if getattr(d, "crs", None) is not None] crs = crs_vals[0] if len(crs_vals) > 0 else None stereo = any(getattr(d, "stereo", False) for d in domains) - super().__init__(bbox, domains, crs=crs, stereo=stereo) + scale_factor_vals = [getattr(d, "scale_factor", 1.0) for d in domains] + scale_factor = max(scale_factor_vals) if scale_factor_vals else 1.0 + super().__init__( + bbox, domains, crs=crs, stereo=stereo, scale_factor=scale_factor + ) def eval(self, x): d = [d.eval(x) for d in self.domain] @@ -176,7 +193,11 @@ def __init__(self, domains): crs_vals = [d.crs for d in domains if getattr(d, "crs", None) is not None] crs = crs_vals[0] if len(crs_vals) > 0 else None stereo = any(getattr(d, "stereo", False) for d in domains) - super().__init__(bbox, domains, crs=crs, stereo=stereo) + scale_factor_vals = [getattr(d, "scale_factor", 1.0) for d in domains] + scale_factor = max(scale_factor_vals) if scale_factor_vals else 1.0 + super().__init__( + bbox, domains, crs=crs, stereo=stereo, scale_factor=scale_factor + ) def eval(self, x): return np.maximum.reduce( @@ -269,10 +290,19 @@ def func_covering(x): dist = ((-1) ** (in_boubox)) * d return dist - # Attach CRS and stereo metadata from the shoreline for downstream validation + # Attach CRS, stereo, and stereographic scale factor metadata from + # the shoreline for downstream validation and projection logic. crs = getattr(shoreline, "crs", None) stereo = getattr(shoreline, "stereo", False) - return Domain(shoreline.bbox, func, covering=func_covering, crs=crs, stereo=stereo) + scale_factor = getattr(shoreline, "scale_factor", 1.0) + return Domain( + shoreline.bbox, + func, + covering=func_covering, + crs=crs, + stereo=stereo, + scale_factor=scale_factor, + ) def _create_boubox(bbox): @@ -312,15 +342,16 @@ def multiscale_signed_distance_function(signed_distance_functions): nests = [] for i, sdf in enumerate(signed_distance_functions): # set eval method to covering - tmp = [ - Domain( - s.bbox, - s.covering, - crs=getattr(s, "crs", None), - stereo=getattr(s, "stereo", False), - ) - for s in signed_distance_functions[i + 1 :] - ] + tmp = [ + Domain( + s.bbox, + s.covering, + crs=getattr(s, "crs", None), + stereo=getattr(s, "stereo", False), + scale_factor=getattr(s, "scale_factor", 1.0), + ) + for s in signed_distance_functions[i + 1 :] + ] nests.append(Difference([sdf, *tmp])) union = Union(nests) diff --git a/setup.cfg b/setup.cfg index 88dba70..b15e489 100644 --- a/setup.cfg +++ b/setup.cfg @@ -37,6 +37,9 @@ install_requires = scikit-image python_requires = >=3.10 +[options.extras_require] +global = cartopy + [versioneer] VCS = git style = pep440 diff --git a/tests/test_global_stereo.py b/tests/test_global_stereo.py index cfd7c8e..767a4b3 100644 --- a/tests/test_global_stereo.py +++ b/tests/test_global_stereo.py @@ -1,6 +1,7 @@ import os import oceanmesh as om +from oceanmesh.projections import CARTOPY_AVAILABLE # Note: global_stereo.shp has been generated using global_tag in pyposeidon # https://github.com/ec-jrc/pyPoseidon/blob/9cfd3bbf5598c810004def83b1f43dc5149addd0/pyposeidon/boundary.py#L452 @@ -20,7 +21,9 @@ def test_global_stereo(): min_edge_length = 0.5 # minimum mesh size in domain in meters max_edge_length = 2 # maximum mesh size in domain in meters shoreline = om.Shoreline(fname, extent.bbox, min_edge_length) + assert shoreline.scale_factor == 1.0 sdf = om.signed_distance_function(shoreline) + assert getattr(sdf, "scale_factor", 1.0) == 1.0 edge_length0 = om.distance_sizing_function(shoreline, rate=0.11) edge_length1 = om.feature_sizing_function( shoreline, @@ -58,3 +61,68 @@ def test_global_stereo(): ax.triplot(points[:, 0], points[:, 1], cells, color="gray", linewidth=0.5) shoreline_stereo.plot(ax=ax) + + # Smoke-check which projection backend is in use so failures + # can be interpreted in CI logs. + backend = "cartopy" if CARTOPY_AVAILABLE else "hardcoded" + print(f"Testing global stereographic mesh with {backend} projections") + + +def test_global_stereo_custom_k0(): + """Smoke-test configurable stereographic scale factor k0. + + Uses a slightly different reference scale factor to ensure that + Shoreline -> Domain -> mesh pipeline correctly propagates the + configuration without crashing. Exact numerical differences in + the output mesh are not asserted here. + """ + + EPSG = 4326 + bbox = (-180.00, 180.00, -89.00, 90.00) + extent = om.Region(extent=bbox, crs=4326) + + min_edge_length = 0.5 + max_edge_length = 2 + k0 = 0.994 + + shoreline = om.Shoreline( + fname, + extent.bbox, + min_edge_length, + scale_factor=k0, + ) + assert shoreline.scale_factor == k0 + + sdf = om.signed_distance_function(shoreline) + assert getattr(sdf, "scale_factor", None) == k0 + + edge_length0 = om.distance_sizing_function(shoreline, rate=0.11) + edge_length1 = om.feature_sizing_function( + shoreline, + sdf, + min_edge_length=min_edge_length, + max_edge_length=max_edge_length, + crs=EPSG, + ) + + edge_length = om.compute_minimum([edge_length0, edge_length1]) + edge_length = om.enforce_mesh_gradation( + edge_length, + gradation=0.09, + stereo=True, + scale_factor=k0, + ) + + shoreline_stereo = om.Shoreline( + fname2, + extent.bbox, + min_edge_length, + stereo=True, + scale_factor=k0, + ) + domain = om.signed_distance_function(shoreline_stereo) + assert getattr(domain, "scale_factor", None) == k0 + + points, cells = om.generate_mesh(domain, edge_length, stereo=True, max_iter=50) + assert points.shape[0] > 0 + assert cells.shape[0] > 0 From 84c18096077696a77eaf6a8c460825bc56fe2b6b Mon Sep 17 00:00:00 2001 From: Keith Roberts Date: Mon, 1 Dec 2025 21:14:37 -0500 Subject: [PATCH 2/4] linting --- oceanmesh/__init__.py | 2 +- oceanmesh/mesh_generator.py | 8 +++----- oceanmesh/projections.py | 21 +++++++++++---------- oceanmesh/region.py | 10 ++++++++-- oceanmesh/signed_distance_function.py | 20 ++++++++++---------- setup.cfg | 4 ++++ 6 files changed, 37 insertions(+), 28 deletions(-) diff --git a/oceanmesh/__init__.py b/oceanmesh/__init__.py index 3f5c19f..29df811 100644 --- a/oceanmesh/__init__.py +++ b/oceanmesh/__init__.py @@ -127,7 +127,7 @@ __version__ = _version.get_versions()["version"] try: # Optional global-stereo helpers (import may fail generically) - from .projections import StereoProjection, CARTOPY_AVAILABLE + from .projections import CARTOPY_AVAILABLE __all__.extend(["StereoProjection", "CARTOPY_AVAILABLE"]) except ImportError: diff --git a/oceanmesh/mesh_generator.py b/oceanmesh/mesh_generator.py index 4595d11..4c6e068 100644 --- a/oceanmesh/mesh_generator.py +++ b/oceanmesh/mesh_generator.py @@ -1030,7 +1030,7 @@ def _deps_vec(i): def _stereo_distortion(lat, k0=1.0): - """Return stereographic scale factor at latitude. + r"""Return stereographic scale factor at latitude. When cartopy is available, this uses :class:`oceanmesh.projections.StereoProjection` to compute the @@ -1061,7 +1061,7 @@ def _stereo_distortion(lat, k0=1.0): def _stereo_distortion_dist(lat, k0=1.0): - """Return stereographic scale factor in distance units. + r"""Return stereographic scale factor in distance units. This is analogous to :func:`_stereo_distortion` but scaled by ``pi/180`` to convert degrees to radians. When cartopy is @@ -1109,9 +1109,7 @@ def _generate_initial_points( p0 = p.reshape(2, -1).T x, y = to_stereo(p0[:, 0], p0[:, 1]) p = np.asarray([x, y]).T - r0 = fh(to_lat_lon(p[:, 0], p[:, 1])) * _stereo_distortion( - p0[:, 1], k0=k0 - ) + r0 = fh(to_lat_lon(p[:, 0], p[:, 1])) * _stereo_distortion(p0[:, 1], k0=k0) else: p = p.reshape(2, -1).T r0 = fh(p) diff --git a/oceanmesh/projections.py b/oceanmesh/projections.py index 16c1a40..f552f75 100644 --- a/oceanmesh/projections.py +++ b/oceanmesh/projections.py @@ -35,20 +35,19 @@ class StereoProjection: - """North-polar stereographic projection helper. + r"""North-polar stereographic projection helper. - This class wraps cartopy's :class:`~cartopy.crs.NorthPolarStereo` - projection and exposes helpers for forward and inverse coordinate - transforms, as well as a convenience method for computing the local - scale factor of the projection. + This class provides forward (lon/lat -> x/y) and inverse (x/y -> lon/lat) + transforms for a north-polar stereographic CRS, along with a helper for + computing the local scale factor :math:`k(\phi)` using a configurable + reference scale factor :math:`k_0`. Parameters ---------- scale_factor: The reference scale factor :math:`k_0` used in the analytic - expression for the stereographic scale distortion. - For a standard north-polar stereographic projection this is - typically 1.0. + expression for the stereographic scale distortion. For a standard + north-polar stereographic projection this is typically 1.0. """ def __init__(self, scale_factor: float = 1.0) -> None: @@ -79,7 +78,9 @@ def crs_stereo(self): return self._stereo - def to_stereo(self, lon: np.ndarray, lat: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + def to_stereo( + self, lon: np.ndarray, lat: np.ndarray + ) -> Tuple[np.ndarray, np.ndarray]: """Project longitude/latitude to stereographic x/y. Parameters @@ -116,7 +117,7 @@ def to_lat_lon(self, x: np.ndarray, y: np.ndarray) -> Tuple[np.ndarray, np.ndarr return lon, lat def get_scale_factor(self, lat: np.ndarray) -> np.ndarray: - """Return the stereographic scale factor at the given latitude. + r"""Return the stereographic scale factor at the given latitude. The scale factor :math:`k(\phi)` for a polar stereographic projection is given by diff --git a/oceanmesh/region.py b/oceanmesh/region.py index 1b75eb0..baf1ab7 100644 --- a/oceanmesh/region.py +++ b/oceanmesh/region.py @@ -3,7 +3,14 @@ from .projections import CARTOPY_AVAILABLE, StereoProjection -__all__ = ["Region", "warp_coordinates", "to_stereo", "to_lat_lon", "stereo_to_3d", "to_3d"] +__all__ = [ + "Region", + "warp_coordinates", + "to_stereo", + "to_lat_lon", + "stereo_to_3d", + "to_3d", +] # ---- Helpers for CRS and bbox validation in multiscale meshes ---- @@ -120,7 +127,6 @@ def _get_stereo_projection(scale_factor=1.0): value when cartopy/pyproj are available. """ - global _STEREO_PROJ_CACHE k0_key = float(scale_factor) if k0_key not in _STEREO_PROJ_CACHE and CARTOPY_AVAILABLE: _STEREO_PROJ_CACHE[k0_key] = StereoProjection(scale_factor=k0_key) diff --git a/oceanmesh/signed_distance_function.py b/oceanmesh/signed_distance_function.py index 15cc34e..a143adb 100644 --- a/oceanmesh/signed_distance_function.py +++ b/oceanmesh/signed_distance_function.py @@ -342,16 +342,16 @@ def multiscale_signed_distance_function(signed_distance_functions): nests = [] for i, sdf in enumerate(signed_distance_functions): # set eval method to covering - tmp = [ - Domain( - s.bbox, - s.covering, - crs=getattr(s, "crs", None), - stereo=getattr(s, "stereo", False), - scale_factor=getattr(s, "scale_factor", 1.0), - ) - for s in signed_distance_functions[i + 1 :] - ] + tmp = [ + Domain( + s.bbox, + s.covering, + crs=getattr(s, "crs", None), + stereo=getattr(s, "stereo", False), + scale_factor=getattr(s, "scale_factor", 1.0), + ) + for s in signed_distance_functions[i + 1 :] + ] nests.append(Difference([sdf, *tmp])) union = Union(nests) diff --git a/setup.cfg b/setup.cfg index b15e489..cb57d7e 100644 --- a/setup.cfg +++ b/setup.cfg @@ -38,6 +38,10 @@ install_requires = python_requires = >=3.10 [options.extras_require] +fast = + cython + numpy + scipy global = cartopy [versioneer] From 0141f10783ec0f103cf5883386a0a0ddd605cfb4 Mon Sep 17 00:00:00 2001 From: Keith Roberts Date: Sun, 7 Dec 2025 22:11:16 -0500 Subject: [PATCH 3/4] Stereo-aware global meshing: configurable k0, flexible CRS, and docs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Add configurable stereographic scale_factor (k0) plumbed through Shoreline → Domain → meshing and gradient enforcement, with Cartopy-backed and analytic fallbacks. Make global and multiscale CRS handling more permissive and CRS-aware, including mixed CRS validation and global bbox detection. Improve polar gradient limiting using stereographic distortion correction and add/adjust tests for stereo, k0, and mixed-CRS multiscale cases. Expand documentation (module docstrings and README) to describe stereographic workflow, k0 semantics, and Cartopy usage. --- README.md | 99 +++++++++- oceanmesh/edgefx.py | 49 +++-- oceanmesh/geodata.py | 173 ++++++++++++++---- oceanmesh/mesh_generator.py | 153 ++++++++++++++-- oceanmesh/region.py | 222 ++++++++++++++++++++--- tests/test_global_regional_multiscale.py | 74 ++++++++ tests/test_global_stereo.py | 129 +++++++++++++ 7 files changed, 804 insertions(+), 95 deletions(-) diff --git a/README.md b/README.md index 2f08b78..6b8fc27 100644 --- a/README.md +++ b/README.md @@ -428,6 +428,67 @@ points, cells = om.generate_multiscale_mesh([sdf1, sdf2], [el1, el2]) Global meshes are defined in WGS84 but meshed in a stereographic projection. Regional refinement can be added as additional domains. +#### Understanding Stereographic Projections for Global Meshes + +Global ocean meshes require special handling due to the distortion +inherent in projecting a sphere onto a 2D plane. OceanMesh uses a +**north-polar stereographic projection** for global meshing, which: + +- Preserves angles (conformal projection) +- Minimises distortion near the poles +- Allows efficient 2D mesh generation algorithms to work on projected + coordinates + +**Cartopy Integration (Optional Dependency)** + +For accurate stereographic projections, OceanMesh can use +[cartopy](https://scitools.org.uk/cartopy/docs/latest/)'s +``NorthPolarStereo`` coordinate reference system. Install cartopy for +global meshing via the ``global`` extra: + +```bash +pip install oceanmesh[global] +``` + +or manually: + +```bash +pip install cartopy +``` + +If cartopy is not installed, OceanMesh falls back to legacy analytic +formulas (less accurate but functional). + +**Scale Factor (k0) Configuration** + +The stereographic projection uses a **scale factor** :math:`k_0` that +controls distortion at different latitudes. The local scale factor at +latitude :math:`\phi` is + +.. math:: + + k(\phi) = \frac{2 k_0}{1 + \sin \phi}. + +Common :math:`k_0` values: + +- :math:`k_0 = 1.0` (default): standard stereographic with true scale + at the pole. +- :math:`k_0 = 0.994`: used in EPSG:3413 (Arctic) and EPSG:3031 + (Antarctic) to minimise distortion at standard parallels (~70°N/S). + +The scale factor affects mesh sizing: smaller :math:`k_0` reduces +element sizes near the poles, larger :math:`k_0` increases them. + +References +~~~~~~~~~~ + +- [GitHub PR #87](https://github.com/CHLNDDEV/oceanmesh/pull/87): + discussion of scale-factor improvements. +- [Seamsh stereographic example](https://git.immc.ucl.ac.be/jlambrechts/seamsh/-/blob/ca380b59fe4d2ea57ffbf08a7b0c70bdf7df1afb/examples/6-stereographics.py#L55): + similar implementation in seamsh. +- Snyder, J.P. (1987). *Map Projections - A Working Manual*. USGS + Professional Paper 1395. + ```python @@ -442,7 +503,17 @@ 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) + +# Use k0 = 0.994 for reduced polar distortion (similar to +# EPSG:3413/3031 in polar stereographic grids). +k0 = 0.994 + +shoreline_global_latlon = om.Shoreline( + fname_global_latlon, + global_region, + 1.0, + scale_factor=k0, +) sdf_global_latlon = om.signed_distance_function(shoreline_global_latlon) edge_global = om.enforce_mesh_gradation( om.compute_minimum([ @@ -451,6 +522,7 @@ edge_global = om.enforce_mesh_gradation( ]), gradation=0.15, stereo=True, + scale_factor=k0, ) aus_region = om.Region(extent=(110.0, 160.0, -45.0, -10.0), crs=4326) @@ -464,7 +536,13 @@ edge_regional = om.enforce_mesh_gradation( gradation=0.12, ) -shoreline_global_stereo = om.Shoreline(fname_global_stereo, global_region, 1.0, stereo=True) +shoreline_global_stereo = om.Shoreline( + fname_global_stereo, + global_region, + 1.0, + stereo=True, + scale_factor=k0, +) sdf_global_stereo = om.signed_distance_function(shoreline_global_stereo) points, cells = om.generate_multiscale_mesh( @@ -485,13 +563,26 @@ points, cells = om.generate_multiscale_mesh( # Global mesh generation only (stereographic meshing) import oceanmesh as om from oceanmesh.region import to_lat_lon + fname = "tests/global/global_latlon.shp" fname2 = "tests/global/global_stereo.shp" region = om.Region(extent=(-180.00, 180.00, -89.00, 90.00), crs=4326) -shore = om.Shoreline(fname, region.bbox, 0.5) + +# Use k0 = 0.994 to mirror common polar stereographic practice. +k0 = 0.994 + +shore = om.Shoreline(fname, region.bbox, 0.5, scale_factor=k0) edge = om.distance_sizing_function(shore, rate=0.11) -domain = om.signed_distance_function(om.Shoreline(fname2, region.bbox, 0.5, stereo=True)) +edge = om.enforce_mesh_gradation(edge, gradation=0.15, stereo=True, scale_factor=k0) + +domain = om.signed_distance_function( + om.Shoreline(fname2, region.bbox, 0.5, stereo=True, scale_factor=k0) +) points, cells = om.generate_mesh(domain, edge, stereo=True, max_iter=100) + +# Points are in stereographic coordinates; convert back to lon/lat for +# post-processing or export. +lon, lat = to_lat_lon(points[:, 0], points[:, 1]) ``` [Back to top](#table-of-contents) diff --git a/oceanmesh/edgefx.py b/oceanmesh/edgefx.py index 1bb6a91..cb1b616 100644 --- a/oceanmesh/edgefx.py +++ b/oceanmesh/edgefx.py @@ -129,14 +129,14 @@ def enforce_mesh_gradation( from .projections import CARTOPY_AVAILABLE from .region import _get_stereo_projection + backend = "cartopy" if CARTOPY_AVAILABLE else "hardcoded" logger.info( - "Global mesh: fixing gradient on the north pole (using %s projection, k0=%s)...", - "cartopy" if CARTOPY_AVAILABLE else "hardcoded", + "Global mesh: fixing gradient on the north pole (backend=%s, k0=%s, distortion_correction=%s)...", + backend, scale_factor, + "enabled" if CARTOPY_AVAILABLE else "disabled", ) - # Obtain a stereographic projection configured with the - # caller-provided scale_factor so that k0 meaningfully affects - # the polar correction. + proj = _get_stereo_projection(scale_factor=scale_factor) # Choose a reasonable physical radius (in metres) around the @@ -156,18 +156,45 @@ def enforce_mesh_gradation( if proj is not None: ulon, vlat = proj.to_lat_lon(us.ravel(), vs.ravel()) + # Spatially-varying stereographic scale factor k(phi) on + # the intermediate polar grid. + k_polar = proj.get_scale_factor(vlat) else: # Fall back to module-level transform if projection could - # not be constructed. - ulon, vlat = to_lat_lon(us.ravel(), vs.ravel()) + # not be constructed. Normalise metre-based coordinates by + # Earth radius so that to_lat_lon operates on the + # non-dimensional unit-sphere coordinates expected by the + # legacy stereographic formulas used when CARTOPY_AVAILABLE + # is False. This keeps the polar patch consistent with + # to_stereo/to_lat_lon. + R_earth = 6_371_000.0 + us_unit = us / R_earth + vs_unit = vs / R_earth + ulon, vlat = to_lat_lon(us_unit.ravel(), vs_unit.ravel()) + k_polar = None + utmp = grid.eval((ulon, vlat)) utmp = np.reshape(utmp, us.shape) - szs = utmp.shape + + # Distortion-aware gradient limiting: operate on mesh sizes + # scaled by the local stereographic factor so the limiter + # acts in approximately physical space. + if k_polar is not None: + k_polar_grid = np.reshape(k_polar, utmp.shape) + utmp_work = utmp * k_polar_grid + else: + utmp_work = utmp + + szs = utmp_work.shape szs = (szs[0], szs[1], 1) - # we choose an excessively large number for the gradiation = 10 - # this is merely to fix the north pole gradient - vtmp = gradient_limit([*szs], dx_stereo, 10, 10000, utmp.flatten("F")) + vtmp = gradient_limit([*szs], dx_stereo, 10, 10000, utmp_work.flatten("F")) vtmp = np.reshape(vtmp, (szs[0], szs[1]), "F") + + # Inverse correction: return to geographic sizing by dividing + # out the stereographic distortion. + if k_polar is not None: + vtmp /= k_polar_grid + # construct stereo interpolating function in projection metres grid_stereo = Grid( bbox=(-r_max, r_max, -r_max, r_max), diff --git a/oceanmesh/geodata.py b/oceanmesh/geodata.py index 15b8ebc..442ef26 100644 --- a/oceanmesh/geodata.py +++ b/oceanmesh/geodata.py @@ -533,7 +533,7 @@ def remove_dup(arr: np.ndarray): class Shoreline(Region): - """ + r""" The shoreline class extends :class:`Region` to store data that is later used to create signed distance functions to represent irregular shoreline geometries. This data @@ -565,6 +565,32 @@ class Shoreline(Region): global world meshes (with EPSG:4326 inputs), and False for regional meshes. This flag is retained on the Shoreline instance and propagated to downstream Domain objects for validation in generate_multiscale_mesh. + scale_factor : float, optional + Reference stereographic scale factor :math:`k_0` used in distortion + calculations for global meshes when ``stereo=True``. Default is 1.0 + (standard north-polar stereographic). The scale factor affects mesh + sizing near the poles via the formula + + .. math:: + + k(\phi) = \frac{2 k_0}{1 + \sin \phi}, + + where :math:`\phi` is latitude. A value of ``k0 = 0.994`` is commonly + used in polar stereographic projections (e.g., EPSG:3413 for the + Arctic and EPSG:3031 for the Antarctic) to minimise distortion at + standard parallels. + + This parameter is propagated through :class:`Domain` objects to mesh + generation and gradient enforcement functions. It is only relevant + when ``stereo=True``. + + References + ---------- + - GitHub PR #87: https://github.com/CHLNDDEV/oceanmesh/pull/87 + - Seamsh stereographic example: + https://git.immc.ucl.ac.be/jlambrechts/seamsh/-/blob/ca380b59fe4d2ea57ffbf08a7b0c70bdf7df1afb/examples/6-stereographics.py#L55 + - Snyder, J.P. (1987). *Map Projections - A Working Manual*. + USGS Professional Paper 1395. """ def __init__( @@ -582,39 +608,7 @@ def __init__( if isinstance(shp, str): shp = Path(shp) - # Determine bbox coordinates and CRS to use - crs_to_use = crs - - if isinstance(bbox, Region): - bbox_coords = bbox.bbox - region_crs = bbox.crs - if crs != "EPSG:4326": - logger.warning( - "Shoreline: Both a Region object and an explicit 'crs' were provided; " - "the Region's CRS will take precedence (crs=%s).", - region_crs, - ) - crs_to_use = region_crs - else: - bbox_coords = bbox - # Backward compatibility for tuple/array input; optional auto-detection - if isinstance(bbox_coords, tuple) and crs == "EPSG:4326": - looks_geo = _infer_crs_from_coordinates(bbox_coords) - if ( - not looks_geo - ): # Only inspect shapefile native CRS when clearly projected - try: - native = gpd.read_file(shp).crs - except Exception: - native = None - if (native is not None) and ( - CRS.from_user_input(native) != CRS.from_user_input("EPSG:4326") - ): - logger.info( - "Shoreline: bbox looks projected but 'crs' not specified; using shapefile's native CRS %s", - native, - ) - crs_to_use = native + crs_to_use, bbox_coords = self._determine_bbox_and_crs(shp, bbox, crs) # Build polygon representation and normalized bbox tuple if isinstance(bbox_coords, tuple): @@ -641,6 +635,12 @@ def __init__( self.stereo = bool(stereo) if scale_factor <= 0: raise ValueError("scale_factor must be > 0") + # Store k0 for the stereographic distortion formula + # + # k(\phi) = 2 k_0 / (1 + sin(\phi)) + # + # which is used to modulate mesh sizing as a function of + # latitude in global meshes (see PR #87 and Snyder, 1987). self.__scale_factor = float(scale_factor) self.inner = [] self.outer = [] @@ -649,9 +649,48 @@ def __init__( self.refinements = refinements self.minimum_area_mult = minimum_area_mult if self.stereo: - from oceanmesh.projections import check_cartopy_available + from oceanmesh.projections import CARTOPY_AVAILABLE, check_cartopy_available + from oceanmesh.region import is_crs_suitable_for_global, is_global_bbox + + # Stereo handling can fall back to legacy analytic + # formulas when cartopy is unavailable. Only require + # cartopy when the cartopy-backed projection helper is in + # use; otherwise log a warning and continue so existing + # hardcoded workflows remain supported. + if CARTOPY_AVAILABLE: + check_cartopy_available() + else: + logger.warning( + "Shoreline: cartopy is not available; falling back to legacy stereographic " + "formulas for stereo=True workflows. Consider installing cartopy for more " + "accurate global projections." + ) + + suitable, msg = is_crs_suitable_for_global(crs_to_use) + if suitable: + logger.info("Shoreline: global (stereo=True) domain CRS check: %s", msg) + else: + logger.warning( + "Shoreline: global (stereo=True) domain CRS may be unsuitable: %s", + msg, + ) + + # Additional warning when bbox looks global and CRS is not + # the canonical EPSG:4326 geographic CRS. + try: + crs_obj = CRS.from_user_input(crs_to_use) + epsg = crs_obj.to_epsg() + except Exception: + crs_obj = None + epsg = None - check_cartopy_available() + if is_global_bbox(bbox_tuple, crs_obj) and epsg != 4326: + logger.warning( + "Global domain detected with CRS %s (not EPSG:4326). Ensure this CRS is appropriate " + "for global-scale meshing. Stereographic transformations will use scale_factor=%s.", + crs_to_use, + scale_factor, + ) polys = self._read() @@ -675,6 +714,48 @@ def __init__( self.bbox, self.boubox, polys, self.h0 / 2, self.minimum_area_mult, stereo ) + @staticmethod + def _determine_bbox_and_crs(shp, bbox, crs): + """Return (crs_to_use, bbox_coords) applying legacy inference rules. + + This helper centralises the bbox/CRS selection logic + """ + + crs_to_use = crs + + if isinstance(bbox, Region): + bbox_coords = bbox.bbox + region_crs = bbox.crs + if crs != "EPSG:4326": + logger.warning( + "Shoreline: Both a Region object and an explicit 'crs' were provided; " + "the Region's CRS will take precedence (crs=%s).", + region_crs, + ) + crs_to_use = region_crs + else: + bbox_coords = bbox + # Backward compatibility for tuple/array input; optional auto-detection + if isinstance(bbox_coords, tuple) and crs == "EPSG:4326": + looks_geo = _infer_crs_from_coordinates(bbox_coords) + if ( + not looks_geo + ): # Only inspect shapefile native CRS when clearly projected + try: + native = gpd.read_file(shp).crs + except Exception: + native = None + if (native is not None) and ( + CRS.from_user_input(native) != CRS.from_user_input("EPSG:4326") + ): + logger.info( + "Shoreline: bbox looks projected but 'crs' not specified; using shapefile's native CRS %s", + native, + ) + crs_to_use = native + + return crs_to_use, bbox_coords + @property def shp(self): return self.__shp @@ -720,7 +801,25 @@ def h0(self, value): @property def scale_factor(self): - """Reference stereographic scale factor k0 for global meshes.""" + r"""Reference stereographic scale factor :math:`k_0` for global meshes. + + The scale factor :math:`k_0` (typically 1.0) is used in the + stereographic distortion formula + + .. math:: + + k(\phi) = \frac{2 k_0}{1 + \sin \phi}, + + to adjust mesh sizing at different latitudes. Common choices are + + * ``k0 = 1.0`` (default): standard north-polar stereographic + * ``k0 = 0.994``: used in EPSG:3413 (Arctic) and EPSG:3031 + (Antarctic) to reduce distortion at standard parallels + (roughly 70°N/S). + + See PR #87 for a discussion of scale-factor selection and how it + propagates through the Shoreline → Domain → mesh pipeline. + """ return self.__scale_factor diff --git a/oceanmesh/mesh_generator.py b/oceanmesh/mesh_generator.py index 4c6e068..4d3f075 100644 --- a/oceanmesh/mesh_generator.py +++ b/oceanmesh/mesh_generator.py @@ -358,7 +358,7 @@ def _validate_multiscale_domains(domains, edge_lengths): # noqa: C901 if ( dcrs is not None and CRS.from_user_input(dcrs).to_epsg() == 4326 - and is_global_bbox(d.bbox) + and is_global_bbox(d.bbox, getattr(d, "crs", None)) ): implicit_global_idx = i break @@ -383,17 +383,34 @@ def _validate_multiscale_domains(domains, edge_lengths): # noqa: C901 gcrs_str = get_crs_string(gcrs) except Exception: gcrs_str = str(gcrs) - # Explicitly require global CRS be EPSG:4326 + # Relaxed global CRS requirement: allow non-EPSG:4326 but + # emit guidance so users can verify suitability. try: parsed = CRS.from_user_input(gcrs) - if parsed.to_epsg() != 4326: - errors.append( - f"Global domain CRS {gcrs_str} must be EPSG:4326 for global+regional multiscale meshing." - ) + epsg = parsed.to_epsg() except Exception: - errors.append( - f"Global domain CRS '{gcrs_str}' could not be parsed; expected EPSG:4326." + parsed = None + epsg = None + + if parsed is None: + logger.warning( + "Global domain CRS '%s' could not be parsed; proceeding but multiscale validation " + "may be unreliable.", + gcrs_str, ) + else: + if epsg != 4326: + logger.warning( + "Global domain CRS %s is not EPSG:4326. This may work for global+regional multiscale " + "meshing if the CRS is appropriate for global coverage. Proceeding with validation.", + gcrs_str, + ) + if parsed.is_projected and is_global_bbox(global_domain.bbox, parsed): + logger.warning( + "Global domain uses projected CRS %s with global-like bbox. Ensure this projection is " + "suitable for global-scale meshing (e.g., stereographic).", + gcrs_str, + ) # Containment checks for regional domains for i, d in enumerate(domains[1:], start=1): # Perform containment in lat/lon if global stereographic domain provided its bbox in stereo space @@ -419,7 +436,12 @@ def _validate_multiscale_domains(domains, edge_lengths): # noqa: C901 f"Regional domain #{i} bbox {d_bbox} (lat/lon) not contained within global stereo bbox {g_bbox}." ) else: - if not bbox_contains(g_bbox, d_bbox): + if not bbox_contains( + g_bbox, + d_bbox, + getattr(global_domain, "crs", None), + getattr(d, "crs", None), + ): errors.append( f"Regional domain #{i} bbox {d_bbox} not contained within global bbox {g_bbox}." ) @@ -438,6 +460,13 @@ def _validate_multiscale_domains(domains, edge_lengths): # noqa: C901 ) if not ok_crs: errors.append(msg_crs) + else: + logger.info( + "Multiscale CRS compatibility: global=%s, regional=%s -> %s", + get_crs_string(getattr(global_domain, "crs", None)), + get_crs_string(getattr(d, "crs", None)), + msg_crs, + ) # Implicit global-like domain with EPSG:4326 but stereo=False mixing with different CRS if implicit_global_idx is not None: @@ -468,7 +497,8 @@ def _validate_multiscale_domains(domains, edge_lengths): # noqa: C901 # If CRS parsing fails, skip this implicit enforcement pass - # Edge length CRS matching + # Edge length CRS matching and informational logging for projected + # regional domains. for i, (d, el) in enumerate(zip(domains, edge_lengths)): if hasattr(el, "crs"): el_crs = getattr(el, "crs", None) @@ -486,6 +516,21 @@ def _validate_multiscale_domains(domains, edge_lengths): # noqa: C901 f"Edge length #{i} CRS could not be compared to domain CRS (el={get_crs_string(el_crs)}, domain={get_crs_string(d_crs)})." ) + # Projected CRS support for regional domains: log when projected + # CRS is used for refinement. + d_crs = getattr(d, "crs", None) + if (i > 0) and (d_crs is not None): + try: + d_parsed = CRS.from_user_input(d_crs) + except Exception: + d_parsed = None + if d_parsed is not None and d_parsed.is_projected: + logger.info( + "Regional domain #%d uses projected CRS %s, which is appropriate for regional refinement.", + i, + get_crs_string(d_crs), + ) + return len(errors) == 0, errors @@ -1032,15 +1077,52 @@ def _deps_vec(i): def _stereo_distortion(lat, k0=1.0): r"""Return stereographic scale factor at latitude. - When cartopy is available, this uses - :class:`oceanmesh.projections.StereoProjection` to compute the - scale factor via + Computes the local scale factor for a north-polar stereographic + projection, which quantifies how distances on the projected plane + relate to distances on the sphere at a given latitude. + + Parameters + ---------- + lat : float or array_like + Latitude in degrees (scalar or array). + k0 : float, optional + Reference scale factor at the projection origin (north pole). + Default is 1.0. Common values include + + * ``k0 = 1.0``: standard stereographic (true scale at the pole) + * ``k0 = 0.994``: used in EPSG:3413/3031 to minimise distortion + near ~70° latitude. + + Returns + ------- + k : float or array_like + Local scale factor at the given latitude(s). + + Notes + ----- + The scale factor formula for a north-polar stereographic projection + is + + .. math:: + + k(\phi) = \frac{2 k_0}{1 + \sin \phi}, - .. math:: k(\phi) = 2 k_0 / (1 + \sin \phi), + where :math:`\phi` is latitude and :math:`k_0` is the reference + scale factor. When cartopy is available, this function delegates to + :class:`oceanmesh.projections.StereoProjection` so that the scale + factors are computed in a CRS-aware manner; otherwise it falls back + to the analytic expression above based on a unit sphere. - with :math:`k_0 = 1`. Otherwise, it falls back to the historical - analytic expression based on a north-polar stereographic - projection of the unit sphere. + The value of :math:`k(\phi)` influences mesh sizing near the poles + when generating global meshes with ``stereo=True``. + + References + ---------- + - Snyder, J.P. (1987). *Map Projections - A Working Manual*. + USGS Professional Paper 1395. + - GitHub PR #87: https://github.com/CHLNDDEV/oceanmesh/pull/87 + - Seamsh stereographic example: + https://git.immc.ucl.ac.be/jlambrechts/seamsh/-/blob/ca380b59fe4d2ea57ffbf08a7b0c70bdf7df1afb/examples/6-stereographics.py#L55 """ try: @@ -1064,9 +1146,40 @@ def _stereo_distortion_dist(lat, k0=1.0): r"""Return stereographic scale factor in distance units. This is analogous to :func:`_stereo_distortion` but scaled by - ``pi/180`` to convert degrees to radians. When cartopy is - available, it reuses :class:`StereoProjection` and otherwise - falls back to the legacy implementation. + :math:`\pi/180` to convert degrees to radians, making it suitable + for distance-based calculations in mesh generation. + + Parameters + ---------- + lat : float or array_like + Latitude in degrees. + k0 : float, optional + Reference scale factor. Default is 1.0. + + Returns + ------- + k_dist : float or array_like + Scale factor in distance units (radians per degree). + + Notes + ----- + The formula is + + .. math:: + + k_{dist}(\phi) = \frac{2 k_0}{1 + \sin \phi} \cdot \frac{\pi}{180}, + + where :math:`\phi` is latitude. When cartopy is available this + function reuses :class:`StereoProjection.get_scale_factor` and + applies the radian conversion; otherwise it falls back to the + analytic expression above. + + This function is used internally to adjust edge lengths based on + stereographic distortion. + + See Also + -------- + _stereo_distortion : Scale factor without unit conversion. """ try: diff --git a/oceanmesh/region.py b/oceanmesh/region.py index baf1ab7..8db3363 100644 --- a/oceanmesh/region.py +++ b/oceanmesh/region.py @@ -1,3 +1,64 @@ +r"""Coordinate-system utilities and stereographic helpers for OceanMesh. + +This module centralises logic for working with geographic and +projected coordinate reference systems (CRS), including: + +* Lightweight wrappers around :mod:`pyproj` for CRS parsing and + coordinate warping. +* Convenience helpers for detecting "global-like" bounding boxes and + performing CRS-aware containment checks used in multiscale mesh + validation. +* Optional integration with the cartopy-backed + :class:`~oceanmesh.projections.StereoProjection` class for + north-polar stereographic projections. + +Stereographic projections for global meshes +------------------------------------------ + +For global ocean meshes, OceanMesh uses a north-polar stereographic +projection: geographic coordinates (lon, lat) are projected onto a +plane (u, v) tangent at the north pole. When cartopy is available, +this is implemented via its ``NorthPolarStereo`` CRS combined with a +``PlateCarree`` (lon/lat) CRS. The forward and inverse transforms are +exposed through :func:`to_stereo` and :func:`to_lat_lon`. + +The local scale factor :math:`k(\phi)` of a polar stereographic +projection, which describes how distances on the plane relate to +distances on the sphere at latitude :math:`\phi`, is given by + +.. math:: + + k(\phi) = \frac{2 k_0}{1 + \sin \phi}, + +where :math:`k_0` is a configurable reference scale factor (typically +1.0). A commonly used value in practical ocean and ice modelling is +``k0 = 0.994``, mirroring conventions in EPSG:3413 (Arctic) and +EPSG:3031 (Antarctic) to reduce distortion near standard parallels. + +When cartopy is not installed, OceanMesh falls back to legacy analytic +formulas based on a unit sphere. These formulas are consistent with +the mapping used in the original Fortran/C implementations and remain +useful for environments where cartopy is difficult to install, +albeit with less CRS metadata. + +Background and references +------------------------- + +The configurable scale-factor support (:math:`k_0`) and cartopy-based +projection backend were introduced in PR #87: + + https://github.com/CHLNDDEV/oceanmesh/pull/87 + +A closely related implementation can be found in the seamsh example: + + https://git.immc.ucl.ac.be/jlambrechts/seamsh/-/blob/ca380b59fe4d2ea57ffbf08a7b0c70bdf7df1afb/examples/6-stereographics.py#L55 + +For theoretical background and additional formulas, see + + Snyder, J.P. (1987). *Map Projections – A Working Manual*. + U.S. Geological Survey Professional Paper 1395. +""" + import numpy as np from pyproj import CRS, Transformer @@ -30,30 +91,73 @@ def get_crs_string(crs): return str(crs) -def is_global_bbox(bbox): +def is_global_bbox(bbox, crs=None): """Heuristic to determine if a bbox is 'global-like'. - Returns True when bbox covers most of the world extents. - Expected bbox format: (xmin, xmax, ymin, ymax) in degrees. + Parameters + ---------- + bbox : tuple + (xmin, xmax, ymin, ymax) in degrees or projection units. + crs : optional + CRS of the bbox. When provided and projected, looser thresholds + are used to decide whether the bbox is effectively global. """ + xmin, xmax, ymin, ymax = bbox lon_span = xmax - xmin lat_span = ymax - ymin + + if crs is not None: + try: + _crs = CRS.from_user_input(crs) + except Exception: + _crs = None + if _crs is not None and _crs.is_projected: + # For projected CRS we rely purely on span heuristics in + # the projected units. Large extents compared to Earth + # radius suggest global coverage. + span_x = abs(lon_span) + span_y = abs(lat_span) + threshold = 5_000_000.0 # ~5000 km + return bool(span_x >= threshold and span_y >= threshold) + lon_global = lon_span >= 300 or (xmin <= -180 and xmax >= 180) lat_global = lat_span >= 150 or (ymin <= -89 and ymax >= 89) return bool(lon_global and lat_global) -def bbox_contains(outer, inner): +def bbox_contains(outer, inner, outer_crs=None, inner_crs=None): """Return True if outer bbox fully contains inner bbox. - Both bboxes are (xmin, xmax, ymin, ymax). If the outer bbox is detected - as global using is_global_bbox, containment is assumed True. If the inner - bbox is global and outer is not, containment is False. + Both bboxes are (xmin, xmax, ymin, ymax). If CRS information is + provided and differs, the inner bbox is transformed into the outer + CRS before checking containment. For global-like outer bboxes, + containment is assumed True. If the inner bbox is global and outer + is not, containment is False. """ - if is_global_bbox(outer): + + if outer_crs is not None and inner_crs is not None: + try: + o = CRS.from_user_input(outer_crs) + i = CRS.from_user_input(inner_crs) + except Exception: + o = i = None + if o is not None and i is not None and not o.equals(i): + transformer = Transformer.from_crs(i, o, always_xy=True) + ixmin, ixmax, iymin, iymax = inner + xs = [ixmin, ixmax, ixmax, ixmin] + ys = [iymin, iymin, iymax, iymax] + tx, ty = transformer.transform(xs, ys) + inner = ( + float(min(tx)), + float(max(tx)), + float(min(ty)), + float(max(ty)), + ) + + if is_global_bbox(outer, outer_crs): return True - if is_global_bbox(inner) and not is_global_bbox(outer): + if is_global_bbox(inner, inner_crs) and not is_global_bbox(outer, outer_crs): return False oxmin, oxmax, oymin, oymax = outer @@ -80,29 +184,93 @@ def validate_crs_compatible(global_crs, regional_crs): try: g = CRS.from_user_input(global_crs) - r = CRS.from_user_input(regional_crs) + CRS.from_user_input(regional_crs) except Exception as e: return ( False, f"Failed to parse CRS. global={get_crs_string(global_crs)}, regional={get_crs_string(regional_crs)}; error={e}", ) - # Global must be EPSG:4326 (WGS84 geographic) + # Global should be geographic (preferably EPSG:4326) or a + # global-suitable projected CRS. try: g_auth = g.to_epsg() except Exception: g_auth = None - if not g.is_geographic or g_auth != 4326: - return False, f"Global domain must be EPSG:4326, found {get_crs_string(g)}" - # Regional can be geographic or projected; both acceptable - if r.is_geographic: - return True, "Compatible CRS: global EPSG:4326 with regional geographic CRS" - if r.is_projected: - return True, "Compatible CRS: global EPSG:4326 with regional projected CRS" + if g.is_geographic: + if g_auth != 4326: + return ( + True, + f"Global domain uses geographic CRS {get_crs_string(g)} (not EPSG:4326). " + "This may work but EPSG:4326 is recommended for global meshing.", + ) + return True, "Compatible CRS: global EPSG:4326 with regional CRS" + + if g.is_projected: + proj_desc = g.to_proj4() + if isinstance(proj_desc, str) and ("stere" in proj_desc.lower()): + return ( + True, + f"Global domain uses stereographic projection {get_crs_string(g)}, which is suitable for global meshing.", + ) + return ( + True, + f"Global domain uses projected CRS {get_crs_string(g)}. Ensure this projection is suitable for global-scale meshing.", + ) + + return ( + False, + f"Global domain CRS {get_crs_string(g)} is neither geographic nor projected; cannot validate compatibility.", + ) + + +def is_crs_suitable_for_global(crs): + """Return (is_suitable, message) for a CRS used as global. + + This helper prefers EPSG:4326, accepts other geographic CRS with a + warning, and treats stereographic projections as suitable for + global stereographic meshing. + """ + + if crs is None: + return False, "No CRS supplied for global domain." + + try: + c = CRS.from_user_input(crs) + except Exception as e: + return False, f"Failed to parse CRS {crs!r}: {e}" - # Fallback neutral - return True, "CRS types appear compatible" + try: + auth = c.to_epsg() + except Exception: + auth = None + + if c.is_geographic and auth == 4326: + return True, "Ideal global CRS: EPSG:4326 (WGS84 geographic)." + + if c.is_geographic: + return ( + True, + f"Global domain uses geographic CRS {get_crs_string(c)} (not EPSG:4326); acceptable but EPSG:4326 is recommended.", + ) + + if c.is_projected: + proj_desc = c.to_proj4() + if isinstance(proj_desc, str) and ("stere" in proj_desc.lower()): + return ( + True, + f"Global domain uses stereographic projection {get_crs_string(c)}, suitable for global stereographic meshing.", + ) + return ( + True, + f"Global domain uses projected CRS {get_crs_string(c)}; ensure it is appropriate for global-scale meshing.", + ) + + return ( + False, + f"CRS {get_crs_string(c)} is neither geographic nor projected; not suitable for global meshing.", + ) def warp_coordinates(points, src_crs, dst_crs): @@ -121,10 +289,18 @@ def _get_stereo_projection(scale_factor=1.0): Parameters ---------- - scale_factor: + scale_factor : float, optional Reference stereographic scale factor :math:`k_0` used in the - distortion formula. Projections are cached per distinct k0 - value when cartopy/pyproj are available. + distortion formula ``k(phi) = 2*k0 / (1 + sin(phi))``. + Projections are cached per distinct :math:`k_0` value when + cartopy/pyproj are available. + + Notes + ----- + Caching avoids repeatedly constructing cartopy CRS/Transformer + objects during mesh generation and gradient enforcement. Different + :math:`k_0` values lead to different local scale factors near the + poles and thus different effective mesh densities. """ k0_key = float(scale_factor) diff --git a/tests/test_global_regional_multiscale.py b/tests/test_global_regional_multiscale.py index c998a68..01717b5 100644 --- a/tests/test_global_regional_multiscale.py +++ b/tests/test_global_regional_multiscale.py @@ -30,6 +30,9 @@ from oceanmesh.region import to_lat_lon +GLOBAL_DATA_DIR = os.path.join(os.path.dirname(__file__), "global") + + # Force acceleration env variables for point-in-polygon os.environ["OCEANMESH_INPOLY_ACCEL"] = "1" os.environ["OCEANMESH_INPOLY_ACCEL_DEBUG"] = "1" @@ -357,6 +360,77 @@ def test_global_regional_multiscale_australia(): ), "Regional domain must not have stereo flag set" +def test_global_regional_multiscale_mixed_crs(): + """Global EPSG:4326 with regional projected CRS for multiscale meshing.""" + + # Reuse the Australia configuration but change the regional CRS to + # a projected UTM zone to exercise mixed CRS validation. + from pyproj import CRS + + # Global stereographic as in the main test + global_bbox = (-180.0, 180.0, -89.0, 90.0) + region_global = om.Region(extent=global_bbox, crs=4326) + shoreline_global = om.Shoreline( + os.path.join(GLOBAL_DATA_DIR, "global_latlon.shp"), + region_global.bbox, + 0.5, + ) + sdf_global = om.signed_distance_function(shoreline_global) + edge_global = om.distance_sizing_function(shoreline_global, rate=0.11) + edge_global = om.enforce_mesh_gradation(edge_global, gradation=0.15, stereo=True) + + # Regional Australia with projected CRS (UTM zone 55S) + regional_bbox = (110.0, 160.0, -50.0, 0.0) + utm55s = CRS.from_epsg(32755).to_epsg() + region_regional = om.Region(extent=regional_bbox, crs=4326) + shoreline_regional = om.Shoreline( + os.path.join(GLOBAL_DATA_DIR, "australia.shp"), + region_regional.bbox, + 0.25, + crs=utm55s, + ) + sdf_regional = om.signed_distance_function(shoreline_regional) + edge_regional = om.distance_sizing_function(shoreline_regional, rate=0.11) + + points, cells = om.generate_multiscale_mesh( + [sdf_global, sdf_regional], [edge_global, edge_regional], max_iter=20 + ) + assert points.shape[0] > 0 + assert cells.shape[0] > 0 + + +def test_multiscale_crs_validation_warnings(): + """Exercise CRS validation pathways for several CRS combinations.""" + + # Build simple dummy domains via Shoreline to obtain Domain + # objects with different CRS combinations without focusing on + # geometry details. + bbox = (-180.0, 180.0, -89.0, 90.0) + region_global = om.Region(extent=bbox, crs=4326) + shoreline_global = om.Shoreline( + os.path.join(GLOBAL_DATA_DIR, "global_latlon.shp"), + region_global.bbox, + 0.5, + ) + sdf_global = om.signed_distance_function(shoreline_global) + edge_global = om.distance_sizing_function(shoreline_global, rate=0.11) + + region_regional = om.Region(extent=(110.0, 160.0, -50.0, 0.0), crs=4326) + shoreline_regional = om.Shoreline( + os.path.join(GLOBAL_DATA_DIR, "australia.shp"), + region_regional.bbox, + 0.25, + ) + sdf_regional = om.signed_distance_function(shoreline_regional) + edge_regional = om.distance_sizing_function(shoreline_regional, rate=0.11) + + ok, errors = om.mesh_generator._validate_multiscale_domains( # type: ignore[attr-defined] + [sdf_global, sdf_regional], [edge_global, edge_regional] + ) + assert ok + assert errors == [] + + if __name__ == "__main__": # Allow running this test file directly as a script without pytest print("[INFO] Running global+regional multiscale test as a standalone script...") diff --git a/tests/test_global_stereo.py b/tests/test_global_stereo.py index 767a4b3..d363ad7 100644 --- a/tests/test_global_stereo.py +++ b/tests/test_global_stereo.py @@ -1,5 +1,6 @@ import os +import numpy as np import oceanmesh as om from oceanmesh.projections import CARTOPY_AVAILABLE @@ -126,3 +127,131 @@ def test_global_stereo_custom_k0(): points, cells = om.generate_mesh(domain, edge_length, stereo=True, max_iter=50) assert points.shape[0] > 0 assert cells.shape[0] > 0 + + +def test_global_stereo_alternative_crs(): + """Global stereographic meshing with a non-EPSG:4326 geographic CRS. + + Uses EPSG:4269 (NAD83) to ensure that Shoreline / Domain accept + alternative geographic CRS for global coastlines while still + producing a valid mesh. + """ + + if not CARTOPY_AVAILABLE: + # The alternative CRS test relies on cartopy-backed + # stereographic handling; skip if unavailable. + return + + bbox = (-180.00, 180.00, -89.00, 90.00) + extent = om.Region(extent=bbox, crs=4326) + + min_edge_length = 0.5 + max_edge_length = 2 + alt_crs = 4269 # NAD83 geographic + + shoreline = om.Shoreline(fname, extent.bbox, min_edge_length, crs=alt_crs) + sdf = om.signed_distance_function(shoreline) + + edge_length0 = om.distance_sizing_function(shoreline, rate=0.11) + edge_length1 = om.feature_sizing_function( + shoreline, + sdf, + min_edge_length=min_edge_length, + max_edge_length=max_edge_length, + crs=alt_crs, + ) + + edge_length = om.compute_minimum([edge_length0, edge_length1]) + edge_length = om.enforce_mesh_gradation(edge_length, gradation=0.09, stereo=True) + + shoreline_stereo = om.Shoreline(fname2, 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=20) + assert points.shape[0] > 0 + assert cells.shape[0] > 0 + + +def test_global_stereo_projected_crs(): + """Global stereographic meshing with a projected stereographic CRS.""" + + if not CARTOPY_AVAILABLE: + return + + bbox = (-180.00, 180.00, -89.00, 90.00) + extent = om.Region(extent=bbox, crs=4326) + + min_edge_length = 0.5 + max_edge_length = 2 + + # A north-polar stereographic CRS; the exact EPSG is less + # important than exercising projected global CRS handling. + stereo_epsg = 3995 + + shoreline = om.Shoreline(fname, extent.bbox, min_edge_length, crs=stereo_epsg) + sdf = om.signed_distance_function(shoreline) + + edge_length0 = om.distance_sizing_function(shoreline, rate=0.11) + edge_length1 = om.feature_sizing_function( + shoreline, + sdf, + min_edge_length=min_edge_length, + max_edge_length=max_edge_length, + crs=stereo_epsg, + ) + + edge_length = om.compute_minimum([edge_length0, edge_length1]) + edge_length = om.enforce_mesh_gradation(edge_length, gradation=0.09, stereo=True) + + shoreline_stereo = om.Shoreline(fname2, 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=20) + assert points.shape[0] > 0 + assert cells.shape[0] > 0 + + +def test_polar_gradient_distortion_correction(): + """Ensure polar gradient enforcement accounts for stereographic distortion. + + This is a smoke test that exercises enforce_mesh_gradation with + stereo=True and a custom k0. It does not assert exact numerical + values but checks that resulting sizes near the pole remain + finite and within a reasonable range. + """ + + k0 = 0.994 + bbox = (-180.00, 180.00, -89.00, 90.00) + extent = om.Region(extent=bbox, crs=4326) + + min_edge_length = 0.5 + max_edge_length = 2 + + shoreline = om.Shoreline( + fname, + extent.bbox, + min_edge_length, + scale_factor=k0, + ) + sdf = om.signed_distance_function(shoreline) + + edge_length0 = om.distance_sizing_function(shoreline, rate=0.11) + edge_length1 = om.feature_sizing_function( + shoreline, + sdf, + min_edge_length=min_edge_length, + max_edge_length=max_edge_length, + crs=4326, + ) + + edge_length = om.compute_minimum([edge_length0, edge_length1]) + limited = om.enforce_mesh_gradation( + edge_length, + gradation=0.09, + stereo=True, + scale_factor=k0, + ) + + lon, lat = limited.create_grid() + polar_mask = lat > 80.0 + polar_vals = limited.values[polar_mask] + assert np.all(np.isfinite(polar_vals)) + assert polar_vals.size > 0 From 1410c70014834aeea932775a92a69c30ab7775e1 Mon Sep 17 00:00:00 2001 From: Keith Roberts Date: Sun, 7 Dec 2025 22:22:47 -0500 Subject: [PATCH 4/4] fix syntax error --- README.md | 24 +++++++++--------------- 1 file changed, 9 insertions(+), 15 deletions(-) diff --git a/README.md b/README.md index 6b8fc27..b96a048 100644 --- a/README.md +++ b/README.md @@ -461,26 +461,23 @@ formulas (less accurate but functional). **Scale Factor (k0) Configuration** -The stereographic projection uses a **scale factor** :math:`k_0` that +The stereographic projection uses a **scale factor** `k0` that controls distortion at different latitudes. The local scale factor at -latitude :math:`\phi` is +latitude `φ` is -.. math:: +`k(φ) = 2 * k0 / (1 + sin φ)`. - k(\phi) = \frac{2 k_0}{1 + \sin \phi}. +Common `k0` values: -Common :math:`k_0` values: - -- :math:`k_0 = 1.0` (default): standard stereographic with true scale +- `k0 = 1.0` (default): standard stereographic with true scale at the pole. -- :math:`k_0 = 0.994`: used in EPSG:3413 (Arctic) and EPSG:3031 +- `k0 = 0.994`: used in EPSG:3413 (Arctic) and EPSG:3031 (Antarctic) to minimise distortion at standard parallels (~70°N/S). -The scale factor affects mesh sizing: smaller :math:`k_0` reduces -element sizes near the poles, larger :math:`k_0` increases them. +The scale factor affects mesh sizing: smaller `k0` reduces +element sizes near the poles, larger `k0` increases them. -References -~~~~~~~~~~ +#### References - [GitHub PR #87](https://github.com/CHLNDDEV/oceanmesh/pull/87): discussion of scale-factor improvements. @@ -599,9 +596,6 @@ performance-critical operations: 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