From f4198c672c173f42899d8c27eb6b8b3a290ddee2 Mon Sep 17 00:00:00 2001
From: jayh <jayaram.hariharan@utexas.edu>
Date: Fri, 14 Jan 2022 22:01:11 -0600
Subject: [PATCH 1/7] revise input checking for mobility functions

This commit revises the input checking for the mobility functions.
Now the input checker can digest time-coordinate inputs or
time-index values, and it spits out dimensioned xarray objects
as well as converted (if needed) indices for the time values.

Revises the testing to get things working. Still need to do
the doc-strings and change the actual mobility functions to
stop returning arrays and return dimensioned xarrays...
---
 deltametrics/mobility.py | 265 +++++++++++++++++++++++++--------------
 tests/test_mobility.py   |  82 ++++++------
 2 files changed, 217 insertions(+), 130 deletions(-)

diff --git a/deltametrics/mobility.py b/deltametrics/mobility.py
index 37138c83..3bbc9f6a 100644
--- a/deltametrics/mobility.py
+++ b/deltametrics/mobility.py
@@ -15,7 +15,8 @@
 import xarray as xr
 
 
-def check_inputs(chmap, basevalues, time_window, landmap=None):
+def check_inputs(chmap, basevalues=None, basevalues_idx=None, window=None,
+                 window_idx=None, landmap=None):
     """
     Check the input variable values.
 
@@ -24,101 +25,175 @@ def check_inputs(chmap, basevalues, time_window, landmap=None):
 
     Parameters
     ----------
-    chmap : ndarray
-        A t-x-y shaped array with binary channel maps.
+    chmap : list, xarray.DataArray, numpy.ndarray
+        Either a list of 2-D deltametrics.mask, xarray.DataArray, or
+        numpy.ndarray objects, or a t-x-y 3-D xarray.DataArray or numpy.ndarray
+        with channel mask values.
 
-    basevalues : list
-        List of t indices to use as the base channel maps.
+    basevalues : list, int, float, optional
+        List of time values to use as the base channel map. (or single value)
 
-    time_window : int
-        Number of time slices (t indices) to use as the time lag
+    basevalues_idx : list, optional
+        List of time indices to use as the base channel map. (or single value)
 
-    landmap : ndarray, optional
-        A t-x-y shaped array with binary land maps. Or a x-y 2D array of a
-        binary base mask that contains 1s in the region to be analyzed and
-        0s elsewhere.
+    window : int, float, optional
+        Duration of time to use as the time lag (aka how far from the basemap
+        will be analyzed).
 
-    """
-    # check binary map types - pull out ndarray from mask or xarray if needed
-    if isinstance(chmap, np.ndarray) is True:
-        pass
-    elif issubclass(type(chmap), mask.BaseMask) is True:
-        raise TypeError(
-            'Cannot input a Mask directly to mobility metrics. '
-            'Use a list-of-masks instead.')
-        chmap = chmap.mask
-    elif isinstance(chmap, xr.core.dataarray.DataArray) is True:
-        chmap = chmap.values
-    elif isinstance(chmap, list):
-        # assume it is a timeseries of masks set into a list
-        _arrs = [msk._mask.astype(int) for msk in chmap]
-        chmap = np.array(_arrs)
-    else:
-        raise TypeError('chmap data type not understood.')
+    window_idx : int, float, optional
+        Duration of time in terms of indices (# of save states) to use as the
+        time lag.
 
-    if ((chmap == 0) | (chmap == 1)).all():
-        pass
-    else:
-        raise ValueError('chmap was not binary')
-
-    if landmap is not None:
-        if isinstance(landmap, np.ndarray) is True:
-            pass
-        elif isinstance(landmap, mask.BaseMask) is True:
-            landmap = landmap.mask
-        elif isinstance(landmap, xr.core.dataarray.DataArray) is True:
-            landmap = landmap.values
-        elif isinstance(landmap, list):
-            # assume it is a timeseries of masks set into a list
-            _arrs = [msk._mask.astype(int) for msk in landmap]
-            landmap = np.array(_arrs)
-        else:
-            raise TypeError('landmap data type not understood.')
-        if ((landmap == 0) | (landmap == 1)).all():
-            pass
+    landmap : list, xarray.DataArray, numpy.ndarray, optional
+        Either a list of 2-D deltametrics.mask, xarray.DataArray, or
+        numpy.ndarray objects, or a t-x-y 3-D xarray.DataArray or numpy.ndarray
+        with land mask values.
+
+    """
+    # handle input maps - try to convert some expected types
+    in_maps = {'chmap': chmap, 'landmap': landmap}
+    out_maps = {'chmap': None, 'landmap': None}
+    _skip = False
+    for key in in_maps.keys():
+        if in_maps[key] is None:
+            _skip = True
         else:
-            raise ValueError('landmap was not binary')
+            inmap = in_maps[key]
+            _skip = False
+        # first expect a list of objects and coerce them into xarray.dataarrays
+        if (isinstance(inmap, list)) and (_skip is False):
+            # depending on type convert to xarray.dataarray
+            if isinstance(inmap[0], np.ndarray) is True:
+                dims = ('time', 'x', 'y')  # assumes an ultimate t-x-y shape
+                if len(np.shape(inmap[0])) != 2:
+                    raise ValueError('Expected list to be of 2-D ndarrays.')
+                coords = {'time': np.arange(1),
+                          'x': np.arange(inmap[0].shape[0]),
+                          'y': np.arange(inmap[0].shape[1])}
+                _converted = [xr.DataArray(
+                    data=np.reshape(inmap[i],
+                                    (1, inmap[i].shape[0], inmap[i].shape[1])),
+                    coords=coords, dims=dims)
+                              for i in inmap]
+            elif issubclass(type(inmap[0]), mask.BaseMask) is True:
+                _converted = [i.mask for i in inmap]
+            elif isinstance(inmap[0], xr.DataArray) is True:
+                _converted = inmap
+            else:
+                raise TypeError('Type of objects in the input list is not '
+                                'a recognized type.')
+            # convert the list of xr.DataArrays into a single 3-D one
+            out_maps[key] = _converted[0]  # init the final 3-D DataArray
+            for j in range(1, len(_converted)):
+                # stack them up along the time array into a 3-D dataarray
+                out_maps[key] = xr.concat(
+                    (out_maps[key], _converted[j]), dim='time')
+
+        elif (isinstance(inmap, np.ndarray) is True) and \
+          (len(inmap.shape) == 3) and (_skip is False):
+            dims = ('time', 'x', 'y')  # assumes t-x-y orientation of array
+            coords = {'time': np.arange(inmap.shape[0]),
+                      'x': np.arange(inmap.shape[1]),
+                      'y': np.arange(inmap.shape[2])}
+            out_maps[key] = xr.DataArray(data=inmap, coords=coords, dims=dims)
+
+        elif (issubclass(type(inmap), mask.BaseMask) is True) and \
+          (_skip is False):
+            raise TypeError(
+                'Cannot input a Mask directly to mobility metrics. '
+                'Use a list-of-masks instead.')
+
+        elif (isinstance(inmap, xr.core.dataarray.DataArray) is True) and \
+          (len(inmap.shape) == 3) and (_skip is False):
+            out_maps[key] = inmap
+
+        elif _skip is False:
+            raise TypeError('Input mask data type or format not understood.')
+
+    # can't do this binary check for a list
+    # if ((chmap == 0) | (chmap == 1)).all():
+    #     pass
+    # else:
+    #     raise ValueError('chmap was not binary')
 
     # check basevalues and time_window types
-    if isinstance(basevalues, list) is False:
+    if (basevalues is not None) and \
+      (isinstance(basevalues, list) is False) and \
+      (isinstance(basevalues, int) is False) and \
+      (isinstance(basevalues, float) is False):
         try:
             basevalues = list(basevalues)
+            # convert to indices of the time dimension
+            basevalues = [
+                out_maps['chmap'].time[
+                    np.abs(out_maps['chmap'].time - i).argmin()]
+                for i in basevalues]
         except Exception:
-            raise TypeError('basevalues was not a list or list()-able obj.')
+            raise TypeError('basevalues was not a list or list-able obj.')
 
-    if isinstance(time_window, int) is False:
+    if (basevalues_idx is not None) and \
+      (isinstance(basevalues_idx, list) is False) and \
+      (isinstance(basevalues_idx, int) is False) and \
+      (isinstance(basevalues_idx, float) is False):
         try:
-            time_window = int(time_window)
+            basevalues_idx = list(basevalues_idx)
         except Exception:
-            raise TypeError('time_window was not an int or int()-able obj.')
-
-    # check map shapes (expect 3D t-x-y arrays of same size)
-    if len(chmap.shape) != 3:
-        raise ValueError('Shape of chmap not 3D (expect t-x-y).')
-    if landmap is not None:
-        if len(landmap.shape) != 3:
-            try:
-                tmp_landmap = np.empty(chmap.shape)
-                for i in range(0, tmp_landmap.shape[0]):
-                    tmp_landmap[i, :, :] = landmap
-                landmap = tmp_landmap
-            except Exception:
-                raise ValueError('Landmp does not match chmap, nor could it be'
-                                 ' cast into an array that would match.')
-
-        if np.shape(chmap) != np.shape(landmap):
+            raise TypeError('basevalues_idx was not a list or list-able obj.')
+
+    if (basevalues is not None) and (basevalues_idx is not None):
+        raise Warning(
+            'basevalues and basevalues_idx supplied, using `basevalues`.')
+        base_out = basevalues
+    elif (basevalues is None) and (basevalues_idx is not None):
+        base_out = basevalues_idx
+    elif (basevalues is not None) and (basevalues_idx is None):
+        base_out = basevalues
+    else:
+        raise ValueError('No basevalue or basevalue_idx supplied!')
+
+    if (window is not None) and \
+      (isinstance(window, int) is False) and \
+      (isinstance(window, float) is False):
+        raise TypeError('Input window type was not an integer or float.')
+    elif (window is not None):
+        # convert to index of the time dimension
+        window = out_maps['chmap'].time[
+            np.abs(out_maps['chmap'].time - window).argmin()]
+
+    if (window_idx is not None) and \
+      (isinstance(window_idx, int) is False) and \
+      (isinstance(window_idx, float) is False):
+        raise TypeError(
+            'Input window_idx type was not an integer or float.')
+
+    if (window is not None) and (window_idx is not None):
+        raise Warning(
+            'window and window_idx supplied, using `window`.')
+        win_out = window
+    elif (window is None) and (window_idx is not None):
+        win_out = window_idx
+    elif (window is not None) and (window_idx is None):
+        win_out = window
+    else:
+        raise ValueError('No window or window_idx supplied!')
+
+    # check map shapes align
+    if out_maps['landmap'] is not None:
+        if np.shape(out_maps['chmap']) != np.shape(out_maps['landmap']):
             raise ValueError('Shapes of chmap and landmap do not match.')
 
     # check that the combined basemap + timewindow does not exceed max t-index
-    Kmax = np.max(basevalues) + time_window
-    if Kmax > chmap.shape[0]:
+    Kmax = np.max(base_out) + win_out
+    if Kmax > out_maps['chmap'].shape[0]:
         raise ValueError('Largest basevalue + time_window exceeds max time.')
 
     # return the sanitized variables
-    return chmap, landmap, basevalues, time_window
+    return out_maps['chmap'], out_maps['landmap'], base_out, win_out
 
 
-def calculate_channel_decay(chmap, landmap, basevalues, time_window):
+def calculate_channel_decay(chmap, landmap,
+                            basevalues=None, basevalues_idx=None,
+                            window=None, window_idx=None):
     """
     Calculate channel decay (reduction in dry fraction).
 
@@ -152,10 +227,8 @@ def calculate_channel_decay(chmap, landmap, basevalues, time_window):
 
     """
     # sanitize the inputs first
-    chmap, landmap, basevalues, time_window = check_inputs(chmap,
-                                                           basevalues,
-                                                           time_window,
-                                                           landmap)
+    chmap, landmap, basevalues, time_window = check_inputs(
+        chmap, basevalues, basevalues_idx, window, window_idx, landmap)
 
     # initialize dry fraction array
     dryfrac = np.zeros((len(basevalues), time_window))
@@ -180,14 +253,16 @@ def calculate_channel_decay(chmap, landmap, basevalues, time_window):
             # subtract incremental map from dry map to get the new dry fraction
             base_dry -= chA_step
             # just want binary (1 = never channlized, 0 = channel visited)
-            base_dry[base_dry < 0] = 0  # no need to have negative values
+            base_dry.values[base_dry.values < 0] = 0  # no need to have negative values
             # store remaining dry fraction
             dryfrac[i, Nstep] = np.sum(base_dry) / np.sum(base_map)
 
     return dryfrac
 
 
-def calculate_planform_overlap(chmap, landmap, basevalues, time_window):
+def calculate_planform_overlap(chmap, landmap,
+                               basevalues=None, basevalues_idx=None,
+                               window=None, window_idx=None):
     """
     Calculate channel planform overlap.
 
@@ -222,10 +297,8 @@ def calculate_planform_overlap(chmap, landmap, basevalues, time_window):
 
     """
     # sanitize the inputs first
-    chmap, landmap, basevalues, time_window = check_inputs(chmap,
-                                                           basevalues,
-                                                           time_window,
-                                                           landmap)
+    chmap, landmap, basevalues, time_window = check_inputs(
+        chmap, basevalues, basevalues_idx, window, window_idx, landmap)
 
     # initialize D, phi and Ophi
     D = np.zeros((len(basevalues), time_window))
@@ -257,7 +330,9 @@ def calculate_planform_overlap(chmap, landmap, basevalues, time_window):
     return Ophi
 
 
-def calculate_reworking_fraction(chmap, landmap, basevalues, time_window):
+def calculate_reworking_fraction(chmap, landmap,
+                                 basevalues=None, basevalues_idx=None,
+                                 window=None, window_idx=None):
     """
     Calculate the reworking fraction.
 
@@ -292,10 +367,8 @@ def calculate_reworking_fraction(chmap, landmap, basevalues, time_window):
 
     """
     # sanitize the inputs first
-    chmap, landmap, basevalues, time_window = check_inputs(chmap,
-                                                           basevalues,
-                                                           time_window,
-                                                           landmap)
+    chmap, landmap, basevalues, time_window = check_inputs(
+        chmap, basevalues, basevalues_idx, window, window_idx, landmap)
 
     # initialize unreworked pixels (Nbt) and reworked fraction (fr)
     Nbt = np.zeros((len(basevalues), time_window))
@@ -341,7 +414,8 @@ def calculate_reworking_fraction(chmap, landmap, basevalues, time_window):
     return fr
 
 
-def calculate_channel_abandonment(chmap, basevalues, time_window):
+def calculate_channel_abandonment(chmap, basevalues=None, basevalues_idx=None,
+                                  window=None, window_idx=None):
     """
     Calculate channel abandonment.
 
@@ -372,9 +446,8 @@ def calculate_channel_abandonment(chmap, basevalues, time_window):
 
     """
     # sanitize the inputs first
-    chmap, landmap, basevalues, time_window = check_inputs(chmap,
-                                                           basevalues,
-                                                           time_window)
+    chmap, landmap, basevalues, time_window = check_inputs(
+        chmap, basevalues, basevalues_idx, window, window_idx)
     # initialize values
     PwetA = np.zeros((len(basevalues), time_window))
     chA_base = np.zeros_like(chmap[0, :, :])
@@ -392,7 +465,9 @@ def calculate_channel_abandonment(chmap, basevalues, time_window):
             # get the incremental map
             chA_step = chmap[Nbase+Nstep, :, :]
             # get the number of pixels that were abandonded
-            stepA = len(np.where(chA_base.flatten() > chA_step.flatten())[0])
+            stepA = len(
+                np.where(chA_base.values.flatten() >
+                         chA_step.values.flatten())[0])
             # store this number in the PwetA array for each transient map
             PwetA[i, Nstep] = stepA/baseA
 
@@ -422,6 +497,8 @@ def channel_presence(chmap):
         chans = chmap._mask
     elif isinstance(chmap, np.ndarray) is True:
         chans = chmap
+    elif isinstance(chmap, xr.DataArray) is True:
+        chans = chmap
     else:
         raise TypeError('chmap data type not understood.')
     channel_presence = np.sum(chans, axis=0) / chans.shape[0]
diff --git a/tests/test_mobility.py b/tests/test_mobility.py
index 20f4a3b2..eb45e6e1 100644
--- a/tests/test_mobility.py
+++ b/tests/test_mobility.py
@@ -3,6 +3,7 @@
 import sys
 import os
 import numpy as np
+import xarray as xr
 
 import deltametrics as dm
 from deltametrics import cube
@@ -32,12 +33,13 @@ def test_check_input_list_of_mask():
     """Test that a deltametrics.mask.BaseMask type can be used."""
     # call checker function
     assert isinstance(chmask, list)
-    chmap, landmap, basevalues, time_window = mob.check_inputs(chmask,
-                                                               [0], 1,
-                                                               landmask)
+    chmap, landmap, basevalues, time_window = mob.check_inputs(
+        chmask, basevalues_idx=[0], window_idx=1, landmap=landmask)
     # assert types
-    assert isinstance(chmap, np.ndarray) is True
-    assert isinstance(landmap, np.ndarray) is True
+    assert isinstance(chmap, xr.DataArray) is True
+    assert isinstance(landmap, xr.DataArray) is True
+    assert len(np.shape(chmap)) == 3
+    assert len(np.shape(landmap)) == 3
     assert isinstance(basevalues, list) is True
     assert isinstance(time_window, int) is True
 
@@ -58,8 +60,8 @@ def test_check_xarrays():
                                                                [0], 1,
                                                                land_xarr)
     # assert types
-    assert isinstance(chmap, np.ndarray) is True
-    assert isinstance(landmap, np.ndarray) is True
+    assert isinstance(chmap, xr.DataArray) is True
+    assert isinstance(landmap, xr.DataArray) is True
     assert isinstance(basevalues, list) is True
     assert isinstance(time_window, int) is True
 
@@ -67,24 +69,26 @@ def test_check_xarrays():
 def test_check_input_nolandmask():
     """Test that the check input can run without a landmap."""
     # call checker function
-    chmap, landmap, basevalues, time_window = mob.check_inputs(chmask,
-                                                               [0], 1)
+    chmap, landmap, basevalues, time_window = mob.check_inputs(
+        chmask, basevalues_idx=[0], window_idx=1)
     # assert types
-    assert isinstance(chmap, np.ndarray) is True
+    assert isinstance(chmap, xr.DataArray) is True
     assert landmap is None
     assert isinstance(basevalues, list) is True
     assert isinstance(time_window, int) is True
 
 
+@pytest.mark.xfail(reason='Removed binary check - to be added back later.')
 def test_check_input_notbinary_chmap():
     """Test that nonbinary channel map raises error."""
     ch_nonbin = np.zeros((3, 3, 3))
     ch_nonbin[0, 1, 1] = 1
     ch_nonbin[0, 1, 2] = 2
     with pytest.raises(ValueError):
-        mob.check_inputs(ch_nonbin, [0], 1)
+        mob.check_inputs(ch_nonbin, basevalues_idx=[0], window_idx=1)
 
 
+@pytest.mark.xfail(reason='Removed binary check - to be added back later.')
 def test_check_input_notbinary_landmap():
     """Test that nonbinary land map raises error."""
     land_nonbin = np.zeros((3, 3, 3))
@@ -92,61 +96,63 @@ def test_check_input_notbinary_landmap():
     land_nonbin[0, 1, 2] = 2
     ch_bin = np.zeros_like(land_nonbin)
     with pytest.raises(ValueError):
-        mob.check_inputs(ch_bin, [0], 1, land_nonbin)
+        mob.check_inputs(ch_bin, basevalues_idx=[0], window_idx=1,
+                         landmap=land_nonbin)
 
 
 def test_check_input_invalid_chmap():
     """Test that an invalid channel map input will throw an error."""
     with pytest.raises(TypeError):
-        mob.check_inputs('invalid', [0], 1, landmask)
+        mob.check_inputs('invalid', basevalues_idx=[0], window_idx=1,
+                         landmap=landmask)
 
 
 def test_check_input_invalid_landmap():
     """Test that an invalid landmap will throw an error."""
     with pytest.raises(TypeError):
-        mob.check_inputs(chmask, [0], 1, 'invalid')
-
-
-def test_check_inputs_basevalues_nonlist():
-    """Test that a non-list input as basevalues throws an error."""
-    with pytest.raises(TypeError):
-        mob.check_inputs(chmask, 0, 1)
+        mob.check_inputs(chmask, basevalues_idx=[0], window_idx=1,
+                         landmap='invalid')
 
 
 def test_check_input_invalid_time_window():
     """Test that a non-valid time_window throws an error."""
     with pytest.raises(TypeError):
-        mob.check_inputs(chmask, [0], 'invalid')
+        mob.check_inputs(chmask, basevalues_idx=[0], window_idx='invalid')
 
 
 def test_check_input_2dchanmask():
     """Test that an unexpected channel mask shape throws an error."""
-    with pytest.raises(ValueError):
-        mob.check_inputs(np.ones((5,)), [0], 1)
+    with pytest.raises(TypeError):
+        mob.check_inputs(np.ones((5,)), basevalues_idx=[0], window_idx=1)
 
 
 def test_check_input_diff_shapes():
     """Test that differently shaped channel and land masks throw an error."""
     with pytest.raises(ValueError):
-        mob.check_inputs(chmask, [0], 1, np.ones((3, 3, 3)))
+        mob.check_inputs(chmask, basevalues_idx=[0], window_idx=1,
+                         landmap=np.ones((3, 3, 3)))
 
 
 def test_check_input_1dlandmask():
     """Test a 1d landmask that will throw an error."""
-    with pytest.raises(ValueError):
-        mob.check_inputs(chmask, [0], 1, np.ones((10, 1)))
+    with pytest.raises(TypeError):
+        mob.check_inputs(chmask, basevalues_idx=[0], window_idx=1,
+                         landmap=np.ones((10, 1)))
 
 
 def test_check_input_exceedmaxvals():
     """Test a basevalue + time window combo that exceeds time indices."""
     with pytest.raises(ValueError):
-        mob.check_inputs(chmask, [0], 100)
+        mob.check_inputs(chmask, basevalues_idx=[0], window_idx=100)
 
 
+@pytest.mark.xfail(
+    reason='Removed this functionality - do we want to blindly expand dims?')
 def test_check_input_castlandmap():
-    """Test ability to case a 2D landmask to match 3D channelmap."""
-    chmap, landmap, bv, tw = mob.check_inputs(chmask, [0], 1,
-                                              landmask[0].mask[:, :])
+    """Test ability to cast a 2D landmask to match 3D channelmap."""
+    chmap, landmap, bv, tw = mob.check_inputs(chmask, basevalues_idx=[0],
+                                              window_idx=1,
+                                              landmap=landmask[0].mask[:, :])
     assert np.shape(chmap) == np.shape(landmap)
 
 
@@ -165,8 +171,8 @@ def test_check_input_castlandmap():
     chmap[i, :, :] = chmap[i-1, :, :].copy()
     chmap[i, -1*i, 1] = 0
     chmap[i, -1*i, 2] = 1
-# define the fluvial surface - entire 4x4 area
-fsurf = np.ones((4, 4))
+# define the fluvial surface - entire 4x4 area through all time
+fsurf = np.ones((5, 4, 4))
 # define the index corresponding to the basemap at time 0
 basevalue = [0]
 # define the size of the time window to use
@@ -175,27 +181,31 @@ def test_check_input_castlandmap():
 
 def test_dry_decay():
     """Test dry fraction decay."""
-    dryfrac = mob.calculate_channel_decay(chmap, fsurf, basevalue, time_window)
+    dryfrac = mob.calculate_channel_decay(
+        chmap, fsurf, basevalues_idx=basevalue, window_idx=time_window)
     assert np.all(dryfrac == np.array([[0.75, 0.6875, 0.625, 0.5625, 0.5]]))
 
 
 def test_planform_olap():
     """Test channel planform overlap."""
-    ophi = mob.calculate_planform_overlap(chmap, fsurf, basevalue, time_window)
+    ophi = mob.calculate_planform_overlap(
+        chmap, fsurf, basevalues_idx=basevalue, window_idx=time_window)
     assert pytest.approx(ophi == np.array([[1., 0.66666667, 0.33333333,
                                             0., -0.33333333]]))
 
 
 def test_reworking():
     """Test reworking index."""
-    fr = mob.calculate_reworking_fraction(chmap, fsurf, basevalue, time_window)
+    fr = mob.calculate_reworking_fraction(
+        chmap, fsurf, basevalues_idx=basevalue, window_idx=time_window)
     assert pytest.approx(fr == np.array([[0., 0.08333333, 0.16666667,
                                           0.25, 0.33333333]]))
 
 
 def test_channel_abandon():
     """Test channel abandonment function."""
-    ch_abandon = mob.calculate_channel_abandonment(chmap, basevalue, time_window)
+    ch_abandon = mob.calculate_channel_abandonment(
+        chmap, basevalues_idx=basevalue, window_idx=time_window)
     assert np.all(ch_abandon == np.array([[0., 0.25, 0.5, 0.75, 1.]]))
 
 

From 71b868314c4bf86f56f758e60de254f71bcecbd0 Mon Sep 17 00:00:00 2001
From: jayh <jayaram.hariharan@utexas.edu>
Date: Sat, 15 Jan 2022 11:42:58 -0600
Subject: [PATCH 2/7] update tests

---
 tests/test_mobility.py |  8 ++++----
 tests/test_utils.py    | 12 ++++++------
 2 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/tests/test_mobility.py b/tests/test_mobility.py
index eb45e6e1..2450274b 100644
--- a/tests/test_mobility.py
+++ b/tests/test_mobility.py
@@ -190,16 +190,16 @@ def test_planform_olap():
     """Test channel planform overlap."""
     ophi = mob.calculate_planform_overlap(
         chmap, fsurf, basevalues_idx=basevalue, window_idx=time_window)
-    assert pytest.approx(ophi == np.array([[1., 0.66666667, 0.33333333,
-                                            0., -0.33333333]]))
+    assert pytest.approx(ophi.values == np.array([[1., 0.66666667, 0.33333333,
+                                                   0., -0.33333333]]))
 
 
 def test_reworking():
     """Test reworking index."""
     fr = mob.calculate_reworking_fraction(
         chmap, fsurf, basevalues_idx=basevalue, window_idx=time_window)
-    assert pytest.approx(fr == np.array([[0., 0.08333333, 0.16666667,
-                                          0.25, 0.33333333]]))
+    assert pytest.approx(fr.values == np.array([[0., 0.08333333, 0.16666667,
+                                                 0.25, 0.33333333]]))
 
 
 def test_channel_abandon():
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 5d9f3566..b8e88d2a 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -123,8 +123,8 @@ def test_bad_inputs(self):
 
 def test_linear_fit():
     """Test linear curve fitting."""
-    ch_abandon = mob.calculate_channel_abandonment(chmap, basevalue,
-                                                   time_window)
+    ch_abandon = mob.calculate_channel_abandonment(
+        chmap, basevalues_idx=basevalue, window_idx=time_window)
     yfit, popts, cov, err = utils.curve_fit(ch_abandon, fit='linear')
     assert pytest.approx(yfit == np.array([4.76315477e-24, 2.50000000e-01,
                                            5.00000000e-01, 7.50000000e-01,
@@ -137,8 +137,8 @@ def test_linear_fit():
 
 def test_harmonic_fit():
     """Test harmonic curve fitting."""
-    ch_abandon = mob.calculate_channel_abandonment(chmap, basevalue,
-                                                   time_window)
+    ch_abandon = mob.calculate_channel_abandonment(
+        chmap, basevalues_idx=basevalue, window_idx=time_window)
     yfit, popts, cov, err = utils.curve_fit(ch_abandon, fit='harmonic')
     assert pytest.approx(yfit == np.array([-0.25986438, 0.41294455,
                                            0.11505591, 0.06683947,
@@ -151,8 +151,8 @@ def test_harmonic_fit():
 
 def test_invalid_fit():
     """Test invalid fit parameter."""
-    ch_abandon = mob.calculate_channel_abandonment(chmap, basevalue,
-                                                   time_window)
+    ch_abandon = mob.calculate_channel_abandonment(
+        chmap, basevalues_idx=basevalue, window_idx=time_window)
     with pytest.raises(ValueError):
         utils.curve_fit(ch_abandon, fit='invalid')
 

From 7d9449e70b1e96ef14a9db94f419d4f8e8cc7b9b Mon Sep 17 00:00:00 2001
From: jayh <jayaram.hariharan@utexas.edu>
Date: Sat, 15 Jan 2022 11:43:12 -0600
Subject: [PATCH 3/7] integrate time dimension into mobility functions

---
 deltametrics/mobility.py | 38 +++++++++++++++++++++++++++++---------
 1 file changed, 29 insertions(+), 9 deletions(-)

diff --git a/deltametrics/mobility.py b/deltametrics/mobility.py
index 3bbc9f6a..32e090bc 100644
--- a/deltametrics/mobility.py
+++ b/deltametrics/mobility.py
@@ -231,7 +231,12 @@ def calculate_channel_decay(chmap, landmap,
         chmap, basevalues, basevalues_idx, window, window_idx, landmap)
 
     # initialize dry fraction array
-    dryfrac = np.zeros((len(basevalues), time_window))
+    dims = ('base', 'time')  # base and time-lag dimensions
+    coords = {'base': np.arange(len(basevalues)),
+              'time': chmap.time[:time_window].values}
+    dryfrac = xr.DataArray(
+        data=np.zeros((len(basevalues), time_window)),
+        coords=coords, dims=dims)
 
     # loop through basevalues
     for i in range(0, len(basevalues)):
@@ -301,9 +306,14 @@ def calculate_planform_overlap(chmap, landmap,
         chmap, basevalues, basevalues_idx, window, window_idx, landmap)
 
     # initialize D, phi and Ophi
-    D = np.zeros((len(basevalues), time_window))
-    phi = np.zeros_like(D)
-    Ophi = np.zeros_like(D)
+    dims = ('base', 'time')  # base and time-lag dimensions
+    coords = {'base': np.arange(len(basevalues)),
+              'time': chmap.time[:time_window].values}
+    D = xr.DataArray(
+        data=np.zeros((len(basevalues), time_window)),
+        coords=coords, dims=dims)
+    phi = xr.zeros_like(D)
+    Ophi = xr.zeros_like(D)
 
     # loop through the base maps
     for j in range(0, len(basevalues)):
@@ -371,8 +381,13 @@ def calculate_reworking_fraction(chmap, landmap,
         chmap, basevalues, basevalues_idx, window, window_idx, landmap)
 
     # initialize unreworked pixels (Nbt) and reworked fraction (fr)
-    Nbt = np.zeros((len(basevalues), time_window))
-    fr = np.zeros_like(Nbt)
+    dims = ('base', 'time')  # base and time-lag dimensions
+    coords = {'base': np.arange(len(basevalues)),
+              'time': chmap.time[:time_window].values}
+    Nbt = xr.DataArray(
+        data=np.zeros((len(basevalues), time_window)),
+        coords=coords, dims=dims)
+    fr = xr.zeros_like(Nbt)
 
     # loop through the base maps
     for j in range(0, len(basevalues)):
@@ -449,9 +464,14 @@ def calculate_channel_abandonment(chmap, basevalues=None, basevalues_idx=None,
     chmap, landmap, basevalues, time_window = check_inputs(
         chmap, basevalues, basevalues_idx, window, window_idx)
     # initialize values
-    PwetA = np.zeros((len(basevalues), time_window))
-    chA_base = np.zeros_like(chmap[0, :, :])
-    chA_step = np.zeros_like(chmap[0, :, :])
+    dims = ('base', 'time')  # base and time-lag dimensions
+    coords = {'base': np.arange(len(basevalues)),
+              'time': chmap.time[:time_window].values}
+    PwetA = xr.DataArray(
+        data=np.zeros((len(basevalues), time_window)),
+        coords=coords, dims=dims)
+    chA_base = xr.zeros_like(chmap[0, :, :])
+    chA_step = xr.zeros_like(chmap[0, :, :])
 
     # loop through the basevalues
     for i in range(0, len(basevalues)):

From 6cb7c36889d9d23f34a29597a2d272f099b0415b Mon Sep 17 00:00:00 2001
From: jayh <jayaram.hariharan@utexas.edu>
Date: Sat, 15 Jan 2022 18:06:29 -0600
Subject: [PATCH 4/7] integrate time dimension for mobility

---
 deltametrics/mobility.py                 | 157 +++++++++++++++--------
 docs/source/reference/mobility/index.rst |   1 +
 2 files changed, 107 insertions(+), 51 deletions(-)

diff --git a/deltametrics/mobility.py b/deltametrics/mobility.py
index 32e090bc..a855aba2 100644
--- a/deltametrics/mobility.py
+++ b/deltametrics/mobility.py
@@ -117,24 +117,17 @@ def check_inputs(chmap, basevalues=None, basevalues_idx=None, window=None,
     #     raise ValueError('chmap was not binary')
 
     # check basevalues and time_window types
-    if (basevalues is not None) and \
-      (isinstance(basevalues, list) is False) and \
-      (isinstance(basevalues, int) is False) and \
-      (isinstance(basevalues, float) is False):
+    if (basevalues is not None):
         try:
-            basevalues = list(basevalues)
+            baselist = list(basevalues)
             # convert to indices of the time dimension
-            basevalues = [
-                out_maps['chmap'].time[
-                    np.abs(out_maps['chmap'].time - i).argmin()]
-                for i in basevalues]
+            basevalues = [np.argmin(
+                np.abs(out_maps['chmap'].time.data - i))
+                for i in baselist]
         except Exception:
             raise TypeError('basevalues was not a list or list-able obj.')
 
-    if (basevalues_idx is not None) and \
-      (isinstance(basevalues_idx, list) is False) and \
-      (isinstance(basevalues_idx, int) is False) and \
-      (isinstance(basevalues_idx, float) is False):
+    if (basevalues_idx is not None):
         try:
             basevalues_idx = list(basevalues_idx)
         except Exception:
@@ -157,8 +150,9 @@ def check_inputs(chmap, basevalues=None, basevalues_idx=None, window=None,
         raise TypeError('Input window type was not an integer or float.')
     elif (window is not None):
         # convert to index of the time dimension
-        window = out_maps['chmap'].time[
-            np.abs(out_maps['chmap'].time - window).argmin()]
+        _basetime = np.min(out_maps['chmap'].time.data)  # baseline time
+        _reltime = out_maps['chmap'].time.data - _basetime  # relative time
+        window = np.argmin(np.abs(_reltime - window)) + 1
 
     if (window_idx is not None) and \
       (isinstance(window_idx, int) is False) and \
@@ -204,19 +198,29 @@ def calculate_channel_decay(chmap, landmap,
 
     Parameters
     ----------
-    chmap : ndarray
-        A t-x-y shaped array with binary channel maps.
+    chmap : list, xarray.DataArray, numpy.ndarray
+        Either a list of 2-D deltametrics.mask, xarray.DataArray, or
+        numpy.ndarray objects, or a t-x-y 3-D xarray.DataArray or numpy.ndarray
+        with channel mask values.
 
-    landmap : ndarray
-        A t-x-y shaped array with binary land maps. Or an x-y shaped array
-        with the binary region representing the fluvial surface over which
-        the mobility metric should be computed.
+    landmap : list, xarray.DataArray, numpy.ndarray
+        Either a list of 2-D deltametrics.mask, xarray.DataArray, or
+        numpy.ndarray objects, or a t-x-y 3-D xarray.DataArray or numpy.ndarray
+        with land mask values.
 
-    basevalues : list
-        List of t indices to use as the base channel maps.
+    basevalues : list, int, float, optional
+        List of time values to use as the base channel map. (or single value)
 
-    time_window : int
-        Number of time slices (t indices) to use as the time lag
+    basevalues_idx : list, optional
+        List of time indices to use as the base channel map. (or single value)
+
+    window : int, float, optional
+        Duration of time to use as the time lag (aka how far from the basemap
+        will be analyzed).
+
+    window_idx : int, float, optional
+        Duration of time in terms of indices (# of save states) to use as the
+        time lag.
 
     Returns
     -------
@@ -278,19 +282,29 @@ def calculate_planform_overlap(chmap, landmap,
 
     Parameters
     ----------
-    chmap : ndarray
-        A t-x-y shaped array with binary channel maps
+    chmap : list, xarray.DataArray, numpy.ndarray
+        Either a list of 2-D deltametrics.mask, xarray.DataArray, or
+        numpy.ndarray objects, or a t-x-y 3-D xarray.DataArray or numpy.ndarray
+        with channel mask values.
+
+    landmap : list, xarray.DataArray, numpy.ndarray
+        Either a list of 2-D deltametrics.mask, xarray.DataArray, or
+        numpy.ndarray objects, or a t-x-y 3-D xarray.DataArray or numpy.ndarray
+        with land mask values.
+
+    basevalues : list, int, float, optional
+        List of time values to use as the base channel map. (or single value)
 
-    landmap : ndarray
-        A t-x-y shaped array with binary land maps. Or an x-y shaped array
-        with the binary region representing the fluvial surface over which
-        the mobility metric should be computed.
+    basevalues_idx : list, optional
+        List of time indices to use as the base channel map. (or single value)
 
-    basevalues : list
-        A list of values (t indices) to use for the base channel maps
+    window : int, float, optional
+        Duration of time to use as the time lag (aka how far from the basemap
+        will be analyzed).
 
-    time_window : int
-        Number of time slices (t values) to use for the transient maps
+    window_idx : int, float, optional
+        Duration of time in terms of indices (# of save states) to use as the
+        time lag.
 
     Returns
     -------
@@ -353,19 +367,29 @@ def calculate_reworking_fraction(chmap, landmap,
 
     Parameters
     ----------
-    chmap : ndarray
-        A t-x-y shaped array with binary channel maps
+    chmap : list, xarray.DataArray, numpy.ndarray
+        Either a list of 2-D deltametrics.mask, xarray.DataArray, or
+        numpy.ndarray objects, or a t-x-y 3-D xarray.DataArray or numpy.ndarray
+        with channel mask values.
+
+    landmap : list, xarray.DataArray, numpy.ndarray
+        Either a list of 2-D deltametrics.mask, xarray.DataArray, or
+        numpy.ndarray objects, or a t-x-y 3-D xarray.DataArray or numpy.ndarray
+        with land mask values.
 
-    landmap : ndarray
-        A t-x-y shaped array with binary land maps. Or an x-y shaped array
-        with the binary region representing the fluvial surface over which
-        the mobility metric should be computed.
+    basevalues : list, int, float, optional
+        List of time values to use as the base channel map. (or single value)
+
+    basevalues_idx : list, optional
+        List of time indices to use as the base channel map. (or single value)
 
-    basevalues : list
-        A list of values (t indices) to use for the base channel maps
+    window : int, float, optional
+        Duration of time to use as the time lag (aka how far from the basemap
+        will be analyzed).
 
-    time_window : int
-        Number of time slices (t values) to use for the transient maps
+    window_idx : int, float, optional
+        Duration of time in terms of indices (# of save states) to use as the
+        time lag.
 
     Returns
     -------
@@ -442,14 +466,24 @@ def calculate_channel_abandonment(chmap, basevalues=None, basevalues_idx=None,
 
     Parameters
     ----------
-    chmap : ndarray
-        A t-x-y shaped array with binary channel maps
+    chmap : list, xarray.DataArray, numpy.ndarray
+        Either a list of 2-D deltametrics.mask, xarray.DataArray, or
+        numpy.ndarray objects, or a t-x-y 3-D xarray.DataArray or numpy.ndarray
+        with channel mask values.
+
+    basevalues : list, int, float, optional
+        List of time values to use as the base channel map. (or single value)
 
-    basevalues : list
-        A list of values (t indices) to use for the base channel maps
+    basevalues_idx : list, optional
+        List of time indices to use as the base channel map. (or single value)
+
+    window : int, float, optional
+        Duration of time to use as the time lag (aka how far from the basemap
+        will be analyzed).
 
-    time_window : int
-        Number of time slices (t values) to use for the transient maps
+    window_idx : int, float, optional
+        Duration of time in terms of indices (# of save states) to use as the
+        time lag.
 
     Returns
     -------
@@ -512,6 +546,27 @@ def channel_presence(chmap):
     channel_presence : ndarray
         A x-y shaped array with the normalized channel presence values.
 
+    Examples
+    --------
+
+    .. plot::
+        :include-source:
+
+        >>> golfcube = dm.sample_data.golf()
+        >>> (x, y) = np.shape(golfcube['eta'][-1, ...])
+        >>> # calculate channel masks/presence over final 5 timesteps
+        >>> chmap = np.zeros((5, x, y))  # initialize channel map
+        >>> for i in np.arange(-5, 0):
+        ...     chmap[i, ...] = dm.mask.ChannelMask(
+        ...         golfcube['eta'][i, ...], golfcube['velocity'][i, ...],
+        ...         elevation_threshold=0, flow_threshold=0).mask
+        >>>
+        >>> fig, ax = plt.subplots(1, 2)
+        >>> golfcube.quick_show('eta', ax=ax[0])  # final delta
+        >>> p = ax[1].imshow(dm.mobility.channel_presence(chmap), cmap='Blues')
+        >>> dm.plot.append_colorbar(p, ax[1], label='Channelized Time')
+        >>> plt.show()
+
     """
     if isinstance(chmap, mask.ChannelMask) is True:
         chans = chmap._mask
diff --git a/docs/source/reference/mobility/index.rst b/docs/source/reference/mobility/index.rst
index dadc56a4..9e19ee96 100644
--- a/docs/source/reference/mobility/index.rst
+++ b/docs/source/reference/mobility/index.rst
@@ -27,3 +27,4 @@ The functions are defined in ``deltametrics.mobility``.
   calculate_planform_overlap
   calculate_reworking_fraction
   calculate_channel_abandonment
+  channel_presence

From daba9ba74caa38d306e625b674e1066a6494f998 Mon Sep 17 00:00:00 2001
From: jayh <jayaram.hariharan@utexas.edu>
Date: Sun, 16 Jan 2022 19:28:42 -0600
Subject: [PATCH 5/7] flexible dims, update tests

---
 deltametrics/mobility.py |  67 ++++++++++----
 tests/test_mobility.py   | 194 ++++++++++++++++++++++++++++++++++++---
 2 files changed, 233 insertions(+), 28 deletions(-)

diff --git a/deltametrics/mobility.py b/deltametrics/mobility.py
index a855aba2..ccd7fb43 100644
--- a/deltametrics/mobility.py
+++ b/deltametrics/mobility.py
@@ -74,7 +74,7 @@ def check_inputs(chmap, basevalues=None, basevalues_idx=None, window=None,
                     data=np.reshape(inmap[i],
                                     (1, inmap[i].shape[0], inmap[i].shape[1])),
                     coords=coords, dims=dims)
-                              for i in inmap]
+                              for i in range(len(inmap))]
             elif issubclass(type(inmap[0]), mask.BaseMask) is True:
                 _converted = [i.mask for i in inmap]
             elif isinstance(inmap[0], xr.DataArray) is True:
@@ -152,7 +152,7 @@ def check_inputs(chmap, basevalues=None, basevalues_idx=None, window=None,
         # convert to index of the time dimension
         _basetime = np.min(out_maps['chmap'].time.data)  # baseline time
         _reltime = out_maps['chmap'].time.data - _basetime  # relative time
-        window = np.argmin(np.abs(_reltime - window)) + 1
+        window = int(np.argmin(np.abs(_reltime - window)) + 1)
 
     if (window_idx is not None) and \
       (isinstance(window_idx, int) is False) and \
@@ -181,8 +181,11 @@ def check_inputs(chmap, basevalues=None, basevalues_idx=None, window=None,
     if Kmax > out_maps['chmap'].shape[0]:
         raise ValueError('Largest basevalue + time_window exceeds max time.')
 
+    # collect name of the first dimenstion (should be time assuming t-x-y)
+    dim0 = out_maps['chmap'].dims[0]
+
     # return the sanitized variables
-    return out_maps['chmap'], out_maps['landmap'], base_out, win_out
+    return out_maps['chmap'], out_maps['landmap'], base_out, win_out, dim0
 
 
 def calculate_channel_decay(chmap, landmap,
@@ -231,13 +234,13 @@ def calculate_channel_decay(chmap, landmap,
 
     """
     # sanitize the inputs first
-    chmap, landmap, basevalues, time_window = check_inputs(
+    chmap, landmap, basevalues, time_window, dim0 = check_inputs(
         chmap, basevalues, basevalues_idx, window, window_idx, landmap)
 
     # initialize dry fraction array
-    dims = ('base', 'time')  # base and time-lag dimensions
+    dims = ('base', dim0)  # base and time-lag dimensions
     coords = {'base': np.arange(len(basevalues)),
-              'time': chmap.time[:time_window].values}
+              dim0: chmap[dim0][:time_window].values}
     dryfrac = xr.DataArray(
         data=np.zeros((len(basevalues), time_window)),
         coords=coords, dims=dims)
@@ -316,13 +319,13 @@ def calculate_planform_overlap(chmap, landmap,
 
     """
     # sanitize the inputs first
-    chmap, landmap, basevalues, time_window = check_inputs(
+    chmap, landmap, basevalues, time_window, dim0 = check_inputs(
         chmap, basevalues, basevalues_idx, window, window_idx, landmap)
 
     # initialize D, phi and Ophi
-    dims = ('base', 'time')  # base and time-lag dimensions
+    dims = ('base', dim0)  # base and time-lag dimensions
     coords = {'base': np.arange(len(basevalues)),
-              'time': chmap.time[:time_window].values}
+              dim0: chmap[dim0][:time_window].values}
     D = xr.DataArray(
         data=np.zeros((len(basevalues), time_window)),
         coords=coords, dims=dims)
@@ -401,13 +404,13 @@ def calculate_reworking_fraction(chmap, landmap,
 
     """
     # sanitize the inputs first
-    chmap, landmap, basevalues, time_window = check_inputs(
+    chmap, landmap, basevalues, time_window, dim0 = check_inputs(
         chmap, basevalues, basevalues_idx, window, window_idx, landmap)
 
     # initialize unreworked pixels (Nbt) and reworked fraction (fr)
-    dims = ('base', 'time')  # base and time-lag dimensions
+    dims = ('base', dim0)  # base and time-lag dimensions
     coords = {'base': np.arange(len(basevalues)),
-              'time': chmap.time[:time_window].values}
+              dim0: chmap[dim0][:time_window].values}
     Nbt = xr.DataArray(
         data=np.zeros((len(basevalues), time_window)),
         coords=coords, dims=dims)
@@ -495,12 +498,12 @@ def calculate_channel_abandonment(chmap, basevalues=None, basevalues_idx=None,
 
     """
     # sanitize the inputs first
-    chmap, landmap, basevalues, time_window = check_inputs(
+    chmap, landmap, basevalues, time_window, dim0 = check_inputs(
         chmap, basevalues, basevalues_idx, window, window_idx)
     # initialize values
-    dims = ('base', 'time')  # base and time-lag dimensions
+    dims = ('base', dim0)  # base and time-lag dimensions
     coords = {'base': np.arange(len(basevalues)),
-              'time': chmap.time[:time_window].values}
+              dim0: chmap[dim0][:time_window].values}
     PwetA = xr.DataArray(
         data=np.zeros((len(basevalues), time_window)),
         coords=coords, dims=dims)
@@ -568,13 +571,45 @@ def channel_presence(chmap):
         >>> plt.show()
 
     """
+    tmp_chans = None  # instantiate
     if isinstance(chmap, mask.ChannelMask) is True:
         chans = chmap._mask
     elif isinstance(chmap, np.ndarray) is True:
-        chans = chmap
+        tmp_chans = chmap
     elif isinstance(chmap, xr.DataArray) is True:
         chans = chmap
+    elif isinstance(chmap, list) is True:
+        # convert to numpy.ndarray if possible
+        if (isinstance(chmap[0], np.ndarray) is True) \
+          or (isinstance(chmap[0], xr.DataArray) is True):
+            # init empty array
+            tmp_chans = np.zeros(
+                (len(chmap), chmap[0].squeeze().shape[0],
+                 chmap[0].squeeze().shape[1]))
+            # populate it
+            for i in range(len(chmap)):
+                if isinstance(chmap[0], xr.DataArray) is True:
+                    tmp_chans[i, ...] = chmap[i].data.squeeze()
+                else:
+                    tmp_chans[i, ...] = chmap[i].squeeze()
+        elif issubclass(type(chmap[0]), mask.BaseMask) is True:
+            tmp_chans = [i.mask for i in chmap]
+            # convert list to ndarray
+            chans = np.zeros(
+                (len(tmp_chans), tmp_chans[0].shape[1], tmp_chans[0].shape[2]))
+            for i in range(chans.shape[0]):
+                chans[i, ...] = tmp_chans[i]
+        else:
+            raise ValueError('Invalid values in the supplied list.')
     else:
         raise TypeError('chmap data type not understood.')
+    # if tmp_chans is a numpy.ndarray, dimensions are not known
+    if isinstance(tmp_chans, np.ndarray):
+        dims = ('time', 'x', 'y')  # assumes an ultimate t-x-y shape
+        coords = {'time': np.arange(tmp_chans.shape[0]),
+                  'x': np.arange(tmp_chans.shape[1]),
+                  'y': np.arange(tmp_chans.shape[2])}
+        chans = xr.DataArray(data=tmp_chans, coords=coords, dims=dims)
+    # calculation of channel presence is actually very simple
     channel_presence = np.sum(chans, axis=0) / chans.shape[0]
     return channel_presence
diff --git a/tests/test_mobility.py b/tests/test_mobility.py
index 2450274b..2c52c6f7 100644
--- a/tests/test_mobility.py
+++ b/tests/test_mobility.py
@@ -28,14 +28,44 @@
         dm.mask.LandMask(rcm8cube['eta'][i, :, :],
                          elevation_threshold=0))
 
+# make them into xarrays (list of xarrays)
+dims = ('time', 'x', 'y')  # assumes an ultimate t-x-y shape
+coords = {'time': np.arange(1),
+          'x': np.arange(chmask[0].mask.shape[0]),
+          'y': np.arange(chmask[0].mask.shape[1])}
+# channel masks
+ch_xarr = [xr.DataArray(
+    data=np.reshape(chmask[i].mask.data,
+                    (1, chmask[i].mask.shape[0], chmask[i].mask.shape[1])),
+    coords=coords, dims=dims)
+            for i in range(len(chmask))]
+# land masks
+land_xarr = [xr.DataArray(
+    data=np.reshape(landmask[i].mask.data,
+                    (1, landmask[i].mask.shape[0], landmask[i].mask.shape[1])),
+    coords=coords, dims=dims)
+            for i in range(len(landmask))]
+
+# convert them to ndarrays
+ch_arr = np.zeros((3, chmask[0].shape[0], chmask[0].shape[1]))
+land_arr = np.zeros((3, chmask[0].shape[0], chmask[0].shape[1]))
+ch_arr_list = []
+land_arr_list = []
+for i in range(3):
+    ch_arr[i, ...] = ch_xarr[i].data
+    land_arr[i, ...] = land_xarr[i].data
+    ch_arr_list.append(ch_xarr[i].data.squeeze())
+    land_arr_list.append(land_xarr[i].data.squeeze())
+
 
 def test_check_input_list_of_mask():
     """Test that a deltametrics.mask.BaseMask type can be used."""
     # call checker function
     assert isinstance(chmask, list)
-    chmap, landmap, basevalues, time_window = mob.check_inputs(
+    chmap, landmap, basevalues, time_window, dim0 = mob.check_inputs(
         chmask, basevalues_idx=[0], window_idx=1, landmap=landmask)
     # assert types
+    assert dim0 == 'time'
     assert isinstance(chmap, xr.DataArray) is True
     assert isinstance(landmap, xr.DataArray) is True
     assert len(np.shape(chmap)) == 3
@@ -48,18 +78,60 @@ def test_check_input_single_mask_error():
     """Test that a deltametrics.mask.BaseMask type can be used."""
     # call checker function
     with pytest.raises(TypeError, match=r'Cannot input a Mask .*'):
-        chmap, landmap, basevalues, time_window = mob.check_inputs(
-            chmask[0], [0], 1, landmask[0])
+        chmap, landmap, basevalues, time_window, dim0 = mob.check_inputs(
+            chmask[0], basevalues_idx=[0], window_idx=1,
+            landmap=landmask[0])
 
 
-@pytest.mark.xfail()
 def test_check_xarrays():
     """Test that an xarray.DataArray can be used as an input."""
     # call checker function
-    chmap, landmap, basevalues, time_window = mob.check_inputs(ch_xarr,
-                                                               [0], 1,
-                                                               land_xarr)
+    chmap, landmap, basevalues, time_window, dim0 = mob.check_inputs(
+        ch_xarr, basevalues_idx=[0], window_idx=1,
+        landmap=land_xarr)
+    # assert types
+    assert dim0 == 'time'
+    assert isinstance(chmap, xr.DataArray) is True
+    assert isinstance(landmap, xr.DataArray) is True
+    assert isinstance(basevalues, list) is True
+    assert isinstance(time_window, int) is True
+
+
+def test_check_list_ndarrays():
+    """Test that a list of numpy.ndarray can be used as an input."""
+    # call checker function
+    chmap, landmap, basevalues, time_window, dim0 = mob.check_inputs(
+        ch_arr_list, basevalues_idx=[0], window_idx=1,
+        landmap=land_arr_list)
+    # assert types
+    assert dim0 == 'time'
+    assert isinstance(chmap, xr.DataArray) is True
+    assert isinstance(landmap, xr.DataArray) is True
+    assert isinstance(basevalues, list) is True
+    assert isinstance(time_window, int) is True
+
+
+def test_check_ndarrays():
+    """Test that a numpy.ndarray can be used as an input."""
+    # call checker function
+    chmap, landmap, basevalues, time_window, dim0 = mob.check_inputs(
+        ch_arr, basevalues_idx=[0], window_idx=1,
+        landmap=land_arr)
+    # assert types
+    assert dim0 == 'time'
+    assert isinstance(chmap, xr.DataArray) is True
+    assert isinstance(landmap, xr.DataArray) is True
+    assert isinstance(basevalues, list) is True
+    assert isinstance(time_window, int) is True
+
+
+def test_check_basevalues_window():
+    """Test that basevalues and window inputs work."""
+    # call checker function
+    chmap, landmap, basevalues, time_window, dim0 = mob.check_inputs(
+        ch_arr, basevalues=[0], window=1, landmap=land_arr)
     # assert types
+    assert dim0 == 'time'
     assert isinstance(chmap, xr.DataArray) is True
     assert isinstance(landmap, xr.DataArray) is True
     assert isinstance(basevalues, list) is True
@@ -69,9 +141,10 @@ def test_check_xarrays():
 def test_check_input_nolandmask():
     """Test that the check input can run without a landmap."""
     # call checker function
-    chmap, landmap, basevalues, time_window = mob.check_inputs(
+    chmap, landmap, basevalues, time_window, dim0 = mob.check_inputs(
         chmask, basevalues_idx=[0], window_idx=1)
     # assert types
+    assert dim0 == 'time'
     assert isinstance(chmap, xr.DataArray) is True
     assert landmap is None
     assert isinstance(basevalues, list) is True
@@ -114,12 +187,42 @@ def test_check_input_invalid_landmap():
                          landmap='invalid')
 
 
+def test_check_input_invalid_basevalues():
+    """Test that a non-listable basevalues throws an error."""
+    with pytest.raises(TypeError):
+        mob.check_inputs(chmask, basevalues=0, window_idx='invalid')
+
+
+def test_check_input_invalid_basevalues_idx():
+    """Test that a non-listable basevalues_idx throws an error."""
+    with pytest.raises(TypeError):
+        mob.check_inputs(chmask, basevalues_idx=0, window_idx='invalid')
+
+
+def test_check_no_basevalues_error():
+    """No basevalues will throw an error."""
+    with pytest.raises(ValueError):
+        mob.check_inputs(chmask, window_idx='invalid')
+
+
 def test_check_input_invalid_time_window():
     """Test that a non-valid time_window throws an error."""
+    with pytest.raises(TypeError):
+        mob.check_inputs(chmask, basevalues_idx=[0], window='invalid')
+
+
+def test_check_input_invalid_time_window_idx():
+    """Test that a non-valid time_window_idx throws an error."""
     with pytest.raises(TypeError):
         mob.check_inputs(chmask, basevalues_idx=[0], window_idx='invalid')
 
 
+def test_check_no_time_window():
+    """Test that no time_window throws an error."""
+    with pytest.raises(ValueError):
+        mob.check_inputs(chmask, basevalues_idx=[0])
+
+
 def test_check_input_2dchanmask():
     """Test that an unexpected channel mask shape throws an error."""
     with pytest.raises(TypeError):
@@ -146,13 +249,26 @@ def test_check_input_exceedmaxvals():
         mob.check_inputs(chmask, basevalues_idx=[0], window_idx=100)
 
 
+def test_check_input_invalid_list():
+    """Test a wrong list."""
+    with pytest.raises(TypeError):
+        mob.check_inputs(['str', 5, 1.], basevalues_idx=[0], window_idx=1)
+
+
+def test_check_input_list_wrong_shape():
+    """Test list with wrongly shaped arrays."""
+    in_list = [np.zeros((5, 2, 1, 1)), np.zeros((5, 2, 2, 2))]
+    with pytest.raises(ValueError):
+        mob.check_inputs(in_list, basevalues_idx=[0], window_idx=100)
+
+
 @pytest.mark.xfail(
     reason='Removed this functionality - do we want to blindly expand dims?')
 def test_check_input_castlandmap():
     """Test ability to cast a 2D landmask to match 3D channelmap."""
-    chmap, landmap, bv, tw = mob.check_inputs(chmask, basevalues_idx=[0],
-                                              window_idx=1,
-                                              landmap=landmask[0].mask[:, :])
+    chmap, landmap, bv, tw, dim0 = mob.check_inputs(
+        chmask, basevalues_idx=[0], window_idx=1,
+        landmap=landmask[0].mask[:, :])
     assert np.shape(chmap) == np.shape(landmap)
 
 
@@ -171,6 +287,21 @@ def test_check_input_castlandmap():
     chmap[i, :, :] = chmap[i-1, :, :].copy()
     chmap[i, -1*i, 1] = 0
     chmap[i, -1*i, 2] = 1
+# alt typing for the input map
+# chmap xarrays
+dims = ('time', 'x', 'y')  # assumes an ultimate t-x-y shape
+coords = {'time': np.arange(5),
+          'x': np.arange(chmap.shape[1]),
+          'y': np.arange(chmap.shape[2])}
+chmap_xr = xr.DataArray(data=chmap, coords=coords, dims=dims)
+# lists
+chmap_xr_list = []
+chmap_list = []
+for i in range(5):
+    chmap_list.append(chmap[i, ...])
+    chmap_xr_list.append(chmap_xr[i, ...].squeeze())
+
+
 # define the fluvial surface - entire 4x4 area through all time
 fsurf = np.ones((5, 4, 4))
 # define the index corresponding to the basemap at time 0
@@ -210,9 +341,48 @@ def test_channel_abandon():
 
 
 def test_channel_presence():
-    """Test channel presence."""
+    """Test channel presence with a regular array."""
     chan_presence = mob.channel_presence(chmap)
     assert np.all(chan_presence == np.array([[0., 0.8, 0.2, 0.],
                                              [0., 0.6, 0.4, 0.],
                                              [0., 0.4, 0.6, 0.],
                                              [0., 0.2, 0.8, 0.]]))
+
+
+def test_channel_presence_xarray():
+    """Test channel presence with an xarray."""
+    chan_presence = mob.channel_presence(chmap_xr)
+    assert np.all(chan_presence == np.array([[0., 0.8, 0.2, 0.],
+                                             [0., 0.6, 0.4, 0.],
+                                             [0., 0.4, 0.6, 0.],
+                                             [0., 0.2, 0.8, 0.]]))
+
+
+def test_channel_presence_xarray_list():
+    """Test channel presence with a list of xarrays."""
+    chan_presence = mob.channel_presence(chmap_xr_list)
+    assert np.all(chan_presence == np.array([[0., 0.8, 0.2, 0.],
+                                             [0., 0.6, 0.4, 0.],
+                                             [0., 0.4, 0.6, 0.],
+                                             [0., 0.2, 0.8, 0.]]))
+
+
+def test_channel_presence_array_list():
+    """Test channel presence with a list of arrays."""
+    chan_presence = mob.channel_presence(chmap_list)
+    assert np.all(chan_presence == np.array([[0., 0.8, 0.2, 0.],
+                                             [0., 0.6, 0.4, 0.],
+                                             [0., 0.4, 0.6, 0.],
+                                             [0., 0.2, 0.8, 0.]]))
+
+
+def test_invalid_list_channel_presence():
+    """Test an invalid list."""
+    with pytest.raises(ValueError):
+        mob.channel_presence(['in', 'valid', 'list'])
+
+
+def test_invalid_type_channel_presence():
+    """Test an invalid input typing."""
+    with pytest.raises(TypeError):
+        mob.channel_presence('invalid type input')

From e3f485939f4e758a4ec6b388a6fab5ec9382d00c Mon Sep 17 00:00:00 2001
From: jayh <jayaram.hariharan@utexas.edu>
Date: Sun, 16 Jan 2022 19:58:21 -0600
Subject: [PATCH 6/7] float typing to be flexible

---
 deltametrics/mobility.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/deltametrics/mobility.py b/deltametrics/mobility.py
index ccd7fb43..72fc32e8 100644
--- a/deltametrics/mobility.py
+++ b/deltametrics/mobility.py
@@ -87,7 +87,7 @@ def check_inputs(chmap, basevalues=None, basevalues_idx=None, window=None,
             for j in range(1, len(_converted)):
                 # stack them up along the time array into a 3-D dataarray
                 out_maps[key] = xr.concat(
-                    (out_maps[key], _converted[j]), dim='time')
+                    (out_maps[key], _converted[j]), dim='time').astype(float)
 
         elif (isinstance(inmap, np.ndarray) is True) and \
           (len(inmap.shape) == 3) and (_skip is False):
@@ -95,7 +95,9 @@ def check_inputs(chmap, basevalues=None, basevalues_idx=None, window=None,
             coords = {'time': np.arange(inmap.shape[0]),
                       'x': np.arange(inmap.shape[1]),
                       'y': np.arange(inmap.shape[2])}
-            out_maps[key] = xr.DataArray(data=inmap, coords=coords, dims=dims)
+            out_maps[key] = \
+                xr.DataArray(
+                    data=inmap, coords=coords, dims=dims).astype(float)
 
         elif (issubclass(type(inmap), mask.BaseMask) is True) and \
           (_skip is False):
@@ -105,7 +107,7 @@ def check_inputs(chmap, basevalues=None, basevalues_idx=None, window=None,
 
         elif (isinstance(inmap, xr.core.dataarray.DataArray) is True) and \
           (len(inmap.shape) == 3) and (_skip is False):
-            out_maps[key] = inmap
+            out_maps[key] = inmap.astype(float)
 
         elif _skip is False:
             raise TypeError('Input mask data type or format not understood.')

From 34fb8d08eea407b51d8b521f970c9b8e5dbbea42 Mon Sep 17 00:00:00 2001
From: jayh <jayaram.hariharan@utexas.edu>
Date: Sun, 16 Jan 2022 19:58:33 -0600
Subject: [PATCH 7/7] documentation for mobility

---
 docs/source/reference/mobility/index.rst      |  2 +
 docs/source/reference/mobility/mobilityex.rst | 53 +++++++++++++++++++
 2 files changed, 55 insertions(+)
 create mode 100644 docs/source/reference/mobility/mobilityex.rst

diff --git a/docs/source/reference/mobility/index.rst b/docs/source/reference/mobility/index.rst
index 9e19ee96..010cf8c9 100644
--- a/docs/source/reference/mobility/index.rst
+++ b/docs/source/reference/mobility/index.rst
@@ -13,6 +13,8 @@ to the results of these functions is provided as well.
 
 The functions are defined in ``deltametrics.mobility``.
 
+.. include:: mobilityex.rst
+
 
 .. mobility functions
 .. ===================
diff --git a/docs/source/reference/mobility/mobilityex.rst b/docs/source/reference/mobility/mobilityex.rst
new file mode 100644
index 00000000..15070862
--- /dev/null
+++ b/docs/source/reference/mobility/mobilityex.rst
@@ -0,0 +1,53 @@
+
+Using the Mobility Functions
+----------------------------
+
+To use the mobility functions, you need a set of masks covering various time
+points from the model output. These can be in the form of a list of mask
+objects, numpy ndarrays, or xarrays. Or these can be 3-D arrays or xarrays.
+
+For this example a few masks will be generated and put into lists. Then these
+lists of masks will be used to compute metrics of channel mobility,
+which will then be plotted.
+
+.. note:: needs to be expanded!
+
+.. plot::
+    :include-source:
+    :context: reset
+
+    >>> golfcube = dm.sample_data.golf()
+    >>> channelmask_list = []
+    >>> landmask_list = []
+
+    >>> for i in range(50, 60):
+    ...     landmask_list.append(dm.mask.LandMask(
+    ...         golfcube['eta'][i, ...], elevation_threshold=0))
+    ...     channelmask_list.append(dm.mask.ChannelMask(
+    ...         golfcube['eta'][i, ...], golfcube['velocity'][i, ...],
+    ...         elevation_threshold=0, flow_threshold=0.3))
+
+    >>> dryfrac = dm.mobility.calculate_channel_decay(
+    ...     channelmask_list, landmask_list,
+    ...     basevalues_idx=[0, 1, 2], window_idx=5)
+    >>> Ophi = dm.mobility.calculate_planform_overlap(
+    ...     channelmask_list, landmask_list,
+    ...     basevalues_idx=[0, 1, 2], window_idx=5)
+    >>> fr = dm.mobility.calculate_reworking_fraction(
+    ...     channelmask_list, landmask_list,
+    ...     basevalues_idx=[0, 1, 2], window_idx=5)
+    >>> PwetA = dm.mobility.calculate_channel_abandonment(
+    ...     channelmask_list,
+    ...     basevalues_idx=[0, 1, 2], window_idx=5)
+
+    >>> fig, ax = plt.subplots(2, 2)
+    >>> dryfrac.plot.line(x='time', ax=ax[0, 0])
+    >>> ax[0, 0].set_title('Dry Fraction')
+    >>> Ophi.plot.line(x='time', ax=ax[0, 1])
+    >>> ax[0, 1].set_title('Overlap Values')
+    >>> fr.plot.line(x='time', ax=ax[1, 0])
+    >>> ax[1, 0].set_title('Reworked Fraction')
+    >>> PwetA.plot.line(x='time', ax=ax[1, 1])
+    >>> ax[1, 1].set_title('Abandoned Fraction')
+    >>> plt.tight_layout()
+    >>> plt.show()