Adding figures of Bock and Lauer (2024)

doc/sphinx/source/recipes/index.rst

@@ -94,6 +94,7 @@ Future projections
    recipe_tcr
    recipe_tebaldi21esd
    recipe_climate_change_hotspot
+   recipe_bock24acp
 
 IPCC
 ^^^^

doc/sphinx/source/recipes/recipe_bock24acp.rst

@@ -0,0 +1,142 @@
+.. _recipes_bock24acp:
+
+Cloud properties and their projected changes in CMIP models with low to high climate sensitivity
+================================================================================================
+
+Overview
+--------
+
+The recipes recipe_bock24acp_*.yml reproduce figures (Fig. 3, 4, 6 and 7) from the publication Bock and Lauer, 2024.
+
+
+Available recipes and diagnostics
+---------------------------------
+
+Recipes are stored in recipes/clouds
+
+    * recipe_bock24acp_fig3-4_maps.yml
+    * recipe_bock24acp_fig6_zonal.yml
+    * recipe_bock24acp_fig7_boxplots.yml
+
+Diagnostics are stored in diag_scripts/
+
+    Fig. 3 and 4:
+
+    * clouds/clouds_ecs_groups_maps.py: Geographical maps of the multi-year annual means for group means of historical CMIP simulations from all three ECS groups.
+
+    Fig. 6:
+
+    * clouds/clouds_ecs_groups_zonal.py: Zonally averaged group means.
+
+    Fig. 7:
+
+    * clouds/clouds_ecs_groups_boxplots.py: Boxplots of relative changes for all groups.
+
+
+User settings in recipe
+-----------------------
+
+#. Script clouds_ecs_groups_maps.py 
+
+   *Required settings (scripts)*
+
+   reference: if true, a reference dataset is given within 'variable_group' equal 'OBS'
+
+   *Optional settings (scripts)*
+
+   plot_each_model: one figure for each single model
+
+
+#. Script clouds/clouds_ecs_groups_zonal.py
+
+   *Required settings (scripts)*
+
+   group_by: list of 'variable_group's to have the order
+   plot_type: 'zonal' and 'height' plots are available 
+
+   *Optional settings (scripts)*
+
+   filename_attach: attachment to the output files
+
+
+#. Script clouds/clouds_ecs_groups_boxplots.py
+
+   *Required settings (scripts)*
+
+   exclude_datasets: list of datasets which are not used for the statistics, default is ['MultiModelMean', 'MultiModelP5', 'MultiModelP95']
+   group_by: list of 'variable_group's to have the order
+   plot_type: 'zonal' and 'height' plots are available 
+
+   *Optional settings (scripts)*
+
+   filename_attach: attachment to the output files
+   title: set title of figure
+   y_range: set range of the y-axes
+
+
+Variables
+---------
+
+* clt (atmos, monthly, longitude latitude time)
+* iwp (atmos, monthly, longitude latitude time)
+* lwp (atmos, monthly, longitude latitude time)
+* rlut (atmos, monthly, longitude latitude time)
+* rsut (atmos, monthly, longitude latitude time)
+* rlutcs (atmos, monthly, longitude latitude time)
+* rsutcs (atmos, monthly, longitude latitude time)
+* tas (atmos, monthly, longitude latitude time)
+
+
+Observations and reformat scripts
+---------------------------------
+
+* CERES-EBAF - CERES TOA radiation fluxes (used for calculation of
+  cloud forcing)
+
+  *Reformat script:* cmorizers/data/formatters/datasets/ceres_ebaf.py
+
+
+References
+----------
+
+* Bock, L. and Lauer, A.: Cloud properties and their projected changes in CMIP
+  models with low to high climate sensitivity, Atmos. Chem. Phys., 24, 1587–1605,
+  https://doi.org/10.5194/acp-24-1587-2024, 2024.
+
+
+Example plots
+-------------
+
+.. _fig_bock24acp_1:
+.. figure::  /recipes/figures/bock24acp/map_netcre.png
+   :align:   center
+
+   Geographical map of the multi-year annual mean net cloud radiative effect from 
+   (a) CERES–EBAF Ed4.2 (OBS) and (b–d) group means of historical CMIP simulations
+   from all three ECS groups (Fig. 4).
+
+.. _fig_bock24acp_2:
+.. figure::  /recipes/figures/bock24acp/zonal_diff_clt_ssp585.png
+   :align:   center
+
+   The upper panel show the zonally averaged group means of total cloud
+   fraction from historical simulations (solid lines)
+   and RCP8.5/SSP5-8.5 scenarios (dashed lines) for the three different ECS groups.
+   Lower panels show the corresponding relative differences of all zonally
+   averaged group means between the RCP8.5/SSP5-8.5 scenarios and the corresponding
+   historical simulations. Shading indicates the 5 % and 95 % quantiles of the single
+   model results (Fig. 6a).
+
+.. _fig_bock24acp_3:
+.. figure::  /recipes/figures/bock24acp/boxplot_ssp585_south_oc.png
+   :align:   center
+
+   Relative change (calculated as the difference between the scenario value and the
+   historical value divided by the historical value) of total cloud fraction (clt),
+   ice water path (iwp), liquid water path (lwp), and net cloud radiative effect
+   (netcre) per degree of warming averaged over the Southern Ocean (30–65∘ S). In the
+   box plot, each box indicates the range from the first
+   quartile to the third quartile, the vertical line shows the median, and the
+   whiskers the minimum and maximum values, excluding the outliers. Outliers are
+   defined as being outside 1.5 times the interquartile range (Fig. 7b).
+

esmvaltool/diag_scripts/clouds/clouds_ecs_groups_boxplots.py

@@ -0,0 +1,245 @@
+"""Python diagnostic for plotting boxplots for different regions."""
+import logging
+from pathlib import Path
+
+import iris
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import seaborn as sns
+
+from esmvaltool.diag_scripts.shared import (
+    ProvenanceLogger,
+    group_metadata,
+    run_diagnostic,
+    get_diagnostic_filename,
+    get_plot_filename,
+    select_metadata,
+)
+
+logger = logging.getLogger(Path(__file__).stem)
+
+VAR_NAMES = {
+    'cl': 'cloud_fraction',
+    'cli': 'ice_water_content',
+    'clw': 'liquid_water_content',
+}
+PALETTE = {
+    'high ECS': 'royalblue',
+    'med ECS': 'green',
+    'low ECS': 'orange',
+}
+
+
+def get_provenance_record(ancestor_files):
+    """Create a provenance record describing the diagnostic data and plot."""
+    caption = ("Relative change per degree of warming averaged over the"
+               "chosen region.")
+
+    record = {
+        'caption': caption,
+        'statistics': ['mean'],
+        'domains': ['global'],
+        'plot_types': ['zonal'],
+        'authors': [
+            'bock_lisa',
+        ],
+        'references': [
+            'bock24acp',
+        ],
+        'ancestors': ancestor_files,
+    }
+    return record
+
+
+def _get_multi_model_mean(cubes, var):
+    """Compute multi-model mean."""
+
+    logger.debug("Calculating multi-model mean")
+    dataset_names = []
+    mmm = {}
+    for (dataset, cube) in cubes.items():
+        dataset_names.append(dataset)
+        mmm['dataset'] = cube.data
+    mmm = np.ma.masked_invalid(list(mmm.values()))
+    mmm_cube = cube.copy(data=np.ma.mean(mmm, axis=0))
+    attributes = {
+        'dataset': 'MultiModelMean',
+        'short_name': var,
+        'datasets': '|'.join(dataset_names),
+    }
+    mmm_cube.attributes = attributes
+    return mmm_cube
+
+
+def read_data(filename):
+    """Compute an example diagnostic."""
+    logger.debug("Loading %s", filename)
+    cube = iris.load_cube(filename)
+
+    if cube.var_name == 'cli':
+        cube.convert_units('g/kg')
+    elif cube.var_name == 'clw':
+        cube.convert_units('g/kg')
+
+    cube = iris.util.squeeze(cube)
+    return cube
+
+
+def compute_diff(filename1, filename2):
+    """Compute difference between two cubes."""
+    logger.debug("Loading %s", filename1)
+    cube1 = iris.load_cube(filename1)
+    cube2 = iris.load_cube(filename2)
+
+    if cube1.var_name == 'cli':
+        cube1.convert_units('g/kg')
+        cube2.convert_units('g/kg')
+    elif cube1.var_name == 'clw':
+        cube1.convert_units('g/kg')
+        cube2.convert_units('g/kg')
+
+    cube = cube2 - cube1
+    cube.metadata = cube1.metadata
+    return cube
+
+
+def compute_diff_temp(input_data, group, var, dataset):
+    """Compute relative change per temperture change."""
+
+    dataset_name = dataset['dataset']
+    var = dataset['short_name']
+
+    input_file_1 = dataset['filename']
+
+    var_data_2 = select_metadata(input_data,
+                                 short_name=var,
+                                 dataset=dataset_name,
+                                 variable_group=var+"_"+group[1])
+    if not var_data_2:
+        raise ValueError(
+            f"No '{var}' data for '{dataset_name}' in '{group[1]}' available")
+
+    input_file_2 = var_data_2[0]['filename']
+
+    tas_data_1 = select_metadata(input_data,
+                                 short_name='tas',
+                                 dataset=dataset_name,
+                                 variable_group='tas_'+group[0])
+    tas_data_2 = select_metadata(input_data,
+                                 short_name='tas',
+                                 dataset=dataset_name,
+                                 variable_group='tas_'+group[1])
+    if not tas_data_1:
+        raise ValueError(
+            f"No 'tas' data for '{dataset_name}' in '{group[0]}' available")
+    if not tas_data_2:
+        raise ValueError(
+            f"No 'tas' data for '{dataset_name}' in '{group[1]}' available")
+    input_file_tas_1 = tas_data_1[0]['filename']
+    input_file_tas_2 = tas_data_2[0]['filename']
+
+    cube = read_data(input_file_1)
+
+    cube_diff = compute_diff(input_file_1, input_file_2)
+    cube_tas_diff = compute_diff(input_file_tas_1, input_file_tas_2)
+
+    cube_diff = (100. * (cube_diff / iris.analysis.maths.abs(cube))
+                 / cube_tas_diff)
+
+    return cube_diff
+
+
+def create_data_frame(input_data, cfg):
+    """Create data frame."""
+    data_frame = pd.DataFrame(columns=['Variable', 'Group', 'Dataset', 'Data'])
+
+    ifile = 0
+
+    all_vars = group_metadata(input_data, 'short_name')
+    groups = group_metadata(input_data, 'variable_group', sort='dataset')
+
+    for var in all_vars:
+        if var != 'tas':
+            logger.info("Processing variable %s", var)
+
+            if var == 'clivi':
+                varname = 'iwp'
+            else:
+                varname = var
+
+            for group_names in cfg['group_by']:
+                logger.info("Processing group %s of variable %s",
+                            group_names[0], var)
+
+                for dataset in groups[var + "_" + group_names[0]]:
+                    dataset_name = dataset['dataset']
+
+                    if dataset_name not in cfg['exclude_datasets']:
+                        cube_diff = compute_diff_temp(input_data, group_names,
+                                                      var, dataset)
+
+                        group_name = group_names[0].split('_')[1] + " ECS"
+
+                        data_frame.loc[ifile] = [varname, group_name,
+                                                 dataset_name, cube_diff.data]
+                        ifile = ifile + 1
+
+    data_frame['Data'] = data_frame['Data'].astype(str).astype(float)
+
+    return data_frame
+
+
+def plot_boxplot(data_frame, cfg):
+
+    sns.set_style('darkgrid')
+    sns.set(font_scale=2)
+    sns.boxplot(data=data_frame, x='Variable', y='Data', hue='Group',
+                palette=PALETTE)
+    plt.ylabel('Relative change (%/K)')
+    if 'y_range' in cfg:
+        plt.ylim(cfg.get('y_range'))
+    plt.title(cfg['title'])
+
+    provenance_record = get_provenance_record(
+            ancestor_files=cfg['input_files'])
+
+    # Save plot
+    plot_path = get_plot_filename('boxplot' + '_' + cfg['filename_attach'],
+                                  cfg)
+    plt.savefig(plot_path)
+    logger.info("Wrote %s", plot_path)
+    plt.close()
+
+    with ProvenanceLogger(cfg) as provenance_logger:
+        provenance_logger.log(plot_path, provenance_record)
+
+
+def main(cfg):
+    """Run diagnostic."""
+    cfg.setdefault('exclude_datasets', ['MultiModelMean', 'MultiModelP5',
+                                        'MultiModelP95'])
+    cfg.setdefault('title', 'Test')
+
+    plt.figure(constrained_layout=True, figsize=(12, 8))
+
+    # Get input data
+    input_data = list(cfg['input_data'].values())
+
+    # Create data frame
+    data_frame = create_data_frame(input_data, cfg)
+
+    # Create plot
+    plot_boxplot(data_frame, cfg)
+
+    # Save file
+    basename = "boxplot_region_" + cfg['filename_attach']
+    csv_path = get_diagnostic_filename(basename, cfg).replace('.nc', '.csv')
+    data_frame.to_csv(csv_path)
+    logger.info("Wrote %s", csv_path)
+
+
+if __name__ == '__main__':
+
+    with run_diagnostic() as config:
+        main(config)