Merge pull request #69 from noaa-ocs-modeling/bugfix/hfuntype

Hfun Collector Valid Input Types

ocsmesh/hfun/collector.py

@@ -50,6 +50,7 @@
 from ocsmesh.hfun.raster import HfunRaster
 from ocsmesh.hfun.mesh import HfunMesh
 from ocsmesh.mesh.mesh import Mesh, EuclideanMesh2D
+from ocsmesh.mesh.base import BaseMesh
 from ocsmesh.raster import Raster, get_iter_windows
 from ocsmesh.features.contour import Contour
 from ocsmesh.features.patch import Patch
@@ -722,7 +723,8 @@ def __init__(
             elif isinstance(in_item, EuclideanMesh2D):
                 hfun = HfunMesh(in_item)
 
-            elif isinstance(in_item, str):
+            elif isinstance(in_item, (str, Path)):
+                in_item = str(in_item)
                 if in_item.endswith('.tif'):
                     raster = Raster(in_item)
                     if self._base_shape:
@@ -1359,7 +1361,7 @@ def _type_chk(input_list: List[Any]) -> None:
             If the any of the inputs are not supported
         """
 
-        valid_types = (str, Raster, Mesh, HfunRaster, HfunMesh)
+        valid_types = (str, Path, Raster, BaseMesh, HfunRaster, HfunMesh)
         if not all(isinstance(item, valid_types) for item in input_list):
             raise TypeError(
                 f'Input list items must be of type'

tests/api/__init__.py

@@ -1 +1,7 @@
-
+import os
+os.system("""
+    ls /tmp/test_dem.tif > /dev/null
+    if [ $? -eq 0 ]; then exit; fi
+    wget https://coast.noaa.gov/htdata/raster2/elevation/NCEI_ninth_Topobathy_2014_8483/northeast_sandy/ncei19_n40x75_w073x75_2015v1.tif -O /tmp/fullsize_dem.tif
+    gdalwarp -tr 0.0005 0.0005 -r average -overwrite /tmp/fullsize_dem.tif /tmp/test_dem.tif
+""")

tests/api/common.py

@@ -0,0 +1,123 @@
+import rasterio as rio
+import numpy as np
+from pyproj import CRS
+from jigsawpy import jigsaw_msh_t
+
+
+def create_rectangle_mesh(nx, ny, holes, x_extent=None, y_extent=None):
+    """
+    Note:
+        x = x-index
+        y = y-index
+
+        node-index(node-id)
+
+          25(26)             29(30)
+      5     *---*---*---*---*
+            | \ | \ | \ | \ |
+      4     *---*---*---*---*
+            | \ | \ | \ | \ |
+      3     *---*---*---*---*
+            | \ |   | \ | \ |
+      2     *---*---*---*---*
+            | \ | \ | \ | \ |
+      1     *---*---*---*---*
+            | \ | \ | \ | \ |
+      0     *---*---*---*---*
+          0(1)               4(5)
+
+            0   1   2   3   4
+    """
+
+    if x_extent is None:
+        x_range = range(nx)
+    else:
+        x_range = np.linspace(x_extent[0], x_extent[1], nx)
+
+    if y_extent is None:
+        y_range = range(ny)
+    else:
+        y_range = np.linspace(y_extent[0], y_extent[1], nx)
+
+    X, Y = np.meshgrid(x_range, y_range)
+    verts = np.array(list(zip(X.ravel(), Y.ravel())))
+    cells = []
+    for j in range(ny - 1):
+        for i in range(nx - 1):
+            if (i + 1) + ((nx-1) * j) in holes:
+                continue
+            cells.append([j * nx + i, j * nx + (i + 1), (j + 1) * nx + i])
+            cells.append([j * nx + (i + 1), (j + 1) * nx + (i + 1), (j + 1) * nx + i])
+    # NOTE: Everywhere is above 0 (auto: land) unless modified later
+    vals = np.ones((len(verts), 1)) * 10
+
+    # TODO: Replace with "make_mesh" util function
+    mesh_msht = jigsaw_msh_t()
+    mesh_msht.ndims = +2
+    mesh_msht.mshID = 'euclidean-mesh'
+    mesh_msht.tria3 = np.array(
+        [(c, 0) for c in cells], dtype=jigsaw_msh_t.TRIA3_t
+    )
+    mesh_msht.vert2 = np.array(
+        [(v, 0) for v in verts], dtype=jigsaw_msh_t.VERT2_t
+    )
+    mesh_msht.value = np.array(
+        vals, dtype=jigsaw_msh_t.REALS_t
+    )
+    return mesh_msht
+
+
+# TODO: Move these helper functions to `utils`
+def raster_from_numpy(
+    filename,
+    data,
+    mgrid,
+    crs=CRS.from_epsg(4326)
+):
+    x = mgrid[0][:, 0]
+    y = mgrid[1][0, :]
+    res_x = (x[-1] - x[0]) / data.shape[1]
+    res_y = (y[-1] - y[0]) / data.shape[0]
+    transform = rio.transform.Affine.translation(
+        x[0], y[0]
+    ) * rio.transform.Affine.scale(res_x, res_y)
+    if not isinstance(crs, CRS):
+        crs = CRS.from_user_input(crs)
+
+    with rio.open(
+        filename,
+        'w',
+        driver='GTiff',
+        height=data.shape[0],
+        width=data.shape[1],
+        count=1,
+        dtype=data.dtype,
+        crs=crs,
+        transform=transform,
+    ) as dst:
+        dst.write(data, 1)
+
+
+def msht_from_numpy(
+    coordinates,
+    triangles,
+    crs=CRS.from_epsg(4326)
+):
+    if not isinstance(crs, CRS):
+        crs = CRS.from_user_input(crs)
+    mesh = jigsaw_msh_t()
+    mesh.mshID = 'euclidean-mesh'
+    mesh.ndims = +2
+    mesh.crs = crs
+    mesh.vert2 = np.array(
+        [(coord, 0) for coord in coordinates],
+        dtype=jigsaw_msh_t.VERT2_t
+        )
+    mesh.tria3 = np.array(
+        [(index, 0) for index in triangles],
+        dtype=jigsaw_msh_t.TRIA3_t
+        )
+
+    return mesh
+
+

tests/api/hfun.py

@@ -4,6 +4,7 @@
 from pathlib import Path
 import shutil
 import tempfile
+import warnings
 
 from jigsawpy import jigsaw_msh_t
 import geopandas as gpd
@@ -15,59 +16,139 @@
 
 import ocsmesh
 
+from tests.api.common import raster_from_numpy, msht_from_numpy, create_rectangle_mesh
 
-# TODO: Move these helper functions to `utils`
-def _helper_raster_from_numpy(
-    filename,
-    data,
-    mgrid,
-    crs=CRS.from_epsg(4326)
-):
-    x = mgrid[0][:, 0]
-    y = mgrid[1][0, :]
-    res_x = (x[-1] - x[0]) / data.shape[1]
-    res_y = (y[-1] - y[0]) / data.shape[0]
-    transform = rio.transform.Affine.translation(
-        x[0], y[0]
-    ) * rio.transform.Affine.scale(res_x, res_y)
-    if not isinstance(crs, CRS):
-        crs = CRS.from_user_input(crs)
-
-    with rio.open(
-        filename,
-        'w',
-        driver='GTiff',
-        height=data.shape[0],
-        width=data.shape[1],
-        count=1,
-        dtype=data.dtype,
-        crs=crs,
-        transform=transform,
-    ) as dst:
-        dst.write(data, 1)
-
-
-def _helper_msh_t_from_numpy(
-    coordinates,
-    triangles,
-    crs=CRS.from_epsg(4326)
-):
-    if not isinstance(crs, CRS):
-        crs = CRS.from_user_input(crs)
-    mesh = jigsaw_msh_t()
-    mesh.mshID = 'euclidean-mesh'
-    mesh.ndims = +2
-    mesh.crs = crs
-    mesh.vert2 = np.array(
-        [(coord, 0) for coord in coordinates],
-        dtype=jigsaw_msh_t.VERT2_t
-        )
-    mesh.tria3 = np.array(
-        [(index, 0) for index in triangles],
-        dtype=jigsaw_msh_t.TRIA3_t
-        )
-
-    return mesh
+
+class SizeFunctionType(unittest.TestCase):
+    def setUp(self):
+        self.tdir = Path(tempfile.mkdtemp())
+
+        self.rast = self.tdir / 'rast_1.tif'
+        self.mesh = self.tdir / 'mesh_1.gr3'
+
+        rast_xy = np.mgrid[-1:0.1:0.1, -0.7:0.1:0.1]
+        rast_xy = np.mgrid[0:1.1:0.1, -0.7:0.1:0.1]
+        rast_z = np.ones_like(rast_xy[0])
+
+        raster_from_numpy(
+            self.rast, rast_z, rast_xy, 4326
+        )
+
+        msh_t = create_rectangle_mesh(
+            nx=17, ny=7, holes=[40, 41], x_extent=(-1, 1), y_extent=(0, 1))
+
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                'ignore', category=UserWarning, message='Input mesh has no CRS information'
+            )
+            mesh = ocsmesh.Mesh(msh_t)
+            mesh.write(str(self.mesh), format='grd', overwrite=False)
+
+
+    def test_create_raster_hfun(self):
+        hfun = ocsmesh.Hfun(
+            ocsmesh.Raster(self.rast),
+            hmin=500,
+            hmax=10000
+        )
+        self.assertTrue(isinstance(hfun, ocsmesh.hfun.raster.HfunRaster))
+
+    def test_mesh_raster_hfun(self):
+        hfun = ocsmesh.Hfun(
+            ocsmesh.Mesh.open(self.mesh, crs=4326),
+        )
+        self.assertTrue(isinstance(hfun, ocsmesh.hfun.mesh.HfunMesh))
+
+    def test_collector_raster_hfun(self):
+        hfun = ocsmesh.Hfun(
+            [self.rast],
+            hmin=500,
+            hmax=10000
+        )
+        self.assertTrue(isinstance(hfun, ocsmesh.hfun.collector.HfunCollector))
+
+
+class SizeFunctionCollector(unittest.TestCase):
+
+    def setUp(self):
+        self.tdir = Path(tempfile.mkdtemp())
+
+        self.rast1 = self.tdir / 'rast_1.tif'
+        self.rast2 = self.tdir / 'rast_2.tif'
+        self.mesh1 = self.tdir / 'mesh_1.gr3'
+
+        rast_xy_1 = np.mgrid[-1:0.1:0.1, -0.7:0.1:0.1]
+        rast_xy_2 = np.mgrid[0:1.1:0.1, -0.7:0.1:0.1]
+        rast_z_1 = np.ones_like(rast_xy_1[0])
+
+        raster_from_numpy(
+            self.rast1, rast_z_1, rast_xy_1, 4326
+        )
+        raster_from_numpy(
+            self.rast2, rast_z_1, rast_xy_2, 4326
+        )
+
+        msh_t = create_rectangle_mesh(
+            nx=17, ny=7, holes=[40, 41], x_extent=(-1, 1), y_extent=(0, 1))
+
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                'ignore', category=UserWarning, message='Input mesh has no CRS information'
+            )
+            mesh = ocsmesh.Mesh(msh_t)
+            mesh.write(str(self.mesh1), format='grd', overwrite=False)
+
+
+    def test_multi_path_input(self):
+        hfun_coll = ocsmesh.Hfun(
+            [self.rast1, self.rast2, self.mesh1],
+            hmin=500,
+            hmax=10000
+        )
+        hfun_msht = hfun_coll.msh_t()
+        self.assertTrue(isinstance(hfun_msht, jigsaw_msh_t))
+
+    def test_multi_str_input(self):
+        hfun_coll = ocsmesh.Hfun(
+            [str(i) for i in [self.rast1, self.rast2, self.mesh1]],
+            hmin=500,
+            hmax=10000
+        )
+        hfun_msht = hfun_coll.msh_t()
+        self.assertTrue(isinstance(hfun_msht, jigsaw_msh_t))
+
+    def test_multi_raster_input(self):
+
+        rast1 = ocsmesh.Raster(self.rast1)
+        rast2 = ocsmesh.Raster(self.rast2)
+        hfun_coll = ocsmesh.Hfun(
+            [rast1, rast2],
+            hmin=500,
+            hmax=10000
+        )
+        hfun_msht = hfun_coll.msh_t()
+        self.assertTrue(isinstance(hfun_msht, jigsaw_msh_t))
+
+    def test_multi_mix_input(self):
+        rast1 = ocsmesh.Raster(self.rast1)
+        mesh1 = ocsmesh.Mesh.open(self.mesh1, crs=4326)
+        hfun_coll = ocsmesh.Hfun(
+            [rast1, self.rast2, mesh1],
+            hmin=500,
+            hmax=10000
+        )
+        hfun_msht = hfun_coll.msh_t()
+        self.assertTrue(isinstance(hfun_msht, jigsaw_msh_t))
+
+    def test_mesh_input(self):
+        mesh1 = ocsmesh.Mesh.open(self.mesh1, crs=4326)
+        hfun_coll = ocsmesh.Hfun(
+            [mesh1],
+            hmin=500,
+            hmax=10000
+        )
+        hfun_msht = hfun_coll.msh_t()
+        self.assertTrue(isinstance(hfun_msht, jigsaw_msh_t))
 
 
 
@@ -482,10 +563,10 @@ def setUp(self):
         rast_xy_2 = np.mgrid[0:1.1:0.1, -0.7:0.1:0.1]
         rast_z_1 = np.ones_like(rast_xy_1[0])
 
-        _helper_raster_from_numpy(
+        raster_from_numpy(
             self.rast1, rast_z_1, rast_xy_1, 4326
         )
-        _helper_raster_from_numpy(
+        raster_from_numpy(
             self.rast2, rast_z_1, rast_xy_2, 4326
         )
 
@@ -509,7 +590,7 @@ def setUp(self):
             [0, 2, 6],
             [5, 4, 7],
         ])
-        msh_t = _helper_msh_t_from_numpy(crd, tria, 4326)
+        msh_t = msht_from_numpy(crd, tria, 4326)
         mesh = ocsmesh.Mesh(msh_t)
         mesh.write(str(self.mesh1), format='grd', overwrite=False)