diff --git a/.gitignore b/.gitignore index cef52ef05..a9e772152 100644 --- a/.gitignore +++ b/.gitignore @@ -122,3 +122,6 @@ docs/variables.json # dask dask-worker-space + +# Apple +.DS_Store diff --git a/AUTHORS.rst b/AUTHORS.rst index 8fbf8baad..14cb18707 100644 --- a/AUTHORS.rst +++ b/AUTHORS.rst @@ -46,4 +46,5 @@ Contributors * Javier Diez-Sierra `@JavierDiezSierra `_ * Hui-Min Wang `@Hem-W `_ * Adrien Lamarche `@LamAdr `_ +* Faisal Mahmood `@faimahsho `_ * Sebastian Lehner `@seblehner `_ diff --git a/docs/references.bib b/docs/references.bib index 52a89e295..bc4ed54fb 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -2153,7 +2153,39 @@ @article{droogers2002 type = {Article} } -@Article{robin_2019, +@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} +} + +@article{robin_2019, author = {Robin, Y. and Vrac, M. and Naveau, P. and Yiou, P.}, title = {Multivariate stochastic bias corrections with optimal transport}, journal = {Hydrology and Earth System Sciences}, diff --git a/tests/test_hydrology.py b/tests/test_hydrology.py index 0a967a42a..0572d12e1 100644 --- a/tests/test_hydrology.py +++ b/tests/test_hydrology.py @@ -70,3 +70,33 @@ 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 * 2) * 10 + a[10:50] = 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 * 2) + 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, 0]) + + +class TestLowflowfrequency: + def test_simple(self, q_series): + a = np.ones(365 * 2) * 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, 0]) diff --git a/tests/test_land.py b/tests/test_land.py index a1782f1be..aedc1563b 100644 --- a/tests/test_land.py +++ b/tests/test_land.py @@ -66,3 +66,34 @@ def test_snw_storm_days(snw_series): snw = snw_series(a) out = land.snw_storm_days(snw, thresh="0.5 kg m-2") np.testing.assert_array_equal(out, [9, np.nan]) + + +def test_flow_index(q_series): + a = np.ones(365 * 2) * 10 + a[10:50] = 50 + q = q_series(a) + + out = land.flow_index(q, p=0.95) + np.testing.assert_array_equal(out, 5) + + +def test_high_flow_frequency(q_series): + a = np.zeros(366 * 2) * 10 + a[50:60] = 10 + a[200:210] = 20 + q = q_series(a) + out = land.high_flow_frequency( + q, + threshold_factor=9, + freq="YS", + ) + np.testing.assert_array_equal(out, [20, 0, np.nan]) + + +def test_low_flow_frequency(q_series): + a = np.ones(366 * 2) * 10 + a[50:60] = 1 + a[200:210] = 1 + q = q_series(a) + out = land.low_flow_frequency(q, threshold_factor=0.2, freq="YS") + np.testing.assert_array_equal(out, [20, 0, np.nan]) diff --git a/xclim/data/fr.json b/xclim/data/fr.json index 24cf47b42..64bae13e1 100644 --- a/xclim/data/fr.json +++ b/xclim/data/fr.json @@ -1392,6 +1392,24 @@ "title": "Rayonnement de grandes longueurs d'onde ascendant", "abstract": "Calcul du rayonnement de grandes longueurs d'onde ascendant à partir du rayonnement de grandes longueurs d'onde net et descendant." }, + "FLOW_INDEX": { + "long_name": "Indice de débit", + "description": "Le {p}e percentile du débit normalisé par le débit médian.", + "title": "Indice de débit", + "abstract": "Indice de débit, calculé comme le rapport entre un percentile donné du débit et la médiane." + }, + "HIGH_FLOW_FREQUENCY": { + "long_name": "Fréquence des débits élevés", + "description": "Nombre de jours où le débit est {threshold_factor} fois supérieur à la médiane.", + "title": "Fréquence des débits élevés", + "abstract": "Fréquence des débits élevés, calculée comme le nombre de jours où le débit est supérieur à un multiple de la médiane." + }, + "LOW_FLOW_FREQUENCY": { + "long_name": "Fréquence des débits faibles", + "description": "Nombre de jours où le débit est inférieur à une fraction ({threshold_factor}) de la moyenne.", + "title": "Fréquence des débits faibles", + "abstract": "Fréquence des débits faibles, calculée comme le nombre de jours où le débit est inférieur à une fraction de la moyenne." + }, "CP": { "title": "Portions de froid (Chill Portions) selon le modèle Dynamique.", "abstract": "Les portions de froid (chill portions) sont une mesure du potentiel de débourrement de différentes récoltes basée sur le modèle Dynamique.", diff --git a/xclim/indicators/land/_streamflow.py b/xclim/indicators/land/_streamflow.py index 63a9ba39f..fc73699af 100644 --- a/xclim/indicators/land/_streamflow.py +++ b/xclim/indicators/land/_streamflow.py @@ -3,14 +3,27 @@ from __future__ import annotations from xclim.core.cfchecks import check_valid -from xclim.core.indicator import ResamplingIndicator +from xclim.core.indicator import ( + ReducingIndicator, + ResamplingIndicator, +) from xclim.core.units import declare_units -from xclim.indices import base_flow_index, generic, rb_flashiness_index +from xclim.indices import ( + base_flow_index, + flow_index, + generic, + high_flow_frequency, + low_flow_frequency, + rb_flashiness_index, +) __all__ = [ "base_flow_index", "doy_qmax", "doy_qmin", + "flow_index", + "high_flow_frequency", + "low_flow_frequency", "rb_flashiness_index", ] @@ -74,3 +87,37 @@ def cfcheck(q): compute=declare_units(da="[discharge]")(generic.select_resample_op), parameters={"op": generic.doymin, "out_units": None}, ) + +flow_index = ReducingIndicator( + realm="land", + context="hydro", + title="Flow index", + identifier="flow_index", + var_name="q_flow_index", + long_name="Flow index", + description="{p}th percentile normalized by the median flow.", + units="1", + compute=flow_index, +) + + +high_flow_frequency = Streamflow( + title="High flow frequency", + identifier="high_flow_frequency", + var_name="q_high_flow_frequency", + long_name="High flow frequency", + description="{freq} frequency of flows greater than {threshold_factor} times the median flow.", + units="days", + compute=high_flow_frequency, +) + + +low_flow_frequency = Streamflow( + title="Low flow frequency", + identifier="low_flow_frequency", + var_name="q_low_flow_frequency", + long_name="Low flow frequency", + description="{freq} frequency of flows smaller than a fraction ({threshold_factor}) of the mean flow.", + units="days", + compute=low_flow_frequency, +) diff --git a/xclim/indices/_hydrology.py b/xclim/indices/_hydrology.py index 3785157af..23d320cbc 100644 --- a/xclim/indices/_hydrology.py +++ b/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 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,103 @@ 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:`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" +) -> xr.DataArray: + """ + High flow frequency. + + Calculate the number of days in a given period with flows greater than a specified threshold, given as a + multiple of the median flow. 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. + + Returns + ------- + xr.DataArray + Number of high flow days. + + References + ---------- + :cite:cts:`addor2018,Clausen2000` + """ + median_flow = q.median(dim="time") + threshold = threshold_factor * median_flow + out = threshold_count(q, ">", 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" +) -> xr.DataArray: + """ + Low flow frequency. + + Calculate the number of days in a given period with flows lower than a specified threshold, given by a fraction + of the mean flow. 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. + + Returns + ------- + xr.DataArray + Number of low flow days. + + References + ---------- + :cite:cts:`Olden2003` + """ + mean_flow = q.mean(dim="time") + threshold = threshold_factor * mean_flow + out = threshold_count(q, "<", threshold, freq=freq) + return to_agg_units(out, q, "count")