velocity to strain rate with gauge length support and differentiate s…

…tep support (#399)



---------

Co-authored-by: derrick chambers <[email protected]>
Co-authored-by: Derrick Chambers <[email protected]>

dascore/core/patch.py

@@ -333,6 +333,7 @@ def iresample(self, *args, **kwargs):
     integrate = transform.integrate
     spectrogram = transform.spectrogram
     velocity_to_strain_rate = transform.velocity_to_strain_rate
+    velocity_to_strain_rate_edgeless = transform.velocity_to_strain_rate_edgeless
     dispersion_phase_shift = transform.dispersion_phase_shift
 
     # --- Method Namespaces

dascore/data_registry.txt

@@ -22,3 +22,4 @@ whale_1.hdf5 a09922969e740307bf26dc6ffa7fb9fbb834dc7cd7d4ced02c66b159fb1ce0cd ht
 febus_1.h5 73eba2b6e183b3bca51f8a1448c3b423c979d05ce6c18bfd7fb76b4f9bda5c0b https://github.com/dasdae/test_data/raw/master/das/febus_1.h5
 ap_sensing_1.hdf5 322429f2c44bed5dc72fb9a02f79bb0d3cb71048e93d906d3d24b0605b431b12 https://github.com/dasdae/test_data/raw/master/das/ap_sensing_1.hdf5
 silixa_h5_1.hdf5 d3f1b92b17ae2d00f900426e80d48964fb5a33b9480ef9805721ac756acd4a21 https://github.com/dasdae/test_data/raw/master/das/silixa_h5_1.hdf5
+deformation_rate_event_1.hdf5 be8574ae523de9b17d2a0f9e847f30301e1607e2076366e005cdb3a46b79f172 https://github.com/dasdae/test_data/raw/master/das/deformation_rate_event_1.hdf5

dascore/examples.py

@@ -240,6 +240,15 @@ def example_event_2():
     return out
 
 
+@register_func(EXAMPLE_PATCHES, key="deformation_rate_event_1")
+def deformation_rate_event_1():
+    """
+    An event recorded in an underground mine by a Terra15 unit.
+    """
+    path = fetch("deformation_rate_event_1.hdf5")
+    return dc.spool(path)[0]
+
+
 @register_func(EXAMPLE_PATCHES, key="ricker_moveout")
 def ricker_moveout(
     frequency=15,

dascore/transform/__init__.py

@@ -10,5 +10,5 @@
 from .fourier import dft, idft
 from .integrate import integrate
 from .spectro import spectrogram
-from .strain import velocity_to_strain_rate
+from .strain import velocity_to_strain_rate, velocity_to_strain_rate_edgeless
 from .dispersion import dispersion_phase_shift

dascore/transform/differentiate.py

@@ -2,12 +2,15 @@
 
 from __future__ import annotations
 
+from collections.abc import Sequence
 from operator import truediv
 
 import numpy as np
 
+from dascore.compat import is_array
 from dascore.constants import PatchType
-from dascore.utils.misc import iterate, optional_import
+from dascore.exceptions import ParameterError
+from dascore.utils.misc import broadcast_for_index, iterate, optional_import
 from dascore.utils.patch import (
     _get_data_units_from_dims,
     _get_dx_or_spacing_and_axes,
@@ -37,13 +40,44 @@ def _findiff_diff(data, axis, order, dx):
     return out
 
 
+def _get_diff(order, data, axes, dx_or_spacing):
+    """Get the diff of an array along specified axes."""
+    if order == 2:
+        new_data = _grad_diff(data, axes, dx_or_spacing)
+    else:
+        new_data = _findiff_diff(data, axes, order, dx_or_spacing)
+    return new_data
+
+
+def _strided_diff(order, patch, axes, dx_or_spacing, step):
+    """Calculate a strided differentiation along specified axes."""
+    if len(axes) > 1:
+        msg = "Step in patch.differentiate can only be used along one axis."
+        raise ParameterError(msg)
+    new_data = np.empty_like(patch.data)
+    dx_or_space = dx_or_spacing[0]
+    for step_ in range(step):
+        current_slice = slice(step_, None, step)
+        slicer = broadcast_for_index(patch.ndim, axes[0], current_slice)
+        # Need to either double DX or slice the coordinate spacing to
+        # account for the fact we are skipping some columns/rows.
+        if is_array(dx_or_space):
+            _dx_or_space = dx_or_space[current_slice]
+        else:
+            _dx_or_space = dx_or_space * step
+        sub = _get_diff(order, patch.data[slicer], axes, [_dx_or_space])
+        new_data[slicer] = sub
+    return new_data
+
+
 @patch_function()
 def differentiate(
     patch: PatchType,
-    dim: str | None,
-    order=2,
+    dim: str | Sequence[str] | None,
+    order: int = 2,
+    step: int = 1,
 ) -> PatchType:
-    """
+    r"""
     Calculate first derivative along dimension(s) using centeral diferences.
 
     The shape of the output patch is the same as the input patch. Derivative
@@ -60,29 +94,54 @@ def differentiate(
     order
         The order of the differentiation operator. Must be a possitive, even
         integar.
+    step
+        The number of columns/rows to skip for differention.
+        eg: an array of [a b c d e] uses b and d to calculate diff of c when
+        step = 1 and order = 2. When step = 2, a and e are used to calcualte
+        diff at c.
 
     Notes
     -----
     For order=2 (the default) numpy's gradient function is used. When
-    order != the optional package findiff must be installed in which case
+    order != 2, the optional package findiff must be installed in which case
     order is interpreted as accuracy ("order" means order of differention
     in that package).
 
+    The second order first derivate, for an evenly spaced coordiante,
+    is defined as:
+
+    $$
+    \hat{f}(x) = \frac{f(x + dx) - f(x - dx)}{2dx} + O({dx}^2)
+    $$
+
+    Where $\hat{f}(x)$ is the estiamted derivative of $f$ at $x$, $dx$ is
+    the sample spacing, and $O$ is the error term.
+
     Examples
     --------
     >>> import dascore as dc
     >>> patch = dc.get_example_patch()
+    >>>
+    >>> # Example 1
     >>> # 2nd order stencil for 1st derivative along time dimension
     >>> patch_diff_1 = patch.differentiate(dim='time', order=2)
+    >>>
+    >>> # Example 2
     >>> # 1st derivative along all dimensions
     >>> patch_diff_2 = patch.differentiate(dim=None)
+    >>>
+    >>> # Example 3
+    >>> # 1st derivative using a step size of 3. This spaces out the columns
+    >>> # or rows used for estimating the derivative.
+    >>> patch_diff_3 = patch.differentiate(dim="distance", step=3, order=2)
     """
     dims = iterate(dim if dim is not None else patch.dims)
     dx_or_spacing, axes = _get_dx_or_spacing_and_axes(patch, dims)
-    if order == 2:
-        new_data = _grad_diff(patch.data, axes, dx_or_spacing)
+    if step > 1:
+        new_data = _strided_diff(order, patch, axes, dx_or_spacing, step)
+    # This avoids an extra copy of the array so probably merits its own case.
     else:
-        new_data = _findiff_diff(patch.data, axes, order, dx_or_spacing)
+        new_data = _get_diff(order, patch.data, axes, dx_or_spacing)
     # update units
     data_units = _get_data_units_from_dims(patch, dims, truediv)
     attrs = patch.attrs.update(data_units=data_units)

dascore/transform/strain.py

@@ -1,8 +1,12 @@
-"""Transformations to strain rates."""
+"""Transformations related to strain."""
 
 from __future__ import annotations
 
+import warnings
+
+import dascore as dc
 from dascore.constants import PatchType
+from dascore.exceptions import ParameterError
 from dascore.transform.differentiate import differentiate
 from dascore.utils.patch import patch_function
 
@@ -13,31 +17,183 @@
 )
 def velocity_to_strain_rate(
     patch: PatchType,
-    gauge_multiple: int = 1,
+    step_multiple: int = 2,
+    gauge_multiple: None | int = None,
     order: int = 2,
 ) -> PatchType:
-    """
-    Convert velocity das data to strain rate.
+    r"""
+    Convert velocity DAS data to strain rate using central differences.
 
-    This is primarily used with Terra15 data and simply uses
-    [patch.differentiate](`dascore.transform.differentiate.differentiate`)
-    under the hood.
+    When order=2 and step_multiple=2 the derivative for non-edge values
+    is estimated by:
+
+    $$
+    \hat{f}(x) = \frac{f(x + (n/2)dx) - f(x - (n/2)dx)}{n dx}
+    $$
+
+    Where $dx$ is the distance step and $n$ is the step_multiple. Values for
+    edges are estimate with the appropriate forward/backward stencils so that
+    the shape of the output data match the input data. The equation
+    becomes more complicated for higher order stencils.
 
     Parameters
     ----------
     patch
-        A patch object containing das data. Note: attrs['data_type'] should be
+        A patch object containing DAS data. Note: attrs['data_type'] should be
         velocity.
+    step_multiple
+        The multiples of spatial sampling for the central averaging stencil.
+        Must be even as odd values result in a staggered grid.
     gauge_multiple
-        The multiples of spatial sampling to make the simulated gauge length.
+        Deprecated name for step_multiple. Use that instead.
     order
         The order for the finite difference 1st derivative stencil (accuracy).
+        It must be a multiple of 2
+
+    Examples
+    --------
+    >>> from contextlib import suppress
+    >>> import dascore as dc
+    >>> from dascore.exceptions import MissingOptionalDependencyError
+    >>> patch = dc.get_example_patch("deformation_rate_event_1")
+    >>>
+    >>> # Example 1
+    >>> # Estimate the strain rate with a gauge length twice the distance step.
+    >>> patch_strain = patch.velocity_to_strain_rate(step_multiple=2)
+    >>>
+    >>> # Example 2
+    >>> # Estimate the strain rate with a 10th order filter. This will raise
+    >>> # an exception if the package findiff is not installed.
+    >>> with suppress(MissingOptionalDependencyError):
+    ...     patch_strain = patch.velocity_to_strain_rate(order=10)
+    >>>
+    >>> # Example 3
+    >>> # Estimate strain rate with a 4th order filter and gauge length 4 times
+    >>> # the distance step.
+    >>> with suppress(MissingOptionalDependencyError):
+    ...     patch_strain = patch.velocity_to_strain_rate(step_multiple=4, order=4)
+
+    Notes
+    -----
+    This is primarily used with Terra15 data and simply uses
+    [patch.differentiate](`dascore.transform.differentiate.differentiate`)
+    under the hood to calculate spatial derivatives.
+
+    The output gauge length is equal to the step_multiple multuplied by the
+    spacing along the distance coordinate, although the concept of
+    gauge_length is more complex with higher oder filters. See
+    @yang2022filtering for more info.
+
+    See the [`velocity_to_strain_rate` note](docs/notes/velocity_to_strain_rate.qmd)
+    for more details on step_multiple and order effects.
+
+    The [edgeless](`dascore.Patch.velocity_to_strain_rate_edgeless`) version
+    of this function removes potential edge effects and supports even and odd
+    `step_multiple` values.
     """
-    assert gauge_multiple == 1, "only supporting 1 for now."
+    if gauge_multiple is not None:
+        msg = "gauge_multiple will be removed in the future. Use step_multiple."
+        warnings.warn(msg, DeprecationWarning)
+        step_multiple = gauge_multiple * 2
+
+    if step_multiple % 2 != 0:
+        msg = (
+            "Step_multiple must be even. Use velocity_to_strain_rate_edgeless "
+            "if odd step multiples are required."
+        )
+        raise ParameterError(msg)
+
     coord = patch.get_coord("distance", require_evenly_sampled=True)
     step = coord.step
-    patch = differentiate.func(patch, dim="distance", order=order)
+    patch = differentiate.func(
+        patch, dim="distance", order=order, step=step_multiple // 2
+    )
     new_attrs = patch.attrs.update(
-        data_type="strain_rate", gauge_length=step * gauge_multiple
+        data_type="strain_rate", gauge_length=step * step_multiple
     )
     return patch.update(attrs=new_attrs)
+
+
+@patch_function(
+    required_dims=("distance",),
+    required_attrs={"data_type": "velocity"},
+)
+def velocity_to_strain_rate_edgeless(
+    patch: PatchType,
+    step_multiple: int = 1,
+) -> PatchType:
+    r"""
+    Estimate strain-rate using central differences.
+
+    For odd step_multiple values this function estimates strain by taking a
+    staggered central difference according to:
+
+    $$
+    \hat{f} = \frac{f(x + n * dx/2) - f(x - n * dx/2)}{dx}
+    $$
+
+    Where $dx$ is the spatial sampling and $n$ is the step_multiple. As a result
+    the strain-rate between existing samples is estimated when $n$ is odd. Edges
+    (points where full central differences are not possible) are discarded in
+    the output.
+
+    Parameters
+    ----------
+    patch
+        A patch object containing DAS data. Note: attrs['data_type'] should be
+        velocity.
+    step_multiple
+        The number of spatial sampling steps to use in the central averaging.
+
+    Examples
+    --------
+    >>> import dascore as dc
+    >>> patch = dc.get_example_patch("deformation_rate_event_1")
+    >>>
+    >>> # Example 1
+    >>> # Estimate strain rate with a gauge length equal to distance step.
+    >>> patch_strain = patch.velocity_to_strain_rate_edgeless(step_multiple=1)
+    >>>
+    >>> # Example 2
+    >>> # Estimate strain rate with a gauge length 5 times the distance step.
+    >>> patch_strain = patch.velocity_to_strain_rate_edgeless(step_multiple=5)
+
+    Notes
+    -----
+    See [velocity_to_strain_rate](`dascore.Patch.velocity_to_strain_rate`)
+    for a similar function which does not change the shape of the patch.
+
+    The resulting gauge length is equal to the step_multiple multiplied by
+    the sampling along the distance dimension.
+
+    See the
+    [`velocity_to_strain_rate` note](docs/notes/velocity_to_strain_rate.qmd)
+    for more details on step_multiple and order effects.
+    """
+    coord = patch.get_coord("distance", require_evenly_sampled=True)
+    distance_step = coord.step
+    gauge_length = step_multiple * distance_step
+
+    data_1 = patch.select(distance=(step_multiple, None), samples=True).data
+    data_2 = patch.select(distance=(None, -step_multiple), samples=True).data
+    strain_rate = (data_1 - data_2) / gauge_length
+
+    # Need to get distance values between current ones.
+    dists = patch.get_array("distance")
+    new_dist = (dists[step_multiple:] + dists[:-step_multiple]) / 2
+    new_coords = patch.coords.update(distance=new_dist)
+
+    # Handle unit conversions.
+    new_data_units = None
+    data_units = dc.get_quantity(patch.attrs.data_units)
+    dist_units = dc.get_quantity(patch.get_coord("distance").units)
+    if data_units and dist_units:
+        new_data_units = data_units / dist_units
+
+    new_attrs = patch.attrs.update(
+        data_type="strain_rate",
+        gauge_length=distance_step * step_multiple,
+        data_units=new_data_units,
+    )
+
+    return dc.Patch(data=strain_rate, coords=new_coords, attrs=new_attrs)