Degree day calculation unit handling errors

[tlogan2000]

Setup Information

• Xclim version: 0.47.0
• Python version:3.11
• Operating System:limux

Description

Calculation of growing degree days using different units does not yield the same results :

I believe that in cumulative_difference

xclim/xclim/indices/generic.py

Line 862 in a1a1a25

def cumulative_difference(

we should maybe convert the data to the threshold units instead of the opposite
doing this threshold = convert_units_to(threshold, data) means that output units should be set to 'F days' not 'K days' that I am getting in my output

Steps To Reproduce

import xarray as xr



import matplotlib.pyplot as plt
from clisops.core import subset
import xclim

ds = xr.tutorial.load_dataset("air_temperature", chunks={})
ds = ds.rename({'air':'tas'}).resample(time='D').mean(dim='time')
ds['tas_F'] = xclim.units.convert_units_to(ds.tas, "degF")
# Single point for the example
ds_pt = ds.sel(lon=-73, lat=44, method='nearest')

out1 = xclim.atmos.growing_degree_days(tas=ds_pt.tas, thresh="5 degC", freq="MS")

out2 = xclim.atmos.growing_degree_days(tas=ds_pt.tas_F, thresh="5 degC", freq="MS")

out3 = xclim.atmos.growing_degree_days(tas=ds_pt.tas_F, thresh="278.15 K", freq="MS")

# Plot figure & tadaa! Everything is the same
fig, ax = plt.subplots(figsize=(12, 5))
fig = out1.plot(label="K and degC", linestyle="-", ax=ax)
fig = out2.plot(label="degF and degC", marker="s", markersize=10, linestyle="none", ax=ax)
fig = out3.plot(label="degF and K", marker="o", linestyle="none", ax=ax)
fig = ax.legend()
ax.set_title("Different input units");

Additional context

No response

Contribution

• I would be willing/able to open a Pull Request to address this bug.

Code of Conduct

• I agree to follow this project's Code of Conduct

[Zeitsperre]

This merits a patch release once the fix is merged.

[tlogan2000]

@aulemahal The problem seems to be only with fahrenheit. and I think the fix needs to possibly take place in to_agg_units in the end:

The output gets assigned a cf units attr here

xclim/xclim/core/units.py

Line 555 in a1a1a25

out.attrs["units"] = pint2cfunits(orig_u * freq_u)

but from what I see the cf attrs are assigned directly without actually applying a conversion factor i,e,
pint2cfunits(1 <Unit('degree_Fahrenheit')> * 1 <Unit('day')>) gives 'K d' as an attr but values are still 'F d'

[aulemahal]

Indeed

import xclim as xc
a = xc.core.units.str2pint('1 °F')
b = xc.core.units.str2pint('1 d')

a * b
# 255.92777777777778 <Unit('day * kelvin')>

And this on pint 0.23, so not an issue of the latest version. I think the fix goes in to_agg_units.

EDIT: Actually that doesnt work at all. 1 °F d should be 0.5555 K d, not 255.

[aulemahal]

AH, I see the issue now. Summoning @coxipi with whom we discussed this thing before.

To avoid treating the "delta temperature" issue we decided to ignore the magnitude difference of the integral. 1 °C * 1 d = 274.15 K d which has an incorrect magnitude in the context of degree days, but the correct units, so we simply muted the magnitude part.

But 1 °F * 1 d = 256 K d for pint. It would work if it returned 459.67 R d, but it doesn't...

Previously we did (orig_u - orig_u) * freq_u to circumvent this problem. I guess we'll need something alike when a non-absolute scale is detected.

[tlogan2000]

Would it help to use pint.delta_{orig_unit} if thresh and data are temperature? https://pint.readthedocs.io/en/0.4.2/nonmult.html

[aulemahal]

Yes it would! I am currently looking for a way to elegantly convert non-mult temperatures to deltas in to_agg_units.