From 2ed4770930737e69324703049f7953669dcd89f8 Mon Sep 17 00:00:00 2001 From: Thomas Wilder Date: Wed, 19 Nov 2025 17:43:39 +0200 Subject: [PATCH 01/11] First commits: a blank recipe but named for purpose; a copy of my_little_diagnostic.py into diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py --- .../southern_ocean/diagnostic_acc_sigma_tr.py | 137 ++++++++++++++++++ esmvaltool/recipes/recipe_so_acc_sigma_tr.yml | 0 2 files changed, 137 insertions(+) create mode 100644 esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py create mode 100644 esmvaltool/recipes/recipe_so_acc_sigma_tr.yml diff --git a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py new file mode 100644 index 0000000000..32dcaa0f91 --- /dev/null +++ b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py @@ -0,0 +1,137 @@ +""" +Look at this module for guidance how to write your own. + +Read the README_PERSONAL_DIAGNOSTIC file associated with this example; + +Module for personal diagnostics (example). +Internal imports from exmvaltool work e.g.: + +from esmvalcore.preprocessor import regrid +from esmvaltool.diag_scripts.shared.supermeans import get_supermean + +Pipe output through logger; + +Please consult the documentation for help with esmvaltool's functionalities +and best coding practices. +""" +# place your module imports here: + +# operating system manipulations (e.g. path constructions) +import os + +# to manipulate iris cubes +import iris +import matplotlib.pyplot as plt +from esmvalcore.preprocessor import area_statistics + +# import internal esmvaltool modules here +from esmvaltool.diag_scripts.shared import group_metadata, run_diagnostic + + +def _plot_time_series(cfg, cube, dataset): + """ + Example of personal diagnostic plotting function. + + Arguments: + cfg - nested dictionary of metadata + cube - the cube to plot + dataset - name of the dataset to plot + + Returns: + string; makes some time-series plots + + Note: this function is private; remove the '_' + so you can make it public. + """ + # custom local paths for e.g. plots are supported - + # here is an example + # root_dir = '/group_workspaces/jasmin2/cmip6_prep/' # edit as per need + # out_path = 'esmvaltool_users/valeriu/' # edit as per need + # local_path = os.path.join(root_dir, out_path) + # but one can use the already defined esmvaltool output paths + local_path = cfg["plot_dir"] + + # do the plotting dance + plt.plot(cube.data, label=dataset) + plt.xlabel("Time (months)") + plt.ylabel("Area average") + plt.title("Time series at (ground level - first level)") + plt.tight_layout() + plt.grid() + plt.legend() + png_name = "Time_series_" + dataset + ".png" + plt.savefig(os.path.join(local_path, png_name)) + plt.close() + + # no need to brag :) + return "I made some plots!" + + +def run_my_diagnostic(cfg): + """ + Simple example of a diagnostic. + + This is a basic (and rather esotherical) diagnostic that firstly + loads the needed model data as iris cubes, performs a difference between + values at ground level and first vertical level, then squares the + result. + + Before plotting, we grab the squared result (not all operations on cubes) + and apply an area average on it. This is a useful example of how to use + standard esmvalcore.preprocessor functionality within a diagnostic, and + especially after a certain (custom) diagnostic has been run and the user + needs to perform an operation that is already part of the preprocessor + standard library of functions. + + The user will implement their own (custom) diagnostics, but this + example shows that once the preprocessor has finished a whole lot of + user-specific metrics can be computed as part of the diagnostic, + and then plotted in various manners. + + Arguments: + cfg - nested dictionary of metadata + + Returns: + string; runs the user diagnostic + + """ + # assemble the data dictionary keyed by dataset name + # this makes use of the handy group_metadata function that + # orders the data by 'dataset'; the resulting dictionary is + # keyed on datasets e.g. dict = {'MPI-ESM-LR': [var1, var2...]} + # where var1, var2 are dicts holding all needed information per variable + my_files_dict = group_metadata(cfg["input_data"].values(), "dataset") + + # iterate over key(dataset) and values(list of vars) + for key, value in my_files_dict.items(): + # load the cube from data files only + # using a single variable here so just grab the first (and only) + # list element + cube = iris.load_cube(value[0]["filename"]) + + # the first data analysis bit: simple cube difference: + # perform a difference between ground and first levels + diff_cube = cube[:, 0, :, :] - cube[:, 1, :, :] + # square the difference'd cube just for fun + squared_cube = diff_cube**2.0 + + # the second data analysis bit (slightly more advanced): + # compute an area average over the squared cube + # to apply the area average use a preprocessor function + # rather than writing your own function + area_avg_cube = area_statistics(squared_cube, "mean") + + # finalize your analysis by plotting a time series of the + # diffed, squared and area averaged cube; call the plot function: + _plot_time_series(cfg, area_avg_cube, key) + + # that's it, we're done! + return "I am done with my first ESMValTool diagnostic!" + + +if __name__ == "__main__": + # always use run_diagnostic() to get the config (the preprocessor + # nested dictionary holding all the needed information) + with run_diagnostic() as config: + # list here the functions that need to run + run_my_diagnostic(config) diff --git a/esmvaltool/recipes/recipe_so_acc_sigma_tr.yml b/esmvaltool/recipes/recipe_so_acc_sigma_tr.yml new file mode 100644 index 0000000000..e69de29bb2 From ad5b9453fdc35fba8db0175a1482819083496769 Mon Sep 17 00:00:00 2001 From: Thomas Wilder Date: Tue, 25 Nov 2025 15:56:16 +0200 Subject: [PATCH 02/11] added some modules and defined some functions --- .../southern_ocean/diagnostic_acc_sigma_tr.py | 23 +++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py index 32dcaa0f91..8162a91187 100644 --- a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py +++ b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py @@ -15,6 +15,11 @@ and best coding practices. """ # place your module imports here: +import gsw +import xarray as xr +import numpy as np + +import cmocean # operating system manipulations (e.g. path constructions) import os @@ -27,7 +32,25 @@ # import internal esmvaltool modules here from esmvaltool.diag_scripts.shared import group_metadata, run_diagnostic +# my functions +def _plot_transect(): + + pass + +def _compute_sigma(): + pass + +def _extract_transect(): + + pass + +def main(): + pass + + + +# example functions def _plot_time_series(cfg, cube, dataset): """ Example of personal diagnostic plotting function. From 959033823915cf7a640e52f631db198c06f72fa7 Mon Sep 17 00:00:00 2001 From: Thomas Wilder Date: Thu, 27 Nov 2025 18:50:06 +0200 Subject: [PATCH 03/11] added code to def_main. --- .../southern_ocean/diagnostic_acc_sigma_tr.py | 219 ++++++++++-------- 1 file changed, 117 insertions(+), 102 deletions(-) diff --git a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py index 8162a91187..f9239bd09c 100644 --- a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py +++ b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py @@ -21,6 +21,8 @@ import cmocean +import logging + # operating system manipulations (e.g. path constructions) import os @@ -32,6 +34,10 @@ # import internal esmvaltool modules here from esmvaltool.diag_scripts.shared import group_metadata, run_diagnostic +# This part sends debug statements to stdout +logger = logging.getLogger(os.path.basename(__file__)) +logging.getLogger().addHandler(logging.StreamHandler(sys.stdout)) + # my functions def _plot_transect(): @@ -44,112 +50,121 @@ def _extract_transect(): pass -def main(): +def main(cfg): + + # assemble the data dictionary keyed by dataset name + input_data = cfg["input_data"].values() + my_files_dict = group_metadata(input_data, "dataset") + + # load xarray dataset + ds = xr.open_dataset(value[0]["filename"]) + + pass -# example functions -def _plot_time_series(cfg, cube, dataset): - """ - Example of personal diagnostic plotting function. - - Arguments: - cfg - nested dictionary of metadata - cube - the cube to plot - dataset - name of the dataset to plot - - Returns: - string; makes some time-series plots - - Note: this function is private; remove the '_' - so you can make it public. - """ - # custom local paths for e.g. plots are supported - - # here is an example - # root_dir = '/group_workspaces/jasmin2/cmip6_prep/' # edit as per need - # out_path = 'esmvaltool_users/valeriu/' # edit as per need - # local_path = os.path.join(root_dir, out_path) - # but one can use the already defined esmvaltool output paths - local_path = cfg["plot_dir"] - - # do the plotting dance - plt.plot(cube.data, label=dataset) - plt.xlabel("Time (months)") - plt.ylabel("Area average") - plt.title("Time series at (ground level - first level)") - plt.tight_layout() - plt.grid() - plt.legend() - png_name = "Time_series_" + dataset + ".png" - plt.savefig(os.path.join(local_path, png_name)) - plt.close() - - # no need to brag :) - return "I made some plots!" - - -def run_my_diagnostic(cfg): - """ - Simple example of a diagnostic. - - This is a basic (and rather esotherical) diagnostic that firstly - loads the needed model data as iris cubes, performs a difference between - values at ground level and first vertical level, then squares the - result. - - Before plotting, we grab the squared result (not all operations on cubes) - and apply an area average on it. This is a useful example of how to use - standard esmvalcore.preprocessor functionality within a diagnostic, and - especially after a certain (custom) diagnostic has been run and the user - needs to perform an operation that is already part of the preprocessor - standard library of functions. - - The user will implement their own (custom) diagnostics, but this - example shows that once the preprocessor has finished a whole lot of - user-specific metrics can be computed as part of the diagnostic, - and then plotted in various manners. - - Arguments: - cfg - nested dictionary of metadata - - Returns: - string; runs the user diagnostic - - """ - # assemble the data dictionary keyed by dataset name - # this makes use of the handy group_metadata function that - # orders the data by 'dataset'; the resulting dictionary is - # keyed on datasets e.g. dict = {'MPI-ESM-LR': [var1, var2...]} - # where var1, var2 are dicts holding all needed information per variable - my_files_dict = group_metadata(cfg["input_data"].values(), "dataset") - - # iterate over key(dataset) and values(list of vars) - for key, value in my_files_dict.items(): - # load the cube from data files only - # using a single variable here so just grab the first (and only) - # list element - cube = iris.load_cube(value[0]["filename"]) - - # the first data analysis bit: simple cube difference: - # perform a difference between ground and first levels - diff_cube = cube[:, 0, :, :] - cube[:, 1, :, :] - # square the difference'd cube just for fun - squared_cube = diff_cube**2.0 - - # the second data analysis bit (slightly more advanced): - # compute an area average over the squared cube - # to apply the area average use a preprocessor function - # rather than writing your own function - area_avg_cube = area_statistics(squared_cube, "mean") - - # finalize your analysis by plotting a time series of the - # diffed, squared and area averaged cube; call the plot function: - _plot_time_series(cfg, area_avg_cube, key) - - # that's it, we're done! - return "I am done with my first ESMValTool diagnostic!" +# # example functions +# def _plot_time_series(cfg, cube, dataset): +# """ +# Example of personal diagnostic plotting function. + +# Arguments: +# cfg - nested dictionary of metadata +# cube - the cube to plot +# dataset - name of the dataset to plot + +# Returns: +# string; makes some time-series plots + +# Note: this function is private; remove the '_' +# so you can make it public. +# """ +# # custom local paths for e.g. plots are supported - +# # here is an example +# # root_dir = '/group_workspaces/jasmin2/cmip6_prep/' # edit as per need +# # out_path = 'esmvaltool_users/valeriu/' # edit as per need +# # local_path = os.path.join(root_dir, out_path) +# # but one can use the already defined esmvaltool output paths +# local_path = cfg["plot_dir"] + +# # do the plotting dance +# plt.plot(cube.data, label=dataset) +# plt.xlabel("Time (months)") +# plt.ylabel("Area average") +# plt.title("Time series at (ground level - first level)") +# plt.tight_layout() +# plt.grid() +# plt.legend() +# png_name = "Time_series_" + dataset + ".png" +# plt.savefig(os.path.join(local_path, png_name)) +# plt.close() + +# # no need to brag :) +# return "I made some plots!" + + +# def run_my_diagnostic(cfg): +# """ +# Simple example of a diagnostic. + +# This is a basic (and rather esotherical) diagnostic that firstly +# loads the needed model data as iris cubes, performs a difference between +# values at ground level and first vertical level, then squares the +# result. + +# Before plotting, we grab the squared result (not all operations on cubes) +# and apply an area average on it. This is a useful example of how to use +# standard esmvalcore.preprocessor functionality within a diagnostic, and +# especially after a certain (custom) diagnostic has been run and the user +# needs to perform an operation that is already part of the preprocessor +# standard library of functions. + +# The user will implement their own (custom) diagnostics, but this +# example shows that once the preprocessor has finished a whole lot of +# user-specific metrics can be computed as part of the diagnostic, +# and then plotted in various manners. + +# Arguments: +# cfg - nested dictionary of metadata + +# Returns: +# string; runs the user diagnostic + +# """ +# # assemble the data dictionary keyed by dataset name +# # this makes use of the handy group_metadata function that +# # orders the data by 'dataset'; the resulting dictionary is +# # keyed on datasets e.g. dict = {'MPI-ESM-LR': [var1, var2...]} +# # where var1, var2 are dicts holding all needed information per variable +# my_files_dict = group_metadata(cfg["input_data"].values(), "dataset") + +# # iterate over key(dataset) and values(list of vars) +# for key, value in my_files_dict.items(): +# # load the cube from data files only +# # using a single variable here so just grab the first (and only) +# # list element +# cube = iris.load_cube(value[0]["filename"]) + +# # the first data analysis bit: simple cube difference: +# # perform a difference between ground and first levels +# diff_cube = cube[:, 0, :, :] - cube[:, 1, :, :] +# # square the difference'd cube just for fun +# squared_cube = diff_cube**2.0 + +# # the second data analysis bit (slightly more advanced): +# # compute an area average over the squared cube +# # to apply the area average use a preprocessor function +# # rather than writing your own function +# area_avg_cube = area_statistics(squared_cube, "mean") + +# # finalize your analysis by plotting a time series of the +# # diffed, squared and area averaged cube; call the plot function: +# _plot_time_series(cfg, area_avg_cube, key) + +# # that's it, we're done! +# return "I am done with my first ESMValTool diagnostic!" if __name__ == "__main__": @@ -157,4 +172,4 @@ def run_my_diagnostic(cfg): # nested dictionary holding all the needed information) with run_diagnostic() as config: # list here the functions that need to run - run_my_diagnostic(config) + main(config) From 51db6469c71c052581cda621868adc3f55f2b480 Mon Sep 17 00:00:00 2001 From: Thomas Wilder Date: Thu, 27 Nov 2025 18:50:55 +0200 Subject: [PATCH 04/11] a recipe with details added: description, dataset, variables. need to add preprocessors. --- esmvaltool/recipes/recipe_so_drakepassage.yml | 38 +++++++++++++++++++ 1 file changed, 38 insertions(+) create mode 100644 esmvaltool/recipes/recipe_so_drakepassage.yml diff --git a/esmvaltool/recipes/recipe_so_drakepassage.yml b/esmvaltool/recipes/recipe_so_drakepassage.yml new file mode 100644 index 0000000000..065a83d1c4 --- /dev/null +++ b/esmvaltool/recipes/recipe_so_drakepassage.yml @@ -0,0 +1,38 @@ +# ESMValTool +# recipe_so_acc_sigma_tr.yml +--- +documentation: + description: + Creates a transect through Drake Passage of zonal velocity with + the sigma contour. + title: + Drake Passage zonal velocity transect with sigma contour + + authors: + - wilder_thomas + +datasets: + - dataset: HadGEM3-GC31-LL + +diagnostics: + # ------------------------------------------------------------------------ + # Drake Passage zonal velocity transect with sigma contour + # -------------------------------------------------------------------------------- + drake_passage_transect: + description: Transect of zonal velocity and sigma + variables: + thetao: # Temperature ocean + ensemble: r1i1p1f3 + exp: historical + grid: gn + mip: Omon + project: CMIP6 + timerange: 2000/2001 + so: # Salinity ocean + ensemble: r1i1p1f3 + exp: historical + grid: gn + mip: Omon + project: CMIP6 + timerange: 2000/2001 + scripts: From 3c721cd108320535aa0adc1ac66db9f2ab93d04c Mon Sep 17 00:00:00 2001 From: Thomas Wilder Date: Fri, 28 Nov 2025 18:42:22 +0200 Subject: [PATCH 05/11] - renamed the recipe file. - current recipe working with basic python diagnostic. --- esmvaltool/recipes/recipe_so_acc_sigma_tr.yml | 0 esmvaltool/recipes/recipe_so_drakepassage.yml | 32 +++++++++++++------ 2 files changed, 23 insertions(+), 9 deletions(-) delete mode 100644 esmvaltool/recipes/recipe_so_acc_sigma_tr.yml diff --git a/esmvaltool/recipes/recipe_so_acc_sigma_tr.yml b/esmvaltool/recipes/recipe_so_acc_sigma_tr.yml deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/esmvaltool/recipes/recipe_so_drakepassage.yml b/esmvaltool/recipes/recipe_so_drakepassage.yml index 065a83d1c4..886f81ac5e 100644 --- a/esmvaltool/recipes/recipe_so_drakepassage.yml +++ b/esmvaltool/recipes/recipe_so_drakepassage.yml @@ -9,16 +9,25 @@ documentation: Drake Passage zonal velocity transect with sigma contour authors: - - wilder_thomas + - demora_lee datasets: - dataset: HadGEM3-GC31-LL +preprocessors: + # ------------------------------------------------------------------------ + # Preprocessor settings for Drake Passage transect + # ------------------------------------------------------------------------ + time_mean: + climate_statistics: + operator: 'mean' + period: 'full' + diagnostics: # ------------------------------------------------------------------------ # Drake Passage zonal velocity transect with sigma contour # -------------------------------------------------------------------------------- - drake_passage_transect: + transect: description: Transect of zonal velocity and sigma variables: thetao: # Temperature ocean @@ -28,11 +37,16 @@ diagnostics: mip: Omon project: CMIP6 timerange: 2000/2001 - so: # Salinity ocean - ensemble: r1i1p1f3 - exp: historical - grid: gn - mip: Omon - project: CMIP6 - timerange: 2000/2001 + preprocessor: time_mean + # so: # Salinity ocean + # ensemble: r1i1p1f3 + # exp: historical + # grid: gn + # mip: Omon + # project: CMIP6 + # timerange: 2000/2001 + # preprocessor: time_mean + scripts: + script1: + script: /scratch/project_465001823/twilder/ESMValTool/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py From 777a843e36cd8dc8c8dfd0a530e2865eb8a847c9 Mon Sep 17 00:00:00 2001 From: Thomas Wilder Date: Fri, 28 Nov 2025 18:44:07 +0200 Subject: [PATCH 06/11] - removed some imports and added a log message to indicate it runs without error. --- .../southern_ocean/diagnostic_acc_sigma_tr.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py index f9239bd09c..149a70b96e 100644 --- a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py +++ b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py @@ -15,16 +15,17 @@ and best coding practices. """ # place your module imports here: -import gsw +# import gsw import xarray as xr import numpy as np -import cmocean +# import cmocean import logging # operating system manipulations (e.g. path constructions) import os +import sys # to manipulate iris cubes import iris @@ -52,12 +53,16 @@ def _extract_transect(): def main(cfg): - # assemble the data dictionary keyed by dataset name - input_data = cfg["input_data"].values() - my_files_dict = group_metadata(input_data, "dataset") + # # assemble the data dictionary keyed by dataset name + # input_data = cfg["input_data"].values() + # my_files_dict = group_metadata(input_data, "dataset") - # load xarray dataset - ds = xr.open_dataset(value[0]["filename"]) + # logger.info('my_files_dict: %s', my_files_dict) + + logger.info('we have lift off!') + + # # load xarray dataset + # ds = xr.open_dataset(value[0]["filename"]) pass From 6c7ad5568a3d8869474f706e355436929517f6f7 Mon Sep 17 00:00:00 2001 From: thomaswilder Date: Tue, 31 Mar 2026 21:12:45 +0100 Subject: [PATCH 07/11] dev: updated recipe with correct preprocessor, and begun laying out diagnostic code in a class oop format. --- .../southern_ocean/diagnostic_acc_sigma_tr.py | 24 ++++++++-- esmvaltool/recipes/recipe_so_drakepassage.yml | 46 ++++++++++++++----- 2 files changed, 55 insertions(+), 15 deletions(-) diff --git a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py index 149a70b96e..d0e4dc5d4b 100644 --- a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py +++ b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py @@ -33,12 +33,30 @@ from esmvalcore.preprocessor import area_statistics # import internal esmvaltool modules here -from esmvaltool.diag_scripts.shared import group_metadata, run_diagnostic +from esmvaltool.diag_scripts.shared import group_metadata, run_diagnostic, ProvenanceLogger + +# reuse tools already developed for ocean diagnostics +# from esmvaltool.diag_scripts.ocean import diagnostic_tools as diagtools # This part sends debug statements to stdout logger = logging.getLogger(os.path.basename(__file__)) logging.getLogger().addHandler(logging.StreamHandler(sys.stdout)) + +class TransectDiagnostic: + ''' + Class used to compute a transect of a velocity with + a sigma2 contour overlaid. + ''' + def __init__(self, config): + self.config = config + + def compute_sigma(self, ds): + ''' + Compute the sigma2 contour from the input dataset. + ''' + + # my functions def _plot_transect(): @@ -54,8 +72,8 @@ def _extract_transect(): def main(cfg): # # assemble the data dictionary keyed by dataset name - # input_data = cfg["input_data"].values() - # my_files_dict = group_metadata(input_data, "dataset") + input_data = cfg["input_data"] + grouped = group_metadata(input_data, "dataset") # logger.info('my_files_dict: %s', my_files_dict) diff --git a/esmvaltool/recipes/recipe_so_drakepassage.yml b/esmvaltool/recipes/recipe_so_drakepassage.yml index 886f81ac5e..37a773a8c4 100644 --- a/esmvaltool/recipes/recipe_so_drakepassage.yml +++ b/esmvaltool/recipes/recipe_so_drakepassage.yml @@ -4,7 +4,7 @@ documentation: description: Creates a transect through Drake Passage of zonal velocity with - the sigma contour. + sigma contours overlaid. title: Drake Passage zonal velocity transect with sigma contour @@ -18,10 +18,20 @@ preprocessors: # ------------------------------------------------------------------------ # Preprocessor settings for Drake Passage transect # ------------------------------------------------------------------------ - time_mean: + # time_mean: + # climate_statistics: + # operator: 'mean' + # period: 'full' + + drakepassage_transect: # does this order matter? No. climate_statistics: operator: 'mean' period: 'full' + extract_region: # probably needs changing + start_longitude: 293 + end_longitude: 300 + start_latitude: -65 + end_latitude: -54 diagnostics: # ------------------------------------------------------------------------ @@ -29,6 +39,10 @@ diagnostics: # -------------------------------------------------------------------------------- transect: description: Transect of zonal velocity and sigma + themes: + - phys + realms: + - ocean variables: thetao: # Temperature ocean ensemble: r1i1p1f3 @@ -37,16 +51,24 @@ diagnostics: mip: Omon project: CMIP6 timerange: 2000/2001 - preprocessor: time_mean - # so: # Salinity ocean - # ensemble: r1i1p1f3 - # exp: historical - # grid: gn - # mip: Omon - # project: CMIP6 - # timerange: 2000/2001 - # preprocessor: time_mean + preprocessor: drakepassage_transect + so: # Salinity ocean + ensemble: r1i1p1f3 + exp: historical + grid: gn + mip: Omon + project: CMIP6 + timerange: 2000/2001 + preprocessor: drakepassage_transect + uo: # Zonal velocity ocean + ensemble: r1i1p1f3 + exp: historical + grid: gn + mip: Omon + project: CMIP6 + timerange: 2000/2001 + preprocessor: drakepassage_transect scripts: script1: - script: /scratch/project_465001823/twilder/ESMValTool/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py + script: southern_ocean/diagnostic_acc_sigma_tr.py From 32a73e4c7b794c85eed6188d4b43dfde80cd762e Mon Sep 17 00:00:00 2001 From: thomaswilder Date: Tue, 31 Mar 2026 21:13:17 +0100 Subject: [PATCH 08/11] comment --- .../diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py index d0e4dc5d4b..c1986123a9 100644 --- a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py +++ b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py @@ -42,6 +42,8 @@ logger = logging.getLogger(os.path.basename(__file__)) logging.getLogger().addHandler(logging.StreamHandler(sys.stdout)) +# uses a class style following the example of eady_growth_rate.py. + class TransectDiagnostic: ''' From 894b97580f8cbde24cc0de06784ceebf4dacfb17 Mon Sep 17 00:00:00 2001 From: thomaswilder Date: Wed, 1 Apr 2026 21:38:17 +0100 Subject: [PATCH 09/11] dev: setting up the class structure, exploring input_data dict, and writing logs to log.txt. --- .../southern_ocean/diagnostic_acc_sigma_tr.py | 170 ++++-------------- esmvaltool/recipes/recipe_so_drakepassage.yml | 16 +- 2 files changed, 45 insertions(+), 141 deletions(-) diff --git a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py index c1986123a9..f20da8c11a 100644 --- a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py +++ b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py @@ -23,9 +23,11 @@ import logging + # operating system manipulations (e.g. path constructions) import os import sys +from pathlib import Path # to manipulate iris cubes import iris @@ -39,13 +41,16 @@ # from esmvaltool.diag_scripts.ocean import diagnostic_tools as diagtools # This part sends debug statements to stdout -logger = logging.getLogger(os.path.basename(__file__)) -logging.getLogger().addHandler(logging.StreamHandler(sys.stdout)) +# logger = logging.getLogger(os.path.basename(__file__)) +# logging.getLogger().addHandler(logging.StreamHandler(sys.stdout)) + +# This outputs to the run/transect/script1/log.txt file +logger = logging.getLogger(Path(__file__).stem) # uses a class style following the example of eady_growth_rate.py. -class TransectDiagnostic: +class TransectDiagnostic(): ''' Class used to compute a transect of a velocity with a sigma2 contour overlaid. @@ -53,148 +58,47 @@ class TransectDiagnostic: def __init__(self, config): self.config = config - def compute_sigma(self, ds): - ''' - Compute the sigma2 contour from the input dataset. - ''' - - -# my functions -def _plot_transect(): - - pass - -def _compute_sigma(): - pass - -def _extract_transect(): - - pass - -def main(cfg): - - # # assemble the data dictionary keyed by dataset name - input_data = cfg["input_data"] - grouped = group_metadata(input_data, "dataset") - - # logger.info('my_files_dict: %s', my_files_dict) - - logger.info('we have lift off!') - - # # load xarray dataset - # ds = xr.open_dataset(value[0]["filename"]) + def _get_data(self): + input_data = self.config["input_data"] - pass + logger.info("Input data: %s", input_data) + filenames = input_data.keys() + logger.info("Filenames: %s", filenames) + + # groups = group_metadata(input_data, "variable_group", sort="dataset") + # for group_name in groups: + # logger.info("Processing variable %s", group_name) + # for attributes in groups[group_name]: + # logger.info("Processing dataset %s", attributes["dataset"]) + # input_file = attributes["filename"] + # logger.info("Loading data from %s", input_file) + return input_data -# # example functions -# def _plot_time_series(cfg, cube, dataset): -# """ -# Example of personal diagnostic plotting function. - -# Arguments: -# cfg - nested dictionary of metadata -# cube - the cube to plot -# dataset - name of the dataset to plot - -# Returns: -# string; makes some time-series plots - -# Note: this function is private; remove the '_' -# so you can make it public. -# """ -# # custom local paths for e.g. plots are supported - -# # here is an example -# # root_dir = '/group_workspaces/jasmin2/cmip6_prep/' # edit as per need -# # out_path = 'esmvaltool_users/valeriu/' # edit as per need -# # local_path = os.path.join(root_dir, out_path) -# # but one can use the already defined esmvaltool output paths -# local_path = cfg["plot_dir"] - -# # do the plotting dance -# plt.plot(cube.data, label=dataset) -# plt.xlabel("Time (months)") -# plt.ylabel("Area average") -# plt.title("Time series at (ground level - first level)") -# plt.tight_layout() -# plt.grid() -# plt.legend() -# png_name = "Time_series_" + dataset + ".png" -# plt.savefig(os.path.join(local_path, png_name)) -# plt.close() - -# # no need to brag :) -# return "I made some plots!" - - -# def run_my_diagnostic(cfg): -# """ -# Simple example of a diagnostic. - -# This is a basic (and rather esotherical) diagnostic that firstly -# loads the needed model data as iris cubes, performs a difference between -# values at ground level and first vertical level, then squares the -# result. - -# Before plotting, we grab the squared result (not all operations on cubes) -# and apply an area average on it. This is a useful example of how to use -# standard esmvalcore.preprocessor functionality within a diagnostic, and -# especially after a certain (custom) diagnostic has been run and the user -# needs to perform an operation that is already part of the preprocessor -# standard library of functions. - -# The user will implement their own (custom) diagnostics, but this -# example shows that once the preprocessor has finished a whole lot of -# user-specific metrics can be computed as part of the diagnostic, -# and then plotted in various manners. + def _compute_sigma(self, ds): + ''' + Compute the sigma2 contour from the input dataset. + ''' -# Arguments: -# cfg - nested dictionary of metadata + pass -# Returns: -# string; runs the user diagnostic + def _plot_transect(self): -# """ -# # assemble the data dictionary keyed by dataset name -# # this makes use of the handy group_metadata function that -# # orders the data by 'dataset'; the resulting dictionary is -# # keyed on datasets e.g. dict = {'MPI-ESM-LR': [var1, var2...]} -# # where var1, var2 are dicts holding all needed information per variable -# my_files_dict = group_metadata(cfg["input_data"].values(), "dataset") + pass -# # iterate over key(dataset) and values(list of vars) -# for key, value in my_files_dict.items(): -# # load the cube from data files only -# # using a single variable here so just grab the first (and only) -# # list element -# cube = iris.load_cube(value[0]["filename"]) + def call(self): + input_data = self._get_data() + return input_data -# # the first data analysis bit: simple cube difference: -# # perform a difference between ground and first levels -# diff_cube = cube[:, 0, :, :] - cube[:, 1, :, :] -# # square the difference'd cube just for fun -# squared_cube = diff_cube**2.0 +def main(): + """Run Eady Growth Rate diagnostic.""" + with run_diagnostic() as config: + diagnostic = TransectDiagnostic(config).call() -# # the second data analysis bit (slightly more advanced): -# # compute an area average over the squared cube -# # to apply the area average use a preprocessor function -# # rather than writing your own function -# area_avg_cube = area_statistics(squared_cube, "mean") - -# # finalize your analysis by plotting a time series of the -# # diffed, squared and area averaged cube; call the plot function: -# _plot_time_series(cfg, area_avg_cube, key) - -# # that's it, we're done! -# return "I am done with my first ESMValTool diagnostic!" if __name__ == "__main__": - # always use run_diagnostic() to get the config (the preprocessor - # nested dictionary holding all the needed information) - with run_diagnostic() as config: - # list here the functions that need to run - main(config) + main() diff --git a/esmvaltool/recipes/recipe_so_drakepassage.yml b/esmvaltool/recipes/recipe_so_drakepassage.yml index 37a773a8c4..b8b0aaa9ff 100644 --- a/esmvaltool/recipes/recipe_so_drakepassage.yml +++ b/esmvaltool/recipes/recipe_so_drakepassage.yml @@ -60,14 +60,14 @@ diagnostics: project: CMIP6 timerange: 2000/2001 preprocessor: drakepassage_transect - uo: # Zonal velocity ocean - ensemble: r1i1p1f3 - exp: historical - grid: gn - mip: Omon - project: CMIP6 - timerange: 2000/2001 - preprocessor: drakepassage_transect + # uo: # Zonal velocity ocean + # ensemble: r1i1p1f3 + # exp: historical + # grid: gn + # mip: Omon + # project: CMIP6 + # timerange: 2000/2001 + # preprocessor: drakepassage_transect scripts: script1: From 84f9bd73a76d7cd4dcacfb75c3c80d927c083a55 Mon Sep 17 00:00:00 2001 From: thomaswilder Date: Mon, 1 Jun 2026 15:52:25 +0100 Subject: [PATCH 10/11] refactor,dev: functional approach, loading in single dataset, computing density, working end-to-end. --- .../diag_scripts/southern_ocean/__init__.py | 0 .../southern_ocean/diagnostic_acc_sigma_tr.py | 110 ++++++++++++------ esmvaltool/recipes/recipe_so_drakepassage.yml | 16 +-- 3 files changed, 78 insertions(+), 48 deletions(-) create mode 100644 esmvaltool/diag_scripts/southern_ocean/__init__.py diff --git a/esmvaltool/diag_scripts/southern_ocean/__init__.py b/esmvaltool/diag_scripts/southern_ocean/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py index f20da8c11a..a71ec47f8d 100644 --- a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py +++ b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py @@ -15,7 +15,7 @@ and best coding practices. """ # place your module imports here: -# import gsw +import gsw import xarray as xr import numpy as np @@ -23,7 +23,6 @@ import logging - # operating system manipulations (e.g. path constructions) import os import sys @@ -38,7 +37,7 @@ from esmvaltool.diag_scripts.shared import group_metadata, run_diagnostic, ProvenanceLogger # reuse tools already developed for ocean diagnostics -# from esmvaltool.diag_scripts.ocean import diagnostic_tools as diagtools +from esmvaltool.diag_scripts.ocean import diagnostic_tools as diagtools # This part sends debug statements to stdout # logger = logging.getLogger(os.path.basename(__file__)) @@ -47,58 +46,93 @@ # This outputs to the run/transect/script1/log.txt file logger = logging.getLogger(Path(__file__).stem) -# uses a class style following the example of eady_growth_rate.py. +# module checks +logger.info("GSW version: %s", gsw.__version__) + +def _get_data(cfg): -class TransectDiagnostic(): ''' - Class used to compute a transect of a velocity with - a sigma2 contour overlaid. + Get the data from the input files specified in the recipe. + Open the files and return a dictionary of xarray datasets. ''' - def __init__(self, config): - self.config = config - def _get_data(self): + #TODO needs adjusting for multi dataset recipes - input_data = self.config["input_data"] + input_data = cfg["input_data"] - logger.info("Input data: %s", input_data) + ds = {} + for filename, attributes in input_data.items(): + logger.info("Loading %s", filename) + variable_name = attributes["short_name"] + ds[variable_name] = xr.open_dataset(filename) - filenames = input_data.keys() - logger.info("Filenames: %s", filenames) - - # groups = group_metadata(input_data, "variable_group", sort="dataset") - # for group_name in groups: - # logger.info("Processing variable %s", group_name) - # for attributes in groups[group_name]: - # logger.info("Processing dataset %s", attributes["dataset"]) - # input_file = attributes["filename"] - # logger.info("Loading data from %s", input_file) + return ds +def _compute_sigma(ds): + ''' + Compute the sigma2 variable from the input dataset. + Uses gsw.density.sigma2 from the gsw toolbox + https://teos-10.github.io/GSW-Python/_modules/gsw/_wrapped_ufuncs.html#sigma2 + ''' - return input_data + # logger.info("Numpy so is: %s", ds["so"]["so"].to_numpy()) + # logger.info("Numpy thetao is: %s", ds["thetao"]["thetao"].to_numpy()) - def _compute_sigma(self, ds): - ''' - Compute the sigma2 contour from the input dataset. - ''' + #TODO need to convert to absolute salinity and conservative temperature + p = gsw.p_from_z(-ds["lev"].to_numpy(), np.mean(ds["lat"].to_numpy())) + logger.info("Pressure levels: %s", p) - pass + CT = gsw.CT_from_pt(ds["so"]["so"].to_numpy(), ds["thetao"]["thetao"].to_numpy()) + logger.info("Conservative temperature: %s", CT) - def _plot_transect(self): + SA = gsw.SA_from_SP(ds["so"]["so"].to_numpy(), + p, + ds["lon"].to_numpy(), + ds["lat"].to_numpy()) - pass + # compute sigma2 + sigma2 = gsw.sigma2(SA, CT) - def call(self): - input_data = self._get_data() - return input_data + # create a new xarray dataset with the same coordinates as the input datasets + sigma2_ds = xr.Dataset( + { + "sigma2": (["lev", "lat"], sigma2) + }, + coords={ + "time": ds["so"]["time"], + "lev": ds["so"]["lev"], + "lat": ds["so"]["lat"], + "lon": ds["so"]["lon"] + }, + attrs={ + "short_name": "sigma2", + "long_name": "Potential density anomaly referenced to 2000 dbar", + "units": "kg/m^3", + } + ) + + return sigma2_ds + +def _plot_transect(ds): + + pass -def main(): - """Run Eady Growth Rate diagnostic.""" - with run_diagnostic() as config: - diagnostic = TransectDiagnostic(config).call() + +def main(cfg): + # get the data + datasets = _get_data(cfg) + + logger.info("Datasets: %s", datasets) + + sigma2_ds = _compute_sigma(datasets) + + logger.info("Sigma2 dataset: %s", sigma2_ds) + + # _plot_transect(sigma2_ds) if __name__ == "__main__": - main() + with run_diagnostic() as config: + main(config) diff --git a/esmvaltool/recipes/recipe_so_drakepassage.yml b/esmvaltool/recipes/recipe_so_drakepassage.yml index b8b0aaa9ff..a1493f5559 100644 --- a/esmvaltool/recipes/recipe_so_drakepassage.yml +++ b/esmvaltool/recipes/recipe_so_drakepassage.yml @@ -18,20 +18,16 @@ preprocessors: # ------------------------------------------------------------------------ # Preprocessor settings for Drake Passage transect # ------------------------------------------------------------------------ - # time_mean: - # climate_statistics: - # operator: 'mean' - # period: 'full' - drakepassage_transect: # does this order matter? No. climate_statistics: operator: 'mean' period: 'full' - extract_region: # probably needs changing - start_longitude: 293 - end_longitude: 300 - start_latitude: -65 - end_latitude: -54 + regrid: + target_grid: 1x1 + scheme: area_weighted + extract_transect: # probably needs changing + longitude: 293.0 + latitude: [-65., -54.] diagnostics: # ------------------------------------------------------------------------ From 2c51ee984714a3bb40064a8d97be6035563474fe Mon Sep 17 00:00:00 2001 From: thomaswilder Date: Mon, 1 Jun 2026 16:23:39 +0100 Subject: [PATCH 11/11] dev: sigma2 correctly calculated, ready to implement plot. --- .../southern_ocean/diagnostic_acc_sigma_tr.py | 36 +++++++++++++------ 1 file changed, 25 insertions(+), 11 deletions(-) diff --git a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py index a71ec47f8d..cde30487a4 100644 --- a/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py +++ b/esmvaltool/diag_scripts/southern_ocean/diagnostic_acc_sigma_tr.py @@ -79,17 +79,30 @@ def _compute_sigma(ds): # logger.info("Numpy so is: %s", ds["so"]["so"].to_numpy()) # logger.info("Numpy thetao is: %s", ds["thetao"]["thetao"].to_numpy()) - #TODO need to convert to absolute salinity and conservative temperature - p = gsw.p_from_z(-ds["lev"].to_numpy(), np.mean(ds["lat"].to_numpy())) - logger.info("Pressure levels: %s", p) - - CT = gsw.CT_from_pt(ds["so"]["so"].to_numpy(), ds["thetao"]["thetao"].to_numpy()) - logger.info("Conservative temperature: %s", CT) - - SA = gsw.SA_from_SP(ds["so"]["so"].to_numpy(), - p, - ds["lon"].to_numpy(), - ds["lat"].to_numpy()) + # pull the data and coords as numpy arrays + so = ds["so"]["so"].to_numpy() # (lev, lat) + thetao = ds["thetao"]["thetao"].to_numpy() # (lev, lat) + depth = ds["so"]["lev"].to_numpy() # (lev,) + lat = ds["so"]["lat"].to_numpy() # (lat,) + lon = float(ds["so"]["lon"]) # scalar transect longitude + + # broadcast depth and lat onto the full (lev, lat) transect grid so all + # operands match the 2D salinity/temperature fields (gsw broadcasts + # every argument together element-wise). + depth2d, lat2d = np.meshgrid(depth, lat, indexing="ij") # both (lev, lat) + + # calculate pressure from depth (latitude-dependent) + p = gsw.p_from_z(-depth2d, lat2d) + logger.info("Pressure shape: %s", p.shape) + + # TEOS-10 chain: SA (absolute salinity) -> CT (conservative temp) -> sigma2. + # calculate absolute salinity from practical salinity + SA = gsw.SA_from_SP(so, p, lon, lat2d) + logger.info("Absolute salinity shape: %s", SA.shape) + + # calculate conservative temperature from potential temperature + CT = gsw.CT_from_pt(SA, thetao) + logger.info("Conservative temperature shape: %s", CT.shape) # compute sigma2 sigma2 = gsw.sigma2(SA, CT) @@ -128,6 +141,7 @@ def main(cfg): sigma2_ds = _compute_sigma(datasets) logger.info("Sigma2 dataset: %s", sigma2_ds) + # _plot_transect(sigma2_ds)