From cde1ce3ed44c6d5b8652f741bec484ef4f3df035 Mon Sep 17 00:00:00 2001
From: Coat <alain.coat@actimar.fr>
Date: Fri, 24 Jun 2022 16:10:17 +0200
Subject: [PATCH 01/34] remove nan nodes of poly from edges

---
 oceanmesh/edges.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/oceanmesh/edges.py b/oceanmesh/edges.py
index e6b4137..563495d 100644
--- a/oceanmesh/edges.py
+++ b/oceanmesh/edges.py
@@ -22,13 +22,13 @@ def get_poly_edges(poly):
         A 2D array of integers containing indexes into the `poly` array.
 
     """
-    ix = np.argwhere(np.isnan(poly[:, 0]))
-    ix = np.insert(ix, 0, 0)
+    ix = np.argwhere(np.isnan(poly[:, 0])).ravel()
+    ix = np.insert(ix, 0, -1)
 
     edges = []
     for s in range(len(ix) - 1):
-        col1 = np.arange(ix[s], ix[s + 1] - 1)
-        col2 = np.arange(ix[s] + 1, ix[s + 1])
+        col1 = np.arange(ix[s] + 1, ix[s + 1] - 2)
+        col2 = np.arange(ix[s] + 2, ix[s + 1] - 1)
         tmp = np.vstack((col1, col2)).T
         tmp = np.append(tmp, [[ix[s + 1] - 1, ix[s] + 1]], axis=0)
         edges.append(tmp)

From cc1e46571f2838282e90c4650b5542bc0e16d4c7 Mon Sep 17 00:00:00 2001
From: Coat <alain.coat@actimar.fr>
Date: Fri, 24 Jun 2022 17:09:29 +0200
Subject: [PATCH 02/34] remove nan from nodes for in_poly2

---
 oceanmesh/edgefx.py                   | 9 +++++----
 oceanmesh/signed_distance_function.py | 4 ++--
 2 files changed, 7 insertions(+), 6 deletions(-)

diff --git a/oceanmesh/edgefx.py b/oceanmesh/edgefx.py
index 2d40271..bfee0e7 100644
--- a/oceanmesh/edgefx.py
+++ b/oceanmesh/edgefx.py
@@ -168,8 +168,8 @@ def distance_sizing_function(
     lon, lat = grid.create_grid()
     points = np.vstack((shoreline.inner, shoreline.mainland))
     # remove shoreline components outside the shoreline.boubox
-    boubox = shoreline.boubox
-    e_box = edges.get_poly_edges(boubox)
+    boubox = np.nan_to_num(shoreline.boubox) # remove nan for inpoly2
+    e_box = edges.get_poly_edges(shoreline.boubox)
     mask = np.ones((grid.nx, grid.ny), dtype=bool)
     if len(points) > 0:
         try:
@@ -177,10 +177,11 @@ def distance_sizing_function(
             points = points[in_boubox]
 
             qpts = np.column_stack((lon.flatten(), lat.flatten()))
-            in_boubox, _ = inpoly2(qpts, shoreline.boubox, e_box)
+            in_boubox, _ = inpoly2(qpts, boubox, e_box)
             mask_indices = grid.find_indices(qpts[in_boubox, :], lon, lat)
             mask[mask_indices] = False
-        except (Exception,):
+        except Exception as e:
+            logger.error(e)
             ...
 
     # find location of points on grid
diff --git a/oceanmesh/signed_distance_function.py b/oceanmesh/signed_distance_function.py
index 74460f8..b0bcae0 100644
--- a/oceanmesh/signed_distance_function.py
+++ b/oceanmesh/signed_distance_function.py
@@ -206,9 +206,9 @@ def func(x):
         # Initialize d with some positive number larger than geps
         dist = np.zeros(len(x)) + 1.0
         # are points inside the boubox?
-        in_boubox, _ = inpoly2(x, boubox, e_box)
+        in_boubox, _ = inpoly2(x, np.nan_to_num(boubox), e_box)
         # are points inside the shoreline?
-        in_shoreline, _ = inpoly2(x, poly, e)
+        in_shoreline, _ = inpoly2(x, np.nan_to_num(poly), e)
         # compute dist to shoreline
         try:
             d, _ = tree.query(x, k=1, workers=-1)

From 6ec738b5c8626f150068a77e66722972271d9fd9 Mon Sep 17 00:00:00 2001
From: Coat <alain.coat@actimar.fr>
Date: Fri, 24 Jun 2022 18:56:58 +0200
Subject: [PATCH 03/34] fix E261

---
 oceanmesh/edgefx.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/oceanmesh/edgefx.py b/oceanmesh/edgefx.py
index bfee0e7..c7f1894 100644
--- a/oceanmesh/edgefx.py
+++ b/oceanmesh/edgefx.py
@@ -168,7 +168,7 @@ def distance_sizing_function(
     lon, lat = grid.create_grid()
     points = np.vstack((shoreline.inner, shoreline.mainland))
     # remove shoreline components outside the shoreline.boubox
-    boubox = np.nan_to_num(shoreline.boubox) # remove nan for inpoly2
+    boubox = np.nan_to_num(shoreline.boubox)  # remove nan for inpoly2
     e_box = edges.get_poly_edges(shoreline.boubox)
     mask = np.ones((grid.nx, grid.ny), dtype=bool)
     if len(points) > 0:

From f3c15e90abaeafee9426c5de0f40bf333327379f Mon Sep 17 00:00:00 2001
From: Coat <alain.coat@actimar.fr>
Date: Fri, 24 Jun 2022 23:41:08 +0200
Subject: [PATCH 04/34] reformat with black

---
 oceanmesh/__init__.py                       | 13 ++++++----
 oceanmesh/edgefx.py                         | 17 ++++---------
 oceanmesh/fix_mesh.py                       |  2 +-
 oceanmesh/geodata.py                        | 13 +++-------
 oceanmesh/idw.py                            |  2 +-
 oceanmesh/mesh_generator.py                 | 27 ++++++---------------
 oceanmesh/signed_distance_function.py       | 16 ++----------
 setup.py                                    | 18 +++++---------
 tests/test_bathymetric_gradient_function.py |  5 +---
 tests/test_circ_rect.py                     |  4 ++-
 tests/test_geodata.py                       | 16 +++---------
 tests/test_grade.py                         |  2 ++
 tests/test_irregular_domain.py              |  4 ++-
 tests/test_mesh_generator_simple.py         |  1 +
 tests/test_multiscale.py                    |  1 +
 15 files changed, 49 insertions(+), 92 deletions(-)

diff --git a/oceanmesh/__init__.py b/oceanmesh/__init__.py
index cd2827c..11667e1 100644
--- a/oceanmesh/__init__.py
+++ b/oceanmesh/__init__.py
@@ -2,9 +2,11 @@
 
 # DPZ patches for miniconda on windows using vcpkg to install cgal and dependencies
 #
-if os.name == 'nt':
-    assert os.environ.get('CGAL_BIN', None), "The environment variable CGAL_BIN must be set."
-    os.add_dll_directory(os.environ['CGAL_BIN'])
+if os.name == "nt":
+    assert os.environ.get(
+        "CGAL_BIN", None
+    ), "The environment variable CGAL_BIN must be set."
+    os.add_dll_directory(os.environ["CGAL_BIN"])
 
 from oceanmesh.clean import (
     delete_boundary_faces,
@@ -24,6 +26,7 @@
     wavelength_sizing_function,
 )
 from oceanmesh.edges import draw_edges, get_poly_edges
+from oceanmesh.filterfx import filt2
 from oceanmesh.geodata import DEM, Shoreline
 from oceanmesh.grid import Grid, compute_minimum
 from oceanmesh.region import Region, warp_coordinates
@@ -32,12 +35,12 @@
     Domain,
     Intersection,
     Union,
-    create_circle,
     create_bbox,
+    create_circle,
     multiscale_signed_distance_function,
     signed_distance_function,
 )
-from oceanmesh.filterfx import filt2
+
 from .fix_mesh import fix_mesh, simp_vol
 from .mesh_generator import generate_mesh, generate_multiscale_mesh, plot_mesh
 
diff --git a/oceanmesh/edgefx.py b/oceanmesh/edgefx.py
index c7f1894..0b4b766 100644
--- a/oceanmesh/edgefx.py
+++ b/oceanmesh/edgefx.py
@@ -3,10 +3,10 @@
 import numpy as np
 import scipy.spatial
 import skfmm
-from _HamiltonJacobi import gradient_limit
 from inpoly import inpoly2
 from skimage.morphology import medial_axis
 
+from _HamiltonJacobi import gradient_limit
 from oceanmesh.filterfx import filt2
 
 from . import edges
@@ -128,11 +128,7 @@ def enforce_mesh_gradation(grid, gradation=0.15, crs=4326):
 
 
 def distance_sizing_function(
-    shoreline,
-    rate=0.15,
-    max_edge_length=None,
-    coarsen=1,
-    crs=4326,
+    shoreline, rate=0.15, max_edge_length=None, coarsen=1, crs=4326,
 ):
     """Mesh sizes that vary linearly at `rate` from coordinates in `obj`:Shoreline
     Parameters
@@ -288,7 +284,7 @@ def bathymetric_gradient_sizing_function(
         msg = f"The type_of_filter {type_of_filter} is not known and remains off"
         logger.info(msg)
         by, bx = _earth_gradient(tmpz, dy, dx)  # get slope in x and y directions
-        bs = np.sqrt(bx**2 + by**2)  # get overall slope
+        bs = np.sqrt(bx ** 2 + by ** 2)  # get overall slope
 
     # Calculating the slope function
     eps = 1e-10  # small number to approximate derivative
@@ -443,7 +439,7 @@ def rossby_radius_filter(tmpz, bbox, grid_details, coords, rbfilt, barot):
             by, bx = _earth_gradient(
                 tmpz_ft, dy, dx
             )  # [n2s:n2e]) # get slope in x and y directions
-            tempbs = np.sqrt(bx**2 + by**2)  # get overall slope
+            tempbs = np.sqrt(bx ** 2 + by ** 2)  # get overall slope
 
             bst[rosb == edges[i]] = tempbs[rosb == edges[i]]
 
@@ -632,10 +628,7 @@ def wavelength_sizing_function(
 
 
 def multiscale_sizing_function(
-    list_of_grids,
-    p=3,
-    nnear=28,
-    blend_width=1000,
+    list_of_grids, p=3, nnear=28, blend_width=1000,
 ):
     """Given a list of mesh size functions in a hierarchy
     w.r.t. to minimum mesh size (largest -> smallest),
diff --git a/oceanmesh/fix_mesh.py b/oceanmesh/fix_mesh.py
index 6b3c36d..3f45044 100644
--- a/oceanmesh/fix_mesh.py
+++ b/oceanmesh/fix_mesh.py
@@ -13,7 +13,7 @@ def simp_qual(p, t):
     assert p.ndim == 2 and t.ndim == 2 and p.shape[1] + 1 == t.shape[1]
 
     def length(p1):
-        return np.sqrt((p1**2).sum(1))
+        return np.sqrt((p1 ** 2).sum(1))
 
     a = length(p[t[:, 1]] - p[t[:, 0]])
     b = length(p[t[:, 2]] - p[t[:, 0]])
diff --git a/oceanmesh/geodata.py b/oceanmesh/geodata.py
index 5d30d21..2df1723 100644
--- a/oceanmesh/geodata.py
+++ b/oceanmesh/geodata.py
@@ -96,9 +96,7 @@ def _densify(poly, maxdiff, bbox, radius=0):
         else:
             ni = int(ni)
             icoords = _create_ranges(
-                np.array([lat[i], lon[i]]),
-                np.array([lat[i + 1], lon[i + 1]]),
-                ni + 2,
+                np.array([lat[i], lon[i]]), np.array([lat[i + 1], lon[i + 1]]), ni + 2,
             )
             latout[n : n + ni + 1] = icoords[0, : ni + 1]
             lonout[n : n + ni + 1] = icoords[1, : ni + 1]
@@ -137,7 +135,7 @@ def _classify_shoreline(bbox, boubox, polys, h0, minimum_area_mult):
     """
     logger.debug("Entering:_classify_shoreline")
 
-    _AREAMIN = minimum_area_mult * h0**2
+    _AREAMIN = minimum_area_mult * h0 ** 2
 
     if len(boubox) == 0:
         boubox = _create_boubox(bbox)
@@ -606,12 +604,7 @@ def plot(
 
         xmin, xmax, ymin, ymax = self.bbox
         rect = plt.Rectangle(
-            (xmin, ymin),
-            xmax - xmin,
-            ymax - ymin,
-            fill=None,
-            hatch="////",
-            alpha=0.2,
+            (xmin, ymin), xmax - xmin, ymax - ymin, fill=None, hatch="////", alpha=0.2,
         )
 
         border = 0.10 * (xmax - xmin)
diff --git a/oceanmesh/idw.py b/oceanmesh/idw.py
index f71b8fa..9705f9a 100644
--- a/oceanmesh/idw.py
+++ b/oceanmesh/idw.py
@@ -92,7 +92,7 @@ def __call__(self, q, nnear=6, eps=0, p=1, weights=None):
             elif dist[0] < 1e-10:
                 wz = self.z[ix[0]]
             else:  # weight z s by 1/dist --
-                w = 1 / dist**p
+                w = 1 / dist ** p
                 if weights is not None:
                     w *= weights[ix]  # >= 0
                 w /= np.sum(w)
diff --git a/oceanmesh/mesh_generator.py b/oceanmesh/mesh_generator.py
index 65f3b24..da3f07c 100644
--- a/oceanmesh/mesh_generator.py
+++ b/oceanmesh/mesh_generator.py
@@ -5,6 +5,7 @@
 import matplotlib.tri as tri
 import numpy as np
 import scipy.sparse as spsparse
+
 from _delaunay_class import DelaunayTriangulation as DT
 from _fast_geometry import unique_edges
 
@@ -213,14 +214,7 @@ def generate_mesh(domain, edge_length, **kwargs):
     lock_boundary = opts["lock_boundary"]
 
     if opts["points"] is None:
-        p = _generate_initial_points(
-            min_edge_length,
-            geps,
-            bbox,
-            fh,
-            fd,
-            pfix,
-        )
+        p = _generate_initial_points(min_edge_length, geps, bbox, fh, fd, pfix,)
     else:
         p = opts["points"]
 
@@ -281,7 +275,7 @@ def generate_mesh(domain, edge_length, **kwargs):
         p = _project_points_back(p, fd, deps)
 
         # Show the user some progress so they know something is happening
-        maxdp = delta_t * np.sqrt((Ftot**2).sum(1)).max()
+        maxdp = delta_t * np.sqrt((Ftot ** 2).sum(1)).max()
 
         logger.info(
             f"Iteration #{count+1}, max movement is {maxdp}, there are {len(p)} vertices and {len(t)}"
@@ -331,10 +325,10 @@ def _compute_forces(p, t, fh, min_edge_length, L0mult):
     N = p.shape[0]
     bars = _get_bars(t)
     barvec = p[bars[:, 0]] - p[bars[:, 1]]  # List of bar vectors
-    L = np.sqrt((barvec**2).sum(1))  # L = Bar lengths
+    L = np.sqrt((barvec ** 2).sum(1))  # L = Bar lengths
     L[L == 0] = np.finfo(float).eps
     hbars = fh(p[bars].sum(1) / 2)
-    L0 = hbars * L0mult * ((L**2).sum() / (hbars**2).sum()) ** (1.0 / 2)
+    L0 = hbars * L0mult * ((L ** 2).sum() / (hbars ** 2).sum()) ** (1.0 / 2)
     F = L0 - L
     F[F < 0] = 0  # Bar forces (scalars)
     Fvec = (
@@ -422,7 +416,7 @@ def _deps_vec(i):
             # take the solely ones outside domain
             dgrads = [(fd(p + _deps_vec(i)) - d) / deps for i in range(2)]
             dgrads = list(np.array(dgrads)[:, ix])
-        dgrad2 = sum(dgrad**2 for dgrad in dgrads)
+        dgrad2 = sum(dgrad ** 2 for dgrad in dgrads)
         dgrad2 = np.where(dgrad2 < deps, deps, dgrad2)
         p[ix] -= (d[ix] * np.vstack(dgrads) / dgrad2).T  # Project
     return p
@@ -436,14 +430,9 @@ def _generate_initial_points(min_edge_length, geps, bbox, fh, fd, pfix):
     p = p.reshape(2, -1).T
     r0 = fh(p)
     r0m = np.min(r0[r0 >= min_edge_length])
-    p = p[np.random.rand(p.shape[0]) < r0m**2 / r0**2]
+    p = p[np.random.rand(p.shape[0]) < r0m ** 2 / r0 ** 2]
     p = p[fd(p) < geps]  # Keep only d<0 points
-    return np.vstack(
-        (
-            pfix,
-            p,
-        )
-    )
+    return np.vstack((pfix, p,))
 
 
 def _dist(p1, p2):
diff --git a/oceanmesh/signed_distance_function.py b/oceanmesh/signed_distance_function.py
index b0bcae0..ebeb1c4 100644
--- a/oceanmesh/signed_distance_function.py
+++ b/oceanmesh/signed_distance_function.py
@@ -69,13 +69,7 @@ def _generate_samples(bbox, dim, N):
     _yrange = (bbox[2] - 0.01, bbox[3] + 0.01)
     if dim == 2:
         points.append(
-            [
-                (
-                    random.uniform(*_xrange),
-                    random.uniform(*_yrange),
-                )
-                for i in range(N)
-            ]
+            [(random.uniform(*_xrange), random.uniform(*_yrange),) for i in range(N)]
         )
     elif dim == 3:
         _zrange = (bbox[4] - 0.01, bbox[5] + 0.01)
@@ -249,13 +243,7 @@ def _create_boubox(bbox):
     """Create a bounding box from domain extents `bbox`. Path orientation will be CCW."""
     xmin, xmax, ymin, ymax = bbox
     return np.array(
-        [
-            [xmin, ymin],
-            [xmax, ymin],
-            [xmax, ymax],
-            [xmin, ymax],
-            [xmin, ymin],
-        ],
+        [[xmin, ymin], [xmax, ymin], [xmax, ymax], [xmin, ymax], [xmin, ymin]],
         dtype=float,
     )
 
diff --git a/setup.py b/setup.py
index aa2be55..3c46783 100644
--- a/setup.py
+++ b/setup.py
@@ -2,11 +2,7 @@
 import sys
 
 from pybind11.setup_helpers import Pybind11Extension, build_ext
-from setuptools import setup, find_packages
-
-import versioneer
-
-sys.path.append(os.path.dirname(__file__))
+from setuptools import setup  # , find_packages
 
 import versioneer
 
@@ -25,16 +21,15 @@
     "oceanmesh/cpp/fast_geometry.cpp",
 ]
 
-import os
 if os.name == 'nt':
     home = os.environ['USERPROFILE'].replace('\\', '/')
     vcpkg = f'{home}/OceanMesh/vcpkg/installed/x64-windows'
     ext_modules = [
-        Pybind11Extension(loc, [fi], 
-            include_dirs=[f'{vcpkg}/include'], 
-            extra_link_args=[f'/LIBPATH:{vcpkg}/lib'], 
-            libraries=["gmp", "mpfr"]
-        )
+        Pybind11Extension(loc, [fi],
+                          include_dirs=[f'{vcpkg}/include'],
+                          extra_link_args=[f'/LIBPATH:{vcpkg}/lib'],
+                          libraries=["gmp", "mpfr"]
+                          )
         for fi, loc in zip(files, is_called)
     ]
 else:
@@ -54,4 +49,3 @@
         ext_modules=ext_modules,
         zip_safe=False,
     )
-
diff --git a/tests/test_bathymetric_gradient_function.py b/tests/test_bathymetric_gradient_function.py
index 8eb16a5..32225d6 100644
--- a/tests/test_bathymetric_gradient_function.py
+++ b/tests/test_bathymetric_gradient_function.py
@@ -61,10 +61,7 @@ def test_bathymetric_gradient_function():
     sdf = om.signed_distance_function(shoreline)
 
     edge_length1 = om.feature_sizing_function(
-        shoreline,
-        sdf,
-        max_edge_length=max_edge_length,
-        crs=EPSG,
+        shoreline, sdf, max_edge_length=max_edge_length, crs=EPSG,
     )
     edge_length2 = om.bathymetric_gradient_sizing_function(
         dem,
diff --git a/tests/test_circ_rect.py b/tests/test_circ_rect.py
index a0eee15..1fab6d4 100644
--- a/tests/test_circ_rect.py
+++ b/tests/test_circ_rect.py
@@ -1,9 +1,11 @@
-import os
 import logging
+import os
+
 import fiona
 import geopandas as gpd
 import matplotlib.pyplot as pyplot
 import numpy
+
 import oceanmesh
 
 shp0 = os.path.join(os.path.dirname(__file__), "ocean.shp")
diff --git a/tests/test_geodata.py b/tests/test_geodata.py
index 617add8..22ea5a5 100644
--- a/tests/test_geodata.py
+++ b/tests/test_geodata.py
@@ -1,6 +1,7 @@
 import os
 
 import pytest
+
 from oceanmesh import DEM, Region, Shoreline, edges
 
 fname = os.path.join(os.path.dirname(__file__), "GSHHS_i_L1.shp")
@@ -10,10 +11,7 @@
 
 @pytest.mark.parametrize(
     "boxes_h0",
-    [
-        ((166.0, 176.0, -48.0, -40.0), 0.01),
-        ((-74.0, -70.0, 35.0, 42.0), 0.005),
-    ],
+    [((166.0, 176.0, -48.0, -40.0), 0.01), ((-74.0, -70.0, 35.0, 42.0), 0.005)],
 )
 def test_shoreline(boxes_h0):
     """Read in a shapefile at different scales h0
@@ -30,14 +28,8 @@ def test_shoreline(boxes_h0):
 @pytest.mark.parametrize(
     "files_bboxes",
     [
-        (
-            dfname,
-            (-95.24, -95.21, 28.95, 29.00),
-        ),
-        (
-            tfname,
-            (-95.24, -95.21, 28.95, 29.00),
-        ),
+        (dfname, (-95.24, -95.21, 28.95, 29.00),),
+        (tfname, (-95.24, -95.21, 28.95, 29.00),),
     ],
 )
 def test_geodata(files_bboxes):
diff --git a/tests/test_grade.py b/tests/test_grade.py
index ec5deff..bc9756c 100644
--- a/tests/test_grade.py
+++ b/tests/test_grade.py
@@ -1,6 +1,8 @@
 import os
+
 import matplotlib.pyplot as plt
 import matplotlib.tri as tri
+
 import oceanmesh as om
 
 
diff --git a/tests/test_irregular_domain.py b/tests/test_irregular_domain.py
index abc0250..77384ad 100644
--- a/tests/test_irregular_domain.py
+++ b/tests/test_irregular_domain.py
@@ -1,7 +1,9 @@
 import os
-import matplotlib.tri as tri
+
 import matplotlib.pyplot as plt
+import matplotlib.tri as tri
 import numpy as np
+
 import oceanmesh as om
 
 
diff --git a/tests/test_mesh_generator_simple.py b/tests/test_mesh_generator_simple.py
index abee4dd..1122022 100644
--- a/tests/test_mesh_generator_simple.py
+++ b/tests/test_mesh_generator_simple.py
@@ -1,4 +1,5 @@
 import numpy as np
+
 from oceanmesh import generate_mesh, simp_vol
 
 
diff --git a/tests/test_multiscale.py b/tests/test_multiscale.py
index 66dfa36..b3ff24f 100644
--- a/tests/test_multiscale.py
+++ b/tests/test_multiscale.py
@@ -1,4 +1,5 @@
 import os
+
 import matplotlib.gridspec as gridspec
 import matplotlib.pyplot as plt
 import matplotlib.tri as tri

From 2a4543d5c1992bde32aa93a6dbdadac52262f627 Mon Sep 17 00:00:00 2001
From: Coat <alain.coat@actimar.fr>
Date: Sat, 25 Jun 2022 00:20:52 +0200
Subject: [PATCH 05/34] doc after black

---
 README.md | 42 +++++++++++++++++++++---------------------
 1 file changed, 21 insertions(+), 21 deletions(-)

diff --git a/README.md b/README.md
index c93c097..03969c0 100644
--- a/README.md
+++ b/README.md
@@ -139,11 +139,20 @@ Setting the region
 ```python
 import oceanmesh as om
 
-EPSG = 32619 # A Python int, dict, or str containing the CRS information (in this case UTM19N)
-bbox = (-70.29637, -43.56508, -69.65537, 43.88338) # the extent of the domain (can also be a multi-polygon delimited by rows of np.nan)
-extent = om.Region(extent=bbox, crs=4326) # set the region (the bbox is given here in EPSG:4326 or WGS84)
-extent = extent.transform_to(EPSG) # Now I transform to the desired EPSG (UTM19N)
-print(extent.bbox) # now the extents are in the desired CRS and can be passed to various functions later on
+EPSG = 32619  # A Python int, dict, or str containing the CRS information (in this case UTM19N)
+bbox = (
+    -70.29637,
+    -43.56508,
+    -69.65537,
+    43.88338,
+)  # the extent of the domain (can also be a multi-polygon delimited by rows of np.nan)
+extent = om.Region(
+    extent=bbox, crs=4326
+)  # set the region (the bbox is given here in EPSG:4326 or WGS84)
+extent = extent.transform_to(EPSG)  # Now I transform to the desired EPSG (UTM19N)
+print(
+    extent.bbox
+)  # now the extents are in the desired CRS and can be passed to various functions later on
 ```
 
 
@@ -209,9 +218,7 @@ fdem = "datasets/EastCoast.nc"
 # DEM is read into memory.
 # Note: the DEM will be projected to the desired CRS automatically.
 EPSG = 4326
-dem = om.DEM(
-    fdem, crs=EPSG
-)
+dem = om.DEM(fdem, crs=EPSG)
 dem.plot(
     xlabel="longitude (WGS84 degrees)",
     ylabel="latitude (WGS84 degrees)",
@@ -354,7 +361,9 @@ dem = om.DEM(fdem, 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)
-edge_length2 = om.wavelength_sizing_function(dem, wl=100, period=12.42*3600) # use the M2-tide period (in seconds)
+edge_length2 = om.wavelength_sizing_function(
+    dem, wl=100, period=12.42 * 3600
+)  # use the M2-tide period (in seconds)
 # 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)
@@ -392,10 +401,7 @@ shoreline = om.Shoreline(fname, extent.bbox, min_edge_length)
 sdf = om.signed_distance_function(shoreline)
 
 edge_length1 = om.feature_sizing_function(
-    shoreline,
-    sdf,
-    max_edge_length=max_edge_length,
-    crs=EPSG,
+    shoreline, sdf, max_edge_length=max_edge_length, crs=EPSG,
 )
 edge_length2 = om.bathymetric_gradient_sizing_function(
     dem,
@@ -478,10 +484,7 @@ points, cells = om.laplacian2(points, cells)
 
 # write the mesh with meshio
 meshio.write_points_cells(
-    "new_york.vtk",
-    points,
-    [("triangle", cells)],
-    file_format="vtk",
+    "new_york.vtk", points, [("triangle", cells)], file_format="vtk",
 )
 ```
 
@@ -528,10 +531,7 @@ sdf2 = om.signed_distance_function(s2)
 el2 = om.distance_sizing_function(s2)
 # Control the element size transition
 # from coarse to fine with the kwargs prefixed with `blend`
-points, cells = om.generate_multiscale_mesh(
-    [sdf1, sdf2],
-    [el1, el2],
-)
+points, cells = om.generate_multiscale_mesh([sdf1, sdf2], [el1, el2],)
 # remove degenerate mesh faces and other common problems in the mesh
 points, cells = om.make_mesh_boundaries_traversable(points, cells)
 # remove singly connected elements (elements connected to only one other element)

From 7e66133fb9dffce55d1530d6e883106d0697a925 Mon Sep 17 00:00:00 2001
From: alcoat <97431609+alcoat@users.noreply.github.com>
Date: Sat, 25 Jun 2022 14:48:33 +0200
Subject: [PATCH 06/34] Update testing.yml

Add python 3.10
---
 .github/workflows/testing.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml
index fafa637..35e9538 100644
--- a/.github/workflows/testing.yml
+++ b/.github/workflows/testing.yml
@@ -27,7 +27,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: [3.7, 3.8, 3.9]
+        python-version: [3.7, 3.8, 3.9, 3.10]
     steps:
     - uses: actions/checkout@v2
     - name: Set up Python ${{ matrix.python-version }}

From 21981ce52592e7e1037d67c9e25a785e10502bd8 Mon Sep 17 00:00:00 2001
From: alcoat <97431609+alcoat@users.noreply.github.com>
Date: Sat, 25 Jun 2022 14:55:37 +0200
Subject: [PATCH 07/34] Update testing.yml

Quoted version
---
 .github/workflows/testing.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml
index 35e9538..9ea35ff 100644
--- a/.github/workflows/testing.yml
+++ b/.github/workflows/testing.yml
@@ -27,7 +27,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: [3.7, 3.8, 3.9, 3.10]
+        python-version: [3.7, 3.8, 3.9, '3.10' ]
     steps:
     - uses: actions/checkout@v2
     - name: Set up Python ${{ matrix.python-version }}

From a5018b88319c28ca8b87fcc1eca7767763f78183 Mon Sep 17 00:00:00 2001
From: Coat <alain.coat@actimar.fr>
Date: Sat, 25 Jun 2022 16:27:18 +0200
Subject: [PATCH 08/34] upgrade black 22.3.0

---
 README.md                                   | 15 +++++++++---
 oceanmesh/edgefx.py                         | 15 ++++++++----
 oceanmesh/fix_mesh.py                       |  2 +-
 oceanmesh/geodata.py                        | 13 ++++++++---
 oceanmesh/idw.py                            |  2 +-
 oceanmesh/mesh_generator.py                 | 26 +++++++++++++++------
 oceanmesh/signed_distance_function.py       |  8 ++++++-
 tests/test_bathymetric_gradient_function.py |  5 +++-
 tests/test_geodata.py                       | 10 ++++++--
 9 files changed, 73 insertions(+), 23 deletions(-)

diff --git a/README.md b/README.md
index 03969c0..883f650 100644
--- a/README.md
+++ b/README.md
@@ -401,7 +401,10 @@ shoreline = om.Shoreline(fname, extent.bbox, min_edge_length)
 sdf = om.signed_distance_function(shoreline)
 
 edge_length1 = om.feature_sizing_function(
-    shoreline, sdf, max_edge_length=max_edge_length, crs=EPSG,
+    shoreline,
+    sdf,
+    max_edge_length=max_edge_length,
+    crs=EPSG,
 )
 edge_length2 = om.bathymetric_gradient_sizing_function(
     dem,
@@ -484,7 +487,10 @@ points, cells = om.laplacian2(points, cells)
 
 # write the mesh with meshio
 meshio.write_points_cells(
-    "new_york.vtk", points, [("triangle", cells)], file_format="vtk",
+    "new_york.vtk",
+    points,
+    [("triangle", cells)],
+    file_format="vtk",
 )
 ```
 
@@ -531,7 +537,10 @@ sdf2 = om.signed_distance_function(s2)
 el2 = om.distance_sizing_function(s2)
 # Control the element size transition
 # from coarse to fine with the kwargs prefixed with `blend`
-points, cells = om.generate_multiscale_mesh([sdf1, sdf2], [el1, el2],)
+points, cells = om.generate_multiscale_mesh(
+    [sdf1, sdf2],
+    [el1, el2],
+)
 # remove degenerate mesh faces and other common problems in the mesh
 points, cells = om.make_mesh_boundaries_traversable(points, cells)
 # remove singly connected elements (elements connected to only one other element)
diff --git a/oceanmesh/edgefx.py b/oceanmesh/edgefx.py
index 0b4b766..dd94353 100644
--- a/oceanmesh/edgefx.py
+++ b/oceanmesh/edgefx.py
@@ -128,7 +128,11 @@ def enforce_mesh_gradation(grid, gradation=0.15, crs=4326):
 
 
 def distance_sizing_function(
-    shoreline, rate=0.15, max_edge_length=None, coarsen=1, crs=4326,
+    shoreline,
+    rate=0.15,
+    max_edge_length=None,
+    coarsen=1,
+    crs=4326,
 ):
     """Mesh sizes that vary linearly at `rate` from coordinates in `obj`:Shoreline
     Parameters
@@ -284,7 +288,7 @@ def bathymetric_gradient_sizing_function(
         msg = f"The type_of_filter {type_of_filter} is not known and remains off"
         logger.info(msg)
         by, bx = _earth_gradient(tmpz, dy, dx)  # get slope in x and y directions
-        bs = np.sqrt(bx ** 2 + by ** 2)  # get overall slope
+        bs = np.sqrt(bx**2 + by**2)  # get overall slope
 
     # Calculating the slope function
     eps = 1e-10  # small number to approximate derivative
@@ -439,7 +443,7 @@ def rossby_radius_filter(tmpz, bbox, grid_details, coords, rbfilt, barot):
             by, bx = _earth_gradient(
                 tmpz_ft, dy, dx
             )  # [n2s:n2e]) # get slope in x and y directions
-            tempbs = np.sqrt(bx ** 2 + by ** 2)  # get overall slope
+            tempbs = np.sqrt(bx**2 + by**2)  # get overall slope
 
             bst[rosb == edges[i]] = tempbs[rosb == edges[i]]
 
@@ -628,7 +632,10 @@ def wavelength_sizing_function(
 
 
 def multiscale_sizing_function(
-    list_of_grids, p=3, nnear=28, blend_width=1000,
+    list_of_grids,
+    p=3,
+    nnear=28,
+    blend_width=1000,
 ):
     """Given a list of mesh size functions in a hierarchy
     w.r.t. to minimum mesh size (largest -> smallest),
diff --git a/oceanmesh/fix_mesh.py b/oceanmesh/fix_mesh.py
index 3f45044..6b3c36d 100644
--- a/oceanmesh/fix_mesh.py
+++ b/oceanmesh/fix_mesh.py
@@ -13,7 +13,7 @@ def simp_qual(p, t):
     assert p.ndim == 2 and t.ndim == 2 and p.shape[1] + 1 == t.shape[1]
 
     def length(p1):
-        return np.sqrt((p1 ** 2).sum(1))
+        return np.sqrt((p1**2).sum(1))
 
     a = length(p[t[:, 1]] - p[t[:, 0]])
     b = length(p[t[:, 2]] - p[t[:, 0]])
diff --git a/oceanmesh/geodata.py b/oceanmesh/geodata.py
index 2df1723..5d30d21 100644
--- a/oceanmesh/geodata.py
+++ b/oceanmesh/geodata.py
@@ -96,7 +96,9 @@ def _densify(poly, maxdiff, bbox, radius=0):
         else:
             ni = int(ni)
             icoords = _create_ranges(
-                np.array([lat[i], lon[i]]), np.array([lat[i + 1], lon[i + 1]]), ni + 2,
+                np.array([lat[i], lon[i]]),
+                np.array([lat[i + 1], lon[i + 1]]),
+                ni + 2,
             )
             latout[n : n + ni + 1] = icoords[0, : ni + 1]
             lonout[n : n + ni + 1] = icoords[1, : ni + 1]
@@ -135,7 +137,7 @@ def _classify_shoreline(bbox, boubox, polys, h0, minimum_area_mult):
     """
     logger.debug("Entering:_classify_shoreline")
 
-    _AREAMIN = minimum_area_mult * h0 ** 2
+    _AREAMIN = minimum_area_mult * h0**2
 
     if len(boubox) == 0:
         boubox = _create_boubox(bbox)
@@ -604,7 +606,12 @@ def plot(
 
         xmin, xmax, ymin, ymax = self.bbox
         rect = plt.Rectangle(
-            (xmin, ymin), xmax - xmin, ymax - ymin, fill=None, hatch="////", alpha=0.2,
+            (xmin, ymin),
+            xmax - xmin,
+            ymax - ymin,
+            fill=None,
+            hatch="////",
+            alpha=0.2,
         )
 
         border = 0.10 * (xmax - xmin)
diff --git a/oceanmesh/idw.py b/oceanmesh/idw.py
index 9705f9a..f71b8fa 100644
--- a/oceanmesh/idw.py
+++ b/oceanmesh/idw.py
@@ -92,7 +92,7 @@ def __call__(self, q, nnear=6, eps=0, p=1, weights=None):
             elif dist[0] < 1e-10:
                 wz = self.z[ix[0]]
             else:  # weight z s by 1/dist --
-                w = 1 / dist ** p
+                w = 1 / dist**p
                 if weights is not None:
                     w *= weights[ix]  # >= 0
                 w /= np.sum(w)
diff --git a/oceanmesh/mesh_generator.py b/oceanmesh/mesh_generator.py
index da3f07c..1b7f042 100644
--- a/oceanmesh/mesh_generator.py
+++ b/oceanmesh/mesh_generator.py
@@ -214,7 +214,14 @@ def generate_mesh(domain, edge_length, **kwargs):
     lock_boundary = opts["lock_boundary"]
 
     if opts["points"] is None:
-        p = _generate_initial_points(min_edge_length, geps, bbox, fh, fd, pfix,)
+        p = _generate_initial_points(
+            min_edge_length,
+            geps,
+            bbox,
+            fh,
+            fd,
+            pfix,
+        )
     else:
         p = opts["points"]
 
@@ -275,7 +282,7 @@ def generate_mesh(domain, edge_length, **kwargs):
         p = _project_points_back(p, fd, deps)
 
         # Show the user some progress so they know something is happening
-        maxdp = delta_t * np.sqrt((Ftot ** 2).sum(1)).max()
+        maxdp = delta_t * np.sqrt((Ftot**2).sum(1)).max()
 
         logger.info(
             f"Iteration #{count+1}, max movement is {maxdp}, there are {len(p)} vertices and {len(t)}"
@@ -325,10 +332,10 @@ def _compute_forces(p, t, fh, min_edge_length, L0mult):
     N = p.shape[0]
     bars = _get_bars(t)
     barvec = p[bars[:, 0]] - p[bars[:, 1]]  # List of bar vectors
-    L = np.sqrt((barvec ** 2).sum(1))  # L = Bar lengths
+    L = np.sqrt((barvec**2).sum(1))  # L = Bar lengths
     L[L == 0] = np.finfo(float).eps
     hbars = fh(p[bars].sum(1) / 2)
-    L0 = hbars * L0mult * ((L ** 2).sum() / (hbars ** 2).sum()) ** (1.0 / 2)
+    L0 = hbars * L0mult * ((L**2).sum() / (hbars**2).sum()) ** (1.0 / 2)
     F = L0 - L
     F[F < 0] = 0  # Bar forces (scalars)
     Fvec = (
@@ -416,7 +423,7 @@ def _deps_vec(i):
             # take the solely ones outside domain
             dgrads = [(fd(p + _deps_vec(i)) - d) / deps for i in range(2)]
             dgrads = list(np.array(dgrads)[:, ix])
-        dgrad2 = sum(dgrad ** 2 for dgrad in dgrads)
+        dgrad2 = sum(dgrad**2 for dgrad in dgrads)
         dgrad2 = np.where(dgrad2 < deps, deps, dgrad2)
         p[ix] -= (d[ix] * np.vstack(dgrads) / dgrad2).T  # Project
     return p
@@ -430,9 +437,14 @@ def _generate_initial_points(min_edge_length, geps, bbox, fh, fd, pfix):
     p = p.reshape(2, -1).T
     r0 = fh(p)
     r0m = np.min(r0[r0 >= min_edge_length])
-    p = p[np.random.rand(p.shape[0]) < r0m ** 2 / r0 ** 2]
+    p = p[np.random.rand(p.shape[0]) < r0m**2 / r0**2]
     p = p[fd(p) < geps]  # Keep only d<0 points
-    return np.vstack((pfix, p,))
+    return np.vstack(
+        (
+            pfix,
+            p,
+        )
+    )
 
 
 def _dist(p1, p2):
diff --git a/oceanmesh/signed_distance_function.py b/oceanmesh/signed_distance_function.py
index ebeb1c4..e58682d 100644
--- a/oceanmesh/signed_distance_function.py
+++ b/oceanmesh/signed_distance_function.py
@@ -69,7 +69,13 @@ def _generate_samples(bbox, dim, N):
     _yrange = (bbox[2] - 0.01, bbox[3] + 0.01)
     if dim == 2:
         points.append(
-            [(random.uniform(*_xrange), random.uniform(*_yrange),) for i in range(N)]
+            [
+                (
+                    random.uniform(*_xrange),
+                    random.uniform(*_yrange),
+                )
+                for i in range(N)
+            ]
         )
     elif dim == 3:
         _zrange = (bbox[4] - 0.01, bbox[5] + 0.01)
diff --git a/tests/test_bathymetric_gradient_function.py b/tests/test_bathymetric_gradient_function.py
index 32225d6..8eb16a5 100644
--- a/tests/test_bathymetric_gradient_function.py
+++ b/tests/test_bathymetric_gradient_function.py
@@ -61,7 +61,10 @@ def test_bathymetric_gradient_function():
     sdf = om.signed_distance_function(shoreline)
 
     edge_length1 = om.feature_sizing_function(
-        shoreline, sdf, max_edge_length=max_edge_length, crs=EPSG,
+        shoreline,
+        sdf,
+        max_edge_length=max_edge_length,
+        crs=EPSG,
     )
     edge_length2 = om.bathymetric_gradient_sizing_function(
         dem,
diff --git a/tests/test_geodata.py b/tests/test_geodata.py
index 22ea5a5..f8f2113 100644
--- a/tests/test_geodata.py
+++ b/tests/test_geodata.py
@@ -28,8 +28,14 @@ def test_shoreline(boxes_h0):
 @pytest.mark.parametrize(
     "files_bboxes",
     [
-        (dfname, (-95.24, -95.21, 28.95, 29.00),),
-        (tfname, (-95.24, -95.21, 28.95, 29.00),),
+        (
+            dfname,
+            (-95.24, -95.21, 28.95, 29.00),
+        ),
+        (
+            tfname,
+            (-95.24, -95.21, 28.95, 29.00),
+        ),
     ],
 )
 def test_geodata(files_bboxes):

From a7d997724f0bbe60361622d9b04d1b83a54c3087 Mon Sep 17 00:00:00 2001
From: Coat <alain.coat@actimar.fr>
Date: Sun, 26 Jun 2022 13:29:13 +0200
Subject: [PATCH 09/34] try rasterio>=1.3a1 for python 3.10

---
 setup.cfg | 1 +
 1 file changed, 1 insertion(+)

diff --git a/setup.cfg b/setup.cfg
index 79ea8d8..7fa8341 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -30,6 +30,7 @@ install_requires =
         geopandas
         shapely
         matplotlib
+        rasterio>=1.3a1
         rioxarray
         scikit-fmm
         scikit-image

From 5a3d00913f9843e2730345d1c0c1a9e527b9f6c7 Mon Sep 17 00:00:00 2001
From: Coat <alain.coat@actimar.fr>
Date: Sun, 26 Jun 2022 13:39:19 +0200
Subject: [PATCH 10/34] apply black to setup.py

---
 Makefile |  2 +-
 setup.py | 18 ++++++++++--------
 2 files changed, 11 insertions(+), 9 deletions(-)

diff --git a/Makefile b/Makefile
index e5be404..4504206 100644
--- a/Makefile
+++ b/Makefile
@@ -11,5 +11,5 @@ clean:
 
 format:
 	isort -rc oceanmesh/ tests/*.py
-	black oceanmesh/ tests/*.py
+	black setup.py oceanmesh/ tests/*.py
 	blacken-docs README.md
diff --git a/setup.py b/setup.py
index 3c46783..a19dafa 100644
--- a/setup.py
+++ b/setup.py
@@ -21,15 +21,17 @@
     "oceanmesh/cpp/fast_geometry.cpp",
 ]
 
-if os.name == 'nt':
-    home = os.environ['USERPROFILE'].replace('\\', '/')
-    vcpkg = f'{home}/OceanMesh/vcpkg/installed/x64-windows'
+if os.name == "nt":
+    home = os.environ["USERPROFILE"].replace("\\", "/")
+    vcpkg = f"{home}/OceanMesh/vcpkg/installed/x64-windows"
     ext_modules = [
-        Pybind11Extension(loc, [fi],
-                          include_dirs=[f'{vcpkg}/include'],
-                          extra_link_args=[f'/LIBPATH:{vcpkg}/lib'],
-                          libraries=["gmp", "mpfr"]
-                          )
+        Pybind11Extension(
+            loc,
+            [fi],
+            include_dirs=[f"{vcpkg}/include"],
+            extra_link_args=[f"/LIBPATH:{vcpkg}/lib"],
+            libraries=["gmp", "mpfr"],
+        )
         for fi, loc in zip(files, is_called)
     ]
 else:

From f41fd71528dfd76c338c052bc4d4a5a722165790 Mon Sep 17 00:00:00 2001
From: Keith Roberts <kroberts@baird.com>
Date: Sun, 26 Jun 2022 13:24:16 -0400
Subject: [PATCH 11/34] Update README.md

* Update CI badge to point to GitHub workflows.
---
 README.md | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/README.md b/README.md
index 883f650..1882184 100644
--- a/README.md
+++ b/README.md
@@ -2,10 +2,11 @@ oceanmesh: Automatic coastal ocean mesh generation
 =====================================================
 :ocean: :cyclone:
 
-[![CircleCI](https://circleci.com/gh/circleci/circleci-docs.svg?style=svg)](https://circleci.com/gh/CHLNDDEV/oceanmesh)
+[![Tests](https://github.com/CHLNDDEV/oceanmesh/actions/workflows/testing.yml/badge.svg)](https://github.com/CHLNDDEV/oceanmesh/actions/workflows/testing.yml)
+
 [![CodeCov](https://codecov.io/gh/CHLNDDEV/oceanmesh/branch/master/graph/badge.svg)](https://codecov.io/gh/CHLNDDEV/oceanmesh)
 
-Coastal ocean mesh generation from ESRI Shapefiles and digital elevation models.
+Coastal ocean mesh generation from vector and raster GIS data.
 
 Table of contents
 =================

From b4c1c8133983a40e99c68a926bec85e701c7d654 Mon Sep 17 00:00:00 2001
From: Coat <alain.coat@actimar.fr>
Date: Tue, 28 Jun 2022 22:04:50 +0200
Subject: [PATCH 12/34] remove duplicate points from Shoreline

---
 oceanmesh/geodata.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/oceanmesh/geodata.py b/oceanmesh/geodata.py
index 5d30d21..5ac7008 100644
--- a/oceanmesh/geodata.py
+++ b/oceanmesh/geodata.py
@@ -424,6 +424,13 @@ def _is_overlapping(bbox1, bbox2):
     return x1min < x2max and x2min < x1max and y1min < y2max and y2min < y1max
 
 
+def remove_dup(arr: np.ndarray):
+    """Remove duplicate element from np.ndarray"""
+    result = np.concatenate((arr[np.nonzero(np.diff(arr))[0]], [arr[-1]]))
+
+    return result
+
+
 class Shoreline(Region):
     """
     The shoreline class extends :class:`Region` to store data
@@ -567,6 +574,7 @@ def _read(self):
                     poly = np.asarray(g.coords)
                 else:  # a polygon
                     poly = np.asarray(g.exterior.coords.xy).T
+                poly = remove_dup(poly)
                 polys.append(np.row_stack((poly, delimiter)))
 
         if len(polys) == 0:

From 169930e4273e12f0e3d365991c6b3b1a2df334b8 Mon Sep 17 00:00:00 2001
From: Coat <alain.coat@actimar.fr>
Date: Wed, 29 Jun 2022 10:21:19 +0200
Subject: [PATCH 13/34] better nan removal

---
 oceanmesh/signed_distance_function.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/oceanmesh/signed_distance_function.py b/oceanmesh/signed_distance_function.py
index e58682d..89be469 100644
--- a/oceanmesh/signed_distance_function.py
+++ b/oceanmesh/signed_distance_function.py
@@ -199,14 +199,14 @@ def signed_distance_function(shoreline, invert=False):
     )
     e = edges.get_poly_edges(poly)
 
-    boubox = shoreline.boubox
-    e_box = edges.get_poly_edges(boubox)
+    boubox = np.nan_to_num(shoreline.boubox)
+    e_box = edges.get_poly_edges(shoreline.boubox)
 
     def func(x):
         # Initialize d with some positive number larger than geps
         dist = np.zeros(len(x)) + 1.0
         # are points inside the boubox?
-        in_boubox, _ = inpoly2(x, np.nan_to_num(boubox), e_box)
+        in_boubox, _ = inpoly2(x, boubox, e_box)
         # are points inside the shoreline?
         in_shoreline, _ = inpoly2(x, np.nan_to_num(poly), e)
         # compute dist to shoreline

From 78b11938b6c587757c6bbe933d2482cba415b64a Mon Sep 17 00:00:00 2001
From: Coat <alain.coat@actimar.fr>
Date: Wed, 29 Jun 2022 10:22:05 +0200
Subject: [PATCH 14/34] more explict edges buildind

---
 oceanmesh/edges.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/oceanmesh/edges.py b/oceanmesh/edges.py
index 563495d..e49bc4d 100644
--- a/oceanmesh/edges.py
+++ b/oceanmesh/edges.py
@@ -27,10 +27,12 @@ def get_poly_edges(poly):
 
     edges = []
     for s in range(len(ix) - 1):
-        col1 = np.arange(ix[s] + 1, ix[s + 1] - 2)
-        col2 = np.arange(ix[s] + 2, ix[s + 1] - 1)
+        ix_start = ix[s] + 1
+        ix_end = ix[s + 1] - 1
+        col1 = np.arange(ix_start, ix_end - 1)
+        col2 = np.arange(ix_start + 1, ix_end)
         tmp = np.vstack((col1, col2)).T
-        tmp = np.append(tmp, [[ix[s + 1] - 1, ix[s] + 1]], axis=0)
+        tmp = np.append(tmp, [[ix_end, ix_start]], axis=0)
         edges.append(tmp)
     return np.concatenate(edges, axis=0)
 

From 1afbfa6c82510361141665a3f80d140709ca42c7 Mon Sep 17 00:00:00 2001
From: Keith Roberts <krober@usp.br>
Date: Wed, 29 Jun 2022 08:12:39 -0400
Subject: [PATCH 15/34] remove stray print statement

---
 oceanmesh/mesh_generator.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/oceanmesh/mesh_generator.py b/oceanmesh/mesh_generator.py
index 1b7f042..d350dcd 100644
--- a/oceanmesh/mesh_generator.py
+++ b/oceanmesh/mesh_generator.py
@@ -197,7 +197,6 @@ def generate_mesh(domain, edge_length, **kwargs):
     _check_bbox(bbox)
     bbox = np.array(bbox).reshape(-1, 2)
 
-    print(min_edge_length)
     assert min_edge_length > 0, "`min_edge_length` must be > 0"
 
     assert opts["max_iter"] > 0, "`max_iter` must be > 0"

From 12ce85f043e471765facb94cee1ac67c914ceee7 Mon Sep 17 00:00:00 2001
From: Alain Coat <97431609+alcoat@users.noreply.github.com>
Date: Wed, 8 Jan 2025 19:38:08 +0100
Subject: [PATCH 16/34] update python version

remove 3.7 and 3.8, add 3.11 and 3.12
---
 .github/workflows/testing.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml
index 27f13ea..df8f3d8 100644
--- a/.github/workflows/testing.yml
+++ b/.github/workflows/testing.yml
@@ -27,7 +27,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: [3.7, 3.8, 3.9, '3.10' ]
+        python-version: [ 3.9, '3.10', '3.11', '3.12']
     steps:
     - uses: actions/checkout@v2
     - name: Set up Python ${{ matrix.python-version }}

From 27511a1e44a1dc73577e2eb5ae1c7c04bfdd9893 Mon Sep 17 00:00:00 2001
From: Alain Coat <97431609+alcoat@users.noreply.github.com>
Date: Wed, 8 Jan 2025 19:39:52 +0100
Subject: [PATCH 17/34] Update testing.yml

---
 .github/workflows/testing.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml
index df8f3d8..f61ab81 100644
--- a/.github/workflows/testing.yml
+++ b/.github/workflows/testing.yml
@@ -12,7 +12,7 @@ jobs:
     steps:
     - uses: actions/setup-python@v2
       with:
-        python-version: "3.7"
+        python-version: "3.9"
     - uses: actions/checkout@v2
     - name: Lint with flake8
       run: |

From 3f02702334d29294f1da61e1c4531edd487c5974 Mon Sep 17 00:00:00 2001
From: Alain Coat <97431609+alcoat@users.noreply.github.com>
Date: Wed, 8 Jan 2025 19:43:48 +0100
Subject: [PATCH 18/34] Add fiona

Fiona needed as geopandas extension
---
 setup.cfg | 1 +
 1 file changed, 1 insertion(+)

diff --git a/setup.cfg b/setup.cfg
index 7fa8341..6d2d3ce 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -28,6 +28,7 @@ install_requires =
         pyproj
         scipy
         geopandas
+        fiona
         shapely
         matplotlib
         rasterio>=1.3a1

From 82a28691d5308e51b5af6c913aa72e0f88cba876 Mon Sep 17 00:00:00 2001
From: Alain Coat <97431609+alcoat@users.noreply.github.com>
Date: Wed, 8 Jan 2025 19:56:02 +0100
Subject: [PATCH 19/34] keep python 3.8

---
 .github/workflows/testing.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml
index f61ab81..fc707fc 100644
--- a/.github/workflows/testing.yml
+++ b/.github/workflows/testing.yml
@@ -27,7 +27,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: [ 3.9, '3.10', '3.11', '3.12']
+        python-version: [ 3.8, 3.9, '3.10', '3.11', '3.12']
     steps:
     - uses: actions/checkout@v2
     - name: Set up Python ${{ matrix.python-version }}

From 1dc56c634fcf0b675623f77fd7535858f0222c07 Mon Sep 17 00:00:00 2001
From: Coat <alain.coat@actimar.fr>
Date: Wed, 8 Jan 2025 22:52:24 +0100
Subject: [PATCH 20/34] fix black format

---
 oceanmesh/geodata.py | 6 +++---
 oceanmesh/idw.py     | 1 +
 2 files changed, 4 insertions(+), 3 deletions(-)

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

From 8122f62ebdf3f984831a5b546407c8e65216434f Mon Sep 17 00:00:00 2001
From: Coat <alain.coat@actimar.fr>
Date: Wed, 8 Jan 2025 23:20:57 +0100
Subject: [PATCH 21/34] numpy deprecated

---
 oceanmesh/fix_mesh.py | 2 +-
 oceanmesh/geodata.py  | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/oceanmesh/fix_mesh.py b/oceanmesh/fix_mesh.py
index 8e3be27..a7984d1 100644
--- a/oceanmesh/fix_mesh.py
+++ b/oceanmesh/fix_mesh.py
@@ -96,7 +96,7 @@ def unique_rows(A, return_index=False, return_inverse=False):
 
     orig_dtype = A.dtype
     ncolumns = A.shape[1]
-    dtype = np.dtype((np.character, orig_dtype.itemsize * ncolumns))
+    dtype = np.dtype(('S1', orig_dtype.itemsize * ncolumns))
     B, I, J = np.unique(A.view(dtype), return_index=True, return_inverse=True)
 
     B = B.view(orig_dtype).reshape((-1, ncolumns), order="C")
diff --git a/oceanmesh/geodata.py b/oceanmesh/geodata.py
index 9f0a514..00af47d 100644
--- a/oceanmesh/geodata.py
+++ b/oceanmesh/geodata.py
@@ -669,7 +669,7 @@ def _read(self):
                     raise ValueError(f"Unsupported geometry type: {g.geom_type}")
 
                 poly = remove_dup(poly)
-                polys.append(np.row_stack((poly, delimiter)))
+                polys.append(np.vstack((poly, delimiter)))
 
         if len(polys) == 0:
             raise ValueError("Shoreline data does not intersect with bbox")

From dbe55ccd2c6cc3cb5a99b7323c2b0142b0879b64 Mon Sep 17 00:00:00 2001
From: Coat <alain.coat@actimar.fr>
Date: Wed, 8 Jan 2025 23:22:35 +0100
Subject: [PATCH 22/34] fix black

---
 oceanmesh/fix_mesh.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/oceanmesh/fix_mesh.py b/oceanmesh/fix_mesh.py
index a7984d1..3dba736 100644
--- a/oceanmesh/fix_mesh.py
+++ b/oceanmesh/fix_mesh.py
@@ -96,7 +96,7 @@ def unique_rows(A, return_index=False, return_inverse=False):
 
     orig_dtype = A.dtype
     ncolumns = A.shape[1]
-    dtype = np.dtype(('S1', orig_dtype.itemsize * ncolumns))
+    dtype = np.dtype(("S1", orig_dtype.itemsize * ncolumns))
     B, I, J = np.unique(A.view(dtype), return_index=True, return_inverse=True)
 
     B = B.view(orig_dtype).reshape((-1, ncolumns), order="C")

From 9e8734c253b705de5cad29cdb46e5f915932e62a Mon Sep 17 00:00:00 2001
From: alcoat <alain.coat@actimar.fr>
Date: Sat, 11 Jan 2025 23:34:12 +0100
Subject: [PATCH 23/34] back to np.character

---
 oceanmesh/edgefx.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/oceanmesh/edgefx.py b/oceanmesh/edgefx.py
index d2ebba7..65a330a 100644
--- a/oceanmesh/edgefx.py
+++ b/oceanmesh/edgefx.py
@@ -7,11 +7,11 @@
 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 0dc55d092c1a4df4d716a6327781a7496574ee1e Mon Sep 17 00:00:00 2001
From: alcoat <alain.coat@actimar.fr>
Date: Sat, 11 Jan 2025 23:34:47 +0100
Subject: [PATCH 24/34] back to np.character

---
 oceanmesh/fix_mesh.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/oceanmesh/fix_mesh.py b/oceanmesh/fix_mesh.py
index 3dba736..8e3be27 100644
--- a/oceanmesh/fix_mesh.py
+++ b/oceanmesh/fix_mesh.py
@@ -96,7 +96,7 @@ def unique_rows(A, return_index=False, return_inverse=False):
 
     orig_dtype = A.dtype
     ncolumns = A.shape[1]
-    dtype = np.dtype(("S1", orig_dtype.itemsize * ncolumns))
+    dtype = np.dtype((np.character, orig_dtype.itemsize * ncolumns))
     B, I, J = np.unique(A.view(dtype), return_index=True, return_inverse=True)
 
     B = B.view(orig_dtype).reshape((-1, ncolumns), order="C")

From 145c99c4fd85b816d383f27cb91f7befc123bef7 Mon Sep 17 00:00:00 2001
From: alcoat <alain.coat@actimar.fr>
Date: Sat, 11 Jan 2025 23:35:52 +0100
Subject: [PATCH 25/34] black

---
 oceanmesh/mesh_generator.py                 | 1 -
 tests/test_bathymetric_gradient_function.py | 1 +
 2 files changed, 1 insertion(+), 1 deletion(-)

diff --git a/oceanmesh/mesh_generator.py b/oceanmesh/mesh_generator.py
index 06ae670..889cbbf 100644
--- a/oceanmesh/mesh_generator.py
+++ b/oceanmesh/mesh_generator.py
@@ -7,7 +7,6 @@
 import matplotlib.tri as tri
 import numpy as np
 import scipy.sparse as spsparse
-
 from _delaunay_class import DelaunayTriangulation as DT
 from _fast_geometry import unique_edges
 
diff --git a/tests/test_bathymetric_gradient_function.py b/tests/test_bathymetric_gradient_function.py
index 0bb0d72..c91c23d 100644
--- a/tests/test_bathymetric_gradient_function.py
+++ b/tests/test_bathymetric_gradient_function.py
@@ -1,6 +1,7 @@
 import os
 
 import pytest
+
 import oceanmesh as om
 
 fname = os.path.join(os.path.dirname(__file__), "GSHHS_i_L1.shp")

From ca9726d90041aea450053caec7ec86d29e0a8963 Mon Sep 17 00:00:00 2001
From: alcoat <alain.coat@actimar.fr>
Date: Sun, 12 Jan 2025 12:06:07 +0100
Subject: [PATCH 26/34] change dtype

---
 oceanmesh/fix_mesh.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/oceanmesh/fix_mesh.py b/oceanmesh/fix_mesh.py
index 8e3be27..df0c5ef 100644
--- a/oceanmesh/fix_mesh.py
+++ b/oceanmesh/fix_mesh.py
@@ -96,7 +96,7 @@ def unique_rows(A, return_index=False, return_inverse=False):
 
     orig_dtype = A.dtype
     ncolumns = A.shape[1]
-    dtype = np.dtype((np.character, orig_dtype.itemsize * ncolumns))
+    dtype = np.dtype((f"S{orig_dtype.itemsize * ncolumns}"))
     B, I, J = np.unique(A.view(dtype), return_index=True, return_inverse=True)
 
     B = B.view(orig_dtype).reshape((-1, ncolumns), order="C")

From 9d6c361704794a0a362daaafb9d2a47afbf51039 Mon Sep 17 00:00:00 2001
From: alcoat <alain.coat@actimar.fr>
Date: Sun, 12 Jan 2025 12:06:28 +0100
Subject: [PATCH 27/34] upgrade target version

---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index bf5a947..f4c6b52 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -3,4 +3,4 @@ requires = ["setuptools>=42", "wheel", "pybind11>=2.6.0", "versioneer-518"]
 build-backend = "setuptools.build_meta"
 
 [tool.black]
-target-version = ['py37']
\ No newline at end of file
+target-version = ['py39']

From 3f35bb6f8dce6bada182331f53fc8039551db9e7 Mon Sep 17 00:00:00 2001
From: alcoat <alain.coat@actimar.fr>
Date: Sun, 12 Jan 2025 12:57:08 +0100
Subject: [PATCH 28/34] remove depracated np.matrix

---
 oceanmesh/clean.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/oceanmesh/clean.py b/oceanmesh/clean.py
index d1d4f98..ec31ad0 100644
--- a/oceanmesh/clean.py
+++ b/oceanmesh/clean.py
@@ -471,7 +471,7 @@ def _closest_node(node, nodes):
     L[L < eps] = eps
     L = L[:, None]
     for it in range(max_iter):
-        pnew = np.divide(S * np.matrix(vertices), np.hstack((W, W)))
+        pnew = np.divide(S @ np.array(vertices), np.hstack((W, W)))
         pnew[bnd, :] = vertices[bnd, :]
         vertices = pnew
         Lnew = np.sqrt(

From a49642093904e8e1932b34c04d265b5bfe8ae639 Mon Sep 17 00:00:00 2001
From: alcoat <alain.coat@actimar.fr>
Date: Sun, 12 Jan 2025 13:51:13 +0100
Subject: [PATCH 29/34] numpy 2 compat

---
 oceanmesh/fix_mesh.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/oceanmesh/fix_mesh.py b/oceanmesh/fix_mesh.py
index df0c5ef..da39ea5 100644
--- a/oceanmesh/fix_mesh.py
+++ b/oceanmesh/fix_mesh.py
@@ -98,6 +98,8 @@ def unique_rows(A, return_index=False, return_inverse=False):
     ncolumns = A.shape[1]
     dtype = np.dtype((f"S{orig_dtype.itemsize * ncolumns}"))
     B, I, J = np.unique(A.view(dtype), return_index=True, return_inverse=True)
+    # NUMPY 2 compatibility
+    J = J.reshape(-1)
 
     B = B.view(orig_dtype).reshape((-1, ncolumns), order="C")
 

From add31275a1585814678e9df253fba901d375d5ae Mon Sep 17 00:00:00 2001
From: alcoat <alain.coat@actimar.fr>
Date: Sun, 12 Jan 2025 14:15:56 +0100
Subject: [PATCH 30/34] adjust python3.8 requirements

---
 setup.py | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/setup.py b/setup.py
index a19dafa..f694ecd 100644
--- a/setup.py
+++ b/setup.py
@@ -1,5 +1,6 @@
 import os
 import sys
+import configparser
 
 from pybind11.setup_helpers import Pybind11Extension, build_ext
 from setuptools import setup  # , find_packages
@@ -44,6 +45,19 @@
 cmdclass = versioneer.get_cmdclass()
 cmdclass.update({"build_ext": build_ext})
 
+def get_requirements():
+    """
+    Fix
+    """
+
+    config = configparser.ConfigParser()
+    config.read("setup.cfg")
+    requirements = config["options"]["install_requires"].split()
+
+    if sys.version_info <= (3,8):
+        requirements.append("fiona<=1.9")
+    return requirements
+
 if __name__ == "__main__":
     setup(
         cmdclass=cmdclass,

From 28e2009795227edcfac776c02149dfe4642bdce9 Mon Sep 17 00:00:00 2001
From: alcoat <alain.coat@actimar.fr>
Date: Sun, 12 Jan 2025 14:17:59 +0100
Subject: [PATCH 31/34] black

---
 setup.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/setup.py b/setup.py
index f694ecd..d546b70 100644
--- a/setup.py
+++ b/setup.py
@@ -45,6 +45,7 @@
 cmdclass = versioneer.get_cmdclass()
 cmdclass.update({"build_ext": build_ext})
 
+
 def get_requirements():
     """
     Fix
@@ -54,10 +55,11 @@ def get_requirements():
     config.read("setup.cfg")
     requirements = config["options"]["install_requires"].split()
 
-    if sys.version_info <= (3,8):
+    if sys.version_info <= (3, 8):
         requirements.append("fiona<=1.9")
     return requirements
 
+
 if __name__ == "__main__":
     setup(
         cmdclass=cmdclass,

From 68a379b76dd7a3c5a397f9de5e4e6de84d5a27c0 Mon Sep 17 00:00:00 2001
From: alcoat <alain.coat@actimar.fr>
Date: Sun, 12 Jan 2025 14:22:12 +0100
Subject: [PATCH 32/34] adjust requirements

---
 setup.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/setup.py b/setup.py
index d546b70..b7b1f3b 100644
--- a/setup.py
+++ b/setup.py
@@ -62,6 +62,7 @@ def get_requirements():
 
 if __name__ == "__main__":
     setup(
+        install_requires=get_requirements(),
         cmdclass=cmdclass,
         version=versioneer.get_version(),
         ext_modules=ext_modules,

From 7b91a5353d0ed1385de7aaecc9f123a30cde203e Mon Sep 17 00:00:00 2001
From: alcoat <alain.coat@actimar.fr>
Date: Sun, 12 Jan 2025 14:39:26 +0100
Subject: [PATCH 33/34] fix requirements

---
 setup.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/setup.py b/setup.py
index b7b1f3b..c271861 100644
--- a/setup.py
+++ b/setup.py
@@ -55,8 +55,10 @@ def get_requirements():
     config.read("setup.cfg")
     requirements = config["options"]["install_requires"].split()
 
-    if sys.version_info <= (3, 8):
-        requirements.append("fiona<=1.9")
+    if sys.version_info < (3, 9):
+        requirements.remove("fiona")
+        requirements.append("fiona<1.10")
+
     return requirements
 
 

From 214f8cdef043e9008e829a144c1171a22bfe69a7 Mon Sep 17 00:00:00 2001
From: Alain Coat <97431609+alcoat@users.noreply.github.com>
Date: Sun, 12 Jan 2025 19:42:56 +0100
Subject: [PATCH 34/34] Remove python 3.12

Action failed with python 3.12, however it works
---
 .github/workflows/testing.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml
index fc707fc..055cf65 100644
--- a/.github/workflows/testing.yml
+++ b/.github/workflows/testing.yml
@@ -27,7 +27,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: [ 3.8, 3.9, '3.10', '3.11', '3.12']
+        python-version: [ 3.8, 3.9, '3.10', '3.11']
     steps:
     - uses: actions/checkout@v2
     - name: Set up Python ${{ matrix.python-version }}