diff --git a/deltametrics/mobility.py b/deltametrics/mobility.py index 37138c83..72fc32e8 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,174 @@ 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 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: + _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').astype(float) + + 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).astype(float) + + 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.astype(float) + + 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): try: - basevalues = list(basevalues) + baselist = list(basevalues) + # convert to indices of the time dimension + 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.') + raise TypeError('basevalues was not a list or list-able obj.') - if isinstance(time_window, int) is False: + if (basevalues_idx is not None): 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 + _basetime = np.min(out_maps['chmap'].time.data) # baseline time + _reltime = out_maps['chmap'].time.data - _basetime # relative time + window = int(np.argmin(np.abs(_reltime - window)) + 1) + + 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.') + # collect name of the first dimenstion (should be time assuming t-x-y) + dim0 = out_maps['chmap'].dims[0] + # return the sanitized variables - return chmap, landmap, basevalues, time_window + return out_maps['chmap'], out_maps['landmap'], base_out, win_out, dim0 -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). @@ -129,19 +203,29 @@ def calculate_channel_decay(chmap, landmap, basevalues, time_window): 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 ------- @@ -152,13 +236,16 @@ 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, dim0 = check_inputs( + chmap, basevalues, basevalues_idx, window, window_idx, landmap) # initialize dry fraction array - dryfrac = np.zeros((len(basevalues), time_window)) + dims = ('base', dim0) # base and time-lag dimensions + coords = {'base': np.arange(len(basevalues)), + dim0: chmap[dim0][: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)): @@ -180,14 +267,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. @@ -198,19 +287,29 @@ def calculate_planform_overlap(chmap, landmap, basevalues, time_window): 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 - A list of values (t indices) to use for 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 values) to use for the transient 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). + + window_idx : int, float, optional + Duration of time in terms of indices (# of save states) to use as the + time lag. Returns ------- @@ -222,15 +321,18 @@ 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, dim0 = check_inputs( + 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', dim0) # base and time-lag dimensions + coords = {'base': np.arange(len(basevalues)), + dim0: chmap[dim0][: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)): @@ -257,7 +359,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. @@ -268,19 +372,29 @@ def calculate_reworking_fraction(chmap, landmap, basevalues, time_window): 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 - A list of values (t indices) to use for 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 values) to use for the transient 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). + + window_idx : int, float, optional + Duration of time in terms of indices (# of save states) to use as the + time lag. Returns ------- @@ -292,14 +406,17 @@ 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, dim0 = 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)) - fr = np.zeros_like(Nbt) + dims = ('base', dim0) # base and time-lag dimensions + coords = {'base': np.arange(len(basevalues)), + dim0: chmap[dim0][: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)): @@ -341,7 +458,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. @@ -353,14 +471,24 @@ def calculate_channel_abandonment(chmap, basevalues, time_window): 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 - A list of values (t indices) to use for 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 values) to use for the transient 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). + + window_idx : int, float, optional + Duration of time in terms of indices (# of save states) to use as the + time lag. Returns ------- @@ -372,13 +500,17 @@ 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, dim0 = 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', dim0) # base and time-lag dimensions + coords = {'base': np.arange(len(basevalues)), + dim0: chmap[dim0][: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)): @@ -392,7 +524,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 @@ -417,12 +551,67 @@ 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() + """ + tmp_chans = None # instantiate if isinstance(chmap, mask.ChannelMask) is True: chans = chmap._mask elif isinstance(chmap, np.ndarray) is True: + 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/docs/source/reference/mobility/index.rst b/docs/source/reference/mobility/index.rst index dadc56a4..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 .. =================== @@ -27,3 +29,4 @@ The functions are defined in ``deltametrics.mobility``. calculate_planform_overlap calculate_reworking_fraction calculate_channel_abandonment + channel_presence 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() diff --git a/tests/test_mobility.py b/tests/test_mobility.py index 20f4a3b2..2c52c6f7 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 @@ -27,17 +28,48 @@ 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(chmask, - [0], 1, - landmask) + chmap, landmap, basevalues, time_window, dim0 = 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 dim0 == 'time' + 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 @@ -46,20 +78,62 @@ 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 isinstance(chmap, np.ndarray) is True - assert isinstance(landmap, np.ndarray) is True + 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 assert isinstance(time_window, int) is True @@ -67,24 +141,27 @@ 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, dim0 = mob.check_inputs( + chmask, basevalues_idx=[0], window_idx=1) # assert types - assert isinstance(chmap, np.ndarray) is True + assert dim0 == 'time' + 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 +169,106 @@ 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') + mob.check_inputs(chmask, basevalues_idx=[0], window_idx=1, + 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_inputs_basevalues_nonlist(): - """Test that a non-list input as basevalues throws an error.""" +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, 0, 1) + 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, [0], 'invalid') + 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(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) + + +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 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, dim0 = mob.check_inputs( + chmask, basevalues_idx=[0], window_idx=1, + landmap=landmask[0].mask[:, :]) assert np.shape(chmap) == np.shape(landmap) @@ -165,8 +287,23 @@ 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)) +# 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 basevalue = [0] # define the size of the time window to use @@ -175,34 +312,77 @@ 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) - assert pytest.approx(ophi == np.array([[1., 0.66666667, 0.33333333, - 0., -0.33333333]])) + ophi = mob.calculate_planform_overlap( + chmap, fsurf, basevalues_idx=basevalue, window_idx=time_window) + 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, basevalue, time_window) - assert pytest.approx(fr == np.array([[0., 0.08333333, 0.16666667, - 0.25, 0.33333333]])) + fr = mob.calculate_reworking_fraction( + chmap, fsurf, basevalues_idx=basevalue, window_idx=time_window) + assert pytest.approx(fr.values == 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.]])) 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') 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')