Added support for global meshes (stereographic projection)

README.md

@@ -275,12 +275,12 @@ extent = om.Region(extent=(-75.00, -70.001, 40.0001, 41.9000), crs=EPSG)
 min_edge_length = 0.01  # minimum mesh size in domain in projection
 shoreline = om.Shoreline(fname, extent.bbox, min_edge_length)
 edge_length = om.distance_sizing_function(shoreline, rate=0.15)
-ax = edge_length.plot(
+fig, ax, pc = edge_length.plot(
     xlabel="longitude (WGS84 degrees)",
     ylabel="latitude (WGS84 degrees)",
     title="Distance sizing function",
     cbarlabel="mesh size (degrees)",
-    hold=True,
+    holding=True,
 )
 shoreline.plot(ax=ax)
 ```
@@ -303,12 +303,12 @@ sdf = om.signed_distance_function(shoreline)
 edge_length = om.feature_sizing_function(
     shoreline, sdf, max_edge_length=0.05, plot=True
 )
-ax = edge_length.plot(
+fig, ax, pc = edge_length.plot(
     xlabel="longitude (WGS84 degrees)",
     ylabel="latitude (WGS84 degrees)",
     title="Feature sizing function",
     cbarlabel="mesh size (degrees)",
-    hold=True,
+    holding=True,
     xlim=[-74.3, -73.8],
     ylim=[40.3, 40.8],
 )
@@ -332,12 +332,12 @@ shoreline = om.Shoreline(fname, extent.bbox, min_edge_length)
 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)
-ax = edge_length.plot(
+fig, ax, pc = edge_length.plot(
     xlabel="longitude (WGS84 degrees)",
     ylabel="latitude (WGS84 degrees)",
     title="Feature sizing function with gradation bound",
     cbarlabel="mesh size (degrees)",
-    hold=True,
+    holding=True,
     xlim=[-74.3, -73.8],
     ylim=[40.3, 40.8],
 )
@@ -358,7 +358,8 @@ fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
 
 min_edge_length = 0.01
 
-dem = om.DEM(fdem, crs=4326)
+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, min_edge_length)
 sdf = om.signed_distance_function(shoreline)
 edge_length1 = om.feature_sizing_function(shoreline, sdf, max_edge_length=0.05)
@@ -368,12 +369,12 @@ edge_length2 = om.wavelength_sizing_function(
 # Compute the minimum of the sizing functions
 edge_length = om.compute_minimum([edge_length1, edge_length2])
 edge_length = om.enforce_mesh_gradation(edge_length, gradation=0.15)
-ax = edge_length.plot(
+fig, ax, pc = edge_length.plot(
     xlabel="longitude (WGS84 degrees)",
     ylabel="latitude (WGS84 degrees)",
     title="Feature sizing function + wavelength + gradation bound",
     cbarlabel="mesh size (degrees)",
-    hold=True,
+    holding=True,
     xlim=[-74.3, -73.8],
     ylim=[40.3, 40.8],
 )
@@ -385,6 +386,7 @@ shoreline.plot(ax=ax)
 
 The distance, feature size, and/or wavelength mesh size functions can lead to coarse mesh resolution in deeper waters that under-resolve and smooth over the sharp topographic gradients that characterize the continental shelf break. These slope features can be important for coastal ocean models in order to capture dissipative effects driven by the internal tides, transmissional reflection at the shelf break that controls the astronomical tides, and trapped shelf waves. The bathymetry field contains often excessive details that are not relevant for most flows, thus the bathymetry can be smoothed by a variety of filters (e.g., lowpass, bandpass, and highpass filters) before calculating the mesh sizes.
 
+<!--pytest-codeblocks:skip-->
 ```python
 import oceanmesh as om
 
@@ -394,7 +396,7 @@ fname = "gshhg-shp-2.3.7/GSHHS_shp/f/GSHHS_f_L1.shp"
 EPSG = 4326  # EPSG:4326 or WGS84
 bbox = (-74.4, -73.4, 40.2, 41.2)
 extent = om.Region(extent=bbox, crs=EPSG)
-dem = om.DEM(fdem, crs=4326)
+dem = om.DEM(fdem, crs=EPSG)
 
 min_edge_length = 0.0025  # minimum mesh size in domain in projection
 max_edge_length = 0.10  # maximum mesh size in domain in projection
@@ -414,7 +416,7 @@ edge_length2 = om.bathymetric_gradient_sizing_function(
     min_edge_length=min_edge_length,
     max_edge_length=max_edge_length,
     crs=EPSG,
-)
+)  # will be reactivated
 edge_length3 = om.compute_minimum([edge_length1, edge_length2])
 edge_length3 = om.enforce_mesh_gradation(edge_length3, gradation=0.15)
 ```

oceanmesh/__init__.py

@@ -38,7 +38,14 @@
     get_polygon_coordinates,
 )
 from oceanmesh.grid import Grid, compute_minimum
-from oceanmesh.region import Region, warp_coordinates
+from oceanmesh.region import (
+    Region,
+    stereo_to_3d,
+    to_3d,
+    to_lat_lon,
+    to_stereo,
+    warp_coordinates,
+)
 from oceanmesh.signed_distance_function import (
     Difference,
     Domain,
@@ -63,6 +70,10 @@
 __all__ = [
     "create_bbox",
     "Region",
+    "stereo_to_3d",
+    "to_lat_lon",
+    "to_3d",
+    "to_stereo",
     "compute_minimum",
     "create_circle_coords",
     "bathymetric_gradient_sizing_function",

oceanmesh/boundary.py

@@ -1,7 +1,7 @@
+import matplotlib.pyplot as plt
 import numpy as np
 
 from .edges import get_winded_boundary_edges
-import matplotlib.pyplot as plt
 
 __all__ = ["identify_ocean_boundary_sections"]
 

oceanmesh/clean.py

@@ -256,7 +256,7 @@ def delete_exterior_faces(vertices, faces, min_disconnected_area):
         # Delete where nflag == 1 from tmp t1 mesh
         t1 = np.delete(t1, nflag == 1, axis=0)
         logger.info(
-            f"ACCEPTED: Deleting {int(np.sum(nflag==0))} faces outside the main mesh"
+            f"ACCEPTED: Deleting {int(np.sum(nflag == 0))} faces outside the main mesh"
         )
 
         # Calculate the remaining area

oceanmesh/edgefx.py

@@ -7,15 +7,16 @@
 import numpy as np
 import scipy.spatial
 import skfmm
-from _HamiltonJacobi import gradient_limit
 from inpoly import inpoly2
 from shapely.geometry import LineString
 from skimage.morphology import medial_axis
 
+from _HamiltonJacobi import gradient_limit
 from oceanmesh.filterfx import filt2
 
 from . import edges
 from .grid import Grid
+from .region import to_lat_lon, to_stereo
 
 logger = logging.getLogger(__name__)
 
@@ -87,7 +88,7 @@ def enforce_mesh_size_bounds_elevation(grid, dem, bounds):
     return grid
 
 
-def enforce_mesh_gradation(grid, gradation=0.15, crs="EPSG:4326"):
+def enforce_mesh_gradation(grid, gradation=0.15, crs="EPSG:4326", stereo=False):
     """Enforce a mesh size gradation bound `gradation` on a :class:`grid`
 
     Parameters
@@ -122,6 +123,46 @@ def enforce_mesh_gradation(grid, gradation=0.15, crs="EPSG:4326"):
     cell_size = cell_size.flatten("F")
     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
+
+        us, vs = np.meshgrid(
+            np.arange(-1, 1, dx_stereo), np.arange(-1, 1, dx_stereo), indexing="ij"
+        )
+        ulon, vlat = to_lat_lon(us.ravel(), vs.ravel())
+        utmp = grid.eval((ulon, vlat))
+        utmp = np.reshape(utmp, us.shape)
+        szs = utmp.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 = np.reshape(vtmp, (szs[0], szs[1]), "F")
+        # construct stereo interpolating function
+        grid_stereo = Grid(
+            bbox=(-1, 1, -1, 1),
+            dx=dx_stereo,
+            values=vtmp,
+            hmin=grid.hmin,
+            extrapolate=grid.extrapolate,
+            crs=crs,
+        )
+        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]))
+        logger.info(
+            "Global mesh: reinject back stereographic gradient and recomputing gradient..."
+        )
+        cell_size = tmp.flatten("F")
+        tmp = gradient_limit([*sz], elen, gradation, 10000, cell_size)
+        tmp = np.reshape(tmp, (sz[0], sz[1]), "F")
+
     grid_limited = Grid(
         bbox=grid.bbox,
         dx=grid.dx,
@@ -453,7 +494,7 @@ def bathymetric_gradient_sizing_function(
     nx, ny = dem.nx, dem.ny
     coords = (xg, yg)
 
-    if crs == "EPSG:4326":
+    if crs == "EPSG:4326" or crs == 4326:
         mean_latitude = np.mean(dem.bbox[2:])
         meters_per_degree = (
             111132.92
@@ -463,7 +504,6 @@ def bathymetric_gradient_sizing_function(
         )
         dy *= meters_per_degree
         dx *= meters_per_degree
-
     grid_details = (nx, ny, dx, dy)
 
     if type_of_filter == "barotropic" and filter_quotient > 0:
@@ -500,7 +540,7 @@ def bathymetric_gradient_sizing_function(
     grid.values = (2 * np.pi / slope_parameter) * np.abs(dp) / (bs + eps)
 
     # Convert back to degrees from meters (if geographic)
-    if crs == "EPSG:4326":
+    if crs == "EPSG:4326" or crs == 4326:
         grid.values /= meters_per_degree
 
     if max_edge_length is not None:
@@ -818,23 +858,26 @@ def wavelength_sizing_function(
     lon, lat = dem.create_grid()
     tmpz = dem.eval((lon, lat))
 
-    if crs == "EPSG:4326":
+    dx, dy = dem.dx, dem.dy  # for gradient function
+
+    if crs == "EPSG:4326" or crs == 4326:
         mean_latitude = np.mean(dem.bbox[2:])
         meters_per_degree = (
             111132.92
             - 559.82 * np.cos(2 * mean_latitude)
             + 1.175 * np.cos(4 * mean_latitude)
             - 0.0023 * np.cos(6 * mean_latitude)
         )
-
+        dy *= meters_per_degree
+        dx *= meters_per_degree
     grid = Grid(
         bbox=dem.bbox, dx=dem.dx, dy=dem.dy, extrapolate=True, values=0.0, crs=crs
     )
     tmpz[np.abs(tmpz) < 1] = 1
     grid.values = period * np.sqrt(gravity * np.abs(tmpz)) / wl
 
     # Convert back to degrees from meters (if geographic)
-    if crs == "EPSG:4326":
+    if crs == "EPSG:4326" or crs == 4326:
         grid.values /= meters_per_degree
 
     if min_edgelength is None:
@@ -892,7 +935,7 @@ def multiscale_sizing_function(
         # interpolate all finer nests onto coarse func and enforce gradation rate
         for k, finer in enumerate(list_of_grids[idx1 + 1 :]):
             logger.info(
-                f"  Interpolating sizing function #{idx1+1 + k} onto sizing function #{idx1}"
+                f"  Interpolating sizing function #{idx1 + 1 + k} onto sizing function #{idx1}"
             )
             _wkt = finer.crs.to_dict()
             if "units" in _wkt:

oceanmesh/filterfx.py

@@ -26,7 +26,7 @@ def filt2(Z, res, wl, filtertype, truncate=2.6):
 highpass, bandpass or bandstop"
         )
 
-    if hasattr(wl, "__len__") and type(wl) != np.ndarray:
+    if hasattr(wl, "__len__") and ~isinstance(type(wl), np.ndarray):
         wl = np.array(wl)
 
     if np.any(wl <= 2 * res):
@@ -37,12 +37,12 @@ def filt2(Z, res, wl, filtertype, truncate=2.6):
 
     if filtertype in ["bandpass", "bandstop"]:
         if hasattr(wl, "__len__"):
-            if len(wl) != 2 or type(wl) == str:
+            if len(wl) != 2 or isinstance(wl, str):
                 raise TypeError(
                     "Wavelength lambda must be a two-element array for a bandpass filter."
                 )
 
-            if type(wl) != np.array:
+            if ~isinstance(wl, np.ndarray):
                 wl = np.array(list(wl))
 
         else:

oceanmesh/fix_mesh.py

@@ -1,6 +1,7 @@
-import numpy as np
 import warnings
 
+import numpy as np
+
 
 def simp_qual(p, t):
     """Simplex quality radius-to-edge ratio

oceanmesh/geodata.py

@@ -18,7 +18,7 @@
 from rasterio.windows import from_bounds
 
 from .grid import Grid
-from .region import Region
+from .region import Region, to_lat_lon, to_stereo
 
 nan = np.nan
 fiona_version = fiona.__version__
@@ -174,7 +174,7 @@ def _poly_length(coords):
     return np.sum(np.sqrt(np.sum(np.diff(c, axis=0) ** 2, axis=1)))
 
 
-def _classify_shoreline(bbox, boubox, polys, h0, minimum_area_mult):
+def _classify_shoreline(bbox, boubox, polys, h0, minimum_area_mult, stereo=False):
     """Classify segments in numpy.array `polys` as either `inner` or `mainland`.
     (1) The `mainland` category contains segments that are not totally enclosed inside the `bbox`.
     (2) The `inner` (i.e., islands) category contains segments totally enclosed inside the `bbox`.
@@ -210,13 +210,21 @@ def _classify_shoreline(bbox, boubox, polys, h0, minimum_area_mult):
     for poly in polyL:
         pSGP = shapely.geometry.Polygon(poly[:-2, :])
         if bSGP.contains(pSGP):
-            if pSGP.area >= _AREAMIN:
+            if stereo:
+                # convert back to Lat/Lon coordinates for the area testing
+                area = _poly_area(*to_lat_lon(*np.asarray(pSGP.exterior.xy)))
+            else:
+                area = pSGP.area
+            if area >= _AREAMIN:
                 inner = np.append(inner, poly, axis=0)
         elif pSGP.overlaps(bSGP):
-            # Append polygon segment to mainland
-            mainland = np.vstack((mainland, poly))
-            # Clip polygon segment from boubox and regenerate path
-            bSGP = bSGP.difference(pSGP)
+            if stereo:
+                bSGP = pSGP
+            else:
+                bSGP = bSGP.difference(pSGP)
+                # Append polygon segment to mainland
+                mainland = np.vstack((mainland, poly))
+                # Clip polygon segment from boubox and regenerate path
 
     out = np.empty(shape=(0, 2))
 
@@ -515,6 +523,7 @@ def __init__(
         refinements=1,
         minimum_area_mult=4.0,
         smooth_shoreline=True,
+        stereo=False,
     ):
         if isinstance(shp, str):
             shp = Path(shp)
@@ -545,6 +554,20 @@ def __init__(
 
         polys = self._read()
 
+        if stereo:
+            u, v = to_stereo(polys[:, 0], polys[:, 1])
+            polys = np.array([u, v]).T
+            polys = polys[
+                ~(np.isinf(u) | np.isinf(v))
+            ]  # remove eventual inf values (south pole)
+            self.bbox = (
+                np.nanmin(polys[:, 0] * 0.99),
+                np.nanmax(polys[:, 0] * 0.99),
+                np.nanmin(polys[:, 1] * 0.99),
+                np.nanmax(polys[:, 1] * 0.99),
+            )  # so that bbox overlaps with antarctica > and becomes the outer boundary
+            self.boubox = np.asarray(_create_boubox(self.bbox))
+
         if smooth_shoreline:
             polys = _smooth_shoreline(polys, self.refinements)
 
@@ -553,7 +576,7 @@ def __init__(
         polys = _clip_polys(polys, self.bbox)
 
         self.inner, self.mainland, self.boubox = _classify_shoreline(
-            self.bbox, self.boubox, polys, self.h0 / 2, self.minimum_area_mult
+            self.bbox, self.boubox, polys, self.h0 / 2, self.minimum_area_mult, stereo
         )
 
     @property