Recipe simple diagnostics for Southern Ocean/Antarctica

doc/sphinx/source/recipes/recipe_ocean_south.rst

@@ -0,0 +1,70 @@
+.. _recipes_ocean_south:
+
+Simple southern ocean diagnostics
+=====
+
+Overview
+--------
+
+A few simple recipes recreated in ESMValTool for running on bulk datasets e.g. CMIP6, from the `COSIMA cookbook <https://cosima-recipes.readthedocs.io/en/latest/index.html>`_.
+
+
+Available recipes and diagnostics
+---------------------------------
+
+Recipes are stored in esmvaltool/recipes/
+
+* recipe_ocean_south.yml
+
+Diagnostics are stored in esmvaltool/diag_scripts/ocean_south/
+
+* plot_barotropic.py: for plotting barotropic streamfunction
+* plot_hovmoller.py: for plotting Hovmoller diagrams of ocean temperature and salinity
+* plot_kinetic_energy.py: for plotting ocean kinetic energy
+
+
+User settings in recipe
+-----------------------
+
+*None*
+
+
+Variables
+---------
+
+* umo (ocean, monthly mean, longitude latitude time)
+* so (ocean, monthly mean, longitude latitude time depth)
+* thetao (ocean, monthly mean, longitude latitude time depth)
+* uo (ocean, monthly mean, longitude latitude time depth)
+* vo (ocean, monthly mean, longitude latitude time depth)
+
+
+
+References
+----------
+
+* https://cosima-recipes.readthedocs.io/en/latest/02-Appetisers/Barotropic_Streamfunction.html
+* https://cosima-recipes.readthedocs.io/en/latest/02-Appetisers/Hovmoller_Temperature_Depth.html
+* https://cosima-recipes.readthedocs.io/en/latest/03-Mains/Eddy-Mean_Kinetic_Energy_Decomposition.html
+
+
+Example plots
+-------------
+
+.. _fig_1:
+.. figure::  /recipes/figures/ocean_south/hovmoller_ACCESS-OM2.png
+   :align:   center
+
+   Hovmoller diagram of ocean temperature and salinity for the ACCESS-OM2 model.
+
+.. _fig_2:
+.. figure::  /recipes/figures/ocean_south/kinetic_energy_ACCESS-OM2-025.png
+   :align:   center
+
+   Kinetic energy around Antarctica for the ACCESS-OM2-025 model.
+
+.. _fig_3:
+.. figure::  /recipes/figures/ocean_south/barotropic_streamfunction_ACCESS-OM2.png
+   :align:   center
+
+   Barotropic streamfunction around Antarctica for the ACCESS-OM2 model.

esmvaltool/diag_scripts/ocean_south/plot_barotropic.py

@@ -0,0 +1,110 @@
+"""diagnostic script to plot"""
+
+import iris
+import os
+import logging
+import numpy as np
+import iris.plot as iplt
+import matplotlib.path as mpath
+import matplotlib.pyplot as plt
+import cartopy.crs as ccrs
+import cmocean
+from esmvaltool.diag_scripts.shared import (
+    run_diagnostic,
+    save_figure,
+    group_metadata,
+    select_metadata,
+)
+from esmvalcore.preprocessor import (
+    climate_statistics,
+    regrid
+)
+
+
+# This part sends debug statements to stdout
+logger = logging.getLogger(os.path.basename(__file__))
+
+
+def circumpolar_map():
+    fig = plt.figure(figsize = (12, 8))
+    ax = plt.axes(projection = ccrs.SouthPolarStereo())
+    ax.set_extent([-180, 180, -80, -40], crs = ccrs.PlateCarree())
+    ax.set_facecolor('lightgrey')
+    # Map the plot boundaries to a circle
+    theta = np.linspace(0, 2 * np.pi, 100)
+    center, radius = [0.5, 0.5], 0.5
+    verts = np.vstack([np.sin(theta), np.cos(theta)]).T
+    circle = mpath.Path(verts * radius + center)
+    ax.set_boundary(circle, transform = ax.transAxes)
+
+    return fig, ax
+
+
+def get_provenance_record(caption, ancestor_files):
+    """Create a provenance record describing the diagnostic data and plot."""
+
+    record = {
+        "caption": caption,
+        "statistics": ["other"],
+        "domains": ["polar"],
+        "plot_types": ["map"],
+        "authors": [
+            "chun_felicity",
+        ],
+        "references": [
+            "access-nri",
+        ],
+        "ancestors": ancestor_files,
+    }
+    return record
+
+
+def main(cfg):
+    """# in script cumsum, units handling, climate mean"""
+
+    input_data = cfg["input_data"].values()
+
+    for dataset in input_data:
+        # Load the data
+        input_file = dataset['filename']
+        name = dataset['dataset']
+        cube = iris.load_cube(input_file)
+
+        # units cumulative sum and climate mean
+        p0 = 1035.0
+        cube = cube.collapsed('depth', iris.analysis.SUM)
+        cube.data = cube.data / p0  # divide by density for volume
+        cube.data = cube.data / 1e6  # then divide by 10^6 for Sv
+        cube = regrid(cube, target_grid='0.5x0.5', scheme="linear")
+        cube.data = cube.data.cumsum(1)  #latitude index 1
+        logger.info(cube.summary(shorten=True))
+        cube.units = 'Sv'
+        cube = climate_statistics(cube, period='full', operator='mean')
+
+        # plot
+        fig, ax = circumpolar_map()
+
+        p1 = iplt.contourf(cube,
+                    levels = np.arange(-50, 150, 10),
+                    extend = 'both',
+                    cmap = cmocean.cm.speed,
+                    )
+        plt.colorbar(p1, orientation='vertical', label='$\psi$ (Sv)')
+
+        plt.title(f'Barotropic streamfunction from {name}')
+        # Save output
+        prov_record = get_provenance_record(
+            f'Barotropic streamfunction from {name}.',
+            [input_file],
+        )
+        save_figure(
+            f"barotropic_streamfunction_{name}",
+            prov_record,
+            cfg,
+            figure=fig,
+            dpi=300,
+            )
+
+if __name__ == "__main__":
+    with run_diagnostic() as config:
+        main(config)

esmvaltool/diag_scripts/ocean_south/plot_hovmoller.py

@@ -0,0 +1,177 @@
+"""diagnostic script to plot"""
+
+import iris
+import os
+import logging
+import numpy as np
+import iris.plot as iplt
+from matplotlib.gridspec import GridSpec
+from matplotlib import ticker
+import matplotlib.pyplot as plt
+
+import cmocean
+from esmvaltool.diag_scripts.shared import (
+    run_diagnostic,
+    save_figure,
+    group_metadata,
+    select_metadata,
+)
+from esmvalcore.preprocessor import (
+    area_statistics,
+    anomalies
+)
+
+
+# This part sends debug statements to stdout
+logger = logging.getLogger(os.path.basename(__file__))
+
+
+def plot_hovmoller(fsize = 14):
+    """Setup of figures properties."""
+    plt.rcParams['font.size'] = fsize
+    plt.rcParams['xtick.labelsize'] = fsize-2
+    plt.rcParams['ytick.labelsize'] = fsize-2
+
+    fig = plt.figure(figsize = (8, 6))
+    grid = GridSpec(100, 100)
+
+    ax = [fig.add_subplot(grid[:30, :45]), fig.add_subplot(grid[:30, 48:93]),
+          fig.add_subplot(grid[32:, :45]), fig.add_subplot(grid[32:, 48:93])]
+
+    for i in range(len(ax)):
+        # ax[i].xaxis.set_major_formatter(date_format)
+        ax[i].tick_params(axis='x', labelrotation=45)
+        if i != 0 and i != 2:
+            ax[i].set_yticklabels([]) 
+
+    return fig, ax
+
+
+def hovmoller(ds_cube, startyr):
+    """Calculate anomalies summed over area for depth vs time plot."""
+    cube = anomalies(ds_cube, period='full', reference={'start_year':startyr, 'start_month':1,'start_day':1,
+                                                            'end_year':startyr, 'end_month':12, 'end_day':31})
+    # anomalies doesn't retain cell measure required for area statistics
+    cube.add_cell_measure(ds_cube.cell_measure('cell_area'), data_dims=[2,3])
+    # sum anomalies
+    cube = area_statistics(cube, operator='sum')
+
+    #divide by total area cube (anomalies are summed.. weighted by cell_area(multiplied)
+    tarea = totalarea_depth(ds_cube)
+    return cube / tarea
+
+
+def totalarea_depth(cube):
+    """Calculate total area at each depth level."""
+    tmp1 = cube[0] #get first time step
+    # divide by itself for 1s mask 3-dimensional
+    tmp1.data = tmp1.data / tmp1.data
+    # check cell_measure exists?
+    if cube.cell_measure('cell_area'):
+        return area_statistics(tmp1, operator='sum')
+
+
+def plot_details(ax, cf_temp, cf_salt):
+    for i in range(len(ax)):
+        if i < 2:
+            ax[i].set_ylim(500, 0)
+            ax[i].set_xticklabels([])
+        else:
+            ax[i].set_xlabel("")
+            ax[i].set_ylim(5000, 500)
+    ax[0].set_ylabel("Depth [m]")
+    ax[1].set_ylabel("")
+    ax[2].set_ylabel("Depth [m]")
+    ax[3].set_ylabel("")
+
+    # Colorbars
+    bar = plt.axes([0.11, 0.97, 0.35, 0.02])
+    cbar_1 = plt.colorbar(cf_temp, cax = bar, orientation = 'horizontal', extend='both', format= '%.2f')  
+    cbar_1.set_label("Temperature [$\degree$C]")
+
+    bar = plt.axes([0.50, 0.97, 0.35, 0.02])
+    cbar_2 = plt.colorbar(cf_salt, cax = bar, orientation = 'horizontal', extend='both', format= '%.2f')  
+    cbar_2.set_label("Salinity [psu]")
+
+    for cbar in [cbar_1, cbar_2]:
+        tick_locator = ticker.MaxNLocator(nbins=3, prune='both') ## The ticker needs to called within the loop
+        cbar.locator = tick_locator
+        cbar.update_ticks()
+
+def get_provenance_record(caption, ancestor_files):
+    """Create a provenance record describing the diagnostic data and plot."""
+
+    record = {
+        "caption": caption,
+        "statistics": ["other"],
+        "domains": ["global"],
+        "plot_types": ["vert"],
+        "authors": [
+            "chun_felicity",
+        ],
+        "references": [
+            "access-nri",
+        ],
+        "ancestors": ancestor_files,
+    }
+    return record
+
+
+def main(cfg):
+    """Create Hovmoller diagrams."""
+
+    input_data = cfg["input_data"].values()
+
+    # group by dataset
+    ds_groups = group_metadata(
+        input_data,
+        "dataset",
+    )
+
+    levels_tempsal = {'thetao':np.arange(-0.3, 0.31, 0.01), 'so':np.arange(-0.03, 0.031, 0.001)}
+    cmapcm = {'thetao':cmocean.cm.balance, 'so':cmocean.cm.curl}
+    axi = {'thetao':0, 'so':1}
+    # for each ds plt ts, sal 
+    for grp, var_attr in ds_groups.items():
+        fig, ax = plot_hovmoller(fsize = 14)
+        files=[]
+        for metadata in var_attr:  #[thetao, so]
+            name = metadata["dataset"]
+            shortname = metadata["short_name"]
+            input_file = metadata['filename']
+            cube = iris.load_cube(metadata['filename'])
+            hov_cube = hovmoller(cube, metadata['start_year'])
+
+            # plot twice
+            for i in [axi[shortname], axi[shortname]+2]:
+                p1 = iplt.contourf(hov_cube,
+                            levels = levels_tempsal[shortname],
+                            extend = 'both',
+                            cmap = cmapcm[shortname],
+                            axes = ax[i], #0,2 temp
+                            )
+            files.append(input_file)
+
+            if shortname == 'so':
+                cf_salt = p1
+            else:
+                cf_temp = p1
+
+        plot_details(ax, cf_temp, cf_salt)
+        # Save output
+        prov_record = get_provenance_record(
+            f'Depth-Time Temperature and Salinity from {name}.',
+            files,
+        )
+        save_figure(
+            f"hovmoller_{name}",
+            prov_record,
+            cfg,
+            figure=fig,
+            dpi=300,
+            )
+
+
+if __name__ == "__main__":
+    with run_diagnostic() as config:
+        main(config)