diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 868e9acc7..c529fcc05 100644 --- a/CHANGELOG.rst +++ b/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 ^^^^^^^^^ diff --git a/docs/references.bib b/docs/references.bib index 9f2419dd9..0cbc779e7 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" +} diff --git a/tests/test_atmos.py b/tests/test_atmos.py index 23929550d..80757e5f7 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,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) 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/tests/test_indices.py b/tests/test_indices.py index d1da9c54c..7269931c4 100644 --- a/tests/test_indices.py +++ b/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")) diff --git a/xclim/data/fr.json b/xclim/data/fr.json index d11f72253..ba7d3be7b 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": { + "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." } } diff --git a/xclim/indicators/atmos/_temperature.py b/xclim/indicators/atmos/_temperature.py index d1efd311a..46d6903ee 100644 --- a/xclim/indicators/atmos/_temperature.py +++ b/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,6 +107,12 @@ 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.""" @@ -107,6 +120,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", @@ -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, +) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 38f939c29..d26b0a918 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -35,6 +35,8 @@ __all__ = [ "biologically_effective_degree_days", + "chill_portions", + "chill_units", "cool_night_index", "corn_heat_units", "dryness_index", @@ -1524,3 +1526,159 @@ def hardiness_zones( zones = zones.assign_attrs(units="") return zones + + +def _accumulate_intermediate(prev_E, prev_xi, curr_xs, curr_ak1): + """Accumulate the intermediate product based on the previous concentration and the current temperature.""" + curr_S = np.where(prev_E < 1, prev_E, prev_E - prev_E * prev_xi) + return curr_xs - (curr_xs - curr_S) * np.exp(-curr_ak1) + + +def _chill_portion_one_season(tas_K): + """Computes the chill portion for a single season based on the dynamic model on a numpy array.""" + # Constants as described in Luedeling et al. (2009) + E0 = 4153.5 + E1 = 12888.8 + A0 = 139500 + A1 = 2.567e18 + SLP = 1.6 + TETMLT = 277 + AA = A0 / A1 + EE = E1 - E0 + + ftmprt = SLP * TETMLT * (tas_K - TETMLT) / tas_K + sr = np.exp(ftmprt) + xi = sr / (1 + sr) + xs = AA * np.exp(EE / tas_K) + ak1 = A1 * np.exp(-E1 / tas_K) + + inter_E = np.zeros_like(tas_K) + for i in range(1, tas_K.shape[-1]): + inter_E[..., i] = _accumulate_intermediate( + inter_E[..., i - 1], xi[..., i - 1], xs[..., i], ak1[..., i] + ) + delta = np.where(inter_E >= 1, inter_E * xi, 0) + + return delta + + +def _apply_chill_portion_one_season(tas_K): + """Apply the chill portion function on to an xarray DataArray.""" + tas_K = tas_K.chunk(time=-1) + return xarray.apply_ufunc( + _chill_portion_one_season, + tas_K, + input_core_dims=[["time"]], + output_core_dims=[["time"]], + output_dtypes=[tas_K.dtype], + dask="parallelized", + ).sum("time") + + +@declare_units(tas="[temperature]") +def chill_portions( + tas: xarray.DataArray, freq: str = "YS", **indexer +) -> xarray.DataArray: + """Chill portion based on the dynamic model + + Chill portions are a measure to estimate the bud breaking potential of different crop. + The constants and functions are taken from Luedeling et al. (2009) which formalises + the method described in Fishman et al. (1987). The model computes the accumulation of + an intermediate product that is transformed to the final product once it exceeds a + certain concentration. The intermediate product can be broken down at higher temperatures + but the final product is stable even at higher temperature. Thus the dynamic model is + more accurate than the Utah model especially in moderate climates like Israel, + California or Spain. + + Parameters + ---------- + tas : xr.DataArray + Hourly temperature. + + Returns + ------- + xr.DataArray, [unitless] + Chill portions after the Dynamic Model + + 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 + ) + # TODO: use resample_map once #1848 is merged + return ( + tas_K.resample(time=freq) + .map(_apply_chill_portion_one_season) + .assign_attrs(units="") + ) + + +@declare_units(tas="[temperature]") +def chill_units(tas: xarray.DataArray, freq: str = "YS") -> xarray.DataArray: + """Chill units using the Utah model + + Chill units are a measure to estimate the bud breaking potential of different crop based on Richardson et al. (1974). + The Utah model assigns a weight to each hour depending on the temperature recognising that high temperatures can actual decrease, + the potential for bud breaking. + + Parameters + ---------- + tas : xr.DataArray + Hourly temperature. + + Returns + ------- + xr.DataArray, [unitless] + Chill units using the Utah model + + 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( + (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), + ), + ), + ) + cu = cu.where(tas.notnull()) + return cu.resample(time=freq).sum().assign_attrs(units="") diff --git a/xclim/indices/helpers.py b/xclim/indices/helpers.py index 2a172212c..f4943d97b 100644 --- a/xclim/indices/helpers.py +++ b/xclim/indices/helpers.py @@ -8,6 +8,7 @@ from __future__ import annotations from collections.abc import Mapping +from datetime import timedelta from inspect import stack from typing import Any, cast @@ -559,3 +560,153 @@ 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( + (np.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 _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 to 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. + + 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 it's used for daily aggregation, we assume that sunrise globally happens at midnight + and the sunsets after `daylength` hours computed via the :py:func:`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=_add_one_day(data.isel(time=-1).time) + ), + ], + 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), + ) + # 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 = nighttime_hours = ( + hourly.time.dt.hour + 1 - hourly.daylength + ).clip(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, + ), + ).assign_attrs(units=tasmin.units)