Refactor robustness methods

CHANGES.rst

@@ -2,6 +2,14 @@
 Changelog
 =========
 
+v0.47.0 (unreleased)
+--------------------
+Contributors to this version: Pascal Bourgault (:user:`aulemahal`).
+
+New features and enhancements
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+* New ``ipcc-advanced`` test added to ``xclim.ensembles.change_significance``.
+
 v0.46.0 (2023-10-24)
 --------------------
 Contributors to this version: David Huard (:user:`huard`), Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Éric Dupuis (:user:`coxipi`).

docs/notebooks/ensembles.ipynb

@@ -275,6 +275,188 @@
     "ax.legend()\n",
     "plt.show()"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Change significance and model agreement\n",
+    "\n",
+    "Aside from arithmetic reductions of the ensemble like shown above, xclim also implements helper functions for characterizing the significance of simulated climatic change and whether the different members agree or not over the sign of the change. See [xc.ensembles.change_significance](../apidoc/xclim.ensembles.rstl#xclim.ensembles._robustness.change_significance). Here we will use that function to generate masks replicating the ones in the [IPCC Atlas](https://interactive-atlas.ipcc.ch/), based on the approaches outlined in Cross-Chapter Box 1 of the [IPCC Atlas (AR6, WG1)](https://doi.org/10.1017/9781009157896.021). Approach A concerns observational trend and is out-of-scope for this the ensembles module.\n",
+    "\n",
+    "First we generate some fake annual mean temperature data. `ref` is the data on the reference period and `fut` is a future projection. In our example, there will be 5 different members in the ensemble."
+   ]
+  },
+  {
+   "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": [
+    "#### Approach B: The simple method.\n",
+    "\n",
+    "Zones where there is no agreement on the sign of the change are to be hatched. We call the function without any change significance test."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "chg_frac, pos_frac = ensembles.change_significance(fut, ref, test=None)\n",
+    "agreement = (pos_frac > 0.8) | (pos_frac < 0.2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here, `chg_frac` is the fraction of models showing significant change. Since `test = None`, this array is all ones, significance was not tested. `pos_frac` is the fraction of members showing positive change. Usually, when less than 80% of the members agree on the sign of the change, the region is hatched."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(6, 4))\n",
+    "mean_delta.plot.contourf(ax=ax)\n",
+    "agreement.plot.contourf(\n",
+    "    ax=ax,\n",
+    "    levels=[-0.5, 0.5, 1.5],\n",
+    "    cmap=\"none\",\n",
+    "    hatches=[\"///\", None],\n",
+    "    add_colorbar=False,\n",
+    ")\n",
+    "ax.legend(\n",
+    "    handles=[\n",
+    "        Rectangle((0, 0), 2, 2, fill=False, hatch=\"///\", label=\"Low model agreement\")\n",
+    "    ],\n",
+    "    bbox_to_anchor=(0.0, 1.1),\n",
+    "    loc=\"upper left\",\n",
+    ");"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Approach C: The advanced alternative.\n",
+    "\n",
+    "Here, we define three categories : robust change, conflicting change and no significant change or robust change. Change is said to be significant if the delta exceed an internal variability. The IPCC document suggests two different ways for computing this variability threshold:\n",
+    "\n",
+    "- $\\gamma = \\sqrt{2}\\cdot1.645\\cdot\\sigma_{20yr}$ where $\\sigma_{20yr}$ is the standard deviation of 20-year means, computed from non-overlapping periods in the pre-industrical control runs, after detrending with a quadratic fit.\n",
+    "- $\\gamma = \\sqrt{\\frac{2}{20}}\\cdot1.645\\cdot\\sigma_{1yr}$ where $\\sigma_{1yr}$ is the interanual standard deviation computed over the linearly detrended reference period.\n",
+    "\n",
+    "The first method involves data that is not always easily available. When such data is at hand, the xclim function could be used with a custom computation of the $\\gamma$ (`gamma`) threshold:\n",
+    "\n",
+    "```python\n",
+    "chg_frac, pos_frac = ensembles.change_significance(fut, ref, test='threshold', abs_thresh=gamma)\n",
+    "```\n",
+    "\n",
+    "The second method is more accessible and xclim implements it directly. However, `pos_frac` output of xclim is the fraction of the members showing _significant_ positive change, while the IPCC defines the agreement separately from the change significance. So we run the function twice."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "chg_frac, _ = ensembles.change_significance(fut, ref, test=\"ipcc-advanced\")\n",
+    "_, pos_frac = ensembles.change_significance(fut, ref, test=None)\n",
+    "\n",
+    "# Fraction of agreement, i.e.: max of positive or negative fraction\n",
+    "agree_frac = xr.concat((chg_frac * pos_frac, chg_frac * (1 - pos_frac)), \"sign\").max(\n",
+    "    \"sign\"\n",
+    ")\n",
+    "\n",
+    "robust = xr.where((chg_frac >= 0.66) & (agree_frac >= 0.8), 1, 0)\n",
+    "nochange_notrobust = xr.where((chg_frac < 0.66), 2, 0)\n",
+    "conflicting = xr.where((chg_frac >= 0.66) & (agree_frac < 0.8), 3, 0)\n",
+    "\n",
+    "mask = robust + nochange_notrobust + conflicting"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(6, 4))\n",
+    "mean_delta.plot.contourf(ax=ax)\n",
+    "mask.plot.contourf(\n",
+    "    ax=ax,\n",
+    "    levels=[0.5, 1.5, 2.5, 3.5],\n",
+    "    cmap=\"none\",\n",
+    "    hatches=[None, \"\\\\\\\\\\\\\", \"xxx\"],\n",
+    "    add_colorbar=False,\n",
+    ")\n",
+    "ax.legend(\n",
+    "    handles=[\n",
+    "        Rectangle((0, 0), 2, 2, fill=False, hatch=h, label=lbl)\n",
+    "        for h, lbl in (\n",
+    "            [\"\\\\\\\\\\\\\", \"No significant change\"],\n",
+    "            [\"xxx\", \"Conflicting signals\"],\n",
+    "        )\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}
+}

xclim/ensembles/_robustness.py

@@ -17,6 +17,7 @@
 from packaging.version import Version
 
 from xclim.core.formatting import update_history
+from xclim.indices.generic import detrend
 
 
 def change_significance(  # noqa: C901
@@ -49,7 +50,7 @@ def change_significance(  # noqa: C901
         a distribution across the future period.
         `fut` and `ref` must be of the same type (Dataset or DataArray). If they are
         Dataset, they must have the same variables (name and coords).
-    test : {'ttest', 'welch-ttest', 'mannwhitney-utest', 'brownforsythe-test', 'threshold', None}
+    test : {'ttest', 'welch-ttest', 'mannwhitney-utest', 'brownforsythe-test', 'ipcc-advanced', 'threshold', None}
         Name of the statistical test used to determine if there was significant change. See notes.
     weights : xr.DataArray
         Weights to apply along the 'realization' dimension. This array cannot contain missing values.
@@ -111,13 +112,19 @@ def change_significance(  # noqa: C901
         Brown-Forsythe test assuming skewed, non-normal distributions. Same significance criterion as 'ttest'.
       'threshold' :
         Change is considered significative if the absolute delta exceeds a given threshold (absolute or relative).
+      'ipcc-advanced' :
+        An approximation of the "advanced approach" used in the IPCC Atlas and described in :cite:t:`ipccatlas_ar6wg1`.
+        Change is considered significant if the delta exceeds :math:`1.645\sqrt{\frac{2}{10}}\sigma_{ref'}$ where $\sigma_{ref'}` is
+        the interannual standard deviation of the linearly detrended reference (`ref`). See notebook :ref:`notebooks/ensembles:Ensembles`
+        for more details.
       None :
-        Significant change is not tested and, thus, members showing no change are
-        included in the `sign_frac` output.
+        Significant change is not tested. Members showing any non-zero positive change are
+        included in the `pos_frac` output.
 
     References
     ----------
     :cite:cts:`tebaldi_mapping_2011`
+    :cite:cts:`ipccatlas_ar6wg1`
 
     Example
     -------
@@ -167,6 +174,7 @@ def change_significance(  # noqa: C901
         "mannwhitney-utest": ["p_change"],
         "brownforsythe-test": ["p_change"],
         "threshold": ["abs_thresh", "rel_thresh"],
+        "ipcc-advanced": [],
     }
 
     # Get delta, either from fut or from fut - ref
@@ -323,7 +331,12 @@ def mwu_wrapper(f, r):  # This specific test can't manage an all-NaN slice
             changed = abs(delta / ref.mean("time")) > kwargs["rel_thresh"]
         else:
             raise ValueError("Invalid argument combination for test='threshold'.")
-
+    elif test == "ipcc-advanced":
+        # Detrend ref
+        refy = ref.resample(time="YS").mean()
+        ref_detrended = detrend(refy, dim="time", deg=1)
+        gamma = np.sqrt(2 / 20) * 1.645 * ref_detrended.std("time")
+        changed = abs(delta) > gamma
     elif test is not None:
         raise ValueError(
             f"Statistical test {test} must be one of {', '.join(test_params.keys())}."

xclim/indices/generic.py

@@ -40,6 +40,7 @@
     "count_occurrences",
     "cumulative_difference",
     "default_freq",
+    "detrend",
     "diurnal_temperature_range",
     "domain_count",
     "doymax",
@@ -1004,3 +1005,30 @@ def _get_zone(da):
         )
 
     return zones
+
+
+def detrend(ds, dim="time", deg=1):
+    """Detrend data along a given dimension computing a polynomial trend of a given order.
+
+    Parameters
+    ----------
+    ds : xr.Dataset or xr.DataArray
+      The data to detrend. If a Dataset, detrending is done on all data variables.
+    dim : str
+      Dimension along which to compute the trend.
+    deg : int
+      Degree of the polynomial to fit.
+
+    Returns
+    -------
+    detrended : xr.Dataset or xr.DataArray
+      Same as `ds`, but with its trend removed (subtracted).
+    """
+    if isinstance(ds, xr.Dataset):
+        return ds.map(detrend, keep_attrs=False, dim=dim, deg=deg)
+    # is a DataArray
+    # detrend along a single dimension
+    coeff = ds.polyfit(dim=dim, deg=deg)
+    trend = xr.polyval(ds[dim], coeff.polyfit_coefficients)
+    with xr.set_options(keep_attrs=True):
+        return ds - trend