Repository: DeltaRCM/DeltaMetrics

Issue:
Implement options for `z`, `dz`, and `nz` at every entry point of the api, make sure consistent
implement options for `z`, `dz`, and `nz` at every entry point to the api, need to ensure the behavior is consistent, and a single default is set everywhere. Preferably a default would be set in only one location. 

Could make default `None` everywhere, and then if all `None`s into the `_determine_strat_coordinates()` function, set the default `dz` there?

Create some tests!

Discussion:

Patch:
diff --git a/deltametrics/strat.py b/deltametrics/strat.py
index f2bf71c7..ce9f7ba7 100644
--- a/deltametrics/strat.py
+++ b/deltametrics/strat.py
@@ -35,7 +35,7 @@ def compute_compensation(line1, line2):
     pass
 
 
-def compute_boxy_stratigraphy_volume(elev, prop, dz=None, z=None,
+def compute_boxy_stratigraphy_volume(elev, prop, z=None, dz=None, nz=None,
                                      return_cube=False):
     """Process t-x-y data volume to boxy stratigraphy volume.
 
@@ -57,20 +57,40 @@ def compute_boxy_stratigraphy_volume(elev, prop, dz=None, z=None,
     prop : :obj:`ndarray`
         The `t-x-y` ndarry of property data to process into the stratigraphy.
 
-    dz : :obj:`float`
-        Vertical resolution of stratigraphy, in meters.
+    z : :obj:`ndarray`, optional
+        Vertical coordinates for stratigraphy, in meters. Optional, and
+        mutually exclusive with :obj:`dz` and :obj:`nz`,
+        see :obj:`_determine_strat_coordinates` for complete description.
+
+    dz : :obj:`float`, optional
+        Vertical resolution of stratigraphy, in meters. Optional, and mutually
+        exclusive with :obj:`z` and :obj:`nz`,
+        see :obj:`_determine_strat_coordinates` for complete description.
+
+    nz : :obj:`int`, optional
+        Number of intervals for vertical coordinates of stratigraphy.
+        Optional, and mutually exclusive with :obj:`z` and :obj:`dz`,
+        see :obj:`_determine_strat_coordinates` for complete description.
 
     return_cube : :obj:`boolean`, optional
         Whether to return the stratigraphy as a
         :obj:`~deltametrics.cube.FrozenStratigraphyCube` instance. Default is
         to return an `ndarray` and :obj:`elevations` `ndarray`.
 
+        .. warning:: not implemented!
+
     Returns
     -------
-    stratigraphy :
+    stratigraphy_cube : :obj:`~deltametrics.cube.StratigraphyCube`
+        Not Implemented.
 
-    elevations :
+    stratigraphy : :obj:`ndarray`
+        A `z-x-y` `ndarray` with data from `prop` placed into voxels to fill
+        stratigraphy.
 
+    elevations : :obj:`ndarray`
+        A `z-x-y` `ndarray` with elevations for each voxel in the stratigraphy
+        array.
     """
     # verify dimensions
     if elev.shape != prop.shape:
@@ -80,7 +100,7 @@ def compute_boxy_stratigraphy_volume(elev, prop, dz=None, z=None,
 
     # compute preservation from low-level funcs
     strata, _ = _compute_elevation_to_preservation(elev)
-    z = _determine_strat_coordinates(elev, dz=dz, z=z)
+    z = _determine_strat_coordinates(elev, z=z, dz=dz, nz=nz)
     strata_coords, data_coords = _compute_preservation_to_cube(strata, z=z)
 
     # copy data out and into the stratigraphy based on coordinates
@@ -100,7 +120,7 @@ def compute_boxy_stratigraphy_volume(elev, prop, dz=None, z=None,
         return stratigraphy, elevations
 
 
-def compute_boxy_stratigraphy_coordinates(elev, dz=None, z=None,
+def compute_boxy_stratigraphy_coordinates(elev, z=None, dz=None, nz=None,
                                           return_cube=False,
                                           return_strata=False):
     """Process t-x-y data volume to boxy stratigraphy coordinates.
@@ -116,19 +136,57 @@ def compute_boxy_stratigraphy_coordinates(elev, dz=None, z=None,
     elev : :obj:`ndarray`
         The `t-x-y` ndarry of elevation data to determine stratigraphy.
 
-    prop : :obj:`ndarray`
-        The `t-x-y` ndarry of property data to process into the stratigraphy.
+    z : :obj:`ndarray`, optional
+        Vertical coordinates for stratigraphy, in meters. Optional, and
+        mutually exclusive with :obj:`dz` and :obj:`nz`,
+        see :obj:`_determine_strat_coordinates` for complete description.
+
+    dz : :obj:`float`, optional
+        Vertical resolution of stratigraphy, in meters. Optional, and mutually
+        exclusive with :obj:`z` and :obj:`nz`,
+        see :obj:`_determine_strat_coordinates` for complete description.
 
-    dz : :obj:`float`
-        Vertical resolution of stratigraphy, in meters.
+    nz : :obj:`int`, optional
+        Number of intervals for vertical coordinates of stratigraphy.
+        Optional, and mutually exclusive with :obj:`z` and :obj:`dz`,
+        see :obj:`_determine_strat_coordinates` for complete description.
+
+    return_cube : :obj:`boolean`, optional
+        Whether to return the stratigraphy as a
+        :obj:`~deltametrics.cube.FrozenStratigraphyCube` instance. Default is
+        `False`, do not return a cube.
+
+        .. warning:: not implemented!
+
+    return_strata : :obj:`boolean`, optional
+        Whether to return the stratigraphy elevations, as returned from
+        internal computations in :obj:`_compute_elevation_to_preservation`.
+        Default is `False`, do not return strata. 
 
     Returns
     -------
-    stratigraphy_cube :
+    stratigraphy_cube : :obj:`~deltametrics.cube.StratigraphyCube`
+        Not Implemented.
+
+    strata_coords : :obj:`ndarray`
+        An `N x 3` array of `z-x-y` coordinates where information is preserved
+        in the boxy stratigraphy. Rows in `strat_coords` correspond with rows
+        in `data_coords`. See :obj:`_compute_preservation_to_cube` for
+        computation implementation.
+
+    data_coords : :obj:`ndarray`
+        An `N x 3` array of `t-x-y` coordinates where information is to be
+        extracted from the data array. Rows in `data_coords` correspond with
+        rows in `strat_coords`.  See :obj:`_compute_preservation_to_cube` for
+        computation implementation.
+
+    strata : :obj:`ndarray`
+        A `t-x-y` `ndarry` of stratal surface elevations. Returned as third
+        argument only if `return_strata=True`.
     """
     # compute preservation from low-level funcs
     strata, _ = _compute_elevation_to_preservation(elev)
-    z = _determine_strat_coordinates(elev, dz=dz, z=z)
+    z = _determine_strat_coordinates(elev, z=z, dz=dz, nz=nz)
     strata_coords, data_coords = _compute_preservation_to_cube(strata, z=z)
 
     if return_cube:
@@ -540,15 +598,17 @@ def _determine_strat_coordinates(elev, z=None, dz=None, nz=None):
         where elevation is expected to along the zeroth axis.
 
     z : :obj:`ndarray`, optional
-        Array of Z values to use, returned unchanged if supplied.
+        Array of Z values to use for bounding intervals (i.e., points in `z`),
+        returned unchanged if supplied.
 
     dz : :obj:`float`, optional
         Interval in created Z array. Z array is created as
         ``np.arange(np.min(elev), np.max(elev)+dz, step=dz)``.
 
     nz : :obj:`int`, optional
-        Number of intervals in `z`. Z array is created as
-        ``np.linspace(np.min(elev), np.max(elev), num=nz, endpoint=True)``.
+        Number of *intervals* in `z`, that is, the number of points in `z` is
+        `nz+1`. Z array is created as ``np.linspace(np.min(elev), np.max
+        (elev), num=nz, endpoint=True)``.
     """
     # if nothing is supplied
     if (dz is None) and (z is None) and (nz is None):
@@ -577,6 +637,6 @@ def _determine_strat_coordinates(elev, z=None, dz=None, nz=None):
             raise _valerr
         max_dos = np.max(elev.data)
         min_dos = np.min(elev.data)
-        return np.linspace(min_dos, max_dos, num=nz, endpoint=True)
+        return np.linspace(min_dos, max_dos, num=nz+1, endpoint=True)
     else:
         raise RuntimeError('No coordinates determined. Check inputs.')


Test patch:
diff --git a/tests/test_strat.py b/tests/test_strat.py
index eb377cff..c9ac2a22 100644
--- a/tests/test_strat.py
+++ b/tests/test_strat.py
@@ -18,7 +18,7 @@ class TestComputeBoxyStratigraphyVolume:
     elev = golfcube['eta']
     time = golfcube['time']
 
-    def test_returns_volume_and_elevations(self):
+    def test_returns_volume_and_elevations_given_dz(self):
         s, e = strat.compute_boxy_stratigraphy_volume(
             self.elev, self.time, dz=0.05)
         assert s.ndim == 3
@@ -33,6 +33,13 @@ class TestComputeBoxyStratigraphyVolume:
         assert s.shape == e.shape
         assert np.all(e[:, 0, 0] == z)
 
+    def test_returns_volume_and_elevations_given_nz(self):
+        s, e = strat.compute_boxy_stratigraphy_volume(
+            self.elev, self.time, nz=33)
+        assert s.ndim == 3
+        assert s.shape == e.shape
+        assert s.shape[0] == 33 + 1
+
     @pytest.mark.xfail(raises=NotImplementedError,
                        strict=True, reason='Not yet developed.')
     def test_return_cube(self):
@@ -72,6 +79,18 @@ class TestComputeBoxyStratigraphyCoordinates:
     elev = golfcube['eta']
     time = golfcube['time']
 
+    def test_returns_sc_dc_given_dz(self):
+        # check for the warning when no option passed
+        with pytest.warns(UserWarning):
+            sc, dc = strat.compute_boxy_stratigraphy_coordinates(
+                self.elev)
+        # now specify dz
+        sc, dc = strat.compute_boxy_stratigraphy_coordinates(
+            self.elev, dz=0.05)
+        assert sc.shape == dc.shape
+        assert sc.shape[1] == 3
+        assert sc.shape[0] > 1  # don't know how big it will be
+
     def test_returns_sc_dc(self):
         sc, dc = strat.compute_boxy_stratigraphy_coordinates(
             self.elev, dz=0.05)
@@ -91,6 +110,14 @@ class TestComputeBoxyStratigraphyCoordinates:
         sc, dc = strat.compute_boxy_stratigraphy_coordinates(self.elev, z=z)
         assert np.min(sc[:, 0]) == 0
         assert np.max(sc[:, 0]) == 6
+        # check that the number of z values matches len(z)
+        assert np.unique(sc[:, 0]).shape[0] == len(z)
+
+    def test_returns_sc_dc_given_nz(self):
+        sc, dc = strat.compute_boxy_stratigraphy_coordinates(self.elev, nz=13)
+        assert np.min(sc[:, 0]) == 0
+        # check that the number of z values matches nz
+        assert np.unique(sc[:, 0]).shape[0] == 13
 
     @pytest.mark.xfail(raises=NotImplementedError,
                        strict=True, reason='Not yet developed.')
@@ -347,18 +374,18 @@ class TestDetermineStratCoordinates:
     def test_given_nz(self):
         e = np.array([0, 1, 1, 2, 1])
         z = strat._determine_strat_coordinates(e, nz=50)
-        assert len(z) == 50
+        assert len(z) == 50 + 1
         assert z[-1] == 2.00
         assert z[0] == 0
-        assert z[1] - z[0] == pytest.approx(2 / 49)
+        assert z[1] - z[0] == pytest.approx(2 / 50)  # delta_Z / nz
 
     def test_given_nz_negative_endpoint(self):
         e = np.array([0, 1, 1, 50, -1])
         z = strat._determine_strat_coordinates(e, nz=50)
-        assert len(z) == 50
+        assert len(z) == 50 + 1
         assert z[-1] == pytest.approx(50)
         assert z[0] == -1
-        assert z[1] - z[0] == pytest.approx(51 / 49)
+        assert z[1] - z[0] == pytest.approx(51 / 50)  # delta_Z / nz
 
     def test_given_nz_zero(self):
         e = np.array([0, 1, 1, 2, 1])