Repository: cerfacs-globc/icclim

Issue:
BUG: Unit handling for temperature difference indices
<!--
If you have encountered a problem while using icclim or you would like to suggest a change,
please follow the structure below to describe your request.
-->
<!-- import icclim; print(icclim.__version__) -->
<!-- import sys; print(sys.version) -->
* icclim version: 6.1.3
* Python version: 3.10.8


### Description
It appears that for indices that involve calculating the difference between two temperatures (e.g. DTR, ETR, VDTR), unit conversion is performed after the difference is taken. This is problematic if the input units are Kelvin, because the difference between two temperatures in Kelvin is equivalent to the difference in Celsius and hence no unit conversion is required.


### Minimal reproducible example
[Example 2: Index ETR](https://icclim.readthedocs.io/en/stable/how_to/recipes_ecad.html#example-2-index-etr)
```python
files_tasmax = [
    "tasmax_day_CNRM-CM5_historical_r1i1p1_19300101-19341231.nc",
    "tasmax_day_CNRM-CM5_historical_r1i1p1_19350101-19391231.nc",
]
files_tasmin = [
    "tasmin_day_CNRM-CM5_historical_r1i1p1_19300101-19341231.nc",
    "tasmin_day_CNRM-CM5_historical_r1i1p1_19350101-19391231.nc",
]

out_f = "ETR_year_CNRM-CM5_historical_r1i1p1_1930-1939.nc"  # OUTPUT FILE: annual values of ETR

icclim.index(
    index_name="ETR",
    in_files=[files_tasmax, files_tasmin],
    var_name=["tasmax", "tasmin"],
    slice_mode="year",
    out_file=out_f,
)
```

### Output received
The output (in `out_f`) is a data array (with units = C) whose values range from -270 to -182 because unit conversion (i.e. subtracting 273.15 to move from K to C) appears to have been performed after the difference in temperatures (in Kelvin) was calculated.

```python
import xarray as xr

ds = xr.open_dataset('ETR_year_CNRM-CM5_histrorical_r1i1p1_1930-1939.nc')

data_min = ds['ETR'].values.min()
data_max = ds['ETR'].values.max()
data_units = ds['ETR'].attrs['units']

print('Data minimum:', data_min)
print('Data maximum:', data_max)
print('Data units:', data_units)
```

```
Data minimum: -270.46213
Data maximum: -182.72186
Data units: C
```

Discussion:

Patch:
diff --git a/doc/source/_static/logo_icclim_colored__displayed.svg b/doc/source/_static/logo_icclim_colored__displayed.svg
index c3d52e8..de59981 100644
--- a/doc/source/_static/logo_icclim_colored__displayed.svg
+++ b/doc/source/_static/logo_icclim_colored__displayed.svg
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:cbb3b0dc6978cb3a9b3d72fc9f27d13ceb2f593ac081bb9ff8f687535dd923e5
+oid sha256:02be51c905baba036d5fbc1d622eb8509ee6fa11ac40209e30f4fa28ca72add6
 size 66296
diff --git a/doc/source/_static/logo_icclim_grey__displayed.svg b/doc/source/_static/logo_icclim_grey__displayed.svg
index ffff686..ff05fb3 100644
--- a/doc/source/_static/logo_icclim_grey__displayed.svg
+++ b/doc/source/_static/logo_icclim_grey__displayed.svg
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:12c01ff2f5def2cc8c950f4f822b869adceadec85f1011ae5769dff33429094a
+oid sha256:4e786640ea61f8e0353a78c0db676e891ad081c4074451951e396ab48234b1f6
 size 63678
diff --git a/doc/source/_static/logo_icclim_white__displayed.svg b/doc/source/_static/logo_icclim_white__displayed.svg
index 1b05f43..732124a 100644
--- a/doc/source/_static/logo_icclim_white__displayed.svg
+++ b/doc/source/_static/logo_icclim_white__displayed.svg
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:c78c30348dce8a5c484e4fcd697fc9b4c00816daa80789e42fb93416bc47a1b4
+oid sha256:041391f619b6a855034ff99d0659d172bbb8f626c2cec7e23351dc7344c2ed59
 size 39021
diff --git a/doc/source/references/release_notes.rst b/doc/source/references/release_notes.rst
index 2e949f7..45b6cbf 100644
--- a/doc/source/references/release_notes.rst
+++ b/doc/source/references/release_notes.rst
@@ -1,6 +1,12 @@
 Release history
 ===============
 
+6.3.0 (unreleased)
+------------------
+* [fix] **BREAKING CHANGE** The indicators based on the difference between two variables (ecad: DTR, ETR, vDTR and anomaly) gave wrong values due to a bad unit conversion of the output.
+  This was for example the case when the input variables are in Kelvin, the difference between the two variables is still in Kelvin but it cannot be converted to degree Celsius with the ususal `+273.15`.
+  To workaround this issue, we first convert inputs to the expected output unit and then we compute the index.
+
 6.2.0
 -----
 * [maint] Upgrade and adapt to xclim 0.40.
diff --git a/icclim/_generated_api.py b/icclim/_generated_api.py
index 55a9fd4..5410b1c 100644
--- a/icclim/_generated_api.py
+++ b/icclim/_generated_api.py
@@ -6,7 +6,7 @@ This module exposes each climate index as individual functions for convenience.
 # flake8: noqa E501
 from __future__ import annotations
 
-import datetime
+from datetime import datetime
 from typing import Sequence
 
 from xarray.core.dataset import Dataset
diff --git a/icclim/generic_indices/generic_indicators.py b/icclim/generic_indices/generic_indicators.py
index 4a07351..ed40b63 100644
--- a/icclim/generic_indices/generic_indicators.py
+++ b/icclim/generic_indices/generic_indicators.py
@@ -213,7 +213,6 @@ class GenericIndicator(ResamplingIndicator):
         output_unit: str | None,
         coef: float | None,
         sampling_method: str,
-        is_compared_to_reference: bool,
     ) -> list[ClimateVariable]:
         if not _same_freq_for_all(climate_vars):
             raise InvalidIcclimArgumentError(
@@ -227,14 +226,24 @@ class GenericIndicator(ResamplingIndicator):
                 f"{self.name} can only be computed with the following"
                 f" sampling_method(s): {self.sampling_methods}"
             )
-        if output_unit is not None and _is_amount_unit(output_unit):
-            for climate_var in climate_vars:
-                current_unit = climate_var.studied_data.attrs.get(UNITS_KEY, None)
-                if current_unit is not None and not _is_amount_unit(current_unit):
-                    with xc_units.context("hydro"):
-                        climate_var.studied_data = rate2amount(
-                            climate_var.studied_data, out_units=output_unit
-                        )
+        if output_unit is not None:
+            if _is_amount_unit(output_unit):
+                climate_vars = _convert_rates_to_amounts(
+                    climate_vars=climate_vars, output_unit=output_unit
+                )
+            elif _is_a_diff_indicator(self) and output_unit != "%":
+                # [gh:255] Indicators computing the difference between two
+                # variables must first convert the units of input variables
+                # to the expected output unit in order to avoid converting
+                # the output which would be a relative "temperature".
+                # In other words: a 15 Kelvin difference *is* equivalent
+                # to a 15 degC and *is not* to a -258.15 degC.
+                for climate_var in climate_vars:
+                    climate_var.studied_data = convert_units_to(
+                        climate_var.studied_data, target=output_unit
+                    )
+            else:
+                pass  # nothing to do
         if coef is not None:
             for climate_var in climate_vars:
                 climate_var.studied_data = coef * climate_var.studied_data
@@ -277,7 +286,6 @@ class GenericIndicator(ResamplingIndicator):
             output_unit=config.out_unit,
             coef=config.coef,
             sampling_method=config.sampling_method,
-            is_compared_to_reference=config.is_compared_to_reference,
         )
         result = self.process(
             climate_vars=climate_vars,
@@ -723,9 +731,10 @@ def _check_couple_of_vars(
     if len(climate_vars) != 2:
         raise InvalidIcclimArgumentError(
             f"{indicator.name} can only be computed on two variables sharing the same"
-            f" unit (e.g. 2 temperatures). Either provide a `base_period_time_range` to"
-            f" create a reference variable or directly provide a secondary variable"
-            f" with `in_files` or `var_name`."
+            f" unit (e.g. 2 temperatures). You can either provide a secondary variable"
+            f" with `in_files` or `var_name`, or you can let icclim compute this"
+            f" second variable as a subset of the first one using"
+            f" `base_period_time_range`."
         )
 
 
@@ -1164,6 +1173,26 @@ def _is_amount_unit(unit: str) -> bool:
         return False
 
 
+def _is_a_diff_indicator(indicator: Indicator) -> bool:
+    return (
+        indicator == GenericIndicatorRegistry.DifferenceOfExtremes
+        or indicator == GenericIndicatorRegistry.MeanOfDifference
+        or indicator == GenericIndicatorRegistry.MeanOfAbsoluteOneTimeStepDifference
+        or indicator == GenericIndicatorRegistry.DifferenceOfMeans
+    )
+
+
+def _convert_rates_to_amounts(climate_vars: list[ClimateVariable], output_unit: str):
+    for climate_var in climate_vars:
+        current_unit = climate_var.studied_data.attrs.get(UNITS_KEY, None)
+        if current_unit is not None and not _is_amount_unit(current_unit):
+            with xc_units.context("hydro"):
+                climate_var.studied_data = rate2amount(
+                    climate_var.studied_data, out_units=output_unit
+                )
+    return climate_vars
+
+
 def _to_percent(da: DataArray, sampling_freq: Frequency) -> DataArray:
     if sampling_freq == FrequencyRegistry.MONTH:
         da = da / da.time.dt.daysinmonth * 100
diff --git a/icclim/models/climate_variable.py b/icclim/models/climate_variable.py
index a42b048..512f5f1 100644
--- a/icclim/models/climate_variable.py
+++ b/icclim/models/climate_variable.py
@@ -32,9 +32,10 @@ from icclim.pre_processing.input_parsing import (
 
 @dataclass
 class ClimateVariable:
-    """Internal icclim structure. It groups together the input variable (studied_data),
-    its associated metadata (standard_variable) and the threshold it must be compared
-    to.
+    """Internal icclim structure.
+    It groups together the input variable (studied_data),
+    its associated metadata (standard_var) and, if any,
+    the threshold it must be compared to.
 
     Attributes
     ----------
diff --git a/icclim/models/constants.py b/icclim/models/constants.py
index a2d683c..f2de4ff 100644
--- a/icclim/models/constants.py
+++ b/icclim/models/constants.py
@@ -3,7 +3,7 @@ from __future__ import annotations
 # fmt: off
 # flake8: noqa
 
-ICCLIM_VERSION = "6.2.0"
+ICCLIM_VERSION = "6.3.0"
 
 # placeholders for user_index
 PERCENTILE_THRESHOLD_STAMP = "p"
diff --git a/setup.py b/setup.py
index 292b466..ab7a3f3 100644
--- a/setup.py
+++ b/setup.py
@@ -21,7 +21,7 @@ MINIMAL_REQUIREMENTS = [
 
 setup(
     name="icclim",
-    version="6.2.0",
+    version="6.3.0",
     packages=find_packages(),
     author="Christian P.",
     author_email="christian.page@cerfacs.fr",
diff --git a/tools/extract-icclim-funs.py b/tools/extract-icclim-funs.py
index 528b7c4..34b49b2 100644
--- a/tools/extract-icclim-funs.py
+++ b/tools/extract-icclim-funs.py
@@ -87,7 +87,7 @@ This module exposes each climate index as individual functions for convenience.
 # flake8: noqa E501
 from __future__ import annotations
 
-import datetime
+from datetime import datetime
 from typing import Sequence
 
 from xarray.core.dataset import Dataset


Test patch:
diff --git a/icclim/tests/test_main.py b/icclim/tests/test_main.py
index 2e99137..6f590f0 100644
--- a/icclim/tests/test_main.py
+++ b/icclim/tests/test_main.py
@@ -130,7 +130,7 @@ class Test_Integration:
 
     def test_index_DTR(self):
         ds = self.data.to_dataset(name="toto")
-        ds["tutu"] = self.data
+        ds["tutu"] = self.data + 10
         res = icclim.index(
             index_name="DTR",
             in_files=ds,
@@ -138,7 +138,21 @@ class Test_Integration:
             var_name=["toto", "tutu"],
         )
         assert f"icclim version: {ICCLIM_VERSION}" in res.attrs["history"]
-        np.testing.assert_array_equal(0, res.DTR)
+        np.testing.assert_array_equal(-10, res.DTR)
+
+    def test_index_DTR__with_unit_conversion(self):
+        ds = self.data.to_dataset(name="toto")
+        ds["tutu"] = self.data + 10
+        ds["toto"].attrs["units"] = "K"
+        ds["tutu"].attrs["units"] = "K"
+        res = icclim.dtr(
+            in_files=ds,
+            out_file=self.OUTPUT_FILE,
+            var_name=["toto", "tutu"],
+        )
+        assert f"icclim version: {ICCLIM_VERSION}" in res.attrs["history"]
+        np.testing.assert_array_equal(-10, res.DTR)
+        np.testing.assert_array_equal("°C", res.DTR.attrs["units"])
 
     def test_index_CD(self):
         ds = self.data.to_dataset(name="tas")