docs/notebooks/sdba.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Statistical Downscaling and Bias-Adjustment\n",
    "\n",
    "`xclim` provides tools and utilities to ease the bias-adjustement process through its `xclim.sdba` module. Almost all adjustment algorithms conform to the `train` - `adjust` scheme, formalized within `TrainAdjust` classes. 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.\n",
    "\n",
    "This presents examples, while a bit more info and the API are given on [this page](../sdba.rst).\n",
    "\n",
    "A very simple \"Quantile Mapping\" approach is available through the \"Empirical Quantile Mapping\" object. The object is created through the `.train` method of the class, and the simulation is adjusted with `.adjust`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cftime\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "\n",
    "%matplotlib inline\n",
    "plt.style.use(\"seaborn\")\n",
    "plt.rcParams[\"figure.figsize\"] = (11, 5)\n",
    "\n",
    "# Create toy data to explore bias adjustment, here fake temperature timeseries\n",
    "t = xr.cftime_range(\"2000-01-01\", \"2030-12-31\", freq=\"D\", calendar=\"noleap\")\n",
    "ref = xr.DataArray(\n",
    "    (\n",
    "        -20 * np.cos(2 * np.pi * t.dayofyear / 365)\n",
    "        + 2 * np.random.random_sample((t.size,))\n",
    "        + 273.15\n",
    "        + 0.1 * (t - t[0]).days / 365\n",
    "    ),  # \"warming\" of 1K per decade,\n",
    "    dims=(\"time\",),\n",
    "    coords={\"time\": t},\n",
    "    attrs={\"units\": \"K\"},\n",
    ")\n",
    "sim = xr.DataArray(\n",
    "    (\n",
    "        -18 * np.cos(2 * np.pi * t.dayofyear / 365)\n",
    "        + 2 * np.random.random_sample((t.size,))\n",
    "        + 273.15\n",
    "        + 0.11 * (t - t[0]).days / 365\n",
    "    ),  # \"warming\" of 1.1K per decade\n",
    "    dims=(\"time\",),\n",
    "    coords={\"time\": t},\n",
    "    attrs={\"units\": \"K\"},\n",
    ")\n",
    "\n",
    "ref = ref.sel(time=slice(None, \"2015-01-01\"))\n",
    "hist = sim.sel(time=slice(None, \"2015-01-01\"))\n",
    "\n",
    "ref.plot(label=\"Reference\")\n",
    "sim.plot(label=\"Model\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xclim import sdba\n",
    "\n",
    "QM = sdba.EmpiricalQuantileMapping.train(\n",
    "    ref, hist, nquantiles=15, group=\"time\", kind=\"+\"\n",
    ")\n",
    "scen = QM.adjust(sim, extrapolation=\"constant\", interp=\"nearest\")\n",
    "\n",
    "ref.groupby(\"time.dayofyear\").mean().plot(label=\"Reference\")\n",
    "hist.groupby(\"time.dayofyear\").mean().plot(label=\"Model - biased\")\n",
    "scen.sel(time=slice(\"2000\", \"2015\")).groupby(\"time.dayofyear\").mean().plot(\n",
    "    label=\"Model - adjusted - 2000-15\", linestyle=\"--\"\n",
    ")\n",
    "scen.sel(time=slice(\"2015\", \"2030\")).groupby(\"time.dayofyear\").mean().plot(\n",
    "    label=\"Model - adjusted - 2015-30\", linestyle=\"--\"\n",
    ")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous example, a simple Quantile Mapping algorithm was used with 15 quantiles and one group of values. The model performs well, but our toy data is also quite smooth and well-behaved so this is not surprising. A more complex example could have biais distribution varying strongly across months. To perform the adjustment with different factors for each months, one can pass `group='time.month'`. Moreover, to reduce the risk of sharp change in the adjustment at the interface of the months, `interp='linear'` can be passed to `adjust` and the adjustment factors will be interpolated linearly. Ex: the factors for the 1st of May will be the average of those for april and those for may."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "QM_mo = sdba.EmpiricalQuantileMapping.train(\n",
    "    ref, hist, nquantiles=15, group=\"time.month\", kind=\"+\"\n",
    ")\n",
    "scen = QM_mo.adjust(sim, extrapolation=\"constant\", interp=\"linear\")\n",
    "\n",
    "ref.groupby(\"time.dayofyear\").mean().plot(label=\"Reference\")\n",
    "hist.groupby(\"time.dayofyear\").mean().plot(label=\"Model - biased\")\n",
    "scen.sel(time=slice(\"2000\", \"2015\")).groupby(\"time.dayofyear\").mean().plot(\n",
    "    label=\"Model - adjusted - 2000-15\", linestyle=\"--\"\n",
    ")\n",
    "scen.sel(time=slice(\"2015\", \"2030\")).groupby(\"time.dayofyear\").mean().plot(\n",
    "    label=\"Model - adjusted - 2015-30\", linestyle=\"--\"\n",
    ")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The training data (here the adjustment factors) is available for inspection in the `ds` attribute of the adjustment object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "QM_mo.ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "QM_mo.ds.af.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Grouping\n",
    "\n",
    "For basic time period grouping (months, day of year, season), passing a string to the methods needing it is sufficient. Most methods acting on grouped data also accept a `window` int argument to pad the groups with data from adjacent ones. Units of `window` are the sampling frequency of the main grouping dimension (usually `time`). For more complex grouping, or simply for clarity, one can pass a `xclim.sdba.base.Grouper` directly.\n",
    "\n",
    "Example here with another, simpler, adjustment method. Here we want `sim` to be scaled so that its mean fits the one of `ref`. Scaling factors are to be computed separately for each day of the year, but including 15 days on either side of the day. This means that the factor for the 1st of May is computed including all values from the 16th of April to the 15th of May (of all years)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "group = sdba.Grouper(\"time.dayofyear\", window=31)\n",
    "QM_doy = sdba.Scaling.train(ref, hist, group=group, kind=\"+\")\n",
    "scen = QM_doy.adjust(sim)\n",
    "\n",
    "ref.groupby(\"time.dayofyear\").mean().plot(label=\"Reference\")\n",
    "hist.groupby(\"time.dayofyear\").mean().plot(label=\"Model - biased\")\n",
    "scen.sel(time=slice(\"2000\", \"2015\")).groupby(\"time.dayofyear\").mean().plot(\n",
    "    label=\"Model - adjusted - 2000-15\", linestyle=\"--\"\n",
    ")\n",
    "scen.sel(time=slice(\"2015\", \"2030\")).groupby(\"time.dayofyear\").mean().plot(\n",
    "    label=\"Model - adjusted - 2015-30\", linestyle=\"--\"\n",
    ")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "QM_doy.ds.af.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modular approach\n",
    "\n",
    "The `sdba` module adopts a modular approach instead of implementing published and named methods directly.\n",
    "A generic bias adjustment process is laid out as follows:\n",
    "\n",
    "- preprocessing on `ref`, `hist` and `sim` (using methods in `xclim.sdba.processing` or `xclim.sdba.detrending`)\n",
    "- creating and training the adjustment object `Adj = Adjustment.train(obs, hist, **kwargs)` (from `xclim.sdba.adjustment`)\n",
    "- adjustment `scen = Adj.adjust(sim, **kwargs)`\n",
    "- post-processing on `scen` (for example: re-trending)\n",
    "\n",
    "The train-adjust approach allows to inspect the trained adjustment object. The training information is stored in the underlying `Adj.ds` dataset and often has a `af` variable with the adjustment factors. Its layout and the other available variables vary between the different algorithm, refer to their part of the API docs.\n",
    "\n",
    "For heavy processing, this separation allows the computation and writing to disk of the training dataset before performing the adjustment(s). See the [advanced notebook](sdba-advanced.ipynb).\n",
    "\n",
    "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.\n",
    "\n",
    "### First example : pr and frequency adaptation\n",
    "\n",
    "The next example generates fake precipitation data and adjusts the `sim` timeseries but also adds a step where the dry-day frequency of `hist` is adapted so that is fits the one of `ref`. This ensures well-behaved adjustment factors for the smaller quantiles. Note also that we are passing `kind='*'` to use the multiplicative mode. Adjustment factors will be multiplied/divided instead of being added/substracted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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\"))\n",
    "\n",
    "pr_ref.plot(alpha=0.9, label=\"Reference\")\n",
    "pr_sim.plot(alpha=0.7, label=\"Model\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1st try without adapt_freq\n",
    "QM = sdba.EmpiricalQuantileMapping.train(\n",
    "    pr_ref, pr_hist, nquantiles=15, kind=\"*\", group=\"time\"\n",
    ")\n",
    "scen = QM.adjust(pr_sim)\n",
    "\n",
    "pr_ref.sel(time=\"2010\").plot(alpha=0.9, label=\"Reference\")\n",
    "pr_hist.sel(time=\"2010\").plot(alpha=0.7, label=\"Model - biased\")\n",
    "scen.sel(time=\"2010\").plot(alpha=0.6, label=\"Model - adjusted\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the figure above, `scen` has small peaks where `sim` is 0. This problem originates from the fact that there are more \"dry days\" (days with almost no precipitation) in `hist` than in `ref`. The next example works around the problem using frequency-adaptation, as described in [Themeßl et al. (2010)](https://doi.org/10.1007/s10584-011-0224-4)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2nd try with adapt_freq\n",
    "sim_ad, pth, dP0 = sdba.processing.adapt_freq(\n",
    "    pr_ref, pr_sim, thresh=\"0.05 mm d-1\", group=\"time\"\n",
    ")\n",
    "QM_ad = sdba.EmpiricalQuantileMapping.train(\n",
    "    pr_ref, sim_ad, nquantiles=15, kind=\"*\", group=\"time\"\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": [
    "### Second example: tas and detrending\n",
    "\n",
    "The next example reuses the fake temperature timeseries generated at the beginning and applies the same QM adjustment method. However, for a better adjustment, we will scale sim to ref and then detrend the series, assuming the trend is linear. When `sim` (or `sim_scl`) is detrended, its values are now anomalies, so we need to normalize `ref` and `hist` so we can compare similar values.\n",
    "\n",
    "This process is detailed here to show how the sdba module should be used in custom adjustment processes, but this specific method also exists as `sdba.DetrendedQuantileMapping` and is based on [Cannon et al. 2015](https://doi.org/10.1175/JCLI-D-14-00754.1). However, `DetrendedQuantileMapping` normalizes over a `time.dayofyear` group, regardless of what is passed in the  `group` argument. As done here, it is anyway recommended to use `dayofyear` groups when normalizing, especially for variables with strong seasonal variations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "doy_win31 = sdba.Grouper(\"time.dayofyear\", window=15)\n",
    "Sca = sdba.Scaling.train(ref, hist, group=doy_win31, kind=\"+\")\n",
    "sim_scl = Sca.adjust(sim)\n",
    "\n",
    "detrender = sdba.detrending.PolyDetrend(degree=1, group=\"time.dayofyear\", kind=\"+\")\n",
    "sim_fit = detrender.fit(sim_scl)\n",
    "sim_detrended = sim_fit.detrend(sim_scl)\n",
    "\n",
    "ref_n, _ = sdba.processing.normalize(ref, group=doy_win31, kind=\"+\")\n",
    "hist_n, _ = sdba.processing.normalize(hist, group=doy_win31, kind=\"+\")\n",
    "\n",
    "QM = sdba.EmpiricalQuantileMapping.train(\n",
    "    ref_n, hist_n, nquantiles=15, group=\"time.month\", kind=\"+\"\n",
    ")\n",
    "scen_detrended = QM.adjust(sim_detrended, extrapolation=\"constant\", interp=\"nearest\")\n",
    "scen = sim_fit.retrend(scen_detrended)\n",
    "\n",
    "\n",
    "ref.groupby(\"time.dayofyear\").mean().plot(label=\"Reference\")\n",
    "sim.groupby(\"time.dayofyear\").mean().plot(label=\"Model - biased\")\n",
    "scen.sel(time=slice(\"2000\", \"2015\")).groupby(\"time.dayofyear\").mean().plot(\n",
    "    label=\"Model - adjusted - 2000-15\", linestyle=\"--\"\n",
    ")\n",
    "scen.sel(time=slice(\"2015\", \"2030\")).groupby(\"time.dayofyear\").mean().plot(\n",
    "    label=\"Model - adjusted - 2015-30\", linestyle=\"--\"\n",
    ")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Third example : Multi-method protocol - Hnilica et al. 2017\n",
    "In [their paper of 2017](https://doi.org/10.1002/joc.4890), Hnilica, Hanel and Puš present a bias-adjustment method based on the principles of Principal Components Analysis. The idea is simple : use principal components to define coordinates on the reference and on the simulation and then transform the simulation data from the latter to the former. Spatial correlation can thus be conserved by taking different points as the dimensions of the transform space. The method was demonstrated in the article by bias-adjusting precipitation over different drainage basins.\n",
    "\n",
    "The same method could be used for multivariate adjustment. The principle would be the same, concatening the different variables into a single dataset along a new dimension. An example is given in the [advanced notebook](sdba-advanced.ipynb).\n",
    "\n",
    "Here we show how the modularity of `xclim.sdba` can be used to construct a quite complex adjustment protocol involving two adjustment methods : quantile mapping and principal components. Evidently, as this example uses only 2 years of data, it is not complete. It is meant to show how the adjustment functions and how the API can be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We are using xarray's \"air_temperature\" dataset\n",
    "ds = xr.tutorial.open_dataset(\"air_temperature\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To get an exagerated example we select different points\n",
    "# here \"lon\" will be our dimension of two \"spatially correlated\" points\n",
    "reft = ds.air.isel(lat=21, lon=[40, 52]).drop_vars([\"lon\", \"lat\"])\n",
    "simt = ds.air.isel(lat=18, lon=[17, 35]).drop_vars([\"lon\", \"lat\"])\n",
    "\n",
    "# Principal Components Adj, no grouping and use \"lon\" as the space dimensions\n",
    "PCA = sdba.PrincipalComponents.train(reft, simt, group=\"time\", crd_dim=\"lon\")\n",
    "scen1 = PCA.adjust(simt)\n",
    "\n",
    "# QM, no grouping, 20 quantiles and additive adjustment\n",
    "EQM = sdba.EmpiricalQuantileMapping.train(\n",
    "    reft, scen1, group=\"time\", nquantiles=50, kind=\"+\"\n",
    ")\n",
    "scen2 = EQM.adjust(scen1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# some Analysis figures\n",
    "fig = plt.figure(figsize=(12, 16))\n",
    "gs = plt.matplotlib.gridspec.GridSpec(3, 2, fig)\n",
    "\n",
    "axPCA = plt.subplot(gs[0, :])\n",
    "axPCA.scatter(reft.isel(lon=0), reft.isel(lon=1), s=20, label=\"Reference\")\n",
    "axPCA.scatter(simt.isel(lon=0), simt.isel(lon=1), s=10, label=\"Simulation\")\n",
    "axPCA.scatter(scen2.isel(lon=0), scen2.isel(lon=1), s=3, label=\"Adjusted - PCA+EQM\")\n",
    "axPCA.set_xlabel(\"Point 1\")\n",
    "axPCA.set_ylabel(\"Point 2\")\n",
    "axPCA.set_title(\"PC-space\")\n",
    "axPCA.legend()\n",
    "\n",
    "refQ = reft.quantile(EQM.ds.quantiles, dim=\"time\")\n",
    "simQ = simt.quantile(EQM.ds.quantiles, dim=\"time\")\n",
    "scen1Q = scen1.quantile(EQM.ds.quantiles, dim=\"time\")\n",
    "scen2Q = scen2.quantile(EQM.ds.quantiles, dim=\"time\")\n",
    "for i in range(2):\n",
    "    if i == 0:\n",
    "        axQM = plt.subplot(gs[1, 0])\n",
    "    else:\n",
    "        axQM = plt.subplot(gs[1, 1], sharey=axQM)\n",
    "    axQM.plot(refQ.isel(lon=i), simQ.isel(lon=i), label=\"No adj\")\n",
    "    axQM.plot(refQ.isel(lon=i), scen1Q.isel(lon=i), label=\"PCA\")\n",
    "    axQM.plot(refQ.isel(lon=i), scen2Q.isel(lon=i), label=\"PCA+EQM\")\n",
    "    axQM.plot(\n",
    "        refQ.isel(lon=i), refQ.isel(lon=i), color=\"k\", linestyle=\":\", label=\"Ideal\"\n",
    "    )\n",
    "    axQM.set_title(f\"QQ plot - Point {i + 1}\")\n",
    "    axQM.set_xlabel(\"Reference\")\n",
    "    axQM.set_xlabel(\"Model\")\n",
    "    axQM.legend()\n",
    "\n",
    "axT = plt.subplot(gs[2, :])\n",
    "reft.isel(lon=0).plot(ax=axT, label=\"Reference\")\n",
    "simt.isel(lon=0).plot(ax=axT, label=\"Unadjusted sim\")\n",
    "# scen1.isel(lon=0).plot(ax=axT, label='PCA only')\n",
    "scen2.isel(lon=0).plot(ax=axT, label=\"PCA+EQM\")\n",
    "axT.legend()\n",
    "axT.set_title(\"Timeseries - Point 1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fourth example : Multivariate bias-adjustment with multiple steps - Cannon 2018\n",
    "\n",
    "This section replicates the \"MBCn\" algorithm described by [Cannon (2018)](https://doi.org/10.1007/s00382-017-3580-6). The method relies on some univariate algorithm, an adaption of the N-pdf transform of [Pitié et al. (2005)](https://ieeexplore.ieee.org/document/1544887/) and a final reordering step.\n",
    "\n",
    "In the following, we use the AHCCD and CanESM2 data are reference and simulation and we correct both `pr` and `tasmax` together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xclim.core.units import convert_units_to\n",
    "from xclim.testing import open_dataset\n",
    "\n",
    "dref = open_dataset(\n",
    "    \"sdba/ahccd_1950-2013.nc\", chunks={\"location\": 1}, drop_variables=[\"lat\", \"lon\"]\n",
    ").sel(time=slice(\"1981\", \"2010\"))\n",
    "dref = dref.assign(\n",
    "    tasmax=convert_units_to(dref.tasmax, \"K\"),\n",
    "    pr=convert_units_to(dref.pr, \"kg m-2 s-1\"),\n",
    ")\n",
    "dsim = open_dataset(\n",
    "    \"sdba/CanESM2_1950-2100.nc\", chunks={\"location\": 1}, drop_variables=[\"lat\", \"lon\"]\n",
    ")\n",
    "\n",
    "dhist = dsim.sel(time=slice(\"1981\", \"2010\"))\n",
    "dsim = dsim.sel(time=slice(\"2041\", \"2070\"))\n",
    "dref"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Perform an initial univariate adjustment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# additive for tasmax\n",
    "QDMtx = sdba.QuantileDeltaMapping.train(\n",
    "    dref.tasmax, dhist.tasmax, nquantiles=20, kind=\"+\", group=\"time\"\n",
    ")\n",
    "# Adjust both hist and sim, we'll feed both to the Npdf transform.\n",
    "scenh_tx = QDMtx.adjust(dhist.tasmax)\n",
    "scens_tx = QDMtx.adjust(dsim.tasmax)\n",
    "\n",
    "# remove == 0 values in pr:\n",
    "dref[\"pr\"] = sdba.processing.jitter_under_thresh(dref.pr, \"0.01 mm d-1\")\n",
    "dhist[\"pr\"] = sdba.processing.jitter_under_thresh(dhist.pr, \"0.01 mm d-1\")\n",
    "dsim[\"pr\"] = sdba.processing.jitter_under_thresh(dsim.pr, \"0.01 mm d-1\")\n",
    "\n",
    "# multiplicative for pr\n",
    "QDMpr = sdba.QuantileDeltaMapping.train(\n",
    "    dref.pr, dhist.pr, nquantiles=20, kind=\"*\", group=\"time\"\n",
    ")\n",
    "# Adjust both hist and sim, we'll feed both to the Npdf transform.\n",
    "scenh_pr = QDMpr.adjust(dhist.pr)\n",
    "scens_pr = QDMpr.adjust(dsim.pr)\n",
    "\n",
    "scenh = xr.Dataset(dict(tasmax=scenh_tx, pr=scenh_pr))\n",
    "scens = xr.Dataset(dict(tasmax=scens_tx, pr=scens_pr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Stack the variables to multivariate arrays and standardize them\n",
    "The standardization process ensure the mean and standard deviation of each column (variable) is 0 and 1 respectively.\n",
    "\n",
    "`hist` and `sim` are standardized together so the two series are coherent. We keep the mean and standard deviation to be reused when we build the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Stack the variables (tasmax and pr)\n",
    "ref = sdba.processing.stack_variables(dref)\n",
    "scenh = sdba.processing.stack_variables(scenh)\n",
    "scens = sdba.processing.stack_variables(scens)\n",
    "\n",
    "# Standardize\n",
    "ref, _, _ = sdba.processing.standardize(ref)\n",
    "\n",
    "allsim, savg, sstd = sdba.processing.standardize(xr.concat((scenh, scens), \"time\"))\n",
    "hist = allsim.sel(time=scenh.time)\n",
    "sim = allsim.sel(time=scens.time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Perform the N-dimensional probability density function transform\n",
    "\n",
    "The NpdfTransform will iteratively randomly rotate our arrays in the \"variables\" space and apply the univariate adjustment before rotating it back. In Cannon (2018) and Pitié et al. (2005), it can be seen that the source array's joint distribution converges toward the target's joint distribution when a large number of iterations is done."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xclim import set_options\n",
    "\n",
    "# See the advanced notebook for details on how this option work\n",
    "with set_options(sdba_extra_output=True):\n",
    "    out = sdba.adjustment.NpdfTransform.adjust(\n",
    "        ref,\n",
    "        hist,\n",
    "        sim,\n",
    "        base=sdba.QuantileDeltaMapping,  # Use QDM as the univariate adjustment.\n",
    "        base_kws={\"nquantiles\": 20, \"group\": \"time\"},\n",
    "        n_iter=20,  # perform 20 iteration\n",
    "        n_escore=1000,  # only send 1000 points to the escore metric (it is realy slow)\n",
    "    )\n",
    "\n",
    "scenh = out.scenh.rename(time_hist=\"time\")  # Bias-adjusted historical period\n",
    "scens = out.scen  # Bias-adjusted future period\n",
    "extra = out.drop_vars([\"scenh\", \"scen\"])\n",
    "\n",
    "# Un-standardize (add the mean and the std back)\n",
    "scenh = sdba.processing.unstandardize(scenh, savg, sstd)\n",
    "scens = sdba.processing.unstandardize(scens, savg, sstd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Restoring the trend\n",
    "\n",
    "The NpdfT has given us new \"hist\" and \"sim\" arrays with a correct rank structure. However, the trend is lost in this process. We reorder the result of the initial adjustment according to the rank structure of the NpdfT outputs to get our final bias-adjusted series.\n",
    "\n",
    "`sdba.processing.reordering` : 'ref' the argument that provides the order, 'sim' is the argument to reorder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scenh = sdba.processing.reordering(hist, scenh, group=\"time\")\n",
    "scens = sdba.processing.reordering(sim, scens, group=\"time\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scenh = sdba.processing.unstack_variables(scenh)\n",
    "scens = sdba.processing.unstack_variables(scens)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### There we are!\n",
    "\n",
    "Let's trigger all the computations. Here we write the data to disk and use `compute=False` in order to trigger the whole computation tree only once. There seems to be no way in xarray to do the same with a `load` call."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from dask import compute\n",
    "from dask.diagnostics import ProgressBar\n",
    "\n",
    "tasks = [\n",
    "    scenh.isel(location=2).to_netcdf(\"mbcn_scen_hist_loc2.nc\", compute=False),\n",
    "    scens.isel(location=2).to_netcdf(\"mbcn_scen_sim_loc2.nc\", compute=False),\n",
    "    extra.escores.isel(location=2)\n",
    "    .to_dataset()\n",
    "    .to_netcdf(\"mbcn_escores_loc2.nc\", compute=False),\n",
    "]\n",
    "\n",
    "with ProgressBar():\n",
    "    compute(tasks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the series and look at the distance scores to see how well the Npdf transform has converged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scenh = xr.open_dataset(\"mbcn_scen_hist_loc2.nc\")\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "dref.isel(location=2).tasmax.plot(ax=ax, label=\"Reference\")\n",
    "scenh.tasmax.plot(ax=ax, label=\"Adjusted\", alpha=0.65)\n",
    "dhist.isel(location=2).tasmax.plot(ax=ax, label=\"Simulated\")\n",
    "\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "escores = xr.open_dataarray(\"mbcn_escores_loc2.nc\")\n",
    "diff_escore = escores.differentiate(\"iterations\")\n",
    "diff_escore.plot()\n",
    "plt.title(\"Difference of the subsequent e-scores.\")\n",
    "plt.ylabel(\"E-scores difference\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff_escore"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tutorial continues in the [advanced notebook](sdba-advanced.ipynb) with more on optimization with dask, other fancier detrending algorithms and an example pipeline for heavy processing."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}