xclim.convert_units_to issues with accuracy

[Zeitsperre]

Setup Information

• Xclim version: v0.53.2
• Python version: Python 3.11
• Operating System: Linux

Description

On PAVICS I'm seeing a difference in the expected unit outputs when comparing Celsius, Kelvin and Fahrenheit

Steps To Reproduce

ds = xr.tutorial.open_dataset("air_temperature", chunks="auto")

# Xclim likes daily data
ds = ds.resample(time="D").mean(keep_attrs=True)
ds = ds.rename({"air": "tas"})

# A temperature dataset in Fahrenheit
ds["tas_F"] = xclim.units.convert_units_to(ds.tas, "degF")

ds_pt = subset.subset_gridpoint(ds, lon=-73, lat=44)

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 (not) the same(???)
fig, ax = plt.subplots(figsize=(12, 5))
out1.plot(label="K and degC", linestyle="-", ax=ax)
out2.plot(label="degF and degC", marker="s", markersize=10, linestyle="none", ax=ax)
out3.plot(label="degF and K", marker="o", linestyle="none", ax=ax)
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]

Is it possible that the outputs now reflect the units of the input array by default? i.e. degC days vs degF days? If so, this is a breaking behavioural change.

This appears to be the case. For out1, this is fixed by:

out1 = xclim.core.units.convert_units_to(out1, "degF days")

Update: This is very odd behaviour. This should be reporting in "K days" by default.

[tlogan2000]

Indicator definitions seem to clearly indicate output units of K days but seems to not get applied correctly? e.g.

xclim/xclim/indicators/atmos/_temperature.py

Line 814 in e4aca67

growing_degree_days = TempWithIndexing(

[tlogan2000]

I don't think we explicitly test Fahrenheit conversions very often (or even at all?)

[coxipi]

Oups, I realize that what you guys were saying, I'll leave the plot there

out1.plot(label="K and degC", linestyle="-", ax=ax)
f = convert_units_to("degC day", "degF day")
(out2/f).plot(label="degF and degC / conv-factor", marker="s", markersize=10, linestyle="none", ax=ax)
(out3/f).plot(label="degF and K / conv-factor", marker="o", linestyle="none", ax=ax)
ax.legend()
ax.set_title("Different input units");

I have a different gridpoint, so the plot is a bit different, but the idea is the same. By taking the conversion factor between degC d and degF d f = convert_units_to("degC day", "degF day"), the curves match

[Zeitsperre]

@coxipi

My issue here is that the initial calculations report that all three curves report having units K days. After converting the units as you mention, the resulting units for that curve become d degF. The fact that they all match in value is great, but their units don't correspond.

If you try going the other way around (f = convert_units_to("degF day", "degC day") they won't match either. Something is clearly wrong.

When converting out2 and out3 to degF:

[coxipi]

Yes, I understand.

The problem is already present in the index function

print(xclim.indices.growing_degree_days(tas=ds_pt.tas, thresh="5 degC", freq="MS").units)
print(xclim.indices.growing_degree_days(tas=ds_pt.tas_F, thresh="5 degC", freq="MS").units)

>>> 'K d'
>>> 'K d'

both give 'K d', but the magnitudes of the values are clearly 'K d' and 'degF d'

[coxipi]

print(cumulative_difference(ds_pt.tas, "5 degC", ">", "MS").units) # 'd K' expected
print(cumulative_difference(ds_pt.tas_F, "5 degC", ">", "MS").units) # 'd degF' expected?
>>> d K
>>> d K

cumulative_difference uses to_agg_units, which has the following line:

            else:
                out.attrs.update(pint2cfattrs(orig_u * freq_u, is_difference))

[aulemahal]

This is #1789. It was fixed in #1806, then broken again in #1836.

[Zeitsperre]

This was fixed in #2122