From 813b11ce15d851c60839cc544cd1a13a4b201033 Mon Sep 17 00:00:00 2001 From: Elizaveta Malinina Date: Wed, 29 Jan 2025 00:03:58 +0000 Subject: [PATCH 01/12] add style file for malinina24 --- .../shared/plot/styles_python/malinina24.yml | 51 +++++++++++++++++++ 1 file changed, 51 insertions(+) create mode 100644 esmvaltool/diag_scripts/shared/plot/styles_python/malinina24.yml diff --git a/esmvaltool/diag_scripts/shared/plot/styles_python/malinina24.yml b/esmvaltool/diag_scripts/shared/plot/styles_python/malinina24.yml new file mode 100644 index 0000000000..a6631f3a2f --- /dev/null +++ b/esmvaltool/diag_scripts/shared/plot/styles_python/malinina24.yml @@ -0,0 +1,51 @@ +############################################################################### +# PYTHON STYLE FILE +# malinina24.yml +############################################################################### +# This file defines plot attributes for certain datasets and groups. +# Template: +# +# DATASET: +# avgstd: 0 +# color: '#ffffff' +# dash: -- +# facecolor: none +# mark: x +# thick: 1 +################################################################################ + +ERA5: + avgstd: 1 + color: '#cd5c5c' + dash: '-' + facecolor: '#cd5c5c' + mark: 'o' + thick: 1 +nat: + avgstd: 1 + color: '#005000' + dash: '-' + facecolor: '#005000' + mark: 'o' + thick: 1 +factual: + avgstd: 1 + color: '#c57a00' + dash: '-' + facecolor: '#c57a00' + mark: 'o' + thick: 1 +ssp245: + avgstd: 1 + color: '#4677c0' + dash: '-' + facecolor: '#4677c0' + mark: 'o' + thick: 1 +default: + avgstd: 0 + color: '#1A0F50' + dash: '-' + facecolor: '#1A0F50' + mark: o + thick: 1 \ No newline at end of file From 2ad0dc32b50923b75d62c8870fcdd957e1e5f83e Mon Sep 17 00:00:00 2001 From: Elizaveta Malinina Date: Wed, 29 Jan 2025 00:11:33 +0000 Subject: [PATCH 02/12] add reference --- esmvaltool/references/malinina24.bibtex | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 esmvaltool/references/malinina24.bibtex diff --git a/esmvaltool/references/malinina24.bibtex b/esmvaltool/references/malinina24.bibtex new file mode 100644 index 0000000000..f5fd7d7437 --- /dev/null +++ b/esmvaltool/references/malinina24.bibtex @@ -0,0 +1,10 @@ +@article{malinina24, + title={The 2021 heatwave was less rare in Western Canada than previously thought}, + author={Malinina, Elizaveta and Gillett, Nathan P}, + journal={Weather and Climate Extremes}, + volume={43}, + pages={100642}, + year={2024}, + doi={10.1016/j.wace.2024.100642}, + publisher={Elsevier} +} \ No newline at end of file From 4bd1c0b0ed8caedb090722dd7b3a862e11c95d85 Mon Sep 17 00:00:00 2001 From: Elizaveta Malinina Date: Wed, 29 Jan 2025 00:13:26 +0000 Subject: [PATCH 03/12] add mpl-style file for malinina24 --- .../matplotlib/malinina24.mplstyle | 28 +++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100644 esmvaltool/diag_scripts/shared/plot/styles_python/matplotlib/malinina24.mplstyle diff --git a/esmvaltool/diag_scripts/shared/plot/styles_python/matplotlib/malinina24.mplstyle b/esmvaltool/diag_scripts/shared/plot/styles_python/matplotlib/malinina24.mplstyle new file mode 100644 index 0000000000..f262b98295 --- /dev/null +++ b/esmvaltool/diag_scripts/shared/plot/styles_python/matplotlib/malinina24.mplstyle @@ -0,0 +1,28 @@ +############################################################################### +# Malinina&Gillett (2024) +# Style file for matplotlib plots +############################################################################### + +# Font +font.size : 10 + +# Fontsizes +axes.labelsize : medium +axes.titlesize : larger +xtick.labelsize : medium +ytick.labelsize : medium + +# Object sizes +lines.linewidth : 1 +lines.markersize : 6 + +# Axis +axes.spines.top : False +axes.spines.right : False +axes.linewidth: 0.5 +xtick.direction : out +ytick.direction : out +xtick.major.width : 0.5 +ytick.major.width : 0.5 +xtick.minor.width : 0.4 +ytick.minor.width : 0.4 \ No newline at end of file From fe356251f9aece3e7afd11ae1cbf0a4cfe894a4a Mon Sep 17 00:00:00 2001 From: Elizaveta Malinina Date: Wed, 29 Jan 2025 00:21:38 +0000 Subject: [PATCH 04/12] add variability diagnostics --- .../eccc_extremes/extreme_temperature_map.py | 173 +++++++++++++++ .../extremes_temperature_analytics.py | 176 +++++++++++++++ .../eccc_extremes/simple_timeseries.py | 201 +++++++++++++++++ .../eccc_extremes/variability_evaluation.py | 210 ++++++++++++++++++ 4 files changed, 760 insertions(+) create mode 100644 esmvaltool/diag_scripts/eccc_extremes/extreme_temperature_map.py create mode 100644 esmvaltool/diag_scripts/eccc_extremes/extremes_temperature_analytics.py create mode 100644 esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py create mode 100644 esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py diff --git a/esmvaltool/diag_scripts/eccc_extremes/extreme_temperature_map.py b/esmvaltool/diag_scripts/eccc_extremes/extreme_temperature_map.py new file mode 100644 index 0000000000..6993ff89e4 --- /dev/null +++ b/esmvaltool/diag_scripts/eccc_extremes/extreme_temperature_map.py @@ -0,0 +1,173 @@ +import iris +import iris.plot as iplt +import cartopy.crs as ccrs +import cartopy.feature as cf +import logging +import numpy as np +import matplotlib.pyplot as plt +from matplotlib.colors import LinearSegmentedColormap +import os + + +# import internal esmvaltool modules here +from esmvaltool.diag_scripts.shared import run_diagnostic, group_metadata +import esmvaltool.diag_scripts.shared.plot as eplot + + +# # This part sends debug statements to stdout +logger = logging.getLogger(os.path.basename(__file__)) +# logging.getLogger().addHandler(logging.StreamHandler(sys.stdout)) + +def create_map_plot(data_dic, mixns, ano_map_cb, cfg): + + b_min = np.floor(np.min([np.percentile(ano_map_cb.data, 1, method='closest_observation'), mixns['min']])) + b_max = np.ceil(np.max([np.percentile(ano_map_cb.data, 99, method='closest_observation'), mixns['max']])) + prov_borders = cf.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='50m', facecolor='none', edgecolor='k') + + prov_stat_dic = {'AB':{'lat': 54, 'lon': -114.5}, 'BC': {'lat': 54, 'lon': -123}, + 'MB': {'lat': 54, 'lon': -99}, 'NB': {'lat': 46.6, 'lon': -67.5}, + 'NL': {'lat': 53, 'lon': -62}, 'NS': {'lat': 44.25, 'lon': -65.7}, + 'NT': {'lat': 63, 'lon': -121}, 'NU': {'lat': 65, 'lon': -99}, + 'ON': {'lat': 50, 'lon': -86}, 'PE': {'lat': 46.4, 'lon': -63.2}, + 'QC': {'lat': 53, 'lon': -72}, 'SK':{'lat': 54, 'lon': -107}, + 'YT': {'lat': 63, 'lon': -136}, 'ID':{'lat': 44, 'lon': -116}, + 'ME': {'lat': 45, 'lon': -70.2}, 'MI': {'lat': 44.8, 'lon': -85.7}, + 'MN': {'lat': 46, 'lon': -95}, 'MT': {'lat': 46.5, 'lon': -110}, + 'ND': {'lat': 47, 'lon': -101}, 'NY': {'lat': 43, 'lon': -76}, + 'OR': {'lat': 43.5, 'lon': -121}, 'SD': {'lat': 44.5, 'lon': -101}, + 'WA': {'lat': 47, 'lon': -121}, 'WI': {'lat': 44, 'lon': -90}, + 'WY': {'lat': 43, 'lon': -108}} + + if cfg['cbar_positive']: + levels = np.arange(b_min, b_max+0.1, 0.1) + cmap = 'Reds' + else: + abs_max = np.abs([b_min, b_max]).max() + levels = np.arange(-1*abs_max, abs_max+0.1, 0.1) + cmap = 'RdBu_r' + + ratio_colors = [(0, (215/256, 179/256, 119/256)), + (0.4, (245/256, 240/256, 223/256)), + # (0.286, (245/256, 240/256, 223/256)), + (1.0, (66/256, 106/256, 90/256))] + ratio_cmap = LinearSegmentedColormap.from_list('ratio_cmap',ratio_colors) + + for dataset in data_dic.keys(): + fig_map, ax_map = plt.subplots(nrows=2, ncols=1, sharex=False, sharey=False, + subplot_kw={'projection': ccrs.LambertConformal(central_longitude=-91.87, + standard_parallels=(49,77))}) + ax_map = ax_map.flatten() + + fig_map.set_size_inches(7., 9.) + fig_map.set_dpi(200) + + model_data = list() + weights = list() + + for i in range(0,len(data_dic[dataset])): + model_data.append(data_dic[dataset][i].data) + weights.append(data_dic[dataset][i].attributes['ensemble_weight']* + data_dic[dataset][i].attributes['reverse_dtsts_n']) + weights = np.array(weights) + model_data = np.array(model_data) + # if dataset == 'Multi-Model-Mean': + # mean_arr = np.average(model_data, axis=0, weights=weights) + # else: + mean_arr = np.average(model_data, axis=0) + + map_c = ax_map[0].contourf(data_dic[dataset][i].coord('longitude').points, + data_dic[dataset][i].coord('latitude').points, + mean_arr, levels=levels, extend='both', cmap=cmap) + ax_map[0].set_title(dataset+' '+cfg['title_var_label']) + + iplt.contourf(ano_map_cb, axes=ax_map[1], levels=levels, extend='both', cmap=cmap) + ax_map[1].set_title('ERA5 '+ cfg['title_var_label']) + + [a.coastlines(linewidth=0.5) for a in ax_map] + [a.add_feature(cf.BORDERS) for a in ax_map] + [a.add_feature(prov_borders) for a in ax_map] + [a.set_extent([-63, -123,37,75]) for a in ax_map] + fig_map.subplots_adjust(left=0.02, bottom=0.1, right=0.99, top=0.92, wspace=0.02, hspace=0.1) + cax = fig_map.add_axes([0.1,0.07,0.8,0.015]) + fig_map.colorbar(map_c, cax=cax, orientation='horizontal', label=cfg['cbar_label'] +', '+ cfg['cbar_units']) + + fig_map.savefig(os.path.join(cfg['plot_dir'], 'map_anomalies_'+dataset+'.'+cfg['output_file_type'])) + + if cfg['ratio']: + fig_rat, ax_rat = plt.subplots(nrows=1, ncols=1, sharex=False, sharey=False, + subplot_kw={'projection': ccrs.LambertConformal(central_longitude=-91.87, + standard_parallels=(49,77))}) + + fig_rat.set_size_inches(7., 5.) + fig_rat.set_dpi(200) + + map_rat = iplt.contourf(1/(ano_map_cb/mean_arr), axes=ax_rat, levels=np.arange(0.2, 2.1, 0.1), extend='both', cmap=ratio_cmap) + + ax_rat.set_title('Ratio of '+cfg['title_var_label']+' ('+dataset+'/ERA5)') + ax_rat.coastlines(linewidth=0.5) + ax_rat.add_feature(cf.BORDERS, linewidth=0.5) + ax_rat.add_feature(prov_borders, linewidth=0.5) + for pr_st in prov_stat_dic.keys(): + ax_rat.text(prov_stat_dic[pr_st]['lon'],prov_stat_dic[pr_st]['lat'], pr_st, transform=ccrs.PlateCarree()) + + ax_rat.set_extent([-60.5, -121.5,37.5,73.5]) + + fig_rat.subplots_adjust(left=0.02, bottom=0.15, right=0.99, top=0.92, wspace=0.02, hspace=0.1) + cax = fig_rat.add_axes([0.1,0.13,0.8,0.012]) + fig_rat.colorbar(map_rat, cax=cax, orientation='horizontal', label=cfg['cbar_label'] +' ratio') + + fig_rat.savefig(os.path.join(cfg['plot_dir'], 'ratio_anomalies_'+dataset+'.'+cfg['output_file_type'])) + + + return + + +def main(cfg): + + input_data = cfg['input_data'] + + groups = group_metadata(input_data.values(), 'variable_group', sort=True) + obs_info = groups.pop('obs') + ano_map_cb = iris.load_cube(obs_info[0]['filename']) + + groups_l = list(groups.keys()) + + plotting_dic = {} + mixns = {} + + for group in groups_l: + mins = list(); maxs = list() + plotting_dic[group] = {}; mixns[group] = {} + group_data = groups[group] + datasets = group_metadata(group_data, 'dataset') + ens_cubelist = iris.cube.CubeList() + for dataset in datasets.keys(): + filepaths = list(group_metadata(datasets[dataset], 'filename').keys()) + n_real = len(filepaths) + mod_cubelist = iris.cube.CubeList() + for filepath in filepaths: + mod_cb = iris.load_cube(filepath) + mod_cb.attributes['ensemble_weight'] = 1 / n_real + mod_cb.attributes['reverse_dtsts_n'] = 1/ len(datasets) + ens_cubelist.append(mod_cb) + mod_cubelist.append(mod_cb) + mins.append(mod_cb.data.min()); maxs.append(mod_cb.data.max()) + plotting_dic[group][dataset] = mod_cubelist + mixns[group]['max'] = np.asarray(maxs).max() + mixns[group]['min'] = np.asarray(mins).min() + plotting_dic[group]['Multi-Model-Mean'] = ens_cubelist + + st_file = eplot.get_path_to_mpl_style(cfg.get('mpl_style')) + plt.style.use(st_file) + + create_map_plot(plotting_dic['txx_all'], mixns['txx_all'], ano_map_cb, cfg) + + logger.info('Success') + + +if __name__ == '__main__': + # always use run_diagnostic() to get the config (the preprocessor + # nested dictionary holding all the needed information) + with run_diagnostic() as config: + # list here the functions that need to run + main(config) \ No newline at end of file diff --git a/esmvaltool/diag_scripts/eccc_extremes/extremes_temperature_analytics.py b/esmvaltool/diag_scripts/eccc_extremes/extremes_temperature_analytics.py new file mode 100644 index 0000000000..803903737b --- /dev/null +++ b/esmvaltool/diag_scripts/eccc_extremes/extremes_temperature_analytics.py @@ -0,0 +1,176 @@ +import csv +import iris +import climextremes as cex +import pandas as pd +import logging +import numpy as np +import matplotlib.pyplot as plt +import os +from scipy.stats import genextreme as gev + +# import internal esmvaltool modules here +from esmvaltool.diag_scripts.shared import run_diagnostic, select_metadata, group_metadata +import esmvaltool.diag_scripts.shared.plot as eplot +from esmvaltool.diag_scripts.ocean import diagnostic_tools as diagtools +# # This part sends debug statements to stdout +logger = logging.getLogger(os.path.basename(__file__)) +# logging.getLogger().addHandler(logging.StreamHandler(sys.stdout)) + + +def calculate_stds(data_dic, obs_gev_data, cfg): + + era_var = obs_gev_data['ano_obs_cb'].data.std() + + col_mod = (25/255, 14/255, 79/255) + col_obs = 'indianred' + + fig_stds, ax_stds = plt.subplots(nrows=1, ncols=1) + + fig_stds.set_size_inches(8., 9.) + fig_stds.set_dpi(200) + + y_ticks = np.arange(0, len(data_dic.keys())) + y_labs = np.zeros(len(data_dic.keys()), dtype='1: + # clean this noncense before submitting! + ax_ts.fill_between(t_s, min_var_arr, max_var_arr, color=col_mod, lw=0, alpha=0.25) + ax_ts.plot([t_s-0.5, t_s+0.5], [0,0], c='grey') + ax_ts.set_xlim(t_s[0]-0.5, t_s[-1]+0.5) + ax_ts.set_ylim(b_min, b_max) + ax_ts.legend(loc=0, ncols=2, fancybox=False, frameon=False) + ax_ts.set_ylabel(cfg['var_label'] + ' anomaly, C') + if dataset == 'Multi-Model-Mean': + fig_ts.suptitle(cfg.get('litera')+' '+cfg['var_label']+' anomalies in '+cfg['region']+' relative to '+ str(cfg['reference_period'][0]) \ + + '-' + str(cfg['reference_period'][1])+ ' from '+str(len(data_dic.keys()) - 1) +' CMIP6 models' , fontsize = 'x-large') + else: + fig_ts.suptitle(cfg.get('litera')+' '+cfg['var_label']+' anomalies in '+cfg['region']+' relative to '+ str(cfg['reference_period'][0]) \ + + '-' + str(cfg['reference_period'][1])+ ' from '+dataset, fontsize = 'x-large') + + plt.tight_layout() + + fig_ts.savefig(os.path.join(cfg['plot_dir'], 'ts_'+cfg['region'].lower()+'_'+cfg['var_label'].lower()+'_'+dataset+'.'+cfg['output_file_type'])) + + return + + +def main(cfg): + + input_data = cfg['input_data'] + + groups = group_metadata(input_data.values(), 'variable_group', sort=True) + + ano_obs_cb, obs_gev_data = obtain_obs_info(groups, cfg) + + # change later, so far to identify the group name + group = [k for k in groups.keys() if 'txnx' in k][0] + + plotting_dic = {} + mixns = {} + + mins = list(); maxs = list() + group_data = groups[group] + datasets = group_metadata(group_data, 'dataset') + ens_var_cubelist = iris.cube.CubeList() + for dataset in datasets.keys(): + filepaths = list(group_metadata(datasets[dataset], 'filename').keys()) + n_real = len(filepaths) + mod_var_cubelist = iris.cube.CubeList() + for filepath in filepaths: + mod_cb = iris.load_cube(filepath) + # adding weights to the data cubes + mod_cb.attributes['ensemble_weight'] = 1 / n_real + mod_cb.attributes['reverse_dtsts_n'] = 1/ len(datasets) + ens_var_cubelist.append(mod_cb); mod_var_cubelist.append(mod_cb) + mins.append(mod_cb.collapsed('time', iris.analysis.MIN).data) + maxs.append(mod_cb.collapsed('time', iris.analysis.MAX).data) + plotting_dic[dataset] = {'var_data': mod_var_cubelist} + mixns['max'] = np.asarray(maxs).max() + mixns['min'] = np.asarray(mins).min() + plotting_dic['Multi-Model-Mean'] = {'var_data' : ens_var_cubelist} + + st_file = eplot.get_path_to_mpl_style(cfg.get('mpl_style')) + plt.style.use(st_file) + + create_timeseries(plotting_dic, mixns, obs_gev_data, cfg) + + calculate_stds(plotting_dic, obs_gev_data, cfg) + + logger.info('Success') + + +if __name__ == '__main__': + # always use run_diagnostic() to get the config (the preprocessor + # nested dictionary holding all the needed information) + with run_diagnostic() as config: + # list here the functions that need to run + main(config) diff --git a/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py b/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py new file mode 100644 index 0000000000..1e8fd953a6 --- /dev/null +++ b/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py @@ -0,0 +1,201 @@ +import iris +import logging +import numpy as np +import matplotlib.pyplot as plt +import os + +# import internal esmvaltool modules here +from esmvaltool.diag_scripts.shared import run_diagnostic, group_metadata, save_figure +import esmvaltool.diag_scripts.shared.plot as eplot + +# # This part sends debug statements to stdout +logger = logging.getLogger(os.path.basename(__file__)) + + +def obtain_reference(data_group: list): + ''' + This function cretes a dictionary with reference data. + Parameters: + ----------- + data_group: + list with metadata for reference group + + Returns: + -------- + reference_dic: + dictionary with the reference data. Dataset name is a keyword. + ref_max: + maximum value for all reference datasets + ref_min: + minimum value for all reference datasets + ''' + + reference_dic = {} + + ref_fnames = group_metadata(data_group, 'filename') + + mins = list() + maxs = list() + + for dataset_f in ref_fnames.keys(): + dataset_n = ref_fnames[dataset_f][0]['dataset'] + if len(group_metadata(data_group, 'dataset')[dataset_n])>1: + key = ref_fnames[dataset_f][0]['alias'] + else: + key = dataset_n + ref_cb = iris.load_cube(dataset_f) + reference_dic[key] = ref_cb + # make key dataset if there's only one filename of this dataset otherwise alias + mins.append(ref_cb.collapsed('time', iris.analysis.MIN).data) + maxs.append(ref_cb.collapsed('time', iris.analysis.MAX).data) + + ref_max = np.asarray(maxs).max() + ref_min = np.asarray(mins).min() + + return reference_dic, ref_max, ref_min + + +def create_provenance(caption: str): + '''Creates provenance dictionary''' + + provenance_dic = {'authors': ['malinina_elizaveta'], + 'caption': caption, + 'references': ['malinina24']} + + return provenance_dic + + +def plot_timeseries(data_dic: dict, reference_dic: dict, cfg: dict, + min_val: float, max_val: float): + ''' + This function plots timeseries and the max/min values + Parameters: + ----------- + data_dic: + dictionary with the data to plot + reference_dic: + dictionary with the refrence data + cfg: + config dictionary, comes from ESMValCore + min_val: + minimum value for the xlims + max_val: + maximum value for the xlims + ''' + + b_max = np.ceil(max_val*1.05) + b_min = np.floor(min_val*1.05) + + mpl_st_file = eplot.get_path_to_mpl_style(cfg.get('mpl_style')) + plt.style.use(mpl_st_file) + + for dataset in data_dic.keys(): + fig_ts, ax_ts = plt.subplots(nrows=1, ncols=1) + + model_var_data = list() + weights = list() + + for i in range(0,len(data_dic[dataset])): + var_cb = data_dic[dataset][i] + model_var_data.append(var_cb.data) + weights.append(var_cb.attributes['ensemble_weight']* + var_cb.attributes['reverse_dtsts_n']) + weights = np.array(weights) + model_var_data = np.array(model_var_data) + mean_var_arr = np.average(model_var_data, axis=0) + min_var_arr = np.min(model_var_data, axis=0) + max_var_arr = np.max(model_var_data, axis=0) + + years = [var_cb.coord('time').cell(i).point.year for i in range(var_cb.coord('time').shape[0])] + + color_st = eplot.get_dataset_style(dataset, cfg.get('color_style')) + ax_ts.plot(years, mean_var_arr, c=color_st['color'], zorder=3, label=dataset) + if len(data_dic[dataset])>1: + ax_ts.fill_between(years, min_var_arr, max_var_arr, color=color_st['color'], lw=0, alpha=0.25) + + for ref_data in reference_dic.keys(): + ref_cb = reference_dic.get(ref_data) + ref_color_st = eplot.get_dataset_style(ref_data, cfg.get('color_style')) + ref_years = [ref_cb.coord('time').cell(i).point.year for i in range(ref_cb.coord('time').shape[0])] + ax_ts.plot(ref_years,ref_cb.data, label=ref_data, c=ref_color_st['color'], zorder=3) + + ax_ts.axhline(color='grey', zorder=1) + ax_ts.set_ylim(b_min, b_max) + ax_ts.legend(loc=0, ncols=4, fancybox=False, frameon=False) + + variable = cfg['var_label'] if cfg.get('var_label') else var_cb.var_name + exp_variable = variable.replace('_', ' ') + units = cfg['units'] if cfg.get('units') else var_cb.units + region = cfg['region'] if cfg.get('units') else 'region' + + ax_ts.set_ylabel(f'{exp_variable}, {units}') + ax_ts.set_xlabel('year') + + default_caption = f'Timeseries of {exp_variable} in {region}' + + caption = cfg['figure_caption'] if cfg.get('figure_caption') else default_caption + + caption += f' ({dataset})' + + fig_ts.suptitle(caption) + + prov_dic = create_provenance(caption) + + plt.tight_layout() + + fig_path = os.path.join(cfg['plot_dir'], f'timeseries_{variable}_{region}_{dataset}') + + save_figure(fig_path, prov_dic, cfg, fig_ts, close=True) + + return + + +def main(cfg): + + input_data = cfg['input_data'] + + groups = group_metadata(input_data.values(), 'variable_group', sort=True) + + maxs, mins = [], [] + data_dic = {} + + reference_dic, ref_max, ref_min = obtain_reference(groups.pop('reference')) + + maxs.append(ref_max) + mins.append(ref_min) + + remaining_metadata = [] + for k in groups.keys(): + remaining_metadata.extend(groups[k]) + + datasets = group_metadata(remaining_metadata, 'dataset') + ens_var_cubelist = iris.cube.CubeList() + for dataset in datasets.keys(): + filepaths = list(group_metadata(datasets[dataset], 'filename').keys()) + n_real = len(filepaths) + mod_var_cubelist = iris.cube.CubeList() + for filepath in filepaths: + mod_cb = iris.load_cube(filepath) + # adding weights to the data cubes + mod_cb.attributes['ensemble_weight'] = 1 / n_real + mod_cb.attributes['reverse_dtsts_n'] = 1/ len(datasets) + ens_var_cubelist.append(mod_cb) + mod_var_cubelist.append(mod_cb) + mins.append(mod_cb.collapsed('time', iris.analysis.MIN).data) + maxs.append(mod_cb.collapsed('time', iris.analysis.MAX).data) + data_dic[dataset] = mod_var_cubelist + data_dic['Multi-Model'] = ens_var_cubelist + + plot_timeseries(data_dic, reference_dic, cfg, + min_val=np.asarray(mins).min(), + max_val=np.asarray(maxs).max()) + + logger.info('Success') + + +if __name__ == '__main__': + # always use run_diagnostic() to get the config (the preprocessor + # nested dictionary holding all the needed information) + with run_diagnostic() as config: + # list here the functions that need to run + main(config) \ No newline at end of file diff --git a/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py b/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py new file mode 100644 index 0000000000..b7803908b0 --- /dev/null +++ b/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py @@ -0,0 +1,210 @@ +import iris +import iris.cube +import logging +import numpy as np +import matplotlib.pyplot as plt +import os + +# import internal esmvaltool modules here +from esmvaltool.diag_scripts.shared import run_diagnostic, group_metadata, save_figure +import esmvaltool.diag_scripts.shared.plot as eplot +# # This part sends debug statements to stdout +logger = logging.getLogger(os.path.basename(__file__)) +# logging.getLogger().addHandler(logging.StreamHandler(sys.stdout)) + + +def obtain_reference(data_group: list): + ''' + This function cretes a dictionary with reference data. + Parameters: + ----------- + data_group: + list with metadata for reference group + + Returns: + -------- + reference_dic: + dictionary with the reference data. Dataset name is a keyword + ''' + + reference_dic = {} + + ref_fnames = group_metadata(data_group, 'filename') + + for dataset_f in ref_fnames.keys(): + dataset_n = ref_fnames[dataset_f][0]['dataset'] + if len(group_metadata(data_group, 'dataset')[dataset_n])>1: + key = ref_fnames[dataset_f][0]['alias'] + else: + key = dataset_n + ref_cb = iris.load_cube(dataset_f) + reference_dic[key] = ref_cb.data.std() + + return reference_dic + + +def create_provenance(caption: str): + '''Creates provenance dictionary''' + + provenance_dic = {'authors': ['malinina_elizaveta'], + 'caption': caption, + 'references': ['malinina24']} + + return provenance_dic + + +def calculate_ensemble_stdev(ens_cubelist: iris.cube.CubeList, + weight_flag: bool): + ''' + This function calculates standard deviation of multi model ensemble + Parameters: + ----------- + ens_cubelist: + CubeList with the ensemble timeseries + weight_flag: + flag for weighting or not the data + + Returns: + -------- + ens_std: float + standard deviation for the ensemble + ''' + + wghts , wdata = [], [] + for i in range(0, len(ens_cubelist)): + cb = ens_cubelist[i] + cb_wght = cb.attributes['ensemble_weight']*cb.attributes['reverse_dtsts_n'] + wghts.append(np.full(cb.shape, cb_wght)) + wdata.append(ens_cubelist[i].data) + + wghts = np.asarray(wghts).flatten() + wdata = np.asarray(wdata).flatten() + + if weight_flag: + mmm = np.average(wdata, weights=wghts) + ens_std = np.sqrt(np.average((wdata - mmm)**2, weights=wghts)) + else: + ens_std = wdata.std() + + return ens_std + + +def plot_stdevs(data_dic, reference_dic, cfg): + ''' + This function plots timeseries and the max/min values + Parameters: + ----------- + data_dic: + dictionary with the data to plot + reference_dic: + dictionary with the reference data + cfg: + config dictionary, comes from ESMValCore + ''' + + + mpl_st_file = eplot.get_path_to_mpl_style(cfg.get('mpl_style')) + plt.style.use(mpl_st_file) + + fig_stds, ax_stds = plt.subplots(nrows=1, ncols=1) + + fig_stds.set_size_inches(8., 9.) + fig_stds.set_dpi(200) + + y_ticks = np.arange(0, len(data_dic.keys())) + y_labs = np.zeros(len(data_dic.keys()), dtype=' Date: Fri, 31 Jan 2025 22:16:05 +0000 Subject: [PATCH 05/12] remove wrong diags --- .../eccc_extremes/extreme_temperature_map.py | 173 ----------------- .../extremes_temperature_analytics.py | 176 ------------------ 2 files changed, 349 deletions(-) delete mode 100644 esmvaltool/diag_scripts/eccc_extremes/extreme_temperature_map.py delete mode 100644 esmvaltool/diag_scripts/eccc_extremes/extremes_temperature_analytics.py diff --git a/esmvaltool/diag_scripts/eccc_extremes/extreme_temperature_map.py b/esmvaltool/diag_scripts/eccc_extremes/extreme_temperature_map.py deleted file mode 100644 index 6993ff89e4..0000000000 --- a/esmvaltool/diag_scripts/eccc_extremes/extreme_temperature_map.py +++ /dev/null @@ -1,173 +0,0 @@ -import iris -import iris.plot as iplt -import cartopy.crs as ccrs -import cartopy.feature as cf -import logging -import numpy as np -import matplotlib.pyplot as plt -from matplotlib.colors import LinearSegmentedColormap -import os - - -# import internal esmvaltool modules here -from esmvaltool.diag_scripts.shared import run_diagnostic, group_metadata -import esmvaltool.diag_scripts.shared.plot as eplot - - -# # This part sends debug statements to stdout -logger = logging.getLogger(os.path.basename(__file__)) -# logging.getLogger().addHandler(logging.StreamHandler(sys.stdout)) - -def create_map_plot(data_dic, mixns, ano_map_cb, cfg): - - b_min = np.floor(np.min([np.percentile(ano_map_cb.data, 1, method='closest_observation'), mixns['min']])) - b_max = np.ceil(np.max([np.percentile(ano_map_cb.data, 99, method='closest_observation'), mixns['max']])) - prov_borders = cf.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='50m', facecolor='none', edgecolor='k') - - prov_stat_dic = {'AB':{'lat': 54, 'lon': -114.5}, 'BC': {'lat': 54, 'lon': -123}, - 'MB': {'lat': 54, 'lon': -99}, 'NB': {'lat': 46.6, 'lon': -67.5}, - 'NL': {'lat': 53, 'lon': -62}, 'NS': {'lat': 44.25, 'lon': -65.7}, - 'NT': {'lat': 63, 'lon': -121}, 'NU': {'lat': 65, 'lon': -99}, - 'ON': {'lat': 50, 'lon': -86}, 'PE': {'lat': 46.4, 'lon': -63.2}, - 'QC': {'lat': 53, 'lon': -72}, 'SK':{'lat': 54, 'lon': -107}, - 'YT': {'lat': 63, 'lon': -136}, 'ID':{'lat': 44, 'lon': -116}, - 'ME': {'lat': 45, 'lon': -70.2}, 'MI': {'lat': 44.8, 'lon': -85.7}, - 'MN': {'lat': 46, 'lon': -95}, 'MT': {'lat': 46.5, 'lon': -110}, - 'ND': {'lat': 47, 'lon': -101}, 'NY': {'lat': 43, 'lon': -76}, - 'OR': {'lat': 43.5, 'lon': -121}, 'SD': {'lat': 44.5, 'lon': -101}, - 'WA': {'lat': 47, 'lon': -121}, 'WI': {'lat': 44, 'lon': -90}, - 'WY': {'lat': 43, 'lon': -108}} - - if cfg['cbar_positive']: - levels = np.arange(b_min, b_max+0.1, 0.1) - cmap = 'Reds' - else: - abs_max = np.abs([b_min, b_max]).max() - levels = np.arange(-1*abs_max, abs_max+0.1, 0.1) - cmap = 'RdBu_r' - - ratio_colors = [(0, (215/256, 179/256, 119/256)), - (0.4, (245/256, 240/256, 223/256)), - # (0.286, (245/256, 240/256, 223/256)), - (1.0, (66/256, 106/256, 90/256))] - ratio_cmap = LinearSegmentedColormap.from_list('ratio_cmap',ratio_colors) - - for dataset in data_dic.keys(): - fig_map, ax_map = plt.subplots(nrows=2, ncols=1, sharex=False, sharey=False, - subplot_kw={'projection': ccrs.LambertConformal(central_longitude=-91.87, - standard_parallels=(49,77))}) - ax_map = ax_map.flatten() - - fig_map.set_size_inches(7., 9.) - fig_map.set_dpi(200) - - model_data = list() - weights = list() - - for i in range(0,len(data_dic[dataset])): - model_data.append(data_dic[dataset][i].data) - weights.append(data_dic[dataset][i].attributes['ensemble_weight']* - data_dic[dataset][i].attributes['reverse_dtsts_n']) - weights = np.array(weights) - model_data = np.array(model_data) - # if dataset == 'Multi-Model-Mean': - # mean_arr = np.average(model_data, axis=0, weights=weights) - # else: - mean_arr = np.average(model_data, axis=0) - - map_c = ax_map[0].contourf(data_dic[dataset][i].coord('longitude').points, - data_dic[dataset][i].coord('latitude').points, - mean_arr, levels=levels, extend='both', cmap=cmap) - ax_map[0].set_title(dataset+' '+cfg['title_var_label']) - - iplt.contourf(ano_map_cb, axes=ax_map[1], levels=levels, extend='both', cmap=cmap) - ax_map[1].set_title('ERA5 '+ cfg['title_var_label']) - - [a.coastlines(linewidth=0.5) for a in ax_map] - [a.add_feature(cf.BORDERS) for a in ax_map] - [a.add_feature(prov_borders) for a in ax_map] - [a.set_extent([-63, -123,37,75]) for a in ax_map] - fig_map.subplots_adjust(left=0.02, bottom=0.1, right=0.99, top=0.92, wspace=0.02, hspace=0.1) - cax = fig_map.add_axes([0.1,0.07,0.8,0.015]) - fig_map.colorbar(map_c, cax=cax, orientation='horizontal', label=cfg['cbar_label'] +', '+ cfg['cbar_units']) - - fig_map.savefig(os.path.join(cfg['plot_dir'], 'map_anomalies_'+dataset+'.'+cfg['output_file_type'])) - - if cfg['ratio']: - fig_rat, ax_rat = plt.subplots(nrows=1, ncols=1, sharex=False, sharey=False, - subplot_kw={'projection': ccrs.LambertConformal(central_longitude=-91.87, - standard_parallels=(49,77))}) - - fig_rat.set_size_inches(7., 5.) - fig_rat.set_dpi(200) - - map_rat = iplt.contourf(1/(ano_map_cb/mean_arr), axes=ax_rat, levels=np.arange(0.2, 2.1, 0.1), extend='both', cmap=ratio_cmap) - - ax_rat.set_title('Ratio of '+cfg['title_var_label']+' ('+dataset+'/ERA5)') - ax_rat.coastlines(linewidth=0.5) - ax_rat.add_feature(cf.BORDERS, linewidth=0.5) - ax_rat.add_feature(prov_borders, linewidth=0.5) - for pr_st in prov_stat_dic.keys(): - ax_rat.text(prov_stat_dic[pr_st]['lon'],prov_stat_dic[pr_st]['lat'], pr_st, transform=ccrs.PlateCarree()) - - ax_rat.set_extent([-60.5, -121.5,37.5,73.5]) - - fig_rat.subplots_adjust(left=0.02, bottom=0.15, right=0.99, top=0.92, wspace=0.02, hspace=0.1) - cax = fig_rat.add_axes([0.1,0.13,0.8,0.012]) - fig_rat.colorbar(map_rat, cax=cax, orientation='horizontal', label=cfg['cbar_label'] +' ratio') - - fig_rat.savefig(os.path.join(cfg['plot_dir'], 'ratio_anomalies_'+dataset+'.'+cfg['output_file_type'])) - - - return - - -def main(cfg): - - input_data = cfg['input_data'] - - groups = group_metadata(input_data.values(), 'variable_group', sort=True) - obs_info = groups.pop('obs') - ano_map_cb = iris.load_cube(obs_info[0]['filename']) - - groups_l = list(groups.keys()) - - plotting_dic = {} - mixns = {} - - for group in groups_l: - mins = list(); maxs = list() - plotting_dic[group] = {}; mixns[group] = {} - group_data = groups[group] - datasets = group_metadata(group_data, 'dataset') - ens_cubelist = iris.cube.CubeList() - for dataset in datasets.keys(): - filepaths = list(group_metadata(datasets[dataset], 'filename').keys()) - n_real = len(filepaths) - mod_cubelist = iris.cube.CubeList() - for filepath in filepaths: - mod_cb = iris.load_cube(filepath) - mod_cb.attributes['ensemble_weight'] = 1 / n_real - mod_cb.attributes['reverse_dtsts_n'] = 1/ len(datasets) - ens_cubelist.append(mod_cb) - mod_cubelist.append(mod_cb) - mins.append(mod_cb.data.min()); maxs.append(mod_cb.data.max()) - plotting_dic[group][dataset] = mod_cubelist - mixns[group]['max'] = np.asarray(maxs).max() - mixns[group]['min'] = np.asarray(mins).min() - plotting_dic[group]['Multi-Model-Mean'] = ens_cubelist - - st_file = eplot.get_path_to_mpl_style(cfg.get('mpl_style')) - plt.style.use(st_file) - - create_map_plot(plotting_dic['txx_all'], mixns['txx_all'], ano_map_cb, cfg) - - logger.info('Success') - - -if __name__ == '__main__': - # always use run_diagnostic() to get the config (the preprocessor - # nested dictionary holding all the needed information) - with run_diagnostic() as config: - # list here the functions that need to run - main(config) \ No newline at end of file diff --git a/esmvaltool/diag_scripts/eccc_extremes/extremes_temperature_analytics.py b/esmvaltool/diag_scripts/eccc_extremes/extremes_temperature_analytics.py deleted file mode 100644 index 803903737b..0000000000 --- a/esmvaltool/diag_scripts/eccc_extremes/extremes_temperature_analytics.py +++ /dev/null @@ -1,176 +0,0 @@ -import csv -import iris -import climextremes as cex -import pandas as pd -import logging -import numpy as np -import matplotlib.pyplot as plt -import os -from scipy.stats import genextreme as gev - -# import internal esmvaltool modules here -from esmvaltool.diag_scripts.shared import run_diagnostic, select_metadata, group_metadata -import esmvaltool.diag_scripts.shared.plot as eplot -from esmvaltool.diag_scripts.ocean import diagnostic_tools as diagtools -# # This part sends debug statements to stdout -logger = logging.getLogger(os.path.basename(__file__)) -# logging.getLogger().addHandler(logging.StreamHandler(sys.stdout)) - - -def calculate_stds(data_dic, obs_gev_data, cfg): - - era_var = obs_gev_data['ano_obs_cb'].data.std() - - col_mod = (25/255, 14/255, 79/255) - col_obs = 'indianred' - - fig_stds, ax_stds = plt.subplots(nrows=1, ncols=1) - - fig_stds.set_size_inches(8., 9.) - fig_stds.set_dpi(200) - - y_ticks = np.arange(0, len(data_dic.keys())) - y_labs = np.zeros(len(data_dic.keys()), dtype='1: - # clean this noncense before submitting! - ax_ts.fill_between(t_s, min_var_arr, max_var_arr, color=col_mod, lw=0, alpha=0.25) - ax_ts.plot([t_s-0.5, t_s+0.5], [0,0], c='grey') - ax_ts.set_xlim(t_s[0]-0.5, t_s[-1]+0.5) - ax_ts.set_ylim(b_min, b_max) - ax_ts.legend(loc=0, ncols=2, fancybox=False, frameon=False) - ax_ts.set_ylabel(cfg['var_label'] + ' anomaly, C') - if dataset == 'Multi-Model-Mean': - fig_ts.suptitle(cfg.get('litera')+' '+cfg['var_label']+' anomalies in '+cfg['region']+' relative to '+ str(cfg['reference_period'][0]) \ - + '-' + str(cfg['reference_period'][1])+ ' from '+str(len(data_dic.keys()) - 1) +' CMIP6 models' , fontsize = 'x-large') - else: - fig_ts.suptitle(cfg.get('litera')+' '+cfg['var_label']+' anomalies in '+cfg['region']+' relative to '+ str(cfg['reference_period'][0]) \ - + '-' + str(cfg['reference_period'][1])+ ' from '+dataset, fontsize = 'x-large') - - plt.tight_layout() - - fig_ts.savefig(os.path.join(cfg['plot_dir'], 'ts_'+cfg['region'].lower()+'_'+cfg['var_label'].lower()+'_'+dataset+'.'+cfg['output_file_type'])) - - return - - -def main(cfg): - - input_data = cfg['input_data'] - - groups = group_metadata(input_data.values(), 'variable_group', sort=True) - - ano_obs_cb, obs_gev_data = obtain_obs_info(groups, cfg) - - # change later, so far to identify the group name - group = [k for k in groups.keys() if 'txnx' in k][0] - - plotting_dic = {} - mixns = {} - - mins = list(); maxs = list() - group_data = groups[group] - datasets = group_metadata(group_data, 'dataset') - ens_var_cubelist = iris.cube.CubeList() - for dataset in datasets.keys(): - filepaths = list(group_metadata(datasets[dataset], 'filename').keys()) - n_real = len(filepaths) - mod_var_cubelist = iris.cube.CubeList() - for filepath in filepaths: - mod_cb = iris.load_cube(filepath) - # adding weights to the data cubes - mod_cb.attributes['ensemble_weight'] = 1 / n_real - mod_cb.attributes['reverse_dtsts_n'] = 1/ len(datasets) - ens_var_cubelist.append(mod_cb); mod_var_cubelist.append(mod_cb) - mins.append(mod_cb.collapsed('time', iris.analysis.MIN).data) - maxs.append(mod_cb.collapsed('time', iris.analysis.MAX).data) - plotting_dic[dataset] = {'var_data': mod_var_cubelist} - mixns['max'] = np.asarray(maxs).max() - mixns['min'] = np.asarray(mins).min() - plotting_dic['Multi-Model-Mean'] = {'var_data' : ens_var_cubelist} - - st_file = eplot.get_path_to_mpl_style(cfg.get('mpl_style')) - plt.style.use(st_file) - - create_timeseries(plotting_dic, mixns, obs_gev_data, cfg) - - calculate_stds(plotting_dic, obs_gev_data, cfg) - - logger.info('Success') - - -if __name__ == '__main__': - # always use run_diagnostic() to get the config (the preprocessor - # nested dictionary holding all the needed information) - with run_diagnostic() as config: - # list here the functions that need to run - main(config) From 83cb5e720e34f6719949bfba16fe16f8ef60d4a8 Mon Sep 17 00:00:00 2001 From: Elizaveta Malinina Date: Fri, 31 Jan 2025 22:55:12 +0000 Subject: [PATCH 06/12] styling formatting --- .../eccc_extremes/simple_timeseries.py | 118 +++++++++++------ .../eccc_extremes/variability_evaluation.py | 123 +++++++++++------- 2 files changed, 149 insertions(+), 92 deletions(-) diff --git a/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py b/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py index 1e8fd953a6..389e12fccd 100644 --- a/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py +++ b/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py @@ -1,25 +1,30 @@ -import iris import logging -import numpy as np -import matplotlib.pyplot as plt import os -# import internal esmvaltool modules here -from esmvaltool.diag_scripts.shared import run_diagnostic, group_metadata, save_figure +import iris +import matplotlib.pyplot as plt +import numpy as np + import esmvaltool.diag_scripts.shared.plot as eplot +# import internal esmvaltool modules here +from esmvaltool.diag_scripts.shared import ( + group_metadata, + run_diagnostic, + save_figure, +) # # This part sends debug statements to stdout logger = logging.getLogger(os.path.basename(__file__)) def obtain_reference(data_group: list): - ''' - This function cretes a dictionary with reference data. + """This function cretes a dictionary with reference data. + Parameters: ----------- - data_group: + data_group: list with metadata for reference group - + Returns: -------- reference_dic: @@ -28,7 +33,7 @@ def obtain_reference(data_group: list): maximum value for all reference datasets ref_min: minimum value for all reference datasets - ''' + """ reference_dic = {} @@ -39,9 +44,9 @@ def obtain_reference(data_group: list): for dataset_f in ref_fnames.keys(): dataset_n = ref_fnames[dataset_f][0]['dataset'] - if len(group_metadata(data_group, 'dataset')[dataset_n])>1: + if len(group_metadata(data_group, 'dataset')[dataset_n]) > 1: key = ref_fnames[dataset_f][0]['alias'] - else: + else: key = dataset_n ref_cb = iris.load_cube(dataset_f) reference_dic[key] = ref_cb @@ -56,19 +61,21 @@ def obtain_reference(data_group: list): def create_provenance(caption: str): - '''Creates provenance dictionary''' + """Creates provenance dictionary.""" - provenance_dic = {'authors': ['malinina_elizaveta'], - 'caption': caption, - 'references': ['malinina24']} + provenance_dic = { + 'authors': ['malinina_elizaveta'], + 'caption': caption, + 'references': ['malinina24'] + } return provenance_dic -def plot_timeseries(data_dic: dict, reference_dic: dict, cfg: dict, - min_val: float, max_val: float): - ''' - This function plots timeseries and the max/min values +def plot_timeseries(data_dic: dict, reference_dic: dict, cfg: dict, + min_val: float, max_val: float): + """This function plots timeseries and the max/min values. + Parameters: ----------- data_dic: @@ -79,51 +86,72 @@ def plot_timeseries(data_dic: dict, reference_dic: dict, cfg: dict, config dictionary, comes from ESMValCore min_val: minimum value for the xlims - max_val: + max_val: maximum value for the xlims - ''' + """ - b_max = np.ceil(max_val*1.05) - b_min = np.floor(min_val*1.05) + b_max = np.ceil(max_val * 1.05) + b_min = np.floor(min_val * 1.05) mpl_st_file = eplot.get_path_to_mpl_style(cfg.get('mpl_style')) plt.style.use(mpl_st_file) - for dataset in data_dic.keys(): + for dataset in data_dic.keys(): fig_ts, ax_ts = plt.subplots(nrows=1, ncols=1) model_var_data = list() weights = list() - for i in range(0,len(data_dic[dataset])): + for i in range(0, len(data_dic[dataset])): var_cb = data_dic[dataset][i] model_var_data.append(var_cb.data) - weights.append(var_cb.attributes['ensemble_weight']* - var_cb.attributes['reverse_dtsts_n']) + weights.append(var_cb.attributes['ensemble_weight'] * + var_cb.attributes['reverse_dtsts_n']) weights = np.array(weights) model_var_data = np.array(model_var_data) mean_var_arr = np.average(model_var_data, axis=0) min_var_arr = np.min(model_var_data, axis=0) max_var_arr = np.max(model_var_data, axis=0) - years = [var_cb.coord('time').cell(i).point.year for i in range(var_cb.coord('time').shape[0])] + years = [ + var_cb.coord('time').cell(i).point.year + for i in range(var_cb.coord('time').shape[0]) + ] color_st = eplot.get_dataset_style(dataset, cfg.get('color_style')) - ax_ts.plot(years, mean_var_arr, c=color_st['color'], zorder=3, label=dataset) - if len(data_dic[dataset])>1: - ax_ts.fill_between(years, min_var_arr, max_var_arr, color=color_st['color'], lw=0, alpha=0.25) + ax_ts.plot(years, + mean_var_arr, + c=color_st['color'], + zorder=3, + label=dataset) + if len(data_dic[dataset]) > 1: + ax_ts.fill_between(years, + min_var_arr, + max_var_arr, + color=color_st['color'], + lw=0, + alpha=0.25) for ref_data in reference_dic.keys(): ref_cb = reference_dic.get(ref_data) - ref_color_st = eplot.get_dataset_style(ref_data, cfg.get('color_style')) - ref_years = [ref_cb.coord('time').cell(i).point.year for i in range(ref_cb.coord('time').shape[0])] - ax_ts.plot(ref_years,ref_cb.data, label=ref_data, c=ref_color_st['color'], zorder=3) + ref_color_st = eplot.get_dataset_style(ref_data, + cfg.get('color_style')) + ref_years = [ + ref_cb.coord('time').cell(i).point.year + for i in range(ref_cb.coord('time').shape[0]) + ] + ax_ts.plot(ref_years, + ref_cb.data, + label=ref_data, + c=ref_color_st['color'], + zorder=3) ax_ts.axhline(color='grey', zorder=1) ax_ts.set_ylim(b_min, b_max) ax_ts.legend(loc=0, ncols=4, fancybox=False, frameon=False) - variable = cfg['var_label'] if cfg.get('var_label') else var_cb.var_name + variable = cfg['var_label'] if cfg.get( + 'var_label') else var_cb.var_name exp_variable = variable.replace('_', ' ') units = cfg['units'] if cfg.get('units') else var_cb.units region = cfg['region'] if cfg.get('units') else 'region' @@ -133,7 +161,8 @@ def plot_timeseries(data_dic: dict, reference_dic: dict, cfg: dict, default_caption = f'Timeseries of {exp_variable} in {region}' - caption = cfg['figure_caption'] if cfg.get('figure_caption') else default_caption + caption = cfg['figure_caption'] if cfg.get( + 'figure_caption') else default_caption caption += f' ({dataset})' @@ -143,7 +172,8 @@ def plot_timeseries(data_dic: dict, reference_dic: dict, cfg: dict, plt.tight_layout() - fig_path = os.path.join(cfg['plot_dir'], f'timeseries_{variable}_{region}_{dataset}') + fig_path = os.path.join(cfg['plot_dir'], + f'timeseries_{variable}_{region}_{dataset}') save_figure(fig_path, prov_dic, cfg, fig_ts, close=True) @@ -176,9 +206,9 @@ def main(cfg): mod_var_cubelist = iris.cube.CubeList() for filepath in filepaths: mod_cb = iris.load_cube(filepath) - # adding weights to the data cubes + # adding weights to the data cubes mod_cb.attributes['ensemble_weight'] = 1 / n_real - mod_cb.attributes['reverse_dtsts_n'] = 1/ len(datasets) + mod_cb.attributes['reverse_dtsts_n'] = 1 / len(datasets) ens_var_cubelist.append(mod_cb) mod_var_cubelist.append(mod_cb) mins.append(mod_cb.collapsed('time', iris.analysis.MIN).data) @@ -186,8 +216,10 @@ def main(cfg): data_dic[dataset] = mod_var_cubelist data_dic['Multi-Model'] = ens_var_cubelist - plot_timeseries(data_dic, reference_dic, cfg, - min_val=np.asarray(mins).min(), + plot_timeseries(data_dic, + reference_dic, + cfg, + min_val=np.asarray(mins).min(), max_val=np.asarray(maxs).max()) logger.info('Success') @@ -198,4 +230,4 @@ def main(cfg): # nested dictionary holding all the needed information) with run_diagnostic() as config: # list here the functions that need to run - main(config) \ No newline at end of file + main(config) diff --git a/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py b/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py index b7803908b0..def0914ba0 100644 --- a/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py +++ b/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py @@ -1,31 +1,37 @@ +import logging +import os + import iris import iris.cube -import logging -import numpy as np import matplotlib.pyplot as plt -import os +import numpy as np -# import internal esmvaltool modules here -from esmvaltool.diag_scripts.shared import run_diagnostic, group_metadata, save_figure import esmvaltool.diag_scripts.shared.plot as eplot +# import internal esmvaltool modules here +from esmvaltool.diag_scripts.shared import ( + group_metadata, + run_diagnostic, + save_figure, +) + # # This part sends debug statements to stdout logger = logging.getLogger(os.path.basename(__file__)) # logging.getLogger().addHandler(logging.StreamHandler(sys.stdout)) def obtain_reference(data_group: list): - ''' - This function cretes a dictionary with reference data. + """This function cretes a dictionary with reference data. + Parameters: ----------- - data_group: + data_group: list with metadata for reference group - + Returns: -------- reference_dic: dictionary with the reference data. Dataset name is a keyword - ''' + """ reference_dic = {} @@ -33,9 +39,9 @@ def obtain_reference(data_group: list): for dataset_f in ref_fnames.keys(): dataset_n = ref_fnames[dataset_f][0]['dataset'] - if len(group_metadata(data_group, 'dataset')[dataset_n])>1: + if len(group_metadata(data_group, 'dataset')[dataset_n]) > 1: key = ref_fnames[dataset_f][0]['alias'] - else: + else: key = dataset_n ref_cb = iris.load_cube(dataset_f) reference_dic[key] = ref_cb.data.std() @@ -44,36 +50,39 @@ def obtain_reference(data_group: list): def create_provenance(caption: str): - '''Creates provenance dictionary''' + """Creates provenance dictionary.""" - provenance_dic = {'authors': ['malinina_elizaveta'], - 'caption': caption, - 'references': ['malinina24']} + provenance_dic = { + 'authors': ['malinina_elizaveta'], + 'caption': caption, + 'references': ['malinina24'] + } return provenance_dic -def calculate_ensemble_stdev(ens_cubelist: iris.cube.CubeList, - weight_flag: bool): - ''' - This function calculates standard deviation of multi model ensemble +def calculate_ensemble_stdev(ens_cubelist: iris.cube.CubeList, + weight_flag: bool): + """This function calculates standard deviation of multi model ensemble. + Parameters: ----------- - ens_cubelist: + ens_cubelist: CubeList with the ensemble timeseries weight_flag: flag for weighting or not the data - + Returns: -------- ens_std: float standard deviation for the ensemble - ''' + """ - wghts , wdata = [], [] + wghts, wdata = [], [] for i in range(0, len(ens_cubelist)): cb = ens_cubelist[i] - cb_wght = cb.attributes['ensemble_weight']*cb.attributes['reverse_dtsts_n'] + cb_wght = cb.attributes['ensemble_weight'] * cb.attributes[ + 'reverse_dtsts_n'] wghts.append(np.full(cb.shape, cb_wght)) wdata.append(ens_cubelist[i].data) @@ -83,15 +92,15 @@ def calculate_ensemble_stdev(ens_cubelist: iris.cube.CubeList, if weight_flag: mmm = np.average(wdata, weights=wghts) ens_std = np.sqrt(np.average((wdata - mmm)**2, weights=wghts)) - else: + else: ens_std = wdata.std() return ens_std def plot_stdevs(data_dic, reference_dic, cfg): - ''' - This function plots timeseries and the max/min values + """This function plots timeseries and the max/min values. + Parameters: ----------- data_dic: @@ -100,8 +109,7 @@ def plot_stdevs(data_dic, reference_dic, cfg): dictionary with the reference data cfg: config dictionary, comes from ESMValCore - ''' - + """ mpl_st_file = eplot.get_path_to_mpl_style(cfg.get('mpl_style')) plt.style.use(mpl_st_file) @@ -118,23 +126,34 @@ def plot_stdevs(data_dic, reference_dic, cfg): color_st = eplot.get_dataset_style(model, cfg.get('color_style')) if model != 'Multi-Model': for i in range(0, len(data_dic[model])): - single_dot = ax_stds.scatter(data_dic[model][i], nm+1, marker='o', - c=color_st['color'], zorder=2, - label='individual member') - y_labs[nm+1] = model + single_dot = ax_stds.scatter(data_dic[model][i], + nm + 1, + marker='o', + c=color_st['color'], + zorder=2, + label='individual member') + y_labs[nm + 1] = model else: - square = ax_stds.scatter(data_dic[model], 0, c=color_st['color'], s=70, - marker='s', zorder=3, label='full ensemble') + square = ax_stds.scatter(data_dic[model], + 0, + c=color_st['color'], + s=70, + marker='s', + zorder=3, + label='full ensemble') y_labs[0] = 'ALL' legend_handles = [square, single_dot] for ref_dataset in reference_dic.keys(): - ref_color_st = eplot.get_dataset_style(ref_dataset, cfg.get('color_style')) - ref_line = ax_stds.axvline(reference_dic[ref_dataset], c=ref_color_st['color'], zorder=2, - label=ref_dataset) + ref_color_st = eplot.get_dataset_style(ref_dataset, + cfg.get('color_style')) + ref_line = ax_stds.axvline(reference_dic[ref_dataset], + c=ref_color_st['color'], + zorder=2, + label=ref_dataset) legend_handles.append(ref_line) - ax_stds.set_ylim(len(data_dic.keys()) -0.8, -0.2) + ax_stds.set_ylim(len(data_dic.keys()) - 0.8, -0.2) ax_stds.set_yticks(y_ticks, labels=y_labs) ax_stds.grid(which='both', c='silver', zorder=1) @@ -148,19 +167,25 @@ def plot_stdevs(data_dic, reference_dic, cfg): default_caption = f'{variable} varibility in {region}' - caption = cfg['figure_caption'] if cfg.get('figure_caption') else default_caption + caption = cfg['figure_caption'] if cfg.get( + 'figure_caption') else default_caption fig_stds.suptitle(caption) prov_dic = create_provenance(caption) - plt.legend(handles=legend_handles, bbox_to_anchor =(0.5,-0.1), - loc='lower center', ncols=4, fancybox=False, frameon=False) + plt.legend(handles=legend_handles, + bbox_to_anchor=(0.5, -0.1), + loc='lower center', + ncols=4, + fancybox=False, + frameon=False) plt.tight_layout() - fig_path = os.path.join(cfg['plot_dir'], f'variability_{variable}_{region}') + fig_path = os.path.join(cfg['plot_dir'], + f'variability_{variable}_{region}') - save_figure(fig_path, prov_dic, cfg, fig_stds, close=True) + save_figure(fig_path, prov_dic, cfg, fig_stds, close=True) return @@ -187,13 +212,13 @@ def main(cfg): mod_var_list = [] for filepath in filepaths: mod_cb = iris.load_cube(filepath) - # adding weights to the data cubes + # adding weights to the data cubes mod_cb.attributes['ensemble_weight'] = 1 / n_real - mod_cb.attributes['reverse_dtsts_n'] = 1/ len(datasets) + mod_cb.attributes['reverse_dtsts_n'] = 1 / len(datasets) ens_var_cubelist.append(mod_cb) mod_var_list.append(mod_cb.data.std()) data_dic[dataset] = mod_var_list - multi_model_std = calculate_ensemble_stdev(ens_var_cubelist, + multi_model_std = calculate_ensemble_stdev(ens_var_cubelist, cfg.get('model_weighting')) data_dic['Multi-Model'] = multi_model_std @@ -207,4 +232,4 @@ def main(cfg): # nested dictionary holding all the needed information) with run_diagnostic() as config: # list here the functions that need to run - main(config) \ No newline at end of file + main(config) From 801b48238c59849388c54d7641eb21b8dd471389 Mon Sep 17 00:00:00 2001 From: Elizaveta Malinina Date: Fri, 31 Jan 2025 23:23:13 +0000 Subject: [PATCH 07/12] shorten long lines, change docstrings --- .../diag_scripts/eccc_extremes/simple_timeseries.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py b/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py index 389e12fccd..2fd614bb2d 100644 --- a/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py +++ b/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py @@ -49,8 +49,9 @@ def obtain_reference(data_group: list): else: key = dataset_n ref_cb = iris.load_cube(dataset_f) + # make key dataset if there's only one filename + # of this dataset otherwise alias reference_dic[key] = ref_cb - # make key dataset if there's only one filename of this dataset otherwise alias mins.append(ref_cb.collapsed('time', iris.analysis.MIN).data) maxs.append(ref_cb.collapsed('time', iris.analysis.MAX).data) @@ -62,7 +63,6 @@ def obtain_reference(data_group: list): def create_provenance(caption: str): """Creates provenance dictionary.""" - provenance_dic = { 'authors': ['malinina_elizaveta'], 'caption': caption, @@ -74,7 +74,8 @@ def create_provenance(caption: str): def plot_timeseries(data_dic: dict, reference_dic: dict, cfg: dict, min_val: float, max_val: float): - """This function plots timeseries and the max/min values. + """ + This function plots timeseries and the max/min values. Parameters: ----------- @@ -181,7 +182,7 @@ def plot_timeseries(data_dic: dict, reference_dic: dict, cfg: dict, def main(cfg): - + """Initiates diagnostic""" input_data = cfg['input_data'] groups = group_metadata(input_data.values(), 'variable_group', sort=True) From 49ed9d04f0287e28ccfc614e0a105837313d7f77 Mon Sep 17 00:00:00 2001 From: Elizaveta Malinina Date: Sat, 1 Feb 2025 00:15:23 +0000 Subject: [PATCH 08/12] add recipe model evaluation --- .../recipe_extremes_evaluation_SK_BC.yml | 191 ++++++++++++++++++ 1 file changed, 191 insertions(+) create mode 100644 esmvaltool/recipes/eccc_extremes/recipe_extremes_evaluation_SK_BC.yml diff --git a/esmvaltool/recipes/eccc_extremes/recipe_extremes_evaluation_SK_BC.yml b/esmvaltool/recipes/eccc_extremes/recipe_extremes_evaluation_SK_BC.yml new file mode 100644 index 0000000000..735679b6c0 --- /dev/null +++ b/esmvaltool/recipes/eccc_extremes/recipe_extremes_evaluation_SK_BC.yml @@ -0,0 +1,191 @@ +# ESMValTool +# recipe_extremes_evaluation_SK_BC.yml +--- +documentation: + + title: Extremes variability in the CMIP models + + description: This is an analysis to evaluate warm extremes (TXx) variability + in CMIP6 models used in event attribution. This recipe reproduces + Figs. 6 and 8 from Malinina and Gillett (2024). + + authors: + - malinina_elizaveta + +preprocessors: + + preproc_area_average_var: &prep_aa_var + area_statistics: + operator: mean + annual_statistics: + operator: max + anomalies: + period: full + reference: + start_year: 1940 + start_month: 1 + start_day: 1 + end_year: 1959 + end_month: 12 + end_day: 31 + + preproc_model_bc: + custom_order: True + regrid: + target_grid: 0.5x0.5 + lon_offset: 0.25 + lat_offset: 0.25 + scheme: linear + extract_shape: + shapefile: british_columbia.shp + method: contains + crop: True + *prep_aa_var + + preproc_era_bc: + custom_order: True + extract_shape: + shapefile: british_columbia.shp + method: contains + crop: True + *prep_aa_var + + preproc_model_sk: + custom_order: True + regrid: + target_grid: 0.5x0.5 + lon_offset: 0.25 + lat_offset: 0.25 + scheme: linear + extract_shape: + shapefile: saskatchewan.shp + method: contains + crop: True + *prep_aa_var + + preproc_era_sk: + custom_order: True + extract_shape: + shapefile: saskatchewan.shp + method: contains + crop: True + *prep_aa_var + + +datasets_cmip: &datasets_cmip + - {dataset: ACCESS-ESM1-5, ensemble: "r(1:18)i1p1f1", grid: gn} + - {dataset: ACCESS-ESM1-5, ensemble: "r(21:24)i1p1f1", grid: gn} + - {dataset: ACCESS-ESM1-5, ensemble: r26i1p1f1, grid: gn} + - {dataset: ACCESS-ESM1-5, ensemble: "r(28:33)i1p1f1", grid: gn} + - {dataset: ACCESS-ESM1-5, ensemble: "r(35:40)i1p1f1", grid: gn} + - {dataset: ACCESS-CM2, ensemble: "r(1:5)i1p1f1", grid: gn} + - {dataset: AWI-CM-1-1-MR, ensemble: r1i1p1f1, grid: gn} + - {dataset: BCC-CSM2-MR, ensemble: r1i1p1f1, grid: gn} + - {dataset: CanESM5, ensemble: "r(1:25)i1p1f1", grid: gn} + - {dataset: CanESM5, ensemble: "r(1:25)i1p2f1", grid: gn} + - {dataset: CMCC-ESM2, ensemble: r1i1p1f1, grid: gn} + - {dataset: CNRM-CM6-1, ensemble: r1i1p1f2, grid: gr} + - {dataset: CNRM-ESM2-1, ensemble: r1i1p1f2, grid: gr} + - { dataset: EC-Earth3, ensemble: r2i1p1f1, grid: gr } + - { dataset: EC-Earth3-CC, ensemble: r1i1p1f1, grid: gr } + - { dataset: EC-Earth3-Veg, ensemble: r1i1p1f1, grid: gr } + - { dataset: EC-Earth3-Veg, ensemble: r6i1p1f1, grid: gr } + - { dataset: GFDL-ESM4, institute: NOAA-GFDL, ensemble: r1i1p1f1, grid: gr1 } + - { dataset: INM-CM5-0, ensemble: r1i1p1f1, grid: gr1 } + - { dataset: INM-CM4-8, ensemble: r1i1p1f1, grid: gr1 } + - { dataset: IPSL-CM6A-LR, ensemble: "r(1:6)i1p1f1", grid: gr } + - { dataset: IPSL-CM6A-LR, ensemble: "r(10:11)i1p1f1", grid: gr } + - { dataset: IPSL-CM6A-LR, ensemble: r14i1p1f1, grid: gr } + - { dataset: IPSL-CM6A-LR, ensemble: r22i1p1f1, grid: gr } + - { dataset: IPSL-CM6A-LR, ensemble: r25i1p1f1, grid: gr } + - {dataset: MIROC6, ensemble: "r(1:3)i1p1f1", grid: gn} + - {dataset: MIROC6, ensemble: "r(21:29)i1p1f1", grid: gn} + - {dataset: MIROC6, ensemble: "r(31:34)i1p1f1", grid: gn} + - {dataset: MIROC6, ensemble: "r(36:37)i1p1f1", grid: gn} + - {dataset: MIROC6, ensemble: "r(39:40)i1p1f1", grid: gn} + - {dataset: MIROC6, ensemble: "r(42:47)i1p1f1", grid: gn} + - {dataset: MIROC6, ensemble: "r(49:50)i1p1f1", grid: gn} + - {dataset: MPI-ESM1-2-LR, ensemble: "r(1:3)i1p1f1", grid: gn} + - {dataset: MPI-ESM1-2-LR, ensemble: "r(5:10)i1p1f1", grid: gn} + - {dataset: MRI-ESM2-0, ensemble: "r(1:5)i1p1f1", grid: gn} + - {dataset: NESM3, ensemble: r1i1p1f1, grid: gn} + - {dataset: NorESM2-LM, ensemble: "r(1:3)i1p1f1", grid: gn} + - {dataset: NorESM2-MM, ensemble: r2i1p1f1, grid: gn} + - {dataset: TaiESM1, ensemble: r1i1p1f1, grid: gn} + +dataset_obs: &dataset_obs + - {dataset: ERA5, project: OBS6, type: reanaly, version: v1, tier: 3} + +diagnostics: + evaluation_bc: + description: Evaluation of model TXx variability in Bristish Columbia + variables: + models: + short_name: tasmax + exp: [historical, ssp245] + timerange: '19400101/20223112' + mip: day + project: CMIP6 + preprocessor: preproc_model_bc + additional_datasets: *datasets_cmip + reference: + short_name: tasmax + timerange: '19400101/20223112' + mip: day + preprocessor: preproc_era_bc + additional_datasets: *dataset_obs + scripts: + timeseries: + script: eccc_extremes/simple_timeseries.py + mpl_style: malinina24 + color_style: malinina24 + model_weighting: False + units: C + region: BC # needed for the figure title + var_label: TXx # needed for axis name + figure_caption: 'TXx anomalies in BC relative to 1940-1959' # optional + variability: + script: eccc_extremes/variability_evaluation.py + mpl_style: malinina24 + color_style: malinina24 + model_weighting: True + units: C + region: BC # needed for the figure title + var_label: TXx # needed for axis name + figure_caption: 'Standard deviations of TXx anomalies in BC relative to 1940-1959' # optional + evaluation_sk: + description: Evaluation of model TXx variability in Bristish Columbia + variables: + models: + short_name: tasmax + exp: [historical, ssp245] + timerange: '19400101/20223112' + mip: day + project: CMIP6 + preprocessor: preproc_model_sk + additional_datasets: *datasets_cmip + reference: + short_name: tasmax + timerange: '19400101/20223112' + mip: day + preprocessor: preproc_era_sk + additional_datasets: *dataset_obs + scripts: + timeseries: + script: eccc_extremes/simple_timeseries.py + mpl_style: malinina24 + color_style: malinina24 + model_weighting: False + units: C + region: SK # needed for the figure title + var_label: TXx # needed for axis name + figure_caption: 'TXx anomalies in SK relative to 1940-1959' # optional + variability: + script: eccc_extremes/variability_evaluation.py + mpl_style: malinina24 + color_style: malinina24 + model_weighting: True + units: C + region: SK # needed for the figure title + var_label: TXx # needed for axis name + figure_caption: 'Standard deviations of TXx anomalies in SK relative to 1940-1959' # optional \ No newline at end of file From 652bb4af2b5e3809a48e90f4f5c47c693449e46c Mon Sep 17 00:00:00 2001 From: Elizaveta Malinina Date: Fri, 9 May 2025 22:55:00 +0000 Subject: [PATCH 09/12] add docstrings --- .../eccc_extremes/simple_timeseries.py | 180 +++++++++++------- .../eccc_extremes/variability_evaluation.py | 175 ++++++++++------- .../recipe_extremes_evaluation_SK_BC.yml | 29 ++- 3 files changed, 232 insertions(+), 152 deletions(-) diff --git a/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py b/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py index 2fd614bb2d..c23e1fd6df 100644 --- a/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py +++ b/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py @@ -1,3 +1,25 @@ +""" +Diagnostic to reproduce fig. 8 from Malinina&Gillett (2024) + +This diagnostic plots simple timeseries for individual models and +Multi-Model mean, weighted by the number of realizations in each model. +If there is more than one realization the min/maxshading is showed as +as a shading around the line. + +The reqired variable groups are `reference` and other data variable +group. The preprocessed data should be a 1D timeseries cube. + +Optional keywords for diagnostic are: + `units`: unit flag which will be used for plotting + `color_style`: name of teh color style file + `mpl_style`: name of the matplotlib style file + `figure_caption`: figure caption to be used + `region`: name of the region + +Author: Elizaveta Malinina (CCCma) + elizaveta.malinina-rieger@ec.gc.ca +""" + import logging import os @@ -6,6 +28,7 @@ import numpy as np import esmvaltool.diag_scripts.shared.plot as eplot + # import internal esmvaltool modules here from esmvaltool.diag_scripts.shared import ( group_metadata, @@ -37,23 +60,23 @@ def obtain_reference(data_group: list): reference_dic = {} - ref_fnames = group_metadata(data_group, 'filename') + ref_fnames = group_metadata(data_group, "filename") mins = list() maxs = list() for dataset_f in ref_fnames.keys(): - dataset_n = ref_fnames[dataset_f][0]['dataset'] - if len(group_metadata(data_group, 'dataset')[dataset_n]) > 1: - key = ref_fnames[dataset_f][0]['alias'] + dataset_n = ref_fnames[dataset_f][0]["dataset"] + if len(group_metadata(data_group, "dataset")[dataset_n]) > 1: + key = ref_fnames[dataset_f][0]["alias"] else: key = dataset_n ref_cb = iris.load_cube(dataset_f) - # make key dataset if there's only one filename + # make key dataset if there's only one filename # of this dataset otherwise alias reference_dic[key] = ref_cb - mins.append(ref_cb.collapsed('time', iris.analysis.MIN).data) - maxs.append(ref_cb.collapsed('time', iris.analysis.MAX).data) + mins.append(ref_cb.collapsed("time", iris.analysis.MIN).data) + maxs.append(ref_cb.collapsed("time", iris.analysis.MAX).data) ref_max = np.asarray(maxs).max() ref_min = np.asarray(mins).min() @@ -64,16 +87,21 @@ def obtain_reference(data_group: list): def create_provenance(caption: str): """Creates provenance dictionary.""" provenance_dic = { - 'authors': ['malinina_elizaveta'], - 'caption': caption, - 'references': ['malinina24'] + "authors": ["malinina_elizaveta"], + "caption": caption, + "references": ["malinina24"], } return provenance_dic -def plot_timeseries(data_dic: dict, reference_dic: dict, cfg: dict, - min_val: float, max_val: float): +def plot_timeseries( + data_dic: dict, + reference_dic: dict, + cfg: dict, + min_val: float, + max_val: float, +): """ This function plots timeseries and the max/min values. @@ -94,7 +122,7 @@ def plot_timeseries(data_dic: dict, reference_dic: dict, cfg: dict, b_max = np.ceil(max_val * 1.05) b_min = np.floor(min_val * 1.05) - mpl_st_file = eplot.get_path_to_mpl_style(cfg.get('mpl_style')) + mpl_st_file = eplot.get_path_to_mpl_style(cfg.get("mpl_style")) plt.style.use(mpl_st_file) for dataset in data_dic.keys(): @@ -106,8 +134,10 @@ def plot_timeseries(data_dic: dict, reference_dic: dict, cfg: dict, for i in range(0, len(data_dic[dataset])): var_cb = data_dic[dataset][i] model_var_data.append(var_cb.data) - weights.append(var_cb.attributes['ensemble_weight'] * - var_cb.attributes['reverse_dtsts_n']) + weights.append( + var_cb.attributes["ensemble_weight"] + * var_cb.attributes["reverse_dtsts_n"] + ) weights = np.array(weights) model_var_data = np.array(model_var_data) mean_var_arr = np.average(model_var_data, axis=0) @@ -115,57 +145,64 @@ def plot_timeseries(data_dic: dict, reference_dic: dict, cfg: dict, max_var_arr = np.max(model_var_data, axis=0) years = [ - var_cb.coord('time').cell(i).point.year - for i in range(var_cb.coord('time').shape[0]) + var_cb.coord("time").cell(i).point.year + for i in range(var_cb.coord("time").shape[0]) ] - color_st = eplot.get_dataset_style(dataset, cfg.get('color_style')) - ax_ts.plot(years, - mean_var_arr, - c=color_st['color'], - zorder=3, - label=dataset) + color_st = eplot.get_dataset_style(dataset, cfg.get("color_style")) + ax_ts.plot( + years, mean_var_arr, c=color_st["color"], zorder=3, label=dataset + ) if len(data_dic[dataset]) > 1: - ax_ts.fill_between(years, - min_var_arr, - max_var_arr, - color=color_st['color'], - lw=0, - alpha=0.25) + ax_ts.fill_between( + years, + min_var_arr, + max_var_arr, + color=color_st["color"], + lw=0, + alpha=0.25, + ) for ref_data in reference_dic.keys(): ref_cb = reference_dic.get(ref_data) - ref_color_st = eplot.get_dataset_style(ref_data, - cfg.get('color_style')) + ref_color_st = eplot.get_dataset_style( + ref_data, cfg.get("color_style") + ) ref_years = [ - ref_cb.coord('time').cell(i).point.year - for i in range(ref_cb.coord('time').shape[0]) + ref_cb.coord("time").cell(i).point.year + for i in range(ref_cb.coord("time").shape[0]) ] - ax_ts.plot(ref_years, - ref_cb.data, - label=ref_data, - c=ref_color_st['color'], - zorder=3) - - ax_ts.axhline(color='grey', zorder=1) + ax_ts.plot( + ref_years, + ref_cb.data, + label=ref_data, + c=ref_color_st["color"], + zorder=3, + ) + + ax_ts.axhline(color="grey", zorder=1) ax_ts.set_ylim(b_min, b_max) ax_ts.legend(loc=0, ncols=4, fancybox=False, frameon=False) - variable = cfg['var_label'] if cfg.get( - 'var_label') else var_cb.var_name - exp_variable = variable.replace('_', ' ') - units = cfg['units'] if cfg.get('units') else var_cb.units - region = cfg['region'] if cfg.get('units') else 'region' + variable = ( + cfg["var_label"] if cfg.get("var_label") else var_cb.var_name + ) + exp_variable = variable.replace("_", " ") + units = cfg["units"] if cfg.get("units") else var_cb.units + region = cfg["region"] if cfg.get("units") else "region" - ax_ts.set_ylabel(f'{exp_variable}, {units}') - ax_ts.set_xlabel('year') + ax_ts.set_ylabel(f"{exp_variable}, {units}") + ax_ts.set_xlabel("year") - default_caption = f'Timeseries of {exp_variable} in {region}' + default_caption = f"Timeseries of {exp_variable} in {region}" - caption = cfg['figure_caption'] if cfg.get( - 'figure_caption') else default_caption + caption = ( + cfg["figure_caption"] + if cfg.get("figure_caption") + else default_caption + ) - caption += f' ({dataset})' + caption += f" ({dataset})" fig_ts.suptitle(caption) @@ -173,8 +210,9 @@ def plot_timeseries(data_dic: dict, reference_dic: dict, cfg: dict, plt.tight_layout() - fig_path = os.path.join(cfg['plot_dir'], - f'timeseries_{variable}_{region}_{dataset}') + fig_path = os.path.join( + cfg["plot_dir"], f"timeseries_{variable}_{region}_{dataset}" + ) save_figure(fig_path, prov_dic, cfg, fig_ts, close=True) @@ -183,14 +221,14 @@ def plot_timeseries(data_dic: dict, reference_dic: dict, cfg: dict, def main(cfg): """Initiates diagnostic""" - input_data = cfg['input_data'] + input_data = cfg["input_data"] - groups = group_metadata(input_data.values(), 'variable_group', sort=True) + groups = group_metadata(input_data.values(), "variable_group", sort=True) maxs, mins = [], [] data_dic = {} - reference_dic, ref_max, ref_min = obtain_reference(groups.pop('reference')) + reference_dic, ref_max, ref_min = obtain_reference(groups.pop("reference")) maxs.append(ref_max) mins.append(ref_min) @@ -199,34 +237,36 @@ def main(cfg): for k in groups.keys(): remaining_metadata.extend(groups[k]) - datasets = group_metadata(remaining_metadata, 'dataset') + datasets = group_metadata(remaining_metadata, "dataset") ens_var_cubelist = iris.cube.CubeList() for dataset in datasets.keys(): - filepaths = list(group_metadata(datasets[dataset], 'filename').keys()) + filepaths = list(group_metadata(datasets[dataset], "filename").keys()) n_real = len(filepaths) mod_var_cubelist = iris.cube.CubeList() for filepath in filepaths: mod_cb = iris.load_cube(filepath) # adding weights to the data cubes - mod_cb.attributes['ensemble_weight'] = 1 / n_real - mod_cb.attributes['reverse_dtsts_n'] = 1 / len(datasets) + mod_cb.attributes["ensemble_weight"] = 1 / n_real + mod_cb.attributes["reverse_dtsts_n"] = 1 / len(datasets) ens_var_cubelist.append(mod_cb) mod_var_cubelist.append(mod_cb) - mins.append(mod_cb.collapsed('time', iris.analysis.MIN).data) - maxs.append(mod_cb.collapsed('time', iris.analysis.MAX).data) + mins.append(mod_cb.collapsed("time", iris.analysis.MIN).data) + maxs.append(mod_cb.collapsed("time", iris.analysis.MAX).data) data_dic[dataset] = mod_var_cubelist - data_dic['Multi-Model'] = ens_var_cubelist + data_dic["Multi-Model"] = ens_var_cubelist - plot_timeseries(data_dic, - reference_dic, - cfg, - min_val=np.asarray(mins).min(), - max_val=np.asarray(maxs).max()) + plot_timeseries( + data_dic, + reference_dic, + cfg, + min_val=np.asarray(mins).min(), + max_val=np.asarray(maxs).max(), + ) - logger.info('Success') + logger.info("Success") -if __name__ == '__main__': +if __name__ == "__main__": # always use run_diagnostic() to get the config (the preprocessor # nested dictionary holding all the needed information) with run_diagnostic() as config: diff --git a/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py b/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py index def0914ba0..fca218fd2c 100644 --- a/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py +++ b/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py @@ -1,3 +1,30 @@ +""" +Diagnostic to reproduce fig. 6 from Malinina&Gillett (2024) + +This diagnostic calculates a standard deviation for individual model +realizations and if requested weighted by the number of realizations +standard deviation for all the models. The resulting plot shows +indivisual model realizations as dots and all models as the square. +The reference dataset is shown as the vertical line. + +The reqired variable groups are `reference` and other data variable +group. The preprocessed data should be a 1D cube. + +The mandatory keywords for diagnostic are required: + `model_weighting`: boolean flag showing if multi-model std_dev + should be weighted + `units`: unit flag which will be used for plotting + +Optional keywords for diagnostic are: + `color_style`: name of teh color style file + `mpl_style`: name of the matplotlib style file + `figure_caption`: figure caption to be used + `region`: name of the region + +Author: Elizaveta Malinina (CCCma) + elizaveta.malinina-rieger@ec.gc.ca +""" + import logging import os @@ -7,6 +34,7 @@ import numpy as np import esmvaltool.diag_scripts.shared.plot as eplot + # import internal esmvaltool modules here from esmvaltool.diag_scripts.shared import ( group_metadata, @@ -35,12 +63,12 @@ def obtain_reference(data_group: list): reference_dic = {} - ref_fnames = group_metadata(data_group, 'filename') + ref_fnames = group_metadata(data_group, "filename") for dataset_f in ref_fnames.keys(): - dataset_n = ref_fnames[dataset_f][0]['dataset'] - if len(group_metadata(data_group, 'dataset')[dataset_n]) > 1: - key = ref_fnames[dataset_f][0]['alias'] + dataset_n = ref_fnames[dataset_f][0]["dataset"] + if len(group_metadata(data_group, "dataset")[dataset_n]) > 1: + key = ref_fnames[dataset_f][0]["alias"] else: key = dataset_n ref_cb = iris.load_cube(dataset_f) @@ -53,16 +81,17 @@ def create_provenance(caption: str): """Creates provenance dictionary.""" provenance_dic = { - 'authors': ['malinina_elizaveta'], - 'caption': caption, - 'references': ['malinina24'] + "authors": ["malinina_elizaveta"], + "caption": caption, + "references": ["malinina24"], } return provenance_dic -def calculate_ensemble_stdev(ens_cubelist: iris.cube.CubeList, - weight_flag: bool): +def calculate_ensemble_stdev( + ens_cubelist: iris.cube.CubeList, weight_flag: bool +): """This function calculates standard deviation of multi model ensemble. Parameters: @@ -81,8 +110,9 @@ def calculate_ensemble_stdev(ens_cubelist: iris.cube.CubeList, wghts, wdata = [], [] for i in range(0, len(ens_cubelist)): cb = ens_cubelist[i] - cb_wght = cb.attributes['ensemble_weight'] * cb.attributes[ - 'reverse_dtsts_n'] + cb_wght = ( + cb.attributes["ensemble_weight"] * cb.attributes["reverse_dtsts_n"] + ) wghts.append(np.full(cb.shape, cb_wght)) wdata.append(ens_cubelist[i].data) @@ -91,7 +121,7 @@ def calculate_ensemble_stdev(ens_cubelist: iris.cube.CubeList, if weight_flag: mmm = np.average(wdata, weights=wghts) - ens_std = np.sqrt(np.average((wdata - mmm)**2, weights=wghts)) + ens_std = np.sqrt(np.average((wdata - mmm) ** 2, weights=wghts)) else: ens_std = wdata.std() @@ -111,79 +141,90 @@ def plot_stdevs(data_dic, reference_dic, cfg): config dictionary, comes from ESMValCore """ - mpl_st_file = eplot.get_path_to_mpl_style(cfg.get('mpl_style')) + mpl_st_file = eplot.get_path_to_mpl_style(cfg.get("mpl_style")) plt.style.use(mpl_st_file) fig_stds, ax_stds = plt.subplots(nrows=1, ncols=1) - fig_stds.set_size_inches(8., 9.) + fig_stds.set_size_inches(8.0, 9.0) fig_stds.set_dpi(200) y_ticks = np.arange(0, len(data_dic.keys())) - y_labs = np.zeros(len(data_dic.keys()), dtype=' Date: Fri, 9 May 2025 23:05:59 +0000 Subject: [PATCH 10/12] fix whitespace --- .../diag_scripts/shared/plot/styles_python/malinina24.yml | 8 ++++---- .../plot/styles_python/matplotlib/malinina24.mplstyle | 2 +- esmvaltool/references/malinina24.bibtex | 2 +- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/esmvaltool/diag_scripts/shared/plot/styles_python/malinina24.yml b/esmvaltool/diag_scripts/shared/plot/styles_python/malinina24.yml index a6631f3a2f..eded076054 100644 --- a/esmvaltool/diag_scripts/shared/plot/styles_python/malinina24.yml +++ b/esmvaltool/diag_scripts/shared/plot/styles_python/malinina24.yml @@ -21,21 +21,21 @@ ERA5: facecolor: '#cd5c5c' mark: 'o' thick: 1 -nat: +nat: avgstd: 1 color: '#005000' dash: '-' facecolor: '#005000' mark: 'o' thick: 1 -factual: +factual: avgstd: 1 color: '#c57a00' dash: '-' facecolor: '#c57a00' mark: 'o' thick: 1 -ssp245: +ssp245: avgstd: 1 color: '#4677c0' dash: '-' @@ -48,4 +48,4 @@ default: dash: '-' facecolor: '#1A0F50' mark: o - thick: 1 \ No newline at end of file + thick: 1 diff --git a/esmvaltool/diag_scripts/shared/plot/styles_python/matplotlib/malinina24.mplstyle b/esmvaltool/diag_scripts/shared/plot/styles_python/matplotlib/malinina24.mplstyle index f262b98295..810071bf15 100644 --- a/esmvaltool/diag_scripts/shared/plot/styles_python/matplotlib/malinina24.mplstyle +++ b/esmvaltool/diag_scripts/shared/plot/styles_python/matplotlib/malinina24.mplstyle @@ -25,4 +25,4 @@ ytick.direction : out xtick.major.width : 0.5 ytick.major.width : 0.5 xtick.minor.width : 0.4 -ytick.minor.width : 0.4 \ No newline at end of file +ytick.minor.width : 0.4 diff --git a/esmvaltool/references/malinina24.bibtex b/esmvaltool/references/malinina24.bibtex index f5fd7d7437..e68a718b31 100644 --- a/esmvaltool/references/malinina24.bibtex +++ b/esmvaltool/references/malinina24.bibtex @@ -7,4 +7,4 @@ year={2024}, doi={10.1016/j.wace.2024.100642}, publisher={Elsevier} -} \ No newline at end of file +} From 8269ae85d6ed9999a37076617234339733521664 Mon Sep 17 00:00:00 2001 From: Elizaveta Malinina Date: Fri, 9 May 2025 23:39:28 +0000 Subject: [PATCH 11/12] fixing codacy errors --- esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py | 6 ++---- .../diag_scripts/eccc_extremes/variability_evaluation.py | 1 + 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py b/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py index c23e1fd6df..89214b4fcb 100644 --- a/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py +++ b/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py @@ -62,8 +62,7 @@ def obtain_reference(data_group: list): ref_fnames = group_metadata(data_group, "filename") - mins = list() - maxs = list() + maxs, mins = [], [] for dataset_f in ref_fnames.keys(): dataset_n = ref_fnames[dataset_f][0]["dataset"] @@ -128,8 +127,7 @@ def plot_timeseries( for dataset in data_dic.keys(): fig_ts, ax_ts = plt.subplots(nrows=1, ncols=1) - model_var_data = list() - weights = list() + model_var_data, weights = [], [] for i in range(0, len(data_dic[dataset])): var_cb = data_dic[dataset][i] diff --git a/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py b/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py index fca218fd2c..f88869d522 100644 --- a/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py +++ b/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py @@ -232,6 +232,7 @@ def plot_stdevs(data_dic, reference_dic, cfg): def main(cfg): + """Initiates diagnostic""" input_data = cfg["input_data"] groups = group_metadata(input_data.values(), "variable_group", sort=True) From 0d0e01a950007cb3ae81e9fcb8982b444b6546ce Mon Sep 17 00:00:00 2001 From: Elizaveta Malinina Date: Fri, 9 May 2025 23:46:11 +0000 Subject: [PATCH 12/12] fixing codacy errors --- .../eccc_extremes/simple_timeseries.py | 32 ++++++++--------- .../eccc_extremes/variability_evaluation.py | 36 +++++++++---------- 2 files changed, 34 insertions(+), 34 deletions(-) diff --git a/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py b/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py index 89214b4fcb..f76daed30d 100644 --- a/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py +++ b/esmvaltool/diag_scripts/eccc_extremes/simple_timeseries.py @@ -43,18 +43,18 @@ def obtain_reference(data_group: list): """This function cretes a dictionary with reference data. - Parameters: - ----------- - data_group: + Parameters + ---------- + data_group : list with metadata for reference group - Returns: - -------- - reference_dic: + Returns + ------- + reference_dic : dictionary with the reference data. Dataset name is a keyword. - ref_max: + ref_max : maximum value for all reference datasets - ref_min: + ref_min : minimum value for all reference datasets """ @@ -104,17 +104,17 @@ def plot_timeseries( """ This function plots timeseries and the max/min values. - Parameters: - ----------- - data_dic: + Parameters + ---------- + data_dic : dictionary with the data to plot - reference_dic: + reference_dic : dictionary with the refrence data - cfg: + cfg : config dictionary, comes from ESMValCore - min_val: + min_val : minimum value for the xlims - max_val: + max_val : maximum value for the xlims """ @@ -218,7 +218,7 @@ def plot_timeseries( def main(cfg): - """Initiates diagnostic""" + """Initiates diagnostic.""" input_data = cfg["input_data"] groups = group_metadata(input_data.values(), "variable_group", sort=True) diff --git a/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py b/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py index f88869d522..992c12c716 100644 --- a/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py +++ b/esmvaltool/diag_scripts/eccc_extremes/variability_evaluation.py @@ -50,14 +50,14 @@ def obtain_reference(data_group: list): """This function cretes a dictionary with reference data. - Parameters: - ----------- - data_group: + Parameters + ---------- + data_group : list with metadata for reference group - Returns: - -------- - reference_dic: + Returns + ------- + reference_dic : dictionary with the reference data. Dataset name is a keyword """ @@ -94,16 +94,16 @@ def calculate_ensemble_stdev( ): """This function calculates standard deviation of multi model ensemble. - Parameters: - ----------- - ens_cubelist: + Parameters + ---------- + ens_cubelist : CubeList with the ensemble timeseries - weight_flag: + weight_flag : flag for weighting or not the data - Returns: - -------- - ens_std: float + Returns + ------- + ens_std : float standard deviation for the ensemble """ @@ -131,13 +131,13 @@ def calculate_ensemble_stdev( def plot_stdevs(data_dic, reference_dic, cfg): """This function plots timeseries and the max/min values. - Parameters: - ----------- - data_dic: + Parameters + ---------- + data_dic : dictionary with the data to plot - reference_dic: + reference_dic : dictionary with the reference data - cfg: + cfg : config dictionary, comes from ESMValCore """