Linear interpolation for seasonal bias adjustment leads to non-smooth distribution

[saschahofmann]

Setup Information

• Xclim version: 0.52.2
• Python version: 3.11.6
• Operating System: Linux

Description

If I change the grouper from the first example of the docs to use time.season instead of time.month. I get some drastic transitions when the seasons change:

I remember we did additions to enable the interpolation for these string columns of seasons but now I am not sure its working correctly. I can see that it changes the coordinates to integers and adds the cyclic bounds so I am not sure why this would be happening?

Steps To Reproduce

Use the example from the tutorial or paste the following

import matplotlib.pyplot as plt
import nc_time_axis  # noqa
import numpy as np
import xarray as xr

%matplotlib inline
plt.style.use("seaborn-v0_8")
plt.rcParams["figure.figsize"] = (11, 5)

# Create toy data to explore bias adjustment, here fake temperature timeseries
t = xr.cftime_range("2000-01-01", "2030-12-31", freq="D", calendar="noleap")
ref = xr.DataArray(
    (
        -20 * np.cos(2 * np.pi * t.dayofyear / 365)
        + 2 * np.random.random_sample((t.size,))
        + 273.15
        + 0.1 * (t - t[0]).days / 365
    ),  # "warming" of 1K per decade,
    dims=("time",),
    coords={"time": t},
    attrs={"units": "K"},
)
sim = xr.DataArray(
    (
        -18 * np.cos(2 * np.pi * t.dayofyear / 365)
        + 2 * np.random.random_sample((t.size,))
        + 273.15
        + 0.11 * (t - t[0]).days / 365
    ),  # "warming" of 1.1K per decade
    dims=("time",),
    coords={"time": t},
    attrs={"units": "K"},
)

ref = ref.sel(time=slice(None, "2015-01-01"))
hist = sim.sel(time=slice(None, "2015-01-01"))
QM_mo = sdba.QuantileDeltaMapping.train(
    ref, hist, nquantiles=15, group="time.season", kind="+"
)
scen = QM_mo.adjust(sim, extrapolation="constant", interp="linear")

ref.groupby("time.dayofyear").mean().plot(label="Reference")
hist.groupby("time.dayofyear").mean().plot(label="Model - biased")
scen.sel(time=slice("2000", "2015")).groupby("time.dayofyear").mean().plot(
    label="Model - adjusted - 2000-15", linestyle="--"
)
scen.sel(time=slice("2015", "2030")).groupby("time.dayofyear").mean().plot(
    label="Model - adjusted - 2015-30", linestyle="--"
)
plt.legend()

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

[saschahofmann]

Ok I think I am getting closer to the problem. If I understand it correctly, newx in interp_on_quantiles should be a float of the season but instead it still seems to be an int. Looking even closer this is the relevant part of the get_index function of the grouper class (note the line for self.prop=='season'):

ind = da.indexes[self.dim]
        if self.prop == "week":
            i = da[self.dim].copy(data=ind.isocalendar().week).astype(int)
        elif self.prop == "season":
            i = da[self.dim].copy(data=ind.month % 12 // 3)
        else:
            i = getattr(ind, self.prop)

        if not np.issubdtype(i.dtype, np.integer):
            raise ValueError(
                f"Index {self.name} is not of type int (rather {i.dtype}), "
                f"but {self.__class__.__name__} requires integer indexes."
            )

        if interp and self.dim == "time" and self.prop == "month":
            i = ind.month - 0.5 + ind.day / ind.days_in_month

so I think the problem is that at no point are the actual seasonal weights computed.

I'll try to work out how this should like and will report back.