Support for temperature differences in CF attributes (#1836)

Set ``units_metadata`` attribute to ``"temperature: difference"`` for
arrays representing a difference between two temperatures.

<!--Please ensure the PR fulfills the following requirements! -->
<!-- If this is your first PR, make sure to add your details to the
AUTHORS.rst! -->
### Pull Request Checklist:
- [x] This PR addresses an already opened issue (for bug fixes /
features)
    - This PR fixes #1822
- [x] Tests for the changes have been added (for bug fixes / features)
- [x] (If applicable) Documentation has been added / updated (for bug
fixes / features)
- [x] CHANGELOG.rst has been updated (with summary of main changes)
- [x] Link to issue (:issue:`number`) and pull request (:pull:`number`)
has been added

### What kind of change does this PR introduce?

* Add ``units_metadata`` attribute set to ``"temperature: difference"``
whenever values represent a difference between two temperatures (degree
days, diurnal amplitude, bias, etc).
* Adjust the units handling functions to recognize this attribute. 
* Modify indices that set units to `delta_degC`, which is not CF
compliant, by the original input DataArray units and the
``units_metadata`` attribute.
* Create new function ``pint2cfattrs`` to convert delta units into a
dictionary of attributes, including ``units_metadata``.


### Does this PR introduce a breaking change?

Yes, units returned by some indices will change.

### Other information:

CHANGELOG.rst

@@ -22,6 +22,7 @@ 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`).
+* Add attribute ``units_metadata`` to outputs representing a difference between temperatures. This is needed to disambiguate temperature differences from absolute temperature. Changes affect indicators ``daily_temperature_range``, ``daily_temperature_range_variability``, ``extreme_temperature_range``, ``interday_diurnal_temperature_range``, and all degree-day indicators. Implemented using a new ``pint2cfattrs`` function to convert pint units to a dictionary of CF attributes. ``units2pint`` is also modified to support ``units_metadata`` attributes in DataArrays. Some SDBA properties and measures previously returning units of ``delta_degC`` will now return the original input DataArray units accompanied with the ``units_metadata`` attribute. (:issue:`1822`, :pull:`1830`).
 * ``xclim.indices.run_length.windowed_max_run_sum`` accumulates positive values across runs and yields the the maximum valued run. (:pull:`1926`).
 * Helper function ``xclim.indices.helpers.make_hourly_temperature`` to estimate hourly temperatures from daily min and max temperatures. (:pull:`1909`).
 * New global option ``resample_map_blocks`` to wrap all ``resample().map()`` code inside a ``xr.map_blocks`` to lower the number of dask tasks. Uses utility ``xclim.indices.helpers.resample_map`` and requires ``flox`` to ensure the chunking allows such block-mapping. Defaults to False. (:pull:`1848`).

docs/notebooks/units.ipynb

@@ -220,6 +220,48 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
+    "## Temperature differences vs absolute temperature\n",
+    "\n",
+    "Temperature anomalies and biases as well as degree-days indicators are all *differences* between temperatures. If we assign those differences units of degrees Celsius, then converting to Kelvins or Fahrenheits will yield nonsensical values. ``pint`` solves this using *delta* units such as ``delta_degC`` and ``delta_degF``. \n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# If we have a DataArray storing a temperature anomaly of 2°C, converting to Kelvin will yield a nonsensical value 0f 275.15.\n",
+    "# Fortunately, pint has delta units to represent temperature differences.\n",
+    "display(units.convert_units_to(xr.DataArray([2], attrs={\"units\": \"delta_degC\"}), \"K\"))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The issue for ``xclim`` is that there are no equivalent delta units in the CF convention. To resolve this ambiguity, the [CF convention](https://cfconventions.org/Data/cf-conventions/cf-conventions-1.11/cf-conventions.html#temperature-units) recommends including a ``units_metadata`` attribute set to ``\"temperature: difference\"``, and this is supported in ``xclim`` as of version 0.52. The function ``units2pint`` interprets the ``units_metadata`` attribute and returns a ``pint`` delta unit as needed. To convert a ``pint`` delta unit to CF attributes, use the function ``pint2cfattrs``, which returns a dictionary with the ``units`` and ``units_metadata`` attributes (``pint2cfunits`` cannot support the convention because it only returns the unit string).   "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "delta = xr.DataArray(\n",
+    "    [2], attrs={\"units\": \"K\", \"units_metadata\": \"temperature: difference\"}\n",
+    ")\n",
+    "units.convert_units_to(delta, \"delta_degC\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "\n",
+    "\n",
+    "\n",
     "## Other utilities\n",
     "\n",
     "Many helper functions are defined in `xclim.core.units`, see [Unit handling module](../api.rst#units-handling-submodule).\n"

tests/test_generic.py

@@ -334,6 +334,12 @@ def test_forbidden(self, tas_series):
         with pytest.raises(NotImplementedError):
             generic.cumulative_difference(tas, threshold="10 degC", op="!=")
 
+    def test_delta_units(self, tas_series):
+        tas = tas_series(np.array([-10, 15, 20, 3, 10]) + K2C)
+
+        out = generic.cumulative_difference(tas, threshold="10 degC", op=">=")
+        assert "units_metadata" in out.attrs
+
 
 class TestFirstDayThreshold:
     @pytest.mark.parametrize(

tests/test_indices.py

@@ -2115,6 +2115,7 @@ def test_static_daily_temperature_range(self, tasmin_series, tasmax_series):
         tasmin, tasmax = self.static_tmin_tmax_setup(tasmin_series, tasmax_series)
         dtr = xci.daily_temperature_range(tasmin, tasmax, freq="YS")
         assert dtr.units == "K"
+        assert dtr.units_metadata == "temperature: difference"
         output = np.mean(tasmax - tasmin)
 
         np.testing.assert_equal(dtr, output)
@@ -2134,12 +2135,14 @@ def test_static_variable_daily_temperature_range(
         dtr = xci.daily_temperature_range_variability(tasmin, tasmax, freq="YS")
 
         np.testing.assert_almost_equal(dtr, 2.667, decimal=3)
+        assert dtr.units_metadata == "temperature: difference"
 
     def test_static_extreme_temperature_range(self, tasmin_series, tasmax_series):
         tasmin, tasmax = self.static_tmin_tmax_setup(tasmin_series, tasmax_series)
         etr = xci.extreme_temperature_range(tasmin, tasmax)
 
         np.testing.assert_array_almost_equal(etr, 31.7)
+        assert etr.units_metadata == "temperature: difference"
 
     def test_uniform_freeze_thaw_cycles(self, tasmin_series, tasmax_series):
         temp_values = np.zeros(365)

tests/test_sdba/test_properties.py

@@ -287,7 +287,8 @@ def test_annual_cycle(self, open_dataset):
 
         assert amp.long_name.startswith("Absolute amplitude of the annual cycle")
         assert phase.long_name.startswith("Phase of the annual cycle")
-        assert amp.units == "delta_degC"
+        assert amp.units == "K"
+        assert amp.units_metadata == "temperature: difference"
         assert relamp.units == "%"
         assert phase.units == ""
 
@@ -330,7 +331,8 @@ def test_annual_range(self, open_dataset):
 
         assert amp.long_name.startswith("Average annual absolute amplitude")
         assert phase.long_name.startswith("Average annual phase")
-        assert amp.units == "delta_degC"
+        assert amp.units == "K"
+        assert amp.units_metadata == "temperature: difference"
         assert relamp.units == "%"
         assert phase.units == ""
 

tests/test_units.py

@@ -7,8 +7,6 @@
 import pytest
 import xarray as xr
 from dask import array as dsk
-from packaging.version import Version
-from pint import __version__ as __pint_version__
 
 from xclim import indices, set_options
 from xclim.core import Quantified, ValidationError
@@ -21,6 +19,7 @@
     declare_units,
     infer_context,
     lwethickness2amount,
+    pint2cfattrs,
     pint2cfunits,
     pint_multiply,
     rate2amount,
@@ -119,6 +118,14 @@ def test_cf_conversion_amount2lwethickness_amount2rate(self):
         assert out.attrs["units"] == "kg d-1 m-2"  # CF equivalent unit
         assert out.attrs["standard_name"] == "rainfall_flux"
 
+    def test_temperature_difference(self):
+        delta = xr.DataArray(
+            [2], attrs={"units": "K", "units_metadata": "temperature: difference"}
+        )
+        out = convert_units_to(source=delta, target="delta_degC")
+        assert out == 2
+        assert out.attrs["units"] == "degC"
+
 
 class TestUnitConversion:
     def test_pint2cfunits(self):
@@ -343,9 +350,9 @@ def index(
             "sum",
             "integral",
             365,
-            ("K d", "d K"),
+            ("degC d", "d degC"),
         ),  # dependent on numpy/pint version
-        ("°F", "sum", "integral", 365, "d °R"),  # not sure why the order is different
+        ("°F", "sum", "integral", 365, "d degF"),  # not sure why the order is different
     ],
 )
 def test_to_agg_units(in_u, opfunc, op, exp, exp_u):
@@ -355,15 +362,37 @@ def test_to_agg_units(in_u, opfunc, op, exp, exp_u):
         coords={"time": xr.cftime_range("1993-01-01", periods=365, freq="D")},
         attrs={"units": in_u},
     )
+    if units(in_u).dimensionality == "[temperature]":
+        da.attrs["units_metadata"] = "temperature: difference"
 
     # FIXME: This is emitting warnings from deprecated DataArray.argmax() usage.
     out = to_agg_units(getattr(da, opfunc)(), da, op)
     np.testing.assert_allclose(out, exp)
-
     if isinstance(exp_u, tuple):
-        if Version(__pint_version__) < Version("0.24.1"):
-            assert out.attrs["units"] == exp_u[0]
-        else:
-            assert out.attrs["units"] == exp_u[1]
+        assert out.attrs["units"] in exp_u
     else:
         assert out.attrs["units"] == exp_u
+
+
+def test_pint2cfattrs():
+    attrs = pint2cfattrs(units.degK, is_difference=True)
+    assert attrs == {"units": "K", "units_metadata": "temperature: difference"}
+
+    attrs = pint2cfattrs(units.meter, is_difference=True)
+    assert "units_metadata" not in attrs
+
+    attrs = pint2cfattrs(units.delta_degC)
+    assert attrs == {"units": "degC", "units_metadata": "temperature: difference"}
+
+
+def test_temp_difference_rountrip():
+    """Test roundtrip of temperature difference units."""
+    attrs = {"units": "degC", "units_metadata": "temperature: difference"}
+    da = xr.DataArray([1], attrs=attrs)
+    pu = units2pint(da)
+    # Confirm that we get delta pint units
+    assert pu == units.delta_degC
+
+    # and that converting those back to cf attrs gives the same result
+    attrs = pint2cfattrs(pu)
+    assert attrs == {"units": "degC", "units_metadata": "temperature: difference"}

xclim/core/indicator.py

@@ -345,6 +345,7 @@ class Indicator(IndicatorRegistrar):
         "standard_name",
         "long_name",
         "units",
+        "units_metadata",
         "cell_methods",
         "description",
         "comment",
@@ -836,7 +837,10 @@ def __call__(self, *args, **kwds):
         for out, attrs, base_attrs in zip(outs, out_attrs, self.cf_attrs, strict=False):
             if self.n_outs > 1:
                 var_id = base_attrs["var_name"]
+            # Set default units
             attrs.update(units=out.units)
+            if "units_metadata" in out.attrs:
+                attrs["units_metadata"] = out.attrs["units_metadata"]
             attrs.update(
                 self._update_attrs(
                     params.copy(),