Major revision (#63)

* IMPRV: adding fort.14 writer

* IMPRV: adding t3s writer

* IMPRV: adding in t3s writer

* Forgot to commit

* Bugfix: fixing reading of DEM

* IMPRV: clipping DEMs and new method total_bounds for Region

* Major revision of DEM class, slope edgefx, new distance from point edgefx, impr. vis. of edgefx, new helper functions in several places

* Get polygon from various input types

* IMPRV: adding in resolution from line function

* BUGFIX: handle multipolygons. IMPRV: begin work on labeling boundaries

* IMPRV: adding automatic detection of ocean boundaries

* IMPRV: more refinement for the ocean boundary

* IMPRV: general improvements

* BUGFIX: wl mesh size function in geographic coords

* Fix code style issues with Black

* FIX: formatting

* FIX: addressing flake8

* Fix code style issues with Black

---------

Co-authored-by: oceanmesh <[email protected]>
Co-authored-by: Lint Action <[email protected]>

oceanmesh/__init__.py

@@ -8,31 +8,54 @@
     ), "The environment variable CGAL_BIN must be set."
     os.add_dll_directory(os.environ["CGAL_BIN"])
 
-from oceanmesh.clean import (delete_boundary_faces, delete_exterior_faces,
-                             delete_faces_connected_to_one_face,
-                             delete_interior_faces, laplacian2,
-                             make_mesh_boundaries_traversable)
-from oceanmesh.edgefx import (bathymetric_gradient_sizing_function,
-                              distance_sizing_function, enforce_mesh_gradation,
-                              enforce_mesh_size_bounds_elevation,
-                              feature_sizing_function,
-                              multiscale_sizing_function,
-                              wavelength_sizing_function)
+from oceanmesh.boundary import identify_ocean_boundary_sections
+from oceanmesh.clean import (
+    delete_boundary_faces,
+    delete_exterior_faces,
+    delete_faces_connected_to_one_face,
+    delete_interior_faces,
+    laplacian2,
+    make_mesh_boundaries_traversable,
+    mesh_clean,
+)
+from oceanmesh.edgefx import (
+    bathymetric_gradient_sizing_function,
+    distance_sizing_from_line_function,
+    distance_sizing_from_point_function,
+    distance_sizing_function,
+    enforce_mesh_gradation,
+    enforce_mesh_size_bounds_elevation,
+    feature_sizing_function,
+    multiscale_sizing_function,
+    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.geodata import (
+    DEM,
+    Shoreline,
+    create_circle_coords,
+    get_polygon_coordinates,
+)
 from oceanmesh.grid import Grid, compute_minimum
 from oceanmesh.region import Region, warp_coordinates
 from oceanmesh.signed_distance_function import (
-    Difference, Domain, Intersection, Union, create_bbox, create_circle,
-    multiscale_signed_distance_function, signed_distance_function)
+    Difference,
+    Domain,
+    Intersection,
+    Union,
+    create_bbox,
+    create_circle,
+    multiscale_signed_distance_function,
+    signed_distance_function,
+)
 
 from .fix_mesh import fix_mesh, simp_vol
-
 from .mesh_generator import (
     generate_mesh,
     generate_multiscale_mesh,
-    plot_mesh,
+    plot_mesh_bathy,
+    plot_mesh_connectivity,
     write_to_fort14,
     write_to_t3s,
 )
@@ -41,16 +64,21 @@
     "create_bbox",
     "Region",
     "compute_minimum",
+    "create_circle_coords",
     "bathymetric_gradient_sizing_function",
     "multiscale_sizing_function",
     "delete_boundary_faces",
     "delete_faces_connected_to_one_face",
-    "plot_mesh",
+    "distance_sizing_from_point_function",
+    "distance_sizing_from_line_function",
+    "plot_mesh_connectivity",
+    "plot_mesh_bathy",
     "make_mesh_boundaries_traversable",
     "enforce_mesh_size_bounds_elevation",
     "laplacian2",
     "delete_interior_faces",
     "delete_exterior_faces",
+    "mesh_clean",
     "SizeFunction",
     "Grid",
     "DEM",
@@ -60,8 +88,10 @@
     "Union",
     "Difference",
     "Intersection",
+    "identify_ocean_boundary_sections",
     "Shoreline",
     "generate_multiscale_mesh",
+    "get_polygon_coordinates",
     "distance_sizing_function",
     "feature_sizing_function",
     "enforce_mesh_gradation",

oceanmesh/boundary.py

@@ -0,0 +1,106 @@
+import numpy as np
+
+from .edges import get_winded_boundary_edges
+import matplotlib.pyplot as plt
+
+__all__ = ["identify_ocean_boundary_sections"]
+
+
+def identify_ocean_boundary_sections(
+    points,
+    cells,
+    topobathymetry,
+    depth_threshold=-50.0,
+    min_nodes_threshold=10,
+    plot=False,
+):
+    """Identify the contiguous sections on the ocean boundary based on depth
+    that could be forced in a numerical model as ocean-type boundaries (e.g., elevation-specified)
+
+    Parameters
+    ----------
+    points: numpy.ndarray
+        Array of points (x,y)
+    cells : numpy.ndarray
+        Array of cells
+    topobathymetry : numpy.ndarray
+        Array of topobathymetry values (depth below datum negative)
+    depth_threshold : float, optional
+        Depth threshold to identify ocean boundary nodes, by default -50 m below the datum
+    min_nodes_threshold : int, optional
+        Minimum number of nodes to be considered a boundary section, by default 10
+    plot : bool, optional
+        Plot the mesh and the identified boundary sections, by default False
+
+    Returns
+    --------
+    boundary_sections : list
+        List of tuples of the nodes that define the ocean boundary sections
+        Note these map back into the points array.
+
+    """
+    # Identify the nodes on the boundary of the mesh
+    boundary_edges = get_winded_boundary_edges(cells)
+    boundary_edges = boundary_edges.flatten()
+    unique_indexes = np.unique(boundary_edges, return_index=True)[1]
+    boundary_nodes_unmasked = [
+        boundary_edges[unique_index] for unique_index in sorted(unique_indexes)
+    ]
+    # Define a boolean array of valid nodes
+    bathymetry_on_boundary = topobathymetry[boundary_nodes_unmasked]
+    # Append a NaN value to the array to align with the original
+    bathymetry_on_boundary = np.append(bathymetry_on_boundary, np.nan)
+    stops = np.nonzero(bathymetry_on_boundary <= depth_threshold)[0]
+
+    # Plot the mesh
+    if plot:
+        fig, ax = plt.subplots()
+        ax.triplot(points[:, 0], points[:, 1], cells, color="k", lw=0.1)
+
+    first = True
+    boundary_sections = []
+    start_node = None
+    end_node = None
+    for idx, (s1, s2) in enumerate(zip(stops[:-1], stops[1:])):
+        if s2 - s1 < min_nodes_threshold:
+            if first:
+                start_node = s1
+                first = False
+            # We've reached the end of the list
+            elif idx == len(stops) - 2:
+                # Append the start and end nodes to the boundary sections list
+                end_node = s2
+                boundary_sections.append([start_node, end_node])
+            # Its not the end and we haven't found a section yet
+            else:
+                end_node = s2
+        elif s2 - s1 >= min_nodes_threshold and not first:
+            # Append the start and end nodes to the boundary sections list
+            boundary_sections.append([start_node, end_node])
+            # Reset the start node, the last node didn't satisfy the threshold
+            # and it appears we have a new section
+            start_node = s1
+            first = True
+        # We've reached the end of the list
+        elif idx == len(stops) - 2:
+            # Save the end node
+            end_node = s2
+            # Append the start and end nodes to the boundary sections list and finish
+            boundary_sections.append([start_node, end_node])
+    if plot:
+        for s1, s2 in boundary_sections:
+            ax.scatter(
+                points[boundary_nodes_unmasked[s1:s2], 0],
+                points[boundary_nodes_unmasked[s1:s2], 1],
+                5,
+                c="r",
+            )
+            ax.set_title("Identified ocean boundary sections")
+        plt.show()
+
+    # Map back to the original node indices associated with the points array
+    boundary_sections = [
+        (boundary_nodes_unmasked[s1], boundary_nodes_unmasked[s2])
+        for s1, s2 in boundary_sections
+    ]
+    return boundary_sections

oceanmesh/clean.py

@@ -16,9 +16,61 @@
     "delete_faces_connected_to_one_face",
     "delete_boundary_faces",
     "laplacian2",
+    "mesh_clean",
 ]
 
 
+def mesh_clean(
+    points,
+    cells,
+    min_element_qual=0.01,
+    min_percent_disconnected_area=0.05,
+    max_iter=20,
+    tol=0.01,
+    pfix=None,
+):
+    """Clean a mesh by removing bad quality elements and boundary faces.
+
+    Parameters
+    ----------
+    points : array-like
+        Mesh vertices.
+    cells : array-like
+        Mesh connectivity table.
+    min_element_qual : float, optional
+        Enforce Mmnimum quality of elements through deletion. The default is 0.01.
+    min_percent_disconnected_area : float, optional
+        Delete disconnected areas smaller than this threshold. The default is 0.05 x
+        the total area of the mesh.
+    max_iter : int, optional
+        Maximum number of iterations for the Laplacian smoothing. The default is 20.
+    tol : float, optional
+        Tolerance for the Laplacian smoothing. The default is 0.01. Laplacian terminates
+        after max_iter or when the maximum change in any vertex is less than tol.
+    pfix : array-like, optional
+        Indices of points to fix during the Laplacian smoothing. The default is None.
+
+    Returns
+    -------
+    points : array-like
+        Mesh vertices cleaned.
+    cells : array-like
+        Mesh connectivity table cleaned.
+
+    """
+    points, cells = make_mesh_boundaries_traversable(
+        points, cells, min_disconnected_area=min_percent_disconnected_area
+    )
+    points, cells = delete_boundary_faces(points, cells, min_qual=min_element_qual)
+    points, cells = delete_faces_connected_to_one_face(points, cells)
+    points, cells = laplacian2(points, cells, max_iter=max_iter, tol=tol, pfix=pfix)
+    points, cells = make_mesh_boundaries_traversable(
+        points, cells, min_disconnected_area=min_percent_disconnected_area
+    )
+    points, cells, _ = fix_mesh(points, cells, delete_unused=True)
+    return points, cells
+
+
 def _arg_sortrows(arr):
     """Before a multi column sort like MATLAB's sortrows"""
     i = arr[:, 1].argsort()  # First sort doesn't need to be stable.
@@ -293,7 +345,7 @@ def _depth_first_search(faces):
     return nflag
 
 
-def delete_faces_connected_to_one_face(vertices, faces, max_iter=5):
+def delete_faces_connected_to_one_face(vertices, faces, max_iter=np.inf):
     """Iteratively deletes faces connected to one face.
 
     Parameters
@@ -304,7 +356,8 @@ def delete_faces_connected_to_one_face(vertices, faces, max_iter=5):
         The "uncleaned" mesh connectivity.
     max_iter: float, optional
         The number of iterations to repeatedly delete faces connected to one face
-
+        If the mesh is well-formed, this will converge in a finite
+        number of iterations. Default is np.inf.
 
     Returns
     -------