Repository: Ouranosinc/xclim

Issue:
DetrendedQuantileMapping only applies the last correction factor
* xclim version: '0.16.1-beta' (master, as of 2020-05-08)
* Python version: 3.8.2
* Operating System: Linux

### Description
The figure below is the daily average : `plt.plot(data.groupby("time.dayofyear").mean("time").values)`. The black line is my OBS data, the blue line is my RAW data, the green line is the result of EmpiricalQQMapping and the orange line is the result of DetrendedQuantileMapping.

![image](https://user-images.githubusercontent.com/38501935/81433867-41b59400-9133-11ea-8adf-b3be4e38fa67.png)

These are the options I used to create the orange line:
```
group_month_nowindow = xclim.sdba.utils.Grouper('time.month')
QQ_m_nw_ni = xclim.sdba.adjustment.DetrendedQuantileMapping(nquantiles=50, kind='+', group=group_month_nowindow, extrapolation='constant')
QQ_m_nw_ni.train(obs, hist)
qq_m_nw_ni = QQ_m_nw_ni.adjust(hist, degree=3)
```

And this is `QQ_m_nw_ni.ds.af.plot.contourf()`
![image](https://user-images.githubusercontent.com/38501935/81434207-d3250600-9133-11ea-9934-0f0d2e31771a.png)

The huge jump in May leads me to believe that the correction factor associated to the last quantile is what is being applied for the whole distribution.

### Additional Information
In `xclim.sdba.adjustment` and `_adjust` (line 230), I added the line with `+269` (which is roughly the annual average) :

```
sim = apply_correction(sim,
            broadcast(self.ds.scaling, sim, group=self.group, interp=self.interp),
            self.kind,
        )
        sim_fit = PolyDetrend(degree=degree, kind=self.kind).fit(sim)
        sim_detrended = sim_fit.detrend(sim)
        sim_detrended.values = sim_detrended.values + 269
        scen_detrended = super()._adjust(sim_detrended)
        scen = sim_fit.retrend(scen_detrended)
        scen.values = scen.values - 269  # I need to remove that +269
```
With that additional line, I don't have the bug. This leads me to believe that this is a discrepancy between `sim`, which is centered around 269°K, and `sim_detrended`, which is normalized and centered around 0.




Discussion:
coveralls: ## Pull Request Test Coverage Report for [Build 1851](https://coveralls.io/builds/30735548)

* **0** of **18**   **(0.0%)**  changed or added relevant lines in **1** file are covered.
* **2** unchanged lines in **1** file lost coverage.
* Overall coverage decreased (**-0.3%**) to **75.732%**

---

|  Changes Missing Coverage | Covered Lines | Changed/Added Lines | % |
| :-----|--------------|--------|---: |
| [xclim/sdba/adjustment.py](https://coveralls.io/builds/30735548/source?filename=xclim%2Fsdba%2Fadjustment.py#L12) | 0 | 18 | 0.0%
<!-- | **Total:** | **0** | **18** | **0.0%** | -->

|  Files with Coverage Reduction | New Missed Lines | % |
| :-----|--------------|--: |
| [xclim/sdba/adjustment.py](https://coveralls.io/builds/30735548/source?filename=xclim%2Fsdba%2Fadjustment.py#L258) | 2 | 0% |
<!-- | **Total:** | **2** |  | -->

|  Totals | [![Coverage Status](https://coveralls.io/builds/30735548/badge)](https://coveralls.io/builds/30735548) |
| :-- | --: |
| Change from base [Build 1849](https://coveralls.io/builds/30734731): |  -0.3% |
| Covered Lines: | 2381 |
| Relevant Lines: | 3144 |

---
##### 💛  - [Coveralls](https://coveralls.io)

aulemahal: Merged fixes for QDM in this PR to save time. Dask-arrays needed a new version of `rank` so I replicated the xarray one but with an `apply_ufunc` that allows for dask parallelization. 
aulemahal: My bad, the fix to `xclim.sdba.utils.broadcast` did not fix  #449 for all cases, so I am reverting it for now.
Also, the new rank util partially fixes #448. However, it is not a great use as long as 449 is not fixed...

Patch:
diff --git a/docs/sdba.rst b/docs/sdba.rst
index d25e699e..1ac5f555 100644
--- a/docs/sdba.rst
+++ b/docs/sdba.rst
@@ -2,18 +2,18 @@
 Bias adjustment and downscaling algorithms
 ==========================================
 
-`xarray` data structures allow for relatively straightforward implementations of simple bias-adjustment and downscaling algorithms documented in :ref:`bias-adjustment-algos`. Each algorithm is split into `train` and `adjust` components. The `train` function will compare two DataArrays `x` and `y`, and return a `Dataset` object storing the *transfer* information allowing to go from `x` to `y`. This dataset can then be used as an input to the `adjust` function to apply this information to `x`. `x` could be the same `DataArray` used for training, or another `DataArray` with similar characteristics.
+`xarray` data structures allow for relatively straightforward implementations of simple bias-adjustment and downscaling algorithms documented in :ref:`bias-adjustment-algos`. Each algorithm is split into `train` and `adjust` components. The `train` function will compare two DataArrays `x` and `y`, and create a dataset storing the *transfer* information allowing to go from `x` to `y`. This dataset, stored in the adjustment object, can then be used by the `adjust` method to apply this information to `x`. `x` could be the same `DataArray` used for training, or another `DataArray` with similar characteristics.
 
 For example, given a daily time series of observations `ref`, a model simulation over the observational period `hist` and a model simulation over a future period `sim`, we would apply a bias-adjustment method such as *detrended quantile mapping* (DQM) as::
 
   from xclim.sdba.adjustment import DetrendedQuantileMapping
   dqm = DetrendedQuantileMapping()
   dqm.train(ref, hist)
-  scen = dqm.predict(sim)
+  scen = dqm.adjust(sim)
 
 Most method can either be applied additively or multiplicatively. Also, most methods can be applied independently on different time groupings (monthly, seasonally) or according to the day of the year and a rolling window width.
 
-When transfer factors are applied in adjustment, they can be interpolated according to the time grouping. This helps avoid discontinuities in adjustment factors at the beginning of each season or month and is computationaly cheaper than computing adjustment factors for each day of the year.
+When transfer factors are applied in adjustment, they can be interpolated according to the time grouping. This helps avoid discontinuities in adjustment factors at the beginning of each season or month and is computationaly cheaper than computing adjustment factors for each day of the year. (Currently only implemented for monthly grouping)
 
 
 Application in multivariate settings
diff --git a/xclim/sdba/__init__.py b/xclim/sdba/__init__.py
index ac3a24d1..4c1d70d7 100644
--- a/xclim/sdba/__init__.py
+++ b/xclim/sdba/__init__.py
@@ -9,9 +9,9 @@ Adjustment algorithms all conform to the `train` - `adjust` scheme, formalized w
 Given a reference time series (ref), historical simulations (hist) and simulations to be adjusted (sim),
 any bias-adjustment method would be applied by first estimating the adjustment factors between the historical simulation and the observations series, and then applying these factors to `sim`, which could be a future simulation::
 
-  Adj = Adjustment(group="time.month", interp="linear")
+  Adj = Adjustment(group="time.month")
   Adj.train(ref, hist)
-  scen = Adj.adjust(sim)
+  scen = Adj.adjust(sim, interp="linear")
   Adj.ds.af  # adjustment factors.
 
 The `group` argument allows adjustment factors to be estimated independently for different periods: the full
@@ -37,13 +37,17 @@ A generic bias adjustment process is to layed out like:
 - preprocessing on `ref`, `hist` and `sim` (using methods in `xclim.sdba.processing` or `xclim.sdba.detrending`)
 - creating the adjustment object `Adj = Adjustment(**kwargs)` (from `xclim.sdba.adjustment`)
 - training `Adj.train(obs, sim)`
-- adjustment `scen = Adj.adjust(sim)`
+- adjustment with corresponding arguments `scen = Adj.adjust(sim, **kwargs)`
 - post-processing on `scen` (for example: re-trending)
 
-The train-adjust approach allows to inspect the trained adjustment object. The adjustment information is stored in
+The train-adjust approach allows to inspect the trained adjustment object. The trained information is stored in
 the underlying `Adj.ds` dataset and always has a `af` variable with the adjustment factors. Its layout and the
 other available variables vary between the different algorithm, refer to :ref:`bias-adjustment-algos`.
 
+Parameters needed by the training and the adjustment are saved to the `Adj.ds` dataset as a  `adj_params` attribute.
+Other parameters, those only needed by the adjustment are passed in the `adjust` call and written to the history attribute
+in the output scenario dataarray.
+
 Grouping
 ========
 
diff --git a/xclim/sdba/adjustment.py b/xclim/sdba/adjustment.py
index 60e96570..6aa1ccde 100644
--- a/xclim/sdba/adjustment.py
+++ b/xclim/sdba/adjustment.py
@@ -18,8 +18,10 @@ from .utils import equally_spaced_nodes
 from .utils import extrapolate_qm
 from .utils import get_correction
 from .utils import interp_on_quantiles
+from .utils import invert
 from .utils import map_cdf
 from .utils import MULTIPLICATIVE
+from .utils import rank
 from xclim.core.calendar import get_calendar
 from xclim.core.formatting import update_history
 
@@ -73,6 +75,8 @@ class BaseAdjustment(Parametrizable):
         ----------
         sim : DataArray
           Time series to be bias-adjusted, usually a model output.
+        kwargs :
+          Algorithm-specific keyword arguments, see class doc.
         """
         if not self.__trained:
             raise ValueError("train() must be called before adjusting.")
@@ -91,8 +95,9 @@ class BaseAdjustment(Parametrizable):
                 stacklevel=4,
             )
         scen = self._adjust(sim, **kwargs)
+        params = ", ".join([f"{k}={repr(v)}" for k, v in kwargs.items()])
         scen.attrs["history"] = update_history(
-            f"Bias-adjusted with method {str(self)}", sim
+            f"Bias-adjusted with {str(self)}.adjust(sim, {params})", sim
         )
         return scen
 
@@ -127,17 +132,22 @@ class EmpiricalQuantileMapping(BaseAdjustment):
 
     Parameters
     ----------
+    At init:
+
     nquantiles : int
       The number of quantiles to use. Two endpoints at 1e-6 and 1 - 1e-6 will be added.
     kind : {'+', '*'}
       The adjustment kind, either additive or multiplicative.
-    interp : {'nearest', 'linear', 'cubic'}
-      The interpolation method to use then interpolating the adjustment factors.
-    extrapolation : {'constant', 'nan'}
-      The type of extrapolation to use. See :py:func:`xclim.sdba.utils.extrapolate_qm` for details.
     group : Union[str, Grouper]
       The grouping information. See :py:class:`xclim.sdba.base.Grouper` for details.
 
+    In adjust:
+
+    interp : {'nearest', 'linear', 'cubic'}
+      The interpolation method to use when interpolating the adjustment factors. Defaults to "nearset".
+    extrapolation : {'constant', 'nan'}
+      The type of extrapolation to use. See :py:func:`xclim.sdba.utils.extrapolate_qm` for details. Defaults to "constant".
+
     References
     ----------
     Dequé, M. (2007). Frequency of precipitation and temperature extremes over France in an anthropogenic scenario: Model results and statistical correction according to observed values. Global and Planetary Change, 57(1–2), 16–26. https://doi.org/10.1016/j.gloplacha.2006.11.030
@@ -149,16 +159,10 @@ class EmpiricalQuantileMapping(BaseAdjustment):
         *,
         nquantiles: int = 20,
         kind: str = ADDITIVE,
-        interp: str = "nearest",
-        extrapolation: str = "constant",
         group: Union[str, Grouper] = "time",
     ):
         super().__init__(
-            nquantiles=nquantiles,
-            kind=kind,
-            interp=interp,
-            extrapolation=extrapolation,
-            group=group,
+            nquantiles=nquantiles, kind=kind, group=group,
         )
 
     def _train(self, ref, hist):
@@ -182,11 +186,9 @@ class EmpiricalQuantileMapping(BaseAdjustment):
         )
         self._make_dataset(af=af, hist_q=hist_q)
 
-    def _adjust(self, sim):
-        af, hist_q = extrapolate_qm(
-            self.ds.af, self.ds.hist_q, method=self.extrapolation
-        )
-        af = interp_on_quantiles(sim, hist_q, af, group=self.group, method=self.interp)
+    def _adjust(self, sim, interp="nearest", extrapolation="constant"):
+        af, hist_q = extrapolate_qm(self.ds.af, self.ds.hist_q, method=extrapolation)
+        af = interp_on_quantiles(sim, hist_q, af, group=self.group, method=interp)
 
         return apply_correction(sim, af, self.kind)
 
@@ -199,9 +201,12 @@ class DetrendedQuantileMapping(EmpiricalQuantileMapping):
     1. A scaling factor that would make the mean of `hist` match the mean of `ref` is computed.
     2. `ref` and `hist` are normalized by removing the group-wise mean.
     3. Adjustment factors are computed between the quantiles of the normalized `ref` and `hist`.
-    4. `sim` is corrected by the scaling factor then detrended using a linear fit.
+    4. `sim` is corrected by the scaling factor, normalized by the group-wise mean and then detrended using a linear fit.*
     5. Values of detrended `sim` are matched to the corresponding quantiles of normalized `hist` and corrected accordingly.
-    6. The trend is put back on the result.
+    6. The group-wise mean and trend are put back on the result.*
+
+    * Steps 4 and 6 include a group-wise normalization to overcome the fact that detrending is not made group-wise.
+    A future release of xclim will have grouped detrending and not require this extra step any more.
 
     .. math::
 
@@ -210,6 +215,26 @@ class DetrendedQuantileMapping(EmpiricalQuantileMapping):
     where :math:`F` is the cumulative distribution function (CDF) and :math:`\overline{xyz}` is the linear trend of the data.
     This equation is valid for multiplicative adjustment. Based on the DQM method of [Cannon2015]_.
 
+    Parameters
+    ----------
+    At init:
+
+    nquantiles : int
+      The number of quantiles to use. Two endpoints at 1e-6 and 1 - 1e-6 will be added.
+    kind : {'+', '*'}
+      The adjustment kind, either additive or multiplicative.
+    group : Union[str, Grouper]
+      The grouping information. See :py:class:`xclim.sdba.base.Grouper` for details.
+
+    In adjust:
+
+    interp : {'nearest', 'linear', 'cubic'}
+      The interpolation method to use when interpolating the adjustment factors. Defaults to "nearest".
+    detrend : int or BaseDetrend instance
+      The method to use when detrending. If an int is passed, it is understood as a PolyDetrend (polynomial detrending) degree. Defaults to 1 (linear detrending)
+    extrapolation : {'constant', 'nan'}
+      The type of extrapolation to use. See :py:func:`xclim.sdba.utils.extrapolate_qm` for details. Defaults to "constant".
+
     References
     ----------
     .. [Cannon2015] Cannon, A. J., Sobie, S. R., & Murdock, T. Q. (2015). Bias correction of GCM precipitation by quantile mapping: How well do methods preserve changes in quantiles and extremes? Journal of Climate, 28(17), 6938–6959. https://doi.org/10.1175/JCLI-D-14-00754.1
@@ -228,16 +253,46 @@ class DetrendedQuantileMapping(EmpiricalQuantileMapping):
             description="Scaling factor making the mean of hist match the one of hist.",
         )
 
-    def _adjust(self, sim, degree=0):
+    def _adjust(self, sim, interp="nearest", extrapolation="constant", detrend=1):
+        # Apply preliminary scaling from obs to hist
         sim = apply_correction(
             sim,
-            broadcast(self.ds.scaling, sim, group=self.group, interp=self.interp),
+            broadcast(self.ds.scaling, sim, group=self.group, interp=interp),
             self.kind,
         )
-        sim_fit = PolyDetrend(degree=degree, kind=self.kind).fit(sim)
-        sim_detrended = sim_fit.detrend(sim)
-        scen_detrended = super()._adjust(sim_detrended)
-        scen = sim_fit.retrend(scen_detrended)
+
+        # Normalize sim group-wise
+        # This group-wise pre normalization + reapplication further down
+        # is to circumvent #442 (detrending is not groupwise)
+        mu_sim = self.group.apply("mean", sim)
+        sim_norm = apply_correction(
+            sim,
+            broadcast(
+                invert(mu_sim, kind=self.kind), sim, group=self.group, interp=interp,
+            ),
+            kind=self.kind,
+        )
+
+        # Find trend on sim (for debugging purpose, the detrending is overridable)
+        if isinstance(detrend, int):
+            detrend = PolyDetrend(degree=detrend, kind=self.kind)
+
+        sim_fit = detrend.fit(sim_norm)
+        sim_detrended = sim_fit.detrend(sim_norm)
+
+        # Adjust using `EmpiricalQuantileMapping.adjust`
+        scen_detrended = super()._adjust(
+            sim_detrended, extrapolation=extrapolation, interp=interp
+        )
+        # Retrend
+        scen_norm = sim_fit.retrend(scen_detrended)
+
+        # # Reapply-mean
+        scen = apply_correction(
+            scen_norm,
+            broadcast(mu_sim, sim, group=self.group, interp=interp),
+            kind=self.kind,
+        )
         return scen
 
 
@@ -254,17 +309,35 @@ class QuantileDeltaMapping(EmpiricalQuantileMapping):
     where :math:`F` is the cumulative distribution function (CDF). This equation is valid for multiplicative adjustment.
     The algorithm is based on the "QDM" method of [Cannon2015]_.
 
+    Parameters
+    ----------
+    At init:
+
+    nquantiles : int
+      The number of quantiles to use. Two endpoints at 1e-6 and 1 - 1e-6 will be added.
+    kind : {'+', '*'}
+      The adjustment kind, either additive or multiplicative.
+    group : Union[str, Grouper]
+      The grouping information. See :py:class:`xclim.sdba.base.Grouper` for details.
+
+    In adjust:
+
+    interp : {'nearest', 'linear', 'cubic'}
+      The interpolation method to use when interpolating the adjustment factors. Defaults to "nearest".
+    extrapolation : {'constant', 'nan'}
+      The type of extrapolation to use. See :py:func:`xclim.sdba.utils.extrapolate_qm` for details. Defaults to "constant".
+
     References
     ----------
     .. [Cannon2015] Cannon, A. J., Sobie, S. R., & Murdock, T. Q. (2015). Bias correction of GCM precipitation by quantile mapping: How well do methods preserve changes in quantiles and extremes? Journal of Climate, 28(17), 6938–6959. https://doi.org/10.1175/JCLI-D-14-00754.1
     """
 
-    def _adjust(self, sim):
-        af, _ = extrapolate_qm(self.ds.af, self.ds.hist_q, method=self.extrapolation)
+    def _adjust(self, sim, interp="nearest", extrapolation="constant"):
+        af, _ = extrapolate_qm(self.ds.af, self.ds.hist_q, method=extrapolation)
 
-        sim_q = self.group.apply(xr.DataArray.rank, sim, main_only=True, pct=True)
+        sim_q = self.group.apply(rank, sim, main_only=True, pct=True)
         sel = {"quantiles": sim_q}
-        af = broadcast(af, sim, group=self.group, interp=self.interp, sel=sel)
+        af = broadcast(af, sim, group=self.group, interp=interp, sel=sel)
 
         return apply_correction(sim, af, self.kind)
 
@@ -295,12 +368,17 @@ class LOCI(BaseAdjustment):
 
     Parameters
     ----------
+    At init:
+
     group : Union[str, Grouper]
       The grouping information. See :py:class:`xclim.sdba.base.Grouper` for details.
     thresh : float
       The threshold in `ref` above which the values are scaled.
+
+    In adjust:
+
     interp : {'nearest', 'linear', 'cubic'}
-      The interpolation method to use then interpolating the adjustment factors.
+      The interpolation method to use then interpolating the adjustment factors. Defaults to "linear".
 
     References
     ----------
@@ -308,14 +386,8 @@ class LOCI(BaseAdjustment):
     """
 
     @parse_group
-    def __init__(
-        self,
-        *,
-        group: Union[str, Grouper] = "time",
-        thresh: float = None,
-        interp: str = "linear",
-    ):
-        super().__init__(group=group, thresh=thresh, interp=interp)
+    def __init__(self, *, group: Union[str, Grouper] = "time", thresh: float = None):
+        super().__init__(group=group, thresh=thresh)
 
     def _train(self, ref, hist):
         s_thresh = map_cdf(hist, ref, self.thresh, group=self.group).isel(
@@ -335,9 +407,9 @@ class LOCI(BaseAdjustment):
         s_thresh.attrs.update(long_name="Threshold over modeled data")
         self._make_dataset(hist_thresh=s_thresh, ref_thresh=self.thresh, af=af)
 
-    def _adjust(self, sim):
-        sth = broadcast(self.ds.hist_thresh, sim, group=self.group, interp=self.interp)
-        factor = broadcast(self.ds.af, sim, group=self.group, interp=self.interp)
+    def _adjust(self, sim, interp="linear"):
+        sth = broadcast(self.ds.hist_thresh, sim, group=self.group, interp=interp)
+        factor = broadcast(self.ds.af, sim, group=self.group, interp=interp)
         with xr.set_options(keep_attrs=True):
             scen = (factor * (sim - sth) + self.ds.ref_thresh).clip(min=0)
         return scen
@@ -351,17 +423,22 @@ class Scaling(BaseAdjustment):
 
     Parameters
     ----------
+    At init:
+
     group : Union[str, Grouper]
       The grouping information. See :py:class:`xclim.sdba.base.Grouper` for details.
     kind : {'+', '*'}
       The adjustment kind, either additive or multiplicative.
+
+    In adjust:
+
     interp : {'nearest', 'linear', 'cubic'}
-      The interpolation method to use then interpolating the adjustment factors.
+      The interpolation method to use then interpolating the adjustment factors. Defaults to "nearest".
     """
 
     @parse_group
-    def __init__(self, *, group="time", kind=ADDITIVE, interp="nearest"):
-        super().__init__(group=group, kind=kind, interp=interp)
+    def __init__(self, *, group="time", kind=ADDITIVE):
+        super().__init__(group=group, kind=kind)
 
     def _train(self, ref, hist):
         mean_hist = self.group.apply("mean", hist)
@@ -370,6 +447,6 @@ class Scaling(BaseAdjustment):
         af.attrs.update(long_name="Scaling adjustment factors")
         self._make_dataset(af=af)
 
-    def _adjust(self, sim):
-        factor = broadcast(self.ds.af, sim, group=self.group, interp=self.interp)
+    def _adjust(self, sim, interp="nearest"):
+        factor = broadcast(self.ds.af, sim, group=self.group, interp=interp)
         return apply_correction(sim, factor, self.kind)
diff --git a/xclim/sdba/utils.py b/xclim/sdba/utils.py
index ad7d05fd..5bdc3ab7 100644
--- a/xclim/sdba/utils.py
+++ b/xclim/sdba/utils.py
@@ -4,6 +4,7 @@ from typing import Sequence
 from typing import Union
 from warnings import warn
 
+import bottleneck as bn
 import numpy as np
 import xarray as xr
 from boltons.funcutils import wraps
@@ -211,6 +212,7 @@ def broadcast(
                 warn(
                     "Broadcasting operations in multiple dimensions can only be done with linear and nearest-neighbor interpolation, not cubic. Using linear."
                 )
+
             grouped = grouped.interp(sel, method=interp)
 
         for var in sel.keys():
@@ -462,3 +464,47 @@ def interp_on_quantiles(
         dask="parallelized",
         output_dtypes=[np.float],
     )
+
+
+def rank(da, dim="time", pct=False):
+    """Ranks data.
+
+    Replicates `xr.DataArray.rank` but with support for dask-stored data. Xarray's docstring is below:
+
+    Equal values are assigned a rank that is the average of the ranks that
+    would have been otherwise assigned to all of the values within that
+    set.  Ranks begin at 1, not 0. If pct, computes percentage ranks.
+
+    NaNs in the input array are returned as NaNs.
+
+    The `bottleneck` library is required.
+
+    Parameters
+    ----------
+    dim : hashable
+        Dimension over which to compute rank.
+    pct : bool, optional
+        If True, compute percentage ranks, otherwise compute integer ranks.
+
+    Returns
+    -------
+    ranked : DataArray
+        DataArray with the same coordinates and dtype 'float64'.
+    """
+
+    def _nanrank(data):
+        func = bn.nanrankdata if data.dtype.kind == "f" else bn.rankdata
+        ranked = func(data, axis=-1)
+        if pct:
+            count = np.sum(~np.isnan(data), axis=-1, keepdims=True)
+            ranked /= count
+        return ranked
+
+    return xr.apply_ufunc(
+        _nanrank,
+        da,
+        input_core_dims=[[dim]],
+        output_core_dims=[[dim]],
+        dask="parallelized",
+        output_dtypes=[da.dtype],
+    )


Test patch:
diff --git a/tests/test_sdba/test_adjustment.py b/tests/test_sdba/test_adjustment.py
index 2b802763..ac54ae4e 100644
--- a/tests/test_sdba/test_adjustment.py
+++ b/tests/test_sdba/test_adjustment.py
@@ -1,5 +1,6 @@
 import numpy as np
 import pytest
+import xarray as xr
 from scipy.stats import norm
 from scipy.stats import uniform
 
@@ -40,7 +41,7 @@ class TestLoci:
         np.testing.assert_array_almost_equal(p, ref, dec)
 
         assert "history" in p.attrs
-        assert "Bias-adjusted with method" in p.attrs["history"]
+        assert "Bias-adjusted with LOCI(" in p.attrs["history"]
 
 
 class TestScaling:
@@ -108,11 +109,9 @@ class TestDQM:
         hist = sim = series(x, name)
         ref = series(y, name)
 
-        DQM = DetrendedQuantileMapping(
-            kind=kind, group="time", nquantiles=50, interp="linear"
-        )
+        DQM = DetrendedQuantileMapping(kind=kind, group="time", nquantiles=50,)
         DQM.train(ref, hist)
-        p = DQM.adjust(sim)
+        p = DQM.adjust(sim, interp="linear")
 
         q = DQM.ds.quantiles
         ex = apply_correction(xd.ppf(q), invert(xd.mean(), kind), kind)
@@ -132,7 +131,7 @@ class TestDQM:
         sim2 = apply_correction(sim, ff, kind)
         ref2 = apply_correction(ref, ff, kind)
 
-        p2 = DQM.adjust(sim2)
+        p2 = DQM.adjust(sim2, interp="linear")
 
         np.testing.assert_array_almost_equal(p2[middle], ref2[middle], 1)
 
@@ -142,7 +141,7 @@ class TestDQM:
         )
         sim3 = apply_correction(sim, trend, kind)
         ref3 = apply_correction(ref, trend, kind)
-        p3 = DQM.adjust(sim3)
+        p3 = DQM.adjust(sim3, interp="linear")
         np.testing.assert_array_almost_equal(p3[middle], ref3[middle], 1)
 
     @pytest.mark.parametrize("kind,name", [(ADDITIVE, "tas"), (MULTIPLICATIVE, "pr")])
@@ -201,6 +200,17 @@ class TestDQM:
         np.testing.assert_almost_equal(p.mean(), 41.6, 0)
         np.testing.assert_almost_equal(p.std(), 15.0, 0)
 
+    def test_group_norm(self, tas_series):
+        x = np.arange(365)
+        ref = tas_series(273 - 10 * np.cos(2 * np.pi * x / 365), start="2001-01-01")
+        hist = sim = tas_series(
+            273 - 8 * np.cos(2 * np.pi * x / 365), start="2001-01-01"
+        )
+        DQM = DetrendedQuantileMapping(group="time.month")
+        DQM.train(ref, hist)
+        scen = DQM.adjust(sim, interp="linear")
+        xr.testing.assert_allclose(ref, scen, rtol=1.5e-4)
+
 
 class TestQDM:
     @pytest.mark.parametrize("kind,name", [(ADDITIVE, "tas"), (MULTIPLICATIVE, "pr")])
@@ -224,11 +234,9 @@ class TestQDM:
         hist = sim = series(x, name)
         ref = series(y, name)
 
-        QDM = QuantileDeltaMapping(
-            kind=kind, group="time", nquantiles=10, interp="linear"
-        )
+        QDM = QuantileDeltaMapping(kind=kind, group="time", nquantiles=10,)
         QDM.train(ref, hist)
-        p = QDM.adjust(sim)
+        p = QDM.adjust(sim, interp="linear")
 
         q = QDM.ds.coords["quantiles"]
         expected = get_correction(xd.ppf(q), yd.ppf(q), kind)
@@ -241,8 +249,12 @@ class TestQDM:
         middle = (u > 1e-2) * (u < 0.99)
         np.testing.assert_array_almost_equal(p[middle], ref[middle], 1)
 
+    @pytest.mark.parametrize("use_dask", [True, False])
     @pytest.mark.parametrize("kind,name", [(ADDITIVE, "tas"), (MULTIPLICATIVE, "pr")])
-    def test_mon_U(self, mon_series, series, mon_triangular, kind, name):
+    @pytest.mark.parametrize("add_dims", [True, False])
+    def test_mon_U(
+        self, mon_series, series, mon_triangular, add_dims, kind, name, use_dask
+    ):
         """
         Train on
         hist: U
@@ -263,6 +275,15 @@ class TestQDM:
 
         # Test train
         hist = sim = series(x, name)
+        if use_dask:
+            sim = sim.chunk({"time": -1})
+        if add_dims:
+            hist = hist.expand_dims(site=[0, 1, 2, 3, 4])
+            sim = sim.expand_dims(site=[0, 1, 2, 3, 4])
+            sel = {"site": 0}
+        else:
+            sel = {}
+
         ref = mon_series(y, name)
 
         QDM = QuantileDeltaMapping(kind=kind, group="time.month", nquantiles=40)
@@ -275,10 +296,12 @@ class TestQDM:
         expected = apply_correction(
             mon_triangular[:, np.newaxis], expected[np.newaxis, :], kind
         )
-        np.testing.assert_array_almost_equal(QDM.ds.af.sel(quantiles=q), expected, 1)
+        np.testing.assert_array_almost_equal(
+            QDM.ds.af.sel(quantiles=q, **sel), expected, 1
+        )
 
         # Test predict
-        np.testing.assert_array_almost_equal(p, ref, 1)
+        np.testing.assert_array_almost_equal(p.isel(**sel), ref, 1)
 
     def test_cannon(self, cannon_2015_dist, cannon_2015_rvs):
         ref, hist, sim = cannon_2015_rvs(15000, random=False)
@@ -326,11 +349,9 @@ class TestQM:
         # Test train
         hist = sim = series(x, name)
         ref = series(y, name)
-        QM = EmpiricalQuantileMapping(
-            kind=kind, group="time", nquantiles=50, interp="linear"
-        )
+        QM = EmpiricalQuantileMapping(kind=kind, group="time", nquantiles=50,)
         QM.train(ref, hist)
-        p = QM.adjust(sim)
+        p = QM.adjust(sim, interp="linear")
 
         q = QM.ds.coords["quantiles"]
         expected = get_correction(xd.ppf(q), yd.ppf(q), kind)
diff --git a/tests/test_sdba/test_utils.py b/tests/test_sdba/test_utils.py
index 6690fd02..db822a1e 100644
--- a/tests/test_sdba/test_utils.py
+++ b/tests/test_sdba/test_utils.py
@@ -109,3 +109,18 @@ def test_interp_on_quantiles(shape, group, method):
             fut_corr.values, obs.where(fut != 1000).values, rtol=2e-3
         )
         xr.testing.assert_equal(fut_corr.isnull(), fut == 1000)
+
+
+@pytest.mark.parametrize("use_dask", [True, False])
+def test_rank(use_dask):
+    arr = np.random.random_sample(size=(10, 10, 1000))
+    da = xr.DataArray(arr, dims=("x", "y", "time"))
+
+    if use_dask:
+        da = da.chunk({"x": 1})
+
+    ranks = u.rank(da, dim="time", pct=False)
+
+    exp = arr.argsort().argsort() + 1
+
+    np.testing.assert_array_equal(ranks.values, exp)