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
+----------
+
+* COSIMA(Consortium for Ocean-Sea Ice Modelling in Australia) recipe: https://cosima-recipes.readthedocs.io/en/latest/DocumentedExamples/SeaIce_Obs_Model_Compare.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,10 @@ authors:
     name: Stevens, Mark
     institute: NCAR, US
     orcid:
+  steketee_anton:
+    name: Steketee, Anton
+    institute: ACCESS-NRI
+    orcid:
   # Former viewers (not active viewers)
   adeniyi_kemisola:
     name: Adeniyi, Kemisola
@@ -771,6 +780,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,161 @@
+"""Diagnostic script to plot extent and differences.
+
+based on code from Anton Steketee's COSIMA cookbook notebook
+https://cosima-recipes.readthedocs.io/en/latest/DocumentedExamples
+    /SeaIce_Obs_Model_Compare.html
+"""
+
+import logging
+import os
+import calendar
+from cartopy.crs import SouthPolarStereo
+import xarray as xr
+import xesmf
+import matplotlib.pyplot as plt
+import matplotlib.lines as mlines
+import numpy as np
+import pandas as pd
+
+from esmvaltool.diag_scripts.shared import (run_diagnostic, save_figure,
+                                            save_data)
+from esmvaltool.diag_scripts.shared._base import get_plot_filename
+
+# This part sends debug statements to stdout
+logger = logging.getLogger(os.path.basename(__file__))
+
+
+def model_regrid_diff(mod_si, obs_si, cdr, mon, latmax):
+    """Regrid and compute difference grid."""
+    lat_min = mod_si.lat.min().values.item()
+    mod_si = mod_si.where(mod_si.lat < latmax, drop=True)
+    regrd_out = obs_si.where(obs_si.lat > lat_min, drop=True)
+
+    # regrid model data to observations
+    regridder_access_sh = xesmf.Regridder(mod_si.isel(time=0).drop(['i', 'j']),
+                                          regrd_out.isel(time=0),
+                                          'bilinear',
+                                          periodic=True,
+                                          unmapped_to_nan=True)
+
+    model_mean = mod_si.siconc.sel(
+        time=mod_si.siconc.time.dt.month.isin(mon)).mean('time')
+    mod_regrid = regridder_access_sh(model_mean)
+    diff_ds = mod_regrid - cdr
+
+    return diff_ds, mod_regrid
+
+
+def map_diff(mod_si_ls, obs_si, months, cfg, prov):
+    """Create figure mapping extents for models and months."""
+    # get lat max for regridding
+    latmax = obs_si.lat.max().values.item()
+
+    proj = SouthPolarStereo(true_scale_latitude=-70)
+
+    # fig set up, width for 2 models, check len mod_si_ls
+    figure = plt.figure(figsize=(9, len(months) * 4))
+    j = 0  # to iterate through positions on figure
+
+    for mon in months:
+        cdr = obs_si.siconc.sel(
+            time=obs_si.siconc.time.dt.month.isin(mon)).mean('time')
+        i = 1
+        for mod_label, mod_si in mod_si_ls.items():
+
+            diff_ds, mod_regrid = model_regrid_diff(mod_si, obs_si, cdr, mon,
+                                                    latmax)
+            # save plot data
+            save_data(''.join([mod_label, calendar.month_abbr[mon], '_mean']), prov,
+                      cfg, mod_regrid.to_iris())
+            save_data(''.join([mod_label, calendar.month_abbr[mon], '_obs_diff']), prov,
+                      cfg, diff_ds.to_iris())
+
+            axes = plt.subplot(len(months), 3, i + j * 3, projection=proj)
+
+            diffmap = axes.contourf(diff_ds.x,
+                                    diff_ds.y,
+                                    diff_ds,
+                                    levels=np.arange(-90, 91, 20),
+                                    cmap='RdBu')
+            cs_cdr = cdr.plot.contour(levels=[15], ax=axes)
+            cs_mod = mod_regrid.plot.contour(levels=[15],
+                                             ax=axes,
+                                             colors=['black'])
+
+            plt.title(' '.join([calendar.month_abbr[mon], mod_label]))
+
+            i += 1
+        j += 1
+
+    line_cdr = mlines.Line2D([], [],
+                             color=cs_cdr.collections[0].get_edgecolor(),
+                             label="Observed Extent")
+
+    line_mod = mlines.Line2D([], [],
+                             color=cs_mod.collections[0].get_edgecolor(),
+                             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')
+    return figure
+
+
+def main(cfg):
+    """Compute."""
+    # Get a description of the preprocessed data that we will use as input.
+    input_data = cfg['input_data'].values()
+    data = []
+
+    # Find input datasets to use
+    for dataset in input_data:
+
+        input_file = [dataset['filename'], dataset['dataset']]
+        # drop areacello dataset for map
+        if dataset['short_name'] == 'siconc':
+            data.append(input_file)
+
+    inputs_df = pd.DataFrame(data, columns=['filename', 'dataset'])
+
+    logger.info("input siconc data:", inputs_df)
+    mod_si_dict = {}
+    for filepath, data_name in inputs_df.itertuples(index=False):
+        if data_name == 'NSIDC-G02202-sh':
+            # Load the data
+            obs_si = xr.open_dataset(filepath)
+        else:
+            mod_si_dict[data_name] = xr.open_dataset(filepath)
+
+    logger.info("creating map differences")
+    provenance_record = get_provenance_record(inputs_df['filename'].to_list())
+
+    mapfig = map_diff(mod_si_dict, obs_si, cfg['months'], cfg,
+                      provenance_record)
+    # Save output
+    output_path = get_plot_filename('map_difference', cfg)
+
+    save_figure(output_path, provenance_record, cfg, figure=mapfig)
+
+
+def get_provenance_record(ancestor_files):
+    """Build provenance dictionary."""
+    record = {
+        'ancestors': ancestor_files,
+        'authors': ['chun_felicity', 'steketee_anton'],
+        'caption': '',
+        'domains': ['shpolar'],
+        'plot_types': ['polar'],
+        'references': [],
+        'statistics': ['diff'],
+    }
+    return record
+
+
+if __name__ == '__main__':
+
+    with run_diagnostic() as config:
+        main(config)