diff --git a/doc/sphinx/source/recipes/figures/climwip/temperature_change_weighted_map.png b/doc/sphinx/source/recipes/figures/climwip/temperature_change_weighted_map.png new file mode 100644 index 0000000000..be817ad0be Binary files /dev/null and b/doc/sphinx/source/recipes/figures/climwip/temperature_change_weighted_map.png differ diff --git a/doc/sphinx/source/recipes/recipe_climwip.rst b/doc/sphinx/source/recipes/recipe_climwip.rst index 2219eced90..5f29b0a093 100644 --- a/doc/sphinx/source/recipes/recipe_climwip.rst +++ b/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 diff --git a/esmvaltool/diag_scripts/weighting/plot_utilities.py b/esmvaltool/diag_scripts/weighting/plot_utilities.py new file mode 100644 index 0000000000..6a5101a673 --- /dev/null +++ b/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 diff --git a/esmvaltool/diag_scripts/weighting/weighted_temperature_graph.py b/esmvaltool/diag_scripts/weighting/weighted_temperature_graph.py index 7cc493ac21..b66a80f908 100644 --- a/esmvaltool/diag_scripts/weighting/weighted_temperature_graph.py +++ b/esmvaltool/diag_scripts/weighting/weighted_temperature_graph.py @@ -3,7 +3,6 @@ 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 @@ -11,84 +10,21 @@ 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,8 +81,9 @@ 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' @@ -154,39 +91,6 @@ def plot_shaded(xrange, upper, lower, color, **kwargs): 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'] @@ -194,8 +98,7 @@ def main(cfg): 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]) diff --git a/esmvaltool/diag_scripts/weighting/weighted_temperature_map.py b/esmvaltool/diag_scripts/weighting/weighted_temperature_map.py new file mode 100644 index 0000000000..3833b22ef0 --- /dev/null +++ b/esmvaltool/diag_scripts/weighting/weighted_temperature_map.py @@ -0,0 +1,171 @@ +"""Implementation of a mapplot for the climwip weighting scheme. + +Lukas Brunner et al. section 2.4 +https://iopscience.iop.org/article/10.1088/1748-9326/ab492f +""" +import logging +import os +from pathlib import Path + +import cartopy.crs as ccrs +import matplotlib.pyplot as plt +import numpy as np +from cartopy.mpl.ticker import LatitudeFormatter, LongitudeFormatter +from climwip import ( + get_diagnostic_filename, + get_plot_filename, + log_provenance, + read_model_data, +) + +from esmvaltool.diag_scripts.shared import run_diagnostic +from esmvaltool.diag_scripts.weighting.plot_utilities import ( + calculate_percentiles, + read_metadata, + read_weights, +) + +logger = logging.getLogger(os.path.basename(__file__)) + + +def mapplot(dataarray, cfg, title_pattern, filename_part, ancestors, + **colormesh_args): + """Visualize weighted temperature.""" + period = '{start_year}-{end_year}'.format(**read_metadata(cfg)['tas'][0]) + if 'tas_reference' in read_metadata(cfg).keys(): + meta = read_metadata(cfg)['tas_reference'] + period = 'change: {} minus {start_year}-{end_year}'.format( + period, **meta[0]) + metric = cfg['model_aggregation'] + if isinstance(metric, int): + metric = f'{metric}perc' + proj = ccrs.PlateCarree(central_longitude=0) + figure, axes = plt.subplots(subplot_kw={'projection': proj}) + + dataarray.plot.pcolormesh( + ax=axes, + transform=ccrs.PlateCarree(), + levels=9, + robust=True, + extend='both', + **colormesh_args + # colorbar size often does not fit nicely + # https://stackoverflow.com/questions/18195758/set-matplotlib-colorbar-size-to-match-graph + # cbar_kwargs={'fraction': .021} + ) + + lons = dataarray.lon.values + lats = dataarray.lat.values + longitude_formatter = LongitudeFormatter() + latitude_formatter = LatitudeFormatter() + default_xticks = np.arange(np.floor(lons.min()), np.ceil(lons.max()), 10) + default_yticks = np.arange(np.floor(lats.min()), np.ceil(lats.max()), 10) + + axes.coastlines() + axes.set_xticks(cfg.get('xticks', default_xticks), crs=proj) + axes.set_yticks(cfg.get('yticks', default_yticks), crs=proj) + axes.xaxis.set_ticks_position('both') + axes.yaxis.set_ticks_position('both') + axes.xaxis.set_major_formatter(longitude_formatter) + axes.yaxis.set_major_formatter(latitude_formatter) + axes.set_xlabel('') + axes.set_ylabel('') + + title = title_pattern.format(metric=metric, period=period) + axes.set_title(title) + + filename_plot = get_plot_filename(filename_part, cfg) + figure.savefig(filename_plot, dpi=300, bbox_inches='tight') + plt.close(figure) + + filename_data = get_diagnostic_filename(filename_part, cfg, extension='nc') + dataarray.to_netcdf(filename_data) + + log_provenance(title, filename_plot, cfg, ancestors) + log_provenance(title, filename_data, cfg, ancestors) + + +def visualize_and_save_temperature(temperature, cfg: dict, ancestors: list): + """Wrap mapplot: absolute temperature.""" + title_pattern = ('Weighted {metric} temperature\n' + r'{period} ($\degree$C)') + filename_part = 'temperature_change_weighted_map' + mapplot(temperature, + cfg, + title_pattern, + filename_part, + ancestors, + cmap='Reds') + + +def visualize_and_save_difference(temperature_difference, cfg: dict, + ancestors: list): + """Wrap mapplot: temperature difference between weighted and unweighted.""" + title_pattern = '\n'.join([ + 'Difference: weighted minus unweighted {metric} temperature', + r'{period} ($\degree$C)', + ]) + filename_part = 'temperature_change_difference_map' + mapplot(temperature_difference, + cfg, + title_pattern, + filename_part, + ancestors, + center=0) + + +def model_aggregation(dataset, metric, weights=None): + """Call mean or percentile calculation.""" + if isinstance(metric, int): + return calculate_percentiles(dataset, [metric], + weights).squeeze('percentile', drop=True) + if metric.lower() == 'mean': + if weights is not None: + dataset = dataset.weighted(weights) + return dataset.mean('model_ensemble') + + if metric.lower() == 'median': + return calculate_percentiles(dataset, [50], + weights).squeeze('percentile', drop=True) + + errmsg = f'model_aggregation {metric} is not implemented!' + raise NotImplementedError(errmsg) + + +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'] + model_data, model_data_files = read_model_data(models) + + # if a historical period is given calculate the change + if 'tas_reference' in read_metadata(cfg).keys(): + models_hist = read_metadata(cfg)['tas_reference'] + model_data_hist, model_data_files_hist = read_model_data(models_hist) + model_data_files += model_data_files_hist + model_data = model_data - model_data_hist + + metric = cfg.get('model_aggregation', 'mean') + unweighted_mean = model_aggregation(model_data, metric) + weighted_mean = model_aggregation(model_data, metric, weights) + + visualize_and_save_temperature( + weighted_mean, + cfg, + model_data_files, + ) + + visualize_and_save_difference( + weighted_mean - unweighted_mean, + cfg, + model_data_files, + ) + + +if __name__ == '__main__': + with run_diagnostic() as config: + main(config) diff --git a/esmvaltool/recipes/recipe_climwip.yml b/esmvaltool/recipes/recipe_climwip.yml index e7b20f41a4..59f603ff58 100644 --- a/esmvaltool/recipes/recipe_climwip.yml +++ b/esmvaltool/recipes/recipe_climwip.yml @@ -102,3 +102,25 @@ diagnostics: script: weighting/weighted_temperature_graph.py ancestors: [calculate_weights_climwip/climwip, tas] weights: 'weights.nc' + + weighted_temperature_map: + variables: + tas: &map_settings + short_name: tas + start_year: 2081 + end_year: 2100 + mip: Amon + preprocessor: climatological_mean + tas_reference: + <<: *map_settings + start_year: 1995 + end_year: 2014 + scripts: + weighted_temperature_map: + script: weighting/weighted_temperature_map.py + ancestors: [calculate_weights_climwip/climwip, tas, tas_reference] + weights: 'weights.nc' + # optional arguments + model_aggregation: mean # [ mean (default) | median | integer in (0, 100) ] + xticks: [0, 10, 20, 30, 40] # if not given ticks will be set automatically + yticks: [30, 40, 50, 60, 70, 80]