From 9cb18338ab58d51aca481e5e05cebf6cea48194d Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Mon, 9 Sep 2024 14:48:24 +0200 Subject: [PATCH 01/36] Add chill_portion to _agro --- xclim/indices/_agro.py | 73 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 72 insertions(+), 1 deletion(-) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 9bfcb0cb4..4c5d560cf 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1,6 +1,5 @@ # noqa: D100 from __future__ import annotations - import warnings from typing import cast @@ -27,6 +26,7 @@ from xclim.indices.helpers import _gather_lat, day_lengths from xclim.indices.stats import standardized_index + # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start # See https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html @@ -1552,3 +1552,74 @@ def hardiness_zones( zones = zones.assign_attrs(units="") return zones + + +# Chill portion constans Dynamic Model 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 + + +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.""" + 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] + ) + delt = np.where(inter_E >= 1, inter_E * xi, 0) + + return delt.cumsum(axis=-1) + + +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"]], + dask="parallelized", + ) + + +@declare_units(tas="[temperature]") +def chill_portion(tas: xarray.DataArray, time_dim: str = "time") -> 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. + time_dim : str, optional + The name of the time dimension (default: "time"). + """ + tas_K = convert_units_to(tas, "K") + return tas_K.groupby(f"{time_dim}.year").apply(_apply_chill_portion_one_season) From f6846b1e9a4213a3aed00cea4c6cb4f536749797 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 9 Sep 2024 12:54:23 +0000 Subject: [PATCH 02/36] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- xclim/indices/_agro.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index bc6ea1298..cb0120651 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1,5 +1,6 @@ # noqa: D100 from __future__ import annotations + from typing import cast import numpy as np @@ -25,7 +26,6 @@ from xclim.indices.helpers import _gather_lat, day_lengths from xclim.indices.stats import standardized_index - # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start # See https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html From 52819cbcc797f062252dec197a12e848421a4906 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Mon, 9 Sep 2024 15:47:06 +0200 Subject: [PATCH 03/36] Add make hourly data to helper function --- xclim/indices/helpers.py | 132 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 131 insertions(+), 1 deletion(-) diff --git a/xclim/indices/helpers.py b/xclim/indices/helpers.py index 4af47b962..5e0f2add6 100644 --- a/xclim/indices/helpers.py +++ b/xclim/indices/helpers.py @@ -6,9 +6,10 @@ """ from __future__ import annotations - from collections.abc import Mapping +from datetime import timedelta from inspect import stack +from math import pi from typing import Any, cast import cf_xarray # noqa: F401, pylint: disable=unused-import @@ -550,3 +551,132 @@ def _gather_lon(da: xr.DataArray) -> xr.DataArray: "Try passing it explicitly (`lon=ds.lon`)." ) raise ValueError(msg) from err + + +def _compute_daytime_temperature( + hour_after_sunrise: xr.DataArray, + tasmin: xr.DataArray, + tasmax: xr.DataArray, + daylength: xr.DataArray, +) -> xr.DataArray: + """Compute daytime temperature based on a sinusoidal profile. + + Minimum temperature is reached at sunrise and maximum temperature 2h before sunset. + + Parameters + ---------- + hours_after_sunrise : xarray.DataArray + Hours after the last sunrise. + tasmin : xarray.DataArray + Daily minimum temperature. + tasmax : xarray.DataArray + Daily maximum temperature. + daylength : xarray.DataArray + Length of the day in hours. + + Returns + ------- + xarray.DataArray + Hourly daytime temperature. + """ + return (tasmax - tasmin) * np.sin( + (pi * hour_after_sunrise) / (daylength + 4) + ) + tasmin + + +def _compute_nighttime_temperature( + hours_after_sunset: xr.DataArray, + tasmin: xr.DataArray, + tas_sunset: xr.DataArray, + daylength: xr.DataArray, +) -> xr.DataArray: + """Compute nighttime temperature based on a logarithmic profile. + + Temperature at sunset is computed from previous daytime temperature, + minimum temperature is reached at sunrise. + + Parameters + ---------- + hours_after_sunset : xarray.DataArray + Hours after the last sunset. + tasmin : xarray.DataArray + Daily minimum temperature. + tas_sunset : xarray.DataArray + Temperature at last sunset. + daylength : xarray.DataArray + Length of the day in hours. + + Returns + ------- + xarray.DataArray + Hourly nighttime temperature. + """ + return tas_sunset - ((tas_sunset - tasmin) / np.log(24 - daylength)) * np.log( + hours_after_sunset + ) + + +def make_hourly_temperature(tasmin: xr.DataArray, tasmax: xr.DataArray) -> xr.DataArray: + """Compute hourly temperatures from tasmin and tasmax. + + Based on the Linvill et al. "Calculating Chilling Hours and Chill Units from Daily + Maximum and Minimum Temperature Observations", HortScience, 1990 + we assume a sinusoidal temperature profile during daylight and a logarithmic decrease after sunset + with tasmin reached at sunsrise and tasmax reached 2h before sunset. + + For simplicity and because we do daily aggregation, we assume that sunrise globally happens at midnight + and the sunsets after `daylength` hours computed via the day_lengths function. + + Parameters + ---------- + tasmin : xarray.DataArray + Daily minimum temperature. + tasmax : xarray.DataArray + Daily maximum temperature. + + Returns + ------- + xarray.DataArray + Hourly temperature. + """ + data = xr.merge([tasmin, tasmax]) + # We add one more timestamp so the resample function includes the last day + data = xr.concat( + [ + data, + data.isel(time=-1).assign_coords( + time=data.isel(time=-1).time + timedelta(days=1) + ), + ], + dim="time", + ) + + daylength = day_lengths(data.time, data.lat) + # Create daily chunks to avoid memory issues after the resampling + data = data.assign( + daylength=daylength, + sunset_temp=_compute_daytime_temperature( + daylength, data.tasmin, data.tasmax, daylength + ), + next_tasmin=data.tasmin.shift(time=-1).ffill("time"), + ).chunk(time=1) + # Compute hourly data by resampling and remove the last time stamp that was added earlier + hourly = data.resample(time="h").ffill().isel(time=slice(0, -1)) + # To avoid "invalid value encountered in log" warning we set hours before sunset to 1 + nighttime_hours = xr.where( + (hourly.time.dt.hour + 1 - hourly.daylength) > 0, + hourly.time.dt.hour + 1 - hourly.daylength, + 1, + ) + return xr.where( + hourly.time.dt.hour < hourly.daylength, + _compute_daytime_temperature( + hourly.time.dt.hour, hourly.tasmin, hourly.tasmax, hourly.daylength + ), + _compute_nighttime_temperature( + nighttime_hours, + hourly.next_tasmin, + hourly.sunset_temp, + hourly.daylength - 1, + ), + ) From d236173fbcd218613ee94f3a9866156241a2f32e Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 9 Sep 2024 13:48:00 +0000 Subject: [PATCH 04/36] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- xclim/indices/helpers.py | 1 + 1 file changed, 1 insertion(+) diff --git a/xclim/indices/helpers.py b/xclim/indices/helpers.py index 5e0f2add6..b211ff675 100644 --- a/xclim/indices/helpers.py +++ b/xclim/indices/helpers.py @@ -6,6 +6,7 @@ """ from __future__ import annotations + from collections.abc import Mapping from datetime import timedelta from inspect import stack From e244e590bc5a8fddc62b882d56071b56a7d1f718 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Mon, 9 Sep 2024 15:50:54 +0200 Subject: [PATCH 05/36] Fix typos --- xclim/indices/_agro.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index cb0120651..4dec00971 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1,6 +1,5 @@ # noqa: D100 from __future__ import annotations - from typing import cast import numpy as np @@ -26,6 +25,7 @@ from xclim.indices.helpers import _gather_lat, day_lengths from xclim.indices.stats import standardized_index + # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start # See https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html @@ -1526,7 +1526,7 @@ def hardiness_zones( return zones -# Chill portion constans Dynamic Model described in Luedeling et al. (2009) +# Chill portion constants Dynamic Model described in Luedeling et al. (2009) E0 = 4153.5 E1 = 12888.8 A0 = 139500 @@ -1556,9 +1556,9 @@ def _chill_portion_one_season(tas_K): inter_E[..., i] = _accumulate_intermediate( inter_E[..., i - 1], xi[..., i - 1], xs[..., i], ak1[..., i] ) - delt = np.where(inter_E >= 1, inter_E * xi, 0) + delta = np.where(inter_E >= 1, inter_E * xi, 0) - return delt.cumsum(axis=-1) + return delta.cumsum(axis=-1) def _apply_chill_portion_one_season(tas_K): From 90f2bad0576e59d4ff32ead1be628582c8ec42cd Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 9 Sep 2024 13:52:01 +0000 Subject: [PATCH 06/36] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- xclim/indices/_agro.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 4dec00971..ce4213d35 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1,5 +1,6 @@ # noqa: D100 from __future__ import annotations + from typing import cast import numpy as np @@ -25,7 +26,6 @@ from xclim.indices.helpers import _gather_lat, day_lengths from xclim.indices.stats import standardized_index - # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start # See https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html From 7d522f66d790719692eee2104f32723550aa78ab Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Mon, 9 Sep 2024 16:00:53 +0200 Subject: [PATCH 07/36] Add utah model --- xclim/indices/_agro.py | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 4dec00971..524faa1fc 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1592,6 +1592,42 @@ def chill_portion(tas: xarray.DataArray, time_dim: str = "time") -> xarray.DataA Hourly temperature. time_dim : str, optional The name of the time dimension (default: "time"). + + Returns + ------- + xr.DataArray """ tas_K = convert_units_to(tas, "K") return tas_K.groupby(f"{time_dim}.year").apply(_apply_chill_portion_one_season) + + +@declare_units(tas="[temperature]") +def utah_model_chill_unit(tas): + """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 + """ + return 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), + ), + ), + ) From 44ab60b464e62a304ddb75500357fa37b85b5d47 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Mon, 9 Sep 2024 16:10:14 +0200 Subject: [PATCH 08/36] Update changelog --- CHANGELOG.rst | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index b28b0e176..df2d9d881 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -4,16 +4,18 @@ 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`). 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. (: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 ^^^^^^^^^ @@ -22,9 +24,6 @@ Bug fixes Breaking changes ^^^^^^^^^^^^^^^^ * `transform` argument of `OTC/dOTC` classes (and child functions) is changed to `normalization`, and `numIterMax` is changed to `num_iter_max` in `utils.optimal_transport` (:pull:`1896`). - -Breaking changes -^^^^^^^^^^^^^^^^ * `platformdirs` is no longer a direct dependency of `xclim`, but `pooch` is required to use many of the new testing functions (installable via `pip install pooch` or `pip install 'xclim[dev]'`). (:pull:`1889`). * The following previously-deprecated functions have now been removed from `xclim`: ``xclim.core.calendar.convert_calendar``, ``xclim.core.calendar.date_range``, ``xclim.core.calendar.date_range_like``, ``xclim.core.calendar.interp_calendar``, ``xclim.core.calendar.days_in_year``, ``xclim.core.calendar.datetime_to_decimal_year``. For guidance on how to migrate to alternatives, see the `version 0.50.0 Breaking changes <#v0-50-0-2024-06-17>`_. (:issue:`1010`, :pull:`1845`). From 2abdbf7fbb722de79638a71defdbf09580a9048b Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Mon, 9 Sep 2024 17:11:33 +0200 Subject: [PATCH 09/36] Implement Zeitsperre comments to changelog and constants into subfunction --- CHANGELOG.rst | 2 +- xclim/indices/_agro.py | 29 +++++++++++++++-------------- 2 files changed, 16 insertions(+), 15 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index df2d9d881..9c1b42352 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -9,7 +9,7 @@ Contributors to this version: Adrien Lamarche (:user:`LamAdr`), Trevor James Smi 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. (:pull:`1909`). +* 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 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 2ea339890..1c895109b 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1,6 +1,5 @@ # noqa: D100 from __future__ import annotations - from typing import cast import numpy as np @@ -26,6 +25,7 @@ from xclim.indices.helpers import _gather_lat, day_lengths from xclim.indices.stats import standardized_index + # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start # See https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html @@ -35,6 +35,8 @@ __all__ = [ "biologically_effective_degree_days", + "chill_portions", + "chill_units", "cool_night_index", "corn_heat_units", "dryness_index", @@ -1526,17 +1528,6 @@ def hardiness_zones( return zones -# Chill portion constants Dynamic Model 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 - - 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) @@ -1545,6 +1536,16 @@ def _accumulate_intermediate(prev_E, prev_xi, curr_xs, 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) @@ -1574,7 +1575,7 @@ def _apply_chill_portion_one_season(tas_K): @declare_units(tas="[temperature]") -def chill_portion(tas: xarray.DataArray, time_dim: str = "time") -> xarray.DataArray: +def chill_portions(tas: xarray.DataArray, time_dim: str = "time") -> xarray.DataArray: """Chill portion based on the dynamic model Chill portions are a measure to estimate the bud breaking potential of different crop. @@ -1602,7 +1603,7 @@ def chill_portion(tas: xarray.DataArray, time_dim: str = "time") -> xarray.DataA @declare_units(tas="[temperature]") -def utah_model_chill_unit(tas): +def chill_units(tas): """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). From c787fff82a90507800824064641d1fb2dadb4931 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 9 Sep 2024 15:13:27 +0000 Subject: [PATCH 10/36] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- xclim/indices/_agro.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 1c895109b..d03a9631a 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1,5 +1,6 @@ # noqa: D100 from __future__ import annotations + from typing import cast import numpy as np @@ -25,7 +26,6 @@ from xclim.indices.helpers import _gather_lat, day_lengths from xclim.indices.stats import standardized_index - # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start # See https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html From bc382eb1306fab23c5f687e9a85ce3f91e9ad39e Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Mon, 9 Sep 2024 17:17:40 +0200 Subject: [PATCH 11/36] Use numpy pi --- xclim/indices/helpers.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/xclim/indices/helpers.py b/xclim/indices/helpers.py index b211ff675..53ab9fbd3 100644 --- a/xclim/indices/helpers.py +++ b/xclim/indices/helpers.py @@ -6,11 +6,9 @@ """ from __future__ import annotations - from collections.abc import Mapping from datetime import timedelta from inspect import stack -from math import pi from typing import Any, cast import cf_xarray # noqa: F401, pylint: disable=unused-import @@ -581,7 +579,7 @@ def _compute_daytime_temperature( Hourly daytime temperature. """ return (tasmax - tasmin) * np.sin( - (pi * hour_after_sunrise) / (daylength + 4) + (np.pi * hour_after_sunrise) / (daylength + 4) ) + tasmin From 4d1b9d7952472b74a9b1519ff247981ad0165d9d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 9 Sep 2024 15:20:16 +0000 Subject: [PATCH 12/36] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- xclim/indices/helpers.py | 1 + 1 file changed, 1 insertion(+) diff --git a/xclim/indices/helpers.py b/xclim/indices/helpers.py index 53ab9fbd3..f7eb36e2f 100644 --- a/xclim/indices/helpers.py +++ b/xclim/indices/helpers.py @@ -6,6 +6,7 @@ """ from __future__ import annotations + from collections.abc import Mapping from datetime import timedelta from inspect import stack From c79eb834580bc8c863efd282bf68cfd5351fce30 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Tue, 10 Sep 2024 10:55:07 +0200 Subject: [PATCH 13/36] Add chill_units test to test_indices --- tests/test_indices.py | 26 +++++++++++++++++++- xclim/indicators/atmos/_temperature.py | 34 +++++++++++++++++++++++++- xclim/indices/_agro.py | 28 ++++++++++++--------- xclim/testing/helpers.py | 2 +- 4 files changed, 76 insertions(+), 14 deletions(-) diff --git a/tests/test_indices.py b/tests/test_indices.py index d1da9c54c..29f247939 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -12,7 +12,6 @@ # saving the results in a reference netcdf dataset. We could then compare the hailstorm output to this reference as # a first line of defense. from __future__ import annotations - import calendar import numpy as np @@ -29,6 +28,7 @@ from xclim.core.options import set_options from xclim.core.units import convert_units_to, units + K2C = 273.15 @@ -346,6 +346,30 @@ 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): + assert False + + 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")) diff --git a/xclim/indicators/atmos/_temperature.py b/xclim/indicators/atmos/_temperature.py index 661169468..3a982996f 100644 --- a/xclim/indicators/atmos/_temperature.py +++ b/xclim/indicators/atmos/_temperature.py @@ -4,9 +4,14 @@ from xclim import indices from xclim.core import cfchecks -from xclim.core.indicator import Daily, Indicator, ResamplingIndicatorWithIndexing +from xclim.core.indicator import ( + Daily, + Indicator, + ResamplingIndicatorWithIndexing, +) from xclim.core.utils import InputKind + __all__ = [ "australian_hardiness_zones", "biologically_effective_degree_days", @@ -100,6 +105,13 @@ class Temp(Daily): keywords = "temperature" +class TempHourly(Indicator): + """Temperature indicators involving hourly temperature.""" + + src_freq = "H" + keywords = "temperature" + + class TempWithIndexing(ResamplingIndicatorWithIndexing): """Indicators involving daily temperature and adding an indexing possibility.""" @@ -1445,3 +1457,23 @@ def cfcheck(self, tas, snd=None): compute=indices.hardiness_zones, parameters={"method": "usda"}, ) + +chill_portions = TempHourly( + title="Chill portions", + identifier="cp", + units="", + long_name="Chill portions after the Dynamic Model", + allowed_periods=["A"], + varname="cp", + compute=indices.chill_portions, +) + +chill_units = TempHourly( + title="Chill units", + identifier="cu", + units="", + long_name="Chill units after the Utah Model", + allowed_periods=["A"], + varname="cu", + compute=indices.chill_units, +) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index d03a9631a..5a0ea2cfe 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1,6 +1,5 @@ # noqa: D100 from __future__ import annotations - from typing import cast import numpy as np @@ -26,6 +25,7 @@ from xclim.indices.helpers import _gather_lat, day_lengths from xclim.indices.stats import standardized_index + # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start # See https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html @@ -1603,7 +1603,7 @@ def chill_portions(tas: xarray.DataArray, time_dim: str = "time") -> xarray.Data @declare_units(tas="[temperature]") -def chill_units(tas): +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). @@ -1619,16 +1619,22 @@ def chill_units(tas): ------- xr.DataArray """ - return xarray.where( - (tas <= 1.4) | ((tas > 12.4) & (tas <= 15.9)), - 0, + tas = convert_units_to(tas, "degC") + cu = ( xarray.where( - ((tas > 1.4) & (tas <= 2.4)) | ((tas > 9.1) & (tas <= 12.4)), - 0.5, + (tas <= 1.4) | ((tas > 12.4) & (tas <= 15.9)), + 0, xarray.where( - (tas > 2.4) & (tas <= 9.1), - 1, - xarray.where((tas > 15.9) & (tas <= 17.9), -0.5, -1), + ((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), + ), ), - ), + ) + .assign_attrs(units="") + .rename("cu") ) + return cu.resample(time=freq).sum() diff --git a/xclim/testing/helpers.py b/xclim/testing/helpers.py index f715f1586..9d45ee717 100644 --- a/xclim/testing/helpers.py +++ b/xclim/testing/helpers.py @@ -1,7 +1,6 @@ """Module for loading testing data.""" from __future__ import annotations - import logging import os import warnings @@ -22,6 +21,7 @@ shortwave_upwelling_radiation_from_net_downwelling, ) + logger = logging.getLogger("xclim") __all__ = [ From de45168641bb952db64ad62f126c3c41ae47c363 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 10 Sep 2024 08:56:08 +0000 Subject: [PATCH 14/36] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- tests/test_indices.py | 2 +- xclim/indicators/atmos/_temperature.py | 7 +------ xclim/indices/_agro.py | 2 +- xclim/testing/helpers.py | 2 +- 4 files changed, 4 insertions(+), 9 deletions(-) diff --git a/tests/test_indices.py b/tests/test_indices.py index 29f247939..dca5474f6 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -12,6 +12,7 @@ # saving the results in a reference netcdf dataset. We could then compare the hailstorm output to this reference as # a first line of defense. from __future__ import annotations + import calendar import numpy as np @@ -28,7 +29,6 @@ from xclim.core.options import set_options from xclim.core.units import convert_units_to, units - K2C = 273.15 diff --git a/xclim/indicators/atmos/_temperature.py b/xclim/indicators/atmos/_temperature.py index 3a982996f..37ee83f23 100644 --- a/xclim/indicators/atmos/_temperature.py +++ b/xclim/indicators/atmos/_temperature.py @@ -4,14 +4,9 @@ from xclim import indices from xclim.core import cfchecks -from xclim.core.indicator import ( - Daily, - Indicator, - ResamplingIndicatorWithIndexing, -) +from xclim.core.indicator import Daily, Indicator, ResamplingIndicatorWithIndexing from xclim.core.utils import InputKind - __all__ = [ "australian_hardiness_zones", "biologically_effective_degree_days", diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 5a0ea2cfe..e7abc98cf 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1,5 +1,6 @@ # noqa: D100 from __future__ import annotations + from typing import cast import numpy as np @@ -25,7 +26,6 @@ from xclim.indices.helpers import _gather_lat, day_lengths from xclim.indices.stats import standardized_index - # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start # See https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html diff --git a/xclim/testing/helpers.py b/xclim/testing/helpers.py index 9d45ee717..f715f1586 100644 --- a/xclim/testing/helpers.py +++ b/xclim/testing/helpers.py @@ -1,6 +1,7 @@ """Module for loading testing data.""" from __future__ import annotations + import logging import os import warnings @@ -21,7 +22,6 @@ shortwave_upwelling_radiation_from_net_downwelling, ) - logger = logging.getLogger("xclim") __all__ = [ From 177bee01d6d516351ceac51a743f9fb60b90429a Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Tue, 10 Sep 2024 12:20:35 +0200 Subject: [PATCH 15/36] Add make_hourly_temperature test --- tests/test_helpers.py | 42 ++++++++++++++++++++++++++++++++++++++++ xclim/indices/helpers.py | 28 +++++++++++++++++++++++---- 2 files changed, 66 insertions(+), 4 deletions(-) diff --git a/tests/test_helpers.py b/tests/test_helpers.py index 699cb3b6b..c8fff785a 100644 --- a/tests/test_helpers.py +++ b/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) diff --git a/xclim/indices/helpers.py b/xclim/indices/helpers.py index f7eb36e2f..21a6d01d0 100644 --- a/xclim/indices/helpers.py +++ b/xclim/indices/helpers.py @@ -6,7 +6,6 @@ """ from __future__ import annotations - from collections.abc import Mapping from datetime import timedelta from inspect import stack @@ -616,6 +615,27 @@ def _compute_nighttime_temperature( ) +def _add_one_day(time: xr.DataArray) -> xr.DataArray: + """Add one day to a time coordinate. + + Depending on the calendar/dtype of the time array we need can use numpy's or + datetime's (for cftimes) timedelta. + + Parameters + ---------- + time : xr.DataArray + Time coordinate. + + Returns + ------- + xr.DataArray + Next day. + """ + if time.dtype == "O": + return time + timedelta(days=1) + return time + np.timedelta64(1, "D") + + def make_hourly_temperature(tasmin: xr.DataArray, tasmax: xr.DataArray) -> xr.DataArray: """Compute hourly temperatures from tasmin and tasmax. @@ -624,7 +644,7 @@ def make_hourly_temperature(tasmin: xr.DataArray, tasmax: xr.DataArray) -> xr.Da we assume a sinusoidal temperature profile during daylight and a logarithmic decrease after sunset with tasmin reached at sunsrise and tasmax reached 2h before sunset. - For simplicity and because we do daily aggregation, we assume that sunrise globally happens at midnight + For simplicity and because it's used for daily aggregation, we assume that sunrise globally happens at midnight and the sunsets after `daylength` hours computed via the day_lengths function. Parameters @@ -645,7 +665,7 @@ def make_hourly_temperature(tasmin: xr.DataArray, tasmax: xr.DataArray) -> xr.Da [ data, data.isel(time=-1).assign_coords( - time=data.isel(time=-1).time + timedelta(days=1) + time=_add_one_day(data.isel(time=-1).time) ), ], dim="time", @@ -679,4 +699,4 @@ def make_hourly_temperature(tasmin: xr.DataArray, tasmax: xr.DataArray) -> xr.Da hourly.sunset_temp, hourly.daylength - 1, ), - ) + ).assign_attrs(units=tasmin.units) From 7dda632d302455f3d6c155e2ce81a7769bc15148 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 10 Sep 2024 10:21:35 +0000 Subject: [PATCH 16/36] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- xclim/indices/helpers.py | 1 + 1 file changed, 1 insertion(+) diff --git a/xclim/indices/helpers.py b/xclim/indices/helpers.py index 21a6d01d0..6779e12df 100644 --- a/xclim/indices/helpers.py +++ b/xclim/indices/helpers.py @@ -6,6 +6,7 @@ """ from __future__ import annotations + from collections.abc import Mapping from datetime import timedelta from inspect import stack From 3bbd6935479609cdc0d90a66b96cd05bc964edd0 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Mon, 16 Sep 2024 23:13:37 +0200 Subject: [PATCH 17/36] Update descriptions --- xclim/indicators/atmos/_temperature.py | 13 +++++++++++++ xclim/indices/helpers.py | 2 +- 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/xclim/indicators/atmos/_temperature.py b/xclim/indicators/atmos/_temperature.py index 37ee83f23..0fc305d39 100644 --- a/xclim/indicators/atmos/_temperature.py +++ b/xclim/indicators/atmos/_temperature.py @@ -7,6 +7,7 @@ from xclim.core.indicator import Daily, Indicator, ResamplingIndicatorWithIndexing from xclim.core.utils import InputKind + __all__ = [ "australian_hardiness_zones", "biologically_effective_degree_days", @@ -1457,6 +1458,15 @@ def cfcheck(self, tas, snd=None): title="Chill portions", identifier="cp", units="", + 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). " + "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=["A"], varname="cp", @@ -1467,6 +1477,9 @@ def cfcheck(self, tas, snd=None): title="Chill units", identifier="cu", units="", + 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=["A"], varname="cu", diff --git a/xclim/indices/helpers.py b/xclim/indices/helpers.py index 21a6d01d0..0d651fd95 100644 --- a/xclim/indices/helpers.py +++ b/xclim/indices/helpers.py @@ -618,7 +618,7 @@ def _compute_nighttime_temperature( def _add_one_day(time: xr.DataArray) -> xr.DataArray: """Add one day to a time coordinate. - Depending on the calendar/dtype of the time array we need can use numpy's or + Depending on the calendar/dtype of the time array we need to use numpy's or datetime's (for cftimes) timedelta. Parameters From 0fcceb49b1f336e7179e938b8e4eff2ca15f1418 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 16 Sep 2024 21:14:32 +0000 Subject: [PATCH 18/36] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- xclim/indicators/atmos/_temperature.py | 1 - 1 file changed, 1 deletion(-) diff --git a/xclim/indicators/atmos/_temperature.py b/xclim/indicators/atmos/_temperature.py index 0fc305d39..9115b3fcf 100644 --- a/xclim/indicators/atmos/_temperature.py +++ b/xclim/indicators/atmos/_temperature.py @@ -7,7 +7,6 @@ from xclim.core.indicator import Daily, Indicator, ResamplingIndicatorWithIndexing from xclim.core.utils import InputKind - __all__ = [ "australian_hardiness_zones", "biologically_effective_degree_days", From 6585a9c9adc4d2fc7508604f02f43a2a3778cf55 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Mon, 16 Sep 2024 23:48:32 +0200 Subject: [PATCH 19/36] Take maximum along time axis for chill_portions --- xclim/indicators/atmos/_temperature.py | 7 +++++++ xclim/indices/_agro.py | 16 +++++++++------- 2 files changed, 16 insertions(+), 7 deletions(-) diff --git a/xclim/indicators/atmos/_temperature.py b/xclim/indicators/atmos/_temperature.py index 0fc305d39..9a937eed3 100644 --- a/xclim/indicators/atmos/_temperature.py +++ b/xclim/indicators/atmos/_temperature.py @@ -11,6 +11,8 @@ __all__ = [ "australian_hardiness_zones", "biologically_effective_degree_days", + "chill_portions", + "chill_units", "cold_spell_days", "cold_spell_duration_index", "cold_spell_frequency", @@ -1458,6 +1460,9 @@ def cfcheck(self, tas, snd=None): title="Chill portions", identifier="cp", units="", + standard_name="chill_portions", + # TODO: check what this does + cell_methods="", 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). " @@ -1480,6 +1485,8 @@ def cfcheck(self, tas, snd=None): 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.", + # TODO: check what this does + standard_name="chill_units", long_name="Chill units after the Utah Model", allowed_periods=["A"], varname="cu", diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index e7abc98cf..567e840be 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1,6 +1,5 @@ # noqa: D100 from __future__ import annotations - from typing import cast import numpy as np @@ -26,6 +25,7 @@ from xclim.indices.helpers import _gather_lat, day_lengths from xclim.indices.stats import standardized_index + # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start # See https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html @@ -1571,11 +1571,11 @@ def _apply_chill_portion_one_season(tas_K): input_core_dims=[["time"]], output_core_dims=[["time"]], dask="parallelized", - ) + ).max("time") @declare_units(tas="[temperature]") -def chill_portions(tas: xarray.DataArray, time_dim: str = "time") -> xarray.DataArray: +def chill_portions(tas: xarray.DataArray, freq: str = "YS") -> xarray.DataArray: """Chill portion based on the dynamic model Chill portions are a measure to estimate the bud breaking potential of different crop. @@ -1591,15 +1591,17 @@ def chill_portions(tas: xarray.DataArray, time_dim: str = "time") -> xarray.Data ---------- tas : xr.DataArray Hourly temperature. - time_dim : str, optional - The name of the time dimension (default: "time"). Returns ------- xr.DataArray """ - tas_K = convert_units_to(tas, "K") - return tas_K.groupby(f"{time_dim}.year").apply(_apply_chill_portion_one_season) + tas_K: xarray.DataArray = convert_units_to(tas, "K") + return ( + tas_K.resample(time=freq) + .apply(_apply_chill_portion_one_season) + .assign_attrs(units="") + ) @declare_units(tas="[temperature]") From 98900dfa10ade7f3c6d1f5039f6d17a356e7c3b6 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 17 Sep 2024 07:16:34 +0000 Subject: [PATCH 20/36] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- xclim/indices/_agro.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 567e840be..b451a0bb2 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1,5 +1,6 @@ # noqa: D100 from __future__ import annotations + from typing import cast import numpy as np @@ -25,7 +26,6 @@ from xclim.indices.helpers import _gather_lat, day_lengths from xclim.indices.stats import standardized_index - # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start # See https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html From 0a496bea4867dd6eb598d81224ceadf6b4807722 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Tue, 17 Sep 2024 10:03:47 +0200 Subject: [PATCH 21/36] Test chill portions --- tests/test_indices.py | 8 +++++--- xclim/indices/_agro.py | 3 ++- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/tests/test_indices.py b/tests/test_indices.py index dca5474f6..5ef662a6e 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -12,7 +12,6 @@ # saving the results in a reference netcdf dataset. We could then compare the hailstorm output to this reference as # a first line of defense. from __future__ import annotations - import calendar import numpy as np @@ -29,6 +28,7 @@ from xclim.core.options import set_options from xclim.core.units import convert_units_to, units + K2C = 273.15 @@ -347,7 +347,9 @@ def test_bedd(self, method, end_date, deg_days, max_deg_days): np.testing.assert_array_equal(bedd, bedd_high_lat) def test_chill_portions(self, tas_series): - assert False + 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 @@ -365,7 +367,7 @@ def test_chill_units(self, tas_series): + num_cu_min_1 * [20.0] ) + K2C, - freq="H", + 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 diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 567e840be..bbd44757e 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1599,7 +1599,8 @@ def chill_portions(tas: xarray.DataArray, freq: str = "YS") -> xarray.DataArray: tas_K: xarray.DataArray = convert_units_to(tas, "K") return ( tas_K.resample(time=freq) - .apply(_apply_chill_portion_one_season) + .map(_apply_chill_portion_one_season) + .rename("cp") .assign_attrs(units="") ) From 72dc9806bd4a079e0ba3adc1c8c0bc8752d1b565 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 17 Sep 2024 08:06:22 +0000 Subject: [PATCH 22/36] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- tests/test_indices.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_indices.py b/tests/test_indices.py index 5ef662a6e..7269931c4 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -12,6 +12,7 @@ # saving the results in a reference netcdf dataset. We could then compare the hailstorm output to this reference as # a first line of defense. from __future__ import annotations + import calendar import numpy as np @@ -28,7 +29,6 @@ from xclim.core.options import set_options from xclim.core.units import convert_units_to, units - K2C = 273.15 From 14b44bd4dbc0c822230ee71afb1fd8a63f45c592 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Thu, 19 Sep 2024 09:40:24 +0200 Subject: [PATCH 23/36] Apply suggestions from aulemahal Co-authored-by: Pascal Bourgault --- xclim/indicators/atmos/_temperature.py | 5 +---- xclim/indices/_agro.py | 7 +++---- xclim/indices/helpers.py | 2 +- 3 files changed, 5 insertions(+), 9 deletions(-) diff --git a/xclim/indicators/atmos/_temperature.py b/xclim/indicators/atmos/_temperature.py index f898133d9..87e6036e8 100644 --- a/xclim/indicators/atmos/_temperature.py +++ b/xclim/indicators/atmos/_temperature.py @@ -105,7 +105,7 @@ class Temp(Daily): class TempHourly(Indicator): """Temperature indicators involving hourly temperature.""" - src_freq = "H" + src_freq = "h" keywords = "temperature" @@ -1459,7 +1459,6 @@ def cfcheck(self, tas, snd=None): title="Chill portions", identifier="cp", units="", - standard_name="chill_portions", # TODO: check what this does cell_methods="", description="Chill portions are a measure to estimate the bud breaking potential of different crops. " @@ -1484,8 +1483,6 @@ def cfcheck(self, tas, snd=None): 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.", - # TODO: check what this does - standard_name="chill_units", long_name="Chill units after the Utah Model", allowed_periods=["A"], varname="cu", diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 41baba872..bcffc4966 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1594,7 +1594,8 @@ def chill_portions(tas: xarray.DataArray, freq: str = "YS") -> xarray.DataArray: Returns ------- - xr.DataArray + xr.DataArray, [unitless] + Chill portions after the Dynamic Model """ tas_K: xarray.DataArray = convert_units_to(tas, "K") return ( @@ -1637,7 +1638,5 @@ def chill_units(tas: xarray.DataArray, freq: str = "YS") -> xarray.DataArray: ), ), ) - .assign_attrs(units="") - .rename("cu") ) - return cu.resample(time=freq).sum() + return cu.resample(time=freq).sum().assign_attrs(units="") diff --git a/xclim/indices/helpers.py b/xclim/indices/helpers.py index 4c7d921c6..b3056525a 100644 --- a/xclim/indices/helpers.py +++ b/xclim/indices/helpers.py @@ -655,7 +655,7 @@ def make_hourly_temperature(tasmin: xr.DataArray, tasmax: xr.DataArray) -> xr.Da with tasmin reached at sunsrise and tasmax reached 2h before sunset. For simplicity and because it's used for daily aggregation, we assume that sunrise globally happens at midnight - and the sunsets after `daylength` hours computed via the day_lengths function. + and the sunsets after `daylength` hours computed via the :py:func:`day_lengths` function. Parameters ---------- From 9d65ca1501afbaa23f3c1c02cd1b6e29b26a1b19 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 19 Sep 2024 07:41:25 +0000 Subject: [PATCH 24/36] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- xclim/indices/_agro.py | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index bcffc4966..f1a94e65d 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1624,19 +1624,17 @@ def chill_units(tas: xarray.DataArray, freq: str = "YS") -> xarray.DataArray: xr.DataArray """ tas = convert_units_to(tas, "degC") - cu = ( + cu = xarray.where( + (tas <= 1.4) | ((tas > 12.4) & (tas <= 15.9)), + 0, xarray.where( - (tas <= 1.4) | ((tas > 12.4) & (tas <= 15.9)), - 0, + ((tas > 1.4) & (tas <= 2.4)) | ((tas > 9.1) & (tas <= 12.4)), + 0.5, 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), - ), + (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="") From f266153a02329b3c415e8d9b62b4d36e9b42d307 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Thu, 19 Sep 2024 09:45:14 +0200 Subject: [PATCH 25/36] Add unitless to chill_units and abstract to chill_portion --- xclim/indicators/atmos/_temperature.py | 4 ++++ xclim/indices/_agro.py | 5 +++-- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/xclim/indicators/atmos/_temperature.py b/xclim/indicators/atmos/_temperature.py index f898133d9..524b650f0 100644 --- a/xclim/indicators/atmos/_temperature.py +++ b/xclim/indicators/atmos/_temperature.py @@ -7,6 +7,7 @@ from xclim.core.indicator import Daily, Indicator, ResamplingIndicatorWithIndexing from xclim.core.utils import InputKind + __all__ = [ "australian_hardiness_zones", "biologically_effective_degree_days", @@ -1464,6 +1465,9 @@ def cfcheck(self, tas, snd=None): cell_methods="", 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 " diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 098031dae..3e0778c31 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1,6 +1,5 @@ # noqa: D100 from __future__ import annotations - from typing import cast import numpy as np @@ -26,6 +25,7 @@ from xclim.indices.helpers import _gather_lat, day_lengths from xclim.indices.stats import standardized_index + # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start # See https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html @@ -1620,7 +1620,8 @@ def chill_units(tas: xarray.DataArray, freq: str = "YS") -> xarray.DataArray: Returns ------- - xr.DataArray + xr.DataArray, [unitless] + Chill units using the Utah model """ tas = convert_units_to(tas, "degC") cu = ( From d0eea99d61f080a56da365ea841d195a4657298a Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Thu, 19 Sep 2024 10:21:56 +0200 Subject: [PATCH 26/36] Apply more suggestions Remove explicit chunking in make_hourly Remove varname from indicators use sum instead of combination of cumsum and max for cp --- xclim/indicators/atmos/_temperature.py | 3 --- xclim/indices/_agro.py | 7 +++---- xclim/indices/helpers.py | 2 +- 3 files changed, 4 insertions(+), 8 deletions(-) diff --git a/xclim/indicators/atmos/_temperature.py b/xclim/indicators/atmos/_temperature.py index fd2b66e6b..970b0d985 100644 --- a/xclim/indicators/atmos/_temperature.py +++ b/xclim/indicators/atmos/_temperature.py @@ -7,7 +7,6 @@ from xclim.core.indicator import Daily, Indicator, ResamplingIndicatorWithIndexing from xclim.core.utils import InputKind - __all__ = [ "australian_hardiness_zones", "biologically_effective_degree_days", @@ -1476,7 +1475,6 @@ def cfcheck(self, tas, snd=None): "especially in moderate climates like Israel, California or Spain.", long_name="Chill portions after the Dynamic Model", allowed_periods=["A"], - varname="cp", compute=indices.chill_portions, ) @@ -1489,6 +1487,5 @@ def cfcheck(self, tas, snd=None): "actually decrease the potential for bud breaking.", long_name="Chill units after the Utah Model", allowed_periods=["A"], - varname="cu", compute=indices.chill_units, ) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index f8a3b876e..f3e3fa7c7 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1,5 +1,6 @@ # noqa: D100 from __future__ import annotations + from typing import cast import numpy as np @@ -25,7 +26,6 @@ from xclim.indices.helpers import _gather_lat, day_lengths from xclim.indices.stats import standardized_index - # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start # See https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html @@ -1559,7 +1559,7 @@ def _chill_portion_one_season(tas_K): ) delta = np.where(inter_E >= 1, inter_E * xi, 0) - return delta.cumsum(axis=-1) + return delta def _apply_chill_portion_one_season(tas_K): @@ -1571,7 +1571,7 @@ def _apply_chill_portion_one_season(tas_K): input_core_dims=[["time"]], output_core_dims=[["time"]], dask="parallelized", - ).max("time") + ).sum("time") @declare_units(tas="[temperature]") @@ -1601,7 +1601,6 @@ def chill_portions(tas: xarray.DataArray, freq: str = "YS") -> xarray.DataArray: return ( tas_K.resample(time=freq) .map(_apply_chill_portion_one_season) - .rename("cp") .assign_attrs(units="") ) diff --git a/xclim/indices/helpers.py b/xclim/indices/helpers.py index b3056525a..cc1a90dcb 100644 --- a/xclim/indices/helpers.py +++ b/xclim/indices/helpers.py @@ -689,7 +689,7 @@ def make_hourly_temperature(tasmin: xr.DataArray, tasmax: xr.DataArray) -> xr.Da daylength, data.tasmin, data.tasmax, daylength ), next_tasmin=data.tasmin.shift(time=-1).ffill("time"), - ).chunk(time=1) + ) # Compute hourly data by resampling and remove the last time stamp that was added earlier hourly = data.resample(time="h").ffill().isel(time=slice(0, -1)) # To avoid "invalid value encountered in log" warning we set hours before sunset to 1 From 0659e6a7a3439a0d56f0041df8fb392e60564f05 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Thu, 19 Sep 2024 10:51:07 +0200 Subject: [PATCH 27/36] Integrate two more comments --- xclim/indices/helpers.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/xclim/indices/helpers.py b/xclim/indices/helpers.py index cc1a90dcb..f4943d97b 100644 --- a/xclim/indices/helpers.py +++ b/xclim/indices/helpers.py @@ -688,16 +688,16 @@ def make_hourly_temperature(tasmin: xr.DataArray, tasmax: xr.DataArray) -> xr.Da sunset_temp=_compute_daytime_temperature( daylength, data.tasmin, data.tasmax, daylength ), - next_tasmin=data.tasmin.shift(time=-1).ffill("time"), + next_tasmin=data.tasmin.shift(time=-1), ) # Compute hourly data by resampling and remove the last time stamp that was added earlier hourly = data.resample(time="h").ffill().isel(time=slice(0, -1)) + # To avoid "invalid value encountered in log" warning we set hours before sunset to 1 - nighttime_hours = xr.where( - (hourly.time.dt.hour + 1 - hourly.daylength) > 0, - hourly.time.dt.hour + 1 - hourly.daylength, - 1, - ) + nighttime_hours = nighttime_hours = ( + hourly.time.dt.hour + 1 - hourly.daylength + ).clip(1) + return xr.where( hourly.time.dt.hour < hourly.daylength, _compute_daytime_temperature( From ad5a7ec8278eea1de9860d900e83b7d7986468f1 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Thu, 19 Sep 2024 17:59:34 +0200 Subject: [PATCH 28/36] Use ResamplingIndicatorWithIndexing for chill_units --- tests/test_atmos.py | 14 ++++++++++++++ xclim/indicators/atmos/_temperature.py | 15 +++++++++++---- xclim/indices/_agro.py | 1 + 3 files changed, 26 insertions(+), 4 deletions(-) diff --git a/tests/test_atmos.py b/tests/test_atmos.py index 23929550d..f0b86ceee 100644 --- a/tests/test_atmos.py +++ b/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,16 @@ 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) diff --git a/xclim/indicators/atmos/_temperature.py b/xclim/indicators/atmos/_temperature.py index 970b0d985..bdc94548a 100644 --- a/xclim/indicators/atmos/_temperature.py +++ b/xclim/indicators/atmos/_temperature.py @@ -116,6 +116,13 @@ class TempWithIndexing(ResamplingIndicatorWithIndexing): 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", @@ -1459,8 +1466,7 @@ def cfcheck(self, tas, snd=None): title="Chill portions", identifier="cp", units="", - # TODO: check what this does - cell_methods="", + 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). ", @@ -1478,14 +1484,15 @@ def cfcheck(self, tas, snd=None): compute=indices.chill_portions, ) -chill_units = TempHourly( +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=["A"], + allowed_periods=["Y"], compute=indices.chill_units, ) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index f3e3fa7c7..cbff71ca6 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1570,6 +1570,7 @@ def _apply_chill_portion_one_season(tas_K): tas_K, input_core_dims=[["time"]], output_core_dims=[["time"]], + output_dtypes=[tas_K.dtype], dask="parallelized", ).sum("time") From 6dc26baef4c506a2c673e72fc15c0eda59a54442 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Thu, 19 Sep 2024 22:24:49 +0200 Subject: [PATCH 29/36] Add indexing to chill_portions --- tests/test_atmos.py | 21 +++++++++++++++++++++ xclim/indicators/atmos/_temperature.py | 12 ++++++++---- xclim/indices/_agro.py | 8 ++++++-- 3 files changed, 35 insertions(+), 6 deletions(-) diff --git a/tests/test_atmos.py b/tests/test_atmos.py index f0b86ceee..e7e6eacbd 100644 --- a/tests/test_atmos.py +++ b/tests/test_atmos.py @@ -638,3 +638,24 @@ def test_chill_units(atmosds): 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) diff --git a/xclim/indicators/atmos/_temperature.py b/xclim/indicators/atmos/_temperature.py index bdc94548a..be5ee18da 100644 --- a/xclim/indicators/atmos/_temperature.py +++ b/xclim/indicators/atmos/_temperature.py @@ -4,7 +4,12 @@ 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__ = [ @@ -102,10 +107,9 @@ class Temp(Daily): keywords = "temperature" -class TempHourly(Indicator): +class TempHourly(Hourly): """Temperature indicators involving hourly temperature.""" - src_freq = "h" keywords = "temperature" @@ -1480,7 +1484,7 @@ def cfcheck(self, tas, snd=None): "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=["A"], + allowed_periods=["Y"], compute=indices.chill_portions, ) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index cbff71ca6..0120def86 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1576,7 +1576,9 @@ def _apply_chill_portion_one_season(tas_K): @declare_units(tas="[temperature]") -def chill_portions(tas: xarray.DataArray, freq: str = "YS") -> xarray.DataArray: +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. @@ -1598,7 +1600,9 @@ def chill_portions(tas: xarray.DataArray, freq: str = "YS") -> xarray.DataArray: xr.DataArray, [unitless] Chill portions after the Dynamic Model """ - tas_K: xarray.DataArray = convert_units_to(tas, "K") + 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) From f9a392fe6aac4d779298ac9dc5ebb97a31d8aa0f Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Wed, 25 Sep 2024 14:38:51 +0200 Subject: [PATCH 30/36] Update expected test values --- tests/test_atmos.py | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/tests/test_atmos.py b/tests/test_atmos.py index e7e6eacbd..f55e32fbe 100644 --- a/tests/test_atmos.py +++ b/tests/test_atmos.py @@ -650,12 +650,5 @@ def test_chill_portions(atmosds): # 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) + exp = [np.nan, 99.91534493, 92.5473925, 99.03177047, np.nan] + np.testing.assert_allclose(cp.isel(location=0), exp, rtol=1e-03) From 43a9d13d52b060128ee150aca1b2404dd1766b0b Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Wed, 25 Sep 2024 16:23:07 +0200 Subject: [PATCH 31/36] Fix nans in cu computation --- tests/test_atmos.py | 4 ++-- xclim/indices/_agro.py | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/test_atmos.py b/tests/test_atmos.py index f55e32fbe..95dcc56c2 100644 --- a/tests/test_atmos.py +++ b/tests/test_atmos.py @@ -636,8 +636,8 @@ def test_chill_units(atmosds): 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) + 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): diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 0120def86..1c552abeb 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1642,4 +1642,5 @@ def chill_units(tas: xarray.DataArray, freq: str = "YS") -> xarray.DataArray: ), ), ) + cu = cu.where(tas.notnull()) return cu.resample(time=freq).sum().assign_attrs(units="") From 56290de1cca3aa158e57e342d145b86441128cc9 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Thu, 26 Sep 2024 10:10:52 +0200 Subject: [PATCH 32/36] Add references --- docs/references.bib | 42 ++++++++++++++++++++++++++++++++++++++++++ tests/test_atmos.py | 4 ++++ xclim/indices/_agro.py | 37 +++++++++++++++++++++++++++++++++++++ 3 files changed, 83 insertions(+) diff --git a/docs/references.bib b/docs/references.bib index 9f2419dd9..0801cec37 100644 --- a/docs/references.bib +++ b/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" +} \ No newline at end of file diff --git a/tests/test_atmos.py b/tests/test_atmos.py index 95dcc56c2..80757e5f7 100644 --- a/tests/test_atmos.py +++ b/tests/test_atmos.py @@ -636,6 +636,8 @@ def test_chill_units(atmosds): 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) @@ -650,5 +652,7 @@ def test_chill_portions(atmosds): # 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) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 1c552abeb..277d0dfec 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1599,6 +1599,30 @@ def chill_portions( ------- xr.DataArray, [unitless] Chill portions after the Dynamic Model + + Notes + ----- + Typically, this indicator is computed for a period of the year. You can use the `**indexer` arguments + of `select_time` in combination with the `freq` argument to select e.g. a winter period: + + .. code-block:: python + + cp = chill_portions(tas, date_bounds=("09-01", "03-30"), freq="YS-JUL") + + Note that incomplete periods will lead to NaNs. + + Examples + -------- + >>> from xclim.indices import chill_portions + >>> from xclim.indices.helpers import make_hourly_temperature + >>> tasmin = xr.open_dataset(path_to_tasmin_file).tasmin + >>> tasmax = xr.open_dataset(path_to_tasmax_file).tasmax + >>> tas_hourly = make_hourly_temperature(tasmin, tasmax) + >>> cp = chill_portions(tasmin) + + References + ---------- + :cite:cts:`fishman_chill_1987,luedeling_chill_2009` """ tas_K: xarray.DataArray = select_time( convert_units_to(tas, "K"), drop=True, **indexer @@ -1627,6 +1651,19 @@ def chill_units(tas: xarray.DataArray, freq: str = "YS") -> xarray.DataArray: ------- xr.DataArray, [unitless] Chill units using the Utah model + + Examples + -------- + >>> from xclim.indices import chill_units + >>> from xclim.indices.helpers import make_hourly_temperature + >>> tasmin = xr.open_dataset(path_to_tasmin_file).tasmin + >>> tasmax = xr.open_dataset(path_to_tasmax_file).tasmax + >>> tas_hourly = make_hourly_temperature(tasmin, tasmax) + >>> cu = chill_units(tasmin) + + References + ---------- + :cite:cts:`richardson_chill_1974` """ tas = convert_units_to(tas, "degC") cu = xarray.where( From 765d67d954e755fe5d461967cfc6d1c7342afba6 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 26 Sep 2024 08:11:53 +0000 Subject: [PATCH 33/36] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- docs/references.bib | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/references.bib b/docs/references.bib index 0801cec37..0cbc779e7 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -2244,4 +2244,4 @@ @article { richardson_chill_1974 doi = "10.21273/HORTSCI.9.4.331", pages= "331 - 332", url = "https://journals.ashs.org/hortsci/view/journals/hortsci/9/4/article-p331.xml" -} \ No newline at end of file +} From af5d5d16cf6ea3bb3cd03acd134f19cbe1bcc742 Mon Sep 17 00:00:00 2001 From: Pascal Bourgault Date: Tue, 1 Oct 2024 16:17:59 -0400 Subject: [PATCH 34/36] add french translations --- xclim/data/fr.json | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/xclim/data/fr.json b/xclim/data/fr.json index d11f72253..fe19a9cef 100644 --- a/xclim/data/fr.json +++ b/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": { + "long_name": "Chill portions selon le modèle Dynamique", + "description": "Les « chill portions » sont une mesure du potentiel de débourrement de différentes récoltes.", + "abstract": "Les « chill portions » sont une mesure du potentiel de débourrement de différentes récoltes.", + "title": "Chill portions selon le modèle Dynamique." + }, + "CU": { + "title": "Unités de froid (chill) selon le modèle Utah", + "long_name": "Unités de froid selon le modèle Utah", + "description": "Les « chill units » sont une mesure du potentiel de débourrement de différentes récoltes basés sur le modèle Utah.", + "abstract": "Les « chill units » sont une mesure du potentiel de débourrement de différentes récoltes basés sur le modèle Utah." } } From b6eeec1d091b5fd2a09e983b5252b3302f886085 Mon Sep 17 00:00:00 2001 From: saschahofmann Date: Wed, 2 Oct 2024 10:53:43 +0200 Subject: [PATCH 35/36] Add todo for using resample_map --- xclim/indices/_agro.py | 1 + 1 file changed, 1 insertion(+) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 277d0dfec..d26b0a918 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1627,6 +1627,7 @@ def chill_portions( tas_K: xarray.DataArray = select_time( convert_units_to(tas, "K"), drop=True, **indexer ) + # TODO: use resample_map once #1848 is merged return ( tas_K.resample(time=freq) .map(_apply_chill_portion_one_season) From 427a3136ad88ad05e2f93ed4e668eb084f173693 Mon Sep 17 00:00:00 2001 From: Pascal Bourgault Date: Wed, 2 Oct 2024 09:22:25 -0400 Subject: [PATCH 36/36] Meilleures traductions --- xclim/data/fr.json | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/xclim/data/fr.json b/xclim/data/fr.json index fe19a9cef..ba7d3be7b 100644 --- a/xclim/data/fr.json +++ b/xclim/data/fr.json @@ -1387,15 +1387,15 @@ "abstract": "Calcul du rayonnement de grandes longueurs d'onde ascendant à partir du rayonnement de grandes longueurs d'onde net et descendant." }, "CP": { - "long_name": "Chill portions selon le modèle Dynamique", - "description": "Les « chill portions » sont une mesure du potentiel de débourrement de différentes récoltes.", - "abstract": "Les « chill portions » sont une mesure du potentiel de débourrement de différentes récoltes.", - "title": "Chill portions selon le modèle Dynamique." + "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) selon le modèle Utah", + "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 « chill units » sont une mesure du potentiel de débourrement de différentes récoltes basés sur le modèle Utah.", - "abstract": "Les « chill units » sont une mesure du potentiel de débourrement de différentes récoltes basés sur 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." } }