diff --git a/pyposeidon/mjigsaw.py b/pyposeidon/mjigsaw.py
index b19cda4..3009894 100644
--- a/pyposeidon/mjigsaw.py
+++ b/pyposeidon/mjigsaw.py
@@ -243,13 +243,13 @@ def read_msh(filename, **kwargs):
 
     A, idxA = np.unique(nodes["tag"], return_inverse=True)
     B, idxB = np.unique(tria["a"], return_inverse=True)
-    IDX = np.in1d(A, B)
+    IDX = np.isin(A, B)
     tria["a"] = idxA[IDX][idxB]
     B, idxB = np.unique(tria["b"], return_inverse=True)
-    IDX = np.in1d(A, B)
+    IDX = np.isin(A, B)
     tria["b"] = idxA[IDX][idxB]
     B, idxB = np.unique(tria["c"], return_inverse=True)
-    IDX = np.in1d(A, B)
+    IDX = np.isin(A, B)
     tria["c"] = idxA[IDX][idxB]
 
     # Drop invalid edges
@@ -258,10 +258,10 @@ def read_msh(filename, **kwargs):
     ### Re-index edges
     A, idxA = np.unique(nodes["tag"], return_inverse=True)
     B, idxB = np.unique(edges["e1"], return_inverse=True)
-    IDX = np.in1d(A, B)
+    IDX = np.isin(A, B)
     edges["e1"] = idxA[IDX][idxB]
     B, idxB = np.unique(edges["e2"], return_inverse=True)
-    IDX = np.in1d(A, B)
+    IDX = np.isin(A, B)
     edges["e2"] = idxA[IDX][idxB]
     # clean up
     nodes = nodes.drop("tag", axis=1)
@@ -479,20 +479,21 @@ def get(contours, **kwargs):
 
     hfun_scal = kwargs.get("hfun_scal", "ABSOLUTE")
     hfun_min = kwargs.get("hfun_min", 0.0)
-    hfun_max = kwargs.get("hfun_max", 1000.0)
+    hfun_max = kwargs.get("hfun_max", 100.0)
+    mesh_eps1 = kwargs.get("mesh_eps1", 0.3)
 
     with open(fjig, "w") as f:
         f.write("GEOM_FILE ={}\n".format(tag + "-geo.msh"))
         f.write("MESH_FILE ={}\n".format(tag + ".msh"))
         if bgmesh:
             f.write("HFUN_FILE ={}\n".format(tag + "-hfun.msh"))
-        f.write("HFUN_SCAL = {}\n".format(hfun_scal))
-        f.write("HFUN_HMAX = {}\n".format(hfun_max))
-        f.write("HFUN_HMIN = {}\n".format(hfun_min))
+        f.write(f"HFUN_SCAL = {hfun_scal}\n")
+        f.write(f"HFUN_HMAX = {hfun_max}\n")
+        f.write(f"HFUN_HMIN = {hfun_min}\n")
         f.write("MESH_DIMS = 2\n")
         f.write("MESH_TOP1 = TRUE\n")
         #        f.write('MESH_TOP2 = TRUE\n')
-        f.write("MESH_EPS1 = 1.0\n")
+        f.write(f"MESH_EPS1 = {mesh_eps1}\n")
         f.write("MESH_RAD2 = 1\n")
         f.write("GEOM_FEAT = TRUE\n")
         f.write("VERBOSITY = 2")
diff --git a/pyposeidon/tools.py b/pyposeidon/tools.py
index 4435f12..ec91e77 100644
--- a/pyposeidon/tools.py
+++ b/pyposeidon/tools.py
@@ -156,7 +156,6 @@ def setup_logging(
         logger.addHandler(file_handler)
 
 
-# TODO Handle master/develop version
 def parse_schism_version(version_output: str) -> str:
     if "schism develop" in version_output:
         version = "develop"
@@ -530,7 +529,7 @@ def dequantize(
     missing_value: int,
     dtype: npt.DTypeLike,
 ) -> npt.NDArray[np.float_]:
-    array = bn.replace(array.astype(np.float_), missing_value, np.nan)
+    array = bn.replace(array.astype(np.float64), missing_value, np.nan)
     dequantized = (array * scale_factor) + add_offset
     return dequantized
 
diff --git a/pyposeidon/utils/coastfix.py b/pyposeidon/utils/coastfix.py
index 785fca6..1c081cf 100644
--- a/pyposeidon/utils/coastfix.py
+++ b/pyposeidon/utils/coastfix.py
@@ -25,9 +25,9 @@ def simplify(geo):
 
     if (geo.geom_type == "Polygon").all():
         try:
-            geo_ = list(geo.buffer(0).unary_union.geoms)
+            geo_ = list(geo.buffer(0).union_all().geoms)
         except TypeError:
-            geo_ = [geo.buffer(0).unary_union]
+            geo_ = [geo.buffer(0).union_all()]
 
         geo = gp.GeoDataFrame(geometry=geo_)
         geo = gp.GeoDataFrame(geometry=geo.buffer(0))
@@ -40,7 +40,7 @@ def simplify(geo):
             for idx, geom in dg.iterrows():
                 pl = shapely.polygonize_full(dg.loc[[idx]].boundary.explode(index_parts=False).values)
                 df = gp.GeoDataFrame(geometry=[pl[0]]).explode(index_parts=False)
-                dff = gp.GeoDataFrame(geometry=[df.unary_union])
+                dff = gp.GeoDataFrame(geometry=[df.union_all()])
                 geo.loc[idx, "geometry"] = dff.geometry.buffer(0).values[0]
 
             if geo.is_valid.all() and (geo.boundary.geom_type == "LineString").all():
diff --git a/tests/conftest.py b/tests/conftest.py
index 88c0847..671ee8c 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -79,7 +79,8 @@ def save_as(url: str, path: pathlib.Path) -> None:
 
 
 TEST_DATA = {
-    DATA_DIR / "ocean.zip": "https://www.dropbox.com/sh/cgyk6ychwimjlzr/AADlueXH1I1F01mlgerMUSb6a?dl=1",
+    DATA_DIR
+    / "ocean.parquet": "https://www.dropbox.com/scl/fi/msiru06qdh7cq99lsxicx/ocean.parquet?rlkey=ku7afiues47y91a15hdvtcybu&st=nayi0hc2&dl=1",
     DATA_DIR / "bl.zip": "https://www.dropbox.com/sh/9tfdl67sll1ax8c/AACntQeIavCzRfTZZ9Tp8uFda?dl=1",
     DATA_DIR / "dem.nc": "https://www.dropbox.com/s/l16crheqc9d89gy/dem.nc?dl=1",
     DATA_DIR / "dem.tif": "https://www.dropbox.com/s/dgdlnr2p5q66np7/dem.tif?dl=1",
diff --git a/tests/test_boundary.py b/tests/test_boundary.py
index 472d276..7d19d75 100644
--- a/tests/test_boundary.py
+++ b/tests/test_boundary.py
@@ -6,9 +6,9 @@
 
 noaa = DATA_DIR / "bl.zip"
 
-COAST_FILE = (DATA_DIR / "ocean.zip").as_posix()
+COAST_FILE = (DATA_DIR / "ocean.parquet").as_posix()
 
-land = gp.read_file(COAST_FILE).drop("FID", axis=1)
+land = gp.read_file(COAST_FILE)
 coast = gp.GeoDataFrame(geometry=land.boundary)
 
 INPUTS = pytest.mark.parametrize("input", [land, coast])
diff --git a/tests/test_dem_fix.py b/tests/test_dem_fix.py
index 5d401cf..30bdbb5 100644
--- a/tests/test_dem_fix.py
+++ b/tests/test_dem_fix.py
@@ -6,7 +6,7 @@
 
 from . import DATA_DIR
 
-COAST_FILE = (DATA_DIR / "ocean.zip").as_posix()
+COAST_FILE = (DATA_DIR / "ocean.parquet").as_posix()
 
 DEM_SOURCES = pytest.mark.parametrize(
     "dem_source",
@@ -29,7 +29,7 @@
 
 @pytest.fixture(scope="session")
 def coasts():
-    coast = gp.read_file(COAST_FILE).drop("FID", axis=1)
+    coast = gp.read_file(COAST_FILE)
     return coast
 
 
diff --git a/tests/test_dem_tile.py b/tests/test_dem_tile.py
index 5047051..dc0ca13 100644
--- a/tests/test_dem_tile.py
+++ b/tests/test_dem_tile.py
@@ -6,7 +6,7 @@
 
 from . import DATA_DIR
 
-COAST_FILE = (DATA_DIR / "ocean.zip").as_posix()
+COAST_FILE = (DATA_DIR / "ocean.parquet").as_posix()
 
 DEM_SOURCES = pytest.mark.parametrize(
     "dem_source",
@@ -18,7 +18,7 @@
 
 @pytest.fixture(scope="session")
 def coasts():
-    #    coast = gp.read_file(COAST_FILE).drop("FID", axis=1)
+    #    coast = gp.read_file(COAST_FILE)
     cr = "h"
     coast = cf.NaturalEarthFeature(
         category="physical", name="land", scale="{}m".format({"l": 110, "i": 50, "h": 10}[cr])
diff --git a/tests/test_mesh.py b/tests/test_mesh.py
index 9688e68..e97d8ac 100644
--- a/tests/test_mesh.py
+++ b/tests/test_mesh.py
@@ -5,7 +5,7 @@
 
 from . import DATA_DIR
 
-COAST_FILE = (DATA_DIR / "ocean.zip").as_posix()
+COAST_FILE = (DATA_DIR / "ocean.parquet").as_posix()
 
 WINDOWS = pytest.mark.parametrize(
     "window",
diff --git a/tests/test_mesh_global.py b/tests/test_mesh_global.py
index d79dd64..eff2b29 100644
--- a/tests/test_mesh_global.py
+++ b/tests/test_mesh_global.py
@@ -9,13 +9,13 @@
 
 DEM_FILE = (DATA_DIR / "dem.nc").as_posix()
 
-COAST_FILE = (DATA_DIR / "ocean.zip").as_posix()
+COAST_FILE = (DATA_DIR / "ocean.parquet").as_posix()
 
 
 @pytest.mark.parametrize("ggor", ["jigsaw", "gmsh"])
 @pytest.mark.parametrize("bgmesh", [None, DEM_FILE])
 @pytest.mark.parametrize("bindings", [True, False])
-@pytest.mark.parametrize("cbuffer", [None, 0.001])
+@pytest.mark.parametrize("cbuffer", [None, 0.01])
 def test_io(pytestconfig, tmpdir, ggor, bgmesh, bindings, cbuffer):
     # Skip the test unless --runslow has been passed
     if bgmesh is not None:
@@ -61,7 +61,7 @@ def test_io(pytestconfig, tmpdir, ggor, bgmesh, bindings, cbuffer):
 @pytest.mark.parametrize("ggor", ["jigsaw", "gmsh"])
 @pytest.mark.parametrize("bgmesh", [None, DEM_FILE])
 @pytest.mark.parametrize("bindings", [True, False])
-@pytest.mark.parametrize("cbuffer", [None, 0.001])
+@pytest.mark.parametrize("cbuffer", [None, 0.01])
 def test_validate(pytestconfig, tmpdir, ggor, cbuffer, bgmesh, bindings):
     if bgmesh is not None:
         if not pytestconfig.getoption("--runslow"):
@@ -69,6 +69,9 @@ def test_validate(pytestconfig, tmpdir, ggor, cbuffer, bgmesh, bindings):
     if ggor == "jigsaw":
         pytest.xfail("Fixing these is a WIP")
 
+    #    if ggor == "jigsaw":
+    #        pytest.xfail("Fixing these is a WIP")
+
     mesh = pmesh.set(
         type="tri2d",
         geometry="global",
diff --git a/tests/test_schism_total.py b/tests/test_schism_total.py
index d1af198..5a05dbd 100644
--- a/tests/test_schism_total.py
+++ b/tests/test_schism_total.py
@@ -8,7 +8,7 @@
 
 DEM_FILE = DATA_DIR / "dem.nc"
 
-COAST_FILE = (DATA_DIR / "ocean.zip").as_posix()
+COAST_FILE = (DATA_DIR / "ocean.parquet").as_posix()
 
 WINDOWS = pytest.mark.parametrize(
     "window",