Add taper processing method (#154)

* start work on taper

* tests for bad inputs

* work taper

* add documentation for taper

* add test for bad window, don't shadow built-in type

dascore/core/patch.py

@@ -248,6 +248,7 @@ def to_xarray(self):
     interpolate = dascore.proc.interpolate
     normalize = dascore.proc.normalize
     standardize = dascore.proc.standardize
+    taper = dascore.proc.taper
 
     # --- Method Namespaces
     # Note: these can't be cached_property (from functools) or references

dascore/examples.py

@@ -153,7 +153,11 @@ def _random_patch_lat_lon():
 
 @register_func(EXAMPLE_PATCHES, key="example_event_1")
 def _example_event_1():
-    """Returns an example of a passive event recorded by DAS."""
+    """
+    Returns an example of a passive event recorded by DAS.
+
+    This event is from @stanvek2022fracture.
+    """
     path = fetch("example_dasdae_event_1.h5")
     return dc.spool(path)[0]
 

dascore/proc/__init__.py

@@ -7,3 +7,4 @@
 from .filter import median_filter, pass_filter, sobel_filter
 from .resample import decimate, interpolate, iresample, resample
 from .select import select
+from .taper import taper

dascore/proc/taper.py

@@ -0,0 +1,124 @@
+"""
+Processing for applying a taper.
+"""
+
+import numpy as np
+import pandas as pd
+from scipy.signal import windows  # the best operating system?
+
+from dascore.constants import PatchType
+from dascore.exceptions import ParameterError
+from dascore.utils.docs import compose_docstring
+from dascore.utils.misc import broadcast_for_index
+from dascore.utils.patch import get_dim_value_from_kwargs, patch_function
+
+TAPER_FUNCTIONS = dict(
+    barthann=windows.barthann,
+    bartlett=windows.bartlett,
+    blackman=windows.blackman,
+    blackmanharris=windows.blackmanharris,
+    bohman=windows.bohman,
+    hamming=windows.hamming,
+    hann=windows.hann,
+    nuttall=windows.nuttall,
+    parzen=windows.parzen,
+    triang=windows.triang,
+)
+
+
+def _get_taper_slices(patch, kwargs):
+    """Get slice for start/end of patch."""
+    dim, axis, value = get_dim_value_from_kwargs(patch, kwargs)
+    d_len = patch.shape[axis]
+    # TODO add unit support once patch_refactor branch lands
+    start, stop = np.broadcast_to(np.array(value), (2,))
+    slice1, slice2 = slice(None), slice(None)
+    if not pd.isnull(start):
+        start = int(start * d_len)
+        slice1 = slice(None, start)
+    if not pd.isnull(stop):
+        stop = int(d_len - stop * d_len)
+        slice2 = slice(stop, None)
+    return axis, (start, stop), slice1, slice2
+
+
+def _get_window_function(window_type):
+    """Get the window function to use for taper."""
+    # get taper function or raise if it isn't known.
+    if window_type not in TAPER_FUNCTIONS:
+        msg = (
+            f"'{window_type}' is not a known window type. "
+            f"Options are: {sorted(TAPER_FUNCTIONS)}"
+        )
+        raise ParameterError(msg)
+    func = TAPER_FUNCTIONS[window_type]
+    return func
+
+
+def _validate_windows(samps, shape, axis):
+    """Validate the the windows don't overlap or exceed dim len."""
+    samps = np.array([x for x in samps if not pd.isnull(x)])
+    if len(samps) > 1 and samps[0] > samps[1]:
+        msg = "Taper windows cannot overlap"
+        raise ParameterError(msg)
+    elif np.any(samps < 0) or np.any(samps > shape[axis]):
+        msg = "Total taper lengths exceed total dim length"
+        raise ParameterError(msg)
+
+
+@patch_function()
+@compose_docstring(taper_type=sorted(TAPER_FUNCTIONS))
+def taper(
+    patch: PatchType,
+    window_type: str = "hann",
+    **kwargs,
+) -> PatchType:
+    """
+    Taper the ends of the signal.
+
+    Parameters
+    ----------
+    patch
+        The patch instance.
+    window_type
+        The type of window to use For tapering. Supported Options are:
+            {taper_type}.
+    **kwargs
+        Used to specify the dimension along which to taper and the percentage
+        of total length of the dimension. If a single value is passed, the
+        taper will be applied to both ends. A length two tuple can specify
+        different values for each end, or no taper on one end.
+
+    Returns
+    -------
+    The tapered patch.
+
+    Examples
+    --------
+    >>> import dascore as dc
+    >>> patch = dc.get_example_patch() # generate example patch
+    >>> # Apply an Hanning taper to 5% of each end for time dimension.
+    >>> patch_taper1 = patch.taper(time=0.05, window_type="hann")
+    >>> # Apply a triangular taper to 10% of the start of the distance dimension.
+    >>> patch_taper2 = patch.taper(distance=(0.10, None), window_type='triang')
+    """
+    func = _get_window_function(window_type)
+    # get taper values in samples.
+    out = np.array(patch.data)
+    shape = out.shape
+    n_dims = len(out.shape)
+    axis, samps, start_slice, end_slice = _get_taper_slices(patch, kwargs)
+    _validate_windows(samps, shape, axis)
+    if samps[0] is not None:
+        window = func(2 * int(samps[0]))[: samps[0]]
+        # get indices window (which will broadcast) and data
+        data_inds = broadcast_for_index(n_dims, axis, start_slice)
+        window_inds = broadcast_for_index(n_dims, axis, slice(None), fill_none=True)
+        out[data_inds] = out[data_inds] * window[window_inds]
+    if samps[1] is not None:
+        diff = shape[axis] - samps[1]
+        window = func(2 * diff)[diff:]
+        data_inds = broadcast_for_index(n_dims, axis, end_slice)
+        window_inds = broadcast_for_index(n_dims, axis, slice(None), fill_none=True)
+        out[data_inds] = out[data_inds] * window[window_inds]
+    return patch.new(data=out)

dascore/utils/misc.py

@@ -142,6 +142,27 @@ def get_slice(array, cond=Optional[tuple]) -> slice:
     return slice(start, stop)
 
 
+def broadcast_for_index(
+    n_dims: int, axis: int, value: Union[slice, int], fill_none=False
+):
+    """
+    For a given shape of array, return empty slices except for slice axis.
+
+    Parameters
+    ----------
+    n_dims
+        The number of dimensions in the array that will be indexed.
+    axis
+        The axis number.
+    value
+        A slice object.
+    fill_none
+        If True, fill non axis dims with None, else slice(None)
+    """
+    fill = None if fill_none else slice(None)
+    return tuple(fill if x != axis else value for x in range(n_dims))
+
+
 def _get_sampling_rate(sampling_period):
     """
     Get a sampling rate as a float.

docs/index.qmd

@@ -96,10 +96,11 @@ out = (
 import dascore as dc
 patch = (
     dc.get_example_patch('example_event_1')
+    .taper(time=0.05)
     .pass_filter(time=(None, 300))
 )
 
-patch.viz.waterfall(show=True, scale=0.04);
+patch.viz.waterfall(show=True, scale=0.2);
 ```
 
 # Installation

docs/references.bib

@@ -34,3 +34,11 @@ @article{lindsey2021fiber
   year={2021},
   publisher={Annual Reviews}
 }
+
+@article{stanvek2022fracture,
+  title={Fracture Imaging Using DAS-Recorded Microseismic Events},
+  author={Stan{\v{e}}k, Franti{\v{s}}ek and Jin, Ge and Simmons, James},
+  journal={Frontiers in Earth Science},
+  volume={10},
+  year={2022}
+}

docs/tutorial/processing.qmd

@@ -25,7 +25,7 @@ patch = dc.examples.sin_wave_patch(
     frequency=[200, 10],
     channel_count=2,
 )
-_ = patch.viz.wiggle(show=True)
+patch.viz.wiggle(show=True);
 ```
 
 ## IIR filter
@@ -34,7 +34,7 @@ Next we decimate by 10x using IIR filter
 
 ```{python}
 decimated_iir = patch.decimate(time=10, filter_type='iir')
-_ = decimated_iir.viz.wiggle(show=True)
+decimated_iir.viz.wiggle(show=True);
 ```
 
 Notice the lowpass filter removed the 200 Hz signal and only
@@ -46,7 +46,7 @@ Next we decimate by 10x using FIR filter.
 
 ```{python}
 decimated_fir = patch.decimate(time=10, filter_type='fir')
-_ = decimated_fir.viz.wiggle(show=True)
+decimated_fir.viz.wiggle(show=True);
 ```
 
 ## No Filter
@@ -55,5 +55,75 @@ Next, we decimate without a filter to purposely induce aliasing.
 
 ```{python}
 decimated_no_filt = patch.decimate(time=10, filter_type=None)
-_ = decimated_no_filt.viz.wiggle(show=True)
+decimated_no_filt.viz.wiggle(show=True);
+```
+
+# Taper
+
+The [taper function](`dascore.Patch.taper`) is used to gently taper the edges of
+a patch dimension to zero. To see this, let's create a patch of all 1s and apply
+the taper function
+
+```{python}
+import numpy as np
+import dascore as dc
+
+_patch = dc.get_example_patch()
+patch_ones = _patch.new(data=np.ones_like(_patch.data))
+```
+
+The following code will apply a 10% taper, meaning the first and last 10% of the
+samples, along the time dimension.
+
+```{python}
+patch_ones.taper(time=0.1).viz.waterfall();
+```
+
+Of course, any dimensions works, and passing a tuple of values enables
+different amounts of tapering for each end of the dimension.
+
+
+```{python}
+patch_ones.taper(distance=(0.1, 0.3)).viz.waterfall();
+```
+
+Either end can also be `None` to indicated tapering on a single side of the patch
+dimension.
+
+```{python}
+patch_ones.taper(distance=(None, 0.1)).viz.waterfall();
+```
+
+
+## Managing Edge Effects
+One of the main applications of tapering is to mitigate edge artefacts, or "edge effects"
+associated with several other types of processing, including filtering. Consider the
+following:
+
+### Bandpass filtering without tapering
+
+```{python}
+import dascore as dc
+
+patch = (
+    dc.get_example_patch('example_event_1')
+    .pass_filter(time=(None, 300))
+)
+
+patch.viz.waterfall(show=True, scale=0.04);
+
+```
+
+### Bandpass filtering with tapering
+
+```{python}
+import dascore as dc
+
+patch = (
+    dc.get_example_patch('example_event_1')
+    .taper(time=0.05)
+    .pass_filter(time=(None, 300))
+)
+
+patch.viz.waterfall(show=True, scale=0.15);
 ```