implement recession analysis

python/events/src/hydrotools/events/baseflow/eckhardt.py

@@ -23,6 +23,7 @@
 """
 
 import numpy as np
+import numpy.lib as npl
 import numpy.typing as npt
 from numba import jit, float64
 
@@ -45,7 +46,23 @@ def linear_recession_analysis(
         recession_constant: float
             The recession parameter, a, from Eckhardt (2005, 2008).
     """
-    return 0.9
+    # Compute differences from time step to the next
+    differences = np.diff(series) < 0.0
+
+    # Inspect overlapping windows to identify recession periods
+    periods = npl.stride_tricks.sliding_window_view(
+        differences,
+        window_shape=window
+        )
+    recessions = np.all(periods, axis=1)
+
+    # Extract recession period central values
+    indices = np.where(recessions)[0] + window // 2
+    x = series[indices-1]
+    y = series[indices]
+
+    # Compute recession constant using regression through the origin
+    return np.sum(x * y) / np.sum(x * x)
 
 @jit(float64(float64[:], float64), nogil=True)
 def maximum_baseflow_analysis(
@@ -60,7 +77,8 @@ def maximum_baseflow_analysis(
         series: array-like, required
             A numpy array of streamflow values. Assumes last value in series is baseflow.
         recession_constant: float, required
-            Linear reservoir recession constant.
+            Linear reservoir recession constant, a, from Eckhardt (2005, 2008). Must be
+            between 0.0 and 1.0, typically between 0.85 and 0.95 for daily streamflow.
             
         Returns
         -------
@@ -93,9 +111,11 @@ def separate_baseflow(
         series: array-type, required
             A numpy array of streamflow values. Assumes first value in series is baseflow.
         recession_constant: float, required
-            Linear reservoir recession constant, a, from Eckhardt (2005, 2008).
+            Linear reservoir recession constant, a, from Eckhardt (2005, 2008). Must be
+            between 0.0 and 1.0, typically between 0.85 and 0.95 for daily streamflow.
         maximum_baseflow_index: float
-            The maximum baseflow index, $BFI_{max}$, from Eckhardt (2005, 2008).
+            The maximum baseflow index, $BFI_{max}$, from Eckhardt (2005, 2008). Must be
+            between 0.0 and 1.0, typically between 0.25 and 0.8 for daily streamflow.
             
         Returns
         -------

python/events/tests/test_baseflow.py

@@ -5,22 +5,17 @@
 import pandas as pd
 
 def test_linear_recession_analysis():
-    s = np.exp(-0.9 * np.linspace(0.0, 1.0, 100))
+    s = np.exp(-0.8 * np.linspace(0.0, 1.0, 100))
     a = bf.linear_recession_analysis(s)
 
-    assert a == 0.9
+    assert a >= 0.0
+    assert a <= 1.0
 
 def test_maximum_baseflow_analysis():
-    # s = np.exp(-0.9 * np.linspace(0.0, 1.0, 100))
     rng = np.random.default_rng()
     s = rng.normal(100.0, 10.0, 100)
 
-    from time import perf_counter
-    start = perf_counter()
     bfi_max = bf.maximum_baseflow_analysis(s, 0.9)
-    end = perf_counter()
-    print(f"{end-start:.6f} s")
-    print(bfi_max)
 
     assert bfi_max >= 0.0
     assert bfi_max <= 1.0
@@ -29,6 +24,5 @@ def test_separate_baseflow():
     rng = np.random.default_rng()
     s = rng.normal(100.0, 10.0, 100)
 
-    # Test numpy
     b = bf.separate_baseflow(s, 0.9, 0.5)
     assert b[0] == s[0]