feat: add method for returning binned spectrum values (#75)

Author: vegard-solum-4ss

src/waveresponse/_core.py

@@ -1680,14 +1680,77 @@ def grid(self, freq_hz=False, degrees=False):
 
         return freq, dirs, vals
 
+    def bingrid(self, freq_hz=False, degrees=False, complex_convert="rectangular"):
+        """
+        Return a copy of the spectrum's frequency and direction coordinates,
+        along with the corresponding binned spectrum values.
+
+        The spectrum values are interpolated and then integrated over their
+        respective direction bins, resulting in total energy per unit frequency.
+        This differs from the ``grid`` method, which returns the spectral density
+        values directly without bin integration.
+
+        Parameters
+        ----------
+        freq_hz : bool
+            Whether to return frequencies in hertz (Hz). If ``False``, angular
+            frequency in radians per second (rad/s) is used.
+        degrees : bool
+            Whether to return directions in degrees. If ``False``, radians are used.
+        complex_convert : str, optional
+            How to convert complex number grid values before interpolating. Should
+            be 'rectangular' or 'polar'. If 'rectangular' (default), complex values
+            are converted to rectangular form (i.e., real and imaginary part) before
+            interpolating. If 'polar', the values are instead converted to polar
+            form (i.e., amplitude and phase) before interpolating. The values are
+            converted back to complex form after interpolation.
+
+        Returns
+        -------
+        freq : ndarray, shape (N,)
+            1D array of frequency coordinates.
+        dirs : ndarray, shape (M,)
+            1D array of direction coordinates.
+        vals : ndarray, shape (N, M)
+            2D array of binned spectrum energy values (energy per unit frequency).
+            ``N = len(freq)``, ``M = len(dirs)``.
+        """
+
+        dirs_tmp = np.r_[
+            -2.0 * np.pi + self._dirs[-1], self._dirs, 2.0 * np.pi + self._dirs[0]
+        ]
+
+        dirs_bin = np.empty(self._dirs.size * 2 + 1, dtype=self._dirs.dtype)
+        dirs_bin[0::2] = dirs_tmp[:-1] + np.diff(dirs_tmp) / 2.0  # bin boundaries
+        dirs_bin[1::2] = self._dirs  # bin centers
+
+        interp_fun = self._interpolate_function(
+            complex_convert=complex_convert, method="linear", bounds_error=True
+        )
+
+        dirsnew, freqnew = np.meshgrid(dirs_bin, self._freq, indexing="ij", sparse=True)
+        vals_tmp = interp_fun((dirsnew, freqnew)).T
+
+        vals_binned = np.column_stack(
+            [
+                trapezoid(vals_tmp[:, i : i + 3], dirs_bin[i : i + 3], axis=1)
+                for i in range(0, 2 * len(self._dirs), 2)
+            ]
+        )
+
+        if freq_hz:
+            vals_binned *= 2.0 * np.pi
+
+        return self.freq(freq_hz=freq_hz), self.dirs(degrees=degrees), vals_binned
+
     def interpolate(
         self,
         freq,
         dirs,
         freq_hz=False,
         degrees=False,
+        complex_convert="rectangular",
         fill_value=0.0,
-        **kwargs,
     ):
         """
         Interpolate (linear) the spectrum values to match the given frequency and direction
@@ -1706,6 +1769,13 @@ def interpolate(
             If frequency is given in 'Hz'. If ``False``, 'rad/s' is assumed.
         degrees : bool
             If direction is given in 'degrees'. If ``False``, 'radians' is assumed.
+        complex_convert : str, optional
+            How to convert complex number grid values before interpolating. Should
+            be 'rectangular' or 'polar'. If 'rectangular' (default), complex values
+            are converted to rectangular form (i.e., real and imaginary part) before
+            interpolating. If 'polar', the values are instead converted to polar
+            form (i.e., amplitude and phase) before interpolating. The values are
+            converted back to complex form after interpolation.
         fill_value : float or None
             The value used for extrapolation (i.e., `freq` outside the bounds of
             the provided grid). If ``None``, values outside the frequency domain
@@ -1723,6 +1793,7 @@ def interpolate(
             dirs,
             freq_hz=freq_hz,
             degrees=degrees,
+            complex_convert=complex_convert,
             fill_value=fill_value,
         )
 

tests/test_core.py

@@ -3277,6 +3277,245 @@ def test_grid__hz_deg(self):
         np.testing.assert_array_almost_equal(dirs_out, dirs_expect)
         np.testing.assert_array_almost_equal(vals_out, vals_expect)
 
+    def test_bingrid(self):
+        freq_in = np.array([0.2, 0.4, 0.6, 0.8, 1.0])
+        dirs_in = np.array([0.0, 90.0, 180.0])
+        vals_in = np.ones((len(freq_in), len(dirs_in)))
+        spectrum = DirectionalSpectrum(
+            freq_in,
+            dirs_in,
+            vals_in,
+            freq_hz=True,
+            degrees=True,
+            clockwise=True,
+            waves_coming_from=True,
+        )
+
+        freq_out, dirs_out, vals_out = spectrum.bingrid(freq_hz=False, degrees=False)
+
+        freq_expect = 2.0 * np.pi * freq_in
+        dirs_expect = (np.pi / 180.0) * dirs_in
+        vals_expect = 1.0 / (2.0 * np.pi * (np.pi / 180.0)) * vals_in
+
+        vals_expect[:, 0] *= 135.0 * (np.pi / 180.0)
+        vals_expect[:, 1] *= 90.0 * (np.pi / 180.0)
+        vals_expect[:, 2] *= 135.0 * (np.pi / 180.0)
+
+        np.testing.assert_array_almost_equal(freq_out, freq_expect)
+        np.testing.assert_array_almost_equal(dirs_out, dirs_expect)
+        np.testing.assert_array_almost_equal(vals_out, vals_expect)
+
+    def test_bingrid2(self):
+        freq_in = np.array([0.2, 0.4, 0.6, 0.8, 1.0]) * 2.0 * np.pi
+        dirs_in = np.radians(np.array([0.0, 90.0, 180.0]))
+        vals_in = np.ones((len(freq_in), len(dirs_in)))
+        spectrum = DirectionalSpectrum(
+            freq_in,
+            dirs_in,
+            vals_in,
+            freq_hz=False,
+            degrees=False,
+            clockwise=True,
+            waves_coming_from=True,
+        )
+
+        freq_out, dirs_out, vals_out = spectrum.bingrid(freq_hz=False, degrees=False)
+
+        freq_expect = freq_in
+        dirs_expect = dirs_in
+        vals_expect = vals_in
+
+        vals_expect[:, 0] *= 135.0 * (np.pi / 180.0)
+        vals_expect[:, 1] *= 90.0 * (np.pi / 180.0)
+        vals_expect[:, 2] *= 135.0 * (np.pi / 180.0)
+
+        np.testing.assert_array_almost_equal(freq_out, freq_expect)
+        np.testing.assert_array_almost_equal(dirs_out, dirs_expect)
+        np.testing.assert_array_almost_equal(vals_out, vals_expect)
+
+    @pytest.mark.parametrize(
+        "dirs_in",
+        [
+            np.array([0.0, 90.0, 180.0, 270.0]),
+            np.array([10.0, 100.0, 190.0, 280.0]),
+            np.array([85.0, 175.0, 265.0, 355.0]),
+        ],
+    )
+    def test_bingrid_rads_rad(self, dirs_in):
+        freq_in = np.arange(0.0, 1, 0.1)
+        vals_in = np.column_stack(
+            [
+                np.zeros_like(freq_in),
+                np.ones_like(freq_in),
+                4 * np.ones_like(freq_in),
+                np.ones_like(freq_in),
+            ]
+        )
+        spectrum = DirectionalSpectrum(
+            freq_in,
+            dirs_in,
+            vals_in,
+            freq_hz=True,
+            degrees=True,
+            clockwise=True,
+            waves_coming_from=True,
+        )
+
+        freq_out, dirs_out, vals_out = spectrum.bingrid(freq_hz=False, degrees=False)
+
+        freq_expect = 2.0 * np.pi * freq_in
+        dirs_expect = (np.pi / 180.0) * dirs_in
+        vals_expect = np.column_stack(
+            [
+                22.5 * np.ones_like(freq_in),
+                112.5 * np.ones_like(freq_in),
+                292.5 * np.ones_like(freq_in),
+                112.5 * np.ones_like(freq_in),
+            ]
+        ) / (2.0 * np.pi)
+
+        np.testing.assert_array_almost_equal(freq_out, freq_expect)
+        np.testing.assert_array_almost_equal(dirs_out, dirs_expect)
+        np.testing.assert_array_almost_equal(vals_out, vals_expect)
+
+    @pytest.mark.parametrize(
+        "dirs_in",
+        [
+            np.array([0.0, 90.0, 180.0, 270.0]),
+            np.array([10.0, 100.0, 190.0, 280.0]),
+            np.array([85.0, 175.0, 265.0, 355.0]),
+        ],
+    )
+    def test_bingrid_hz_rad(self, dirs_in):
+        freq_in = np.arange(0.0, 1, 0.1)
+        vals_in = np.column_stack(
+            [
+                np.zeros_like(freq_in),
+                np.ones_like(freq_in),
+                4 * np.ones_like(freq_in),
+                np.ones_like(freq_in),
+            ]
+        )
+
+        spectrum = DirectionalSpectrum(
+            freq_in,
+            dirs_in,
+            vals_in,
+            freq_hz=True,
+            degrees=True,
+            clockwise=True,
+            waves_coming_from=True,
+        )
+
+        freq_out, dirs_out, vals_out = spectrum.bingrid(freq_hz=True, degrees=False)
+
+        freq_expect = freq_in
+        dirs_expect = dirs_in * (np.pi / 180.0)
+        vals_expect = np.column_stack(
+            [
+                22.5 * np.ones_like(freq_in),
+                112.5 * np.ones_like(freq_in),
+                292.5 * np.ones_like(freq_in),
+                112.5 * np.ones_like(freq_in),
+            ]
+        )
+
+        np.testing.assert_array_almost_equal(freq_out, freq_expect)
+        np.testing.assert_array_almost_equal(dirs_out, dirs_expect)
+        np.testing.assert_array_almost_equal(vals_out, vals_expect)
+
+    @pytest.mark.parametrize(
+        "dirs_in",
+        [
+            np.array([0.0, 90.0, 180.0, 270.0]),
+            np.array([10.0, 100.0, 190.0, 280.0]),
+            np.array([85.0, 175.0, 265.0, 355.0]),
+        ],
+    )
+    def test_bingrid_rads_deg(self, dirs_in):
+        freq_in = np.arange(0.0, 1, 0.1)
+        vals_in = np.column_stack(
+            [
+                np.zeros_like(freq_in),
+                np.ones_like(freq_in),
+                4 * np.ones_like(freq_in),
+                np.ones_like(freq_in),
+            ]
+        )
+
+        spectrum = DirectionalSpectrum(
+            freq_in,
+            dirs_in,
+            vals_in,
+            freq_hz=True,
+            degrees=True,
+            clockwise=True,
+            waves_coming_from=True,
+        )
+
+        freq_out, dirs_out, vals_out = spectrum.bingrid(freq_hz=False, degrees=True)
+
+        freq_expect = 2.0 * np.pi * freq_in
+        dirs_expect = dirs_in
+        vals_expect = np.column_stack(
+            [
+                22.5 * np.ones_like(freq_in),
+                112.5 * np.ones_like(freq_in),
+                292.5 * np.ones_like(freq_in),
+                112.5 * np.ones_like(freq_in),
+            ]
+        ) / (2.0 * np.pi)
+
+        np.testing.assert_array_almost_equal(freq_out, freq_expect)
+        np.testing.assert_array_almost_equal(dirs_out, dirs_expect)
+        np.testing.assert_array_almost_equal(vals_out, vals_expect)
+
+    @pytest.mark.parametrize(
+        "dirs_in",
+        [
+            np.array([0.0, 90.0, 180.0, 270.0]),
+            np.array([10.0, 100.0, 190.0, 280.0]),
+            np.array([85.0, 175.0, 265.0, 355.0]),
+        ],
+    )
+    def test_bingrid_hz_deg(self, dirs_in):
+        freq_in = np.arange(0.0, 1, 0.1)
+        vals_in = np.column_stack(
+            [
+                np.zeros_like(freq_in),
+                np.ones_like(freq_in),
+                4 * np.ones_like(freq_in),
+                np.ones_like(freq_in),
+            ]
+        )
+        spectrum = DirectionalSpectrum(
+            freq_in,
+            dirs_in,
+            vals_in,
+            freq_hz=True,
+            degrees=True,
+            clockwise=True,
+            waves_coming_from=True,
+        )
+
+        freq_out, dirs_out, vals_out = spectrum.bingrid(freq_hz=True, degrees=True)
+
+        freq_expect = freq_in
+        dirs_expect = dirs_in
+        dirs_expect = dirs_in
+        vals_expect = np.column_stack(
+            [
+                22.5 * np.ones_like(freq_in),
+                112.5 * np.ones_like(freq_in),
+                292.5 * np.ones_like(freq_in),
+                112.5 * np.ones_like(freq_in),
+            ]
+        )
+
+        np.testing.assert_array_almost_equal(freq_out, freq_expect)
+        np.testing.assert_array_almost_equal(dirs_out, dirs_expect)
+        np.testing.assert_array_almost_equal(vals_out, vals_expect)
+
     def test_from_spectrum1d(self):
         freq = np.array([0.0, 0.5, 1.0])
         dirs = np.array([0.0, 90.0, 180.0, 270.0])