From 06404e23c650e0aca98231daf0a276a8b1f756ec Mon Sep 17 00:00:00 2001
From: Pascal Bourgault <bourgault.pascal@ouranos.ca>
Date: Wed, 26 Jun 2024 11:09:06 -0400
Subject: [PATCH] Fix fahrenheit deltas

---
 CHANGES.rst         |  6 +++++-
 tests/test_units.py |  2 ++
 xclim/core/units.py | 30 +++++++++++++++++++++++++++---
 3 files changed, 34 insertions(+), 4 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 42ff6c252..3df3cc150 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -4,7 +4,11 @@ Changelog
 
 v0.51.0 (unreleased)
 --------------------
-Contributors to this version: Trevor James Smith (:user:`Zeitsperre`).
+Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`).
+
+Bug fixes
+^^^^^^^^^
+* Units of degree-days computations with Fahrenheit input fixed to yield "°R d". Added a new ``core.units.ensure_absolute_temperature``  method to convert from delta to absolute temperatures (:issue:`1789`, :pull:`1804`).
 
 Internal changes
 ^^^^^^^^^^^^^^^^
diff --git a/tests/test_units.py b/tests/test_units.py
index e94db9140..b59626fde 100644
--- a/tests/test_units.py
+++ b/tests/test_units.py
@@ -332,6 +332,8 @@ def index(
         ("", "sum", "count", 365, "d"),
         ("kg m-2", "var", "var", 0, "kg2 m-4"),
         ("°C", "argmax", "doymax", 0, ""),
+        ("°C", "sum", "integral", 365, "K d"),
+        ("°F", "sum", "integral", 365, "d °R"),  # not sure why the order is different
     ],
 )
 def test_to_agg_units(in_u, opfunc, op, exp, exp_u):
diff --git a/xclim/core/units.py b/xclim/core/units.py
index f06e00040..38a5b5d76 100644
--- a/xclim/core/units.py
+++ b/xclim/core/units.py
@@ -463,6 +463,25 @@ def infer_sampling_units(
     return out
 
 
+DELTA_ABSOLUTE_TEMP = {
+    units.delta_degC: units.kelvin,
+    units.delta_degF: units.rankine,
+}
+
+
+def ensure_absolute_temperature(units: str):
+    """Convert temperature units to their absolute counterpart, assuming they represented a difference (delta).
+
+    Celsius becomes Kelvin, Fahrenheit becomes Rankine. Does nothing for other units.
+    """
+    a = str2pint(units)
+    # ensure a delta pint unit
+    a = a - 0 * a
+    if a.units in DELTA_ABSOLUTE_TEMP:
+        return pint2cfunits(DELTA_ABSOLUTE_TEMP[a.units])
+    return units
+
+
 def to_agg_units(
     out: xr.DataArray, orig: xr.DataArray, op: str, dim: str = "time"
 ) -> xr.DataArray:
@@ -526,11 +545,16 @@ def to_agg_units(
     >>> degdays.units
     'K d'
     """
-    if op in ["amin", "min", "amax", "max", "mean", "std", "sum"]:
+    if op in ["amin", "min", "amax", "max", "mean", "sum"]:
         out.attrs["units"] = orig.attrs["units"]
 
+    elif op in ["std"]:
+        out.attrs["units"] = ensure_absolute_temperature(orig.attrs["units"])
+
     elif op in ["var"]:
-        out.attrs["units"] = pint2cfunits(str2pint(orig.units) ** 2)
+        out.attrs["units"] = pint2cfunits(
+            str2pint(ensure_absolute_temperature(orig.units)) ** 2
+        )
 
     elif op in ["doymin", "doymax"]:
         out.attrs.update(
@@ -539,7 +563,7 @@ def to_agg_units(
 
     elif op in ["count", "integral"]:
         m, freq_u_raw = infer_sampling_units(orig[dim])
-        orig_u = str2pint(orig.units)
+        orig_u = str2pint(ensure_absolute_temperature(orig.units))
         freq_u = str2pint(freq_u_raw)
         out = out * m