diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index f2de0f723..f5d6ff905 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -61,6 +61,7 @@ repos: hooks: - id: nbstripout files: '.ipynb' + args: [ "--extra-keys", "metadata.kernelspec" ] - repo: https://github.com/pycqa/pydocstyle rev: 6.3.0 hooks: diff --git a/CHANGES.rst b/CHANGES.rst index b530cf4b5..3f600abf9 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -4,7 +4,7 @@ Changelog v0.46.0 (unreleased) -------------------- -Contributors to this version: Éric Dupuis (:user:`coxipi`), Trevor James Smith (:user:`Zeitsperre`), David Huard (:user:`huard`), Pascal Bourgault (:user:`aulemahal`). +Contributors to this version: David Huard (:user:`huard`), Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Éric Dupuis (:user:`coxipi`). New features and enhancements ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -15,6 +15,7 @@ New features and enhancements * Refactored the indicators page with the addition of a search bar. * Indicator ``generic.stats`` now accepts any frequency (previously only daily). (:pull:`1498`). * Added argument ``out_units`` to ``select_resample_op`` to bypass limitations of ``to_agg_units`` in custom indicators. Add "var" to supported operations in ``to_agg_units``. (:pull:`1498`). +* `adapt_freq_thresh` argument added `sdba` training functions, allowing to perform frequency adaptation appropriately in each map block. (:pull:`1407`). Bug fixes ^^^^^^^^^ diff --git a/docs/notebooks/analogs.ipynb b/docs/notebooks/analogs.ipynb index 4c292eea4..1cc7aebaf 100644 --- a/docs/notebooks/analogs.ipynb +++ b/docs/notebooks/analogs.ipynb @@ -248,11 +248,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/cli.ipynb b/docs/notebooks/cli.ipynb index 9027f82e6..aac26b04e 100644 --- a/docs/notebooks/cli.ipynb +++ b/docs/notebooks/cli.ipynb @@ -407,11 +407,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/customize.ipynb b/docs/notebooks/customize.ipynb index 506ae4fb0..a5f8a9b34 100644 --- a/docs/notebooks/customize.ipynb +++ b/docs/notebooks/customize.ipynb @@ -206,11 +206,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/ensembles-advanced.ipynb b/docs/notebooks/ensembles-advanced.ipynb index e8a9c24de..572da5910 100644 --- a/docs/notebooks/ensembles-advanced.ipynb +++ b/docs/notebooks/ensembles-advanced.ipynb @@ -389,11 +389,6 @@ ], "metadata": { "celltoolbar": "Edit Metadata", - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/ensembles.ipynb b/docs/notebooks/ensembles.ipynb index 94ddfbc40..0a48b600c 100644 --- a/docs/notebooks/ensembles.ipynb +++ b/docs/notebooks/ensembles.ipynb @@ -279,11 +279,6 @@ ], "metadata": { "celltoolbar": "Edit Metadata", - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/example.ipynb b/docs/notebooks/example.ipynb index bf3a075b1..26d8f0b9b 100644 --- a/docs/notebooks/example.ipynb +++ b/docs/notebooks/example.ipynb @@ -693,11 +693,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", @@ -709,11 +704,6 @@ "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.10" - }, - "vscode": { - "interpreter": { - "hash": "978613a455e3ebf412cb0984a5d2225f42b37aab91916c56e2a684ef83795a8e" - } } }, "nbformat": 4, diff --git a/docs/notebooks/extendxclim.ipynb b/docs/notebooks/extendxclim.ipynb index e162781a2..d4e9a4c04 100644 --- a/docs/notebooks/extendxclim.ipynb +++ b/docs/notebooks/extendxclim.ipynb @@ -604,11 +604,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/frequency_analysis.ipynb b/docs/notebooks/frequency_analysis.ipynb index 44a16a1bf..0e4c0d086 100644 --- a/docs/notebooks/frequency_analysis.ipynb +++ b/docs/notebooks/frequency_analysis.ipynb @@ -176,11 +176,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/partitioning.ipynb b/docs/notebooks/partitioning.ipynb index eeb7311a3..76ccf61f6 100644 --- a/docs/notebooks/partitioning.ipynb +++ b/docs/notebooks/partitioning.ipynb @@ -228,11 +228,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/sdba-advanced.ipynb b/docs/notebooks/sdba-advanced.ipynb index 308a98b6e..0354a0eb7 100644 --- a/docs/notebooks/sdba-advanced.ipynb +++ b/docs/notebooks/sdba-advanced.ipynb @@ -692,6 +692,118 @@ "plt.legend()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Frequency adaption with a rolling window" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the previous example, we performed bias adjustment with a rolling window. Here we show how to include frequency adaptation (see `sdba.ipynb` for the simple case `group=\"time\"`). We first generate the same precipitation dataset used in `sdba.ipynb`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import xarray as xr\n", + "\n", + "t = xr.cftime_range(\"2000-01-01\", \"2030-12-31\", freq=\"D\", calendar=\"noleap\")\n", + "\n", + "vals = np.random.randint(0, 1000, size=(t.size,)) / 100\n", + "vals_ref = (4 ** np.where(vals < 9, vals / 100, vals)) / 3e6\n", + "vals_sim = (\n", + " (1 + 0.1 * np.random.random_sample((t.size,)))\n", + " * (4 ** np.where(vals < 9.5, vals / 100, vals))\n", + " / 3e6\n", + ")\n", + "\n", + "pr_ref = xr.DataArray(\n", + " vals_ref, coords={\"time\": t}, dims=(\"time\",), attrs={\"units\": \"mm/day\"}\n", + ")\n", + "pr_ref = pr_ref.sel(time=slice(\"2000\", \"2015\"))\n", + "pr_sim = xr.DataArray(\n", + " vals_sim, coords={\"time\": t}, dims=(\"time\",), attrs={\"units\": \"mm/day\"}\n", + ")\n", + "pr_hist = pr_sim.sel(time=slice(\"2000\", \"2015\"))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Bias adjustment on a rolling window can be performed in the same way as shown in `sdba.ipynb`, but instead of being a single string precising the time grouping (e.g. `time.month`), the `group` argument is built with `sdba.Grouper` function" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "# adapt_freq with a sdba.Grouper\n", + "from xclim import sdba\n", + "\n", + "group = sdba.Grouper(\"time.dayofyear\", window=31)\n", + "sim_ad, pth, dP0 = sdba.processing.adapt_freq(\n", + " pr_ref, pr_sim, thresh=\"0.05 mm d-1\", group=group\n", + ")\n", + "QM_ad = sdba.EmpiricalQuantileMapping.train(\n", + " pr_ref, sim_ad, nquantiles=15, kind=\"*\", group=group\n", + ")\n", + "scen_ad = QM_ad.adjust(pr_sim)\n", + "\n", + "pr_ref.sel(time=\"2010\").plot(alpha=0.9, label=\"Reference\")\n", + "pr_sim.sel(time=\"2010\").plot(alpha=0.7, label=\"Model - biased\")\n", + "scen_ad.sel(time=\"2010\").plot(alpha=0.6, label=\"Model - adjusted\")\n", + "plt.legend()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the figure above, `scen` occasionally has small peaks where `sim` is 0, indicating that there are more \"dry days\" (days with almost no precipitation) in `hist` than in `ref`. The frequency-adaptation [Themeßl et al. (2010)](https://doi.org/10.1007/s10584-011-0224-4) performed in the step above only worked partially. \n", + "\n", + "The reason for this is the following. The first step above combines precipitations in 365 overlapping blocks of 31 days * Y years, one block for each day of the year. Each block is adapted, and the 16th day-of-year slice (at the center of the block) is assigned to the corresponding day-of-year in the adapted dataset `sim_ad`. As we proceed to the training, we re-form those 31 days * Y years blocks, but this step does not invert the last one: There can still be more zeroes in the simulation than in the reference. \n", + "\n", + "To alleviate this issue, another way of proceeding is to perform a frequency adaptation on the blocks, and then use the same blocks in the training step, as we show below.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# adapt_freq directly in the training step\n", + "group = sdba.Grouper(\"time.dayofyear\", window=31)\n", + "\n", + "QM_ad = sdba.EmpiricalQuantileMapping.train(\n", + " pr_ref,\n", + " sim_ad,\n", + " nquantiles=15,\n", + " kind=\"*\",\n", + " group=group,\n", + " adapt_freq_thresh=\"0.05 mm d-1\",\n", + ")\n", + "scen_ad = QM_ad.adjust(pr_sim)\n", + "\n", + "pr_ref.sel(time=\"2010\").plot(alpha=0.9, label=\"Reference\")\n", + "pr_sim.sel(time=\"2010\").plot(alpha=0.7, label=\"Model - biased\")\n", + "scen_ad.sel(time=\"2010\").plot(alpha=0.6, label=\"Model - adjusted\")\n", + "plt.legend()" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -834,11 +946,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", @@ -849,7 +956,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.9" + "version": "3.11.4" }, "toc": { "base_numbering": 1, diff --git a/docs/notebooks/sdba.ipynb b/docs/notebooks/sdba.ipynb index 86882a4b6..7f0e75f90 100644 --- a/docs/notebooks/sdba.ipynb +++ b/docs/notebooks/sdba.ipynb @@ -676,11 +676,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", @@ -691,7 +686,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.9" + "version": "3.11.4" }, "toc": { "base_numbering": 1, diff --git a/docs/notebooks/units.ipynb b/docs/notebooks/units.ipynb index 1273da58a..5ff1dbe6b 100644 --- a/docs/notebooks/units.ipynb +++ b/docs/notebooks/units.ipynb @@ -230,11 +230,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/usage.ipynb b/docs/notebooks/usage.ipynb index 0d703c29b..69a01c629 100644 --- a/docs/notebooks/usage.ipynb +++ b/docs/notebooks/usage.ipynb @@ -357,11 +357,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/xclim_training/Exercices.ipynb b/docs/notebooks/xclim_training/Exercices.ipynb index 3eb52a9bf..abdd32d61 100644 --- a/docs/notebooks/xclim_training/Exercices.ipynb +++ b/docs/notebooks/xclim_training/Exercices.ipynb @@ -85,11 +85,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/xclim_training/XARRAY_calcul_moy_saisonniere.ipynb b/docs/notebooks/xclim_training/XARRAY_calcul_moy_saisonniere.ipynb index f52798420..70c21e78b 100644 --- a/docs/notebooks/xclim_training/XARRAY_calcul_moy_saisonniere.ipynb +++ b/docs/notebooks/xclim_training/XARRAY_calcul_moy_saisonniere.ipynb @@ -227,11 +227,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/xclim_training/XCLIM Demo - Ensembles.ipynb b/docs/notebooks/xclim_training/XCLIM Demo - Ensembles.ipynb index 4f4be7967..e8d1bfd98 100644 --- a/docs/notebooks/xclim_training/XCLIM Demo - Ensembles.ipynb +++ b/docs/notebooks/xclim_training/XCLIM Demo - Ensembles.ipynb @@ -468,11 +468,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/xclim_training/XCLIM_calculate_index-Exemple.ipynb b/docs/notebooks/xclim_training/XCLIM_calculate_index-Exemple.ipynb index 7b3caf952..a5e055d0b 100644 --- a/docs/notebooks/xclim_training/XCLIM_calculate_index-Exemple.ipynb +++ b/docs/notebooks/xclim_training/XCLIM_calculate_index-Exemple.ipynb @@ -498,11 +498,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/xclim_training/finch.ipynb b/docs/notebooks/xclim_training/finch.ipynb index 91f1419c9..3b59f3ff8 100644 --- a/docs/notebooks/xclim_training/finch.ipynb +++ b/docs/notebooks/xclim_training/finch.ipynb @@ -145,11 +145,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/xclim_training/intro_xarray.ipynb b/docs/notebooks/xclim_training/intro_xarray.ipynb index e07713293..d566d1d72 100644 --- a/docs/notebooks/xclim_training/intro_xarray.ipynb +++ b/docs/notebooks/xclim_training/intro_xarray.ipynb @@ -423,11 +423,6 @@ ], "metadata": { "celltoolbar": "Slideshow", - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/docs/notebooks/xclim_training/readme.ipynb b/docs/notebooks/xclim_training/readme.ipynb index 9c38191f0..68b657244 100644 --- a/docs/notebooks/xclim_training/readme.ipynb +++ b/docs/notebooks/xclim_training/readme.ipynb @@ -99,11 +99,6 @@ } ], "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", diff --git a/xclim/sdba/_adjustment.py b/xclim/sdba/_adjustment.py index 407b61e1f..a28f9c03a 100644 --- a/xclim/sdba/_adjustment.py +++ b/xclim/sdba/_adjustment.py @@ -9,21 +9,33 @@ import numpy as np import xarray as xr +from xclim.core.units import convert_units_to, infer_context, units from xclim.indices.stats import _fitfunc_1d # noqa from . import nbutils as nbu from . import utils as u +from ._processing import _adapt_freq from .base import Grouper, map_blocks, map_groups from .detrending import PolyDetrend from .processing import escore +def _adapt_freq_hist(ds, adapt_freq_thresh): + """Adapt frequency of null values of `hist` in order to match `ref`.""" + with units.context(infer_context(ds.ref.attrs.get("standard_name"))): + thresh = convert_units_to(adapt_freq_thresh, ds.ref) + dim = ["time"] + ["window"] * ("window" in ds.hist.dims) + return _adapt_freq.func( + xr.Dataset(dict(sim=ds.hist, ref=ds.ref)), thresh=thresh, dim=dim + ).sim_ad + + @map_groups( af=[Grouper.PROP, "quantiles"], hist_q=[Grouper.PROP, "quantiles"], scaling=[Grouper.PROP], ) -def dqm_train(ds, *, dim, kind, quantiles) -> xr.Dataset: +def dqm_train(ds, *, dim, kind, quantiles, adapt_freq_thresh) -> xr.Dataset: """Train step on one group. Notes @@ -31,16 +43,22 @@ def dqm_train(ds, *, dim, kind, quantiles) -> xr.Dataset: Dataset must contain the following variables: ref : training target hist : training data + + adapt_freq_thresh : str | None + Threshold for frequency adaptation. See :py:class:`xclim.sdba.processing.adapt_freq` for details. + Default is None, meaning that frequency adaptation is not performed. """ + hist = _adapt_freq_hist(ds, adapt_freq_thresh) if adapt_freq_thresh else ds.hist + refn = u.apply_correction(ds.ref, u.invert(ds.ref.mean(dim), kind), kind) - histn = u.apply_correction(ds.hist, u.invert(ds.hist.mean(dim), kind), kind) + histn = u.apply_correction(hist, u.invert(hist.mean(dim), kind), kind) ref_q = nbu.quantile(refn, quantiles, dim) hist_q = nbu.quantile(histn, quantiles, dim) af = u.get_correction(hist_q, ref_q, kind) mu_ref = ds.ref.mean(dim) - mu_hist = ds.hist.mean(dim) + mu_hist = hist.mean(dim) scaling = u.get_correction(mu_hist, mu_ref, kind=kind) return xr.Dataset(data_vars=dict(af=af, hist_q=hist_q, scaling=scaling)) @@ -50,15 +68,20 @@ def dqm_train(ds, *, dim, kind, quantiles) -> xr.Dataset: af=[Grouper.PROP, "quantiles"], hist_q=[Grouper.PROP, "quantiles"], ) -def eqm_train(ds, *, dim, kind, quantiles) -> xr.Dataset: +def eqm_train(ds, *, dim, kind, quantiles, adapt_freq_thresh) -> xr.Dataset: """EQM: Train step on one group. Dataset variables: ref : training target hist : training data + + adapt_freq_thresh : str | None + Threshold for frequency adaptation. See :py:class:`xclim.sdba.processing.adapt_freq` for details. + Default is None, meaning that frequency adaptation is not performed. """ + hist = _adapt_freq_hist(ds, adapt_freq_thresh) if adapt_freq_thresh else ds.hist ref_q = nbu.quantile(ds.ref, quantiles, dim) - hist_q = nbu.quantile(ds.hist, quantiles, dim) + hist_q = nbu.quantile(hist, quantiles, dim) af = u.get_correction(hist_q, ref_q, kind) diff --git a/xclim/sdba/_processing.py b/xclim/sdba/_processing.py index b297b31e9..164129abe 100644 --- a/xclim/sdba/_processing.py +++ b/xclim/sdba/_processing.py @@ -11,10 +11,11 @@ import numpy as np import xarray as xr +from xarray.core.utils import get_temp_dimname from . import nbutils as nbu from .base import Grouper, map_groups -from .utils import ADDITIVE, apply_correction, ecdf, invert +from .utils import ADDITIVE, apply_correction, ecdf, invert, rank @map_groups( @@ -74,22 +75,15 @@ def _adapt_freq( pth = nbu.vecquantiles(ds.ref, P0_sim, dim).where(dP0 > 0) # Probabilities and quantiles computed within all dims, but correction along the first one only. - if "window" in dim: - # P0_sim was computed using the window, but only the original time series is corrected. - # Grouper.apply does this step, but if done here it makes the code faster. - sim = ds.sim.isel(window=(ds.sim.window.size - 1) // 2) - else: - sim = ds.sim - dim = dim[0] - + sim = ds.sim # Get the percentile rank of each value in sim. - rank = sim.rank(dim, pct=True) + rnk = rank(sim, dim=dim, pct=True) # Frequency-adapted sim sim_ad = sim.where( dP0 < 0, # dP0 < 0 means no-adaptation. sim.where( - (rank < P0_ref) | (rank > P0_sim), # Preserve current values + (rnk < P0_ref) | (rnk > P0_sim), # Preserve current values # Generate random numbers ~ U[T0, Pth] (pth.broadcast_like(sim) - thresh) * np.random.random_sample(size=sim.shape) diff --git a/xclim/sdba/adjustment.py b/xclim/sdba/adjustment.py index 9b7db46ae..edf166b1d 100644 --- a/xclim/sdba/adjustment.py +++ b/xclim/sdba/adjustment.py @@ -328,6 +328,9 @@ class EmpiricalQuantileMapping(TrainAdjust): group : Union[str, Grouper] The grouping information. See :py:class:`xclim.sdba.base.Grouper` for details. Default is "time", meaning an single adjustment group along dimension "time". + adapt_freq_thresh : str | None + Threshold for frequency adaptation. See :py:class:`xclim.sdba.processing.adapt_freq` for details. + Default is None, meaning that frequency adaptation is not performed. Adjust step: @@ -352,6 +355,7 @@ def _train( nquantiles: int | np.ndarray = 20, kind: str = ADDITIVE, group: str | Grouper = "time", + adapt_freq_thresh: str | None = None, ): if np.isscalar(nquantiles): quantiles = equally_spaced_nodes(nquantiles).astype(ref.dtype) @@ -363,6 +367,7 @@ def _train( group=group, kind=kind, quantiles=quantiles, + adapt_freq_thresh=adapt_freq_thresh, ) ds.af.attrs.update( @@ -417,6 +422,9 @@ class DetrendedQuantileMapping(TrainAdjust): group : Union[str, Grouper] The grouping information. See :py:class:`xclim.sdba.base.Grouper` for details. Default is "time", meaning a single adjustment group along dimension "time". + adapt_freq_thresh : str | None + Threshold for frequency adaptation. See :py:class:`xclim.sdba.processing.adapt_freq` for details. + Default is None, meaning that frequency adaptation is not performed. Adjust step: @@ -445,6 +453,7 @@ def _train( nquantiles: int | np.ndarray = 20, kind: str = ADDITIVE, group: str | Grouper = "time", + adapt_freq_thresh: str | None = None, ): if group.prop not in ["group", "dayofyear"]: warn( @@ -461,6 +470,7 @@ def _train( group=group, quantiles=quantiles, kind=kind, + adapt_freq_thresh=adapt_freq_thresh, ) ds.af.attrs.update( diff --git a/xclim/sdba/utils.py b/xclim/sdba/utils.py index 9eba7ca4f..fbe42d7a0 100644 --- a/xclim/sdba/utils.py +++ b/xclim/sdba/utils.py @@ -14,6 +14,7 @@ from dask import array as dsk from scipy.interpolate import griddata, interp1d from scipy.stats import spearmanr +from xarray.core.utils import get_temp_dimname from xclim.core.calendar import _interpolate_doy_calendar # noqa from xclim.core.utils import ensure_chunk_size @@ -461,7 +462,9 @@ def interp_on_quantiles( ) -def rank(da: xr.DataArray, dim: str = "time", pct: bool = False) -> xr.DataArray: +def rank( + da: xr.DataArray, dim: str | list[str] = "time", pct: bool = False +) -> xr.DataArray: """Ranks data along a dimension. Replicates `xr.DataArray.rank` but as a function usable in a Grouper.apply(). Xarray's docstring is below: @@ -469,12 +472,14 @@ def rank(da: xr.DataArray, dim: str = "time", pct: bool = False) -> xr.DataArray Equal values are assigned a rank that is the average of the ranks that would have been otherwise assigned to all the values within that set. Ranks begin at 1, not 0. If pct, computes percentage ranks, ranging from 0 to 1. + A list of dimensions can be provided and the ranks are then computed separately for each dimension. + Parameters ---------- da: xr.DataArray Source array. - dim : str, hashable - Dimension over which to compute rank. + dim : str | list[str], hashable + Dimension(s) over which to compute rank. pct : bool, optional If True, compute percentage ranks, otherwise compute integer ranks. Percentage ranks range from 0 to 1, in opposition to xarray's implementation, @@ -493,11 +498,26 @@ def rank(da: xr.DataArray, dim: str = "time", pct: bool = False) -> xr.DataArray -------- xarray.DataArray.rank """ - rnk = da.rank(dim, pct=pct) + da_dims, da_coords = da.dims, da.coords + dims = dim if isinstance(dim, list) else [dim] + rnk_dim = dims[0] if len(dims) == 1 else get_temp_dimname(da_dims, "temp") + + # multi-dimensional ranking through stacking + if len(dims) > 1: + da = da.stack(**{rnk_dim: dims}) + rnk = da.rank(rnk_dim, pct=pct) + if pct: - mn = rnk.min(dim) - mx = rnk.max(dim) - return mx * (rnk - mn) / (mx - mn) + mn = rnk.min(rnk_dim) + mx = rnk.max(rnk_dim) + rnk = mx * (rnk - mn) / (mx - mn) + + if len(dims) > 1: + rnk = ( + rnk.unstack(rnk_dim) + .transpose(*da_dims) + .drop_vars([d for d in dims if d not in da_coords]) + ) return rnk