Add chill portion and chill units

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`).
+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`).
 
 Announcements
 ^^^^^^^^^^^^^
@@ -14,11 +14,13 @@ 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 ``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 features and enhancements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 * New generic ``xclim.indices.generic.spell_mask``  that returns a mask of which days are part of a spell. Supports multivariate conditions and weights. Used in new generic index ``xclim.indices.generic.bivariate_spell_length_statistics`` that extends ``spell_length_statistics`` to two variables.  (:pull:`1885`).
 * Indicator parameters can now be assigned a new name, different from the argument name in the compute function. (:pull:`1885`).
+* Helper function ``xclim.indices.helpers.make_hourly_temperature`` to estimate hourly temperatures from daily min and max temperatures. (:pull:`1909`).
 
 Bug fixes
 ^^^^^^^^^

tests/test_atmos.py

@@ -6,6 +6,7 @@
 import xarray as xr
 
 from xclim import atmos, set_options
+from xclim.indices.helpers import make_hourly_temperature
 
 K2C = 273.16
 
@@ -624,3 +625,37 @@ def test_late_frost_days(self, atmosds):
         lfd = atmos.late_frost_days(tasmin, date_bounds=("04-01", "06-30"))
 
         np.testing.assert_allclose(lfd.isel(time=0), exp, rtol=1e-03)
+
+
+def test_chill_units(atmosds):
+    tasmax = atmosds.tasmax
+    tasmin = atmosds.tasmin
+    tas = make_hourly_temperature(tasmin, tasmax)
+    cu = atmos.chill_units(tas, date_bounds=("04-01", "06-30"))
+    assert cu.attrs["units"] == "1"
+    assert cu.name == "cu"
+    assert cu.time.size == 4
+
+    exp = [-5029.5, -6634.5, -5993.0, -6596.0, -5654.0]
+    np.testing.assert_allclose(cu.isel(time=0), exp, rtol=1e-03)
+
+
+def test_chill_portions(atmosds):
+    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")
+    assert cp.attrs["units"] == "1"
+    assert cp.name == "cp"
+    # Although its 4 years of data its 5 seasons starting in July
+    assert cp.time.size == 5
+
+    # First timestamp is all nan so we test on second
+    exp = [
+        99.91534493319304,
+        96.84276565115084,
+        16.720155001987106,
+        77.86606372811245,
+        140.85516682327875,
+    ]
+    np.testing.assert_allclose(cp.isel(time=1), exp, rtol=1e-03)

tests/test_helpers.py

@@ -132,3 +132,45 @@ def test_cosine_of_solar_zenith_angle():
         ]
     )
     np.testing.assert_allclose(cza[:4, :], exp_cza, rtol=1e-3)
+
+
+@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)
+    tasmin = tasmin_series(np.array([0]), units="degC", cftime=cftime).expand_dims(
+        lat=[0]
+    )
+
+    tasmin.lat.attrs["units"] = "degree_north"
+    tas_hourly = helpers.make_hourly_temperature(tasmax, tasmin)
+    assert tas_hourly.attrs["units"] == "degC"
+    assert tas_hourly.time.size == 24
+    expected = np.array(
+        [
+            0.0,
+            3.90180644,
+            7.65366865,
+            11.11140466,
+            14.14213562,
+            16.62939225,
+            18.47759065,
+            19.61570561,
+            20.0,
+            19.61570561,
+            18.47759065,
+            16.62939225,
+            14.14213562,
+            10.32039099,
+            8.0848137,
+            6.49864636,
+            5.26831939,
+            4.26306907,
+            3.41314202,
+            2.67690173,
+            2.02749177,
+            1.44657476,
+            0.92107141,
+            0.44132444,
+        ]
+    )
+    np.testing.assert_allclose(tas_hourly.isel(lat=0).values, expected)

tests/test_indices.py

@@ -346,6 +346,32 @@ def test_bedd(self, method, end_date, deg_days, max_deg_days):
             if method == "icclim":
                 np.testing.assert_array_equal(bedd, bedd_high_lat)
 
+    def test_chill_portions(self, tas_series):
+        tas = tas_series(np.linspace(0, 15, 120 * 24) + K2C, freq="h")
+        out = xci.chill_portions(tas)
+        assert out[0] == 72.24417644977083
+
+    def test_chill_units(self, tas_series):
+        num_cu_0 = 10
+        num_cu_1 = 20
+        num_cu_05 = 15
+        num_cu_min_05 = 10
+        num_cu_min_1 = 5
+
+        tas = tas_series(
+            np.array(
+                num_cu_0 * [1.1]
+                + num_cu_05 * [2.0]
+                + num_cu_1 * [5.6]
+                + num_cu_min_05 * [16.0]
+                + num_cu_min_1 * [20.0]
+            )
+            + K2C,
+            freq="h",
+        )
+        out = xci.chill_units(tas)
+        assert out[0] == 0.5 * num_cu_05 + num_cu_1 - 0.5 * num_cu_min_05 - num_cu_min_1
+
     def test_cool_night_index(self, open_dataset):
         ds = open_dataset("cmip5/tas_Amon_CanESM2_rcp85_r1i1p1_200701-200712.nc")
         ds = ds.rename(dict(tas="tasmin"))

xclim/indicators/atmos/_temperature.py

@@ -4,12 +4,19 @@
 
 from xclim import indices
 from xclim.core import cfchecks
-from xclim.core.indicator import Daily, Indicator, ResamplingIndicatorWithIndexing
+from xclim.core.indicator import (
+    Daily,
+    Hourly,
+    Indicator,
+    ResamplingIndicatorWithIndexing,
+)
 from xclim.core.utils import InputKind
 
 __all__ = [
     "australian_hardiness_zones",
     "biologically_effective_degree_days",
+    "chill_portions",
+    "chill_units",
     "cold_spell_days",
     "cold_spell_duration_index",
     "cold_spell_frequency",
@@ -100,13 +107,26 @@ class Temp(Daily):
     keywords = "temperature"
 
 
+class TempHourly(Hourly):
+    """Temperature indicators involving hourly temperature."""
+
+    keywords = "temperature"
+
+
 class TempWithIndexing(ResamplingIndicatorWithIndexing):
     """Indicators involving daily temperature and adding an indexing possibility."""
 
     src_freq = "D"
     keywords = "temperature"
 
 
+class TempHourlyWithIndexing(ResamplingIndicatorWithIndexing):
+    """Indicators involving hourly temperature and adding an indexing possibility."""
+
+    src_freq = "h"
+    keywords = "temperature"
+
+
 tn_days_above = TempWithIndexing(
     title="Number of days with minimum temperature above a given threshold",
     identifier="tn_days_above",
@@ -1445,3 +1465,38 @@ def cfcheck(self, tas, snd=None):
     compute=indices.hardiness_zones,
     parameters={"method": "usda"},
 )
+
+chill_portions = TempHourly(
+    title="Chill portions",
+    identifier="cp",
+    units="",
+    cell_methods="time: sum",
+    description="Chill portions are a measure to estimate the bud breaking potential of different crops. "
+    "The constants and functions are taken from Luedeling et al. (2009) which formalises "
+    "the method described in Fishman et al. (1987). ",
+    abstract="Chill portions are a measure to estimate the bud breaking potential of different crops. "
+    "The constants and functions are taken from Luedeling et al. (2009) which formalises "
+    "the method described in Fishman et al. (1987). "
+    "The model computes the accumulation of cold temperatures in a two-step process. "
+    "First, cold temperatures contribute to an intermediate product that is transformed to a chill portion "
+    "once it exceeds a certain concentration. The intermediate product can be broken down at higher temperatures "
+    "but the final product is stable even at higher temperature. "
+    "Thus the dynamic model is more accurate than other chill models like the Chilling hours or Utah model, "
+    "especially in moderate climates like Israel, California or Spain.",
+    long_name="Chill portions after the Dynamic Model",
+    allowed_periods=["Y"],
+    compute=indices.chill_portions,
+)
+
+chill_units = TempHourlyWithIndexing(
+    title="Chill units",
+    identifier="cu",
+    units="",
+    cell_methods="time: sum",
+    description="Chill units are a measure to estimate the bud breaking potential of different crops based on the Utah model developed in "
+    "Richardson et al. (1974). The Utah model assigns a weight to each hour depending on the temperature recognising that high temperatures can "
+    "actually decrease the potential for bud breaking.",
+    long_name="Chill units after the Utah Model",
+    allowed_periods=["Y"],
+    compute=indices.chill_units,
+)

xclim/indices/_agro.py

@@ -35,6 +35,8 @@
 
 __all__ = [
     "biologically_effective_degree_days",
+    "chill_portions",
+    "chill_units",
     "cool_night_index",
     "corn_heat_units",
     "dryness_index",
@@ -1524,3 +1526,120 @@ def hardiness_zones(
 
     zones = zones.assign_attrs(units="")
     return zones
+
+
+def _accumulate_intermediate(prev_E, prev_xi, curr_xs, curr_ak1):
+    """Accumulate the intermediate product based on the previous concentration and the current temperature."""
+    curr_S = np.where(prev_E < 1, prev_E, prev_E - prev_E * prev_xi)
+    return curr_xs - (curr_xs - curr_S) * np.exp(-curr_ak1)
+
+
+def _chill_portion_one_season(tas_K):
+    """Computes the chill portion for a single season based on the dynamic model on a numpy array."""
+    # Constants as described in Luedeling et al. (2009)
+    E0 = 4153.5
+    E1 = 12888.8
+    A0 = 139500
+    A1 = 2.567e18
+    SLP = 1.6
+    TETMLT = 277
+    AA = A0 / A1
+    EE = E1 - E0
+
+    ftmprt = SLP * TETMLT * (tas_K - TETMLT) / tas_K
+    sr = np.exp(ftmprt)
+    xi = sr / (1 + sr)
+    xs = AA * np.exp(EE / tas_K)
+    ak1 = A1 * np.exp(-E1 / tas_K)
+
+    inter_E = np.zeros_like(tas_K)
+    for i in range(1, tas_K.shape[-1]):
+        inter_E[..., i] = _accumulate_intermediate(
+            inter_E[..., i - 1], xi[..., i - 1], xs[..., i], ak1[..., i]
+        )
+    delta = np.where(inter_E >= 1, inter_E * xi, 0)
+
+    return delta
+
+
+def _apply_chill_portion_one_season(tas_K):
+    """Apply the chill portion function on to an xarray DataArray."""
+    tas_K = tas_K.chunk(time=-1)
+    return xarray.apply_ufunc(
+        _chill_portion_one_season,
+        tas_K,
+        input_core_dims=[["time"]],
+        output_core_dims=[["time"]],
+        output_dtypes=[tas_K.dtype],
+        dask="parallelized",
+    ).sum("time")
+
+
+@declare_units(tas="[temperature]")
+def chill_portions(
+    tas: xarray.DataArray, freq: str = "YS", **indexer
+) -> xarray.DataArray:
+    """Chill portion based on the dynamic model
+
+    Chill portions are a measure to estimate the bud breaking potential of different crop.
+    The constants and functions are taken from Luedeling et al. (2009) which formalises
+    the method described in Fishman et al. (1987). The model computes the accumulation of
+    an intermediate product that is transformed to the final product once it exceeds a
+    certain concentration. The intermediate product can be broken down at higher temperatures
+    but the final product is stable even at higher temperature. Thus the dynamic model is
+    more accurate than the Utah model especially in moderate climates like Israel,
+    California or Spain.
+
+    Parameters
+    ----------
+    tas : xr.DataArray
+        Hourly temperature.
+
+    Returns
+    -------
+    xr.DataArray, [unitless]
+        Chill portions after the Dynamic Model
+    """
+    tas_K: xarray.DataArray = select_time(
+        convert_units_to(tas, "K"), drop=True, **indexer
+    )
+    return (
+        tas_K.resample(time=freq)
+        .map(_apply_chill_portion_one_season)
+        .assign_attrs(units="")
+    )
+
+
+@declare_units(tas="[temperature]")
+def chill_units(tas: xarray.DataArray, freq: str = "YS") -> xarray.DataArray:
+    """Chill units using the Utah model
+
+    Chill units are a measure to estimate the bud breaking potential of different crop based on Richardson et al. (1974).
+    The Utah model assigns a weight to each hour depending on the temperature recognising that high temperatures can actual decrease,
+    the potential for bud breaking.
+
+    Parameters
+    ----------
+    tas : xr.DataArray
+        Hourly temperature.
+
+    Returns
+    -------
+    xr.DataArray, [unitless]
+        Chill units using the Utah model
+    """
+    tas = convert_units_to(tas, "degC")
+    cu = xarray.where(
+        (tas <= 1.4) | ((tas > 12.4) & (tas <= 15.9)),
+        0,
+        xarray.where(
+            ((tas > 1.4) & (tas <= 2.4)) | ((tas > 9.1) & (tas <= 12.4)),
+            0.5,
+            xarray.where(
+                (tas > 2.4) & (tas <= 9.1),
+                1,
+                xarray.where((tas > 15.9) & (tas <= 17.9), -0.5, -1),
+            ),
+        ),
+    )
+    return cu.resample(time=freq).sum().assign_attrs(units="")