Source code for esmvaltool.diag_scripts.ocean.diagnostic_transects

"""
Transects diagnostics
=====================

Diagnostic to produce images of a transect. These plost show either latitude or
longitude against depth, and the cube value is used as the colour scale.

Note that this diagnostic assumes that the preprocessors do the bulk of the
hard work, and that the cube received by this diagnostic (via the settings.yml
and metadata.yml files) has no time component, and one of the latitude or
longitude coordinates has been reduced to a single value.

An approproate preprocessor for a 3D+time field would be::

  preprocessors:
    prep_transect:
      climate_statistics:
        operator: mean
      extract_transect: # Atlantic Meridional Transect
        latitude: [-50.,50.]
        longitude: 332.

This tool is part of the ocean diagnostic tools package in the ESMValTool.

Author: Lee de Mora (PML)
        ledm@pml.ac.uk

"""
import itertools
import logging
import os
import sys

import iris
import iris.quickplot as qplt
import matplotlib.pyplot as plt
import numpy as np

from esmvaltool.diag_scripts.ocean import diagnostic_tools as diagtools
from esmvaltool.diag_scripts.shared import run_diagnostic
from esmvaltool.diag_scripts.shared._base import ProvenanceLogger

# This part sends debug statements to stdout
logger = logging.getLogger(os.path.basename(__file__))
logging.getLogger().addHandler(logging.StreamHandler(sys.stdout))


[docs] def titlify(title): """ Check whether a title is too long then add it to current figure. Parameters ---------- title: str The title for the figure. """ cutoff = 40 if len(title) > cutoff: # Find good mid point titles = title.split(' ') length = 0 for itr, word in enumerate(titles): length += len(word) if length > cutoff: titles[itr] += '\n' length = 0. title = ' '.join(titles) plt.title(title)
[docs] def determine_transect_str(cube, region=''): """ Determine the Transect String. Takes a guess at a string to describe the transect. Parameters ---------- cube: iris.cube.Cube Input cube to use to determine the transect name. """ if region: return region options = ['latitude', 'longitude'] cube_dims = [c.standard_name for c in cube.coords()] for option in options: if option not in cube_dims: continue coord = cube.coord(option) if len(coord.points) > 1: continue value = coord.points.mean() value = round(value, 2) if option == 'latitude': return str(value) + ' N' if option == 'longitude': if value > 180.: return str(value - 360.) + ' W' return str(value) + ' E' return ''
[docs] def make_depth_safe(cube): """ Make the depth coordinate safe. If the depth coordinate has a value of zero or above, we replace the zero with the average point of the first depth layer. Parameters ---------- cube: iris.cube.Cube Input cube to make the depth coordinate safe Returns ---------- iris.cube.Cube: Output cube with a safe depth coordinate """ depth = cube.coord('depth') # it's fine if depth.points.min() * depth.points.max() > 0.: return cube if depth.attributes['positive'] != 'down': raise Exception('The depth field is not set up correctly') depth_points = [] bad_points = depth.points <= 0. for itr, point in enumerate(depth.points): if bad_points[itr]: depth_points.append(depth.bounds[itr, :].mean()) else: depth_points.append(point) cube.coord('depth').points = depth_points return cube
[docs] def make_cube_region_dict(cube): """ Take a cube and return a dictionairy region: cube. Each item in the dict is a layer with a separate cube for each layer. ie: cubes[region] = cube from specific region Cubes with no region component are returns as: cubes[''] = cube with no region component. This is based on the method diagnostics_tools.make_cube_layer_dict, however, it wouldn't make sense to look for depth layers here. Parameters ---------- cube: iris.cube.Cube the opened dataset as a cube. Returns --------- dict A dictionairy of layer name : layer cube. """ ##### # Check layering: coords = cube.coords() layers = [] for coord in coords: if coord.standard_name in ['region', ]: layers.append(coord) cubes = {} if not layers: cubes[''] = cube return cubes # iris stores coords as a list with one entry: layer_dim = layers[0] if len(layer_dim.points) in [1, ]: cubes[''] = cube return cubes if layer_dim.standard_name == 'region': coord_dim = cube.coord_dims('region')[0] for layer_index, layer in enumerate(layer_dim.points): slices = [slice(None) for index in cube.shape] slices[coord_dim] = layer_index layer = layer.replace('_', ' ').title() cubes[layer] = cube[tuple(slices)] return cubes
[docs] def determine_set_y_logscale(cfg, metadata): """ Determine whether to use a log scale y axis. Parameters ---------- cfg: dict the opened global config dictionairy, passed by ESMValTool. metadata: dict The metadata dictionairy for a specific model. Returns ---------- bool: Boolean to flag whether to plot as a log scale. """ set_y_logscale = True if 'set_y_logscale' in cfg: set_y_logscale = cfg['set_y_logscale'] if 'set_y_logscale' in metadata: set_y_logscale = metadata['set_y_logscale'] return set_y_logscale
[docs] def make_transects_plots( cfg, metadata, filename, ): """ Make a simple plot of the transect for an indivudual model. This tool loads the cube from the file, checks that the units are sensible BGC units, checks for layers, adjusts the titles accordingly, determines the ultimate file name and format, then saves the image. Parameters ---------- cfg: dict the opened global config dictionairy, passed by ESMValTool. metadata: dict The metadata dictionairy for a specific model. filename: str The preprocessed model file. """ # Load cube and set up units cube = iris.load_cube(filename) cube = diagtools.bgc_units(cube, metadata['short_name']) # Is this data is a multi-model dataset? multi_model = metadata['dataset'].find('MultiModel') > -1 cube = make_depth_safe(cube) cubes = make_cube_region_dict(cube) # Determine y log scale. set_y_logscale = determine_set_y_logscale(cfg, metadata) for region, cube in cubes.items(): # Make a dict of cubes for each layer. qplt.contourf(cube, 15, linewidth=0, rasterized=True) if set_y_logscale: plt.gca().set_yscale('log') if region: region_title = region else: region_title = determine_transect_str(cube, region) # Add title to plot title = ' '.join( [metadata['dataset'], metadata['long_name'], region_title]) titlify(title) # Load image format extention image_extention = diagtools.get_image_format(cfg) # Determine image filename: if multi_model: path = diagtools.folder( cfg['plot_dir']) + os.path.basename(filename).replace( '.nc', region + '_transect' + image_extention) else: path = diagtools.get_image_path( cfg, metadata, suffix=region + 'transect' + image_extention, ) # Saving files: logger.info('Saving plots to %s', path) plt.savefig(path) plt.close() provenance_record = diagtools.prepare_provenance_record( cfg, caption=f'Transect of {title}', statistics=['mean'], domain=['reg'], plot_type=['sect', 'zonal', ], ancestors=[filename], ) with ProvenanceLogger(cfg) as provenance_logger: provenance_logger.log(path, provenance_record)
[docs] def add_sea_floor(cube): """ Add a simple sea floor line from the cube mask. Parameters ---------- cube: iris.cube.Cube Input cube to use to produce the sea floor. """ land_cube = cube.copy() land_cube.data = np.ma.array(land_cube.data) mask = 1. * land_cube.data.mask if mask.shape == (): mask = np.zeros_like(land_cube.data) land_cube.data = np.ma.masked_where(mask == 0, mask) land_cube.data.mask = mask qplt.contour(land_cube, 2, cmap='Greys_r', rasterized=True)
[docs] def make_transect_contours( cfg, metadata, filename, ): """ Make a contour plot of the transect for an indivudual model. This tool loads the cube from the file, checks that the units are sensible BGC units, checks for layers, adjusts the titles accordingly, determines the ultimate file name and format, then saves the image. Parameters ---------- cfg: dict the opened global config dictionairy, passed by ESMValTool. metadata: dict The metadata dictionairy for a specific model. filename: str The preprocessed model file. """ # Load cube and set up units cube = iris.load_cube(filename) cube = diagtools.bgc_units(cube, metadata['short_name']) cube = make_depth_safe(cube) # Load threshold/thresholds. plot_details = {} colours = [] thresholds = diagtools.load_thresholds(cfg, metadata) linewidths = [1 for thres in thresholds] linestyles = ['-' for thres in thresholds] cubes = make_cube_region_dict(cube) for region, cube in cubes.items(): for itr, thres in enumerate(thresholds): colour = diagtools.get_colour_from_cmap(itr, len(thresholds)) label = str(thres) + ' ' + str(cube.units) colours.append(colour) plot_details[thres] = { 'c': colour, 'lw': 1, 'ls': '-', 'label': label } qplt.contour( cube, thresholds, colors=colours, linewidths=linewidths, linestyles=linestyles, rasterized=True) # Determine y log scale. Use gca to set scale if determine_set_y_logscale(cfg, metadata): plt.gca().set_yscale('log') add_sea_floor(cube) # Add legend diagtools.add_legend_outside_right( plot_details, plt.gca(), column_width=0.08, loc='below') # Add title to plot title = ' '.join([ metadata['dataset'], metadata['long_name'], determine_transect_str(cube, region) ]) titlify(title) # Load image format extention image_extention = diagtools.get_image_format(cfg) # Determine image filename: if metadata['dataset'].find('MultiModel') > -1: path = diagtools.folder( cfg['plot_dir']) + os.path.basename(filename) path.replace('.nc', region + '_transect_contour' + image_extention) else: path = diagtools.get_image_path( cfg, metadata, suffix=region + 'transect_contour' + image_extention, ) # Saving files: logger.info('Saving plots to %s', path) plt.savefig(path) plt.close() provenance_record = diagtools.prepare_provenance_record( cfg, caption=f'Transect of {title}', statistics=['mean'], domain=['reg'], plot_type=['sect', 'zonal', ], ancestors=[filename], ) with ProvenanceLogger(cfg) as provenance_logger: provenance_logger.log(path, provenance_record)
[docs] def multi_model_contours( cfg, metadatas, ): """ Make a multi model comparison plot showing several transect contour plots. This tool loads several cubes from the files, checks that the units are sensible BGC units, checks for layers, adjusts the titles accordingly, determines the ultimate file name and format, then saves the image. Parameters ---------- cfg: dict the opened global config dictionairy, passed by ESMValTool. metadatas: dict The metadatas dictionairy for a specific model. """ #### # Load the data for each layer as a separate cube model_cubes = {} regions = {} thresholds = {} set_y_logscale = True for filename in sorted(metadatas): cube = iris.load_cube(filename) cube = diagtools.bgc_units(cube, metadatas[filename]['short_name']) cube = make_depth_safe(cube) cubes = make_cube_region_dict(cube) model_cubes[filename] = cubes for region in model_cubes[filename]: regions[region] = True # Determine y log scale. set_y_logscale = determine_set_y_logscale(cfg, metadatas[filename]) # Load threshold/thresholds. tmp_thresholds = diagtools.load_thresholds(cfg, metadatas[filename]) for threshold in tmp_thresholds: thresholds[threshold] = True # Load image format extention image_extention = diagtools.get_image_format(cfg) # Make a plot for each layer and each threshold for region, threshold in itertools.product(regions, thresholds): logger.info('plotting threshold: \t%s', threshold) title = '' plot_details = {} # Plot each file in the group for index, filename in enumerate(sorted(metadatas)): color = diagtools.get_colour_from_cmap(index, len(metadatas)) linewidth = 1. linestyle = '-' # Determine line style for MultiModel statistics: if 'MultiModel' in metadatas[filename]['dataset']: linewidth = 2. linestyle = ':' # Determine line style for Observations if metadatas[filename]['project'] in diagtools.get_obs_projects(): color = 'black' linewidth = 1.7 linestyle = '-' qplt.contour( model_cubes[filename][region], [ threshold, ], colors=[ color, ], linewidths=linewidth, linestyles=linestyle, rasterized=True) plot_details[filename] = { 'c': color, 'ls': linestyle, 'lw': linewidth, 'label': metadatas[filename]['dataset'] } if set_y_logscale: plt.gca().set_yscale('log') title = metadatas[filename]['long_name'] units = str(model_cubes[filename][region].units) add_sea_floor(model_cubes[filename][region]) # Add title, threshold, legend to plots title = ' '.join([ title, str(threshold), units, determine_transect_str(model_cubes[filename][region], region) ]) titlify(title) plt.legend(loc='best') # Saving files: path = diagtools.get_image_path( cfg, metadatas[filename], prefix='MultipleModels', suffix='_'.join([ 'contour_tramsect', region, str(threshold) + image_extention ]), metadata_id_list=[ 'field', 'short_name', 'preprocessor', 'diagnostic', 'start_year', 'end_year' ], ) # Resize and add legend outside thew axes. plt.gcf().set_size_inches(9., 6.) diagtools.add_legend_outside_right( plot_details, plt.gca(), column_width=0.15) logger.info('Saving plots to %s', path) plt.savefig(path) plt.close() provenance_record = diagtools.prepare_provenance_record( cfg, caption=f'Transect of {title}', statistics=['mean'], domain=['reg'], plot_type=['sect', 'zonal', ], ancestors=list(metadatas.keys()), ) with ProvenanceLogger(cfg) as provenance_logger: provenance_logger.log(path, provenance_record)
[docs] def main(cfg): """ Load the config file and some metadata, then pass them the plot making tools. Parameters ---------- cfg: dict the opened global config dictionairy, passed by ESMValTool. """ ##### for index, metadata_filename in enumerate(cfg['input_files']): logger.info( 'metadata filename:\t%s', metadata_filename, ) metadatas = diagtools.get_input_files(cfg, index=index) thresholds = diagtools.load_thresholds(cfg, next(iter(metadatas.values()))) ####### # Multi model contour plots if thresholds: multi_model_contours( cfg, metadatas, ) for filename in sorted(metadatas): logger.info('-----------------') logger.info( 'model filenames:\t%s', filename, ) ###### # Time series of individual model make_transects_plots(cfg, metadatas[filename], filename) ###### # Contour maps of individual model if thresholds: make_transect_contours(cfg, metadatas[filename], filename) logger.info('Success')
if __name__ == '__main__': with run_diagnostic() as config: main(config)