Add mapplot diagnostic to ClimWIP (#1864)

Co-authored-by: Peter Kalverla <peter.kalverla@gmx.com>

doc/sphinx/source/recipes/recipe_climwip.rst

@@ -58,7 +58,7 @@ User settings in recipe
 
   *Required settings for script*
     * ``ancestors``: must include weights from previous diagnostic
-    * ``weights``: the filename of the weights: 'weights_combined.nc'
+    * ``weights``: the filename of the weights: 'weights.nc'
 
   *Required settings for variables*
     * This script only takes temperature (tas) as input
@@ -75,6 +75,28 @@ User settings in recipe
     * Anomalies can be calculated with respect to a custom reference period
     * Monthly, annual or seasonal average/extraction can be used
 
+3. Script ``weighted_temperature_map.py``
+
+   *Required settings for script*
+     * ``ancestors``: must include weights from previous diagnostic
+     * ``weights``: the filename of the weights: 'weights_combined.nc'
+
+   *Optional settings for script*
+     * ``model_aggregation``: how to aggregate the models: mean (default), median, integer between 0 and 100 given a percentile
+     * ``xticks``: positions to draw xticks at
+     * ``yticks``: positions to draw yticks at
+
+   *Required settings for variables*
+     * This script takes temperature (tas) as input
+     * ``start_year``: provide the period for which to plot a temperature change graph.
+     * ``end_year``: provide the period for which to plot a temperature change graph.
+     * ``mip``: typically Amon
+     * ``preprocessor``: temperature_anomalies
+
+   *Optional settings for variables*
+     * A second variable is optional: temperature reference (tas_reference). If given, maps of temperature change to the reference are drawn, otherwise absolute temperature are drawn.
+     * tas_reference takes the same fields as tas
+
 
 Variables
 ---------
@@ -121,3 +143,9 @@ Example plots
    :align:   center
 
    Interquartile range of temperature anomalies relative to 1981-2010, weighted versus non-weighted.
+
+   .. _fig_climwip_5:
+.. figure::  /recipes/figures/climwip/temperature_change_weighted_map.png
+   :align:   center
+
+   Map of weighted mean temperature change 2081-2100 relative to 1995-2014

esmvaltool/diag_scripts/weighting/plot_utilities.py

@@ -0,0 +1,104 @@
+"""A collection of utility functions for dealing with weights."""
+from collections import defaultdict
+
+import numpy as np
+import xarray as xr
+
+
+def read_weights(filename: str) -> dict:
+    """Read a `.nc` file into a weights DataArray."""
+    weights_ds = xr.open_dataset(filename)
+    return weights_ds['weight']
+
+
+def read_metadata(cfg: dict, groupby: str = 'variable_group') -> dict:
+    """Read the metadata from the config file."""
+    datasets = defaultdict(list)
+
+    metadata = cfg['input_data'].values()
+
+    for item in metadata:
+        variable = item[groupby]
+
+        datasets[variable].append(item)
+
+    return datasets
+
+
+def weighted_quantile(values: list,
+                      quantiles: list,
+                      weights: list = None) -> 'np.array':
+    """Calculate weighted quantiles.
+
+    Analogous to np.quantile, but supports weights.
+
+    Based on: https://stackoverflow.com/a/29677616/6012085
+
+    Parameters
+    ----------
+    values: array_like
+        List of input values.
+    quantiles: array_like
+        List of quantiles between 0.0 and 1.0.
+    weights: array_like
+        List with same length as `values` containing the weights.
+
+    Returns
+    -------
+    np.array
+        Numpy array with computed quantiles.
+    """
+    values = np.array(values)
+    quantiles = np.array(quantiles)
+    if weights is None:
+        weights = np.ones(len(values))
+    weights = np.array(weights)
+
+    if not np.all((quantiles >= 0) & (quantiles <= 1)):
+        raise ValueError('Quantiles should be between 0.0 and 1.0')
+
+    idx = np.argsort(values)
+    values = values[idx]
+    weights = weights[idx]
+
+    weighted_quantiles = np.cumsum(weights) - 0.5 * weights
+
+    # Cast weighted quantiles to 0-1 To be consistent with np.quantile
+    min_val = weighted_quantiles.min()
+    max_val = weighted_quantiles.max()
+    weighted_quantiles = (weighted_quantiles - min_val) / max_val
+
+    return np.interp(quantiles, weighted_quantiles, values)
+
+
+def calculate_percentiles(data: 'xr.DataArray',
+                          percentiles: list,
+                          weights: dict = None) -> 'xr.DataArray':
+    """Calculate (weighted) percentiles.
+
+    Calculate the (weighted) percentiles for the given data.
+
+    Percentiles is a list of values between 0 and 100.
+
+    The `model_ensemble` dimension in weights has to contain at
+    least the same elements as in data.
+    If `weights` is not specified, the non-weighted percentiles are calculated.
+
+    Returns a DataArray with 'percentiles' as the dimension.
+    """
+    if weights is not None:
+        weights = weights.sel(model_ensemble=data.model_ensemble)
+
+    output = xr.apply_ufunc(weighted_quantile,
+                            data,
+                            input_core_dims=[['model_ensemble']],
+                            output_core_dims=[['percentiles']],
+                            kwargs={
+                                'weights': weights,
+                                'quantiles': percentiles / 100
+                            },
+                            vectorize=True)
+
+    output['percentiles'] = percentiles
+
+    return output

esmvaltool/diag_scripts/weighting/weighted_temperature_graph.py

@@ -3,92 +3,28 @@
 Lukas Brunner et al. section 2.4
 https://iopscience.iop.org/article/10.1088/1748-9326/ab492f
 """
-from collections import defaultdict
 import logging
 import os
 from pathlib import Path
 
 import matplotlib.pyplot as plt
 import numpy as np
 import xarray as xr
-from climwip import log_provenance, read_model_data
-from climwip import get_diagnostic_filename
 
+from climwip import get_diagnostic_filename, log_provenance, read_model_data
 from esmvaltool.diag_scripts.shared import get_plot_filename, run_diagnostic
+from esmvaltool.diag_scripts.weighting.plot_utilities import (
+    calculate_percentiles,
+    read_metadata,
+    read_weights,
+)
 
 logger = logging.getLogger(os.path.basename(__file__))
 
 
-def read_weights(filename: str) -> dict:
-    """Read a `.nc` file into a weights DataArray."""
-    weights_ds = xr.open_dataset(filename)
-    return weights_ds.to_dataframe().to_dict()['weight']
-
-
-def read_metadata(cfg: dict) -> dict:
-    """Read the metadata from the config file."""
-    datasets = defaultdict(list)
-
-    metadata = cfg['input_data'].values()
-
-    for item in metadata:
-        variable = item['short_name']
-
-        datasets[variable].append(item)
-
-    return datasets
-
-
-def weighted_quantile(values: list,
-                      quantiles: list,
-                      weights: list = None) -> 'np.array':
-    """Calculate weighted quantiles.
-
-    Analogous to np.quantile, but supports weights.
-
-    Based on: https://stackoverflow.com/a/29677616/6012085
-
-    Parameters
-    ----------
-    values: array_like
-        List of input values.
-    quantiles: array_like
-        List of quantiles between 0.0 and 1.0.
-    weights: array_like
-        List with same length as `values` containing the weights.
-
-    Returns
-    -------
-    np.array
-        Numpy array with computed quantiles.
-    """
-    values = np.array(values)
-    quantiles = np.array(quantiles)
-    if weights is None:
-        weights = np.ones(len(values))
-    weights = np.array(weights)
-
-    if not np.all((quantiles >= 0) & (quantiles <= 1)):
-        raise ValueError('Quantiles should be between 0.0 and 1.0')
-
-    idx = np.argsort(values)
-    values = values[idx]
-    weights = weights[idx]
-
-    weighted_quantiles = np.cumsum(weights) - 0.5 * weights
-
-    # Cast weighted quantiles to 0-1 To be consistent with np.quantile
-    min_val = weighted_quantiles.min()
-    max_val = weighted_quantiles.max()
-    weighted_quantiles = (weighted_quantiles - min_val) / max_val
-
-    return np.interp(quantiles, weighted_quantiles, values)
-
-
 def visualize_and_save_temperatures(temperature: 'xr.DataArray',
                                     iqr: 'xr.DataArray',
-                                    iqr_weighted: 'xr.DataArray',
-                                    cfg: dict,
+                                    iqr_weighted: 'xr.DataArray', cfg: dict,
                                     ancestors: list):
     """Visualize weighted temperature."""
     figure, axes = plt.subplots(dpi=300)
@@ -145,57 +81,24 @@ def plot_shaded(xrange, upper, lower, color, **kwargs):
     figure.savefig(filename_plot, dpi=300, bbox_inches='tight')
     plt.close(figure)
 
-    filename_data = get_diagnostic_filename(
-        'temperature_anomalies', cfg, extension='nc')
+    filename_data = get_diagnostic_filename('temperature_anomalies',
+                                            cfg,
+                                            extension='nc')
     temperature.to_netcdf(filename_data)
 
     caption = 'Temperature anomaly relative to 1981-2010'
     log_provenance(caption, filename_plot, cfg, ancestors)
     log_provenance(caption, filename_data, cfg, ancestors)
 
 
-def calculate_percentiles(data: 'xr.DataArray',
-                          percentiles: list,
-                          weights: dict = None):
-    """Calculate (weighted) percentiles.
-
-    Calculate the (weighted) percentiles for the given data.
-
-    Percentiles is a list of values between 0 and 100.
-
-    Weights is a dictionary where the keys represent the model members,
-    and the values the weights.
-    If `weights` is not specified, the non-weighted percentiles are calculated.
-
-    Returns a DataArray with 'percentiles' as the dimension.
-    """
-    if weights:
-        weights = [weights[member] for member in data.model_ensemble.values]
-
-    output = xr.apply_ufunc(weighted_quantile,
-                            data,
-                            input_core_dims=[['model_ensemble']],
-                            output_core_dims=[['percentiles']],
-                            kwargs={
-                                'weights': weights,
-                                'quantiles': percentiles / 100
-                            },
-                            vectorize=True)
-
-    output['percentiles'] = percentiles
-
-    return output
-
-
 def main(cfg):
     """Plot weighted temperature graph."""
     input_files = cfg['input_files']
     filename = cfg['weights']
     weights_path = Path(input_files[0]) / filename
     weights = read_weights(weights_path)
 
-    models = read_metadata(cfg)['tas']
-
+    models = read_metadata(cfg, 'short_name')['tas']
     model_data, model_data_files = read_model_data(models)
 
     percentiles = np.array([25, 75])