Create _streamflow_flow_indices.py

.gitignore

@@ -122,3 +122,6 @@ docs/variables.json
 
 # dask
 dask-worker-space
+
+# Apple
+.DS_Store

AUTHORS.rst

@@ -45,3 +45,4 @@ Contributors
 * Sascha Hofmann <[email protected]> `@saschahofmann <https://github.com/saschahofmann>`_
 * Javier Diez-Sierra <[email protected]> `@JavierDiezSierra <https://github.com/JavierDiezSierra>`_
 * Hui-Min Wang `@Hem-W <https://github.com/Hem-W>`
+* Faisal Mahmood <[email protected]> <[email protected]> `@faimahsho <https://github.com/faimahsho>`_

docs/references.bib

@@ -2152,3 +2152,35 @@ @article{droogers2002
 	url = {https://www.scopus.com/inward/record.uri?eid=2-s2.0-0036464359&doi=10.1023%2fA%3a1015508322413&partnerID=40&md5=7322aaa4c6874878f5b1dab3c73c1718},
 	type = {Article}
 }
+
+@article{addor2018,
+    author = {Addor, Nans and Nearing, Grey and Prieto, Cristina and Newman, A. and Le Vine, Nataliya and Clark, Martyn},
+    year = {2018},
+    month = {11},
+    pages = {},
+    title = {A Ranking of Hydrological Signatures Based on Their Predictability in Space},
+    journal = {Water Resources Research},
+    doi = {10.1029/2018WR022606}
+}
+
+@article{Clausen2000,
+    author = {Clausen, B and Biggs, Barry},
+    year = {2000},
+    month = {11},
+    pages = {184-197},
+    title = {Flow variables for ecological studies in temperate streams: Groupings based on covariance},
+    volume = {237},
+    journal = {Journal of Hydrology},
+    doi = {10.1016/S0022-1694(00)00306-1}
+}
+
+@article{Olden2003,
+    author = {Olden, Julian and Poff, N.},
+    year = {2003},
+    month = {03},
+    pages = {101 - 121},
+    title = {Redundancy and the Choice of Hydrologic Indices for Characterizing Stream Flow Regimes},
+    volume = {19},
+    journal = {River Research and Applications},
+    doi = {10.1002/rra.700}
+}

tests/test_hydrology.py

@@ -70,3 +70,32 @@ def test_simple(self, snw_series, pr_series):
         out = xci.melt_and_precip_max(snw, pr)
         np.testing.assert_array_equal(out, 2)
         assert out.units == "kg m-2"
+
+
+class TestFlowindex:
+    def test_simple(self, q_series):
+        a = np.ones(365) * 10
+        a[10:30] = 50
+        q = q_series(a)
+        out = xci.flow_index(q, 0.95)
+        np.testing.assert_array_equal(out, 5)
+
+
+class TestHighflowfrequency:
+    def test_simple(self, q_series):
+        a = np.zeros(365)
+        a[50:60] = 10
+        a[200:210] = 20
+        q = q_series(a)
+        out = xci.high_flow_frequency(q, 9, freq="YS")
+        np.testing.assert_array_equal(out, 20)
+
+
+class TestLowflowfrequency:
+    def test_simple(self, q_series):
+        a = np.ones(365) * 10
+        a[50:60] = 1
+        a[200:210] = 1
+        q = q_series(a)
+        out = xci.low_flow_frequency(q, 0.2, freq="YS")
+        np.testing.assert_array_equal(out, 20)

xclim/indices/_hydrology.py

@@ -2,16 +2,20 @@
 from __future__ import annotations
 
 import numpy as np
-import xarray
+import xarray as xr
 
 from xclim.core.calendar import get_calendar
 from xclim.core.missing import at_least_n_valid
-from xclim.core.units import declare_units, rate2amount
+from xclim.core.units import declare_units, rate2amount, to_agg_units
+from xclim.indices.generic import select_time, threshold_count
 
 from . import generic
 
 __all__ = [
     "base_flow_index",
+    "flow_index",
+    "high_flow_frequency",
+    "low_flow_frequency",
     "melt_and_precip_max",
     "rb_flashiness_index",
     "snd_max",
@@ -23,7 +27,7 @@
 
 
 @declare_units(q="[discharge]")
-def base_flow_index(q: xarray.DataArray, freq: str = "YS") -> xarray.DataArray:
+def base_flow_index(q: xr.DataArray, freq: str = "YS") -> xr.DataArray:
     r"""Base flow index.
 
     Return the base flow index, defined as the minimum 7-day average flow divided by the mean flow.
@@ -67,7 +71,7 @@ def base_flow_index(q: xarray.DataArray, freq: str = "YS") -> xarray.DataArray:
 
 
 @declare_units(q="[discharge]")
-def rb_flashiness_index(q: xarray.DataArray, freq: str = "YS") -> xarray.DataArray:
+def rb_flashiness_index(q: xr.DataArray, freq: str = "YS") -> xr.DataArray:
     r"""Richards-Baker flashiness index.
 
     Measures oscillations in flow relative to total flow, quantifying the frequency and rapidity of short term changes
@@ -105,7 +109,7 @@ def rb_flashiness_index(q: xarray.DataArray, freq: str = "YS") -> xarray.DataArr
 
 
 @declare_units(snd="[length]")
-def snd_max(snd: xarray.DataArray, freq: str = "YS-JUL") -> xarray.DataArray:
+def snd_max(snd: xr.DataArray, freq: str = "YS-JUL") -> xr.DataArray:
     """Maximum snow depth.
 
     The maximum daily snow depth.
@@ -126,7 +130,7 @@ def snd_max(snd: xarray.DataArray, freq: str = "YS-JUL") -> xarray.DataArray:
 
 
 @declare_units(snd="[length]")
-def snd_max_doy(snd: xarray.DataArray, freq: str = "YS-JUL") -> xarray.DataArray:
+def snd_max_doy(snd: xr.DataArray, freq: str = "YS-JUL") -> xr.DataArray:
     """Maximum snow depth day of year.
 
     Day of year when surface snow reaches its peak value. If snow depth is 0 over entire period, return NaN.
@@ -157,7 +161,7 @@ def snd_max_doy(snd: xarray.DataArray, freq: str = "YS-JUL") -> xarray.DataArray
 
 
 @declare_units(snw="[mass]/[area]")
-def snw_max(snw: xarray.DataArray, freq: str = "YS-JUL") -> xarray.DataArray:
+def snw_max(snw: xr.DataArray, freq: str = "YS-JUL") -> xr.DataArray:
     """Maximum snow amount.
 
     The maximum daily snow amount.
@@ -178,7 +182,7 @@ def snw_max(snw: xarray.DataArray, freq: str = "YS-JUL") -> xarray.DataArray:
 
 
 @declare_units(snw="[mass]/[area]")
-def snw_max_doy(snw: xarray.DataArray, freq: str = "YS-JUL") -> xarray.DataArray:
+def snw_max_doy(snw: xr.DataArray, freq: str = "YS-JUL") -> xr.DataArray:
     """Maximum snow amount day of year.
 
     Day of year when surface snow amount reaches its peak value. If snow amount is 0 over entire period, return NaN.
@@ -210,8 +214,8 @@ def snw_max_doy(snw: xarray.DataArray, freq: str = "YS-JUL") -> xarray.DataArray
 
 @declare_units(snw="[mass]/[area]")
 def snow_melt_we_max(
-    snw: xarray.DataArray, window: int = 3, freq: str = "YS-JUL"
-) -> xarray.DataArray:
+    snw: xr.DataArray, window: int = 3, freq: str = "YS-JUL"
+) -> xr.DataArray:
     """Maximum snow melt.
 
     The maximum snow melt over a given number of days expressed in snow water equivalent.
@@ -244,8 +248,8 @@ def snow_melt_we_max(
 
 @declare_units(snw="[mass]/[area]", pr="[precipitation]")
 def melt_and_precip_max(
-    snw: xarray.DataArray, pr: xarray.DataArray, window: int = 3, freq: str = "YS-JUL"
-) -> xarray.DataArray:
+    snw: xr.DataArray, pr: xr.DataArray, window: int = 3, freq: str = "YS-JUL"
+) -> xr.DataArray:
     """Maximum snow melt and precipitation.
 
     The maximum snow melt plus precipitation over a given number of days expressed in snow water equivalent.
@@ -279,3 +283,106 @@ def melt_and_precip_max(
     out = agg.resample(time=freq).max(dim="time")
     out.attrs["units"] = snw.units
     return out
+
+
+@declare_units(q="[discharge]")
+def flow_index(q: xr.DataArray, p: float = 0.95) -> xr.DataArray:
+    """
+    Flow index
+
+    Calculate the pth percentile of daily streamflow normalized by the median flow.
+
+    Parameters
+    ----------
+    q : xr.DataArray
+        Daily streamflow data.
+    p : float
+        Percentile for calculating the flow index, between 0 and 1. Default of 0.95 is for high flows.
+
+    Returns
+    -------
+    xr.DataArray
+        Normalized Qp, which is the p th percentile of daily streamflow normalized by the median flow.
+
+    References
+    ----------
+    :cite:cts:`addor2018,Clausen2000`
+    """
+    qp = q.quantile(p, dim="time")
+    q_median = q.median(dim="time")
+    out = qp / q_median
+    out.attrs["units"] = "1"
+    return out
+
+
+@declare_units(q="[discharge]")
+def high_flow_frequency(
+    q: xr.DataArray, threshold_factor: int = 9, freq: str = "YS-OCT", **indexer
+) -> xr.DataArray:
+    """
+    High flow frequency.
+
+    Calculate the number of days in a given period with flows greater than a specified threshold. By default, the
+    period is the water year starting on 1st October and ending on 30th September, as commonly defined in North America.
+
+    Parameters
+    ----------
+    q : xr.DataArray
+        Daily streamflow data.
+    threshold_factor : int
+        Factor by which the median flow is multiplied to set the high flow threshold, default is 9.
+    freq : str, optional
+        Resampling frequency, default is 'YS-OCT' for water year starting in October and ending in September.
+    indexer
+        Indexing parameters to perform a temporal subset of the data.
+        It accepts the same arguments as :py:func:`xclim.indices.generic.select_time`.
+
+    Returns
+    -------
+    xr.DataArray
+        Calculated mean of high flow days per water year
+
+    References
+    ----------
+    :cite:cts:`addor2018,Clausen2000`
+    """
+    median_flow = q.median(dim="time")
+    threshold = threshold_factor * median_flow
+    sel = select_time(q, **indexer)
+    out = threshold_count(sel, ">", threshold, freq=freq)
+    return to_agg_units(out, q, "count")
+
+
+@declare_units(q="[discharge]")
+def low_flow_frequency(
+    q: xr.DataArray, threshold_factor: float = 0.2, freq: str = "YS-OCT", **indexer
+) -> xr.DataArray:
+    """
+    Calculate the number of days in a given period with flows lower than a specified threshold. By default, the period is the water year starting on 1st October and ending on 30th September, as commonly defined in North America.
+
+    Parameters
+    ----------
+    q : xr.DataArray
+        Daily streamflow data.
+    threshold_factor : float
+        Factor by which the mean flow is multiplied to set the low flow threshold, default is 0.2.
+    freq : str, optional
+        Resampling frequency, default is 'YS-OCT' for water year starting in October and ending in September.
+    indexer
+        Indexing parameters to perform a temporal subset of the data.
+        It accepts the same arguments as :py:func:`xclim.indices.generic.select_time`.
+
+    Returns
+    -------
+    xr.DataArray
+        Calculated mean of low flow days per water year.
+
+    References
+    ----------
+    :cite:cts:`Olden2003`
+    """
+    mean_flow = q.mean(dim="time")
+    threshold = threshold_factor * mean_flow
+    sel = select_time(q, **indexer)
+    out = threshold_count(q, "<", threshold, freq=freq)
+    return to_agg_units(out, q, "count")