Add chill portion and chill units (#1909)

Closes #1753 

### TODO:
- [x] Tests for the changes have been added (for bug fixes / features)
- [x] (If applicable) Documentation has been added / updated (for bug
fixes / features)
  - [x] Test `make_hourly_temperature`
  - [x] Test `chill_portions`
  - [x] Test `chill_units`
- [x] CHANGELOG.rst has been updated (with summary of main changes)
- [x] Link to issue (:issue:`number`) and pull request (:pull:`number`)
has been added
- [x] Check `make_hourly_temperature` is working with other calendars
than 360_day

### What kind of change does this PR introduce?
* Adds a chill_portion indicator based on the Dynamic model
* Adds a chill_unit indicator based on the Utah model
* Adds a helper function to estimate hourly temperatures from daily
tasmin and tasmax

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
 ^^^^^^^^^

docs/references.bib

@@ -2203,3 +2203,45 @@ @article{knol_1989
 	publisher = "Springer",
 	number = "1",
 }
+
+@article{luedeling_chill_2009,
+	title = {Climate change impacts on winter chill for temperate fruit and nut production: A review},
+	journal = {Scientia Horticulturae},
+	volume = {144},
+	pages = {218-229},
+	year = {2012},
+	issn = {0304-4238},
+	doi = {https://doi.org/10.1016/j.scienta.2012.07.011},
+	url = {https://www.sciencedirect.com/science/article/pii/S0304423812003305},
+	author = {Eike Luedeling},
+	keywords = {Adaptation, Chilling Hours, Chill Portions, Climate analogues, Dynamic Model, Tree dormancy},
+	abstract = {Temperate fruit and nut species require exposure to chilling conditions in winter to break dormancy and produce high yields. Adequate winter chill is an important site characteristic for commercial orchard operations, and quantifying chill is crucial for orchard management. Climate change may impact winter chill. With a view to adapting orchards to climate change, this review assesses the state of knowledge in modelling winter chill and the performance of various modelling approaches. It then goes on to present assessments of past and projected future changes in winter chill for fruit growing regions and discusses potential adaptation strategies. Some of the most common approaches to modelling chill, in particular the Chilling Hours approach, are very sensitive to temperature increases, and have also been found to perform poorly, especially in warm growing regions. The Dynamic Model offers a more complex but also more accurate alternative, and use of this model is recommended. Chill changes projected with the Dynamic Model are typically much less severe than those estimated with other models. Nevertheless, projections of future chill consistently indicate substantial losses for the warmest growing regions, while temperate regions will experience relatively little change, and cold regions may even see chill increases. Growers can adapt to lower chill by introducing low-chill cultivars, by influencing orchard microclimates and by applying rest-breaking chemicals. Given substantial knowledge gaps in tree dormancy, accurate models are still a long way off. Since timely adaptation is essential for growers of long-lived high-value perennials, alternative ways of adaptation planning are needed. Climate analogues, which are present-day manifestations of future projected climates, can be used for identifying and testing future-adapted species and cultivars. Horticultural researchers and practitioners should work towards the development and widespread adoption of better chill accumulation and dormancy models, for facilitating quantitatively appropriate adaptation planning.}
+}
+
+@article{fishman_chill_1987,
+	title = {The temperature dependence of dormancy breaking in plants: Mathematical analysis of a two-step model involving a cooperative transition},
+	journal = {Journal of Theoretical Biology},
+	volume = {124},
+	number = {4},
+	pages = {473-483},
+	year = {1987},
+	issn = {0022-5193},
+	doi = {https://doi.org/10.1016/S0022-5193(87)80221-7},
+	url = {https://www.sciencedirect.com/science/article/pii/S0022519387802217},
+	author = {Svetlana Fishman and A. Erez and G.A. Couvillon},
+	abstract = {A two-step model describing the thermal dependence of the dormancy breaking phenomenon is developed. The model assumes that the level of dormancy completion is proportional to the amount of a certain dormancy breaking factor which accumulates in plants by a two-step process. The first step represents a reversible process of formation of a precursor for the dormancy breaking factor at low temperatures and its destruction at high temperatures. The rate constants of this process are assumed to be dependent upon the temperature according to the Arrhenius law. The second step is an irreversible cooperative transition from the unstable precursor to a stable dormancy breaking factor. The transition is assumed to occur when a critical level of the precursor is accumulated. The two-step scheme is analysed mathematically. This model explains qualitatively the main observations on dormancy completion made under controlled temperature conditions and relates the parameters of the theory to the measurable characteristics of the system.}
+}
+
+@article { richardson_chill_1974,
+	author = "E. Arlo Richardson and Schuyler D. Seeley and David R. Walker",
+	title = "A Model for Estimating the Completion of Rest for ‘Redhaven’ and ‘Elberta’ Peach Trees1",
+	journal = "HortScience",
+	year = "1974",
+	publisher = "American Society for Horticultural Science",
+	address = "Washington, DC",
+	volume = "9",
+	number = "4",
+	doi = "10.21273/HORTSCI.9.4.331",
+	pages=      "331 - 332",
+	url = "https://journals.ashs.org/hortsci/view/journals/hortsci/9/4/article-p331.xml"
+}

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,34 @@ 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
+
+    # Values are confirmed with chillR package although not an exact match
+    # due to implementation details
+    exp = [1546.5, 1344.0, 1162.0, 1457.5]
+    np.testing.assert_allclose(cu.isel(location=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
+
+    # Values are confirmed with chillR package although not an exact match
+    # 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)

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/data/fr.json

@@ -1385,5 +1385,17 @@
     "description": "Calcul du rayonnement de grandes longueurs d'onde ascendant à partir du rayonnement de grandes longueurs d'onde net et descendant.",
     "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."
+  },
+  "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.",
+    "long_name": "Portions de froid selon le modèle Dynamique",
+    "description": "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."
+  },
+  "CU": {
+    "title": "Unités de froid (chill units) selon le modèle Utah",
+    "abstract": "Les unités de froid (chill units) sont une mesure du potentiel de débourrement de différentes récoltes basée sur le modèle Utah.",
+    "long_name": "Unités de froid selon le modèle Utah",
+    "description": "Les unités de froid (chill units) sont une mesure du potentiel de débourrement de différentes récoltes basée sur le modèle Utah."
   }
 }

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,
+)