velocity to strain rate with gauge length support

dascore/core/patch.py

@@ -331,6 +331,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_fd = transform.velocity_to_strain_rate_fd
     dispersion_phase_shift = transform.dispersion_phase_shift
 
     # --- Method Namespaces

dascore/io/silixah5/core.py

@@ -1,6 +1,7 @@
 """
 Core modules for Silixa H5 support.
 """
+
 from __future__ import annotations
 
 import numpy as np
@@ -16,9 +17,9 @@
 class SilixaPatchAttrs(dc.PatchAttrs):
     """Patch Attributes for Silixa hdf5 format."""
 
-    gauge_length: float = np.NaN
+    gauge_length: float = np.nan
     gauge_length_units: str = "m"
-    pulse_width: float = np.NaN
+    pulse_width: float = np.nan
     pulse_width_units: str = "ns"
 
 

dascore/io/silixah5/utils.py

@@ -1,6 +1,7 @@
 """
 Utility functions for AP sensing module.
 """
+
 import pandas as pd
 
 import dascore as dc

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_fd
 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,11 +40,42 @@ 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:
     """
     Calculate first derivative along dimension(s) using centeral diferences.
@@ -60,6 +94,11 @@ 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
     -----
@@ -72,17 +111,27 @@ def differentiate(
     --------
     >>> 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

@@ -2,6 +2,9 @@
 
 from __future__ import annotations
 
+import numpy as np
+
+import dascore as dc
 from dascore.constants import PatchType
 from dascore.transform.differentiate import differentiate
 from dascore.utils.patch import patch_function
@@ -17,27 +20,81 @@ def velocity_to_strain_rate(
     order: int = 2,
 ) -> PatchType:
     """
-    Convert velocity das data to strain rate.
+    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.
+    under the hood which allows higher order differentiating.
+    It doesn't change the shape of the array.
 
     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.
     gauge_multiple
         The multiples of spatial sampling to make the simulated gauge length.
     order
         The order for the finite difference 1st derivative stencil (accuracy).
+
+    See Also
+    --------
+    - [velocity_to_strain_rate_fd](`dascore.Patch.velocity_to_strain_rate_fd`)
     """
-    assert gauge_multiple == 1, "only supporting 1 for now."
     coord = patch.get_coord("distance", require_evenly_sampled=True)
-    step = coord.step
-    patch = differentiate.func(patch, dim="distance", order=order)
+    distance_step = coord.step
+    patch = differentiate.func(patch, dim="distance", order=order, step=gauge_multiple)
     new_attrs = patch.attrs.update(
-        data_type="strain_rate", gauge_length=step * gauge_multiple
+        data_type="strain_rate", gauge_length=distance_step * gauge_multiple
     )
     return patch.update(attrs=new_attrs)
+
+
+@patch_function(
+    required_dims=("distance",),
+    required_attrs={"data_type": "velocity"},
+)
+def velocity_to_strain_rate_fd(
+    patch: PatchType,
+    gauge_multiple: int = 1,
+) -> PatchType:
+    """
+    Convert velocity DAS data to strain rate using forward differences.
+
+    This is primarily used with Terra15 data and forward differences the data.
+    It does change the shape of the array.
+
+    Parameters
+    ----------
+    patch
+        A patch object containing DAS data. Note: attrs['data_type'] should be
+        velocity.
+    gauge_multiple : int, optional
+        The multiples of spatial sampling to make the simulated gauge length.
+
+    See Also
+    --------
+    - [velocity_to_strain_rate](`dascore.Patch.velocity_to_strain_rate`)
+    """
+    assert isinstance(gauge_multiple, int), "gauge_multiple must be an integer."
+    coord = patch.get_coord("distance", require_evenly_sampled=True)
+    distance_step = coord.step
+    gauge_length = gauge_multiple * distance_step
+
+    data_1 = patch.select(distance=(gauge_multiple, None), samples=True).data
+    data_2 = patch.select(distance=(None, -gauge_multiple), samples=True).data
+    strain_rate = (data_1 - data_2) / gauge_length
+
+    dist_1 = int(np.floor(gauge_multiple / 2))
+    dist_2 = -int(np.ceil(gauge_multiple / 2))
+    new_coords, _ = patch.coords.select(distance=(dist_1, dist_2), samples=True)
+
+    dist_unit = dc.get_quantity(patch.get_coord("distance").units)
+    new_data_units = dc.get_quantity(patch.attrs["data_units"]) / dist_unit
+    new_attrs = patch.attrs.update(
+        data_type="strain_rate",
+        gauge_length=distance_step * gauge_multiple,
+        data_units=new_data_units,
+    )
+
+    return dc.Patch(data=strain_rate, coords=new_coords, attrs=new_attrs)

tests/test_transform/test_differentiate.py

@@ -6,6 +6,7 @@
 import pytest
 
 import dascore as dc
+from dascore.exceptions import ParameterError
 from dascore.transform.differentiate import differentiate
 from dascore.units import get_quantity
 from dascore.utils.time import to_float
@@ -41,6 +42,20 @@ def linear_patch(random_patch):
     return random_patch.new(data=new_data)
 
 
+@pytest.fixture(scope="class")
+def linear_patch_monotonic_coord(random_patch):
+    """A patch increasing linearly along unevenly sampled distance dimension."""
+    state = np.random.RandomState(13)
+    x_values = np.cumsum(state.random(50))
+    dist_coord = dc.get_coord(values=x_values)
+    data = np.broadcast_to(x_values[..., None], (len(x_values), 10))
+    time = dc.get_coord(values=dc.to_datetime64(np.arange(10)))
+    dims = ("distance", "time")
+    coords = {"distance": dist_coord, "time": time}
+    patch = dc.Patch(data=data, coords=coords, dims=dims)
+    return patch
+
+
 class TestDifferentiateOrder2:
     """Simple tests for differentiation using order=2 (numpy.grad)."""
 
@@ -88,6 +103,46 @@ def test_uneven_spacing(self, wacky_dim_patch):
         assert np.allclose(expected, out.data, rtol=0.01)
 
 
+class TestStep:
+    """Tests for step != 1."""
+
+    def test_multi_dim_raises(self, random_patch):
+        """Can only use step on a single dimension."""
+        msg = "can only be used along one axis"
+        with pytest.raises(ParameterError, match=msg):
+            random_patch.differentiate(("time", "distance"), step=2)
+
+    def test_different_steps_comparable(self, linear_patch):
+        """Ensure the values are comparable with different step sizes."""
+        out1 = linear_patch.differentiate("distance", step=1)
+        out2 = linear_patch.differentiate("distance", step=2)
+        out3 = linear_patch.differentiate("distance", step=2)
+        assert np.allclose(out1.data, out2.data)
+        assert np.allclose(out2.data, out3.data)
+
+    def test_different_steps_monotonic_coord(self, linear_patch_monotonic_coord):
+        """Ensure using different steps still gives similar results."""
+        patch = linear_patch_monotonic_coord
+        out1 = patch.differentiate("distance", step=1)
+        out2 = patch.differentiate("distance", step=2)
+        out3 = patch.differentiate("distance", step=5)
+        assert np.allclose(out1.data, out2.data)
+        assert np.allclose(out2.data, out3.data)
+
+    def test_steps_with_order(self, linear_patch):
+        """Tests combinations of step and order."""
+        pytest.importorskip("findiff")
+        values = [
+            linear_patch.differentiate("distance"),
+            linear_patch.differentiate("distance", step=2, order=2),
+            linear_patch.differentiate("distance", step=5, order=6),
+            linear_patch.differentiate("distance", step=1, order=4),
+        ]
+        for num, patch in enumerate(values[:-1]):
+            next_patch = values[num + 1]
+            assert np.allclose(patch.data, next_patch.data)
+
+
 class TestCompareOrders:
     """Ensure differentiation with different orders returns similar results."""