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)
 ```
@@ -591,6 +593,69 @@ plt.show()
 ```
 ![Multiscale](https://user-images.githubusercontent.com/18619644/136119785-8746552d-4ff6-44c3-9aa1-3e4981ba3518.png)
 
+Global mesh generation
+----------------------
+Using oceanmesh is now possible for global meshes.
+The process is done in two steps:
+ * first the definition of the sizing functions in WGS84 coordinates,
+ * then the mesh generation is done in the stereographic projection
+
+```python
+import os
+import numpy as np
+import oceanmesh as om
+from oceanmesh.region import to_lat_lon
+import matplotlib.pyplot as plt
+# utilities functions for plotting
+def crosses_dateline(lon1, lon2):
+    return abs(lon1 - lon2) > 180
+
+def filter_triangles(points, cells):
+    filtered_cells = []
+    for cell in cells:
+        p1, p2, p3 = points[cell[0]], points[cell[1]], points[cell[2]]
+        if not (crosses_dateline(p1[0], p2[0]) or crosses_dateline(p2[0], p3[0]) or crosses_dateline(p3[0], p1[0])):
+            filtered_cells.append(cell)
+    return filtered_cells
+
+# Note: global_stereo.shp has been generated using global_tag() function in pyposeidon
+# https://github.com/ec-jrc/pyPoseidon/blob/9cfd3bbf5598c810004def83b1f43dc5149addd0/pyposeidon/boundary.py#L452
+fname = "tests/global/global_latlon.shp"
+fname2 = "tests/global/global_stereo.shp"
+
+EPSG = 4326  # EPSG:4326 or WGS84
+bbox = (-180.00, 180.00, -89.00, 90.00)
+extent = om.Region(extent=bbox, crs=4326)
+
+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)
+sdf = om.signed_distance_function(shoreline)
+edge_length = om.distance_sizing_function(shoreline, rate=0.11)
+
+# once the size functions have been defined, wed need to mesh inside domain in
+# stereographic projections. This is way we use another coastline which is
+# already translated in a sterographic projection
+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=100)
+
+# remove degenerate mesh faces and other common problems in the mesh
+points, cells = om.make_mesh_boundaries_traversable(points, cells)
+points, cells = om.delete_faces_connected_to_one_face(points, cells)
+
+# apply a Laplacian smoother
+points, cells = om.laplacian2(points, cells, max_iter=100)
+lon, lat = to_lat_lon(points[:, 0], points[:, 1])
+trin = filter_triangles(np.array([lon,lat]).T, cells)
+
+fig, ax, pc = edge_length.plot(holding = True, plot_colorbar=True, cbarlabel = "Resolution in °", cmap='magma')
+ax.triplot(lon, lat, trin, color="w", linewidth=0.25)
+plt.tight_layout()
+plt.show()
+```
+![Global](https://github.com/tomsail/oceanmesh/assets/18373442/a9c45416-78d5-4f3b-a0a3-1afd061d8dbd)
 
 See the tests inside the `testing/` folder for more inspiration. Work is ongoing on this package.
 

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