From 84776ee1b0f41a0826386af5c7cd780094b4fb9e Mon Sep 17 00:00:00 2001 From: Eddy Comyn-Platt <53045993+EddyCMWF@users.noreply.github.com> Date: Fri, 24 Nov 2023 17:43:24 +0000 Subject: [PATCH] heuristic methods added (#10) * heuristic methods added * qa --- earthkit/climate/heuristics.py | 149 +++++++++++++++++++++++++++++++++ tests/test_10_heuristics.py | 68 +++++++++++++++ 2 files changed, 217 insertions(+) create mode 100644 earthkit/climate/heuristics.py create mode 100644 tests/test_10_heuristics.py diff --git a/earthkit/climate/heuristics.py b/earthkit/climate/heuristics.py new file mode 100644 index 0000000..f466f8b --- /dev/null +++ b/earthkit/climate/heuristics.py @@ -0,0 +1,149 @@ +"""Heuristics methods for climate indicators.""" +import typing as T + +import numpy as np + + +def growing_degree_days( + tas_mean: T.Union[None, np.ndarray] = None, + tas_max: T.Union[None, np.ndarray] = None, + tas_min: T.Union[None, np.ndarray] = None, + tas_base: float = 0.0, + time_axis: int = 0, + time_start_index: int = 0, + time_stop_index: T.Union[None, int] = None, +) -> np.ndarray: + """Return the Growing Degree Days (GDD) index calculated from daily Near-Surface Air Temperature data. + + Growing degrees days (GDD) is defined as the mean daily temperature (`tas_mean`), or the average of daily + maximum and minimum temperatures (`tas_max` and `tas_min`), above a certain tas_base base temperature + (`tas_base`) accumulated on a daily basis over a period of time. https://doi.org/10.1002/joc.4535 + + Parameters + ---------- + tas_mean : numpy.ndarray + Daily mean Near-Surface Air Temperature. Required if tas_min and tax_max not provided. + tas_max : numpy.ndarray + Daily maximum Near-Surface Air Temperature. Required if tas_mean not provided. + tas_min : numpy.ndarray + Daily minimum Near-Surface Air Temperature. Must be same units and dimensions as tas_max. + Required if tas_mean not provided. + tas_base : float + Threshold temperature for the GDD index calculation, must be in same units as the tas array[s]. + time_axis : int (optional) + axis of the time dimension in the tas array, default is 0 + time_start_index : int (optional) + Start index to start accumulating degree days, default is 0 + time_stop_index : int (optional) + Stop index to start accumulating degree days, default is None, i.e. end of array + + Returns + ------- + numpy.ndarray: Growing Degree Days index. + + """ + # Check we have either tas_mean, or tas_min and tax_max + l_tas_mean = tas_mean is not None + l_tas_minmax = tas_min is not None and tas_max is not None + + if l_tas_mean and not l_tas_minmax: + gdu = tas_mean - tas_base + elif l_tas_minmax and not l_tas_mean: + gdu = (tas_max + tas_min) / 2 - tas_base + elif l_tas_mean and l_tas_minmax: + raise AssertionError("Please provide tas_mean OR tas_min+tas_max.") + else: + raise AssertionError("Insufficient arguments provided.") + + gdu[gdu < 0.0] = 0.0 + _gdu = np.rollaxis(gdu, time_axis) + + gdd = np.zeros_like(gdu) + _gdd = np.rollaxis(gdd, time_axis) + + if time_stop_index is not None: + _gdd[time_start_index:time_stop_index] = np.cumsum(_gdu[time_start_index:time_stop_index], axis=0) + else: + _gdd[time_start_index:] = np.cumsum(_gdu[time_start_index:], axis=0) + + return gdd + + +def heating_degree_days( + tas_max: np.ndarray, + tas_mean: np.ndarray, + tas_min: np.ndarray, + tas_base: float, +) -> np.ndarray: + """Return daily heating degree days. + + The heating degree days are computed using Spinoni's method which is used by the European + Environmental Agency (EEA). The method is described in the following + `paper (Spinoni et al., 2018) `_. + + Parameters + ---------- + tas_max : numpy.ndarray + Daily maximum Near-Surface Air Temperature. Must be same units and dimensions as tas_mean and tas_min. + tas_mean : numpy.ndarray + Daily mean Near-Surface Air Temperature. Must be same units and dimensions as tas_max and tas_min. + tas_min : numpy.ndarray + Daily minimum Near-Surface Air Temperature. Must be same units and dimensions as tas_max and tas_mean. + tas_base: float + Threshold temperature for the HDD index calculation, must be in same units as the tas arrays. + + Returns + ------- + numpy.ndarray: Daily heating degree days. + + + """ + hdd = np.zeros_like(tas_mean) + _mask = tas_base >= tas_max + hdd[_mask] = tas_base - tas_mean[_mask] + _mask = (tas_base < tas_max) & (tas_base >= tas_mean) + hdd[_mask] = (tas_base - tas_min[_mask]) / 2 - (tas_max[_mask] - tas_base) / 4 + _mask = (tas_base < tas_mean) & (tas_base >= tas_min) + hdd[_mask] = (tas_base - tas_min[_mask]) / 4 + + return hdd + + +def cooling_degree_days( + tas_max: np.ndarray, + tas_mean: np.ndarray, + tas_min: np.ndarray, + tas_base: float, +) -> np.ndarray: + """Return daily cooling degree days. + + The cooling degree days are computed using Spinoni's method which is used by the European + Environmental Agency (EEA). The method is described in the following + `paper (Spinoni et al., 2018) `_. + + Parameters + ---------- + tas_max : numpy.ndarray + Daily maximum Near-Surface Air Temperature. Must be same units and dimensions as tas_mean and tas_min. + tas_mean : numpy.ndarray + Daily mean Near-Surface Air Temperature. Must be same units and dimensions as tas_max and tas_min. + tas_min : numpy.ndarray + Daily minimum Near-Surface Air Temperature. Must be same units and dimensions as tas_max and tas_mean. + tas_base: float + Threshold temperature for the CDD index calculation, must be in same units as the tas arrays. + + Returns + ------- + numpy.ndarray: Daily cooling degree days. + + + """ + cdd = np.zeros_like(tas_mean) + _mask = tas_base < tas_max + cdd[_mask] = tas_mean[_mask] - tas_base + _mask = (tas_base < tas_max) & (tas_base >= tas_mean) + cdd[_mask] = (tas_max[_mask] - tas_base) / 4 + _mask = (tas_base < tas_mean) & (tas_base >= tas_min) + cdd[_mask] = (tas_max[_mask] - tas_base) / 2 - (tas_base - tas_min[_mask]) / 4 + + return cdd diff --git a/tests/test_10_heuristics.py b/tests/test_10_heuristics.py new file mode 100644 index 0000000..930cb24 --- /dev/null +++ b/tests/test_10_heuristics.py @@ -0,0 +1,68 @@ +import unittest + +import numpy as np + +from earthkit.climate import heuristics + + +class TestHeuristics(unittest.TestCase): + def test_growing_degree_days_tas_mean(self): + tas_mean = np.array([[20, 25, 30], [22, 27, 32]]) + tas_base = 15 + expected_result = np.array([[5, 10, 15], [12, 22, 32]]) + + result = heuristics.growing_degree_days(tas_mean=tas_mean, tas_base=tas_base) + + np.testing.assert_array_equal(result, expected_result) + + def test_growing_degree_days_tas_min_max(self): + tas_min = np.array([[18, 23, 28], [20, 25, 30]]) + tas_max = np.array([[22, 27, 32], [24, 29, 34]]) + tas_base = 15 + expected_result = np.array([[5.0, 10.0, 15.0], [12.0, 22.0, 32.0]]) + + result = heuristics.growing_degree_days(tas_min=tas_min, tas_max=tas_max, tas_base=tas_base) + + np.testing.assert_array_equal(result, expected_result) + + def test_growing_degree_days_invalid_input(self): + with self.assertRaises(AssertionError): + heuristics.growing_degree_days() + + def test_growing_degree_days_invalid_input_both(self): + with self.assertRaises(AssertionError): + heuristics.growing_degree_days( + tas_mean=np.array([[20, 25, 30]]), + tas_min=np.array([[18, 23, 28]]), + tas_max=np.array([[22, 27, 32]]), + ) + + def test_heating_degree_days(self): + tas_max = np.array([[10, 15, 20], [12, 18, 24]]) + tas_mean = np.array([[5, 10, 15], [8, 15, 20]]) + tas_min = np.array([[2, 8, 12], [6, 12, 18]]) + tas_base = 10 + expected_result = np.array([[5, 0, 0], [1, 0, 0]]) + + result = heuristics.heating_degree_days( + tas_max=tas_max, tas_mean=tas_mean, tas_min=tas_min, tas_base=tas_base + ) + + np.testing.assert_array_equal(result, expected_result) + + def test_cooling_degree_days(self): + tas_max = np.array([[25, 30, 35], [28, 33, 38]]) + tas_mean = np.array([[20, 25, 30], [22, 28, 34]]) + tas_min = np.array([[15, 20, 25], [18, 23, 28]]) + tas_base = 26 + expected_result = np.array([[0, 1, 4], [0, 2, 8]]) + + result = heuristics.cooling_degree_days( + tas_max=tas_max, tas_mean=tas_mean, tas_min=tas_min, tas_base=tas_base + ) + + np.testing.assert_array_equal(result, expected_result) + + +if __name__ == "__main__": + unittest.main()