Index function temperature_sum (closes #81)

climix/etc/climate_indices_CV22.yml

@@ -1466,7 +1466,7 @@ indices:
       standard_name: integral_of_air_temperature_deficit_wrt_time
       proposed_standard_name: _
       long_name: Heating Degree Days (Tmean < {TT}C)
-      units: degree_Celsius
+      units: degree_Celsius days
       cell_methods:
         - time: mean within days
         - time: sum over days

climix/etc/metadata.yml

@@ -88,3 +88,18 @@ index_functions:
         kind: operator
       reducer:
         kind: reducer
+
+  temperature_sum:
+    description: |
+      Calculates the temperature sum above/below a threshold. First, the threshold
+      is transformed to the same standard_name and units as the input data.
+      Then the thresholding is performed as condition(data, threshold), ie
+      if condition is <, data < threshold.
+      Finally, the sum is calculated for those data values that fulfill
+      the condition after subtraction of the threshold value. If the sum is for
+      values below the threshold the result is multiplied by -1.
+    parameters:
+      threshold:
+        kind: quantity
+      condition:
+        kind: operator

climix/index_functions.py

@@ -289,3 +289,37 @@ class ThresholdedStatistics:
         comb = self.lazy_condition(data, self.threshold.points)
         res = self.lazy_reducer(da.ma.masked_where(~comb, data), axis=axis)
         return res.astype('float32')
+
+
+class TemperatureSum:
+    def __init__(self, threshold, condition):
+        self.threshold = threshold
+        if condition in ['>', '>=']:
+            self.fun = lambda d, t: np.maximum(d - t, 0)
+            self.lazy_fun = lambda d, t: da.maximum(d - t, 0)
+        else:
+            self.fun = lambda d, t: np.maximum(t - d, 0)
+            self.lazy_fun = lambda d, t: da.maximum(t - d, 0)
+        self.extra_coords = [threshold]
+
+    def prepare(self, input_cube):
+        change_units(self.threshold,
+                     input_cube.units,
+                     input_cube.standard_name)
+        self.standard_name = input_cube.standard_name
+        if input_cube.units.is_convertible('degC'):
+            self.units = 'degC days'
+        else:
+            raise RuntimeError("Invalid input units")
+
+    def call_func(self, data, axis, **kwargs):
+        axis = normalize_axis(axis, data.ndim)
+        threshold = self.threshold.points[0]
+        res = np.sum(self.fun(data, threshold), axis=axis)
+        return res.astype('float32')
+
+    def lazy_func(self, data, axis, **kwargs):
+        axis = normalize_axis(axis, data.ndim)
+        threshold = self.threshold.points[0]
+        res = da.sum(self.lazy_fun(data, threshold), axis=axis)
+        return res.astype('float32')

setup.py

@@ -57,6 +57,7 @@ setuptools.setup(
             'spell_length=climix.index_functions:SpellLength',
             'statistics=climix.index_functions:Statistics',
             'thresholded_statistics=climix.index_functions:ThresholdedStatistics',
+            'temperature_sum=climix.index_functions:TemperatureSum',
         ],
     },
     project_urls={