Improve global stereographic meshing with configurable k0 and CRS support

README.md

@@ -428,6 +428,64 @@ 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** `k0` that
+controls distortion at different latitudes. The local scale factor at
+latitude `φ` is
+
+`k(φ) = 2 * k0 / (1 + sin φ)`.
+
+Common `k0` values:
+
+- `k0 = 1.0` (default): standard stereographic with true scale
+  at the pole.
+- `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 `k0` reduces
+element sizes near the poles, larger `k0` 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.
+
 <!--pytest-codeblocks:skip-->
 
 ```python
@@ -442,7 +500,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 +519,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 +533,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 +560,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)
@@ -508,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

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

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,28 +126,78 @@ 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
+
+        backend = "cartopy" if CARTOPY_AVAILABLE else "hardcoded"
+        logger.info(
+            "Global mesh: fixing gradient on the north pole (backend=%s, k0=%s, distortion_correction=%s)...",
+            backend,
+            scale_factor,
+            "enabled" if CARTOPY_AVAILABLE else "disabled",
+        )
+
+        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())
+            # 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. 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")
-        # construct stereo interpolating function
+
+        # 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=(-1, 1, -1, 1),
+            bbox=(-r_max, r_max, -r_max, r_max),
             dx=dx_stereo,
             values=vtmp,
             hmin=grid.hmin,
@@ -155,7 +207,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..."
         )
@@ -824,7 +880,19 @@ def feature_sizing_function(
     # find location of points on grid
     indices = grid_calc.find_indices(points, lon, lat)
     phi2[indices] = -1.0
-    dis = np.abs(skfmm.distance(phi2, [grid_calc.dx, grid_calc.dy]))
+    # Ensure the level-set field contains a zero contour before calling skfmm.
+    # skfmm.distance requires that the input field change sign; if it does not,
+    # it raises a ValueError indicating that there is no zero contour.
+    has_neg = np.any(phi2 < 0.0)
+    has_pos = np.any(phi2 > 0.0)
+    if not (has_neg and has_pos):
+        logger.warning(
+            "Level set field for feature sizing has no zero contour; "
+            "skipping skfmm.distance and using approximate distances based on grid spacing."
+        )
+        dis = np.ones_like(phi2) * max(grid_calc.dx, grid_calc.dy)
+    else:
+        dis = np.abs(skfmm.distance(phi2, [grid_calc.dx, grid_calc.dy]))
 
     if plot:
         fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4))