enh: add functionality for computing extreme values to DirectionalSpectrum class (#45)

Author: vegard-solum-4ss

docs/user_guide/concepts_utils/directional_spectrum.rst

@@ -72,3 +72,20 @@ method with the desired order, `n`.
     m2 = spectrum.moment(2)
 
     # Etc.
+
+Calculate the mean zero-crossing period, Tz:
+
+.. code-block:: python
+
+    spectrum.tz
+
+Calculate extreme values using the :meth:`~waveresponse.DirectionalSpectrum.extreme`
+method. The method takes two arguments: the duration of the process (in seconds),
+and a quantile, ``q``.
+
+.. code-block:: python
+
+    duration = 3 * 3600   # 3 hours
+
+    mpm = spectrum.extreme(duration, q=0.37)   # most probable maximum (MPM)
+    q90 = spectrum.extreme(duration, q=0.99)   # 90-th quantile

docs/user_guide/concepts_utils/wave_spectrum.rst

@@ -80,12 +80,6 @@ Calculate the wave peak period, Tp:
 
     wave.tp
 
-Calculate the mean crossing period, Tz:
-
-.. code-block:: python
-
-    wave.tz
-
 Calculate the wave peak direction:
 
 .. code-block:: python

src/waveresponse/_core.py

@@ -1433,6 +1433,56 @@ def moment(self, n, freq_hz=None):
         m_n = trapz((f**n) * spectrum, f)
         return m_n
 
+    @property
+    def tz(self):
+        """
+        Mean zero-crossing period, Tz, in 'seconds'.
+
+        Calculated from the zeroth- and second-order spectral moments according to:
+
+            ``tz = sqrt(m0 / m2)``
+
+        where the spectral moments are calculated by integrating over frequency in Hz.
+        """
+        m0 = self.moment(0, freq_hz=True)
+        m2 = self.moment(2, freq_hz=True)
+        return np.sqrt(m0 / m2)
+
+    def extreme(self, t, q=0.37):
+        """
+        Compute the q-th quantile extreme value (assuming a Gaussian process).
+
+        The extreme value, ``x``, is calculated according to:
+
+        ``x = sigma * sqrt(2 * ln((t / tz) / ln(1 / q)))``
+
+        where ``sigma`` is the standard deviation of the process, ``t`` is the duration
+        of the process, and ``q`` is the quantile. Setting ``q=0.37`` yields the
+        most probable maximum (MPM).
+
+        The extreme value indicates a level which the maximum value of the process
+        amplitudes will be below with a certain probability. Note that the method
+        only computes extreme values for the maxima, and does not consider the minima.
+
+        Parameters
+        ----------
+        t : float
+            Time/duration in seconds for which the of the process is observed.
+        q : float or array-like
+            Quantile or sequence of quantiles to compute. Must be between 0 and 1
+            (inclusive).
+
+        Returns
+        -------
+        x : float or array
+            Extreme value(s). During the given time period, the maximum value of
+            the process amplitudes will be below the returned value with a given
+            probability.
+
+        """
+        q = np.asarray_chkfinite(q)
+        return self.std() * np.sqrt(2.0 * np.log((t / self.tz) / np.log(1.0 / q)))
+
 
 class WaveSpectrum(DirectionalSpectrum):
     def __repr__(self):
@@ -1445,24 +1495,11 @@ def hs(self):
 
         Calculated from the zeroth-order spectral moment according to:
 
-            ``hs = 4.0 * np.sqrt(m0)``
+            ``hs = 4.0 * sqrt(m0)``
         """
         m0 = self.moment(0)
         return 4.0 * np.sqrt(m0)
 
-    @property
-    def tz(self):
-        """
-        Mean crossing period, Tz, (sometimes called the mean wave period) in 'seconds'.
-
-        Calculated from the zeroth- and second-order spectral moments according to:
-
-            ``tz = np.sqrt(m0 / m2)``
-        """
-        m0 = self.moment(0, freq_hz=False)
-        m2 = self.moment(2, freq_hz=True)
-        return np.sqrt(m0 / m2)
-
     @property
     def tp(self):
         """

tests/test_core.py

@@ -3284,6 +3284,24 @@ def test_moment_m2_rads(self):
         # not exactly same due to error in trapz for higher order functions
         assert m_out == pytest.approx(m_expect, rel=0.1)
 
+    def test_tz(self):
+        f0 = 0.0
+        f1 = 2.0
+
+        freq = np.linspace(f0, f1, 20)
+        dirs = np.arange(5, 360, 10)
+        vals = np.ones((len(freq), len(dirs)))
+        spectrum = DirectionalSpectrum(freq, dirs, vals, freq_hz=True, degrees=True)
+
+        tz_out = spectrum.tz
+
+        m0 = (0.0 - 360.0) * (f0 - f1)
+        m2 = (1.0 / 3.0) * (0.0 - 360.0) * (f0**3 - f1**3)
+
+        tz_expect = np.sqrt(m0 / m2)
+
+        assert tz_out == pytest.approx(tz_expect, rel=0.1)
+
     def test_from_grid(self):
         freq = np.linspace(0, 1.0, 10)
         dirs = np.linspace(0, 360.0, 15, endpoint=False)
@@ -3352,6 +3370,69 @@ def test_imag_raises(self, directional_spectrum):
         with pytest.raises(AttributeError):
             directional_spectrum.imag
 
+    def test_extreme_float(self):
+        f0 = 0.0
+        f1 = 2.0
+
+        freq = np.linspace(f0, f1, 20)
+        dirs = np.arange(5, 360, 10)
+        vals = np.ones((len(freq), len(dirs)))
+        spectrum = wr.DirectionalSpectrum(freq, dirs, vals, freq_hz=True, degrees=True)
+
+        sigma = spectrum.std()
+        tz = spectrum.tz
+
+        T = 360 * 24 * 60.0**2
+        q = 0.99
+        extreme_out = spectrum.extreme(T, q=q)
+
+        extreme_expect = sigma * np.sqrt(2.0 * np.log((T / tz) / np.log(1.0 / q)))
+
+        assert extreme_out == pytest.approx(extreme_expect)
+
+    def test_extreme_list(self):
+        f0 = 0.0
+        f1 = 2.0
+
+        freq = np.linspace(f0, f1, 20)
+        dirs = np.arange(5, 360, 10)
+        vals = np.ones((len(freq), len(dirs)))
+        spectrum = wr.DirectionalSpectrum(freq, dirs, vals, freq_hz=True, degrees=True)
+
+        sigma = spectrum.std()
+        tz = spectrum.tz
+
+        T = 360 * 24 * 60.0**2
+        q = [0.1, 0.5, 0.99]
+        extreme_out = spectrum.extreme(T, q=q)
+
+        extreme_expect = [
+            sigma * np.sqrt(2.0 * np.log((T / tz) / np.log(1.0 / q[0]))),
+            sigma * np.sqrt(2.0 * np.log((T / tz) / np.log(1.0 / q[1]))),
+            sigma * np.sqrt(2.0 * np.log((T / tz) / np.log(1.0 / q[2]))),
+        ]
+
+        np.testing.assert_array_almost_equal(extreme_out, extreme_expect)
+
+    def test_extreme_mpm(self):
+        f0 = 0.0
+        f1 = 2.0
+
+        freq = np.linspace(f0, f1, 20)
+        dirs = np.arange(5, 360, 10)
+        vals = np.ones((len(freq), len(dirs)))
+        spectrum = wr.DirectionalSpectrum(freq, dirs, vals, freq_hz=True, degrees=True)
+
+        sigma = spectrum.std()
+        tz = spectrum.tz
+
+        T = 360 * 24 * 60.0**2
+        extreme_out = spectrum.extreme(T, q=0.37)
+
+        extreme_expect = sigma * np.sqrt(2.0 * np.log(T / tz))
+
+        assert extreme_out == pytest.approx(extreme_expect, rel=1e-3)
+
 
 class Test_WaveSpectrum:
     def test__init__(self):
@@ -3428,24 +3509,6 @@ def test_hs(self):
 
         assert hs_out == pytest.approx(hs_expect)
 
-    def test_tz(self):
-        f0 = 0.0
-        f1 = 2.0
-
-        freq = np.linspace(f0, f1, 20)
-        dirs = np.arange(5, 360, 10)
-        vals = np.ones((len(freq), len(dirs)))
-        wave = WaveSpectrum(freq, dirs, vals, freq_hz=True, degrees=True)
-
-        tz_out = wave.tz
-
-        m0 = (0.0 - 360.0) * (f0 - f1)
-        m2 = (1.0 / 3.0) * (0.0 - 360.0) * (f0**3 - f1**3)
-
-        tz_expect = np.sqrt(m0 / m2)
-
-        assert tz_out == pytest.approx(tz_expect, rel=0.1)
-
     def test_tp_hz(self):
         freq = np.linspace(0, 2, 20)
         dirs = np.arange(5, 360, 10)