Add recipe for sea ice area and extents in southern polar region

.zenodo.json

@@ -386,6 +386,11 @@
       "name": "Bonnet, Pauline",
       "orcid": "0000-0003-3780-0784"
     },
+    {
+      "affiliation": "ACCESS-NRI, Australia",
+      "name": "Chun, Felicity",
+      "orcid": "0009-0007-0845-0953"
+    },
     {
       "affiliation": "MetOffice, UK",
       "name": "Munday, Gregory",

CITATION.cff

@@ -391,6 +391,11 @@ authors:
     family-names: Bonnet
     given-names: Pauline
     orcid: "https://orcid.org/0000-0003-3780-0784"
+  -
+    affiliation: "ACCESS-NRI, Australia"
+    family-names: Chun
+    given-names: Felicity
+    orcid: "https://orcid.org/0009-0007-0845-0953"
   - 
     affiliation: "MetOffice, UK"
     family-names: Munday

doc/sphinx/source/recipes/index.rst

@@ -147,6 +147,7 @@ Other
    recipe_rainfarm
    recipe_seaice
    recipe_seaice_drift
+   recipe_seaice_extents_sh
    recipe_seaice_feedback
    recipe_shapeselect
    recipes_testing

doc/sphinx/source/recipes/recipe_seaice_extents_sh.rst

@@ -0,0 +1,72 @@
+.. _recipes_seaice_extents_sh:
+
+Sea Ice area and extents
+========================
+
+Overview
+--------
+
+This recipe plots sea ice concentration from CICE (sea ice model) output 
+around the southern polar region and compares it to the NSIDC CDR 
+(National Snow and Ice Data Centre, Climate Data Record) dataset.
+
+
+Available recipes and diagnostics
+---------------------------------
+
+Recipes are stored in esmvaltool/recipes/
+
+    * recipe_seaice_extents_sh.yml
+
+Diagnostics are stored in esmvaltool/diag_scripts/seaice_area_extents/
+
+    * seaicearea_trends.py: plot minima and maxima sea ice area trends
+    * seaice_mapextents.py: plot sea ice extent and differences to Observations
+
+
+User settings in recipe
+-----------------------
+
+#. Script seaice_mapextents.py
+
+   *Required settings for script*
+
+   * `months`: months by month number which the mean are to be plotted
+
+
+Variables
+---------
+
+* siconc (seaIce, monthly, longitude latitude time)
+* areacello (fx)
+
+
+Observations and reformat scripts
+---------------------------------
+
+*Note: (1) obs4MIPs data can be used directly without any preprocessing;
+(2) see headers of reformat scripts for non-obs4MIPs data for download
+instructions.*
+
+* NSIDC CDR sh (siconc - esmvaltool/cmorizers/data/formatters/datasets/nsidc_g02202_sh.py)
+
+
+References
+----------
+
+https://cosima-recipes.readthedocs.io/en/latest/Examples/Sea_Ice_Area_Concentration_Volume_with_Obs.html
+
+Example plots
+-------------
+
+.. _trends:
+.. figure::  /recipes/figures/seaice_extents_sh/min_trend.png
+   :align:   center
+
+   Minima trends of sea ice area with observation data. ACCESS OM model data years from 0, added 1652 years to model years for comparability.
+
+.. _map extents:
+.. figure::  /recipes/figures/seaice_extents_sh/map_difference.png
+   :align:   center
+
+   Difference and extents of models with observations for selected months.

esmvaltool/config-references.yml

@@ -157,6 +157,11 @@ authors:
     name: Chen, Jack
     institute: NCAR, USA
     orcid:
+  chun_felicity:
+    name: Chun, Felicity
+    institute: ACCESS-NRI, Australia
+    orcid: https://orcid.org/0009-0007-0845-0953
+    github: flicj191
   cionni_irene:
     name: Cionni, Irene
     institute: ENEA, Italy
@@ -722,6 +727,11 @@ authors:
     name: Stevens, Mark
     institute: NCAR, US
     orcid:
+  steketee_anton:
+    name: Steketee, Anton
+    institute: ACCESS-NRI
+    orcid: https://orcid.org/0009-0002-9081-4106
+    github: anton-seaice
   # Former viewers (not active viewers)
   adeniyi_kemisola:
     name: Adeniyi, Kemisola
@@ -771,6 +781,7 @@ authors:
 
 projects:
   4c: EU H2020 project 4C
+  access-nri: The Australian Earth System Simulator (ACCESS-NRI)
   applicate: EU Horizon 2020 Advanced prediction in polar regions and beyond
   c3s-magic: Copernicus Climate Change Service 34a Lot 2 (MAGIC) project
   climval: BMBF MiKlip Project ClimVal

esmvaltool/diag_scripts/seaice_area_extents/seaice_mapextents.py

@@ -0,0 +1,125 @@
+"""Diagnostic script to plot extent and differences.
+
+based on code from Anton Steketee's COSIMA recipes notebook
+https://cosima-recipes.readthedocs.io/en/latest/Examples/Sea_Ice_Area_Concentration_Volume_with_Obs.html
+"""
+
+import calendar
+import logging
+import os
+
+import iris
+import iris.plot as iplt
+import matplotlib.lines as mlines
+import matplotlib.pyplot as plt
+import numpy as np
+from cartopy.crs import SouthPolarStereo
+from esmvalcore.preprocessor import extract_month
+
+from esmvaltool.diag_scripts.shared import (
+    group_metadata,
+    select_metadata,
+    run_diagnostic,
+    save_figure,
+)
+
+# This part sends debug statements to stdout
+logger = logging.getLogger(os.path.basename(__file__))
+
+
+def map_fig_diff(model_dict, obs_si, months):
+    """Create map with model dictionary: labels, cubes."""
+    # figure set up, width 9 for 2 models
+    figure = plt.figure(figsize=(9, len(months) * 3.5))
+    j = 0  # to iterate through positions on figure
+
+    for mon in months:  # eg.[2,9]
+        i = 1
+        obs_cube = extract_month(obs_si, mon)
+        for mod_label, mod_si in model_dict.items():
+            mod_cube = extract_month(mod_si, mon)
+
+            out = mod_cube.copy()
+            diff = mod_cube.data - obs_cube.data
+            out.data = diff
+
+            plt.subplot(len(months), 3, i + j * 3,
+                        projection=SouthPolarStereo(true_scale_latitude=-70))
+
+            diffmap = iplt.contourf(out, levels=np.arange(-90, 91, 20),
+                                    cmap='RdBu')
+
+            iplt.contour(obs_cube, levels=[15], colors=['yellow'])
+            iplt.contour(mod_cube, levels=[15],
+                         linewidths=1.0, colors=['black'])
+
+            plt.title(calendar.month_abbr[mon] + ' ' + mod_label)
+
+            i += 1
+        j += 1
+
+    line_cdr = mlines.Line2D([], [], color='yellow', label="Observed Extent")
+    line_mod = mlines.Line2D([], [], color='black', label="Modelled Extent")
+
+    plt.legend(handles=[line_cdr, line_mod], loc='center left',
+               bbox_to_anchor=(1.2, 0.5))
+    cax = plt.axes([0.7, 0.55, 0.04, 0.3])
+    _ = plt.colorbar(diffmap, cax=cax,
+                     label='Difference in \nSea Ice Concentration')
+
+    plt.subplots_adjust(left=0.05, bottom=0.05,
+                        right=0.95, top=0.95,
+                        wspace=0.05, hspace=0.05)
+
+    return figure
+
+
+def main(cfg):
+    """Compute."""
+    # Get the preprocessed data that we will use as input.
+    input_data = cfg['input_data'].values()
+
+    groups = group_metadata(input_data, 'variable_group')
+    # assign obs_si - select
+    selection = select_metadata(input_data, project='OBS6')
+    obs_si = iris.load_cube(selection[0]['filename'])
+
+    for group_name in groups.keys():
+        mod_dict = {}
+        ancestor_filels = []
+        logger.info("Processing variable %s", group_name)
+
+        for attr in groups[group_name]:
+            if attr['project'].startswith('OBS'):
+                logger.info("Load OBS dataset %s", attr['dataset'])
+                obs_si = iris.load_cube(attr['filename'])
+            else:
+                logger.info("load model dataset %s", attr['dataset'])
+                mod_dict[attr['dataset']] = iris.load_cube(attr['filename'])
+            ancestor_filels.append(attr['filename'])
+
+        logger.info("creating map differences")
+        mapfig = map_fig_diff(mod_dict, obs_si, cfg['months'])
+
+        provenance_record = get_provenance_record(ancestor_filels)
+        save_figure(group_name, provenance_record, cfg, figure=mapfig)
+
+
+def get_provenance_record(ancestor_files):
+    """Build provenance dictionary."""
+    record = {
+        'ancestors': ancestor_files,
+        'authors': ['chun_felicity', 'steketee_anton'],
+        'caption': 'siconc observations difference from model',
+        'domains': ['shpolar'],
+        'plot_types': ['polar'],
+        'references': [],
+        'statistics': ['diff'],
+    }
+    return record
+
+
+if __name__ == '__main__':
+
+    with run_diagnostic() as config:
+        main(config)

esmvaltool/diag_scripts/seaice_area_extents/seaicearea_trends.py

@@ -0,0 +1,78 @@
+"""Diagnostic script to plot minima and maxima trends.
+
+based on code from Anton Steketee's COSIMA cookbook notebook
+https://cosima-recipes.readthedocs.io/en/latest/Examples/Sea_Ice_Area_Concentration_Volume_with_Obs.html
+"""
+
+import logging
+import os
+
+import iris
+import matplotlib.pyplot as plt
+from iris import quickplot
+
+from esmvaltool.diag_scripts.shared import run_diagnostic, save_figure
+
+# This part sends debug statements to stdout
+logger = logging.getLogger(os.path.basename(__file__))
+
+
+def _prom_dim_coord(cube, _field, _filename):
+    iris.util.promote_aux_coord_to_dim_coord(cube, 'year')
+
+
+def plot_trends(datagroup, provenance_record, cfg):
+    """Create plot for min and max groups."""
+    for variable_group, attributes in datagroup.items():
+        plt.clf()
+        for (filep, vname, datalabel, offset) in attributes:
+            cube = iris.load_cube(filep, vname, _prom_dim_coord)
+            if offset:
+                cube.coord('year').points = [y + offset for y in
+                                             cube.coord('year').points]
+            quickplot.plot(cube, label=datalabel)
+
+        plt.title(f"Trends in Sea-Ice {variable_group.split('_')[1]}ima")
+        plt.legend(loc='upper left')
+        plt.ylabel('Sea-Ice Area (km2)')
+
+        save_figure(variable_group, provenance_record, cfg, dpi=300)
+
+
+def main(cfg):
+    """Compute sea ice area for each input dataset."""
+    provenance_record = {
+        'caption': "sea ice trends southern hemisphere",
+        'authors': [
+            'chun_felicity',
+            'steketee_anton'
+        ],
+        'references': [''],
+        'ancestors': list(cfg['input_data'].keys()),
+    }
+    input_data = cfg['input_data'].values()
+
+    datagroup = {}  # for each variable min and max
+
+    for dataset in input_data:
+        # Load the data
+        if 'offset_years' in dataset:
+            input_file = (dataset['filename'], dataset['short_name'],
+                          dataset['dataset'], dataset['offset_years'])
+        else:
+            input_file = (dataset['filename'], dataset['short_name'],
+                          dataset['dataset'], None)
+        # key for different models
+        logger.info("Dataset: %s", {dataset['long_name']})
+        if dataset['variable_group'] not in datagroup:
+            datagroup[dataset['variable_group']] = []
+        datagroup[dataset['variable_group']].append(input_file)
+
+    logger.info(datagroup)
+    plot_trends(datagroup, provenance_record, cfg)
+
+
+if __name__ == '__main__':
+
+    with run_diagnostic() as config:
+        main(config)