Bias removal with XClim methods

[gmfarache]

Hello,

I'm trying to execute the bias removal with Climpred using Quantile Delta Mapping algorithm (wrapped from XClim package) but the result is not being what I expected.

I have a model dataset which has two spatial coordinates (latitude, longitude) and two temporal coordinates (init, lead). Furthermore, I also have an observed dataset which has two spatial coordinates (latitude, longitude) and one temporal coordinate (time). The spatial and temporal coordinates of both datasets are aligned correctly.
Model dataset:

Observed dataset:

My objective using this bias removal component of Climpred is to input the whole historical data that I have of a forecast model and observed data and them get from the output, in other words the bias corrected dataset, only the current forecast init, which is the last init of the model dataset (in this case init='2024-06-13').

When I execute the bias removal with an algorithm implemented in Climpred such as additive_mean the result is ok, the resulting dataset has all the init and leads from the model dataset with its values corrected, as you can see in the following figure

In other words, executing this process with a additive_mean algorithm the bias corrected dataset has the last init (2024-06-13) which is of interest to me. But when I execute with an algorithm wrapped from XClim or from bias_correction packages such as QuantileDeltaMapping the bias removed dataset doesn't have the last init, it was just dropped from the dataset for a reason that I still don't know why.

Bias removed dataset:

In this case I want the result of bias corrected model from init='2024-06-13' but it is missing for some reason.

Can someone help me on that?

Obs: I'm using alignment='maximize' and train_test_split='unfair'.

[aaronspring]

Thanks for using climpred.

Really unsure why this is happening. I would have expected the same dimensions for xclim wrapped as for the additive mean.

Can you provide sample code to reproduce based on the sample datasets from the docs?

[gmfarache]

Yes, I was expecting the same dimensions for both xclim and additive mean results too.
I wrote a code to test this problem with the sample datasets from the docs, and it doesn't occur. Here it is:

import climpred
from climpred import HindcastEnsemble

initialized = climpred.tutorial.load_dataset("NMME_hindcast_Nino34_sst")
obs = climpred.tutorial.load_dataset("NMME_OIv2_Nino34_sst")
v = "sst"

initialized = initialized.sel(S=slice("1982-01-01", "2020-12-01"))
obs = obs.sel(time=slice("1982-01-01", "2020-12-01"))

hindcast = HindcastEnsemble(
    initialized.sel(model="GFDL-CM2p5-FLOR-A06", M=1.0).drop_vars(["model", "M"]),
).add_observations(obs)

am_res = hindcast.remove_bias(
    how="additive_mean",
    alignment="maximize",
    train_test_split="unfair",
).get_initialized().sortby('init')

obs.sst.attrs = {'units': 'degC'}
initialized.sst.attrs = {'units': 'degC'}

hindcast = HindcastEnsemble(
    initialized.sel(model="GFDL-CM2p5-FLOR-A06", M=1.0).drop_vars(["model", "M"]),
).add_observations(obs)

qdm_res = hindcast.remove_bias(
    how="QuantileDeltaMapping",
    alignment="maximize",
    train_test_split="unfair",
).get_initialized().sortby('init')

print(am_res, qdm_res)

In this case, both qdm_res and am_res datasets has 468 time values, which is the same dimensions as the input datasets. Do you have any clue what could be happening with the datasets that I'm trying to remove the bias?

[aaronspring]

Thanks for the excellent example code despite this not reproducing your error. The code for remove_bias isn't very readable unfortunately.

Things to check for debugging:

• Is every init unique in your dataset?
• try again but only with 100 inits and check. Likely now the latest init is close or equal to the latest time and maybe in the case something gets dropped.
• try the same only with 1-2 leads

[gmfarache]

Thank you for the recomendations, here it goes.
Yes, every init is unique in the dataset. If I try to remove the bias only with 1-2 leads, the problem is the same. And, if I try with only 100 inits it returns the following error:

>>> hindcast.remove_bias(how='QuantileDeltaMapping', alignment='maximize', train_test_split="unfair").get_initialized().sortby('init')
Traceback (most recent call last):
  File "/usr/local/lib/python3.8/site-packages/xclim/sdba/base.py", line 675, in _call_and_transpose_on_exit
    out = func(dsblock, **kwargs).transpose(*all_dims)
  File "/usr/local/lib/python3.8/site-packages/xclim/sdba/_adjustment.py", line 162, in qdm_adjust
    af = u.interp_on_quantiles(
  File "<boltons.funcutils.FunctionBuilder-263>", line 2, in interp_on_quantiles
  File "/usr/local/lib/python3.8/site-packages/xclim/sdba/base.py", line 476, in _parse_group
    return func(*args, **kwargs)
  File "/usr/local/lib/python3.8/site-packages/xclim/sdba/utils.py", line 443, in interp_on_quantiles
    return xr.apply_ufunc(
  File "/usr/local/lib/python3.8/site-packages/xarray/core/computation.py", line 1208, in apply_ufunc
    return apply_dataarray_vfunc(
  File "/usr/local/lib/python3.8/site-packages/xarray/core/computation.py", line 299, in apply_dataarray_vfunc
    args = deep_align(
  File "/usr/local/lib/python3.8/site-packages/xarray/core/alignment.py", line 860, in deep_align
    aligned = align(
  File "/usr/local/lib/python3.8/site-packages/xarray/core/alignment.py", line 797, in align
    aligner.align()
  File "/usr/local/lib/python3.8/site-packages/xarray/core/alignment.py", line 584, in align
    self.align_indexes()
  File "/usr/local/lib/python3.8/site-packages/xarray/core/alignment.py", line 429, in align_indexes
    raise ValueError(
ValueError: cannot align objects with join='exact' where index/labels/sizes are not equal along these coordinates (dimensions): 'month' ('month',)

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/local/lib/python3.8/site-packages/climpred/classes.py", line 2829, in remove_bias
    self = func(
  File "/usr/local/lib/python3.8/site-packages/climpred/bias_removal.py", line 731, in xclim_sdba
    bias_removed_hind = hindcast.verify(
  File "/usr/local/lib/python3.8/site-packages/climpred/classes.py", line 2377, in verify
    res = self._apply_climpred_function(
  File "/usr/local/lib/python3.8/site-packages/climpred/classes.py", line 1845, in _apply_climpred_function
    return func(hind.drop_vars(drop_init), verif.drop_vars(drop_obs), **kwargs)
  File "/usr/local/lib/python3.8/site-packages/climpred/classes.py", line 2292, in _verify
    metric_over_leads = [
  File "/usr/local/lib/python3.8/site-packages/climpred/classes.py", line 2293, in <listcomp>
    _apply_metric_at_given_lead(
  File "/usr/local/lib/python3.8/site-packages/climpred/prediction.py", line 100, in _apply_metric_at_given_lead
    result = metric.function(lforecast, lverif, dim=dim, **metric_kwargs)
  File "/usr/local/lib/python3.8/site-packages/climpred/bias_removal.py", line 714, in bc_func
    c[v] = adjustment(reference[v], model[v], data_to_be_corrected[v])
  File "/usr/local/lib/python3.8/site-packages/climpred/bias_removal.py", line 707, in adjustment
    data_to_be_corrected = dqm.adjust(data_to_be_corrected, **adjust_kwargs)
  File "/usr/local/lib/python3.8/site-packages/xclim/sdba/adjustment.py", line 216, in adjust
    out = self._adjust(sim, *args, **kwargs)
  File "/usr/local/lib/python3.8/site-packages/xclim/sdba/adjustment.py", line 557, in _adjust
    out = qdm_adjust(
  File "<boltons.funcutils.FunctionBuilder-285>", line 2, in _map_blocks
  File "/usr/local/lib/python3.8/site-packages/xclim/sdba/base.py", line 476, in _parse_group
    return func(*args, **kwargs)
  File "/usr/local/lib/python3.8/site-packages/xclim/sdba/base.py", line 694, in _map_blocks
    out = ds.map_blocks(
  File "/usr/local/lib/python3.8/site-packages/xarray/core/dataset.py", line 7731, in map_blocks
    return map_blocks(func, self, args, kwargs, template)
  File "/usr/local/lib/python3.8/site-packages/xarray/core/parallel.py", line 327, in map_blocks
    return func(obj, *args, **kwargs)
  File "/usr/local/lib/python3.8/site-packages/xclim/sdba/base.py", line 677, in _call_and_transpose_on_exit
    raise ValueError(
ValueError: qdm_adjust failed on block with coords : Coordinates:
  * latitude   (latitude) float64 4.75 4.5 4.25 4.0 ... -59.5 -59.75 -60.0
  * longitude  (longitude) float64 270.0 270.2 270.5 270.8 ... 329.5 329.8 330.0
  * time       (time) object 2023-06-04 00:00:00 ... 2023-09-11 00:00:00
  * quantiles  (quantiles) float32 0.025 0.075 0.125 0.175 ... 0.875 0.925 0.975
  * month      (month) int64 6 7 8 9.