diff --git a/oceanmesh/geodata.py b/oceanmesh/geodata.py
index d7268d6..9f0a514 100644
--- a/oceanmesh/geodata.py
+++ b/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.")
 
diff --git a/oceanmesh/idw.py b/oceanmesh/idw.py
index f71b8fa..a165118 100644
--- a/oceanmesh/idw.py
+++ b/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
diff --git a/oceanmesh/mesh_generator.py b/oceanmesh/mesh_generator.py
index 06ae670..9ea56ef 100644
--- a/oceanmesh/mesh_generator.py
+++ b/oceanmesh/mesh_generator.py
@@ -678,17 +678,6 @@ def _deps_vec(i):
     return p
 
 
-def _stereo_distortion(lat):
-    # we use here Stereographic projection of the sphere
-    # from the north pole onto the plane
-    # https://en.wikipedia.org/wiki/Stereographic_projection
-    lat0 = 90
-    ll = lat + lat0
-    lrad = ll / 180 * np.pi
-    res = 2 / (1 + np.sin(lrad))
-    return res
-
-
 def _stereo_distortion_dist(lat):
     lrad = np.radians(lat)
     # Calculate the scale factor for the stereographic projection
@@ -696,21 +685,29 @@ def _stereo_distortion_dist(lat):
     return res
 
 
+def _edge_length_wgs84(lat):
+    lrad = np.radians(lat)
+    return 1 / np.cos(lrad) ** (1 / 2)
+
+
 def _generate_initial_points(min_edge_length, geps, bbox, fh, fd, pfix, stereo=False):
     """Create initial distribution in bounding box (equilateral triangles)"""
     if stereo:
-        bbox = np.array([[-180, 180], [-89, 89]])
-    p = np.mgrid[
-        tuple(slice(min, max + min_edge_length, min_edge_length) for min, max in bbox)
-    ].astype(float)
+        bbox = np.array([[-180, 180], [-90, 90]])
+    p = np.mgrid[tuple(slice(min, max, min_edge_length) for min, max in bbox)].astype(
+        float
+    )
     if stereo:
         # for global meshes in stereographic projections,
         # we need to reproject the points from lon/lat to stereo projection
         # then, we need to rectify their coordinates to lat/lon for the sizing function
+        p += (
+            np.random.rand(*p.shape) * min_edge_length / 2
+        )  # randomise the distribution
         p0 = p.reshape(2, -1).T
         x, y = to_stereo(p0[:, 0], p0[:, 1])
         p = np.asarray([x, y]).T
-        r0 = fh(to_lat_lon(p[:, 0], p[:, 1])) * _stereo_distortion(p0[:, 1])
+        r0 = fh(to_lat_lon(x, y)) * _edge_length_wgs84(p0[:, 1])
     else:
         p = p.reshape(2, -1).T
         r0 = fh(p)