added parametric law for stereo distortion (north pole)

oceanmesh/geodata.py

@@ -794,9 +794,9 @@ def __init__(self, dem, crs="EPSG:4326", bbox=None, extrapolate=False):
                     topobathy = np.transpose(topobathy, (1, 0))
             # Ensure its a floating point array
             topobathy = topobathy.astype(np.float64)
-            topobathy[
-                topobathy == nodata_value
-            ] = np.nan  # set the no-data value to nan
+            topobathy[topobathy == nodata_value] = (
+                np.nan
+            )  # set the no-data value to nan
         elif not dem.exists():
             raise FileNotFoundError(f"File {dem} could not be located.")
 

oceanmesh/idw.py

@@ -1,6 +1,7 @@
 """ invdisttree.py: inverse-distance-weighted interpolation using KDTree
     fast, solid, local
 """
+
 from __future__ import division
 
 import numpy as np

oceanmesh/mesh_generator.py

@@ -581,7 +581,7 @@ def _compute_forces(p, t, fh, min_edge_length, L0mult, opts):
         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) * _parametric(y)
     else:
         hbars = fh(p[bars].sum(1) / 2)
     L0 = hbars * L0mult * (np.nanmedian(L) / np.nanmedian(hbars))
@@ -696,6 +696,12 @@ def _stereo_distortion_dist(lat):
     return res
 
 
+def _parametric(lat):
+    ones = np.ones(lat.shape)
+    res = ((90 - lat) * 2 + 18) / 180 * np.pi
+    return np.minimum(res, ones)
+
+
 def _generate_initial_points(min_edge_length, geps, bbox, fh, fd, pfix, stereo=False):
     """Create initial distribution in bounding box (equilateral triangles)"""
     if stereo: