Merge branch 'main' into fix-1822

CHANGELOG.rst

@@ -4,7 +4,7 @@ Changelog
 
 v0.53.0 (unreleased)
 --------------------
-Contributors to this version: Adrien Lamarche (:user:`LamAdr`), Trevor James Smith (:user:`Zeitsperre`),  Éric Dupuis (:user:`coxipi`), Pascal Bourgault (:user:`aulemahal`), Sascha Hofmann (:user:`saschahofmann`).
+Contributors to this version: Adrien Lamarche (:user:`LamAdr`), Trevor James Smith (:user:`Zeitsperre`),  Éric Dupuis (:user:`coxipi`), Pascal Bourgault (:user:`aulemahal`), Sascha Hofmann (:user:`saschahofmann`), David Huard (:user:`huard`).
 
 Announcements
 ^^^^^^^^^^^^^
@@ -13,9 +13,10 @@ Announcements
 
 New indicators
 ^^^^^^^^^^^^^^
-* New ``heat_spell_frequency``, ``heat_spell_max_length`` and ``heat_spell_total_length`` : spell length statistics on a bivariate condition that uses the average over a window by default. (:pull:`1885`).
+* New ``heat_spell_frequency``, ``heat_spell_max_length`` and ``heat_spell_total_length`` : spell length statistics on a bivariate condition that uses the average over a window by default. (:pull:`1885`, :pull:`1778`).
 * New ``hot_spell_max_magnitude`` : yields the magnitude of the most intensive heat wave. (:pull:`1926`).
 * New ``chill_portion`` and ``chill_unit``: chill portion based on the Dynamic Model and chill unit based on the Utah model indicators. (:issue:`1753`, :pull:`1909`).
+* New ``water_cycle_intensity``: yields the sum of precipitation and actual evapotranspiration. (:issue:`410`, :pull:`1947`).
 
 New features and enhancements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -24,6 +25,10 @@ New features and enhancements
 * Add attribute ``units_metadata`` to outputs representing a difference between temperatures. This is needed to disambiguate temperature differences from absolute temperature. Changes affect indicators ``daily_temperature_range``, ``daily_temperature_range_variability``, ``extreme_temperature_range``, ``interday_diurnal_temperature_range``, and all degree-day indicators. Implemented using a new ``pint2cfattrs`` function to convert pint units to a dictionary of CF attributes. ``units2pint`` is also modified to support ``units_metadata`` attributes in DataArrays. Some SDBA properties and measures previously returning units of ``delta_degC`` will now return the original input DataArray units accompanied with the ``units_metadata`` attribute. (:issue:`1822`, :pull:`1830`).
 * ``xclim.indices.run_length.windowed_max_run_sum`` accumulates positive values across runs and yields the the maximum valued run. (:pull:`1926`).
 * Helper function ``xclim.indices.helpers.make_hourly_temperature`` to estimate hourly temperatures from daily min and max temperatures. (:pull:`1909`).
+* New global option ``resample_map_blocks`` to wrap all ``resample().map()`` code inside a ``xr.map_blocks`` to lower the number of dask tasks. Uses utility ``xclim.indices.helpers.resample_map`` and requires ``flox`` to ensure the chunking allows such block-mapping. Defaults to False. (:pull:`1848`).
+* ``xclim.indices.run_length.runs_with_holes`` allows to input a condition that must be met for a run to start and a second condition that must be met for the run to stop. (:pull:`1778`).
+* New generic compute function ``xclim.indices.generic.thresholded_events`` that finds events based on a threshold condition and returns basic stats for each. See also ``xclim.indices.run_length.find_events``. (:pull:`1778`).
+* ``xclim.core.units.rate2amount`` and ``xclim.core.units.amount2rate`` can now also accept quantities (pint objects or strings), in which case the ``dim`` argument must be the ``time`` coordinate through which we can find the sampling rate. (:pull:`1778`).
 
 Bug fixes
 ^^^^^^^^^

docs/references.bib

@@ -2307,3 +2307,16 @@ @article{zhang_high_2022
   urldate = {2024-10-03},
   langid = {english}
 }
+
+@article{huntington_2018,
+  title = {A new indicator framework for quantifying the intensity of the terrestrial water cycle},
+  journal = {Journal of Hydrology},
+  volume = {559},
+  pages = {361-372},
+  year = {2018},
+  issn = {0022-1694},
+  doi = {https://doi.org/10.1016/j.jhydrol.2018.02.048},
+  url = {https://www.sciencedirect.com/science/article/pii/S0022169418301276},
+  author = {Thomas G. Huntington and Peter K. Weiskel and David M. Wolock and Gregory J. McCabe},
+  keywords = {Intensification of the water cycle, Water cycle intensity, Water-balance model, Evapotranspiration, Aridification},
+}

environment.yml

@@ -11,6 +11,7 @@ dependencies:
   - click >=8.1
   - dask >=2.6.0
   - filelock >=3.14.0
+  - flox >= 0.9
   - jsonpickle >=3.1.0
   - numba >=0.54.1
   - numpy >=1.23.0

pyproject.toml

@@ -111,7 +111,7 @@ docs = [
   "sphinxcontrib-bibtex",
   "sphinxcontrib-svg2pdfconverter[Cairosvg]"
 ]
-extras = ["fastnanquantile >=0.0.2", "POT >=0.9.4"]
+extras = ["fastnanquantile >=0.0.2", "flox >=0.9", "POT >=0.9.4"]
 all = ["xclim[dev]", "xclim[docs]", "xclim[extras]"]
 
 [project.scripts]
@@ -135,7 +135,7 @@ target-version = [
 ]
 
 [tool.bumpversion]
-current_version = "0.52.3-dev.8"
+current_version = "0.52.3-dev.11"
 commit = true
 commit_args = "--no-verify"
 tag = false

tests/conftest.py

@@ -81,6 +81,12 @@ def pr_series():
     return _pr_series
 
 
+@pytest.fixture
+def evspsbl_series():
+    """Return precipitation time series."""
+    return partial(test_timeseries, variable="evspsbl")
+
+
 @pytest.fixture
 def prc_series():
     """Return convective precipitation time series."""

tests/test_atmos.py

@@ -3,6 +3,7 @@
 from __future__ import annotations
 
 import numpy as np
+import pytest
 import xarray as xr
 
 from xclim import atmos, set_options
@@ -642,11 +643,18 @@ def test_chill_units(atmosds):
     np.testing.assert_allclose(cu.isel(location=0), exp, rtol=1e-03)
 
 
-def test_chill_portions(atmosds):
+@pytest.mark.parametrize("use_dask", [True, False])
+def test_chill_portions(atmosds, use_dask):
+    pytest.importorskip("flox")
     tasmax = atmosds.tasmax
     tasmin = atmosds.tasmin
     tas = make_hourly_temperature(tasmin, tasmax)
-    cp = atmos.chill_portions(tas, date_bounds=("09-01", "03-30"), freq="YS-JUL")
+    if use_dask:
+        tas = tas.chunk(time=tas.time.size // 2, location=1)
+
+    with set_options(resample_map_blocks=True):
+        cp = atmos.chill_portions(tas, date_bounds=("09-01", "03-30"), freq="YS-JUL")
+
     assert cp.attrs["units"] == "1"
     assert cp.name == "cp"
     # Although its 4 years of data its 5 seasons starting in July
@@ -656,3 +664,11 @@ def test_chill_portions(atmosds):
     # due to implementation details
     exp = [np.nan, 99.91534493, 92.5473925, 99.03177047, np.nan]
     np.testing.assert_allclose(cp.isel(location=0), exp, rtol=1e-03)
+
+
+def test_water_cycle_intensity(pr_series, evspsbl_series):
+    pr = pr_series(np.ones(31))
+    evspsbl = evspsbl_series(np.ones(31))
+
+    wci = atmos.water_cycle_intensity(pr=pr, evspsbl=evspsbl, freq="MS")
+    np.testing.assert_allclose(wci, 2 * 60 * 60 * 24 * 31)

tests/test_generic.py

@@ -8,6 +8,7 @@
 
 from xclim.core.calendar import doy_to_days_since, select_time
 from xclim.indices import generic
+from xclim.testing.helpers import assert_lazy
 
 K2C = 273.15
 
@@ -774,3 +775,126 @@ def test_spell_length_statistics_multi(tasmin_series, tasmax_series):
     )
     xr.testing.assert_equal(outs, outm)
     np.testing.assert_allclose(outc, 1)
+
+
+class TestThresholdedEvents:
+
+    @pytest.mark.parametrize("use_dask", [True, False])
+    def test_simple(self, pr_series, use_dask):
+        arr = np.array([0, 0, 0, 1, 2, 3, 0, 3, 3, 10, 0, 0, 0, 0, 0, 1, 2, 2, 2, 0, 0, 0, 0, 0, 0, 1, 3, 3, 2, 0, 0, 0, 2, 0, 0, 0, 0])  # fmt: skip
+        pr = pr_series(arr, start="2000-01-01", units="mm")
+        if use_dask:
+            pr = pr.chunk(-1)
+
+        with assert_lazy:
+            out = generic.thresholded_events(
+                pr,
+                thresh="1 mm",
+                op=">=",
+                window=3,
+            )
+
+        assert out.event.size == np.ceil(arr.size / (3 + 1))
+        out = out.load().dropna("event", how="all")
+
+        np.testing.assert_array_equal(out.event_length, [3, 3, 4, 4])
+        np.testing.assert_array_equal(out.event_effective_length, [3, 3, 4, 4])
+        np.testing.assert_array_equal(out.event_sum, [6, 16, 7, 9])
+        np.testing.assert_array_equal(
+            out.event_start,
+            np.array(
+                ["2000-01-04", "2000-01-08", "2000-01-16", "2000-01-26"],
+                dtype="datetime64[ns]",
+            ),
+        )
+
+    @pytest.mark.parametrize("use_dask", [True, False])
+    def test_diff_windows(self, pr_series, use_dask):
+        arr = np.array([0, 0, 0, 1, 2, 3, 0, 3, 3, 10, 0, 0, 0, 0, 0, 1, 2, 2, 2, 0, 0, 0, 0, 0, 0, 1, 3, 3, 2, 0, 0, 0, 2, 0, 0, 0, 0])  # fmt: skip
+        pr = pr_series(arr, start="2000-01-01", units="mm")
+        if use_dask:
+            pr = pr.chunk(-1)
+
+        # different window stop
+        out = (
+            generic.thresholded_events(
+                pr, thresh="2 mm", op=">=", window=3, window_stop=4
+            )
+            .load()
+            .dropna("event", how="all")
+        )
+
+        np.testing.assert_array_equal(out.event_length, [3, 3, 7])
+        np.testing.assert_array_equal(out.event_effective_length, [3, 3, 4])
+        np.testing.assert_array_equal(out.event_sum, [16, 6, 10])
+        np.testing.assert_array_equal(
+            out.event_start,
+            np.array(
+                ["2000-01-08", "2000-01-17", "2000-01-27"], dtype="datetime64[ns]"
+            ),
+        )
+
+    @pytest.mark.parametrize("use_dask", [True, False])
+    def test_cftime(self, pr_series, use_dask):
+        arr = np.array([0, 0, 0, 1, 2, 3, 0, 3, 3, 10, 0, 0, 0, 0, 0, 1, 2, 2, 2, 0, 0, 0, 0, 0, 0, 1, 3, 3, 2, 0, 0, 0, 2, 0, 0, 0, 0])  # fmt: skip
+        pr = pr_series(arr, start="2000-01-01", units="mm").convert_calendar("noleap")
+        if use_dask:
+            pr = pr.chunk(-1)
+
+        # cftime
+        with assert_lazy:
+            out = generic.thresholded_events(
+                pr,
+                thresh="1 mm",
+                op=">=",
+                window=3,
+                window_stop=3,
+            )
+        out = out.load().dropna("event", how="all")
+
+        np.testing.assert_array_equal(out.event_length, [7, 4, 4])
+        np.testing.assert_array_equal(out.event_effective_length, [6, 4, 4])
+        np.testing.assert_array_equal(out.event_sum, [22, 7, 9])
+        exp = xr.DataArray(
+            [1, 2, 3],
+            dims=("time",),
+            coords={
+                "time": np.array(
+                    ["2000-01-04", "2000-01-16", "2000-01-26"], dtype="datetime64[ns]"
+                )
+            },
+        )
+        np.testing.assert_array_equal(
+            out.event_start, exp.convert_calendar("noleap").time
+        )
+
+    @pytest.mark.parametrize("use_dask", [True, False])
+    def test_freq(self, pr_series, use_dask):
+        jan = [0, 0, 0, 1, 2, 3, 0, 3, 3, 10, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 2, 0, 0, 0, 0, 0, 3, 2, 3, 2]  # fmt: skip
+        fev = [2, 2, 1, 0, 0, 0, 3, 3, 4, 5, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]  # fmt: skip
+        pr = pr_series(np.array(jan + fev), start="2000-01-01", units="mm")
+        if use_dask:
+            pr = pr.chunk(-1)
+
+        with assert_lazy:
+            out = generic.thresholded_events(
+                pr, thresh="1 mm", op=">=", window=3, freq="MS", window_stop=3
+            )
+        assert out.event_length.shape == (2, 6)
+        out = out.load().dropna("event", how="all")
+
+        np.testing.assert_array_equal(out.event_length, [[7, 6, 4], [3, 5, np.nan]])
+        np.testing.assert_array_equal(
+            out.event_effective_length, [[6, 6, 4], [3, 5, np.nan]]
+        )
+        np.testing.assert_array_equal(out.event_sum, [[22, 12, 10], [5, 17, np.nan]])
+        np.testing.assert_array_equal(
+            out.event_start,
+            np.array(
+                [
+                    ["2000-01-04", "2000-01-17", "2000-01-28"],
+                    ["2000-02-01", "2000-02-07", "NaT"],
+                ],
+                dtype="datetime64[ns]",
+            ),
+        )

tests/test_helpers.py

@@ -5,8 +5,11 @@
 import pytest
 import xarray as xr
 
+from xclim.core.options import set_options
 from xclim.core.units import convert_units_to
+from xclim.core.utils import uses_dask
 from xclim.indices import helpers
+from xclim.testing.helpers import assert_lazy
 
 
 @pytest.mark.parametrize("method,rtol", [("spencer", 5e3), ("simple", 1e2)])
@@ -134,6 +137,48 @@ def test_cosine_of_solar_zenith_angle():
     np.testing.assert_allclose(cza[:4, :], exp_cza, rtol=1e-3)
 
 
+def _test_function(da, op, dim):
+    return getattr(da, op)(dim)
+
+
+@pytest.mark.parametrize(
+    ["in_chunks", "exp_chunks"], [(60, 6 * (2,)), (30, 12 * (1,)), (-1, (12,))]
+)
+def test_resample_map(tas_series, in_chunks, exp_chunks):
+    pytest.importorskip("flox")
+    tas = tas_series(365 * [1]).chunk(time=in_chunks)
+    with assert_lazy:
+        out = helpers.resample_map(
+            tas, "time", "MS", lambda da: da.mean("time"), map_blocks=True
+        )
+    assert out.chunks[0] == exp_chunks
+    out.load()  # Trigger compute to see if it actually works
+
+
+def test_resample_map_dataset(tas_series, pr_series):
+    pytest.importorskip("flox")
+    tas = tas_series(3 * 365 * [1], start="2000-01-01").chunk(time=365)
+    pr = pr_series(3 * 365 * [1], start="2000-01-01").chunk(time=365)
+    ds = xr.Dataset({"pr": pr, "tas": tas})
+    with set_options(resample_map_blocks=True):
+        with assert_lazy:
+            out = helpers.resample_map(
+                ds,
+                "time",
+                "YS",
+                lambda da: da.mean("time"),
+            )
+    assert out.chunks["time"] == (1, 1, 1)
+    out.load()
+
+
+def test_resample_map_passthrough(tas_series):
+    tas = tas_series(365 * [1])
+    with assert_lazy:
+        out = helpers.resample_map(tas, "time", "MS", lambda da: da.mean("time"))
+    assert not uses_dask(out)
+
+
 @pytest.mark.parametrize("cftime", [False, True])
 def test_make_hourly_temperature(tasmax_series, tasmin_series, cftime):
     tasmax = tasmax_series(np.array([20]), units="degC", cftime=cftime)

tests/test_indices.py

@@ -1454,13 +1454,24 @@ def test_1d(self, tasmax_series, thresh, window, op, expected):
     def test_resampling_order(self, tasmax_series, resample_before_rl, expected):
         a = np.zeros(365)
         a[5:35] = 31
-        tx = tasmax_series(a + K2C)
+        tx = tasmax_series(a + K2C).chunk()
 
         hsf = xci.hot_spell_frequency(
             tx, resample_before_rl=resample_before_rl, freq="MS"
-        )
+        ).load()
         assert hsf[1] == expected
 
+    @pytest.mark.parametrize("resample_map", [True, False])
+    def test_resampling_map(self, tasmax_series, resample_map):
+        pytest.importorskip("flox")
+        a = np.zeros(365)
+        a[5:35] = 31
+        tx = tasmax_series(a + K2C).chunk()
+
+        with set_options(resample_map_blocks=resample_map):
+            hsf = xci.hot_spell_frequency(tx, resample_before_rl=True, freq="MS").load()
+        assert hsf[1] == 1
+
 
 class TestHotSpellMaxLength:
     @pytest.mark.parametrize(
@@ -1746,10 +1757,10 @@ def test_run_start_at_0(self, pr_series):
     def test_resampling_order(self, pr_series, resample_before_rl, expected):
         a = np.zeros(365) + 10
         a[5:35] = 0
-        pr = pr_series(a)
+        pr = pr_series(a).chunk()
         out = xci.maximum_consecutive_dry_days(
             pr, freq="ME", resample_before_rl=resample_before_rl
-        )
+        ).load()
         assert out[0] == expected
 
 
@@ -3986,3 +3997,12 @@ def test_simple(self, sfcWind_series):
         pb = xci.wind_power_potential(b)
 
         np.testing.assert_array_almost_equal(pa, pb, decimal=6)
+
+
+class TestWaterCycleIntensity:
+    def test_simple(self, pr_series, evspsbl_series):
+        pr = pr_series(np.ones(31))
+        evspsbl = evspsbl_series(np.ones(31))
+
+        wci = xci.water_cycle_intensity(pr=pr, evspsbl=evspsbl, freq="MS")
+        np.testing.assert_allclose(wci, 2 * 60 * 60 * 24 * 31)