Refactor robustness methods (#1514)

<!--Please ensure the PR fulfills the following requirements! -->
<!-- If this is your first PR, make sure to add your details to the
AUTHORS.rst! -->
### Pull Request Checklist:
- [ ] This PR addresses an already opened issue (for bug fixes /
features)
    - This PR fixes #xyz
- [ ] Tests for the changes have been added (for bug fixes / features)
- [x] (If applicable) Documentation has been added / updated (for bug
fixes / features)
- [x] CHANGES.rst has been updated (with summary of main changes)
- [ ] Link to issue (:issue:`number`) and pull request (:pull:`number`)
has been added

### What kind of change does this PR introduce?

* Add `ipcc-advanced` as a test to `ensembles.change_significance`. This
corresponds to the simpler alternative of approach C as described in the
Cross-Chapter Box 1 of the [IPCC Atlas (AR6,
WG1)](https://www.cambridge.org/core/books/climate-change-2021-the-physical-science-basis/atlas/24E1C016DBBE4725BDFBC343695DE7DB).
* Add two examples in the "Ensembles" notebook to show how xclim can be
used to get the IPCC-recommended hatching masks.
* Add a very simple `detrend` function.

### Does this PR introduce a breaking change?
No, but it might, see below.

### Other information:
As you can see in the examples, the IPCC decouples the "change
significance" from the "sign agreement", while xclim does not. This make
the second example more complex than expected.

I thought of :
- Changing the output of `change_significance` so that `pos_frac` is
independent of change significance
- Adding a third output so we have: `significant_change_frac,
significant_positive_change_frac, positive_change_frac`.
- Adding a keyword argument to switch between the two modes.
First two are breaking changes while the thirds feels overly heavy for
something this simple...

CHANGES.rst

@@ -4,7 +4,11 @@ Changelog
 
 v0.47.0 (unreleased)
 --------------------
-Contributors to this version: Juliette Lavoie (:user:`juliettelavoie`), Trevor James Smith (:user:`Zeitsperre`).
+Contributors to this version: Juliette Lavoie (:user:`juliettelavoie`), Pascal Bourgault (:user:`aulemahal`), Trevor James Smith (:user:`Zeitsperre`).
+
+New features and enhancements
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+* New functions ``xclim.ensembles.robustness_fractions`` and ``xclim.ensembles.robustness_categories``. The former will replace ``xclim.ensembles.change_significance`` which is now deprecated and will be removed in xclim 0.49 (:pull:`1514`).
 
 Bug fixes
 ^^^^^^^^^

docs/notebooks/ensembles.ipynb

@@ -164,6 +164,7 @@
     "# Set display to HTML style (for fancy output)\n",
     "xr.set_options(display_style=\"html\", display_width=50)\n",
     "\n",
+    "import matplotlib as mpl\n",
     "import matplotlib.pyplot as plt\n",
     "\n",
     "%matplotlib inline\n",
@@ -275,6 +276,187 @@
     "ax.legend()\n",
     "plt.show()"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Change significance and model agreement\n",
+    "\n",
+    "When communicating climate change through plots of projected change, it is often useful to add information on the statistical significance of the values. A common way to represent this information without overloading the figures is through hatching patterns superimposed on the primary data. Two aspects are usually shown : \n",
+    "\n",
+    "- change significance : whether most of the ensemble members project a statistically significant climate change signal, in comparison to their internal variability.\n",
+    "- model agreement : whether the different ensemble members agree on the sign of the change.\n",
+    "\n",
+    "We can then divide the plotted points into categories each with its own hatching pattern, usually leaving the robust data (models agree and enough show a significant change) without hatching. \n",
+    "\n",
+    "Xclim provides some tools to help in generating these hatching masks. First is [xc.ensembles.robustness_fractions](../apidoc/xclim.ensembles.rst#xclim.ensembles._robustness.robustness_fractions) that can characterize the change significance and sign agreement accross ensemble members. To demonstrate its usage, we'll first generate some fake annual mean temperature data. Here, `ref` is the data on the reference period and `fut` is a future projection. There are 5 different members in the ensemble. We tweaked the generation so that all models agree on significant change in the \"south\" while agreement and signifiance of change decreases as we go north and east."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import xarray as xr\n",
+    "from matplotlib.patches import Rectangle\n",
+    "\n",
+    "xr.set_options(keep_attrs=True)\n",
+    "\n",
+    "# Reference period\n",
+    "ref = xr.DataArray(\n",
+    "    20 * np.random.random_sample((5, 30, 10, 10)) + 275,\n",
+    "    dims=(\"realization\", \"time\", \"lat\", \"lon\"),\n",
+    "    coords={\n",
+    "        \"time\": xr.date_range(\"1990\", periods=30, freq=\"YS\"),\n",
+    "        \"lat\": np.arange(40, 50),\n",
+    "        \"lon\": np.arange(-70, -60),\n",
+    "    },\n",
+    "    attrs={\"units\": \"K\"},\n",
+    ")\n",
+    "\n",
+    "# Future\n",
+    "fut = xr.DataArray(\n",
+    "    20 * np.random.random_sample((5, 30, 10, 10)) + 275,\n",
+    "    dims=(\"realization\", \"time\", \"lat\", \"lon\"),\n",
+    "    coords={\n",
+    "        \"time\": xr.date_range(\"2070\", periods=30, freq=\"YS\"),\n",
+    "        \"lat\": np.arange(40, 50),\n",
+    "        \"lon\": np.arange(-70, -60),\n",
+    "    },\n",
+    "    attrs={\"units\": \"K\"},\n",
+    ")\n",
+    "# Add change.\n",
+    "fut = fut + xr.concat(\n",
+    "    [\n",
+    "        xr.DataArray(np.linspace(15, north_delta, num=10), dims=(\"lat\",))\n",
+    "        for north_delta in [15, 10, 0, -7, -10]\n",
+    "    ],\n",
+    "    \"realization\",\n",
+    ")\n",
+    "\n",
+    "deltas = (fut.mean(\"time\") - ref.mean(\"time\")).assign_attrs(\n",
+    "    long_name=\"Temperature change\"\n",
+    ")\n",
+    "mean_delta = deltas.mean(\"realization\")\n",
+    "deltas.plot(col=\"realization\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Change significance can be determined in a lot of different ways. Xclim provides some simple and some more complicated statistical test in `robustness_fractions`. In this example, we'll follow the suggestions found in the Cross-Chapter Box 1 of the [IPCC Atlas chapter (AR6, WG1)](https://doi.org/10.1017/9781009157896.021). Specifically, we are following Approach C, using the alternative for when pre-industrial control data is not available.\n",
+    "\n",
+    "We first compute the different fractions for each robustness aspect."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fractions = ensembles.robustness_fractions(fut, ref, test=\"ipcc-ar6-c\")\n",
+    "fractions"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this output we have:\n",
+    "\n",
+    "- `changed` : The fraction of members showing significant change.\n",
+    "- `positive` : The fraction of members showing positive change, no matter if it is significant or not.\n",
+    "- `changed_positive` : The fraction of members showing significant AND positive change.\n",
+    "- `agree` : The fraction of members agreeing on the sign of change. This is the maximum between `positive` and `1 - positive`.\n",
+    "- `valid` : The fraction of \"valid\" members. A member is valid is there are no NaNs along the time axes of `fut` and  `ref`. In our case, it is 1 everywhere.\n",
+    "\n",
+    "For example, here's the plot of the fraction of members showing significant change."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fractions.changed.plot(figsize=(6, 4))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Xclim provides all this so that one can construct their own robustness maps the way they want. Often, hatching overlays are based on categories defined by some thresholds on the significant change and agreement fractions. The [`xclim.ensembles.robustness_categories`](../apidoc/xclim.ensembles.rst#xclim.ensembles._robustness.robustness_categories) function helps for that common case and defaults to the categories and thresholds used by the IPCC in its Atlas."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "robustness = ensembles.robustness_categories(fractions)\n",
+    "robustness"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The output is a categorical map following the \"flag variables\" CF conventions. Parameters needed for plotting are found in the attributes. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "robustness.plot(figsize=(6, 4))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Matplotlib doesn't provide an easy way of plotting categorial data with a proper legend, so our real plotting script is a bit more complicated, but xclim's output makes it easier."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "cmap = mpl.colors.ListedColormap([\"none\"])  # So we can deactivate pcolor's colormapping\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(6, 4))\n",
+    "mean_delta.plot(ax=ax)\n",
+    "# For each flag value plot the corresponding hatch.\n",
+    "for val, ha in zip(robustness.flag_values, [None, \"\\\\\\\\\\\\\", \"xxx\"]):\n",
+    "    ax.pcolor(\n",
+    "        robustness.lon,\n",
+    "        robustness.lat,\n",
+    "        robustness.where(robustness == val),\n",
+    "        hatch=ha,\n",
+    "        cmap=cmap,\n",
+    "    )\n",
+    "\n",
+    "ax.legend(\n",
+    "    handles=[\n",
+    "        Rectangle((0, 0), 2, 2, fill=False, hatch=h, label=lbl)\n",
+    "        for h, lbl in zip([\"\\\\\\\\\\\\\", \"xxx\"], robustness.flag_descriptions[1:])\n",
+    "    ],\n",
+    "    bbox_to_anchor=(0.0, 1.1),\n",
+    "    loc=\"upper left\",\n",
+    "    ncols=2,\n",
+    ");"
+   ]
   }
  ],
  "metadata": {

docs/references.bib

@@ -2073,3 +2073,16 @@ @article{tobin_2018
 title = {Vulnerabilities and resilience of European power generation to 1.5°C, 2°C and 3°C warming},
 journal = {Environmental Research Letters}
 }
+
+
+@inbook{
+ipccatlas_ar6wg1,
+place={Cambridge},
+title={Atlas},
+DOI={10.1017/9781009157896.021},
+booktitle={Climate Change 2021 – The Physical Science Basis: Working Group I Contribution to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change},
+publisher={Cambridge University Press},
+author={Intergovernmental Panel on Climate Change (IPCC)},
+year={2023},
+pages={1927–2058}
+}