diff --git a/doc/sphinx/source/recipes/figures/sea_surface_salinity/radar.png b/doc/sphinx/source/recipes/figures/sea_surface_salinity/radar.png deleted file mode 100644 index dbecac5107..0000000000 Binary files a/doc/sphinx/source/recipes/figures/sea_surface_salinity/radar.png and /dev/null differ diff --git a/doc/sphinx/source/recipes/figures/sea_surface_salinity/radar_bias.png b/doc/sphinx/source/recipes/figures/sea_surface_salinity/radar_bias.png new file mode 100644 index 0000000000..61b312624e Binary files /dev/null and b/doc/sphinx/source/recipes/figures/sea_surface_salinity/radar_bias.png differ diff --git a/doc/sphinx/source/recipes/figures/sea_surface_salinity/radar_std.png b/doc/sphinx/source/recipes/figures/sea_surface_salinity/radar_std.png new file mode 100644 index 0000000000..4eb8a8e2cc Binary files /dev/null and b/doc/sphinx/source/recipes/figures/sea_surface_salinity/radar_std.png differ diff --git a/doc/sphinx/source/recipes/recipe_sea_surface_salinity.rst b/doc/sphinx/source/recipes/recipe_sea_surface_salinity.rst index 3e9f357c8c..688502bf03 100644 --- a/doc/sphinx/source/recipes/recipe_sea_surface_salinity.rst +++ b/doc/sphinx/source/recipes/recipe_sea_surface_salinity.rst @@ -9,8 +9,11 @@ Overview This recipe compares the regional means of sea surface salinity with a reference dataset (ESACCI-SEA-SURFACE-SALINITY v1 or v2 by default). To do this, the recipe generates plots for the timeseries of each region and -a radar plot showing the correlation between dataset and reference timeseries for -each region during the time they both exist. +a radar plot showing (i) the mean state bias, and (ii) the ratio between the +simulated and observed standard deviations of different regional averages of +sea surface salinity, calculated in the temporal window for which observations +and simulations overlap. + Preprocessor requirements: -------------------------- @@ -202,8 +205,16 @@ References Example plots ------------- -.. figure:: /recipes/figures/sea_surface_salinity/radar.png +.. figure:: /recipes/figures/sea_surface_salinity/radar_bias.png + :align: center + + Radar plot showing the mean state biases (simulation minus observations) + for the regional averages of sea surface salinity in the selected + ocean basins and seas. + +.. figure:: /recipes/figures/sea_surface_salinity/radar_std.png :align: center - Radar plot showing correlation of average sea surface salinity for multiple - regions with the observations + Radar plot showing the ratio between the simulated and observed standard deviations + of the regional averages of sea surface salinity in the selected + ocean basins and seas. diff --git a/esmvaltool/diag_scripts/sea_surface_salinity/compare_salinity.py b/esmvaltool/diag_scripts/sea_surface_salinity/compare_salinity.py index fb308e975c..e1fa384297 100644 --- a/esmvaltool/diag_scripts/sea_surface_salinity/compare_salinity.py +++ b/esmvaltool/diag_scripts/sea_surface_salinity/compare_salinity.py @@ -9,10 +9,9 @@ import iris.quickplot as qplot import matplotlib.pyplot as plt import numpy as np -from esmvalcore.preprocessor._time import regrid_time -from iris.analysis.stats import pearsonr +from esmvalcore.iris_helpers import date2num +from esmvalcore.preprocessor import climate_statistics, regrid_time from iris.coord_categorisation import add_month_number, add_year -from iris.util import unify_time_units from matplotlib.legend import Legend from matplotlib.legend_handler import HandlerBase from matplotlib.text import Text @@ -27,6 +26,12 @@ class CompareSalinity(object): def __init__(self, config): self.cfg = config + self.ticks = { + 'mean': [-0.5, 0.0, 0.5], + 'std_dev': [0.25, 0.5, 1, 2, 4] + } + self.lim = {'mean': [-1, 1], 'std_dev': [0.01, 10]} + self.operation = {'mean': 'bias', 'std_dev': 'std_ratio'} def compute(self): data = group_metadata(self.cfg[names.INPUT_DATA].values(), @@ -50,7 +55,7 @@ def compute(self): time_coord.units.name, calendar='gregorian', ) - unify_time_units((reference, dataset)) + self._unify_time_coordinates([reference, dataset]) logger.debug("Info dataset %s:", alias) logger.debug(dataset) ancestors = (dataset_info[names.FILENAME], reference_ancestor) @@ -92,64 +97,77 @@ def create_radar_plot(self, data_info, data, reference, reference_alias, add_month_number(data, 'time') add_year(data, 'time') - data.remove_coord('time') add_month_number(reference, 'time') add_year(reference, 'time') - reference.remove_coord('time') data_alias = data_info[names.ALIAS] - corr = pearsonr(data, reference, ('month_number', 'year')) - angles = np.linspace(0, 2 * np.pi, corr.shape[0] + 1) - # Initialise the spider plot - ax = plt.subplot(111, polar=True) - for spine in ax.spines.values(): - spine.set_color('grey') - - # Draw one axe per variable + add labels labels yet - letters = [string.ascii_uppercase[i] for i in range(0, corr.shape[0])] - plt.xticks(angles[:-1], - letters, - color='grey', - size=8, - rotation=angles[:-1]) - - # Draw ylabels - ax.set_rlabel_position(0) - plt.yticks([0.25, 0.5, 0.75], ["0.25", "0.5", "0.75"], - color="grey", - size=7) - plt.ylim(0, 1) - - data = np.append(corr.data, corr.data[0]) - - more_angles = np.linspace(0, 2 * np.pi, corr.shape[0] * 20 + 1) - interp_data = np.interp(more_angles, angles, data) - - # Plot data - ax.plot(more_angles, interp_data, linewidth=1, linestyle='solid') - ax.fill(more_angles, interp_data, 'b', alpha=0.1) - ax.legend(letters, - corr.coord('shape_id').points, - loc='upper center', - ncol=2, - frameon=False, - bbox_to_anchor=(0.5, -0.1), - borderaxespad=0.) - plt.title( - f'{data_info[names.SHORT_NAME]} correlation\n' - f'{data_alias} vs {reference_alias}', - pad=20) - plt.tight_layout() - plot_path = os.path.join( - self.cfg[names.PLOT_DIR], - f"{data_info[names.SHORT_NAME]}_comparison_{data_alias}_" - f"{reference_alias}.{self.cfg[names.OUTPUT_FILE_TYPE]}") - plt.savefig(plot_path) - plt.close() - caption = (f"Correlation comparison in diferent regions for " - f"{data_alias} and {reference_alias}") - self._create_prov_record(plot_path, caption, ancestors) + for operator in ['mean', 'std_dev']: + climat_ref = climate_statistics(reference, operator) + climat_data = climate_statistics(data, operator) + if operator == 'mean': + result_data = climat_ref.data - climat_data.data + else: + result_data = climat_ref.data / climat_data.data + + result = climat_ref.copy(result_data) + angles = np.linspace(0, 2 * np.pi, result.shape[0] + 1) + # Initialise the spider plot + ax = plt.subplot(111, polar=True) + for spine in ax.spines.values(): + spine.set_color('grey') + + # Draw one axe per variable + add labels labels yet + letters = [ + string.ascii_uppercase[i] for i in range(0, result.shape[0]) + ] + plt.xticks(angles[:-1], + letters, + color='grey', + size=8, + rotation=angles[:-1]) + + # Draw ylabels + ax.set_rlabel_position(0) + plt.yticks(self.ticks[operator], + list(map(str, self.ticks[operator])), + color="grey", + size=7) + plt.ylim(min(self.lim[operator]), max(self.lim[operator])) + + radar_data = np.append(result.data, result.data[0]) + more_angles = np.linspace(0, 2 * np.pi, result.shape[0] * 20 + 1) + interp_data = np.interp(more_angles, angles, radar_data) + + # Plot data + ax.plot(more_angles, interp_data, linewidth=1, linestyle='solid') + ax.fill(more_angles, interp_data, 'b', alpha=0.1) + ax.legend(letters, + result.coord('shape_id').points, + loc='upper center', + ncol=2, + frameon=False, + bbox_to_anchor=(0.5, -0.1), + borderaxespad=0.) + if operator == 'std_dev': + ax.set_yscale('symlog', linthresh=0.1) + operation = self.operation[operator] + plt.title( + f'{data_info[names.SHORT_NAME]} {operation}\n' + f'{data_alias} vs {reference_alias}', + pad=20) + plt.tight_layout() + plot_path = os.path.join( + self.cfg[names.PLOT_DIR], + f"{data_info[names.SHORT_NAME]}_{operation}" + f"_comparison_{data_alias}_" + f"{reference_alias}.{self.cfg[names.OUTPUT_FILE_TYPE]}") + plt.savefig(plot_path) + plt.close() + caption = ( + f"Absolute {operation} comparison in different regions for " + f"{data_alias} and {reference_alias}") + self._create_prov_record(plot_path, caption, ancestors) def _create_prov_record(self, filepath, caption, ancestors): record = { @@ -172,8 +190,10 @@ def _get_time_offset(self, time_unit): return time_offset def _align_yearly_axes(self, cube): - """ - Perform a time-regridding operation to align time axes for yr data. + """Align years. + + Perform a time-regridding operation to align time axes for yr + data. """ years = [cell.point.year for cell in cube.coord('time').cells()] # be extra sure that the first point is not in the previous year @@ -203,14 +223,12 @@ def _datetime_to_int_days(self, cube): return days def _get_overlap(self, cubes): - """ - Get discrete time overlaps. - This method gets the bounds of coord time - from the cube and assembles a continuous time - axis with smallest unit 1; then it finds the - overlaps by doing a 1-dim intersect; - takes the floor of first date and - ceil of last date. + """Get discrete time overlaps. + + This method gets the bounds of coord time from the cube and + assembles a continuous time axis with smallest unit 1; then it + finds the overlaps by doing a 1-dim intersect; takes the floor + of first date and ceil of last date. """ all_times = [] for cube in cubes: @@ -224,10 +242,9 @@ def _get_overlap(self, cubes): return time_bounds_list def _slice_cube(self, cube, t_1, t_2): - """ - Efficient slicer. - Simple cube data slicer on indices - of common time-data elements. + """Efficient slicer. + + Simple cube data slicer on indices of common time-data elements. """ time_pts = [t for t in cube.coord('time').points] converted_t = self._datetime_to_int_days(cube) @@ -237,6 +254,37 @@ def _slice_cube(self, cube, t_1, t_2): ]) return [idxs[0], idxs[-1]] + @staticmethod + def _get_consistent_time_unit(cubes): + """Fix time units. + + Return cubes' time unit if consistent, standard calendar + otherwise. + """ + t_units = [cube.coord('time').units for cube in cubes] + if len(set(t_units)) == 1: + return t_units[0] + return cf_units.Unit("days since 1850-01-01", calendar="standard") + + def _unify_time_coordinates(self, cubes): + """Make sure all cubes' share the same time coordinate.""" + t_unit = self._get_consistent_time_unit(cubes) + for cube in cubes: + # Extract date info from cube + coord = cube.coord('time') + years = [p.year for p in coord.units.num2date(coord.points)] + months = [p.month for p in coord.units.num2date(coord.points)] + dates = [ + datetime(year, month, 15, 0, 0, 0) + for year, month in zip(years, months) + ] + + # Update the cubes' time coordinate + cube.coord('time').points = date2num(dates, t_unit, coord.dtype) + cube.coord('time').units = t_unit + cube.coord('time').bounds = None + cube.coord('time').guess_bounds() + class TextHandler(HandlerBase): def create_artists(self, legend, text, xdescent, ydescent, width, height, diff --git a/esmvaltool/recipes/recipe_sea_surface_salinity.yml b/esmvaltool/recipes/recipe_sea_surface_salinity.yml index 63c5e30452..f3cc71b355 100644 --- a/esmvaltool/recipes/recipe_sea_surface_salinity.yml +++ b/esmvaltool/recipes/recipe_sea_surface_salinity.yml @@ -32,8 +32,6 @@ preprocessors: ids: - Arctic Ocean - Southern Ocean - - Mediterranean Sea - Eastern Basin - - Mediterranean Sea - Western Basin - North Atlantic Ocean - South Atlantic Ocean - North Pacific Ocean @@ -45,8 +43,16 @@ preprocessors: areacello: datasets: - - {project: CMIP6, exp: historical, dataset: MPI-ESM1-2-HR, ensemble: r1i1p1f1, - start_year: 1950, end_year: 2014, alias: MPI-ESM1-2-HR} + - &cmip6 {project: CMIP6, exp: historical, dataset: ACCESS-CM2, ensemble: r1i1p1f1, + start_year: 1950, end_year: 2014, alias: ACCESS-CM2} + - {<<: *cmip6, dataset: CMCC-CM2-HR4, alias: CMCC-CM2-HR4} + - {<<: *cmip6, dataset: CanESM5, alias: CanESM5} + - {<<: *cmip6, dataset: IPSL-CM6A-LR, alias: IPSL-CM6A-LR} + - {<<: *cmip6, dataset: MIROC6, alias: MIROC6} + - {<<: *cmip6, dataset: MPI-ESM1-2-HR, alias: MPI-ESM1-2-HR} + - {<<: *cmip6, dataset: NorESM2-MM, alias: NorESM2-MM} + - {<<: *cmip6, dataset: GISS-E2-2-H, alias: GISS-E2-2-H, institute: NASA-GISS} + diagnostics: compare_salinity: