Merge 6976fec into 67d806c

.zenodo.json

@@ -28,6 +28,11 @@
       "name": "Braun, Marco",
       "affiliation": "Ouranos",
       "orcid": "0000-0001-5061-3217"
+    },
+    {
+      "name": "Dupuis, Éric",
+      "affiliation": "Ouranos, Montréal, Québec, Canada",
+      "orcid": "0000-0001-7976-4596"
     }
   ],
   "keywords": [

AUTHORS.rst

@@ -20,3 +20,4 @@ Contributors
 
 * Gabriel Rondeau-Genesse <rondeau-genesse.gabriel@ouranos.ca> `@RondeauG <https://github.com/RondeauG>`_
 * Marco Braun <Braun.Marco@ouranos.ca> `@vindelico <https://github.com/vindelico>`_
+* Éric Dupuis <dupuis.eric@ouranos.ca> `@coxipi <https://github.com/coxipi>`_

CHANGELOG.rst

@@ -4,7 +4,7 @@ Changelog
 
 0.4.0 (unreleased)
 ------------------
-Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Marco Braun (:user:`vindelico`), Pascal Bourgault (:user:`aulemahal`), Sarah-Claude Bourdeau-Goulet (:user:`Sarahclaude`)
+Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Marco Braun (:user:`vindelico`), Pascal Bourgault (:user:`aulemahal`), Sarah-Claude Bourdeau-Goulet (:user:`Sarahclaude`), Éric Dupuis (:user:`coxipi`)
 
 New features and enhancements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -13,13 +13,16 @@ New features and enhancements
 * Added style sheet ``transparent.mplstyle`` (:issue:`183`, :pull:`185`)
 * Fix ``NaN`` issues, extreme values in sizes legend and added edgecolors in ``fg.matplotlib.scattermap``  (:pull:`184`).
 * New function ``fg.data`` for fetching package data and defined `matplotlib` style definitions. (:pull:`211`).
+* ``fg.taylordiagram`` can now accept datasets with many dimensions (not only `taylor_params`), provided that they all share the same `ref_std` (e.g. normalized taylor diagrams)  (:pull:`214`).
+* A new optional way to organize points in a `fg.taylordiagram`  with  `colors_key`, `markers_key`  : DataArrays with a common dimension value or a common attrtibute are grouped with the same color/marker (:pull:`214`).
 
 Breaking changes
 ^^^^^^^^^^^^^^^^
 * `figanos` no longer supports Python 3.8. (:pull:`210`).
 * `figanos` now uses a `'src' layout <https://packaging.python.org/en/latest/discussions/src-layout-vs-flat-layout>`_ for the package. (:pull:`210`).
 * `cartopy` has been pinned above v0.23.0 due to a licensing issue. (:pull:`210`).
 * `twine` and `wheel` have been removed from the `dev` requirements. (:pull:`210`).
+* ``fg.taylordiagram`` returns a tuple of `(fig, floating_ax, legend)` instead of only `floating_ax`. (:pull:`214`).
 
 Internal changes
 ^^^^^^^^^^^^^^^^

docs/notebooks/figanos_multiplots.ipynb

@@ -199,6 +199,83 @@
     "im.fig.suptitle(\"Multiple hatchmaps\", y=1.08)\n"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Taylor diagram\n",
+    "\n",
+    "A Taylor diagram displays the correlation along a given dimension between a target and a source DataArray. The standard deviation of each DataArrays are also displayed. In the example below, maximum temperatures `tasmax` are compared along time."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from xclim.testing import open_dataset\n",
+    "from xclim import sdba\n",
+    "from xclim.core.units import convert_units_to\n",
+    "\n",
+    "ds = open_dataset(\"sdba/CanESM2_1950-2100.nc\").sel(time=slice(\"1950\", \"2013\"))\n",
+    "ds1 = open_dataset(\"sdba/ahccd_1950-2013.nc\")\n",
+    "ds2 = open_dataset(\"sdba/nrcan_1950-2013.nc\")\n",
+    "for v in ds.data_vars: \n",
+    "    ds1[v] = convert_units_to(ds1[v], ds[v], context=\"hydro\")\n",
+    "    ds2[v] = convert_units_to(ds2[v], ds[v], context=\"hydro\")\n",
+    "tx_ref = ds.tasmax.isel(location=0)\n",
+    "tx_sim = ds1.tasmax.isel(location=0)\n",
+    "tx_sim2 = ds2.tasmax.isel(location=0)\n",
+    "out = {}\n",
+    "out[\"sim\"] =  sdba.measures.taylordiagram(ref=tx_ref, sim=tx_sim, dim=\"time\")\n",
+    "out[\"sim2\"] = sdba.measures.taylordiagram(ref=tx_ref, sim=tx_sim2, dim=\"time\")\n",
+    "# In the default behaviour of `fg.taylordiagram`, the keys of the dictionary of \n",
+    "# DataArray are used for the legend\n",
+    "# A single DataArray can also be passed.\n",
+    "fg.taylordiagram(out)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Normalized taylor diagram\n",
+    "\n",
+    "If we normalize the standard deviation of our measures, many Taylor diagrams with difference references can be combined in a single plot. In the following example, we keep all three locations in the Datasets from the previous example, and we also include the precipitation instead of only selecting the maximum temperature. For each location (3) and variable (2), a Taylordiagram measure is computed. Each set of correlation and standard deviation is then plotted. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from xclim.testing import open_dataset\n",
+    "from xclim import sdba\n",
+    "from xclim.core.units import convert_units_to\n",
+    "import figanos.matplotlib as fg\n",
+    "\n",
+    "# Here, we have three locations, two variables. We stack variables to convert from\n",
+    "# a Dataset to a DataArray.\n",
+    "da_ref = sdba.stack_variables(ds)\n",
+    "da_sim = sdba.stack_variables(ds2)\n",
+    "\n",
+    "# Each location/variable will have its own set of taylor parameters\n",
+    "out = sdba.measures.taylordiagram(ref=da_ref, sim=da_sim, dim=\"time\")\n",
+    "\n",
+    "# If we normalize the taylor diagrams, they can be compared on the same plot\n",
+    "out[{\"taylor_param\":[0,1]}] = out[{\"taylor_param\":[0,1]}]/out[{\"taylor_param\":0}]\n",
+    "\n",
+    "# The `markers_key` and `colors_key` are used to separate between two different features. \n",
+    "# Here, the type of marker is used to distinguish between locations, and the color \n",
+    "# distinguishes between variables. If those parameters are not specified, then each \n",
+    "# pair (location, multivar)  has simply its own color.\n",
+    "# NB: Negatives correlations are not plotted\n",
+    "fg.taylordiagram(out, markers_key=\"location\", colors_key=\"multivar\")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -367,7 +444,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.0"
+   "version": "3.12.3"
   }
  },
  "nbformat": 4,

src/figanos/matplotlib/plot.py

@@ -1827,6 +1827,8 @@ def taylordiagram(
     legend_kw: dict[str, Any] | None = None,
     std_label: str | None = None,
     corr_label: str | None = None,
+    colors_key: str | None = None,
+    markers_key: str | None = None,
 ):
     """Build a Taylor diagram.
 
@@ -1854,16 +1856,26 @@ def taylordiagram(
         Label for the standard deviation (x and y) axes.
     corr_label : str, optional
         Label for the correlation axis.
+    colors_key : str, optional
+        Attribute or dimension of DataArrays used to separate DataArrays into groups with different colors. If present,
+        it overrides the "color" key in `plot_kw`.
+    markers_key : str, optional
+        Attribute or dimension of DataArrays used to separate DataArrays into groups with different markers. If present,
+        it overrides the "marker" key in `plot_kw`.
 
     Returns
     -------
-    matplotlib.axes.Axes
+    (plt.figure, mpl_toolkits.axisartist.floating_axes.FloatingSubplot, plt.legend)
     """
     plot_kw = empty_dict(plot_kw)
     fig_kw = empty_dict(fig_kw)
     contours_kw = empty_dict(contours_kw)
     legend_kw = empty_dict(legend_kw)
 
+    # preserve order of dimensions if used for marker/color
+    ordered_markers_type = None
+    ordered_colors_type = None
+
     # convert SSP, RCP, CMIP formats in keys
     if isinstance(data, dict):
         data = process_keys(data, convert_scen_name)
@@ -1875,8 +1887,7 @@ def taylordiagram(
         data = {"_no_label": data}  # mpl excludes labels starting with "_" from legend
         plot_kw = {"_no_label": empty_dict(plot_kw)}
     elif not plot_kw:
-        plot_kw = {}
-
+        plot_kw = {k: {} for k in data.keys()}
     # check type
     for key, v in data.items():
         if not isinstance(v, xr.DataArray):
@@ -1886,6 +1897,35 @@ def taylordiagram(
         if key == "reference":
             raise ValueError("'reference' is not allowed as a key in data.")
 
+    # If there are other dimensions than 'taylor_param', create a bigger dict with them
+    data_keys = list(data.keys())
+    for data_key in data_keys:
+        da = data[data_key]
+        dims = list(set(da.dims) - {"taylor_param"})
+        if dims != []:
+            if markers_key in dims:
+                ordered_markers_type = da[markers_key].values
+            if colors_key in dims:
+                ordered_colors_type = da[colors_key].values
+
+            da = da.stack(pl_dims=dims)
+            for i, dim_key in enumerate(da.pl_dims.values):
+                if isinstance(dim_key, list) or isinstance(dim_key, tuple):
+                    dim_key = "-".join([str(k) for k in dim_key])
+                da0 = da.isel(pl_dims=i)
+                # if colors_key/markers_key is a dimension, add it as an attribute for later use
+                if markers_key in dims:
+                    da0.attrs[markers_key] = da0[markers_key].values.item()
+                if colors_key in dims:
+                    da0.attrs[colors_key] = da0[colors_key].values.item()
+                new_data_key = (
+                    f"{data_key}-{dim_key}" if data_key != "_no_label" else dim_key
+                )
+                data[new_data_key] = da0
+                plot_kw[new_data_key] = empty_dict(plot_kw[f"{data_key}"])
+            data.pop(data_key)
+            plot_kw.pop(data_key)
+
     # remove negative correlations
     initial_len = len(data)
     removed = [
@@ -1929,10 +1969,9 @@ def taylordiagram(
     # make labels
     if not std_label:
         try:
-            std_label = (
-                get_localized_term("standard deviation")
-                + f" ({list(data.values())[0].units})"
-            )
+            units = list(data.values())[0].units
+            std_label = get_localized_term("standard deviation")
+            std_label = std_label if units == "" else f"{std_label} ({units})"
         except AttributeError:
             std_label = get_localized_term("standard deviation").capitalize()
 
@@ -1949,12 +1988,11 @@ def taylordiagram(
     transform = PolarAxes.PolarTransform()
 
     # Setup the axis, here we map angles in degrees to angles in radius
-    rlocs = np.array([0, 0.2, 0.4, 0.6, 0.8, 1])
-    rlocs_deg = rlocs * 90
-    tlocs = rlocs_deg * np.pi / 180  # convert degrees to radians
+    # Correlation labels
+    rlocs = np.array([0, 0.2, 0.4, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 1])
+    tlocs = np.arccos(rlocs)  # Conversion to polar angles
     gl1 = gf.FixedLocator(tlocs)  # Positions
-    tf1 = gf.DictFormatter(dict(zip(tlocs, np.flip(rlocs).astype(str))))
-
+    tf1 = gf.DictFormatter(dict(zip(tlocs, map(str, rlocs))))
     # Standard deviation axis extent
     radius_min = std_range[0] * max(max_std)
     radius_max = std_range[1] * max(max_std)
@@ -1968,7 +2006,6 @@ def taylordiagram(
     )
 
     fig = plt.figure(**fig_kw)
-
     floating_ax = FloatingSubplot(fig, 111, grid_helper=ghelper)
     fig.add_subplot(floating_ax)
 
@@ -2013,7 +2050,8 @@ def taylordiagram(
     # rmse contours from reference standard deviation
     if contours:
         radii, angles = np.meshgrid(
-            np.linspace(radius_min, radius_max), np.linspace(0, np.pi / 2)
+            np.linspace(radius_min, radius_max),
+            np.linspace(0, np.pi / 2),
         )
         # Compute centered RMS difference
         rms = np.sqrt(ref_std**2 + radii**2 - 2 * ref_std * radii * np.cos(angles))
@@ -2023,47 +2061,93 @@ def taylordiagram(
 
         ax.clabel(ct, ct.levels, fontsize=8)
 
-        ct_line = Line2D(
+        # points.append(ct_line)
+        ct_line = ax.plot(
             [0],
             [0],
             ls=contours_kw["linestyles"],
             lw=1,
             c="k" if "colors" not in contours_kw else contours_kw["colors"],
             label="rmse",
         )
-        points.append(ct_line)
+        points.append(ct_line[0])
 
     # get color options
     style_colors = matplotlib.rcParams["axes.prop_cycle"].by_key()["color"]
     if len(data) > len(style_colors):
         style_colors = style_colors * math.ceil(len(data) / len(style_colors))
     cat_colors = Path(__file__).parents[1] / "data/ipcc_colors/categorical_colors.json"
+    # get marker options (only used if `markers_key` is set)
+    style_markers = "oDv^<>p*hH+x|_"
+    if len(data) > len(style_markers):
+        style_markers = style_markers * math.ceil(len(data) / len(style_markers))
+
+    # set colors and markers styles based on discrimnating attributes (if specified)
+    if colors_key or markers_key:
+        if colors_key:
+            # get_scen_color : look for SSP, RCP, CMIP model color
+            colors_type = (
+                ordered_colors_type
+                if ordered_colors_type is not None
+                else {da.attrs[colors_key] for da in data.values()}
+            )
+            colorsd = {
+                k: get_scen_color(k, cat_colors) or style_colors[i]
+                for i, k in enumerate(colors_type)
+            }
+        if markers_key:
+            markers_type = (
+                ordered_markers_type
+                if ordered_markers_type is not None
+                else {da.attrs[markers_key] for da in data.values()}
+            )
+            markersd = {k: style_markers[i] for i, k in enumerate(markers_type)}
+
+        for key, da in data.items():
+            if colors_key:
+                plot_kw[key]["color"] = colorsd[da.attrs[colors_key]]
+            if markers_key:
+                plot_kw[key]["marker"] = markersd[da.attrs[markers_key]]
 
     # plot scatter
     for (key, da), i in zip(data.items(), range(len(data))):
         # look for SSP, RCP, CMIP model color
-        if get_scen_color(key, cat_colors):
-            plot_kw[key].setdefault("color", get_scen_color(key, cat_colors))
-        else:
-            plot_kw[key].setdefault("color", style_colors[i])
-
-        # convert corr to polar coordinates
-        plot_corr = (1 - da.sel(taylor_param="corr").values) * 90 * np.pi / 180
-
+        if colors_key is None:
+            plot_kw[key].setdefault(
+                "color", get_scen_color(key, cat_colors) or style_colors[i]
+            )
         # set defaults
         plot_kw[key] = {"label": key} | plot_kw[key]
 
+        # legend will be handled later in this case
+        if markers_key or colors_key:
+            plot_kw[key]["label"] = ""
+
         # plot
         pt = ax.scatter(
-            plot_corr, da.sel(taylor_param="sim_std").values, **plot_kw[key]
+            np.arccos(da.sel(taylor_param="corr").values),
+            da.sel(taylor_param="sim_std").values,
+            **plot_kw[key],
         )
         points.append(pt)
 
     # legend
     legend_kw.setdefault("loc", "upper right")
-    fig.legend(points, [pt.get_label() for pt in points], **legend_kw)
-
-    return floating_ax
+    legend = fig.legend(points, [pt.get_label() for pt in points], **legend_kw)
+
+    # plot new legend if markers/colors represent a certain dimension
+    if colors_key or markers_key:
+        handles = list(floating_ax.get_legend_handles_labels()[0])
+        if markers_key:
+            for k, m in markersd.items():
+                handles.append(Line2D([0], [0], color="k", label=k, marker=m, ls=""))
+        if colors_key:
+            for k, c in colorsd.items():
+                handles.append(Line2D([0], [0], color=c, label=k, ls="-"))
+        legend.remove()
+        legend = fig.legend(handles=handles, **legend_kw)
+
+    return fig, floating_ax, legend
 
 
 def hatchmap(