added example in README + supporting files

README.md

@@ -593,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 = os.path.join(os.path.dirname(__file__), "global", "global_latlon.shp")
+fname2 = os.path.join(os.path.dirname(__file__), "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/geodata.py

@@ -18,7 +18,7 @@
 from rasterio.windows import from_bounds
 
 from .grid import Grid
-from .region import Region, to_lat_lon, to_stereo
+from .region import Region, to_lat_lon
 
 nan = np.nan
 fiona_version = fiona.__version__
@@ -555,11 +555,6 @@ 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),

tests/global/global_latlon.cpg

@@ -0,0 +1 @@
+ISO-8859-1

tests/global/global_latlon.prj

@@ -0,0 +1 @@
+GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

tests/global/global_stereo.cpg

@@ -0,0 +1 @@
+ISO-8859-1

tests/global/global_stereo.prj

@@ -0,0 +1 @@
+GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

tests/test_global_stereo.py

@@ -2,11 +2,10 @@
 
 import oceanmesh as om
 
-# Note: to following file is partof the test data from pyposeidon Tutorial notebooks
-# https://github.com/ec-jrc/pyPoseidon/tree/master/Tutorial
-# You can get them with
-# curl -L url https://www.dropbox.com/sh/nd2b012wrpt6qnq/AAAD7aA_qXztUhlT39YK2yBua?dl=1 > data.zip
-fname = os.path.join(os.path.dirname(__file__), "data", "ocean", "ocean.shp")
+# Note: global_stereo.shp has been generated using global_tag in pyposeidon
+# https://github.com/ec-jrc/pyPoseidon/blob/9cfd3bbf5598c810004def83b1f43dc5149addd0/pyposeidon/boundary.py#L452
+fname = os.path.join(os.path.dirname(__file__), "global", "global_latlon.shp")
+fname2 = os.path.join(os.path.dirname(__file__), "global", "global_stereo.shp")
 
 
 def test_global_stereo():
@@ -34,7 +33,10 @@ def test_global_stereo():
     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(fname, extent.bbox, min_edge_length, stereo=True)
+    # 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)