Repository: DeltaRCM/DeltaMetrics

Issue:
Masks currently locked to sea level == 0
Are masks currently locked to sea level == 0? I think they might be. This should be relaxed, at least with an option to specify during mask creation. 

related to #10. Should consider how we will include boundary conditions in an overhaul of the netcdf file output from pyDeltaRCM too.

Discussion:

Patch:
diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
index 8c71ab35..f18e5f56 100644
--- a/.github/workflows/build.yml
+++ b/.github/workflows/build.yml
@@ -48,10 +48,10 @@ jobs:
       uses: actions/checkout@v2
       with:
         persist-credentials: false
-    - name: Set up Python 3.x
+    - name: Set up Python 3.8
       uses: actions/setup-python@v2
       with:
-        python-version: '3.x'
+        python-version: '3.8'
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
diff --git a/deltametrics/mask.py b/deltametrics/mask.py
index 31e11614..2770a771 100644
--- a/deltametrics/mask.py
+++ b/deltametrics/mask.py
@@ -25,6 +25,16 @@ class BaseMask(object):
         else:
             self.data = data
 
+        # set data to be 3D even if it is not (assuming t-x-y)
+        if len(np.shape(self.data)) == 3:
+            pass
+        elif len(np.shape(self.data)) == 2:
+            self.data = np.reshape(self.data, [1,
+                                               np.shape(self.data)[0],
+                                               np.shape(self.data)[1]])
+        else:
+            raise ValueError('Input data shape was not 2-D nor 3-D')
+
         self._mask = np.zeros(self.data.shape)
 
     @property
@@ -48,15 +58,22 @@ class BaseMask(object):
         """
         return self._mask
 
-    def show(self, **kwargs):
+    def show(self, t=-1, **kwargs):
         """Show the mask.
 
+        Parameters
+        ----------
+        t : int, optional
+            Value of time slice or integer from first dimension in 3D (t-x-y)
+            convention to display the mask from. Default is -1 so the 'final'
+            mask in time is displayed if this argument is not supplied.
+
         Passes `**kwargs` to ``matplotlib.imshow``.
         """
         cmap = kwargs.pop('cmap', 'gray')
         fig, ax = plt.subplots()
         if hasattr(self, 'mask') and np.sum(self.mask) > 0:
-            ax.imshow(self.mask, cmap=cmap, **kwargs)
+            ax.imshow(self.mask[t, :, :], cmap=cmap, **kwargs)
             ax.set_title('A ' + self.mask_type + ' mask')
         else:
             raise AttributeError('No mask computed and/or mask is all 0s.')
@@ -107,49 +124,51 @@ class OAM(object):
         # write the parameters to self
         self.angle_threshold = angle_threshold
 
-        # use the first column to trim the land_width
-        # (or, can we assume access to the DeltaRCM variable `L0`?)
-        i = 0
-        delt = 10
-        while i < self.data.shape[0] and delt != 0:
-            delt = self.data[i, 0] - self.data[i+1, 0]
-            i += 1
-        trim_idx = i - 1  # assign the trimming index
-        data_trim = self.data[trim_idx:, :]
-        # use topo_threshold to identify oceanmap
-        self.oceanmap = (data_trim < self.topo_threshold) * 1.
-        # if all ocean, there is no shore to be found - exit function
-        if (self.oceanmap == 1).all():
-            self.oceanmap = np.zeros_like(self.data)
-            self.shore_image = np.zeros_like(self.data)
-            return
-        # apply seaangles function
-        _, seaangles = self.Seaangles_mod(self.numviews,
-                                          self.oceanmap)
-        # translate flat seaangles values to the shoreline image
-        shore_image = np.zeros_like(data_trim)
-        flat_inds = list(map(lambda x: np.ravel_multi_index(x,
+        # initialize arrays
+        self.oceanmap = np.zeros_like(self.data)
+        self.shore_image = np.zeros_like(self.data)
+
+        # loop through the time dimension
+        for tval in range(0, self.data.shape[0]):
+            # use the first column to trim the land_width
+            # (or, can we assume access to the DeltaRCM variable `L0`?)
+            i = 0
+            delt = 10
+            while i < self.data.shape[1] and delt != 0:
+                delt = self.data[tval, i, 0] - self.data[tval, i+1, 0]
+                i += 1
+            trim_idx = i - 1  # assign the trimming index
+            data_trim = self.data[tval, trim_idx:, :]
+            # use topo_threshold to identify oceanmap
+            omap = (data_trim < self.topo_threshold) * 1.
+            # if all ocean, there is no shore to be found - do next t-slice
+            if (omap == 1).all():
+                pass
+            else:
+                # apply seaangles function
+                _, seaangles = self.Seaangles_mod(self.numviews, omap)
+                # translate flat seaangles values to the shoreline image
+                shore_image = np.zeros_like(data_trim)
+                flat_inds = list(map(lambda x: np.ravel_multi_index(x,
                                                             shore_image.shape),
-                             seaangles[:2, :].T.astype(int)))
-        shore_image.flat[flat_inds] = seaangles[-1, :]
-        # grab contour from seaangles corresponding to angle threshold
-        cs = measure.find_contours(shore_image, np.float(self.angle_threshold))
-        C = cs[0]
-        # convert this extracted contour to the shoreline mask
-        shoremap = np.zeros_like(data_trim)
-        flat_inds = list(map(lambda x: np.ravel_multi_index(x, shoremap.shape),
-                             np.round(C).astype(int)))
-        shoremap.flat[flat_inds] = 1
-
-        # write shoreline map out to data.mask
-        self._mask[trim_idx:, :] = shoremap
-        # assign shore_image to the mask object with proper size
-        self.shore_image = np.zeros(self.data.shape)
-        self.shore_image[trim_idx:, :] = shore_image
-        # resize the oceanmap
-        oceanmap_temp = self.oceanmap.copy()
-        self.oceanmap = np.zeros(self.data.shape)
-        self.oceanmap[trim_idx:, :] = oceanmap_temp
+                                     seaangles[:2, :].T.astype(int)))
+                shore_image.flat[flat_inds] = seaangles[-1, :]
+                # grab contour from seaangles corresponding to angle threshold
+                cs = measure.find_contours(shore_image,
+                                           np.float(self.angle_threshold))
+                C = cs[0]
+                # convert this extracted contour to the shoreline mask
+                shoremap = np.zeros_like(data_trim)
+                flat_inds = list(map(lambda x: np.ravel_multi_index(x,
+                                                            shoremap.shape),
+                                     np.round(C).astype(int)))
+                shoremap.flat[flat_inds] = 1
+                # write shoreline map out to data.mask
+                self._mask[tval, trim_idx:, :] = shoremap
+                # assign shore_image to the mask object with proper size
+                self.shore_image[tval, trim_idx:, :] = shore_image
+                # properly assign the oceanmap to the self.oceanmap
+                self.oceanmap[tval, trim_idx:, :] = omap
 
     def Seaangles_mod(self, numviews, thresholdimg):
         """Extract the opening angle map from an image.
@@ -754,9 +773,14 @@ class EdgeMask(BaseMask, OAM):
             self.landmask = (self.shore_image < self.angle_threshold) * 1
             self.wetmask = self.oceanmap * self.landmask
         # compute the edge mask
-        self._mask = np.maximum(0,
-                                feature.canny(self.wetmask)*1 -
-                                feature.canny(self.landmask)*1)
+        for i in range(0, np.shape(self._mask)[0]):
+            self._mask[i, :, :] = np.maximum(0,
+                                             feature.canny(self.wetmask[i,
+                                                                        :,
+                                                                        :])*1 -
+                                             feature.canny(self.landmask[i,
+                                                                         :,
+                                                                         :])*1)
 
 
 class CenterlineMask(BaseMask):
@@ -865,7 +889,8 @@ class CenterlineMask(BaseMask):
         """
         # skimage.morphology.skeletonize() method
         if self.method == 'skeletonize':
-            self._mask = morphology.skeletonize(self.data)
+            for i in range(0, np.shape(self._mask)[0]):
+                self._mask[i, :, :] = morphology.skeletonize(self.data[i, :, :])
 
         if self.method == 'rivamap':
             # rivamap based method - first check for import error
diff --git a/docs/source/guides/userguide.rst b/docs/source/guides/userguide.rst
index c21621fc..e3d53f57 100644
--- a/docs/source/guides/userguide.rst
+++ b/docs/source/guides/userguide.rst
@@ -2,7 +2,7 @@
 User Guide
 **********
 
-This documentation provides a 
+This documentation provides a
 
 
 
@@ -34,7 +34,7 @@ Setting up your coding environment
 All of the documentation in this package assumes that you have imported the DeltaMetrics package as ``dm``:
 
 .. doctest::
-    
+
     >>> import deltametrics as dm
 
 Additionally, we frequently rely on the `numpy` package, and `matplotlib`. We will assume you have imported these packages by their common shorthand as well; if we import other packages, or other modules from `matplotlib`, these imports will be declared!
@@ -58,7 +58,7 @@ DeltaMetrics centers around the use of “Cubes” in DeltaMetrics language are
 
 Creating the ``rcm8cube`` connects to a dataset, but does not read any of the data into memory, allowing for efficient computation on large datasets. The type of the ``rcm8cube`` is ``DataCube``.
 
-Inspect which variables are available in the ``rcm8cube``. 
+Inspect which variables are available in the ``rcm8cube``.
 
 .. doctest::
 
@@ -97,7 +97,7 @@ Let’s examine the timeseries of bed elevations by taking slices out of the ``'
 
 .. plot:: guides/userguide_bed_timeseries.py
 
-.. note:: 
+.. note::
 
     The 0th dimension of the cube is the *time* dimension, and the 1st and 2nd dimensions are the `y` and `x` dimensions of the model domain, respectively. The `x` dimension is the *cross-channel* dimension, Implementations using non-standard data should permute datasets to match this convention.
 
@@ -144,7 +144,7 @@ Manipulating Section data
 We are often interested in not only the spatiotemporal changes in the planform of the delta, but we want to know what is preserved in the subsurface.
 In DeltaMetrics, we refer to this preserved history as the "stratigraphy", and we provide a number of convenient routines for computing stratigraphy and analyzing the deposits.
 
-Importantly, the stratigraphy (or i.e., which voxels are preserved) is not computed by default when a Cube instance is created. 
+Importantly, the stratigraphy (or i.e., which voxels are preserved) is not computed by default when a Cube instance is created.
 We must directly tell the Cube instance to compute stratigraphy by specifying which variable contains the bed elevation history, because this history dictates preservation.
 
 Mainly, the API works by registering a section of a specified type, and
@@ -157,7 +157,7 @@ For a data cube, sections are most easily instantiated by the :obj:`~deltametric
 
     >>> rcm8cube.register_section('demo', dm.section.StrikeSection(y=10))
 
-which creates a section across a constant y-value ``==10``. 
+which creates a section across a constant y-value ``==10``.
 The path of any `Section` in the ``x-y`` plane can always be accessed via the ``.trace`` attribute.
 We can plot the trace on top the the final bed elevation to see where the section will be located.
 
@@ -225,7 +225,7 @@ We have implemented support for rapid stratigraphy computation for visualization
 These quick stratigraphy computations create a mesh of preserved elevations and fill this matrix with values sliced out of the ``t-x-y`` data.
 
 Notably, the full "boxy" stratigraphy computation is also quite fast.
-More on that below. 
+More on that below.
 Compute the quick stratigraphy as:
 
 .. doctest::
@@ -292,9 +292,9 @@ The following are currently implemented.
     >>> _strike = dm.section.StrikeSection(rcm8cube, y=18)
     >>> _path = dm.section.PathSection(rcm8cube, path=np.column_stack((np.linspace(50, 150, num=4000, dtype=np.int),
     ...                                                                np.linspace(10, 90, num=4000, dtype=np.int))))
-    
 
-The `Section` classes all inherit from the same ``BaseSection`` class, which means they mostly have the same options available to them, and have a common API. 
+
+The `Section` classes all inherit from the same ``BaseSection`` class, which means they mostly have the same options available to them, and have a common API.
 Each has unique instantiation arguments, though, which must be properly specified.
 
 .. doctest::
@@ -321,7 +321,7 @@ Computing and Manipulating Stratigraphy
 1) Does not consider volume of sediment filled by preserved-time indicies, 2) cannot be sliced by planform, 3) irregularity does not lend well to computation and other uses (hydrological studies).
 
 So, we want to be able to create what I refer to as "boxy" stratigraphy.
-This has been done in the past by "placing" values from, e.g., ``strata_sand_frac`` into stratigraphy. 
+This has been done in the past by "placing" values from, e.g., ``strata_sand_frac`` into stratigraphy.
 This requires full computation for any variable you want to examine though.
 Here, we use a method that computes boxy stratigraphy only once, then synthesizes the volume from
 the precomputed sparse indicies.
@@ -331,7 +331,7 @@ Here’s a simple example to demonstrate how we place data into the stratigraphy
 .. doctest::
 
     >>> ets = rcm8cube['eta'][:, 25, 120]  # a "real" slice of the model
-    
+
     >>> fig, ax = plt.subplots(figsize=(8, 4))
     >>> dm.plot.show_one_dimensional_trajectory_to_strata(ets, ax=ax, dz=0.25)
     >>> plt.show() #doctest: +SKIP
@@ -423,7 +423,7 @@ speed up computations if an array is being accessed over and over.
 
     fs = sc8cube.export_frozen_variable('strata_sand_frac')
     fe = sc8cube.Z  # exported volume does not have coordinate information!
-    
+
     fig, ax = plt.subplots(figsize=(10, 2))
     pcm = ax.pcolormesh(np.tile(np.arange(fs.shape[2]), (fs.shape[0], 1)),
        fe[:,10,:], fs[:,10,:], shading='auto',
@@ -468,8 +468,8 @@ methods to create new “variables” in the data cube which hold the binary
 values of the masks.
 
 Currently implemented `Masks`:
-  * ChannelMask 
-  * EdgeMask 
+  * ChannelMask
+  * EdgeMask
   * LandMask
   * ShorelineMask
   * WetMask
@@ -479,7 +479,7 @@ Currently implemented `Masks`:
 
     >>> # use a new cube, code is currently locked to sea_level==0
     >>> maskcube = dm.sample_data.cube.rcm8()
-    
+
     >>> # create the masks from variables in the cube
     >>> land_mask = dm.mask.LandMask(maskcube['eta'][-1, :, :])
     >>> wet_mask = dm.mask.WetMask(maskcube['eta'][-1, :, :])
@@ -497,7 +497,7 @@ Currently implemented `Masks`:
 
     >>> ax0.imshow(maskcube['eta'][-1, :, :]) #doctest: +SKIP
     >>> for i, m in enumerate([land_mask, wet_mask, channel_mask, centerline_mask, edge_mask, shore_mask]):
-    ...     axs[i].imshow(m.mask, cmap='gray') #doctest: +SKIP
+    ...     axs[i].imshow(m.mask[-1, :, :], cmap='gray') #doctest: +SKIP
     ...     axs[i].set_title(m.mask_type) #doctest: +SKIP
     >>> plt.show() #doctest: +SKIP
 
diff --git a/docs/source/pyplots/guides/userguide_masks_all_demo.py b/docs/source/pyplots/guides/userguide_masks_all_demo.py
index ab154a49..6b389178 100644
--- a/docs/source/pyplots/guides/userguide_masks_all_demo.py
+++ b/docs/source/pyplots/guides/userguide_masks_all_demo.py
@@ -18,5 +18,5 @@ ax0.imshow(rcm8cube['eta'][-1, :, :])
 
 for i, m in enumerate([land_mask, wet_mask, channel_mask,
                        centerline_mask, edge_mask, shore_mask]):
-    axs[i].imshow(m.mask, cmap='gray')
+    axs[i].imshow(m.mask[-1, :, :], cmap='gray')
     axs[i].set_title(m.mask_type)
diff --git a/docs/source/reference/mask/index.rst b/docs/source/reference/mask/index.rst
index 140b5aa5..caac6dd5 100644
--- a/docs/source/reference/mask/index.rst
+++ b/docs/source/reference/mask/index.rst
@@ -7,7 +7,13 @@ Masking operations
 The package makes available routines for masking out planview attributes of
 the delta.
 
-This class is defined in ``deltametrics.mask``. 
+This class is defined in ``deltametrics.mask``.
+
+.. note::
+
+   When `Mask` objects are instantiated with single time-slice x-y planform
+   data (i.e., `L x W` data), the returned mask is expanded along the first
+   dimension to return an object with dimensions `1 x L x W`. This ensures that all `Mask` objects will have the same number of dimensions, but may lead to confusion if writing your own functions. 
 
 
 Mask classes
@@ -15,7 +21,7 @@ Mask classes
 
 .. currentmodule:: deltametrics.mask
 
-.. autosummary:: 
+.. autosummary::
 	:toctree: ../../_autosummary
 
 	BaseMask


Test patch:
diff --git a/tests/test_mask.py b/tests/test_mask.py
index 975ce3d1..d0090950 100644
--- a/tests/test_mask.py
+++ b/tests/test_mask.py
@@ -75,6 +75,13 @@ class TestShorelineMask:
         with pytest.raises(TypeError):
             shoremask = mask.ShorelineMask('invalid')
 
+    def test_3d(self):
+        """Test with multiple time slices."""
+        # define the mask
+        shoremask = mask.ShorelineMask(rcm8cube['eta'][30:33, :, :])
+        # assert the shape
+        assert np.shape(shoremask.mask) == (3, 120, 240)
+
 
 class TestLandMask:
     """Tests associated with the mask.LandMask class."""
@@ -164,6 +171,13 @@ class TestLandMask:
         assert np.array_equal(landmask.mask,
                               landmask.mask.astype(bool)) is True
 
+    def test_3d(self):
+        """Test with multiple time slices."""
+        # define the mask
+        landmask = mask.LandMask(rcm8cube['eta'][30:33, :, :])
+        # assert the shape
+        assert np.shape(landmask.mask) == (3, 120, 240)
+
 
 class TestWetMask:
     """Tests associated with the mask.WetMask class."""
@@ -259,6 +273,13 @@ class TestWetMask:
         # assert - expect all values to be 0s
         assert np.all(wetmask.mask == 0)
 
+    def test_3d(self):
+        """Test with multiple time slices."""
+        # define the mask
+        wetmask = mask.WetMask(rcm8cube['eta'][30:33, :, :])
+        # assert the shape
+        assert np.shape(wetmask.mask) == (3, 120, 240)
+
 
 class TestChannelMask:
     """Tests associated with the mask.ChannelMask class."""
@@ -422,6 +443,14 @@ class TestChannelMask:
             channelmask = mask.ChannelMask('bad_velocity',
                                            rcm8cube['eta'][-1, :, :])
 
+    def test_3d(self):
+        """Test with multiple time slices."""
+        # define the mask
+        channelmask = mask.ChannelMask(rcm8cube['velocity'][30:33, :, :],
+                                       rcm8cube['eta'][30:33, :, :])
+        # assert the shape
+        assert np.shape(channelmask.mask) == (3, 120, 240)
+
 
 class TestEdgeMask:
     """Tests associated with the mask.EdgeMask class."""
@@ -564,6 +593,13 @@ class TestEdgeMask:
         # assert - expect all values to be 0s
         assert np.all(edgemask.mask == 0)
 
+    def test_3d(self):
+        """Test with multiple time slices."""
+        # define the mask
+        edgemask = mask.EdgeMask(rcm8cube['eta'][30:33, :, :])
+        # assert the shape
+        assert np.shape(edgemask.mask) == (3, 120, 240)
+
 
 class TestCenterlineMask:
     """Tests associated with the mask.CenterlineMask class."""
@@ -626,6 +662,15 @@ class TestCenterlineMask:
         # assert - expect all values to be 0s
         assert np.all(centerlinemask.mask == 0)
 
+    def test_3d(self):
+        """Test with multiple time slices."""
+        # define the mask
+        channelmask = mask.ChannelMask(rcm8cube['velocity'][30:33, :, :],
+                                       rcm8cube['eta'][30:33, :, :])
+        centerlinemask = mask.CenterlineMask(channelmask)
+        # assert the shape
+        assert np.shape(centerlinemask.mask) == (3, 120, 240)
+
     @pytest.mark.xfail()
     def test_rivamapDefaults(self):
         """Test rivamap extraction of centerlines."""
@@ -690,3 +735,10 @@ def test_init_OAM():
     # assertions
     assert oam.mask_type == 'test'
     assert np.all(oam.data == np.zeros((3, 3)))
+
+
+def test_wrong_size():
+    """Raise error if array with bad dimensions is used."""
+    bad_arr = np.arange(0, 5)
+    with pytest.raises(ValueError):
+        mask.BaseMask('bad1D', bad_arr)