From b3167f82856a7afa58be10a3e5c57928424978c5 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Mon, 13 Oct 2025 16:58:43 +0100 Subject: [PATCH 01/26] In progress: - Removed obs from Roach-style graphic (left description for reference). - Fed in model datasets from recipe differently to allow for obs having a different variable. To do: - Calculate and plot the obs. - Print the obs trend values. --- .../diag_scripts/seaice/seaice_sensitivity.py | 130 +++++++----------- .../recipes/recipe_seaice_sensitivity.yml | 107 +++++++------- 2 files changed, 106 insertions(+), 131 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 9f1c0c41d3..6b3b1053e6 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -6,6 +6,7 @@ import iris import matplotlib.pyplot as plt import pandas as pd +import yaml from matplotlib.colors import Normalize from scipy import stats @@ -37,13 +38,24 @@ def get_provenance_record(cfg, caption): return record -def list_datasets(data): - """Actually returns a set of datatsets, to avoid duplication.""" - logger.debug("listing datasets") - datasets = set() +def create_categorised_dataset_dict(data): + """Create a dictionary of datasets categorised into models and observations.""" + logger.debug("Creating categorised dataset dictionary from %s", data) + # Initialise blank dictionary + dict = { + 'models': {}, + 'observations': {}, + } + + # Iterate over the data for element in data: - datasets.add(element["dataset"]) - return datasets + if 'observation' in element: + if element['observation']: + dict['observations'][element['dataset']] = {} + else: + dict['models'][element['dataset']] = {} + + return dict def extract_cube(data, variable_group): @@ -142,28 +154,6 @@ def write_obs_from_cfg(cfg): obs_dict["notz_style"]["std_dev"] = notz_values["standard deviation"] obs_dict["notz_style"]["plausible"] = notz_values["plausible range"] - # Add a blank dictionary for the Roach-style plot - obs_dict["roach_style"] = {} - - # Add each observation point to the dictionary - roach_values = cfg["observations"]["annual trends (Roach-style plot)"] - for point in roach_values.keys(): - obs_dict["roach_style"][point] = {} - - # Add the individual values for the observation point - obs_dict["roach_style"][point]["annual_tas_trend"] = roach_values[ - point - ]["GMST trend"] - obs_dict["roach_style"][point]["annual_siconc_trend"] = roach_values[ - point - ]["SIA trend"] - obs_dict["roach_style"][point]["r_value"] = roach_values[point][ - "Pearson CC of SIA over GMST" - ] - obs_dict["roach_style"][point]["p_value"] = roach_values[point][ - "significance of SIA over GMST" - ] - return obs_dict @@ -208,7 +198,7 @@ def write_dictionary_to_csv(cfg, model_dict, filename): csv_filepath = f"{cfg['work_dir']}/{filename}.csv" # Write the data to a csv file (via a Pandas DataFrame) - pd.DataFrame.from_dict(model_dict, orient="index").to_csv(csv_filepath) + pd.DataFrame.from_dict(model_dict, orient="index").to_csv(csv_filepath, index_label="Dataset") logger.info("Wrote data to %s", csv_filepath) # Create a provenance record for the csv file @@ -349,54 +339,36 @@ def roach_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): if inner_dict["label"] == "to_label": plt.annotate(dataset, xy=(x, y), xytext=(x + 0.01, y - 0.005)) - # Read from observations dictionary - obs_years = titles_dictionary["obs"]["obs_period"] - obs_dict = titles_dictionary["obs"]["roach_style"] - - # Add the observations - for point in obs_dict.keys(): - # Get the values for the point - x = obs_dict[point]["annual_tas_trend"] - y = obs_dict[point]["annual_siconc_trend"] - r_corr = obs_dict[point]["r_value"] - p_val = obs_dict[point]["p_value"] - - # Provide a default colour for the point if Pearson coefficient is missing - if r_corr is None: - r_corr = 0 - - # Provide a pattern for the point if the p-value is present and sufficiently large - if p_val is not None and p_val >= 0.05: - h = 5 * "/" # This is a hatch pattern - else: - h = None - - # Plot the point only if both x and y values are provided - if x is not None and y is not None: - plt.scatter( - x, - y, - marker="s", - s=150, - c=[r_corr], - hatch=h, - cmap=cmap, - norm=norm, - zorder=0, - edgecolors="black", - ) + # # Provide a default colour for the point if Pearson coefficient is missing + # if r_corr is None: + # r_corr = 0 + # + # # Provide a pattern for the point if the p-value is present and sufficiently large + # if p_val is not None and p_val >= 0.05: + # h = 5 * "/" # This is a hatch pattern + # else: + # h = None + # + # # Plot the point only if both x and y values are provided + # if x is not None and y is not None: + # plt.scatter( + # x, + # y, + # marker="s", + # s=150, + # c=[r_corr], + # hatch=h, + # cmap=cmap, + # norm=norm, + # zorder=0, + # edgecolors="black", + # ) # Add a colour bar plt.colorbar(label="Pearson correlation coefficient") - # Create caption based on whether observational temp trend is present - if obs_dict["first point"]["annual_tas_trend"] is not None: - caption = ( - "Decadal trends of sea ice area and global mean temperature." - f"Observations from {obs_years} are plotted as squares." - ) - else: - caption = "Decadal trends of sea ice area and global mean temperature." + # Create caption + caption = "Decadal trends of sea ice area and global mean temperature." # Save the figure (also closes it) save_figure( @@ -419,22 +391,18 @@ def main(cfg): titles_and_obs_dict = create_titles_dict(input_data, cfg) logger.debug("Titles and observations dictionary: %s", titles_and_obs_dict) - # Initialize blank data dictionary to send to plotting codes later - data_dict = {} - # Get list of datasets from cfg logger.info("Listing datasets in the data") - datasets = list_datasets(input_data) + dataset_dict = create_categorised_dataset_dict(input_data) + data_dict = dataset_dict['models'] + # Iterate over each dataset - for dataset in datasets: + for dataset in data_dict.keys(): # Select only data from that dataset logger.debug("Selecting data from %s", dataset) selection = select_metadata(input_data, dataset=dataset) - # Add the dataset to the dictionary with a blank inner dictionary - data_dict[dataset] = {} - # Add an entry to determine labelling in plots if "label_dataset" in selection[0]: data_dict[dataset]["label"] = "to_label" diff --git a/esmvaltool/recipes/recipe_seaice_sensitivity.yml b/esmvaltool/recipes/recipe_seaice_sensitivity.yml index b64191b41c..1a232805b0 100644 --- a/esmvaltool/recipes/recipe_seaice_sensitivity.yml +++ b/esmvaltool/recipes/recipe_seaice_sensitivity.yml @@ -21,24 +21,33 @@ documentation: maintainer: - parsons_naomi -defaults: &defaults {ensemble: r1i1p1f1, exp: historical, grid: gn, project: CMIP6} -datasets: - - {<<: *defaults, dataset: HadGEM3-GC31-LL, institute: MOHC, ensemble: r1i1p1f3, label_dataset: True} - - {<<: *defaults, dataset: UKESM1-0-LL, institute: MOHC, ensemble: r1i1p1f2, label_dataset: True} - - {<<: *defaults, dataset: ACCESS-CM2, institute: CSIRO-ARCCSS} - - {<<: *defaults, dataset: ACCESS-ESM1-5, institute: CSIRO} - - {<<: *defaults, dataset: BCC-CSM2-MR, institute: BCC} - - {<<: *defaults, dataset: CAMS-CSM1-0, institute: CAMS} - - {<<: *defaults, dataset: CanESM5, institute: CCCma} - - {<<: *defaults, dataset: CESM2, institute: NCAR} - - {<<: *defaults, dataset: CESM2-WACCM, institute: NCAR} - - {<<: *defaults, dataset: CESM2-WACCM-FV2, institute: NCAR} - - {<<: *defaults, dataset: FIO-ESM-2-0, institute: FIO-QLNM} - - {<<: *defaults, dataset: MIROC6, institute: MIROC} - - {<<: *defaults, dataset: MPI-ESM-1-2-HAM, institute: HAMMOZ-Consortium} - - {<<: *defaults, dataset: MPI-ESM1-2-HR, institute: MPI-M} - - {<<: *defaults, dataset: MPI-ESM1-2-LR, institute: MPI-M} - - {<<: *defaults, dataset: MRI-ESM2-0, institute: MRI} +model_defaults: &model_defaults { ensemble: r1i1p1f1, exp: historical, grid: gn, project: CMIP6} +model_datasets: &model_datasets + - {<<: *model_defaults, dataset: HadGEM3-GC31-LL, institute: MOHC, ensemble: r1i1p1f3, label_dataset: True} + - {<<: *model_defaults, dataset: UKESM1-0-LL, institute: MOHC, ensemble: r1i1p1f2, label_dataset: True} + - {<<: *model_defaults, dataset: ACCESS-CM2, institute: CSIRO-ARCCSS} + - {<<: *model_defaults, dataset: ACCESS-ESM1-5, institute: CSIRO} + - {<<: *model_defaults, dataset: BCC-CSM2-MR, institute: BCC} + - {<<: *model_defaults, dataset: CAMS-CSM1-0, institute: CAMS} + - {<<: *model_defaults, dataset: CanESM5, institute: CCCma} + - {<<: *model_defaults, dataset: CESM2, institute: NCAR} + - {<<: *model_defaults, dataset: CESM2-WACCM, institute: NCAR} + - {<<: *model_defaults, dataset: CESM2-WACCM-FV2, institute: NCAR} + - {<<: *model_defaults, dataset: FIO-ESM-2-0, institute: FIO-QLNM} + - {<<: *model_defaults, dataset: MIROC6, institute: MIROC} + - {<<: *model_defaults, dataset: MPI-ESM-1-2-HAM, institute: HAMMOZ-Consortium} + - {<<: *model_defaults, dataset: MPI-ESM1-2-HR, institute: MPI-M} + - {<<: *model_defaults, dataset: MPI-ESM1-2-LR, institute: MPI-M} + - {<<: *model_defaults, dataset: MRI-ESM2-0, institute: MRI} + +obs_defaults: &obs_defaults { project: OBS, observation: True, tier: 2} +tasa_obs: &tasa_obs + - { <<: *obs_defaults, dataset: GISTEMP, type: ground, version: v4 } + - { <<: *obs_defaults, dataset: HadCRUT5, type: ground, version: 5.0.1.0-analysis } + - { <<: *obs_defaults, dataset: NOAAGlobalTemp, type: ground, version: v5.0.0, mip: Amon } + +arctic_si_obs: &arctic_si_obs + - { <<: *obs_defaults, dataset: OSI-450-nh, type: reanaly, version: v3, mip: OImon, supplementary_variables: [{short_name: areacello, mip: fx}]} preprocessors: extract_test_period: @@ -112,9 +121,24 @@ diagnostics: siconc: preprocessor: pp_arctic_sept_sea_ice mip: SImon + additional_datasets: + *model_datasets tas: preprocessor: pp_avg_ann_global_temp mip: Amon + additional_datasets: + *model_datasets + si_obs: + short_name: siconc + preprocessor: pp_arctic_sept_sea_ice + additional_datasets: + *arctic_si_obs + tasa_obs: + short_name: tasa + preprocessor: pp_avg_ann_global_temp + mip: Amon + additional_datasets: + *tasa_obs scripts: sea_ice_sensitivity_script: script: seaice/seaice_sensitivity.py @@ -124,22 +148,6 @@ diagnostics: mean: -4.01 standard deviation: 0.32 plausible range: 1.28 - annual trends (Roach-style plot): - first point: - GMST trend: - SIA trend: - Pearson CC of SIA over GMST: - significance of SIA over GMST: - second point: - GMST trend: - SIA trend: - Pearson CC of SIA over GMST: - significance of SIA over GMST: - third point: - GMST trend: - SIA trend: - Pearson CC of SIA over GMST: - significance of SIA over GMST: antarctic: @@ -148,9 +156,24 @@ diagnostics: siconc: preprocessor: pp_antarctic_avg_ann_sea_ice mip: SImon + additional_datasets: + *model_datasets tas: preprocessor: pp_avg_ann_global_temp mip: Amon + additional_datasets: + *model_datasets +# si_obs: +# short_name: siconc +# preprocessor: pp_antarctic_avg_ann_sea_ice +# additional_datasets: +# *antarctic_si_obs + tasa_obs: + short_name: tasa + preprocessor: pp_avg_ann_global_temp + mip: Amon + additional_datasets: + *tasa_obs scripts: sea_ice_sensitivity_script: script: seaice/seaice_sensitivity.py @@ -160,19 +183,3 @@ diagnostics: mean: standard deviation: plausible range: - annual trends (Roach-style plot): - first point: - GMST trend: - SIA trend: - Pearson CC of SIA over GMST: - significance of SIA over GMST: - second point: - GMST trend: - SIA trend: - Pearson CC of SIA over GMST: - significance of SIA over GMST: - third point: - GMST trend: - SIA trend: - Pearson CC of SIA over GMST: - significance of SIA over GMST: From 9c2103b29af63aae59993535b540765330cd58c7 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Tue, 14 Oct 2025 10:10:15 +0100 Subject: [PATCH 02/26] TEMP --- .../diag_scripts/seaice/seaice_sensitivity.py | 108 +++++++++++++----- 1 file changed, 82 insertions(+), 26 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 6b3b1053e6..76f735db32 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -10,14 +10,14 @@ from matplotlib.colors import Normalize from scipy import stats -from esmvaltool.diag_scripts.shared import ( - ProvenanceLogger, - run_diagnostic, - save_figure, - select_metadata, -) - -logger = logging.getLogger(Path(__file__).stem) +# from esmvaltool.diag_scripts.shared import ( +# ProvenanceLogger, +# run_diagnostic, +# save_figure, +# select_metadata, +# ) +# +# logger = logging.getLogger(Path(__file__).stem) # This is stolen from the example in AutoAssess _plot_mo_metrics.py @@ -44,14 +44,18 @@ def create_categorised_dataset_dict(data): # Initialise blank dictionary dict = { 'models': {}, - 'observations': {}, + 'tasa_obs': {}, + 'siconc_obs': {}, } # Iterate over the data for element in data: if 'observation' in element: - if element['observation']: - dict['observations'][element['dataset']] = {} + if element['observation']: # In case user writes "false" or similar in recipe + if element['variable_group'] == 'tasa_obs': + dict['tasa_obs'][element['dataset']] = {} + elif element['variable_group'] == 'si_obs': + dict['siconc_obs'][element['dataset']] = {} else: dict['models'][element['dataset']] = {} @@ -394,47 +398,99 @@ def main(cfg): # Get list of datasets from cfg logger.info("Listing datasets in the data") dataset_dict = create_categorised_dataset_dict(input_data) - data_dict = dataset_dict['models'] - - # Iterate over each dataset - for dataset in data_dict.keys(): + # Iterate over each model dataset (has both tas and siconc) + logger.info("Calculating from model data") + model_dict = dataset_dict['models'] + for dataset in model_dict.keys(): # Select only data from that dataset logger.debug("Selecting data from %s", dataset) selection = select_metadata(input_data, dataset=dataset) # Add an entry to determine labelling in plots if "label_dataset" in selection[0]: - data_dict[dataset]["label"] = "to_label" + model_dict[dataset]["label"] = "to_label" logger.info("Dataset %s will be labelled", dataset) else: - data_dict[dataset]["label"] = "unlabelled" + model_dict[dataset]["label"] = "unlabelled" logger.info("Not labelling dataset %s in plots", dataset) # Calculations for the Notz-style plot logger.info("Calculating data for Notz-style plot") sensitivity = calculate_direct_sensitivity(selection) # Add to dictionary - data_dict[dataset]["direct_sensitivity_(notz-style)"] = sensitivity + model_dict[dataset]["direct_sensitivity_(notz-style)"] = sensitivity # Calculations for the Roach-style plot logger.info("Calculating data for Roach-style plot") trends = calculate_annual_trends(selection) # Add to dictionary - data_dict[dataset]["annual_siconc_trend"] = trends["si_ann_trend"] - data_dict[dataset]["annual_tas_trend"] = trends["tas_ann_trend"] - data_dict[dataset]["direct_r_val"] = trends["direct_r_val"] - data_dict[dataset]["direct_p_val"] = trends["direct_p_val"] + model_dict[dataset]["annual_siconc_trend"] = trends["si_ann_trend"] + model_dict[dataset]["annual_tas_trend"] = trends["tas_ann_trend"] + model_dict[dataset]["direct_r_val"] = trends["direct_r_val"] + model_dict[dataset]["direct_p_val"] = trends["direct_p_val"] + + # Iterate over each observational dataset (has only one variable) + + # tasa observations + logger.info("Calculating from tasa observational data") + tasa_dict = dataset_dict['tasa_obs'] + for dataset in tasa_dict.keys(): + # Select only data from that dataset + logger.debug("Selecting data from %s", dataset) + selection = select_metadata(input_data, dataset=dataset) + + # Add an entry to determine labelling in plots + if "label_dataset" in selection[0]: + tasa_dict[dataset]["label"] = "to_label" + logger.info("Dataset %s will be labelled", dataset) + else: + tasa_dict[dataset]["label"] = "unlabelled" + logger.info("Not labelling dataset %s in plots", dataset) + + # Calculations for the Roach-style plot + logger.info("Calculating data for Roach-style plot") + tasa_cube = extract_cube(input_data, "tasa") + years = tasa_cube.coord("year").points + tasa_trend = calculate_regression(years, tasa_cube.data) + # Add to dictionary + tasa_dict[dataset]["annual_tasa_trend"] = tasa_trend.slope + + # siconc observations + logger.info("Calculating from siconc observational data") + siconc_dict = dataset_dict['siconc_obs'] + for dataset in siconc_dict.keys(): + # Select only data from that dataset + logger.debug("Selecting data from %s", dataset) + selection = select_metadata(input_data, dataset=dataset) + + # Add an entry to determine labelling in plots + if "label_dataset" in selection[0]: + siconc_dict[dataset]["label"] = "to_label" + logger.info("Dataset %s will be labelled", dataset) + else: + siconc_dict[dataset]["label"] = "unlabelled" + logger.info("Not labelling dataset %s in plots", dataset) + + # Calculations for the Roach-style plot + logger.info("Calculating data for Roach-style plot") + siconc_cube = extract_cube(input_data, "siconc") + years = siconc_cube.coord("year").points + siconc_trend = calculate_regression(years, siconc_cube.data) + # Add to dictionary + siconc_dict[dataset]["annual_siconc_trend"] = siconc_trend.slope # Add the values to plot to a csv file logger.info("Writing values to csv") - write_dictionary_to_csv(cfg, data_dict, "plotted_values") + write_dictionary_to_csv(cfg, dataset_dict, "plotted_values") - # Plot the sensitivities (and save and close the plots) + # Plot the sensitivities (Notz-style, involves model data only) logger.info("Creating Notz-style plot") - notz_style_plot_from_dict(data_dict, titles_and_obs_dict, cfg) + notz_style_plot_from_dict(model_dict, titles_and_obs_dict, cfg) + + # Plot the trends (Roach-style, involves model and observational data ) logger.info("Creating Roach-style plot") - roach_style_plot_from_dict(data_dict, titles_and_obs_dict, cfg) + roach_style_plot_from_dict(dataset_dict, titles_and_obs_dict, cfg) if __name__ == "__main__": From 86c5c84892ec964cae7566e60d5b464e5b9c0885 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Fri, 17 Oct 2025 09:06:20 +0100 Subject: [PATCH 03/26] Interim version --- .../diag_scripts/seaice/seaice_sensitivity.py | 526 ++++-------------- .../recipes/recipe_seaice_sensitivity.yml | 10 +- 2 files changed, 118 insertions(+), 418 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 76f735db32..1ed2666a08 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -6,18 +6,15 @@ import iris import matplotlib.pyplot as plt import pandas as pd -import yaml from matplotlib.colors import Normalize -from scipy import stats +from scipy.stats import linregress -# from esmvaltool.diag_scripts.shared import ( -# ProvenanceLogger, -# run_diagnostic, -# save_figure, -# select_metadata, -# ) -# -# logger = logging.getLogger(Path(__file__).stem) +from esmvaltool.diag_scripts.shared import ( + ProvenanceLogger, + run_diagnostic, +) + +logger = logging.getLogger(Path(__file__).stem) # This is stolen from the example in AutoAssess _plot_mo_metrics.py @@ -38,171 +35,141 @@ def get_provenance_record(cfg, caption): return record -def create_categorised_dataset_dict(data): - """Create a dictionary of datasets categorised into models and observations.""" - logger.debug("Creating categorised dataset dictionary from %s", data) - # Initialise blank dictionary - dict = { +def create_category_dict(cfg): + """Create a structured dictionary for adding values to later on.""" + # Create blank dictionary of correct structure + category_dict = { 'models': {}, 'tasa_obs': {}, 'siconc_obs': {}, } - # Iterate over the data - for element in data: - if 'observation' in element: - if element['observation']: # In case user writes "false" or similar in recipe - if element['variable_group'] == 'tasa_obs': - dict['tasa_obs'][element['dataset']] = {} - elif element['variable_group'] == 'si_obs': - dict['siconc_obs'][element['dataset']] = {} - else: - dict['models'][element['dataset']] = {} - - return dict - - -def extract_cube(data, variable_group): - """Return the variable group's single iris cube from the data.""" - logger.debug("extracting %s cube from %s", variable_group, data) - - # Load any data in the variable_group - selection = select_metadata(data, variable_group=variable_group) + # Initialize the models as a set to avoid duplication + models_set = set() - # Ensure there is only one file in the list - if len(selection) != 1: - raise ValueError( - f"None or too many matching files found for {variable_group}" - ) + # Read the data from the config object + input_data = cfg["input_data"].values() - # Load the cube, [0] is because selection returns a list - cube = iris.load_cube(selection[0]["filename"]) + # Iterate over the datasets to add to the dictionary + for input in input_data: - return cube + # Check for tasa observations + if input['variable_group'] == 'tasa_obs': + category_dict['tasa_obs'][input['dataset']] = {} + # Check for siconc observations + elif input['variable_group'] == 'siconc_obs': + category_dict['siconc_obs'][input['dataset']] = {} -def calculate_regression(independent, dependent): - """Use SciPy stats to calculate the least-squares regression.""" - logger.debug( - "Calculating linear relationship between %s and %s", - dependent, - independent, - ) - - # Use SciPy stats to calculate the regression - # result = slope, intercept, rvalue, pvalue, stderr - result = stats.linregress(independent, dependent) + # Everything else should be a model + else: + models_set.add(input['dataset']) - # Return everything - return result + # Add the models to the dictionary + for model in models_set: + category_dict['models'][model] = {} + return category_dict -def calculate_annual_trends(data): - """ - Calculate annual trends for surface air temperature (tas) and sea ice area (siconc). - Also used for the r and p values from the regression of siconc as a function of tas. - """ - logger.debug("calculating annual trends") +def fetch_cube(dataset, variable, cfg): + """Fetch a data cube for a dataset and variable using info from the config.""" + # Read the data from the config object + input_data = cfg["input_data"].values() - # Load the preprocessed cubes - si_cube = extract_cube(data, "siconc") - tas_cube = extract_cube(data, "tas") + # Find the correct filepath for the dataset + for input in input_data: + if input['dataset'] == dataset: + # Also check the variable matches as models have two entries + # Only matching the first three letters to avoid issues with sic vs siconc + if input['short_name'][:3] == variable[:3]: + filepath = input['filename'] + break + + # Load the cube using iris + cube = iris.load_cube(filepath, variable) + return cube - # Calculate the individual trends over time - years = tas_cube.coord("year").points - si_trend = calculate_regression(years, si_cube.data) - tas_trend = calculate_regression(years, tas_cube.data) - # Calculate the direct regression for r and p values - direct_regression = calculate_regression(tas_cube.data, si_cube.data) +def calculate_annual_trend(cube): + """Calculate the linear trend of a cube over time using scipy.stats.linregres.""" + # Depending on preprocessor, coord may be 'year' or 'time' + if 'year' in cube.coords(): + years = cube.coord("year").points + else: + years = cube.coord("time").points - dictionary = { - "si_ann_trend": si_trend.slope, - "tas_ann_trend": tas_trend.slope, - "direct_r_val": direct_regression.rvalue, - "direct_p_val": direct_regression.pvalue, - } + # slope, intercept, rvalue, pvalue, stderr = linregress(independent, dependent) + trend = linregress(years, cube.data) - return dictionary + # Only the slope is needed in this code + return trend.slope -def calculate_direct_sensitivity(data): - """Calculate slope of sea ice area over global mean temperature.""" - logger.debug("calculating sensitivity") +def calculate_direct_stats(dataset, cfg): + """Calculate the direct sensitivity of siconc to tas for a given dataset.""" + # Fetch the required cubes + siconc_cube = fetch_cube(dataset, 'siconc', cfg) + tas_cube = fetch_cube(dataset, 'tas', cfg) - # Load the preprocessed cubes - si_cube = extract_cube(data, "siconc") - tas_cube = extract_cube(data, "tas") + # Calculate direct regression (tas as independent) + direct_sensitivity = linregress(tas_cube.data, siconc_cube.data) - # Calculate the slope of the direct regression, NOT via time - sensitivity = calculate_regression(tas_cube.data, si_cube.data) + # direct_sensitivity = slope, intercept, rvalue, pvalue, stderr + return direct_sensitivity - return sensitivity.slope +def write_values_to_dict(data_dict, cfg): + """Calculate and write values to the structured dictionary.""" + # Calculate all the values for the models + for model_dataset in data_dict['models'].keys(): + # Calculate annual tas trend + tas_cube = fetch_cube(model_dataset, 'tas', cfg) + ann_tas_trend = calculate_annual_trend(tas_cube) + data_dict['models'][model_dataset]['annual_tas_trend'] = ann_tas_trend -def write_obs_from_cfg(cfg): - """Write the observations from the recipe to a dictionary.""" - # Initialize the dictionary to hold observations - obs_dict = {} + # Calculate annual siconc trend + siconc_cube = fetch_cube(model_dataset, 'siconc', cfg) + ann_siconc_trend = calculate_annual_trend(siconc_cube) + data_dict['models'][model_dataset]['annual_siconc_trend'] = ann_siconc_trend - # Add the observation period to the dictionary - obs_dict["obs_period"] = cfg["observations"]["observation period"] + # Calculate direct sensitivity of siconc to tas + direct_sensitivity = calculate_direct_stats(model_dataset, cfg) + data_dict['models'][model_dataset]['direct_sensitivity'] = direct_sensitivity.slope + data_dict['models'][model_dataset]['direct_r_value'] = direct_sensitivity.rvalue + data_dict['models'][model_dataset]['direct_p_value'] = direct_sensitivity.pvalue - # Add a blank dictionary for the Notz-style plot - obs_dict["notz_style"] = {} + # Calculate just the tasa trend for the tasa observations + for obs_dataset in data_dict['tasa_obs'].keys(): + # Calculate annual tas trend + tasa_cube = fetch_cube(obs_dataset, 'tasa', cfg) + ann_tasa_trend = calculate_annual_trend(tasa_cube) + # Still labelling in final dictionary as tas for consistency with models + data_dict['tasa_obs'][obs_dataset]['annual_tas_trend'] = ann_tasa_trend - # Add the observations values to the dictionary - notz_values = cfg["observations"]["sea ice sensitivity (Notz-style plot)"] - obs_dict["notz_style"]["mean"] = notz_values["mean"] - obs_dict["notz_style"]["std_dev"] = notz_values["standard deviation"] - obs_dict["notz_style"]["plausible"] = notz_values["plausible range"] + # Calculate just the siconc trend for the siconc observations + for obs_dataset in data_dict['siconc_obs'].keys(): + # Calculate annual siconc trend + siconc_cube = fetch_cube(obs_dataset, 'siconc', cfg) + ann_siconc_trend = calculate_annual_trend(siconc_cube) + data_dict['siconc_obs'][obs_dataset]['annual_siconc_trend'] = ann_siconc_trend - return obs_dict + return data_dict -def create_titles_dict(data, cfg): +def write_dictionary_to_csv(sub_dict, filename, cfg): """ - Create a dictionary of appropriate observations and titles. + Output a section of data dictionary to a csv file using Pandas. - Values depend on whether the plot is for the Arctic or Antarctic - and assume the recipe used September Arctic sea ice data or - annually mean averaged Antarctic sea ice data + Only sections of the dictionary should be written at a time as otherwise + the structure is too complex to easily convert to a DataFrame. """ - dictionary = {} - - first_variable = next(iter(data)) - - if first_variable["diagnostic"] == "arctic": - dictionary["titles"] = { - "notz_fig_title": "September Arctic sea-ice area sensitivity to global mean surface temperature", - "notz_ax_title": "dSIA/dGMST", - "notz_plot_filename": "September Arctic sea ice sensitivity", - "roach_fig_title": "Trends in Annual Mean Temperature And September Arctic Sea Ice", - "roach_plot_filename": "September Arctic sea ice trends", - } - - elif first_variable["diagnostic"] == "antarctic": - dictionary["titles"] = { - "notz_fig_title": "Annually Meaned Antarctic Sea Ice Sensitivity", - "notz_ax_title": "dSIA/dGMST", - "notz_plot_filename": "Annual Antarctic sea ice sensitivity", - "roach_fig_title": "Trends in Annual Mean Temperature And Annual Antarctic Sea Ice", - "roach_plot_filename": "Annual Antarctic sea ice trends", - } - - dictionary["obs"] = write_obs_from_cfg(cfg) - - return dictionary - - -def write_dictionary_to_csv(cfg, model_dict, filename): - """Output the model dictionary to a csv file using Pandas.""" - # Read the work directory from the config and create the csv filepath + # Create the csv filepath using info from the config csv_filepath = f"{cfg['work_dir']}/{filename}.csv" # Write the data to a csv file (via a Pandas DataFrame) - pd.DataFrame.from_dict(model_dict, orient="index").to_csv(csv_filepath, index_label="Dataset") + dataframe = pd.DataFrame.from_dict(sub_dict, orient="index") + dataframe.to_csv(csv_filepath) logger.info("Wrote data to %s", csv_filepath) # Create a provenance record for the csv file @@ -213,284 +180,17 @@ def write_dictionary_to_csv(cfg, model_dict, filename): ) -def notz_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): - """Save a plot of sensitivities and observations.""" - # Read from observations dictionary - obs_years = titles_dictionary["obs"]["obs_period"] - obs_dict = titles_dictionary["obs"]["notz_style"] - obs_mean = obs_dict["mean"] - obs_std_dev = obs_dict["std_dev"] - obs_plausible = obs_dict["plausible"] - - # Set up the figure - fig, ax = plt.subplots(figsize=(3.5, 6), layout="constrained") - fig.suptitle(titles_dictionary["titles"]["notz_fig_title"], wrap=True) - ax.set_title( - titles_dictionary["titles"]["notz_ax_title"], wrap=True, fontsize=10 - ) - - # Iterate over the dictionary - for dataset, inner_dict in data_dictionary.items(): - ax.plot( - 0.25, - inner_dict["direct_sensitivity_(notz-style)"], - color="blue", - marker="_", - markersize=20, - ) - - # Label with the dataset if specified - if inner_dict["label"] == "to_label": - plt.annotate( - dataset, - xy=(0.25, inner_dict["direct_sensitivity_(notz-style)"]), - xytext=( - 0.35, - inner_dict["direct_sensitivity_(notz-style)"] - 0.05, - ), - ) - - # Add observations only if obs_mean is a number - if isinstance(obs_mean, int | float): - ax.hlines(obs_mean, 0, 1, linestyle="--", color="black", linewidth=2) - ax.fill_between( - [0, 1], - obs_mean - obs_std_dev, - obs_mean + obs_std_dev, - facecolor="k", - alpha=0.15, - ) - ax.hlines( - obs_mean + obs_plausible, - 0, - 1, - linestyle=":", - color="0.5", - linewidth=1, - ) - ax.hlines( - obs_mean - obs_plausible, - 0, - 1, - linestyle=":", - color="0.5", - linewidth=1, - ) - - # Tidy the figure - ax.set_xlim(0, 1) - ax.set_xticks([]) - ax.set_ylabel(r"dSIA/dGMST ($million \ km^2 \ K^{-1}$)") - - # Create caption based on whether observation mean is presnt - if isinstance(obs_mean, int | float): - caption = ( - "Sensitivity of sea ice area to annual mean global warming." - f"Mean (dashed), standard deviation (shaded) and plausible values from {obs_years}." - ) - else: - caption = "Sensitivity of sea ice area to annual mean global warming." - - # Save the figure (also closes it) - save_figure( - titles_dictionary["titles"]["notz_plot_filename"], - get_provenance_record(cfg, caption), - cfg, - figure=fig, - close=True, - ) - - -def roach_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): - """Save a plot of trend in SIA against trend in GMST to the given filename.""" - # Set up the figure - fig, ax = plt.subplots(figsize=(10, 6), layout="constrained") - fig.suptitle(titles_dictionary["titles"]["roach_fig_title"], wrap=True) - - # Set up for colouring the points - norm = Normalize(vmin=-1, vmax=1) - cmap = plt.get_cmap("PiYG_r") - - # Set up the axes - ax.axhline(color="black", alpha=0.5) - ax.axvline(color="black", alpha=0.5) - ax.set_xlabel(r"Trend in GMST ($K \ decade^{-1}$)") - ax.set_ylabel(r"Trend in SIA ($million \ km^2 \ decade^{-1}$)") - - # Iterate over the dictionary - for dataset, inner_dict in data_dictionary.items(): - # Determine the position of the point - x = 10 * inner_dict["annual_tas_trend"] # for equivalence to decades - y = ( - 10 * inner_dict["annual_siconc_trend"] - ) # for equivalence to decades - - # Determine the colour of the point - r_corr = inner_dict["direct_r_val"] - - # Decide if the point should be hatched - if inner_dict["direct_p_val"] >= 0.05: - h = 5 * "/" # This is a hatch pattern - else: - h = None - - # Plot the point - plt.scatter( - x, y, marker="o", s=150, c=[r_corr], hatch=h, cmap=cmap, norm=norm - ) - - # Label with the dataset if specified - if inner_dict["label"] == "to_label": - plt.annotate(dataset, xy=(x, y), xytext=(x + 0.01, y - 0.005)) - - # # Provide a default colour for the point if Pearson coefficient is missing - # if r_corr is None: - # r_corr = 0 - # - # # Provide a pattern for the point if the p-value is present and sufficiently large - # if p_val is not None and p_val >= 0.05: - # h = 5 * "/" # This is a hatch pattern - # else: - # h = None - # - # # Plot the point only if both x and y values are provided - # if x is not None and y is not None: - # plt.scatter( - # x, - # y, - # marker="s", - # s=150, - # c=[r_corr], - # hatch=h, - # cmap=cmap, - # norm=norm, - # zorder=0, - # edgecolors="black", - # ) - - # Add a colour bar - plt.colorbar(label="Pearson correlation coefficient") - - # Create caption - caption = "Decadal trends of sea ice area and global mean temperature." - - # Save the figure (also closes it) - save_figure( - titles_dictionary["titles"]["roach_plot_filename"], - get_provenance_record(cfg, caption), - cfg, - figure=fig, - close=True, - ) - - def main(cfg): - """Create two plots per diagnostic from preprocessed data.""" - # Get the data from the cfg - logger.info("Getting data from the config") - input_data = cfg["input_data"].values() + # Create the structured dictionary + data_dict = create_category_dict(cfg) - # Titles and observations depend on the diagnostic being plotted - logger.info("Creating titles and observations dictionary") - titles_and_obs_dict = create_titles_dict(input_data, cfg) - logger.debug("Titles and observations dictionary: %s", titles_and_obs_dict) - - # Get list of datasets from cfg - logger.info("Listing datasets in the data") - dataset_dict = create_categorised_dataset_dict(input_data) - - # Iterate over each model dataset (has both tas and siconc) - logger.info("Calculating from model data") - model_dict = dataset_dict['models'] - for dataset in model_dict.keys(): - # Select only data from that dataset - logger.debug("Selecting data from %s", dataset) - selection = select_metadata(input_data, dataset=dataset) - - # Add an entry to determine labelling in plots - if "label_dataset" in selection[0]: - model_dict[dataset]["label"] = "to_label" - logger.info("Dataset %s will be labelled", dataset) - else: - model_dict[dataset]["label"] = "unlabelled" - logger.info("Not labelling dataset %s in plots", dataset) - - # Calculations for the Notz-style plot - logger.info("Calculating data for Notz-style plot") - sensitivity = calculate_direct_sensitivity(selection) - # Add to dictionary - model_dict[dataset]["direct_sensitivity_(notz-style)"] = sensitivity - - # Calculations for the Roach-style plot - logger.info("Calculating data for Roach-style plot") - trends = calculate_annual_trends(selection) - # Add to dictionary - model_dict[dataset]["annual_siconc_trend"] = trends["si_ann_trend"] - model_dict[dataset]["annual_tas_trend"] = trends["tas_ann_trend"] - model_dict[dataset]["direct_r_val"] = trends["direct_r_val"] - model_dict[dataset]["direct_p_val"] = trends["direct_p_val"] - - # Iterate over each observational dataset (has only one variable) - - # tasa observations - logger.info("Calculating from tasa observational data") - tasa_dict = dataset_dict['tasa_obs'] - for dataset in tasa_dict.keys(): - # Select only data from that dataset - logger.debug("Selecting data from %s", dataset) - selection = select_metadata(input_data, dataset=dataset) - - # Add an entry to determine labelling in plots - if "label_dataset" in selection[0]: - tasa_dict[dataset]["label"] = "to_label" - logger.info("Dataset %s will be labelled", dataset) - else: - tasa_dict[dataset]["label"] = "unlabelled" - logger.info("Not labelling dataset %s in plots", dataset) - - # Calculations for the Roach-style plot - logger.info("Calculating data for Roach-style plot") - tasa_cube = extract_cube(input_data, "tasa") - years = tasa_cube.coord("year").points - tasa_trend = calculate_regression(years, tasa_cube.data) - # Add to dictionary - tasa_dict[dataset]["annual_tasa_trend"] = tasa_trend.slope - - # siconc observations - logger.info("Calculating from siconc observational data") - siconc_dict = dataset_dict['siconc_obs'] - for dataset in siconc_dict.keys(): - # Select only data from that dataset - logger.debug("Selecting data from %s", dataset) - selection = select_metadata(input_data, dataset=dataset) - - # Add an entry to determine labelling in plots - if "label_dataset" in selection[0]: - siconc_dict[dataset]["label"] = "to_label" - logger.info("Dataset %s will be labelled", dataset) - else: - siconc_dict[dataset]["label"] = "unlabelled" - logger.info("Not labelling dataset %s in plots", dataset) - - # Calculations for the Roach-style plot - logger.info("Calculating data for Roach-style plot") - siconc_cube = extract_cube(input_data, "siconc") - years = siconc_cube.coord("year").points - siconc_trend = calculate_regression(years, siconc_cube.data) - # Add to dictionary - siconc_dict[dataset]["annual_siconc_trend"] = siconc_trend.slope - - # Add the values to plot to a csv file - logger.info("Writing values to csv") - write_dictionary_to_csv(cfg, dataset_dict, "plotted_values") - - # Plot the sensitivities (Notz-style, involves model data only) - logger.info("Creating Notz-style plot") - notz_style_plot_from_dict(model_dict, titles_and_obs_dict, cfg) - - # Plot the trends (Roach-style, involves model and observational data ) - logger.info("Creating Roach-style plot") - roach_style_plot_from_dict(dataset_dict, titles_and_obs_dict, cfg) + # Calculate and write values to the dictionary + data_dict = write_values_to_dict(data_dict, cfg) + + # Write the model and obs dictionaries to csv files + write_dictionary_to_csv(data_dict['models'], 'models_values', cfg) + write_dictionary_to_csv(data_dict['tasa_obs'], 'tasa_obs_values', cfg) + write_dictionary_to_csv(data_dict['siconc_obs'], 'siconc_obs_values', cfg) if __name__ == "__main__": diff --git a/esmvaltool/recipes/recipe_seaice_sensitivity.yml b/esmvaltool/recipes/recipe_seaice_sensitivity.yml index 1a232805b0..46e19d3bce 100644 --- a/esmvaltool/recipes/recipe_seaice_sensitivity.yml +++ b/esmvaltool/recipes/recipe_seaice_sensitivity.yml @@ -46,7 +46,7 @@ tasa_obs: &tasa_obs - { <<: *obs_defaults, dataset: HadCRUT5, type: ground, version: 5.0.1.0-analysis } - { <<: *obs_defaults, dataset: NOAAGlobalTemp, type: ground, version: v5.0.0, mip: Amon } -arctic_si_obs: &arctic_si_obs +arctic_siconc_obs: &arctic_siconc_obs - { <<: *obs_defaults, dataset: OSI-450-nh, type: reanaly, version: v3, mip: OImon, supplementary_variables: [{short_name: areacello, mip: fx}]} preprocessors: @@ -128,11 +128,11 @@ diagnostics: mip: Amon additional_datasets: *model_datasets - si_obs: + siconc_obs: short_name: siconc preprocessor: pp_arctic_sept_sea_ice additional_datasets: - *arctic_si_obs + *arctic_siconc_obs tasa_obs: short_name: tasa preprocessor: pp_avg_ann_global_temp @@ -163,11 +163,11 @@ diagnostics: mip: Amon additional_datasets: *model_datasets -# si_obs: +# siconc_obs: # short_name: siconc # preprocessor: pp_antarctic_avg_ann_sea_ice # additional_datasets: -# *antarctic_si_obs +# *antarctic_siconc_obs tasa_obs: short_name: tasa preprocessor: pp_avg_ann_global_temp From 58362450a44ca2ef81c1b0680e22b5917c9f786c Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Fri, 21 Nov 2025 10:30:18 +0000 Subject: [PATCH 04/26] printing dictionary --- .../diag_scripts/seaice/seaice_sensitivity.py | 15 +++++++++++++++ .../recipes/recipe_seaice_sensitivity.yml | 17 ++++++++++------- 2 files changed, 25 insertions(+), 7 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 1ed2666a08..d88ff83e28 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -119,6 +119,19 @@ def calculate_direct_stats(dataset, cfg): return direct_sensitivity +def calculate_cross_dataset_stats(tasa_dataset, siconc_dataset, cfg): + """Calculate the sensitivity of siconc to tasa across (obs) datasets.""" + # Fetch the required cubes + siconc_cube = fetch_cube(dataset, 'siconc', cfg) + tas_cubea = fetch_cube(dataset, 'tas', cfg) + + # Calculate direct regression (tas as independent) + direct_sensitivity = linregress(tasa_cube.data, siconc_cube.data) + + # direct_sensitivity = slope, intercept, rvalue, pvalue, stderr + return direct_sensitivity + + def write_values_to_dict(data_dict, cfg): """Calculate and write values to the structured dictionary.""" # Calculate all the values for the models @@ -192,6 +205,8 @@ def main(cfg): write_dictionary_to_csv(data_dict['tasa_obs'], 'tasa_obs_values', cfg) write_dictionary_to_csv(data_dict['siconc_obs'], 'siconc_obs_values', cfg) + # TEMPORARY: Print the dictionary to the log for checking + logger.debug("Final data dictionary:\n%s", data_dict) if __name__ == "__main__": with run_diagnostic() as config: diff --git a/esmvaltool/recipes/recipe_seaice_sensitivity.yml b/esmvaltool/recipes/recipe_seaice_sensitivity.yml index 46e19d3bce..c40fd6f7ea 100644 --- a/esmvaltool/recipes/recipe_seaice_sensitivity.yml +++ b/esmvaltool/recipes/recipe_seaice_sensitivity.yml @@ -42,13 +42,16 @@ model_datasets: &model_datasets obs_defaults: &obs_defaults { project: OBS, observation: True, tier: 2} tasa_obs: &tasa_obs - - { <<: *obs_defaults, dataset: GISTEMP, type: ground, version: v4 } - - { <<: *obs_defaults, dataset: HadCRUT5, type: ground, version: 5.0.1.0-analysis } + - { <<: *obs_defaults, dataset: GISTEMP, type: ground, version: v4 } # Doesn't match Roach obs + - { <<: *obs_defaults, dataset: HadCRUT5, type: ground, version: 5.0.1.0-analysis } # Doesn't match Roach obs - { <<: *obs_defaults, dataset: NOAAGlobalTemp, type: ground, version: v5.0.0, mip: Amon } arctic_siconc_obs: &arctic_siconc_obs - { <<: *obs_defaults, dataset: OSI-450-nh, type: reanaly, version: v3, mip: OImon, supplementary_variables: [{short_name: areacello, mip: fx}]} +antarctic_siconc_obs: &antarctic_siconc_obs + - { <<: *obs_defaults, dataset: OSI-450-sh, type: reanaly, version: v3, mip: OImon, supplementary_variables: [{short_name: areacello, mip: fx}]} + preprocessors: extract_test_period: &extract_test_period @@ -163,11 +166,11 @@ diagnostics: mip: Amon additional_datasets: *model_datasets -# siconc_obs: -# short_name: siconc -# preprocessor: pp_antarctic_avg_ann_sea_ice -# additional_datasets: -# *antarctic_siconc_obs + siconc_obs: + short_name: siconc + preprocessor: pp_antarctic_avg_ann_sea_ice + additional_datasets: + *antarctic_siconc_obs tasa_obs: short_name: tasa preprocessor: pp_avg_ann_global_temp From 72eb2110119b03430289d1549cb5b9252450af0a Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Fri, 21 Nov 2025 10:50:44 +0000 Subject: [PATCH 05/26] Trying to fix trends --- .../diag_scripts/seaice/seaice_sensitivity.py | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index d88ff83e28..25308102fa 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -95,12 +95,12 @@ def calculate_annual_trend(cube): """Calculate the linear trend of a cube over time using scipy.stats.linregres.""" # Depending on preprocessor, coord may be 'year' or 'time' if 'year' in cube.coords(): - years = cube.coord("year").points + no_years = [i for i in range(len(cube.coord("years").points))] else: - years = cube.coord("time").points + no_years = [i for i in range(len(cube.coord("time").points))] # slope, intercept, rvalue, pvalue, stderr = linregress(independent, dependent) - trend = linregress(years, cube.data) + trend = linregress(no_years, cube.data) # Only the slope is needed in this code return trend.slope @@ -119,17 +119,17 @@ def calculate_direct_stats(dataset, cfg): return direct_sensitivity -def calculate_cross_dataset_stats(tasa_dataset, siconc_dataset, cfg): - """Calculate the sensitivity of siconc to tasa across (obs) datasets.""" - # Fetch the required cubes - siconc_cube = fetch_cube(dataset, 'siconc', cfg) - tas_cubea = fetch_cube(dataset, 'tas', cfg) - - # Calculate direct regression (tas as independent) - direct_sensitivity = linregress(tasa_cube.data, siconc_cube.data) - - # direct_sensitivity = slope, intercept, rvalue, pvalue, stderr - return direct_sensitivity +# def calculate_cross_dataset_stats(tasa_dataset, siconc_dataset, cfg): +# """Calculate the sensitivity of siconc to tasa across (obs) datasets.""" +# # Fetch the required cubes +# siconc_cube = fetch_cube(dataset, 'siconc', cfg) +# tas_cubea = fetch_cube(dataset, 'tas', cfg) +# +# # Calculate direct regression (tas as independent) +# direct_sensitivity = linregress(tasa_cube.data, siconc_cube.data) +# +# # direct_sensitivity = slope, intercept, rvalue, pvalue, stderr +# return direct_sensitivity def write_values_to_dict(data_dict, cfg): From b731cb623bd57076b61a71aeda99aee8938ec85d Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Fri, 21 Nov 2025 11:23:01 +0000 Subject: [PATCH 06/26] Cross-obs values --- .../diag_scripts/seaice/seaice_sensitivity.py | 40 ++++++++++++++----- 1 file changed, 29 insertions(+), 11 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 25308102fa..6f586c6df9 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -119,17 +119,17 @@ def calculate_direct_stats(dataset, cfg): return direct_sensitivity -# def calculate_cross_dataset_stats(tasa_dataset, siconc_dataset, cfg): -# """Calculate the sensitivity of siconc to tasa across (obs) datasets.""" -# # Fetch the required cubes -# siconc_cube = fetch_cube(dataset, 'siconc', cfg) -# tas_cubea = fetch_cube(dataset, 'tas', cfg) -# -# # Calculate direct regression (tas as independent) -# direct_sensitivity = linregress(tasa_cube.data, siconc_cube.data) -# -# # direct_sensitivity = slope, intercept, rvalue, pvalue, stderr -# return direct_sensitivity +def calculate_cross_dataset_stats(tasa_dataset, siconc_dataset, cfg): + """Calculate the sensitivity of siconc to tasa across (obs) datasets.""" + # Fetch the required cubes + siconc_cube = fetch_cube(siconc_dataset, 'siconc', cfg) + tasa_cube = fetch_cube(tasa_dataset, 'tasa', cfg) + + # Calculate direct regression (tasa as independent) + direct_sensitivity = linregress(tasa_cube.data, siconc_cube.data) + + # direct_sensitivity = slope, intercept, rvalue, pvalue, stderr + return direct_sensitivity def write_values_to_dict(data_dict, cfg): @@ -167,6 +167,19 @@ def write_values_to_dict(data_dict, cfg): ann_siconc_trend = calculate_annual_trend(siconc_cube) data_dict['siconc_obs'][obs_dataset]['annual_siconc_trend'] = ann_siconc_trend + # Calculate cross-dataset statistics between tasa and siconc observations + for tasa_dataset in data_dict['tasa_obs'].keys(): + for siconc_dataset in data_dict['siconc_obs'].keys(): + # Calculate cross-dataset sensitivity of siconc to tasa + cross_sensitivity = calculate_cross_dataset_stats(tasa_dataset, siconc_dataset, cfg) + # Determine structure of dictionary to store values + key_name = f"{siconc_dataset}_to_{tasa_dataset}" + inner_dict = data_dict['cross-dataset-obs'][key_name] + # Store the values + inner_dict['direct_sensitivity'] = cross_sensitivity.slope + inner_dict['direct_r_value'] = cross_sensitivity.rvalue + inner_dict['direct_p_value'] = cross_sensitivity.pvalue + return data_dict @@ -205,6 +218,11 @@ def main(cfg): write_dictionary_to_csv(data_dict['tasa_obs'], 'tasa_obs_values', cfg) write_dictionary_to_csv(data_dict['siconc_obs'], 'siconc_obs_values', cfg) + # Write the cross-dataset obs dictionary to csv files + for pair in data_dict['cross-dataset-obs'].keys(): + data_dict['cross-dataset-obs'][pair] = data_dict['cross-dataset-obs'][pair] + write_dictionary_to_csv(data_dict['cross-dataset-obs'][pair], f"{pair}", cfg) + # TEMPORARY: Print the dictionary to the log for checking logger.debug("Final data dictionary:\n%s", data_dict) From 65ff7a7945e4b787b44dc223351ba4533dfd23f0 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Fri, 21 Nov 2025 12:16:55 +0000 Subject: [PATCH 07/26] fixing keys --- .../diag_scripts/seaice/seaice_sensitivity.py | 39 ++++++++++--------- 1 file changed, 21 insertions(+), 18 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 6f586c6df9..1a729a0ab1 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -42,6 +42,7 @@ def create_category_dict(cfg): 'models': {}, 'tasa_obs': {}, 'siconc_obs': {}, + 'cross-dataset-obs': {}, } # Initialize the models as a set to avoid duplication @@ -174,6 +175,7 @@ def write_values_to_dict(data_dict, cfg): cross_sensitivity = calculate_cross_dataset_stats(tasa_dataset, siconc_dataset, cfg) # Determine structure of dictionary to store values key_name = f"{siconc_dataset}_to_{tasa_dataset}" + data_dict['cross-dataset-obs'][key_name] = {} inner_dict = data_dict['cross-dataset-obs'][key_name] # Store the values inner_dict['direct_sensitivity'] = cross_sensitivity.slope @@ -207,24 +209,25 @@ def write_dictionary_to_csv(sub_dict, filename, cfg): def main(cfg): - # Create the structured dictionary - data_dict = create_category_dict(cfg) - - # Calculate and write values to the dictionary - data_dict = write_values_to_dict(data_dict, cfg) - - # Write the model and obs dictionaries to csv files - write_dictionary_to_csv(data_dict['models'], 'models_values', cfg) - write_dictionary_to_csv(data_dict['tasa_obs'], 'tasa_obs_values', cfg) - write_dictionary_to_csv(data_dict['siconc_obs'], 'siconc_obs_values', cfg) - - # Write the cross-dataset obs dictionary to csv files - for pair in data_dict['cross-dataset-obs'].keys(): - data_dict['cross-dataset-obs'][pair] = data_dict['cross-dataset-obs'][pair] - write_dictionary_to_csv(data_dict['cross-dataset-obs'][pair], f"{pair}", cfg) - - # TEMPORARY: Print the dictionary to the log for checking - logger.debug("Final data dictionary:\n%s", data_dict) + print(cfg) + # # Create the structured dictionary + # data_dict = create_category_dict(cfg) + # + # # Calculate and write values to the dictionary + # data_dict = write_values_to_dict(data_dict, cfg) + # + # # Write the model and obs dictionaries to csv files + # write_dictionary_to_csv(data_dict['models'], 'models_values', cfg) + # write_dictionary_to_csv(data_dict['tasa_obs'], 'tasa_obs_values', cfg) + # write_dictionary_to_csv(data_dict['siconc_obs'], 'siconc_obs_values', cfg) + # + # # Write the cross-dataset obs dictionary to csv files + # for pair in data_dict['cross-dataset-obs'].keys(): + # data_dict['cross-dataset-obs'][pair] = data_dict['cross-dataset-obs'][pair] + # write_dictionary_to_csv(data_dict['cross-dataset-obs'][pair], f"{pair}", cfg) + # + # # TEMPORARY: Print the dictionary to the log for checking + # logger.debug("Final data dictionary:\n%s", data_dict) if __name__ == "__main__": with run_diagnostic() as config: From c6d6335ded5703c0f1a61acd14135c21b3ef2a1a Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Fri, 21 Nov 2025 12:18:20 +0000 Subject: [PATCH 08/26] temp --- .../diag_scripts/seaice/seaice_sensitivity.py | 38 +++++++++---------- 1 file changed, 19 insertions(+), 19 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 1a729a0ab1..bf6cfcd1fe 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -209,25 +209,25 @@ def write_dictionary_to_csv(sub_dict, filename, cfg): def main(cfg): - print(cfg) - # # Create the structured dictionary - # data_dict = create_category_dict(cfg) - # - # # Calculate and write values to the dictionary - # data_dict = write_values_to_dict(data_dict, cfg) - # - # # Write the model and obs dictionaries to csv files - # write_dictionary_to_csv(data_dict['models'], 'models_values', cfg) - # write_dictionary_to_csv(data_dict['tasa_obs'], 'tasa_obs_values', cfg) - # write_dictionary_to_csv(data_dict['siconc_obs'], 'siconc_obs_values', cfg) - # - # # Write the cross-dataset obs dictionary to csv files - # for pair in data_dict['cross-dataset-obs'].keys(): - # data_dict['cross-dataset-obs'][pair] = data_dict['cross-dataset-obs'][pair] - # write_dictionary_to_csv(data_dict['cross-dataset-obs'][pair], f"{pair}", cfg) - # - # # TEMPORARY: Print the dictionary to the log for checking - # logger.debug("Final data dictionary:\n%s", data_dict) + # Create the structured dictionary + data_dict = create_category_dict(cfg) + + # Calculate and write values to the dictionary + data_dict = write_values_to_dict(data_dict, cfg) + + # Write the model and obs dictionaries to csv files + write_dictionary_to_csv(data_dict['models'], 'models_values', cfg) + write_dictionary_to_csv(data_dict['tasa_obs'], 'tasa_obs_values', cfg) + write_dictionary_to_csv(data_dict['siconc_obs'], 'siconc_obs_values', cfg) + + # Write the cross-dataset obs dictionary to csv files + for pair in data_dict['cross-dataset-obs'].keys(): + data_dict['cross-dataset-obs'][pair] = data_dict['cross-dataset-obs'][pair] + write_dictionary_to_csv(data_dict['cross-dataset-obs'][pair], f"{pair}", cfg) + + # TEMPORARY: Print the dictionary to the log for checking + logger.debug("Final data dictionary:\n%s", data_dict) + if __name__ == "__main__": with run_diagnostic() as config: From 879bf918b7edb3b6dd342c90a0c064bebac31995 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Fri, 21 Nov 2025 13:47:38 +0000 Subject: [PATCH 09/26] adding one plot back in --- .../diag_scripts/seaice/seaice_sensitivity.py | 167 +++++++++++++++++- 1 file changed, 164 insertions(+), 3 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index bf6cfcd1fe..2c771dc7a2 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -12,6 +12,7 @@ from esmvaltool.diag_scripts.shared import ( ProvenanceLogger, run_diagnostic, + save_figure, ) logger = logging.getLogger(Path(__file__).stem) @@ -37,6 +38,7 @@ def get_provenance_record(cfg, caption): def create_category_dict(cfg): """Create a structured dictionary for adding values to later on.""" + logger.debug("Creating blank dictionary.") # Create blank dictionary of correct structure category_dict = { 'models': {}, @@ -75,6 +77,8 @@ def create_category_dict(cfg): def fetch_cube(dataset, variable, cfg): """Fetch a data cube for a dataset and variable using info from the config.""" + logger.debug("Fetching cube for dataset: %s, variable: %s", dataset, variable) + # Read the data from the config object input_data = cfg["input_data"].values() @@ -94,6 +98,8 @@ def fetch_cube(dataset, variable, cfg): def calculate_annual_trend(cube): """Calculate the linear trend of a cube over time using scipy.stats.linregres.""" + logger.debug("Calculating annual trend for cube %s.", cube.name()) + # Depending on preprocessor, coord may be 'year' or 'time' if 'year' in cube.coords(): no_years = [i for i in range(len(cube.coord("years").points))] @@ -109,6 +115,8 @@ def calculate_annual_trend(cube): def calculate_direct_stats(dataset, cfg): """Calculate the direct sensitivity of siconc to tas for a given dataset.""" + logger.debug("Calculating direct sensitivity for dataset %s.", dataset) + # Fetch the required cubes siconc_cube = fetch_cube(dataset, 'siconc', cfg) tas_cube = fetch_cube(dataset, 'tas', cfg) @@ -122,6 +130,8 @@ def calculate_direct_stats(dataset, cfg): def calculate_cross_dataset_stats(tasa_dataset, siconc_dataset, cfg): """Calculate the sensitivity of siconc to tasa across (obs) datasets.""" + logger.debug("Calculating cross sensitivity for datasets %s, %s.", tasa_dataset, siconc_dataset) + # Fetch the required cubes siconc_cube = fetch_cube(siconc_dataset, 'siconc', cfg) tasa_cube = fetch_cube(tasa_dataset, 'tasa', cfg) @@ -135,6 +145,8 @@ def calculate_cross_dataset_stats(tasa_dataset, siconc_dataset, cfg): def write_values_to_dict(data_dict, cfg): """Calculate and write values to the structured dictionary.""" + logger.debug("Writing values to dictionary.") + # Calculate all the values for the models for model_dataset in data_dict['models'].keys(): # Calculate annual tas trend @@ -192,6 +204,8 @@ def write_dictionary_to_csv(sub_dict, filename, cfg): Only sections of the dictionary should be written at a time as otherwise the structure is too complex to easily convert to a DataFrame. """ + logger.debug("Writing dictionary to csv file.") + # Create the csv filepath using info from the config csv_filepath = f"{cfg['work_dir']}/{filename}.csv" @@ -208,25 +222,172 @@ def write_dictionary_to_csv(sub_dict, filename, cfg): ) +def write_obs_from_cfg(cfg): + """Write the Notz-style observations from the recipe to a dictionary.""" + logger.debug("Writing observations from config file.") + + # Initialize the dictionary with observation period + obs_dict = {"obs_period": cfg["observations"]["observation period"]} + + # Add the observations values to the dictionary + notz_values = cfg["observations"]["sea ice sensitivity (Notz-style plot)"] + obs_dict["mean"] = notz_values["mean"] + obs_dict["std_dev"] = notz_values["standard deviation"] + obs_dict["plausible"] = notz_values["plausible range"] + + return obs_dict + + +def create_titles_dict(cfg): + """ + Create a dictionary of appropriate titles and hardcoded observations. + Values depend on whether the plot is for the Arctic or Antarctic + and assume the recipe used September Arctic sea ice data or + annually mean averaged Antarctic sea ice data + """ + logger.debug("Creating titles dictionary.") + dictionary = {} + + data = cfg["input_data"].values() + first_variable = next(iter(data)) + + if first_variable["diagnostic"] == "arctic": + dictionary["titles"] = { + "notz_fig_title": "September Arctic sea-ice area sensitivity to global mean surface temperature", + "notz_ax_title": "dSIA/dGMST", + "notz_plot_filename": "September Arctic sea ice sensitivity", + "roach_fig_title": "Trends in Annual Mean Temperature And September Arctic Sea Ice", + "roach_plot_filename": "September Arctic sea ice trends", + } + + elif first_variable["diagnostic"] == "antarctic": + dictionary["titles"] = { + "notz_fig_title": "Annually Meaned Antarctic Sea Ice Sensitivity", + "notz_ax_title": "dSIA/dGMST", + "notz_plot_filename": "Annual Antarctic sea ice sensitivity", + "roach_fig_title": "Trends in Annual Mean Temperature And Annual Antarctic Sea Ice", + "roach_plot_filename": "Annual Antarctic sea ice trends", + } + + dictionary["obs"] = write_obs_from_cfg(cfg) + + return dictionary + + +def notz_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): + """Save a plot of sensitivities and observations for model datasets.""" + # Read from observations dictionary + obs_years = titles_dictionary["obs"]["obs_period"] + obs_dict = titles_dictionary["obs"] + obs_mean = obs_dict["mean"] + obs_std_dev = obs_dict["std_dev"] + obs_plausible = obs_dict["plausible"] + + # Set up the figure + fig, ax = plt.subplots(figsize=(3.5, 6), layout="constrained") + fig.suptitle(titles_dictionary["titles"]["notz_fig_title"], wrap=True) + ax.set_title( + titles_dictionary["titles"]["notz_ax_title"], wrap=True, fontsize=10 + ) + + # Iterate over the dictionary + for dataset, inner_dict in data_dictionary.items(): + ax.plot( + 0.25, + inner_dict["direct_sensitivity_(notz-style)"], + color="blue", + marker="_", + markersize=20, + ) + + # Label with the dataset if specified + if inner_dict["label"] == "to_label": + plt.annotate( + dataset, + xy=(0.25, inner_dict["direct_sensitivity_(notz-style)"]), + xytext=( + 0.35, + inner_dict["direct_sensitivity_(notz-style)"] - 0.05, + ), + ) + + # Add observations only if obs_mean is a number + if isinstance(obs_mean, int | float): + ax.hlines(obs_mean, 0, 1, linestyle="--", color="black", linewidth=2) + ax.fill_between( + [0, 1], + obs_mean - obs_std_dev, + obs_mean + obs_std_dev, + facecolor="k", + alpha=0.15, + ) + ax.hlines( + obs_mean + obs_plausible, + 0, + 1, + linestyle=":", + color="0.5", + linewidth=1, + ) + ax.hlines( + obs_mean - obs_plausible, + 0, + 1, + linestyle=":", + color="0.5", + linewidth=1, + ) + + # Tidy the figure + ax.set_xlim(0, 1) + ax.set_xticks([]) + ax.set_ylabel(r"dSIA/dGMST ($million \ km^2 \ K^{-1}$)") + + # Create caption based on whether observation mean is presnt + if isinstance(obs_mean, int | float): + caption = ( + "Sensitivity of sea ice area to annual mean global warming." + f"Mean (dashed), standard deviation (shaded) and plausible values from {obs_years}." + ) + else: + caption = "Sensitivity of sea ice area to annual mean global warming." + + # Save the figure (also closes it) + save_figure( + titles_dictionary["titles"]["notz_plot_filename"], + get_provenance_record(cfg, caption), + cfg, + figure=fig, + close=True, + ) + + def main(cfg): # Create the structured dictionary data_dict = create_category_dict(cfg) # Calculate and write values to the dictionary + logger.info(f"Calculating and writing values to dictionary.") data_dict = write_values_to_dict(data_dict, cfg) # Write the model and obs dictionaries to csv files + logger.info(f"Writing dictionaries to csv files.") write_dictionary_to_csv(data_dict['models'], 'models_values', cfg) write_dictionary_to_csv(data_dict['tasa_obs'], 'tasa_obs_values', cfg) write_dictionary_to_csv(data_dict['siconc_obs'], 'siconc_obs_values', cfg) - # Write the cross-dataset obs dictionary to csv files + # Write the cross-dataset obs dictionary to csv files (separately for each pair) for pair in data_dict['cross-dataset-obs'].keys(): data_dict['cross-dataset-obs'][pair] = data_dict['cross-dataset-obs'][pair] write_dictionary_to_csv(data_dict['cross-dataset-obs'][pair], f"{pair}", cfg) - # TEMPORARY: Print the dictionary to the log for checking - logger.debug("Final data dictionary:\n%s", data_dict) + # Titles and observations depend on the diagnostic being plotted + logger.info("Creating titles and observations dictionary") + titles_and_obs_dict = create_titles_dict(cfg) + + # Plot the sensitivities (and save and close the plots) + logger.info("Creating Notz-style plot") + notz_style_plot_from_dict(data_dict, titles_and_obs_dict, cfg) if __name__ == "__main__": From 8deb66a698766655c15ad435d8d500135546db16 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Fri, 21 Nov 2025 14:04:31 +0000 Subject: [PATCH 10/26] error finding --- esmvaltool/diag_scripts/seaice/seaice_sensitivity.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 2c771dc7a2..2e45fe64bd 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -294,7 +294,7 @@ def notz_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): for dataset, inner_dict in data_dictionary.items(): ax.plot( 0.25, - inner_dict["direct_sensitivity_(notz-style)"], + inner_dict["direct_sensitivity"], color="blue", marker="_", markersize=20, @@ -304,10 +304,10 @@ def notz_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): if inner_dict["label"] == "to_label": plt.annotate( dataset, - xy=(0.25, inner_dict["direct_sensitivity_(notz-style)"]), + xy=(0.25, inner_dict["direct_sensitivity"]), xytext=( 0.35, - inner_dict["direct_sensitivity_(notz-style)"] - 0.05, + inner_dict["direct_sensitivity"] - 0.05, ), ) From aa1a3127292fd7e6f828c29bbb90cdf85631a078 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Fri, 21 Nov 2025 14:32:41 +0000 Subject: [PATCH 11/26] error finding --- .../diag_scripts/seaice/seaice_sensitivity.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 2e45fe64bd..9a0d0ccdda 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -47,9 +47,6 @@ def create_category_dict(cfg): 'cross-dataset-obs': {}, } - # Initialize the models as a set to avoid duplication - models_set = set() - # Read the data from the config object input_data = cfg["input_data"].values() @@ -66,11 +63,17 @@ def create_category_dict(cfg): # Everything else should be a model else: - models_set.add(input['dataset']) - - # Add the models to the dictionary - for model in models_set: - category_dict['models'][model] = {} + # Add the model dataset if not already present (appears twice, for tas and siconc) + if input['dataset'] not in category_dict['models']: + category_dict['models'][input['dataset']] = {} + + # Add labelling info + if 'label_dataset' in input and input['label_dataset']: + category_dict['models'][input['dataset']]['label'] = 'to_label' + logger.info("Dataset %s will be labelled", input['dataset']) + else: + category_dict['models'][input['dataset']]['label'] = 'unlabelled' + logger.info("Not labelling dataset %s in plots", input['dataset']) return category_dict From 3347c45112833c30ab0a35d1f52a211822e3500b Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Fri, 21 Nov 2025 14:38:12 +0000 Subject: [PATCH 12/26] error finding --- esmvaltool/diag_scripts/seaice/seaice_sensitivity.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 9a0d0ccdda..4438c2104c 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -388,9 +388,9 @@ def main(cfg): logger.info("Creating titles and observations dictionary") titles_and_obs_dict = create_titles_dict(cfg) - # Plot the sensitivities (and save and close the plots) + # Plot the sensitivities, uses model data only (and obs from recipe) logger.info("Creating Notz-style plot") - notz_style_plot_from_dict(data_dict, titles_and_obs_dict, cfg) + notz_style_plot_from_dict(data_dict['models'], titles_and_obs_dict, cfg) if __name__ == "__main__": From f7468daf385a6bfca1179e52ea400900e5f07f97 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Mon, 24 Nov 2025 08:58:37 +0000 Subject: [PATCH 13/26] retry partial replace --- .../diag_scripts/seaice/seaice_sensitivity.py | 60 +++++++++++++++++++ 1 file changed, 60 insertions(+) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 4438c2104c..29d4d56791 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -365,6 +365,64 @@ def notz_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): ) +def roach_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): + """Save a plot of trend in SIA against trend in GMST to the given filename.""" + # Set up the figure + fig, ax = plt.subplots(figsize=(10, 6), layout="constrained") + fig.suptitle(titles_dictionary["titles"]["roach_fig_title"], wrap=True) + + # Set up for colouring the points + norm = Normalize(vmin=-1, vmax=1) + cmap = plt.get_cmap("PiYG_r") + + # Set up the axes + ax.axhline(color="black", alpha=0.5) + ax.axvline(color="black", alpha=0.5) + ax.set_xlabel(r"Trend in GMST ($K \ decade^{-1}$)") + ax.set_ylabel(r"Trend in SIA ($million \ km^2 \ decade^{-1}$)") + + # Iterate over the models sub-dictionary + for dataset, inner_dict in data_dictionary['models'].items(): + # Determine the position of the point + x = 10 * inner_dict["annual_tas_trend"] # for equivalence to decades + y = ( + 10 * inner_dict["annual_siconc_trend"] + ) # for equivalence to decades + + # Determine the colour of the point + r_corr = inner_dict["direct_r_val"] + + # Decide if the point should be hatched + if inner_dict["direct_p_val"] >= 0.05: + h = 5 * "/" # This is a hatch pattern + else: + h = None + + # Plot the point + plt.scatter( + x, y, marker="o", s=150, c=[r_corr], hatch=h, cmap=cmap, norm=norm + ) + + # Label with the dataset if specified + if inner_dict["label"] == "to_label": + plt.annotate(dataset, xy=(x, y), xytext=(x + 0.01, y - 0.005)) + + # TODO: Obs here + + # Add a colour bar + plt.colorbar(label="Pearson correlation coefficient") + + # Save the figure (also closes it) + caption = "Decadal trends of sea ice area and global mean temperature." + save_figure( + titles_dictionary["titles"]["roach_plot_filename"], + get_provenance_record(cfg, caption), + cfg, + figure=fig, + close=True, + ) + + def main(cfg): # Create the structured dictionary data_dict = create_category_dict(cfg) @@ -391,6 +449,8 @@ def main(cfg): # Plot the sensitivities, uses model data only (and obs from recipe) logger.info("Creating Notz-style plot") notz_style_plot_from_dict(data_dict['models'], titles_and_obs_dict, cfg) + logger.info("Creating Roach-style plot") + roach_style_plot_from_dict(data_dict, titles_and_obs_dict, cfg) if __name__ == "__main__": From 1658389fe1e6097fd126b45447f221cc4004558f Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Mon, 24 Nov 2025 09:06:15 +0000 Subject: [PATCH 14/26] error fixing --- esmvaltool/diag_scripts/seaice/seaice_sensitivity.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 29d4d56791..b2e68e2c83 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -390,10 +390,10 @@ def roach_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): ) # for equivalence to decades # Determine the colour of the point - r_corr = inner_dict["direct_r_val"] + r_corr = inner_dict["direct_r_value"] # Decide if the point should be hatched - if inner_dict["direct_p_val"] >= 0.05: + if inner_dict["direct_p_value"] >= 0.05: h = 5 * "/" # This is a hatch pattern else: h = None From 293720d6d87d1ac91597b05fecb5b99909cb48e1 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Mon, 24 Nov 2025 09:20:32 +0000 Subject: [PATCH 15/26] error with provenance --- esmvaltool/diag_scripts/seaice/seaice_sensitivity.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index b2e68e2c83..0a35a4ac9b 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -414,9 +414,16 @@ def roach_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): # Save the figure (also closes it) caption = "Decadal trends of sea ice area and global mean temperature." + # save_figure( + # titles_dictionary["titles"]["roach_plot_filename"], + # get_provenance_record(cfg, caption), + # cfg, + # figure=fig, + # close=True, + # ) save_figure( titles_dictionary["titles"]["roach_plot_filename"], - get_provenance_record(cfg, caption), + {}, cfg, figure=fig, close=True, From a9a13a68d9a5667519f3b745c9d0bb8f4f7c7dc4 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Mon, 24 Nov 2025 09:39:32 +0000 Subject: [PATCH 16/26] REMOVING provenance --- .../diag_scripts/seaice/seaice_sensitivity.py | 23 ++++++++++++------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 0a35a4ac9b..410238e775 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -217,13 +217,6 @@ def write_dictionary_to_csv(sub_dict, filename, cfg): dataframe.to_csv(csv_filepath) logger.info("Wrote data to %s", csv_filepath) - # Create a provenance record for the csv file - with ProvenanceLogger(cfg) as provenance_logger: - provenance_logger.log( - csv_filepath, - get_provenance_record(cfg, "Annual (not decadal) figures"), - ) - def write_obs_from_cfg(cfg): """Write the Notz-style observations from the recipe to a dictionary.""" @@ -356,9 +349,16 @@ def notz_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): caption = "Sensitivity of sea ice area to annual mean global warming." # Save the figure (also closes it) + # save_figure( + # titles_dictionary["titles"]["notz_plot_filename"], + # get_provenance_record(cfg, caption), + # cfg, + # figure=fig, + # close=True, + # ) save_figure( titles_dictionary["titles"]["notz_plot_filename"], - get_provenance_record(cfg, caption), + {}, cfg, figure=fig, close=True, @@ -444,6 +444,13 @@ def main(cfg): write_dictionary_to_csv(data_dict['tasa_obs'], 'tasa_obs_values', cfg) write_dictionary_to_csv(data_dict['siconc_obs'], 'siconc_obs_values', cfg) + # Create a provenance record for the csv files + with ProvenanceLogger(cfg) as provenance_logger: + provenance_logger.log( + f"{cfg['work_dir']}/figures_as_csv", + get_provenance_record(cfg, "Annual (not decadal) figures"), + ) + # Write the cross-dataset obs dictionary to csv files (separately for each pair) for pair in data_dict['cross-dataset-obs'].keys(): data_dict['cross-dataset-obs'][pair] = data_dict['cross-dataset-obs'][pair] From cc4301eae1a13109c4be5c1bd03642de7fc6a052 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Mon, 24 Nov 2025 09:54:48 +0000 Subject: [PATCH 17/26] error with provenance --- .../diag_scripts/seaice/seaice_sensitivity.py | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 410238e775..c73a88a35e 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -349,16 +349,9 @@ def notz_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): caption = "Sensitivity of sea ice area to annual mean global warming." # Save the figure (also closes it) - # save_figure( - # titles_dictionary["titles"]["notz_plot_filename"], - # get_provenance_record(cfg, caption), - # cfg, - # figure=fig, - # close=True, - # ) save_figure( titles_dictionary["titles"]["notz_plot_filename"], - {}, + get_provenance_record(cfg, caption), cfg, figure=fig, close=True, @@ -421,13 +414,11 @@ def roach_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): # figure=fig, # close=True, # ) - save_figure( - titles_dictionary["titles"]["roach_plot_filename"], - {}, - cfg, - figure=fig, - close=True, + # TODO: Temporary save to try to get around provenance error + fig.savefig( + f"{cfg['plot_dir']}/{titles_dictionary['titles']['roach_plot_filename']}.png" ) + plt.close(fig) def main(cfg): From 17041049bb8123d9fe1141cf86eb02ae01ffc470 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Mon, 24 Nov 2025 09:55:08 +0000 Subject: [PATCH 18/26] error with provenance --- esmvaltool/diag_scripts/seaice/seaice_sensitivity.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index c73a88a35e..60d6895d1c 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -435,7 +435,7 @@ def main(cfg): write_dictionary_to_csv(data_dict['tasa_obs'], 'tasa_obs_values', cfg) write_dictionary_to_csv(data_dict['siconc_obs'], 'siconc_obs_values', cfg) - # Create a provenance record for the csv files + # Create a single provenance record for the csv files with ProvenanceLogger(cfg) as provenance_logger: provenance_logger.log( f"{cfg['work_dir']}/figures_as_csv", From 7fc87d2beadcca344bd79db0df25d644416d8b77 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Mon, 24 Nov 2025 09:57:59 +0000 Subject: [PATCH 19/26] error with provenance --- esmvaltool/diag_scripts/seaice/seaice_sensitivity.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 60d6895d1c..888bfcdb22 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -435,6 +435,11 @@ def main(cfg): write_dictionary_to_csv(data_dict['tasa_obs'], 'tasa_obs_values', cfg) write_dictionary_to_csv(data_dict['siconc_obs'], 'siconc_obs_values', cfg) + # Write the cross-dataset obs dictionary to csv files (separately for each pair) + for pair in data_dict['cross-dataset-obs'].keys(): + data_dict['cross-dataset-obs'][pair] = data_dict['cross-dataset-obs'][pair] + write_dictionary_to_csv(data_dict['cross-dataset-obs'][pair], f"{pair}", cfg) + # Create a single provenance record for the csv files with ProvenanceLogger(cfg) as provenance_logger: provenance_logger.log( @@ -442,11 +447,6 @@ def main(cfg): get_provenance_record(cfg, "Annual (not decadal) figures"), ) - # Write the cross-dataset obs dictionary to csv files (separately for each pair) - for pair in data_dict['cross-dataset-obs'].keys(): - data_dict['cross-dataset-obs'][pair] = data_dict['cross-dataset-obs'][pair] - write_dictionary_to_csv(data_dict['cross-dataset-obs'][pair], f"{pair}", cfg) - # Titles and observations depend on the diagnostic being plotted logger.info("Creating titles and observations dictionary") titles_and_obs_dict = create_titles_dict(cfg) From 2f1df1b8ad861e484a8edd61f845526a629bfd09 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Mon, 24 Nov 2025 10:07:58 +0000 Subject: [PATCH 20/26] nope, error was incorrect indentation --- .../diag_scripts/seaice/seaice_sensitivity.py | 33 ++++++++----------- 1 file changed, 14 insertions(+), 19 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 888bfcdb22..05fd8d128c 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -392,33 +392,28 @@ def roach_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): h = None # Plot the point - plt.scatter( + ax.scatter( x, y, marker="o", s=150, c=[r_corr], hatch=h, cmap=cmap, norm=norm ) # Label with the dataset if specified if inner_dict["label"] == "to_label": - plt.annotate(dataset, xy=(x, y), xytext=(x + 0.01, y - 0.005)) + ax.annotate(dataset, xy=(x, y), xytext=(x + 0.01, y - 0.005)) # TODO: Obs here - # Add a colour bar - plt.colorbar(label="Pearson correlation coefficient") - - # Save the figure (also closes it) - caption = "Decadal trends of sea ice area and global mean temperature." - # save_figure( - # titles_dictionary["titles"]["roach_plot_filename"], - # get_provenance_record(cfg, caption), - # cfg, - # figure=fig, - # close=True, - # ) - # TODO: Temporary save to try to get around provenance error - fig.savefig( - f"{cfg['plot_dir']}/{titles_dictionary['titles']['roach_plot_filename']}.png" - ) - plt.close(fig) + # Add a colour bar + fig.colorbar(label="Pearson correlation coefficient") + + # Save the figure (also closes it) + caption = "Decadal trends of sea ice area and global mean temperature." + save_figure( + titles_dictionary["titles"]["roach_plot_filename"], + get_provenance_record(cfg, caption), + cfg, + figure=fig, + close=True, + ) def main(cfg): From 2be7b2741aa57302c064031d724edfe8daef1c55 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Mon, 24 Nov 2025 10:14:47 +0000 Subject: [PATCH 21/26] revert to less good plt for colourbar --- esmvaltool/diag_scripts/seaice/seaice_sensitivity.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 05fd8d128c..4ca75836bc 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -392,18 +392,18 @@ def roach_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): h = None # Plot the point - ax.scatter( + plt.scatter( x, y, marker="o", s=150, c=[r_corr], hatch=h, cmap=cmap, norm=norm ) # Label with the dataset if specified if inner_dict["label"] == "to_label": - ax.annotate(dataset, xy=(x, y), xytext=(x + 0.01, y - 0.005)) + plt.annotate(dataset, xy=(x, y), xytext=(x + 0.01, y - 0.005)) # TODO: Obs here # Add a colour bar - fig.colorbar(label="Pearson correlation coefficient") + plt.colorbar(label="Pearson correlation coefficient") # Save the figure (also closes it) caption = "Decadal trends of sea ice area and global mean temperature." From 106c8d52d5e28bd792ec269ee78107a8af3b3075 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Mon, 24 Nov 2025 10:34:38 +0000 Subject: [PATCH 22/26] add in obs --- .../diag_scripts/seaice/seaice_sensitivity.py | 27 ++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 4ca75836bc..d5879fa9cb 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -400,7 +400,32 @@ def roach_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): if inner_dict["label"] == "to_label": plt.annotate(dataset, xy=(x, y), xytext=(x + 0.01, y - 0.005)) - # TODO: Obs here + # Add the observations + siconc_dict = data_dictionary['siconc_obs'] + tasa_dict = data_dictionary['tasa_obs'] + + # Iterate over the pairs in cross-dataset-obs + for pair, inner_dict in data_dictionary['cross-dataset-obs'].items(): + # Retrieve the names of the datasets from the pair string + siconc_ds, tasa_ds = pair.split("_to_") + + # Determine the position of the point, from other dictionaries + x = 10 * tasa_dict[tasa_ds]['annual_tas_trend'] # This was labelled as tas, not tasa + y = 10 * siconc_dict[siconc_ds]['annual_siconc_trend'] + + # Determine the colour of the point from the inner dictionary + r_corr = inner_dict["direct_r_value"] + + # Decide if the point should be hatched + if inner_dict["direct_p_value"] >= 0.05: + h = 5 * "/" + else: + h = None + + # Plot the point + plt.scatter( + x, y, marker="s", s=150, c=[r_corr], hatch=h, cmap=cmap, norm=norm, zorder=0, edgecolors="black", + ) # Add a colour bar plt.colorbar(label="Pearson correlation coefficient") From aa74676c21106e7646139157f9d294e8b4117f4e Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Mon, 24 Nov 2025 10:43:55 +0000 Subject: [PATCH 23/26] looking for comparability to paper --- esmvaltool/recipes/recipe_seaice_sensitivity.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/esmvaltool/recipes/recipe_seaice_sensitivity.yml b/esmvaltool/recipes/recipe_seaice_sensitivity.yml index c40fd6f7ea..38560e2cc9 100644 --- a/esmvaltool/recipes/recipe_seaice_sensitivity.yml +++ b/esmvaltool/recipes/recipe_seaice_sensitivity.yml @@ -42,8 +42,8 @@ model_datasets: &model_datasets obs_defaults: &obs_defaults { project: OBS, observation: True, tier: 2} tasa_obs: &tasa_obs - - { <<: *obs_defaults, dataset: GISTEMP, type: ground, version: v4 } # Doesn't match Roach obs - - { <<: *obs_defaults, dataset: HadCRUT5, type: ground, version: 5.0.1.0-analysis } # Doesn't match Roach obs +# - { <<: *obs_defaults, dataset: GISTEMP, type: ground, version: v4 } # Doesn't match Roach obs +# - { <<: *obs_defaults, dataset: HadCRUT5, type: ground, version: 5.0.1.0-analysis } # Doesn't match Roach obs - { <<: *obs_defaults, dataset: NOAAGlobalTemp, type: ground, version: v5.0.0, mip: Amon } arctic_siconc_obs: &arctic_siconc_obs @@ -61,7 +61,7 @@ preprocessors: start_year: 1979 end_day: 31 end_month: 12 - end_year: 2014 + end_year: 2018 # Put this back to 2014 after the test extract_sept: &extract_sept From 76c116a7f60d748267a879df6a1a4703375cdb41 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Mon, 24 Nov 2025 13:30:50 +0000 Subject: [PATCH 24/26] pre-commit changes --- .../diag_scripts/seaice/seaice_sensitivity.py | 164 +++++++++++------- .../recipes/recipe_seaice_sensitivity.yml | 6 +- 2 files changed, 103 insertions(+), 67 deletions(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index d5879fa9cb..280308f501 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -41,10 +41,10 @@ def create_category_dict(cfg): logger.debug("Creating blank dictionary.") # Create blank dictionary of correct structure category_dict = { - 'models': {}, - 'tasa_obs': {}, - 'siconc_obs': {}, - 'cross-dataset-obs': {}, + "models": {}, + "tasa_obs": {}, + "siconc_obs": {}, + "cross-dataset-obs": {}, } # Read the data from the config object @@ -52,46 +52,51 @@ def create_category_dict(cfg): # Iterate over the datasets to add to the dictionary for input in input_data: - # Check for tasa observations - if input['variable_group'] == 'tasa_obs': - category_dict['tasa_obs'][input['dataset']] = {} + if input["variable_group"] == "tasa_obs": + category_dict["tasa_obs"][input["dataset"]] = {} # Check for siconc observations - elif input['variable_group'] == 'siconc_obs': - category_dict['siconc_obs'][input['dataset']] = {} + elif input["variable_group"] == "siconc_obs": + category_dict["siconc_obs"][input["dataset"]] = {} # Everything else should be a model else: # Add the model dataset if not already present (appears twice, for tas and siconc) - if input['dataset'] not in category_dict['models']: - category_dict['models'][input['dataset']] = {} + if input["dataset"] not in category_dict["models"]: + category_dict["models"][input["dataset"]] = {} # Add labelling info - if 'label_dataset' in input and input['label_dataset']: - category_dict['models'][input['dataset']]['label'] = 'to_label' - logger.info("Dataset %s will be labelled", input['dataset']) + if "label_dataset" in input and input["label_dataset"]: + category_dict["models"][input["dataset"]]["label"] = "to_label" + logger.info("Dataset %s will be labelled", input["dataset"]) else: - category_dict['models'][input['dataset']]['label'] = 'unlabelled' - logger.info("Not labelling dataset %s in plots", input['dataset']) + category_dict["models"][input["dataset"]]["label"] = ( + "unlabelled" + ) + logger.info( + "Not labelling dataset %s in plots", input["dataset"] + ) return category_dict def fetch_cube(dataset, variable, cfg): """Fetch a data cube for a dataset and variable using info from the config.""" - logger.debug("Fetching cube for dataset: %s, variable: %s", dataset, variable) + logger.debug( + "Fetching cube for dataset: %s, variable: %s", dataset, variable + ) # Read the data from the config object input_data = cfg["input_data"].values() # Find the correct filepath for the dataset for input in input_data: - if input['dataset'] == dataset: + if input["dataset"] == dataset: # Also check the variable matches as models have two entries # Only matching the first three letters to avoid issues with sic vs siconc - if input['short_name'][:3] == variable[:3]: - filepath = input['filename'] + if input["short_name"][:3] == variable[:3]: + filepath = input["filename"] break # Load the cube using iris @@ -104,7 +109,7 @@ def calculate_annual_trend(cube): logger.debug("Calculating annual trend for cube %s.", cube.name()) # Depending on preprocessor, coord may be 'year' or 'time' - if 'year' in cube.coords(): + if "year" in cube.coords(): no_years = [i for i in range(len(cube.coord("years").points))] else: no_years = [i for i in range(len(cube.coord("time").points))] @@ -121,8 +126,8 @@ def calculate_direct_stats(dataset, cfg): logger.debug("Calculating direct sensitivity for dataset %s.", dataset) # Fetch the required cubes - siconc_cube = fetch_cube(dataset, 'siconc', cfg) - tas_cube = fetch_cube(dataset, 'tas', cfg) + siconc_cube = fetch_cube(dataset, "siconc", cfg) + tas_cube = fetch_cube(dataset, "tas", cfg) # Calculate direct regression (tas as independent) direct_sensitivity = linregress(tas_cube.data, siconc_cube.data) @@ -133,11 +138,15 @@ def calculate_direct_stats(dataset, cfg): def calculate_cross_dataset_stats(tasa_dataset, siconc_dataset, cfg): """Calculate the sensitivity of siconc to tasa across (obs) datasets.""" - logger.debug("Calculating cross sensitivity for datasets %s, %s.", tasa_dataset, siconc_dataset) + logger.debug( + "Calculating cross sensitivity for datasets %s, %s.", + tasa_dataset, + siconc_dataset, + ) # Fetch the required cubes - siconc_cube = fetch_cube(siconc_dataset, 'siconc', cfg) - tasa_cube = fetch_cube(tasa_dataset, 'tasa', cfg) + siconc_cube = fetch_cube(siconc_dataset, "siconc", cfg) + tasa_cube = fetch_cube(tasa_dataset, "tasa", cfg) # Calculate direct regression (tasa as independent) direct_sensitivity = linregress(tasa_cube.data, siconc_cube.data) @@ -151,51 +160,63 @@ def write_values_to_dict(data_dict, cfg): logger.debug("Writing values to dictionary.") # Calculate all the values for the models - for model_dataset in data_dict['models'].keys(): + for model_dataset in data_dict["models"].keys(): # Calculate annual tas trend - tas_cube = fetch_cube(model_dataset, 'tas', cfg) + tas_cube = fetch_cube(model_dataset, "tas", cfg) ann_tas_trend = calculate_annual_trend(tas_cube) - data_dict['models'][model_dataset]['annual_tas_trend'] = ann_tas_trend + data_dict["models"][model_dataset]["annual_tas_trend"] = ann_tas_trend # Calculate annual siconc trend - siconc_cube = fetch_cube(model_dataset, 'siconc', cfg) + siconc_cube = fetch_cube(model_dataset, "siconc", cfg) ann_siconc_trend = calculate_annual_trend(siconc_cube) - data_dict['models'][model_dataset]['annual_siconc_trend'] = ann_siconc_trend + data_dict["models"][model_dataset]["annual_siconc_trend"] = ( + ann_siconc_trend + ) # Calculate direct sensitivity of siconc to tas direct_sensitivity = calculate_direct_stats(model_dataset, cfg) - data_dict['models'][model_dataset]['direct_sensitivity'] = direct_sensitivity.slope - data_dict['models'][model_dataset]['direct_r_value'] = direct_sensitivity.rvalue - data_dict['models'][model_dataset]['direct_p_value'] = direct_sensitivity.pvalue + data_dict["models"][model_dataset]["direct_sensitivity"] = ( + direct_sensitivity.slope + ) + data_dict["models"][model_dataset]["direct_r_value"] = ( + direct_sensitivity.rvalue + ) + data_dict["models"][model_dataset]["direct_p_value"] = ( + direct_sensitivity.pvalue + ) # Calculate just the tasa trend for the tasa observations - for obs_dataset in data_dict['tasa_obs'].keys(): + for obs_dataset in data_dict["tasa_obs"].keys(): # Calculate annual tas trend - tasa_cube = fetch_cube(obs_dataset, 'tasa', cfg) + tasa_cube = fetch_cube(obs_dataset, "tasa", cfg) ann_tasa_trend = calculate_annual_trend(tasa_cube) # Still labelling in final dictionary as tas for consistency with models - data_dict['tasa_obs'][obs_dataset]['annual_tas_trend'] = ann_tasa_trend + data_dict["tasa_obs"][obs_dataset]["annual_tas_trend"] = ann_tasa_trend # Calculate just the siconc trend for the siconc observations - for obs_dataset in data_dict['siconc_obs'].keys(): + for obs_dataset in data_dict["siconc_obs"].keys(): # Calculate annual siconc trend - siconc_cube = fetch_cube(obs_dataset, 'siconc', cfg) + siconc_cube = fetch_cube(obs_dataset, "siconc", cfg) ann_siconc_trend = calculate_annual_trend(siconc_cube) - data_dict['siconc_obs'][obs_dataset]['annual_siconc_trend'] = ann_siconc_trend + data_dict["siconc_obs"][obs_dataset]["annual_siconc_trend"] = ( + ann_siconc_trend + ) # Calculate cross-dataset statistics between tasa and siconc observations - for tasa_dataset in data_dict['tasa_obs'].keys(): - for siconc_dataset in data_dict['siconc_obs'].keys(): + for tasa_dataset in data_dict["tasa_obs"].keys(): + for siconc_dataset in data_dict["siconc_obs"].keys(): # Calculate cross-dataset sensitivity of siconc to tasa - cross_sensitivity = calculate_cross_dataset_stats(tasa_dataset, siconc_dataset, cfg) + cross_sensitivity = calculate_cross_dataset_stats( + tasa_dataset, siconc_dataset, cfg + ) # Determine structure of dictionary to store values key_name = f"{siconc_dataset}_to_{tasa_dataset}" - data_dict['cross-dataset-obs'][key_name] = {} - inner_dict = data_dict['cross-dataset-obs'][key_name] + data_dict["cross-dataset-obs"][key_name] = {} + inner_dict = data_dict["cross-dataset-obs"][key_name] # Store the values - inner_dict['direct_sensitivity'] = cross_sensitivity.slope - inner_dict['direct_r_value'] = cross_sensitivity.rvalue - inner_dict['direct_p_value'] = cross_sensitivity.pvalue + inner_dict["direct_sensitivity"] = cross_sensitivity.slope + inner_dict["direct_r_value"] = cross_sensitivity.rvalue + inner_dict["direct_p_value"] = cross_sensitivity.pvalue return data_dict @@ -375,7 +396,7 @@ def roach_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): ax.set_ylabel(r"Trend in SIA ($million \ km^2 \ decade^{-1}$)") # Iterate over the models sub-dictionary - for dataset, inner_dict in data_dictionary['models'].items(): + for dataset, inner_dict in data_dictionary["models"].items(): # Determine the position of the point x = 10 * inner_dict["annual_tas_trend"] # for equivalence to decades y = ( @@ -401,17 +422,19 @@ def roach_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): plt.annotate(dataset, xy=(x, y), xytext=(x + 0.01, y - 0.005)) # Add the observations - siconc_dict = data_dictionary['siconc_obs'] - tasa_dict = data_dictionary['tasa_obs'] + siconc_dict = data_dictionary["siconc_obs"] + tasa_dict = data_dictionary["tasa_obs"] # Iterate over the pairs in cross-dataset-obs - for pair, inner_dict in data_dictionary['cross-dataset-obs'].items(): + for pair, inner_dict in data_dictionary["cross-dataset-obs"].items(): # Retrieve the names of the datasets from the pair string siconc_ds, tasa_ds = pair.split("_to_") # Determine the position of the point, from other dictionaries - x = 10 * tasa_dict[tasa_ds]['annual_tas_trend'] # This was labelled as tas, not tasa - y = 10 * siconc_dict[siconc_ds]['annual_siconc_trend'] + x = ( + 10 * tasa_dict[tasa_ds]["annual_tas_trend"] + ) # This was labelled as tas, not tasa + y = 10 * siconc_dict[siconc_ds]["annual_siconc_trend"] # Determine the colour of the point from the inner dictionary r_corr = inner_dict["direct_r_value"] @@ -424,7 +447,16 @@ def roach_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): # Plot the point plt.scatter( - x, y, marker="s", s=150, c=[r_corr], hatch=h, cmap=cmap, norm=norm, zorder=0, edgecolors="black", + x, + y, + marker="s", + s=150, + c=[r_corr], + hatch=h, + cmap=cmap, + norm=norm, + zorder=0, + edgecolors="black", ) # Add a colour bar @@ -446,19 +478,23 @@ def main(cfg): data_dict = create_category_dict(cfg) # Calculate and write values to the dictionary - logger.info(f"Calculating and writing values to dictionary.") + logger.info("Calculating and writing values to dictionary.") data_dict = write_values_to_dict(data_dict, cfg) # Write the model and obs dictionaries to csv files - logger.info(f"Writing dictionaries to csv files.") - write_dictionary_to_csv(data_dict['models'], 'models_values', cfg) - write_dictionary_to_csv(data_dict['tasa_obs'], 'tasa_obs_values', cfg) - write_dictionary_to_csv(data_dict['siconc_obs'], 'siconc_obs_values', cfg) + logger.info("Writing dictionaries to csv files.") + write_dictionary_to_csv(data_dict["models"], "models_values", cfg) + write_dictionary_to_csv(data_dict["tasa_obs"], "tasa_obs_values", cfg) + write_dictionary_to_csv(data_dict["siconc_obs"], "siconc_obs_values", cfg) # Write the cross-dataset obs dictionary to csv files (separately for each pair) - for pair in data_dict['cross-dataset-obs'].keys(): - data_dict['cross-dataset-obs'][pair] = data_dict['cross-dataset-obs'][pair] - write_dictionary_to_csv(data_dict['cross-dataset-obs'][pair], f"{pair}", cfg) + for pair in data_dict["cross-dataset-obs"].keys(): + data_dict["cross-dataset-obs"][pair] = data_dict["cross-dataset-obs"][ + pair + ] + write_dictionary_to_csv( + data_dict["cross-dataset-obs"][pair], f"{pair}", cfg + ) # Create a single provenance record for the csv files with ProvenanceLogger(cfg) as provenance_logger: @@ -473,7 +509,7 @@ def main(cfg): # Plot the sensitivities, uses model data only (and obs from recipe) logger.info("Creating Notz-style plot") - notz_style_plot_from_dict(data_dict['models'], titles_and_obs_dict, cfg) + notz_style_plot_from_dict(data_dict["models"], titles_and_obs_dict, cfg) logger.info("Creating Roach-style plot") roach_style_plot_from_dict(data_dict, titles_and_obs_dict, cfg) diff --git a/esmvaltool/recipes/recipe_seaice_sensitivity.yml b/esmvaltool/recipes/recipe_seaice_sensitivity.yml index 38560e2cc9..c40fd6f7ea 100644 --- a/esmvaltool/recipes/recipe_seaice_sensitivity.yml +++ b/esmvaltool/recipes/recipe_seaice_sensitivity.yml @@ -42,8 +42,8 @@ model_datasets: &model_datasets obs_defaults: &obs_defaults { project: OBS, observation: True, tier: 2} tasa_obs: &tasa_obs -# - { <<: *obs_defaults, dataset: GISTEMP, type: ground, version: v4 } # Doesn't match Roach obs -# - { <<: *obs_defaults, dataset: HadCRUT5, type: ground, version: 5.0.1.0-analysis } # Doesn't match Roach obs + - { <<: *obs_defaults, dataset: GISTEMP, type: ground, version: v4 } # Doesn't match Roach obs + - { <<: *obs_defaults, dataset: HadCRUT5, type: ground, version: 5.0.1.0-analysis } # Doesn't match Roach obs - { <<: *obs_defaults, dataset: NOAAGlobalTemp, type: ground, version: v5.0.0, mip: Amon } arctic_siconc_obs: &arctic_siconc_obs @@ -61,7 +61,7 @@ preprocessors: start_year: 1979 end_day: 31 end_month: 12 - end_year: 2018 # Put this back to 2014 after the test + end_year: 2014 extract_sept: &extract_sept From 445c3fd76bde7b3c07789162152b03407dc14adf Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Thu, 27 Nov 2025 15:48:43 +0000 Subject: [PATCH 25/26] Missing line break in caption. --- esmvaltool/diag_scripts/seaice/seaice_sensitivity.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index 280308f501..c8f23dfedf 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -364,7 +364,7 @@ def notz_style_plot_from_dict(data_dictionary, titles_dictionary, cfg): if isinstance(obs_mean, int | float): caption = ( "Sensitivity of sea ice area to annual mean global warming." - f"Mean (dashed), standard deviation (shaded) and plausible values from {obs_years}." + f"\nMean (dashed), standard deviation (shaded) and plausible values from {obs_years}." ) else: caption = "Sensitivity of sea ice area to annual mean global warming." From 578ffe6959930fed54a561b2d50341b1b577d1f2 Mon Sep 17 00:00:00 2001 From: Naomi Parsons Date: Fri, 28 Nov 2025 13:42:41 +0000 Subject: [PATCH 26/26] typo in docstring --- esmvaltool/diag_scripts/seaice/seaice_sensitivity.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py index c8f23dfedf..13417cf0ea 100644 --- a/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py +++ b/esmvaltool/diag_scripts/seaice/seaice_sensitivity.py @@ -105,7 +105,7 @@ def fetch_cube(dataset, variable, cfg): def calculate_annual_trend(cube): - """Calculate the linear trend of a cube over time using scipy.stats.linregres.""" + """Calculate the linear trend of a cube over time using scipy.stats.linregress""" logger.debug("Calculating annual trend for cube %s.", cube.name()) # Depending on preprocessor, coord may be 'year' or 'time'