From 16ba4a2a996a404a80f5dedacb6811e91c49331a Mon Sep 17 00:00:00 2001 From: Franziska Winterstein Date: Thu, 28 May 2026 17:03:22 +0200 Subject: [PATCH 01/33] implemented the option to test for significance when comparing with a reference --- .../config/configurations/defaults/extra_facets_emac.yml | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/esmvalcore/config/configurations/defaults/extra_facets_emac.yml b/esmvalcore/config/configurations/defaults/extra_facets_emac.yml index 169a69860f..8ac3d9d05d 100644 --- a/esmvalcore/config/configurations/defaults/extra_facets_emac.yml +++ b/esmvalcore/config/configurations/defaults/extra_facets_emac.yml @@ -257,9 +257,12 @@ projects: ta: raw_name: [tm1_cav, tm1_ave, tm1] channel: Amon + # tro3: + # raw_name: [O3_cav, O3_ave, O3] + # channel: Amon tro3: - raw_name: [O3_cav, O3_ave, O3] - channel: Amon + raw_name: [O3_p39_cav, O3_cav, O3_ave, O3] + channel: AERmon ua: raw_name: [um1_cav, um1_ave, um1] channel: Amon From 1f407788c833f239c6d073aca4d9c5176f580d65 Mon Sep 17 00:00:00 2001 From: Franziska Winterstein Date: Thu, 11 Jun 2026 14:06:13 +0200 Subject: [PATCH 02/33] correct attaching of ancillary variable --- esmvalcore/preprocessor/_compare_with_refs.py | 36 +++++++++++++++++-- 1 file changed, 33 insertions(+), 3 deletions(-) diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index 4e6fc6c180..b812ce754c 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -8,10 +8,12 @@ import dask import dask.array as da +import dask.array.stats import iris.analysis import iris.analysis.stats import numpy as np from iris.common.metadata import CubeMetadata +from iris.coords import AncillaryVariable from iris.coords import CellMethod from iris.cube import Cube, CubeList from scipy.stats import wasserstein_distance @@ -47,6 +49,7 @@ def bias( products: set[PreprocessorFile] | Iterable[Cube], reference: Cube | None = None, bias_type: BiasType = "absolute", + ttest: bool = False, denominator_mask_threshold: float = 1e-3, keep_reference_dataset: bool = False, ) -> set[PreprocessorFile] | CubeList: @@ -81,6 +84,11 @@ def bias( bias_type: Bias type that is calculated. Must be one of ``'absolute'`` (dataset - ref) or ``'relative'`` ((dataset - ref) / ref). + ttest: + If ``True`` calculate the p-value according to the Welch's Test between + dataset and reference. The p-value is attached as an ancillary variable + to the dataset. A Welch's Test is chosen, since unequal variances and + sample sizes are likely. denominator_mask_threshold: Threshold to mask values close to zero in the denominator (i.e., the reference dataset) during the calculation of relative biases. All @@ -138,7 +146,7 @@ def bias( # If input is an Iterable of Cube objects, calculate bias for each element if all_cubes_given: - cubes = [_calculate_bias(c, reference, bias_type) for c in products] + cubes = [_calculate_bias(c, reference, bias_type, ttest) for c in products] return CubeList(cubes) # Otherwise, iterate over all input products, calculate bias and adapt @@ -150,7 +158,7 @@ def bias( cube = concatenate(product.cubes) # Calculate bias - cube = _calculate_bias(cube, reference, bias_type) + cube = _calculate_bias(cube, reference, bias_type, ttest) # Adapt metadata and provenance information product.attributes["units"] = str(cube.units) @@ -196,10 +204,13 @@ def _get_ref( return (reference, ref_product) -def _calculate_bias(cube: Cube, reference: Cube, bias_type: BiasType) -> Cube: +def _calculate_bias(cube: Cube, reference: Cube, bias_type: BiasType, ttest: bool) -> Cube: """Calculate bias for a single cube relative to a reference cube.""" cube_metadata = cube.metadata + if ttest: + p_value, remaining_axis = _calculate_ttest(cube, reference) + if bias_type == "absolute": cube = cube - reference new_units = cube.units @@ -215,10 +226,29 @@ def _calculate_bias(cube: Cube, reference: Cube, bias_type: BiasType) -> Cube: cube.metadata = cube_metadata cube.units = new_units + # Attach ancillary variable after the calculation + # otherwise, it will be lost + if ttest: + cube.add_ancillary_variable( + AncillaryVariable(p_value, + long_name="pvalue"), + remaining_axis) return cube +def _calculate_ttest(cube: Cube, reference: Cube) -> Cube: + """Calculate the Welch's t-test and attach p_value as an ancillary variable.""" + axis = cube.coord_dims("time")[0] + remaining_axis = sorted(set(range(cube.ndim)) - set([axis])) + + ttest = dask.array.stats.ttest_ind(cube, reference, 0, False) + shape = tuple([list(cube.shape)[dim] for dim in remaining_axis]) + p_value = da.from_delayed(ttest.pvalue, shape, dtype=float) + + return p_value, remaining_axis + + MetricType = Literal[ "rmse", "weighted_rmse", From e22937c8ea5866340f77848ed07b0f3330415267 Mon Sep 17 00:00:00 2001 From: Franziska Winterstein Date: Fri, 12 Jun 2026 11:34:06 +0200 Subject: [PATCH 03/33] prepared ttest as its own preprocessor function, still open how to keep ancillary variables attached through arbitrary preprocessor functions --- esmvalcore/preprocessor/__init__.py | 4 +- esmvalcore/preprocessor/_compare_with_refs.py | 151 ++++++++++++++++-- 2 files changed, 141 insertions(+), 14 deletions(-) diff --git a/esmvalcore/preprocessor/__init__.py b/esmvalcore/preprocessor/__init__.py index 798edb36a8..4a16f06207 100644 --- a/esmvalcore/preprocessor/__init__.py +++ b/esmvalcore/preprocessor/__init__.py @@ -31,7 +31,7 @@ meridional_statistics, zonal_statistics, ) -from esmvalcore.preprocessor._compare_with_refs import bias, distance_metric +from esmvalcore.preprocessor._compare_with_refs import bias, distance_metric, ttest from esmvalcore.preprocessor._concatenate import concatenate from esmvalcore.preprocessor._cycles import amplitude from esmvalcore.preprocessor._dask_progress import _compute_with_progress @@ -219,6 +219,7 @@ # Comparison with reference datasets "bias", "distance_metric", + "ttest", # Remove supplementary variables from cube "remove_supplementary_variables", # Save to file @@ -256,6 +257,7 @@ MULTI_MODEL_FUNCTIONS: set[str] = { "bias", "distance_metric", + "ttest", "ensemble_statistics", "multi_model_statistics", "mask_multimodel", diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index b812ce754c..f7f396a3b1 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -45,6 +45,107 @@ BiasType = Literal["absolute", "relative"] +def ttest( + products: set[PreprocessorFile] | Iterable[Cube], + reference: Cube | None = None, + siglvl: float = 0.05, +) -> set[PreprocessorFile] | CubeList: + """Calculate the ttest between dataset and reference dataset. + + The reference dataset needs to be broadcastable to all input `products`. + This supports `iris' rich broadcasting abilities + `__. To ensure this, the preprocessors + :func:`esmvalcore.preprocessor.regrid` and/or + :func:`esmvalcore.preprocessor.regrid_time` might be helpful. + + Notes + ----- + The reference dataset can be specified with the `reference` argument. If + `reference` is ``None``, exactly one input dataset in the `products` set + needs to have the facet ``reference_for_bias: true`` defined in the recipe. + Please do **not** specify the option `reference` when using this + preprocessor function in a recipe. + + Parameters + ---------- + products: + Input datasets/cubes for which the significance is calculated relative + to a reference dataset/cube. + reference: + Cube which is used as reference for the bttest calculation. If ``None``, + `products` needs to be a :obj:`set` of + `~esmvalcore.preprocessor.PreprocessorFile` objects and exactly one + dataset in `products` needs the facet ``reference_for_bias: true``. Do + not specify this argument in a recipe. + siglvl: + significance level which is used to create the mask of significant and + non-significant data points. + + Returns + ------- + set[PreprocessorFile] | CubeList + Output datasets/cubes. Will be a :obj:`set` of + :class:`~esmvalcore.preprocessor.PreprocessorFile` objects if + `products` is also one, a :class:`~iris.cube.CubeList` otherwise. + The cubes contain an attached ancillary variable which is a mask of the + cube defining significant (1) and non-significant (0) values related + to the bias (in terms of absolute) of cube and reference. + + Raises + ------ + ValueError + Not exactly one input datasets contains the facet ``reference_for_bias: + true`` if ``reference=None``; ``reference=None`` and the input products + are given as iterable of :class:`~iris.cube.Cube` objects; + ``bias_type`` is not one of ``'absolute'`` or ``'relative'``. + """ + ref_product = None + all_cubes_given = all(isinstance(p, Cube) for p in products) + + # Get reference cube if not explicitly given + if reference is None: + if all_cubes_given: + msg = ( + "A list of Cubes is given to this preprocessor; please " + "specify a `reference`" + ) + raise ValueError(msg) + (reference, ref_product) = _get_ref(products, "reference_for_bias") + else: + ref_product = None + + # If input is an Iterable of Cube objects, calculate ttest for each element + if all_cubes_given: + cubes = [_calculate_ttest(c, reference, siglvl) for c in products] + return CubeList(cubes) + + # Otherwise, iterate over all input products, calculate ttest and adapt + # metadata and provenance information accordingly + output_products = set() + for product in products: + if product == ref_product: + continue + cube = concatenate(product.cubes) + + # Calculate ttest + cube = _calculate_ttest(cube, reference, siglvl) + + # Adapt metadata and provenance information + product.attributes["units"] = str(cube.units) + if ref_product is not None: + product.wasderivedfrom(ref_product) + + product.cubes = CubeList([cube]) + output_products.add(product) + print(cube) + # Add reference dataset to output (always) + if ref_product is not None: + output_products.add(ref_product) + + return output_products + + def bias( products: set[PreprocessorFile] | Iterable[Cube], reference: Cube | None = None, @@ -208,8 +309,8 @@ def _calculate_bias(cube: Cube, reference: Cube, bias_type: BiasType, ttest: boo """Calculate bias for a single cube relative to a reference cube.""" cube_metadata = cube.metadata - if ttest: - p_value, remaining_axis = _calculate_ttest(cube, reference) + # if ttest: + # p_value, remaining_axis = _calculate_ttest(cube, reference) if bias_type == "absolute": cube = cube - reference @@ -226,27 +327,51 @@ def _calculate_bias(cube: Cube, reference: Cube, bias_type: BiasType, ttest: boo cube.metadata = cube_metadata cube.units = new_units + # Attach ancillary variable after the calculation # otherwise, it will be lost - if ttest: - cube.add_ancillary_variable( - AncillaryVariable(p_value, - long_name="pvalue"), - remaining_axis) + # if ttest: + # cube.add_ancillary_variable( + # AncillaryVariable(p_value, + # long_name="pvalue"), + # remaining_axis) return cube -def _calculate_ttest(cube: Cube, reference: Cube) -> Cube: - """Calculate the Welch's t-test and attach p_value as an ancillary variable.""" +# def _calculate_ttest(cube: Cube, reference: Cube) -> Cube: + # """Calculate the Welch's t-test and attach p_value as an ancillary variable.""" + # axis = cube.coord_dims("time")[0] + # remaining_axis = sorted(set(range(cube.ndim)) - set([axis])) + + # ttest = dask.array.stats.ttest_ind(cube, reference, 0, False) + # shape = tuple([list(cube.shape)[dim] for dim in remaining_axis]) + # p_value = da.from_delayed(ttest.pvalue, shape, dtype=float) + + # return p_value, remaining_axis + + +def _calculate_ttest(cube: Cube, reference: Cube, siglvl: float) -> Cube: + """Calculate the Welch's t-test and attach a mask of significance.""" + cube_metadata = cube.metadata + axis = cube.coord_dims("time")[0] - remaining_axis = sorted(set(range(cube.ndim)) - set([axis])) + remaining_axes = sorted(set(range(cube.ndim)) - set([axis])) ttest = dask.array.stats.ttest_ind(cube, reference, 0, False) - shape = tuple([list(cube.shape)[dim] for dim in remaining_axis]) + shape = tuple([list(cube.shape)[dim] for dim in remaining_axes]) p_value = da.from_delayed(ttest.pvalue, shape, dtype=float) - return p_value, remaining_axis + significance = da.where(p_value <= siglvl, 1, 0) + + cube.add_ancillary_variable( + AncillaryVariable(significance, long_name="significance"), + remaining_axes + ) + + cube.metadata = cube_metadata + + return cube MetricType = Literal[ @@ -267,7 +392,7 @@ def distance_metric( keep_reference_dataset: bool = True, **kwargs: Any, ) -> set[PreprocessorFile] | CubeList: - r"""Calculate distance metrics. + """Calculate distance metrics. All input datasets need to have identical dimensional coordinates. This can for example be ensured with the preprocessors From fc338d067324766b1d2cb2dcb425aea873592564 Mon Sep 17 00:00:00 2001 From: Franziska Winterstein Date: Tue, 23 Jun 2026 18:06:50 +0200 Subject: [PATCH 04/33] Reverted code s.t. p_value is attached and not a mask 0/1 --- esmvalcore/preprocessor/_compare_with_refs.py | 86 +++++++++---------- 1 file changed, 40 insertions(+), 46 deletions(-) diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index f7f396a3b1..7b21a75887 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -48,7 +48,9 @@ def ttest( products: set[PreprocessorFile] | Iterable[Cube], reference: Cube | None = None, - siglvl: float = 0.05, + coordinate: str | None = 'time', + equal_var: bool = False, + **kwargs: Any, ) -> set[PreprocessorFile] | CubeList: """Calculate the ttest between dataset and reference dataset. @@ -73,24 +75,31 @@ def ttest( Input datasets/cubes for which the significance is calculated relative to a reference dataset/cube. reference: - Cube which is used as reference for the bttest calculation. If ``None``, + Cube which is used as reference for the ttest calculation. If ``None``, `products` needs to be a :obj:`set` of `~esmvalcore.preprocessor.PreprocessorFile` objects and exactly one dataset in `products` needs the facet ``reference_for_bias: true``. Do not specify this argument in a recipe. - siglvl: - significance level which is used to create the mask of significant and - non-significant data points. + coordinate: + Coordinate axis of the input along which to compute the statistic. + Default is 'time'. + equal_var: + If False (default), perform Welch's t-test, which does not assume equal + population variances. Can be set to True, to perform a standard + independent 2 sample test that assumes equal population variances, if + equality is assured + **kwargs: + Any other parameter allowed by dask.array.ttest_ind. Returns ------- set[PreprocessorFile] | CubeList Output datasets/cubes. Will be a :obj:`set` of - :class:`~esmvalcore.preprocessor.PreprocessorFile` objects if - `products` is also one, a :class:`~iris.cube.CubeList` otherwise. - The cubes contain an attached ancillary variable which is a mask of the - cube defining significant (1) and non-significant (0) values related - to the bias (in terms of absolute) of cube and reference. + :class:`~esmvalcore.preprocessor.PreprocessorFile` objects if `products` + is also one, a :class:`~iris.cube.CubeList` otherwise. The cubes + contain an attached ancillary variable which is the p-value w.r.t. the + chosen ttest related to the bias (in terms of absolute) of cube and + reference. Raises ------ @@ -99,6 +108,7 @@ def ttest( true`` if ``reference=None``; ``reference=None`` and the input products are given as iterable of :class:`~iris.cube.Cube` objects; ``bias_type`` is not one of ``'absolute'`` or ``'relative'``. + """ ref_product = None all_cubes_given = all(isinstance(p, Cube) for p in products) @@ -117,7 +127,8 @@ def ttest( # If input is an Iterable of Cube objects, calculate ttest for each element if all_cubes_given: - cubes = [_calculate_ttest(c, reference, siglvl) for c in products] + cubes = [_calculate_ttest(c, reference, coordinate, equal_var, **kwargs) + for c in products] return CubeList(cubes) # Otherwise, iterate over all input products, calculate ttest and adapt @@ -129,7 +140,8 @@ def ttest( cube = concatenate(product.cubes) # Calculate ttest - cube = _calculate_ttest(cube, reference, siglvl) + cube = _calculate_ttest(cube, reference, coordinate, equal_var, + **kwargs) # Adapt metadata and provenance information product.attributes["units"] = str(cube.units) @@ -138,7 +150,7 @@ def ttest( product.cubes = CubeList([cube]) output_products.add(product) - print(cube) + # Add reference dataset to output (always) if ref_product is not None: output_products.add(ref_product) @@ -185,11 +197,6 @@ def bias( bias_type: Bias type that is calculated. Must be one of ``'absolute'`` (dataset - ref) or ``'relative'`` ((dataset - ref) / ref). - ttest: - If ``True`` calculate the p-value according to the Welch's Test between - dataset and reference. The p-value is attached as an ancillary variable - to the dataset. A Welch's Test is chosen, since unequal variances and - sample sizes are likely. denominator_mask_threshold: Threshold to mask values close to zero in the denominator (i.e., the reference dataset) during the calculation of relative biases. All @@ -309,8 +316,9 @@ def _calculate_bias(cube: Cube, reference: Cube, bias_type: BiasType, ttest: boo """Calculate bias for a single cube relative to a reference cube.""" cube_metadata = cube.metadata - # if ttest: - # p_value, remaining_axis = _calculate_ttest(cube, reference) + # save ancillary variables as they get removed in the bias calculation + anc_var = cube.ancillary_variables() + anc_var_dims = [ cube.ancillary_variable_dims(anc) for anc in anc_var ] if bias_type == "absolute": cube = cube - reference @@ -328,44 +336,30 @@ def _calculate_bias(cube: Cube, reference: Cube, bias_type: BiasType, ttest: boo cube.metadata = cube_metadata cube.units = new_units - # Attach ancillary variable after the calculation - # otherwise, it will be lost - # if ttest: - # cube.add_ancillary_variable( - # AncillaryVariable(p_value, - # long_name="pvalue"), - # remaining_axis) + # reattach ancillary variables + for anc, dims in zip(anc_var, anc_var_dims): + cube.add_ancillary_variable(anc, dims) return cube -# def _calculate_ttest(cube: Cube, reference: Cube) -> Cube: - # """Calculate the Welch's t-test and attach p_value as an ancillary variable.""" - # axis = cube.coord_dims("time")[0] - # remaining_axis = sorted(set(range(cube.ndim)) - set([axis])) - - # ttest = dask.array.stats.ttest_ind(cube, reference, 0, False) - # shape = tuple([list(cube.shape)[dim] for dim in remaining_axis]) - # p_value = da.from_delayed(ttest.pvalue, shape, dtype=float) - - # return p_value, remaining_axis - - -def _calculate_ttest(cube: Cube, reference: Cube, siglvl: float) -> Cube: - """Calculate the Welch's t-test and attach a mask of significance.""" +def _calculate_ttest(cube: Cube, reference: Cube, coordinate, equal_var, **kwargs) -> Cube: + """Calculate the Welch's t-test and attach the p-value.""" cube_metadata = cube.metadata - axis = cube.coord_dims("time")[0] + axis = cube.coord_dims(coordinate)[0] remaining_axes = sorted(set(range(cube.ndim)) - set([axis])) - ttest = dask.array.stats.ttest_ind(cube, reference, 0, False) + ttest = dask.array.stats.ttest_ind(cube, + reference, + axis=axis, + equal_var=equal_var) shape = tuple([list(cube.shape)[dim] for dim in remaining_axes]) - p_value = da.from_delayed(ttest.pvalue, shape, dtype=float) - significance = da.where(p_value <= siglvl, 1, 0) + p_value = da.from_delayed(ttest.pvalue, shape, dtype=float) cube.add_ancillary_variable( - AncillaryVariable(significance, long_name="significance"), + AncillaryVariable(p_value, long_name="p_value"), remaining_axes ) From 6cf75acd7790b45a435607f54af82db4778de05c Mon Sep 17 00:00:00 2001 From: Franziska Winterstein Date: Wed, 24 Jun 2026 12:22:14 +0200 Subject: [PATCH 05/33] added example recipe where ttest is used --- test_ttest.yml | 90 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 90 insertions(+) create mode 100644 test_ttest.yml diff --git a/test_ttest.yml b/test_ttest.yml new file mode 100644 index 0000000000..18c20ed49b --- /dev/null +++ b/test_ttest.yml @@ -0,0 +1,90 @@ +# ESMValTool +# This recipe was successfully run with +# - ESMValTool version: +# - ESMValCore version: +--- +documentation: + title: Example recipe for monitor/multi_datasets.py + + description: This recipe includes the zonal mean run by the diagnostic monitor/multi_datasets.py and MESSy data. + + authors: + - winterstein_franziska + + maintainer: + - winterstein_franziska + + references: + - acknow_project + +datasets: + - {alias: ERA5, project: native6, dataset: ERA5, type: reanaly, version: v1, tier: 3, reference_for_bias: True} + - {alias: RD1-PSrad-02, dataset: EMAC, project: EMAC, ensemble: EMAC, exp: RD1-PSrad-02, supplementary_variables: [{short_name: areacella, skip: true}] } + +preprocessors: + + zonal_mean_bias: + custom_order: true + regrid: + scheme: linear + target_grid: ERA5 + lon_offset: false + extract_levels: + levels: {cmor_table: CMIP6, coordinate: plev39} + scheme: linear + coordinate: air_pressure + zonal_statistics: + operator: mean + annual_statistics: + operator: mean + regrid_time: + frequency: 'yr' + ttest: true + climate_statistics: + operator: mean + period: full + bias: + bias_type: absolute + # Note that if bias is used you need to specify 'reference_to_bias: True' + # in that dataset above, which should serve as the reference + # bias_type "absolute" is recommended for temperature + # for other tracer it could make sense to use "relative" + +diagnostics: + + ttest: + description: Zonal mean bias of temperature with ttest + themes: + - phys + realms: + - atmos + variables: + ta: + mip: Amon + channel: Amon + raw_name: tm1_p39_cav + start_year: 2000 + end_year: 2001 + preprocessor: zonal_mean_bias + scripts: + plot_zonal_mean: + script: monitor/multi_datasets.py + facets_used_for_labels: alias + plot_folder: '{plot_dir}' + plot_filename: '{plot_type}_{alias}-reference_{preprocessor}' + plots: + zonal_mean_profile: + siglvl: 0.05 + common_cbar: true + plot_kwargs: + default: + cmap: 'RdBu' + levels: [-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5] + pyplot_kwargs: + title: '{alias} - ERA5' + # this corresponds to the requirements of copernicus + xticks: + ticks: [-60, -30, 0, 30, 60] + labels: ["60° S", "30° S", "0°", "30° N", "60° N"] + xlabel: latitude + ylabel: Pressure [hPa] From 48844bcc2bdb3ec8b78c9d39a3c052384bcaac99 Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 09:07:26 +0200 Subject: [PATCH 06/33] Remove recipe from repo --- test_ttest.yml | 90 -------------------------------------------------- 1 file changed, 90 deletions(-) delete mode 100644 test_ttest.yml diff --git a/test_ttest.yml b/test_ttest.yml deleted file mode 100644 index 18c20ed49b..0000000000 --- a/test_ttest.yml +++ /dev/null @@ -1,90 +0,0 @@ -# ESMValTool -# This recipe was successfully run with -# - ESMValTool version: -# - ESMValCore version: ---- -documentation: - title: Example recipe for monitor/multi_datasets.py - - description: This recipe includes the zonal mean run by the diagnostic monitor/multi_datasets.py and MESSy data. - - authors: - - winterstein_franziska - - maintainer: - - winterstein_franziska - - references: - - acknow_project - -datasets: - - {alias: ERA5, project: native6, dataset: ERA5, type: reanaly, version: v1, tier: 3, reference_for_bias: True} - - {alias: RD1-PSrad-02, dataset: EMAC, project: EMAC, ensemble: EMAC, exp: RD1-PSrad-02, supplementary_variables: [{short_name: areacella, skip: true}] } - -preprocessors: - - zonal_mean_bias: - custom_order: true - regrid: - scheme: linear - target_grid: ERA5 - lon_offset: false - extract_levels: - levels: {cmor_table: CMIP6, coordinate: plev39} - scheme: linear - coordinate: air_pressure - zonal_statistics: - operator: mean - annual_statistics: - operator: mean - regrid_time: - frequency: 'yr' - ttest: true - climate_statistics: - operator: mean - period: full - bias: - bias_type: absolute - # Note that if bias is used you need to specify 'reference_to_bias: True' - # in that dataset above, which should serve as the reference - # bias_type "absolute" is recommended for temperature - # for other tracer it could make sense to use "relative" - -diagnostics: - - ttest: - description: Zonal mean bias of temperature with ttest - themes: - - phys - realms: - - atmos - variables: - ta: - mip: Amon - channel: Amon - raw_name: tm1_p39_cav - start_year: 2000 - end_year: 2001 - preprocessor: zonal_mean_bias - scripts: - plot_zonal_mean: - script: monitor/multi_datasets.py - facets_used_for_labels: alias - plot_folder: '{plot_dir}' - plot_filename: '{plot_type}_{alias}-reference_{preprocessor}' - plots: - zonal_mean_profile: - siglvl: 0.05 - common_cbar: true - plot_kwargs: - default: - cmap: 'RdBu' - levels: [-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5] - pyplot_kwargs: - title: '{alias} - ERA5' - # this corresponds to the requirements of copernicus - xticks: - ticks: [-60, -30, 0, 30, 60] - labels: ["60° S", "30° S", "0°", "30° N", "60° N"] - xlabel: latitude - ylabel: Pressure [hPa] From 04802a9e58f1dbfb5983bc06ae5f7b06c00a2b4e Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 09:24:10 +0200 Subject: [PATCH 07/33] pre-commit --- esmvalcore/preprocessor/__init__.py | 10 ++- esmvalcore/preprocessor/_compare_with_refs.py | 82 ++++++++++--------- 2 files changed, 51 insertions(+), 41 deletions(-) diff --git a/esmvalcore/preprocessor/__init__.py b/esmvalcore/preprocessor/__init__.py index 4a16f06207..6418a0289f 100644 --- a/esmvalcore/preprocessor/__init__.py +++ b/esmvalcore/preprocessor/__init__.py @@ -31,7 +31,11 @@ meridional_statistics, zonal_statistics, ) -from esmvalcore.preprocessor._compare_with_refs import bias, distance_metric, ttest +from esmvalcore.preprocessor._compare_with_refs import ( + bias, + distance_metric, + t_test, +) from esmvalcore.preprocessor._concatenate import concatenate from esmvalcore.preprocessor._cycles import amplitude from esmvalcore.preprocessor._dask_progress import _compute_with_progress @@ -219,7 +223,7 @@ # Comparison with reference datasets "bias", "distance_metric", - "ttest", + "t_test", # Remove supplementary variables from cube "remove_supplementary_variables", # Save to file @@ -257,7 +261,7 @@ MULTI_MODEL_FUNCTIONS: set[str] = { "bias", "distance_metric", - "ttest", + "t_test", "ensemble_statistics", "multi_model_statistics", "mask_multimodel", diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index 7b21a75887..f80b6eb8d4 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -13,8 +13,7 @@ import iris.analysis.stats import numpy as np from iris.common.metadata import CubeMetadata -from iris.coords import AncillaryVariable -from iris.coords import CellMethod +from iris.coords import AncillaryVariable, AuxCoord, CellMethod, DimCoord from iris.cube import Cube, CubeList from scipy.stats import wasserstein_distance @@ -45,14 +44,13 @@ BiasType = Literal["absolute", "relative"] -def ttest( +def t_test( products: set[PreprocessorFile] | Iterable[Cube], reference: Cube | None = None, - coordinate: str | None = 'time', - equal_var: bool = False, + coordinate: str | None = "time", **kwargs: Any, ) -> set[PreprocessorFile] | CubeList: - """Calculate the ttest between dataset and reference dataset. + """Calculate the t-test between dataset and reference dataset. The reference dataset needs to be broadcastable to all input `products`. This supports `iris' rich broadcasting abilities @@ -75,7 +73,7 @@ def ttest( Input datasets/cubes for which the significance is calculated relative to a reference dataset/cube. reference: - Cube which is used as reference for the ttest calculation. If ``None``, + Cube which is used as reference for the t-test calculation. If ``None``, `products` needs to be a :obj:`set` of `~esmvalcore.preprocessor.PreprocessorFile` objects and exactly one dataset in `products` needs the facet ``reference_for_bias: true``. Do @@ -83,13 +81,8 @@ def ttest( coordinate: Coordinate axis of the input along which to compute the statistic. Default is 'time'. - equal_var: - If False (default), perform Welch's t-test, which does not assume equal - population variances. Can be set to True, to perform a standard - independent 2 sample test that assumes equal population variances, if - equality is assured **kwargs: - Any other parameter allowed by dask.array.ttest_ind. + Additional keyword arguments passed to :func:`dask.array.ttest_ind`. Returns ------- @@ -98,7 +91,7 @@ def ttest( :class:`~esmvalcore.preprocessor.PreprocessorFile` objects if `products` is also one, a :class:`~iris.cube.CubeList` otherwise. The cubes contain an attached ancillary variable which is the p-value w.r.t. the - chosen ttest related to the bias (in terms of absolute) of cube and + chosen t-test related to the bias (in terms of absolute) of cube and reference. Raises @@ -125,13 +118,15 @@ def ttest( else: ref_product = None - # If input is an Iterable of Cube objects, calculate ttest for each element + # If input is an Iterable of Cube objects, calculate t-test for each element if all_cubes_given: - cubes = [_calculate_ttest(c, reference, coordinate, equal_var, **kwargs) - for c in products] + cubes = [ + _calculate_t_test(c, reference, coordinate, **kwargs) + for c in products + ] return CubeList(cubes) - # Otherwise, iterate over all input products, calculate ttest and adapt + # Otherwise, iterate over all input products, calculate t-test and adapt # metadata and provenance information accordingly output_products = set() for product in products: @@ -139,9 +134,8 @@ def ttest( continue cube = concatenate(product.cubes) - # Calculate ttest - cube = _calculate_ttest(cube, reference, coordinate, equal_var, - **kwargs) + # Calculate t-test + cube = _calculate_t_test(cube, reference, coordinate, **kwargs) # Adapt metadata and provenance information product.attributes["units"] = str(cube.units) @@ -162,7 +156,6 @@ def bias( products: set[PreprocessorFile] | Iterable[Cube], reference: Cube | None = None, bias_type: BiasType = "absolute", - ttest: bool = False, denominator_mask_threshold: float = 1e-3, keep_reference_dataset: bool = False, ) -> set[PreprocessorFile] | CubeList: @@ -254,7 +247,7 @@ def bias( # If input is an Iterable of Cube objects, calculate bias for each element if all_cubes_given: - cubes = [_calculate_bias(c, reference, bias_type, ttest) for c in products] + cubes = [_calculate_bias(c, reference, bias_type) for c in products] return CubeList(cubes) # Otherwise, iterate over all input products, calculate bias and adapt @@ -266,7 +259,7 @@ def bias( cube = concatenate(product.cubes) # Calculate bias - cube = _calculate_bias(cube, reference, bias_type, ttest) + cube = _calculate_bias(cube, reference, bias_type) # Adapt metadata and provenance information product.attributes["units"] = str(cube.units) @@ -312,13 +305,19 @@ def _get_ref( return (reference, ref_product) -def _calculate_bias(cube: Cube, reference: Cube, bias_type: BiasType, ttest: bool) -> Cube: +def _calculate_bias( + cube: Cube, + reference: Cube, + bias_type: BiasType, +) -> Cube: """Calculate bias for a single cube relative to a reference cube.""" cube_metadata = cube.metadata # save ancillary variables as they get removed in the bias calculation - anc_var = cube.ancillary_variables() - anc_var_dims = [ cube.ancillary_variable_dims(anc) for anc in anc_var ] + ancillary_variables_and_dims = [ + (a, cube.ancillary_variable_dims(a)) + for a in cube.ancillary_variables() + ] if bias_type == "absolute": cube = cube - reference @@ -337,30 +336,37 @@ def _calculate_bias(cube: Cube, reference: Cube, bias_type: BiasType, ttest: boo cube.units = new_units # reattach ancillary variables - for anc, dims in zip(anc_var, anc_var_dims): - cube.add_ancillary_variable(anc, dims) + for ancillary_variable, dims in ancillary_variables_and_dims: + cube.add_ancillary_variable(ancillary_variable, dims) return cube -def _calculate_ttest(cube: Cube, reference: Cube, coordinate, equal_var, **kwargs) -> Cube: +def _calculate_t_test( + cube: Cube, + reference: Cube, + coordinate: str | AuxCoord | DimCoord, + **kwargs: Any, +) -> Cube: """Calculate the Welch's t-test and attach the p-value.""" cube_metadata = cube.metadata axis = cube.coord_dims(coordinate)[0] - remaining_axes = sorted(set(range(cube.ndim)) - set([axis])) + remaining_axes = sorted(set(range(cube.ndim)) - {axis}) - ttest = dask.array.stats.ttest_ind(cube, - reference, - axis=axis, - equal_var=equal_var) + t_test = dask.array.stats.ttest_ind( + cube, + reference, + axis=axis, + **kwargs, + ) shape = tuple([list(cube.shape)[dim] for dim in remaining_axes]) - p_value = da.from_delayed(ttest.pvalue, shape, dtype=float) + p_value = da.from_delayed(t_test.pvalue, shape, dtype=float) cube.add_ancillary_variable( AncillaryVariable(p_value, long_name="p_value"), - remaining_axes + remaining_axes, ) cube.metadata = cube_metadata @@ -386,7 +392,7 @@ def distance_metric( keep_reference_dataset: bool = True, **kwargs: Any, ) -> set[PreprocessorFile] | CubeList: - """Calculate distance metrics. + r"""Calculate distance metrics. All input datasets need to have identical dimensional coordinates. This can for example be ensured with the preprocessors From 1caccba288fedc4d06948b056134560cf8ddb6ed Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 09:24:38 +0200 Subject: [PATCH 08/33] Optimize preprocessor order --- esmvalcore/preprocessor/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/esmvalcore/preprocessor/__init__.py b/esmvalcore/preprocessor/__init__.py index 6418a0289f..520c67ddb0 100644 --- a/esmvalcore/preprocessor/__init__.py +++ b/esmvalcore/preprocessor/__init__.py @@ -221,9 +221,9 @@ # Multi model statistics "multi_model_statistics", # Comparison with reference datasets + "t_test", "bias", "distance_metric", - "t_test", # Remove supplementary variables from cube "remove_supplementary_variables", # Save to file @@ -261,11 +261,11 @@ MULTI_MODEL_FUNCTIONS: set[str] = { "bias", "distance_metric", - "t_test", "ensemble_statistics", "multi_model_statistics", "mask_multimodel", "mask_fillvalues", + "t_test", } From b4e86cf19f643accc3cde9f340e8293cf0ff5723 Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 09:25:57 +0200 Subject: [PATCH 09/33] Optimize function order --- esmvalcore/preprocessor/_compare_with_refs.py | 280 +++++++++--------- 1 file changed, 140 insertions(+), 140 deletions(-) diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index f80b6eb8d4..afc7769414 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -44,114 +44,6 @@ BiasType = Literal["absolute", "relative"] -def t_test( - products: set[PreprocessorFile] | Iterable[Cube], - reference: Cube | None = None, - coordinate: str | None = "time", - **kwargs: Any, -) -> set[PreprocessorFile] | CubeList: - """Calculate the t-test between dataset and reference dataset. - - The reference dataset needs to be broadcastable to all input `products`. - This supports `iris' rich broadcasting abilities - `__. To ensure this, the preprocessors - :func:`esmvalcore.preprocessor.regrid` and/or - :func:`esmvalcore.preprocessor.regrid_time` might be helpful. - - Notes - ----- - The reference dataset can be specified with the `reference` argument. If - `reference` is ``None``, exactly one input dataset in the `products` set - needs to have the facet ``reference_for_bias: true`` defined in the recipe. - Please do **not** specify the option `reference` when using this - preprocessor function in a recipe. - - Parameters - ---------- - products: - Input datasets/cubes for which the significance is calculated relative - to a reference dataset/cube. - reference: - Cube which is used as reference for the t-test calculation. If ``None``, - `products` needs to be a :obj:`set` of - `~esmvalcore.preprocessor.PreprocessorFile` objects and exactly one - dataset in `products` needs the facet ``reference_for_bias: true``. Do - not specify this argument in a recipe. - coordinate: - Coordinate axis of the input along which to compute the statistic. - Default is 'time'. - **kwargs: - Additional keyword arguments passed to :func:`dask.array.ttest_ind`. - - Returns - ------- - set[PreprocessorFile] | CubeList - Output datasets/cubes. Will be a :obj:`set` of - :class:`~esmvalcore.preprocessor.PreprocessorFile` objects if `products` - is also one, a :class:`~iris.cube.CubeList` otherwise. The cubes - contain an attached ancillary variable which is the p-value w.r.t. the - chosen t-test related to the bias (in terms of absolute) of cube and - reference. - - Raises - ------ - ValueError - Not exactly one input datasets contains the facet ``reference_for_bias: - true`` if ``reference=None``; ``reference=None`` and the input products - are given as iterable of :class:`~iris.cube.Cube` objects; - ``bias_type`` is not one of ``'absolute'`` or ``'relative'``. - - """ - ref_product = None - all_cubes_given = all(isinstance(p, Cube) for p in products) - - # Get reference cube if not explicitly given - if reference is None: - if all_cubes_given: - msg = ( - "A list of Cubes is given to this preprocessor; please " - "specify a `reference`" - ) - raise ValueError(msg) - (reference, ref_product) = _get_ref(products, "reference_for_bias") - else: - ref_product = None - - # If input is an Iterable of Cube objects, calculate t-test for each element - if all_cubes_given: - cubes = [ - _calculate_t_test(c, reference, coordinate, **kwargs) - for c in products - ] - return CubeList(cubes) - - # Otherwise, iterate over all input products, calculate t-test and adapt - # metadata and provenance information accordingly - output_products = set() - for product in products: - if product == ref_product: - continue - cube = concatenate(product.cubes) - - # Calculate t-test - cube = _calculate_t_test(cube, reference, coordinate, **kwargs) - - # Adapt metadata and provenance information - product.attributes["units"] = str(cube.units) - if ref_product is not None: - product.wasderivedfrom(ref_product) - - product.cubes = CubeList([cube]) - output_products.add(product) - - # Add reference dataset to output (always) - if ref_product is not None: - output_products.add(ref_product) - - return output_products - - def bias( products: set[PreprocessorFile] | Iterable[Cube], reference: Cube | None = None, @@ -342,38 +234,6 @@ def _calculate_bias( return cube -def _calculate_t_test( - cube: Cube, - reference: Cube, - coordinate: str | AuxCoord | DimCoord, - **kwargs: Any, -) -> Cube: - """Calculate the Welch's t-test and attach the p-value.""" - cube_metadata = cube.metadata - - axis = cube.coord_dims(coordinate)[0] - remaining_axes = sorted(set(range(cube.ndim)) - {axis}) - - t_test = dask.array.stats.ttest_ind( - cube, - reference, - axis=axis, - **kwargs, - ) - shape = tuple([list(cube.shape)[dim] for dim in remaining_axes]) - - p_value = da.from_delayed(t_test.pvalue, shape, dtype=float) - - cube.add_ancillary_variable( - AncillaryVariable(p_value, long_name="p_value"), - remaining_axes, - ) - - cube.metadata = cube_metadata - - return cube - - MetricType = Literal[ "rmse", "weighted_rmse", @@ -752,3 +612,143 @@ def _get_emd( if np.ma.is_masked(arr) or np.ma.is_masked(ref_arr): return np.ma.masked # this is safe because PMFs will be masked arrays return wasserstein_distance(bin_centers, bin_centers, arr, ref_arr) + + +def t_test( + products: set[PreprocessorFile] | Iterable[Cube], + reference: Cube | None = None, + coordinate: str | None = "time", + **kwargs: Any, +) -> set[PreprocessorFile] | CubeList: + """Calculate the t-test between dataset and reference dataset. + + The reference dataset needs to be broadcastable to all input `products`. + This supports `iris' rich broadcasting abilities + `__. To ensure this, the preprocessors + :func:`esmvalcore.preprocessor.regrid` and/or + :func:`esmvalcore.preprocessor.regrid_time` might be helpful. + + Notes + ----- + The reference dataset can be specified with the `reference` argument. If + `reference` is ``None``, exactly one input dataset in the `products` set + needs to have the facet ``reference_for_bias: true`` defined in the recipe. + Please do **not** specify the option `reference` when using this + preprocessor function in a recipe. + + Parameters + ---------- + products: + Input datasets/cubes for which the significance is calculated relative + to a reference dataset/cube. + reference: + Cube which is used as reference for the t-test calculation. If ``None``, + `products` needs to be a :obj:`set` of + `~esmvalcore.preprocessor.PreprocessorFile` objects and exactly one + dataset in `products` needs the facet ``reference_for_bias: true``. Do + not specify this argument in a recipe. + coordinate: + Coordinate axis of the input along which to compute the statistic. + Default is 'time'. + **kwargs: + Additional keyword arguments passed to :func:`dask.array.ttest_ind`. + + Returns + ------- + set[PreprocessorFile] | CubeList + Output datasets/cubes. Will be a :obj:`set` of + :class:`~esmvalcore.preprocessor.PreprocessorFile` objects if `products` + is also one, a :class:`~iris.cube.CubeList` otherwise. The cubes + contain an attached ancillary variable which is the p-value w.r.t. the + chosen t-test related to the bias (in terms of absolute) of cube and + reference. + + Raises + ------ + ValueError + Not exactly one input datasets contains the facet ``reference_for_bias: + true`` if ``reference=None``; ``reference=None`` and the input products + are given as iterable of :class:`~iris.cube.Cube` objects; + ``bias_type`` is not one of ``'absolute'`` or ``'relative'``. + + """ + ref_product = None + all_cubes_given = all(isinstance(p, Cube) for p in products) + + # Get reference cube if not explicitly given + if reference is None: + if all_cubes_given: + msg = ( + "A list of Cubes is given to this preprocessor; please " + "specify a `reference`" + ) + raise ValueError(msg) + (reference, ref_product) = _get_ref(products, "reference_for_bias") + else: + ref_product = None + + # If input is an Iterable of Cube objects, calculate t-test for each element + if all_cubes_given: + cubes = [ + _calculate_t_test(c, reference, coordinate, **kwargs) + for c in products + ] + return CubeList(cubes) + + # Otherwise, iterate over all input products, calculate t-test and adapt + # metadata and provenance information accordingly + output_products = set() + for product in products: + if product == ref_product: + continue + cube = concatenate(product.cubes) + + # Calculate t-test + cube = _calculate_t_test(cube, reference, coordinate, **kwargs) + + # Adapt metadata and provenance information + product.attributes["units"] = str(cube.units) + if ref_product is not None: + product.wasderivedfrom(ref_product) + + product.cubes = CubeList([cube]) + output_products.add(product) + + # Add reference dataset to output (always) + if ref_product is not None: + output_products.add(ref_product) + + return output_products + + +def _calculate_t_test( + cube: Cube, + reference: Cube, + coordinate: str | AuxCoord | DimCoord, + **kwargs: Any, +) -> Cube: + """Calculate the Welch's t-test and attach the p-value.""" + cube_metadata = cube.metadata + + axis = cube.coord_dims(coordinate)[0] + remaining_axes = sorted(set(range(cube.ndim)) - {axis}) + + t_test = dask.array.stats.ttest_ind( + cube, + reference, + axis=axis, + **kwargs, + ) + shape = tuple([list(cube.shape)[dim] for dim in remaining_axes]) + + p_value = da.from_delayed(t_test.pvalue, shape, dtype=float) + + cube.add_ancillary_variable( + AncillaryVariable(p_value, long_name="p_value"), + remaining_axes, + ) + + cube.metadata = cube_metadata + + return cube From 720faab93424ea5f0b28bb26b6e89d8c7b3ea8c5 Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 10:14:20 +0200 Subject: [PATCH 10/33] Add support for lazy and non-lazy data --- esmvalcore/preprocessor/_compare_with_refs.py | 84 +++++++++++-------- 1 file changed, 49 insertions(+), 35 deletions(-) diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index afc7769414..6be0a48d6b 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -8,14 +8,13 @@ import dask import dask.array as da -import dask.array.stats import iris.analysis import iris.analysis.stats import numpy as np from iris.common.metadata import CubeMetadata -from iris.coords import AncillaryVariable, AuxCoord, CellMethod, DimCoord +from iris.coords import AncillaryVariable, CellMethod from iris.cube import Cube, CubeList -from scipy.stats import wasserstein_distance +from scipy.stats import wasserstein_distance, ttest_ind from esmvalcore.iris_helpers import ( ignore_iris_vague_metadata_warnings, @@ -34,7 +33,7 @@ if TYPE_CHECKING: from collections.abc import Iterable - from iris.coords import Coord + from iris.coords import AuxCoord, Coord, DimCoord from esmvalcore.preprocessor import PreprocessorFile @@ -168,7 +167,7 @@ def bias( return output_products -def _get_ref( +def _get_ref( products: Iterable[PreprocessorFile], ref_tag: str, ) -> tuple[Cube, PreprocessorFile]: @@ -617,10 +616,13 @@ def _get_emd( def t_test( products: set[PreprocessorFile] | Iterable[Cube], reference: Cube | None = None, - coordinate: str | None = "time", + coordinate: str | AuxCoord | DimCoord | None = "time", **kwargs: Any, ) -> set[PreprocessorFile] | CubeList: - """Calculate the t-test between dataset and reference dataset. + """Perform t-test between dataset and reference dataset. + + This is a test for the null hypothesis that 2 independent samples (dataset + and reference dataset) have identical average (expected) values. The reference dataset needs to be broadcastable to all input `products`. This supports `iris' rich broadcasting abilities @@ -633,8 +635,8 @@ def t_test( ----- The reference dataset can be specified with the `reference` argument. If `reference` is ``None``, exactly one input dataset in the `products` set - needs to have the facet ``reference_for_bias: true`` defined in the recipe. - Please do **not** specify the option `reference` when using this + needs to have the facet ``reference_for_t_test: true`` defined in the + recipe. Please do **not** specify the option `reference` when using this preprocessor function in a recipe. Parameters @@ -643,34 +645,35 @@ def t_test( Input datasets/cubes for which the significance is calculated relative to a reference dataset/cube. reference: - Cube which is used as reference for the t-test calculation. If ``None``, - `products` needs to be a :obj:`set` of + Cube which is used as reference for the t-test calculation. If + ``None``, `products` needs to be a :obj:`set` of `~esmvalcore.preprocessor.PreprocessorFile` objects and exactly one - dataset in `products` needs the facet ``reference_for_bias: true``. Do - not specify this argument in a recipe. + dataset in `products` needs the facet ``reference_for_t_test: true``. + Do not specify this argument in a recipe. coordinate: - Coordinate axis of the input along which to compute the statistic. - Default is 'time'. + Coordinate axis of the input along which to compute the statistic. If + ``None``, the input will be raveled before computing the statistic. **kwargs: - Additional keyword arguments passed to :func:`dask.array.ttest_ind`. + Additional keyword arguments passed to :func:`scipy.stats.ttest_ind` + (not all input data are lazy) or + :func:`dask.array.stats.ttest_ind.ttest_ind` (all input data are lazy). Returns ------- - set[PreprocessorFile] | CubeList + : Output datasets/cubes. Will be a :obj:`set` of - :class:`~esmvalcore.preprocessor.PreprocessorFile` objects if `products` - is also one, a :class:`~iris.cube.CubeList` otherwise. The cubes - contain an attached ancillary variable which is the p-value w.r.t. the - chosen t-test related to the bias (in terms of absolute) of cube and - reference. + :class:`~esmvalcore.preprocessor.PreprocessorFile` objects if + `products` is also one, a :class:`~iris.cube.CubeList` otherwise. The + cubes contain an attached ancillary variable which is the p-value + w.r.t. the chosen t-test. Raises ------ ValueError - Not exactly one input datasets contains the facet ``reference_for_bias: - true`` if ``reference=None``; ``reference=None`` and the input products - are given as iterable of :class:`~iris.cube.Cube` objects; - ``bias_type`` is not one of ``'absolute'`` or ``'relative'``. + Not exactly one input datasets contains the facet + ``reference_for_t_test: true`` if ``reference=None``; + ``reference=None`` and the input products are given as iterable of + :class:`~iris.cube.Cube` objectss. """ ref_product = None @@ -684,7 +687,7 @@ def t_test( "specify a `reference`" ) raise ValueError(msg) - (reference, ref_product) = _get_ref(products, "reference_for_bias") + (reference, ref_product) = _get_ref(products, "reference_for_t_test") else: ref_product = None @@ -708,7 +711,6 @@ def t_test( cube = _calculate_t_test(cube, reference, coordinate, **kwargs) # Adapt metadata and provenance information - product.attributes["units"] = str(cube.units) if ref_product is not None: product.wasderivedfrom(ref_product) @@ -725,21 +727,33 @@ def t_test( def _calculate_t_test( cube: Cube, reference: Cube, - coordinate: str | AuxCoord | DimCoord, + coordinate: str | AuxCoord | DimCoord | None, **kwargs: Any, ) -> Cube: - """Calculate the Welch's t-test and attach the p-value.""" + """Calculate the t-test and attach the p-value as ancillary variable to cube.""" cube_metadata = cube.metadata - axis = cube.coord_dims(coordinate)[0] - remaining_axes = sorted(set(range(cube.ndim)) - {axis}) + if coordinate is None: + axis = None + remaining_axes = list(range(cube.ndim)) + else: + axis = cube.coord_dims(coordinate)[0] + remaining_axes = sorted(set(range(cube.ndim)) - {axis}) - t_test = dask.array.stats.ttest_ind( - cube, - reference, + if cube.has_lazy_data() and reference.has_lazy_data(): + ttest_func = da.stats.ttest_ind + else: + ttest_func = ttest_ind + + t_test = ttest_func( + cube.core_data(), + reference.core_data(), axis=axis, **kwargs, ) + + print(t_test.pvalue) + assert 0 shape = tuple([list(cube.shape)[dim] for dim in remaining_axes]) p_value = da.from_delayed(t_test.pvalue, shape, dtype=float) From 81928338a19b2fc6d24856ab497a67cd3af076c2 Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 11:07:54 +0200 Subject: [PATCH 11/33] pre-commit --- esmvalcore/preprocessor/_compare_with_refs.py | 31 +++++++++---------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index 6be0a48d6b..2035088734 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -14,7 +14,7 @@ from iris.common.metadata import CubeMetadata from iris.coords import AncillaryVariable, CellMethod from iris.cube import Cube, CubeList -from scipy.stats import wasserstein_distance, ttest_ind +from scipy.stats import ttest_ind, wasserstein_distance from esmvalcore.iris_helpers import ( ignore_iris_vague_metadata_warnings, @@ -167,7 +167,7 @@ def bias( return output_products -def _get_ref( +def _get_ref( products: Iterable[PreprocessorFile], ref_tag: str, ) -> tuple[Cube, PreprocessorFile]: @@ -226,7 +226,7 @@ def _calculate_bias( cube.metadata = cube_metadata cube.units = new_units - # reattach ancillary variables + # Reattach ancillary variables for ancillary_variable, dims in ancillary_variables_and_dims: cube.add_ancillary_variable(ancillary_variable, dims) @@ -673,7 +673,7 @@ def t_test( Not exactly one input datasets contains the facet ``reference_for_t_test: true`` if ``reference=None``; ``reference=None`` and the input products are given as iterable of - :class:`~iris.cube.Cube` objectss. + :class:`~iris.cube.Cube` objects. """ ref_product = None @@ -717,7 +717,7 @@ def t_test( product.cubes = CubeList([cube]) output_products.add(product) - # Add reference dataset to output (always) + # Add reference dataset to output if ref_product is not None: output_products.add(ref_product) @@ -731,20 +731,21 @@ def _calculate_t_test( **kwargs: Any, ) -> Cube: """Calculate the t-test and attach the p-value as ancillary variable to cube.""" + cube = cube.copy() # do not modify input cube cube_metadata = cube.metadata + remaining_axes: tuple[int, ...] if coordinate is None: axis = None - remaining_axes = list(range(cube.ndim)) + remaining_axes = () else: axis = cube.coord_dims(coordinate)[0] - remaining_axes = sorted(set(range(cube.ndim)) - {axis}) + remaining_axes = tuple(sorted(set(range(cube.ndim)) - {axis})) if cube.has_lazy_data() and reference.has_lazy_data(): ttest_func = da.stats.ttest_ind else: ttest_func = ttest_ind - t_test = ttest_func( cube.core_data(), reference.core_data(), @@ -752,17 +753,15 @@ def _calculate_t_test( **kwargs, ) - print(t_test.pvalue) - assert 0 - shape = tuple([list(cube.shape)[dim] for dim in remaining_axes]) - - p_value = da.from_delayed(t_test.pvalue, shape, dtype=float) - cube.add_ancillary_variable( - AncillaryVariable(p_value, long_name="p_value"), + AncillaryVariable( + t_test.pvalue, + long_name="p-value", + var_name="pvalue", + units="1", + ), remaining_axes, ) - cube.metadata = cube_metadata return cube From 16f3c2a0081a2ca9e191360be0609bcc1aacd79a Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 11:08:05 +0200 Subject: [PATCH 12/33] Add tests for t_test with products --- .../test_compare_with_refs.py | 174 ++++++++++++++++-- 1 file changed, 158 insertions(+), 16 deletions(-) diff --git a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py index d9de9c7d4b..5271305e63 100644 --- a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py +++ b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py @@ -1,17 +1,22 @@ """Unit tests for :mod:`esmvalcore.preprocessor._compare_with_refs`.""" import contextlib +import re import dask.array as da import iris import numpy as np import pytest from cf_units import Unit -from iris.coords import CellMeasure, CellMethod +from iris.coords import AncillaryVariable, CellMeasure, CellMethod from iris.cube import Cube, CubeList from iris.exceptions import CoordinateNotFoundError -from esmvalcore.preprocessor._compare_with_refs import bias, distance_metric +from esmvalcore.preprocessor._compare_with_refs import ( + bias, + distance_metric, + t_test, +) from tests import PreprocessorFile @@ -100,20 +105,36 @@ def ref_cubes(): ] +@pytest.mark.parametrize("use_reference_product", [True, False]) @pytest.mark.parametrize(("bias_type", "data", "units"), TEST_BIAS) -def test_bias_products(regular_cubes, ref_cubes, bias_type, data, units): +def test_bias_products( + regular_cubes, + ref_cubes, + bias_type, + data, + units, + use_reference_product, +): """Test calculation of bias with products.""" - ref_product = PreprocessorFile( - ref_cubes, - "REF", - {"reference_for_bias": True}, - ) products = { PreprocessorFile(regular_cubes, "A", {"dataset": "a"}), PreprocessorFile(regular_cubes, "B", {"dataset": "b"}), - ref_product, } - out_products = bias(products, bias_type=bias_type) + + if use_reference_product: + ref_product = PreprocessorFile( + ref_cubes, + "REF", + {"reference_for_bias": True}, + ) + products.add(ref_product) + expected_ancestors = {ref_product} + reference = None + else: + expected_ancestors = set() + reference = ref_cubes[0] + + out_products = bias(products, reference=reference, bias_type=bias_type) assert isinstance(out_products, set) out_dict = products_set_to_dict(out_products) @@ -131,8 +152,8 @@ def test_bias_products(regular_cubes, ref_cubes, bias_type, data, units): assert out_cube.units == units assert out_cube.dim_coords == regular_cubes[0].dim_coords assert out_cube.aux_coords == regular_cubes[0].aux_coords - product_a.wasderivedfrom.assert_called_once() - assert product_a.mock_ancestors == {ref_product} + assert product_a.wasderivedfrom.call_count == len(expected_ancestors) + assert product_a.mock_ancestors == expected_ancestors product_b = out_dict["B"] assert product_b.filename == "B" @@ -146,8 +167,8 @@ def test_bias_products(regular_cubes, ref_cubes, bias_type, data, units): assert out_cube.units == units assert out_cube.dim_coords == regular_cubes[0].dim_coords assert out_cube.aux_coords == regular_cubes[0].aux_coords - product_b.wasderivedfrom.assert_called_once() - assert product_b.mock_ancestors == {ref_product} + assert product_b.wasderivedfrom.call_count == len(expected_ancestors) + assert product_b.mock_ancestors == expected_ancestors @pytest.mark.parametrize(("bias_type", "data", "units"), TEST_BIAS) @@ -382,8 +403,8 @@ def test_invalid_bias_type(regular_cubes, ref_cubes): bias(products, bias_type="invalid_bias_type") -def test_reference_none_cubes(regular_cubes): - """Test reference=None with with cubes.""" +def test_bias_reference_none_cubes(regular_cubes): + """Test bias with reference=None and cubes given.""" msg = ( "A list of Cubes is given to this preprocessor; please specify a " "`reference`" @@ -932,3 +953,124 @@ def test_distance_metric_no_lon_for_area_weights(regular_cubes, metric, error): reference=ref_cube, coords=["time", "latitude"], ) + + +@pytest.mark.parametrize("use_reference_product", [True, False]) +def test_t_test_products(regular_cubes, ref_cubes, use_reference_product): # noqa: PLR0915 + """Test calculation of t_test with products.""" + products = { + PreprocessorFile(regular_cubes, "A", {"dataset": "a"}), + PreprocessorFile(regular_cubes, "B", {"dataset": "b"}), + } + + if use_reference_product: + ref_product = PreprocessorFile( + ref_cubes, + "REF", + {"reference_for_t_test": True}, + ) + products.add(ref_product) + expected_ancestors = {ref_product} + reference = None + else: + expected_ancestors = set() + reference = ref_cubes[0] + + out_products = t_test(products, reference=reference) + + expected_ancillary_variable = AncillaryVariable( + np.array( + [[1.0, 0.6666667], [0.42264974, 0.46547753]], + dtype=np.float32, + ), + long_name="p-value", + var_name="pvalue", + units="1", + ) + + assert isinstance(out_products, set) + out_dict = products_set_to_dict(out_products) + assert len(out_dict) == len(products) + + product_a = out_dict["A"] + assert product_a.filename == "A" + assert product_a.attributes == {"dataset": "a"} + assert len(product_a.cubes) == 1 + out_cube = product_a.cubes[0] + assert out_cube.dtype == np.float32 + assert_allclose(out_cube.data, regular_cubes[0].data) + assert out_cube.var_name == "tas" + assert out_cube.standard_name == "air_temperature" + assert out_cube.units == "K" + assert out_cube.dim_coords == regular_cubes[0].dim_coords + assert out_cube.aux_coords == regular_cubes[0].aux_coords + assert out_cube.ancillary_variables() == [expected_ancillary_variable] + assert product_a.wasderivedfrom.call_count == len(expected_ancestors) + assert product_a.mock_ancestors == expected_ancestors + + product_b = out_dict["B"] + assert product_b.filename == "B" + assert product_b.attributes == {"dataset": "b"} + assert len(product_b.cubes) == 1 + out_cube = product_b.cubes[0] + assert out_cube.dtype == np.float32 + assert_allclose(out_cube.data, regular_cubes[0].data) + assert out_cube.var_name == "tas" + assert out_cube.standard_name == "air_temperature" + assert out_cube.units == "K" + assert out_cube.dim_coords == regular_cubes[0].dim_coords + assert out_cube.aux_coords == regular_cubes[0].aux_coords + assert out_cube.ancillary_variables() == [expected_ancillary_variable] + assert product_b.wasderivedfrom.call_count == len(expected_ancestors) + assert product_b.mock_ancestors == expected_ancestors + + if use_reference_product: + product_ref = out_dict["REF"] + assert product_ref.filename == "REF" + assert product_ref.attributes == {"reference_for_t_test": True} + assert len(product_ref.cubes) == 1 + out_cube = product_ref.cubes[0] + assert out_cube.dtype == np.float32 + assert_allclose(out_cube.data, ref_cubes[0].data) + assert out_cube.var_name == "tas" + assert out_cube.standard_name == "air_temperature" + assert out_cube.units == "K" + assert out_cube.dim_coords == ref_cubes[0].dim_coords + assert out_cube.aux_coords == ref_cubes[0].aux_coords + assert out_cube.ancillary_variables() == [] + product_ref.wasderivedfrom.assert_not_called() + assert product_ref.mock_ancestors == set() + + +def test_t_test_reference_none_cubes(regular_cubes): + """Test t_test with reference=None and cubes given.""" + msg = ( + r"A list of Cubes is given to this preprocessor; please specify a " + r"`reference`" + ) + with pytest.raises(ValueError, match=re.escape(msg)): + t_test(regular_cubes) + + +def test_no_reference_for_t_test(regular_cubes, ref_cubes): + """Test fail when no reference_for_t_test is given.""" + products = { + PreprocessorFile(regular_cubes, "A", {}), + PreprocessorFile(regular_cubes, "B", {}), + PreprocessorFile(ref_cubes, "REF", {}), + } + msg = r"Expected exactly 1 dataset with 'reference_for_t_test: true', found 0" + with pytest.raises(ValueError, match=re.escape(msg)): + t_test(products) + + +def test_two_references_for_t_test(regular_cubes, ref_cubes): + """Test fail when two reference_for_t_test products are given.""" + products = { + PreprocessorFile(regular_cubes, "A", {"reference_for_t_test": False}), + PreprocessorFile(ref_cubes, "REF1", {"reference_for_t_test": True}), + PreprocessorFile(ref_cubes, "REF2", {"reference_for_t_test": True}), + } + msg = r"Expected exactly 1 dataset with 'reference_for_t_test: true', found 2" + with pytest.raises(ValueError, match=re.escape(msg)): + t_test(products) From 93ee1342c18f7ad0e77135ba411d70f5eae5149d Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 11:16:12 +0200 Subject: [PATCH 13/33] Remove duplicate bias test --- .../test_compare_with_refs.py | 40 ++++++------------- 1 file changed, 12 insertions(+), 28 deletions(-) diff --git a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py index 5271305e63..51ab2bdb15 100644 --- a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py +++ b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py @@ -105,36 +105,20 @@ def ref_cubes(): ] -@pytest.mark.parametrize("use_reference_product", [True, False]) @pytest.mark.parametrize(("bias_type", "data", "units"), TEST_BIAS) -def test_bias_products( - regular_cubes, - ref_cubes, - bias_type, - data, - units, - use_reference_product, -): +def test_bias_products(regular_cubes, ref_cubes, bias_type, data, units): """Test calculation of bias with products.""" + ref_product = PreprocessorFile( + ref_cubes, + "REF", + {"reference_for_bias": True}, + ) products = { PreprocessorFile(regular_cubes, "A", {"dataset": "a"}), PreprocessorFile(regular_cubes, "B", {"dataset": "b"}), + ref_product, } - - if use_reference_product: - ref_product = PreprocessorFile( - ref_cubes, - "REF", - {"reference_for_bias": True}, - ) - products.add(ref_product) - expected_ancestors = {ref_product} - reference = None - else: - expected_ancestors = set() - reference = ref_cubes[0] - - out_products = bias(products, reference=reference, bias_type=bias_type) + out_products = bias(products, bias_type=bias_type) assert isinstance(out_products, set) out_dict = products_set_to_dict(out_products) @@ -152,8 +136,8 @@ def test_bias_products( assert out_cube.units == units assert out_cube.dim_coords == regular_cubes[0].dim_coords assert out_cube.aux_coords == regular_cubes[0].aux_coords - assert product_a.wasderivedfrom.call_count == len(expected_ancestors) - assert product_a.mock_ancestors == expected_ancestors + product_a.wasderivedfrom.assert_called_once() + assert product_a.mock_ancestors == {ref_product} product_b = out_dict["B"] assert product_b.filename == "B" @@ -167,8 +151,8 @@ def test_bias_products( assert out_cube.units == units assert out_cube.dim_coords == regular_cubes[0].dim_coords assert out_cube.aux_coords == regular_cubes[0].aux_coords - assert product_b.wasderivedfrom.call_count == len(expected_ancestors) - assert product_b.mock_ancestors == expected_ancestors + product_b.wasderivedfrom.assert_called_once() + assert product_b.mock_ancestors == {ref_product} @pytest.mark.parametrize(("bias_type", "data", "units"), TEST_BIAS) From e5ebef516b41baf1014227080743693ec8689a45 Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 11:54:57 +0200 Subject: [PATCH 14/33] Fix lazy t-test calculation --- esmvalcore/preprocessor/_compare_with_refs.py | 32 +++++++++++++------ 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index 2035088734..90c7764fe1 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -8,6 +8,7 @@ import dask import dask.array as da +import dask.array.stats import iris.analysis import iris.analysis.stats import numpy as np @@ -738,24 +739,37 @@ def _calculate_t_test( if coordinate is None: axis = None remaining_axes = () + shape = () else: axis = cube.coord_dims(coordinate)[0] remaining_axes = tuple(sorted(set(range(cube.ndim)) - {axis})) + shape = tuple(cube.shape[a] for a in remaining_axes) if cube.has_lazy_data() and reference.has_lazy_data(): - ttest_func = da.stats.ttest_ind + t_test = dask.array.stats.ttest_ind( + cube.core_data(), + reference.core_data(), + axis=axis, + **kwargs, + ) + # For some reason, t_test.pvalue is not a da.array, but a dask.Delayed + p_value_arr = da.from_delayed(t_test.pvalue, shape, cube.dtype) else: - ttest_func = ttest_ind - t_test = ttest_func( - cube.core_data(), - reference.core_data(), - axis=axis, - **kwargs, - ) + # Avoid realizing cube.data + cube_data = ( + cube.lazy_data().compute() if cube.has_lazy_data() else cube.data + ) + ref_data = ( + reference.lazy_data().compute() + if reference.has_lazy_data() + else reference.data + ) + t_test = ttest_ind(cube_data, ref_data, axis=axis, **kwargs) + p_value_arr = t_test.pvalue cube.add_ancillary_variable( AncillaryVariable( - t_test.pvalue, + p_value_arr, long_name="p-value", var_name="pvalue", units="1", From 0d1f95cd1a6de336c44a4a7ddc219880b558c9c1 Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 11:55:11 +0200 Subject: [PATCH 15/33] Add tests for t-test calculation of lazy and non-lazy cubes --- .../test_compare_with_refs.py | 97 ++++++++++++++++--- 1 file changed, 81 insertions(+), 16 deletions(-) diff --git a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py index 51ab2bdb15..4bb9185f1a 100644 --- a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py +++ b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py @@ -8,7 +8,7 @@ import numpy as np import pytest from cf_units import Unit -from iris.coords import AncillaryVariable, CellMeasure, CellMethod +from iris.coords import CellMeasure, CellMethod from iris.cube import Cube, CubeList from iris.exceptions import CoordinateNotFoundError @@ -939,6 +939,22 @@ def test_distance_metric_no_lon_for_area_weights(regular_cubes, metric, error): ) +def assert_correct_p_value(cube: Cube, *, is_lazy: bool = False) -> None: + assert len(cube.ancillary_variables()) == 1 + p_value = cube.ancillary_variables()[0] + assert p_value.standard_name is None + assert p_value.long_name == "p-value" + assert p_value.var_name == "pvalue" + assert p_value.units == "1" + assert p_value.attributes == {} + assert p_value.has_lazy_data() is is_lazy + assert p_value.dtype == np.float32 + np.testing.assert_allclose( + p_value.data, + [[1.0, 0.6666667], [0.42264974, 0.46547753]], + ) + + @pytest.mark.parametrize("use_reference_product", [True, False]) def test_t_test_products(regular_cubes, ref_cubes, use_reference_product): # noqa: PLR0915 """Test calculation of t_test with products.""" @@ -962,16 +978,6 @@ def test_t_test_products(regular_cubes, ref_cubes, use_reference_product): # no out_products = t_test(products, reference=reference) - expected_ancillary_variable = AncillaryVariable( - np.array( - [[1.0, 0.6666667], [0.42264974, 0.46547753]], - dtype=np.float32, - ), - long_name="p-value", - var_name="pvalue", - units="1", - ) - assert isinstance(out_products, set) out_dict = products_set_to_dict(out_products) assert len(out_dict) == len(products) @@ -988,7 +994,7 @@ def test_t_test_products(regular_cubes, ref_cubes, use_reference_product): # no assert out_cube.units == "K" assert out_cube.dim_coords == regular_cubes[0].dim_coords assert out_cube.aux_coords == regular_cubes[0].aux_coords - assert out_cube.ancillary_variables() == [expected_ancillary_variable] + assert_correct_p_value(out_cube) assert product_a.wasderivedfrom.call_count == len(expected_ancestors) assert product_a.mock_ancestors == expected_ancestors @@ -1004,7 +1010,7 @@ def test_t_test_products(regular_cubes, ref_cubes, use_reference_product): # no assert out_cube.units == "K" assert out_cube.dim_coords == regular_cubes[0].dim_coords assert out_cube.aux_coords == regular_cubes[0].aux_coords - assert out_cube.ancillary_variables() == [expected_ancillary_variable] + assert_correct_p_value(out_cube) assert product_b.wasderivedfrom.call_count == len(expected_ancestors) assert product_b.mock_ancestors == expected_ancestors @@ -1026,7 +1032,66 @@ def test_t_test_products(regular_cubes, ref_cubes, use_reference_product): # no assert product_ref.mock_ancestors == set() -def test_t_test_reference_none_cubes(regular_cubes): +@pytest.mark.parametrize("lazy", [True, False]) +def test_t_test_cubes(regular_cubes, ref_cubes, lazy): + """Test calculation of t_test with cubes.""" + cube = regular_cubes[0] + ref_cube = ref_cubes[0] + if lazy: + cube.data = cube.lazy_data() + ref_cube.data = ref_cube.lazy_data() + + out_cubes = t_test([cube], ref_cube) + + assert cube.has_lazy_data() is lazy + assert ref_cube.has_lazy_data() is lazy + + assert isinstance(out_cubes, CubeList) + assert len(out_cubes) == 1 + out_cube = out_cubes[0] + + assert out_cube.has_lazy_data() is lazy + assert out_cube.dtype == np.float32 + assert_allclose(out_cube.data, cube.data) + assert out_cube.var_name == "tas" + assert out_cube.standard_name == "air_temperature" + assert out_cube.units == "K" + assert out_cube.dim_coords == regular_cubes[0].dim_coords + assert out_cube.aux_coords == regular_cubes[0].aux_coords + assert_correct_p_value(out_cube, is_lazy=lazy) + + +@pytest.mark.parametrize("ref_lazy", [True, False]) +def test_t_test_cubes_partly_lazy(regular_cubes, ref_cubes, ref_lazy): + """Test calculation of t_test with cubes.""" + cube = regular_cubes[0] + ref_cube = ref_cubes[0] + if ref_lazy: + ref_cube.data = ref_cube.lazy_data() + else: + cube.data = cube.lazy_data() + + out_cubes = t_test([cube], ref_cube) + + assert cube.has_lazy_data() is not ref_lazy + assert ref_cube.has_lazy_data() is ref_lazy + + assert isinstance(out_cubes, CubeList) + assert len(out_cubes) == 1 + out_cube = out_cubes[0] + + assert out_cube.has_lazy_data() is not ref_lazy + assert out_cube.dtype == np.float32 + assert_allclose(out_cube.data, cube.data) + assert out_cube.var_name == "tas" + assert out_cube.standard_name == "air_temperature" + assert out_cube.units == "K" + assert out_cube.dim_coords == regular_cubes[0].dim_coords + assert out_cube.aux_coords == regular_cubes[0].aux_coords + assert_correct_p_value(out_cube) + + +def test_t_test_reference_none_cubes_fail(regular_cubes): """Test t_test with reference=None and cubes given.""" msg = ( r"A list of Cubes is given to this preprocessor; please specify a " @@ -1036,7 +1101,7 @@ def test_t_test_reference_none_cubes(regular_cubes): t_test(regular_cubes) -def test_no_reference_for_t_test(regular_cubes, ref_cubes): +def test_no_reference_for_t_test_fail(regular_cubes, ref_cubes): """Test fail when no reference_for_t_test is given.""" products = { PreprocessorFile(regular_cubes, "A", {}), @@ -1048,7 +1113,7 @@ def test_no_reference_for_t_test(regular_cubes, ref_cubes): t_test(products) -def test_two_references_for_t_test(regular_cubes, ref_cubes): +def test_two_references_for_t_test_fail(regular_cubes, ref_cubes): """Test fail when two reference_for_t_test products are given.""" products = { PreprocessorFile(regular_cubes, "A", {"reference_for_t_test": False}), From e70a9438c015fdd67aedd11c53f6f3681afbcead Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 11:58:30 +0200 Subject: [PATCH 16/33] t_test intrinsically preserves metadata --- esmvalcore/preprocessor/_compare_with_refs.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index 90c7764fe1..1ba3069509 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -733,7 +733,6 @@ def _calculate_t_test( ) -> Cube: """Calculate the t-test and attach the p-value as ancillary variable to cube.""" cube = cube.copy() # do not modify input cube - cube_metadata = cube.metadata remaining_axes: tuple[int, ...] if coordinate is None: @@ -776,6 +775,5 @@ def _calculate_t_test( ), remaining_axes, ) - cube.metadata = cube_metadata return cube From aedd6a621dfe7bdedff4b87605bf6944a49ec68f Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 13:39:28 +0200 Subject: [PATCH 17/33] Support masked data --- esmvalcore/preprocessor/_compare_with_refs.py | 25 ++++-- .../test_compare_with_refs.py | 85 +++++++++++++++++-- 2 files changed, 94 insertions(+), 16 deletions(-) diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index 1ba3069509..b6735fa43a 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -746,25 +746,36 @@ def _calculate_t_test( if cube.has_lazy_data() and reference.has_lazy_data(): t_test = dask.array.stats.ttest_ind( - cube.core_data(), - reference.core_data(), + da.ma.filled(cube.core_data(), np.nan), + da.ma.filled(reference.core_data(), np.nan), axis=axis, **kwargs, ) # For some reason, t_test.pvalue is not a da.array, but a dask.Delayed p_value_arr = da.from_delayed(t_test.pvalue, shape, cube.dtype) + p_value_arr = da.ma.masked_invalid(p_value_arr) else: # Avoid realizing cube.data - cube_data = ( - cube.lazy_data().compute() if cube.has_lazy_data() else cube.data + cube_data = np.ma.filled( + cube.lazy_data().compute() if cube.has_lazy_data() else cube.data, + np.nan, ) - ref_data = ( + ref_data = np.ma.filled( reference.lazy_data().compute() if reference.has_lazy_data() - else reference.data + else reference.data, + np.nan, ) t_test = ttest_ind(cube_data, ref_data, axis=axis, **kwargs) - p_value_arr = t_test.pvalue + p_value_arr = np.ma.masked_invalid(t_test.pvalue) + + # Ensure that float dtype from parents is inherited + if ( + np.issubdtype(cube.dtype, np.floating) + and np.issubdtype(reference.dtype, np.floating) + and cube.dtype == reference.dtype + ): + p_value_arr = p_value_arr.astype(cube.dtype) cube.add_ancillary_variable( AncillaryVariable( diff --git a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py index 4bb9185f1a..49660d4038 100644 --- a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py +++ b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py @@ -1,7 +1,10 @@ """Unit tests for :mod:`esmvalcore.preprocessor._compare_with_refs`.""" +from __future__ import annotations + import contextlib import re +from typing import Any import dask.array as da import iris @@ -939,7 +942,12 @@ def test_distance_metric_no_lon_for_area_weights(regular_cubes, metric, error): ) -def assert_correct_p_value(cube: Cube, *, is_lazy: bool = False) -> None: +def assert_correct_p_value( + cube: Cube, + data: Any, # noqa: ANN401 + *, + is_lazy: bool = False, +) -> None: assert len(cube.ancillary_variables()) == 1 p_value = cube.ancillary_variables()[0] assert p_value.standard_name is None @@ -949,10 +957,9 @@ def assert_correct_p_value(cube: Cube, *, is_lazy: bool = False) -> None: assert p_value.attributes == {} assert p_value.has_lazy_data() is is_lazy assert p_value.dtype == np.float32 - np.testing.assert_allclose( - p_value.data, - [[1.0, 0.6666667], [0.42264974, 0.46547753]], - ) + if np.ma.is_masked(data): + np.testing.assert_equal(p_value.data.mask, data.mask) + np.testing.assert_allclose(p_value.data, data) @pytest.mark.parametrize("use_reference_product", [True, False]) @@ -994,7 +1001,10 @@ def test_t_test_products(regular_cubes, ref_cubes, use_reference_product): # no assert out_cube.units == "K" assert out_cube.dim_coords == regular_cubes[0].dim_coords assert out_cube.aux_coords == regular_cubes[0].aux_coords - assert_correct_p_value(out_cube) + assert_correct_p_value( + out_cube, + [[1.0, 0.6666667], [0.42264974, 0.46547753]], + ) assert product_a.wasderivedfrom.call_count == len(expected_ancestors) assert product_a.mock_ancestors == expected_ancestors @@ -1010,7 +1020,10 @@ def test_t_test_products(regular_cubes, ref_cubes, use_reference_product): # no assert out_cube.units == "K" assert out_cube.dim_coords == regular_cubes[0].dim_coords assert out_cube.aux_coords == regular_cubes[0].aux_coords - assert_correct_p_value(out_cube) + assert_correct_p_value( + out_cube, + [[1.0, 0.6666667], [0.42264974, 0.46547753]], + ) assert product_b.wasderivedfrom.call_count == len(expected_ancestors) assert product_b.mock_ancestors == expected_ancestors @@ -1058,7 +1071,11 @@ def test_t_test_cubes(regular_cubes, ref_cubes, lazy): assert out_cube.units == "K" assert out_cube.dim_coords == regular_cubes[0].dim_coords assert out_cube.aux_coords == regular_cubes[0].aux_coords - assert_correct_p_value(out_cube, is_lazy=lazy) + assert_correct_p_value( + out_cube, + [[1.0, 0.6666667], [0.42264974, 0.46547753]], + is_lazy=lazy, + ) @pytest.mark.parametrize("ref_lazy", [True, False]) @@ -1088,7 +1105,51 @@ def test_t_test_cubes_partly_lazy(regular_cubes, ref_cubes, ref_lazy): assert out_cube.units == "K" assert out_cube.dim_coords == regular_cubes[0].dim_coords assert out_cube.aux_coords == regular_cubes[0].aux_coords - assert_correct_p_value(out_cube) + assert_correct_p_value( + out_cube, + [[1.0, 0.6666667], [0.42264974, 0.46547753]], + ) + + +@pytest.mark.parametrize("lazy", [True, False]) +def test_t_test_cubes_masked(regular_cubes, ref_cubes, lazy): + """Test calculation of t_test with cubes (masked data).""" + cube = regular_cubes[0] + cube.data = np.ma.masked_inside(cube.data, 1.5, 3.5) + cube.data = np.ma.masked_greater(cube.data, 6.5) + ref_cube = ref_cubes[0] + ref_cube.data = np.ma.masked_invalid( + np.array( + [[[2.0, np.nan], [2.0, 2.0]], [[2.0, 2.0], [2.0, 2.0]]], + dtype=np.float32, + ), + ) + if lazy: + cube.data = cube.lazy_data() + ref_cube.data = ref_cube.lazy_data() + + out_cubes = t_test([cube], ref_cube) + + assert cube.has_lazy_data() is lazy + assert ref_cube.has_lazy_data() is lazy + + assert isinstance(out_cubes, CubeList) + assert len(out_cubes) == 1 + out_cube = out_cubes[0] + + assert out_cube.has_lazy_data() is lazy + assert out_cube.dtype == np.float32 + assert_allclose(out_cube.data, cube.data) + assert out_cube.var_name == "tas" + assert out_cube.standard_name == "air_temperature" + assert out_cube.units == "K" + assert out_cube.dim_coords == regular_cubes[0].dim_coords + assert out_cube.aux_coords == regular_cubes[0].aux_coords + assert_correct_p_value( + out_cube, + np.ma.masked_invalid([[1.0, np.nan], [np.nan, np.nan]]), + is_lazy=lazy, + ) def test_t_test_reference_none_cubes_fail(regular_cubes): @@ -1123,3 +1184,9 @@ def test_two_references_for_t_test_fail(regular_cubes, ref_cubes): msg = r"Expected exactly 1 dataset with 'reference_for_t_test: true', found 2" with pytest.raises(ValueError, match=re.escape(msg)): t_test(products) + + +# ADD: +# - masked data +# - coordinate=None +# - **kwargs From 2066b5e94b102c34d8d26be94ca49ea9cbc2b56e Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 14:29:51 +0200 Subject: [PATCH 18/33] dask ttest_ind does not support axis=None, so we will not support it either --- esmvalcore/preprocessor/_compare_with_refs.py | 24 +++++++++---------- .../test_compare_with_refs.py | 18 ++++++++++---- 2 files changed, 25 insertions(+), 17 deletions(-) diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index b6735fa43a..69d931727a 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -15,6 +15,7 @@ from iris.common.metadata import CubeMetadata from iris.coords import AncillaryVariable, CellMethod from iris.cube import Cube, CubeList +from iris.exceptions import CoordinateMultiDimError from scipy.stats import ttest_ind, wasserstein_distance from esmvalcore.iris_helpers import ( @@ -617,7 +618,7 @@ def _get_emd( def t_test( products: set[PreprocessorFile] | Iterable[Cube], reference: Cube | None = None, - coordinate: str | AuxCoord | DimCoord | None = "time", + coordinate: str | AuxCoord | DimCoord = "time", **kwargs: Any, ) -> set[PreprocessorFile] | CubeList: """Perform t-test between dataset and reference dataset. @@ -652,8 +653,7 @@ def t_test( dataset in `products` needs the facet ``reference_for_t_test: true``. Do not specify this argument in a recipe. coordinate: - Coordinate axis of the input along which to compute the statistic. If - ``None``, the input will be raveled before computing the statistic. + Coordinate axis of the input along which to compute the statistic. **kwargs: Additional keyword arguments passed to :func:`scipy.stats.ttest_ind` (not all input data are lazy) or @@ -675,6 +675,8 @@ def t_test( ``reference_for_t_test: true`` if ``reference=None``; ``reference=None`` and the input products are given as iterable of :class:`~iris.cube.Cube` objects. + CoordinateMultiDimError + ``coordinate`` is not 1D. """ ref_product = None @@ -728,21 +730,17 @@ def t_test( def _calculate_t_test( cube: Cube, reference: Cube, - coordinate: str | AuxCoord | DimCoord | None, + coordinate: str | AuxCoord | DimCoord, **kwargs: Any, ) -> Cube: """Calculate the t-test and attach the p-value as ancillary variable to cube.""" cube = cube.copy() # do not modify input cube - remaining_axes: tuple[int, ...] - if coordinate is None: - axis = None - remaining_axes = () - shape = () - else: - axis = cube.coord_dims(coordinate)[0] - remaining_axes = tuple(sorted(set(range(cube.ndim)) - {axis})) - shape = tuple(cube.shape[a] for a in remaining_axes) + if cube.coord(coordinate).ndim > 1: + raise CoordinateMultiDimError(cube.coord(coordinate)) + axis = cube.coord_dims(coordinate)[0] + remaining_axes = tuple(sorted(set(range(cube.ndim)) - {axis})) + shape = tuple(cube.shape[a] for a in remaining_axes) if cube.has_lazy_data() and reference.has_lazy_data(): t_test = dask.array.stats.ttest_ind( diff --git a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py index 49660d4038..2c27580bfb 100644 --- a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py +++ b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py @@ -11,9 +11,9 @@ import numpy as np import pytest from cf_units import Unit -from iris.coords import CellMeasure, CellMethod +from iris.coords import AuxCoord, CellMeasure, CellMethod from iris.cube import Cube, CubeList -from iris.exceptions import CoordinateNotFoundError +from iris.exceptions import CoordinateMultiDimError, CoordinateNotFoundError from esmvalcore.preprocessor._compare_with_refs import ( bias, @@ -957,6 +957,7 @@ def assert_correct_p_value( assert p_value.attributes == {} assert p_value.has_lazy_data() is is_lazy assert p_value.dtype == np.float32 + print(p_value.data) if np.ma.is_masked(data): np.testing.assert_equal(p_value.data.mask, data.mask) np.testing.assert_allclose(p_value.data, data) @@ -1152,6 +1153,16 @@ def test_t_test_cubes_masked(regular_cubes, ref_cubes, lazy): ) +def test_t_test_multidim_coord_fail(regular_cubes): + """Test t_test with multidimensional coordinate.""" + cube = regular_cubes[0] + x_coord = AuxCoord([[0, 0], [0, 0]], var_name="x") + cube.add_aux_coord(x_coord, (0, 1)) + msg = r"Multi-dimensional coordinate not supported: 'x'" + with pytest.raises(CoordinateMultiDimError, match=re.escape(msg)): + t_test([cube], reference=cube, coordinate=x_coord) + + def test_t_test_reference_none_cubes_fail(regular_cubes): """Test t_test with reference=None and cubes given.""" msg = ( @@ -1187,6 +1198,5 @@ def test_two_references_for_t_test_fail(regular_cubes, ref_cubes): # ADD: -# - masked data -# - coordinate=None +# - 1D case # - **kwargs From 1a1bae800f21ec7ca799144828a77ca9ae5d2bc4 Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 14:44:56 +0200 Subject: [PATCH 19/33] Add test case for 1D data --- .../test_compare_with_refs.py | 40 ++++++++++++++++--- 1 file changed, 34 insertions(+), 6 deletions(-) diff --git a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py index 2c27580bfb..1ebaf70994 100644 --- a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py +++ b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py @@ -957,7 +957,6 @@ def assert_correct_p_value( assert p_value.attributes == {} assert p_value.has_lazy_data() is is_lazy assert p_value.dtype == np.float32 - print(p_value.data) if np.ma.is_masked(data): np.testing.assert_equal(p_value.data.mask, data.mask) np.testing.assert_allclose(p_value.data, data) @@ -1079,6 +1078,40 @@ def test_t_test_cubes(regular_cubes, ref_cubes, lazy): ) +@pytest.mark.parametrize("lazy", [True, False]) +def test_t_test_cubes_1d(regular_cubes, ref_cubes, lazy): + """Test calculation of t_test with 1D cubes.""" + cube = regular_cubes[0][0, :, 0] + ref_cube = ref_cubes[0][0, :, 0] + if lazy: + cube.data = cube.lazy_data() + ref_cube.data = ref_cube.lazy_data() + + out_cubes = t_test( + [cube], + ref_cube, + coordinate="latitude", + equal_var=False, + ) + + assert cube.has_lazy_data() is lazy + assert ref_cube.has_lazy_data() is lazy + + assert isinstance(out_cubes, CubeList) + assert len(out_cubes) == 1 + out_cube = out_cubes[0] + + assert out_cube.has_lazy_data() is lazy + assert out_cube.dtype == np.float32 + assert_allclose(out_cube.data, cube.data) + assert out_cube.var_name == "tas" + assert out_cube.standard_name == "air_temperature" + assert out_cube.units == "K" + assert out_cube.dim_coords == cube.dim_coords + assert out_cube.aux_coords == cube.aux_coords + assert_correct_p_value(out_cube, [0.4999999701976776], is_lazy=lazy) + + @pytest.mark.parametrize("ref_lazy", [True, False]) def test_t_test_cubes_partly_lazy(regular_cubes, ref_cubes, ref_lazy): """Test calculation of t_test with cubes.""" @@ -1195,8 +1228,3 @@ def test_two_references_for_t_test_fail(regular_cubes, ref_cubes): msg = r"Expected exactly 1 dataset with 'reference_for_t_test: true', found 2" with pytest.raises(ValueError, match=re.escape(msg)): t_test(products) - - -# ADD: -# - 1D case -# - **kwargs From bba79ac2a8bd19ec2cc3b0cc3baf48bb6655a6de Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 14:57:29 +0200 Subject: [PATCH 20/33] Make sure that bias also preserves cell measures and add corresponding test --- esmvalcore/preprocessor/_compare_with_refs.py | 16 ++++---- .../test_compare_with_refs.py | 37 ++++++++++++++++++- 2 files changed, 45 insertions(+), 8 deletions(-) diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index 69d931727a..a94fe47402 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -198,19 +198,19 @@ def _get_ref( return (reference, ref_product) -def _calculate_bias( - cube: Cube, - reference: Cube, - bias_type: BiasType, -) -> Cube: +def _calculate_bias(cube: Cube, reference: Cube, bias_type: BiasType) -> Cube: """Calculate bias for a single cube relative to a reference cube.""" cube_metadata = cube.metadata - # save ancillary variables as they get removed in the bias calculation + # Save ancillary variables and cell measures as they get removed in the + # bias calculation ancillary_variables_and_dims = [ (a, cube.ancillary_variable_dims(a)) for a in cube.ancillary_variables() ] + cell_measures_and_dims = [ + (c, cube.cell_measure_dims(c)) for c in cube.cell_measures() + ] if bias_type == "absolute": cube = cube - reference @@ -228,9 +228,11 @@ def _calculate_bias( cube.metadata = cube_metadata cube.units = new_units - # Reattach ancillary variables + # Reattach ancillary variables and cell measures for ancillary_variable, dims in ancillary_variables_and_dims: cube.add_ancillary_variable(ancillary_variable, dims) + for cell_measure, dims in cell_measures_and_dims: + cube.add_cell_measure(cell_measure, dims) return cube diff --git a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py index 1ebaf70994..736058febc 100644 --- a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py +++ b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py @@ -11,7 +11,7 @@ import numpy as np import pytest from cf_units import Unit -from iris.coords import AuxCoord, CellMeasure, CellMethod +from iris.coords import AncillaryVariable, AuxCoord, CellMeasure, CellMethod from iris.cube import Cube, CubeList from iris.exceptions import CoordinateMultiDimError, CoordinateNotFoundError @@ -351,6 +351,41 @@ def test_bias_products_and_ref_cube( assert product_a.mock_ancestors == set() +@pytest.mark.parametrize(("bias_type", "data", "units"), TEST_BIAS) +def test_bias_cubes_preserve_metadata( + regular_cubes, + ref_cubes, + bias_type, + data, + units, +): + """Test calculation of bias with cubes.""" + cube = regular_cubes[0] + ancillary_var = AncillaryVariable([0, 1], var_name="a") + cell_measure = CellMeasure([0, 1], var_name="c") + cube.add_ancillary_variable(ancillary_var, 0) + cube.add_cell_measure(cell_measure, 1) + ref_cube = ref_cubes[0] + + out_cubes = bias(regular_cubes, ref_cube, bias_type=bias_type) + + assert isinstance(out_cubes, CubeList) + assert len(out_cubes) == 1 + out_cube = out_cubes[0] + + assert out_cube.dtype == np.float32 + assert_allclose(out_cube.data, data) + assert out_cube.var_name == "tas" + assert out_cube.standard_name == "air_temperature" + assert out_cube.units == units + assert out_cube.dim_coords == cube.dim_coords + assert out_cube.aux_coords == cube.aux_coords + assert out_cube.ancillary_variables() == [ancillary_var] + assert out_cube.ancillary_variable_dims(ancillary_var) == (0,) + assert out_cube.cell_measures() == [cell_measure] + assert out_cube.cell_measure_dims(cell_measure) == (1,) + + def test_no_reference_for_bias(regular_cubes, ref_cubes): """Test fail when no reference_for_bias is given.""" products = { From 04006d715b8cda4e4115e57d9455bae635df32e3 Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Wed, 1 Jul 2026 15:08:27 +0200 Subject: [PATCH 21/33] Improve docs --- esmvalcore/preprocessor/_compare_with_refs.py | 36 +++++++++---------- 1 file changed, 16 insertions(+), 20 deletions(-) diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index a94fe47402..3583a97961 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -628,45 +628,41 @@ def t_test( This is a test for the null hypothesis that 2 independent samples (dataset and reference dataset) have identical average (expected) values. - The reference dataset needs to be broadcastable to all input `products`. - This supports `iris' rich broadcasting abilities - `__. To ensure this, the preprocessors - :func:`esmvalcore.preprocessor.regrid` and/or - :func:`esmvalcore.preprocessor.regrid_time` might be helpful. + All input datasets need to have the same shape. To ensure this, the + preprocessors :func:`esmvalcore.preprocessor.regrid` might be helpful. Notes ----- - The reference dataset can be specified with the `reference` argument. If - `reference` is ``None``, exactly one input dataset in the `products` set - needs to have the facet ``reference_for_t_test: true`` defined in the - recipe. Please do **not** specify the option `reference` when using this + The reference dataset can be specified with the ``reference`` argument. If + ``reference`` is ``None``, exactly one input dataset in the ``products`` + set needs to have the facet ``reference_for_t_test: true`` defined in the + recipe. Please do **not** specify the option ``reference`` when using this preprocessor function in a recipe. Parameters ---------- products: - Input datasets/cubes for which the significance is calculated relative - to a reference dataset/cube. + Input datasets/cubes for which the t-test is performed relative to a + reference dataset/cube. reference: Cube which is used as reference for the t-test calculation. If ``None``, `products` needs to be a :obj:`set` of - `~esmvalcore.preprocessor.PreprocessorFile` objects and exactly one - dataset in `products` needs the facet ``reference_for_t_test: true``. - Do not specify this argument in a recipe. + :class:`~esmvalcore.preprocessor.PreprocessorFile` objects and exactly + one dataset in ``products`` needs the facet ``reference_for_t_test: + true``. Do not specify this argument in a recipe. coordinate: Coordinate axis of the input along which to compute the statistic. **kwargs: Additional keyword arguments passed to :func:`scipy.stats.ttest_ind` - (not all input data are lazy) or - :func:`dask.array.stats.ttest_ind.ttest_ind` (all input data are lazy). + (not all input data are lazy) or :func:`dask.array.stats.ttest_ind` + (all input data are lazy). Returns ------- : Output datasets/cubes. Will be a :obj:`set` of :class:`~esmvalcore.preprocessor.PreprocessorFile` objects if - `products` is also one, a :class:`~iris.cube.CubeList` otherwise. The + ``products`` is also one, a :class:`~iris.cube.CubeList` otherwise. The cubes contain an attached ancillary variable which is the p-value w.r.t. the chosen t-test. @@ -742,7 +738,7 @@ def _calculate_t_test( raise CoordinateMultiDimError(cube.coord(coordinate)) axis = cube.coord_dims(coordinate)[0] remaining_axes = tuple(sorted(set(range(cube.ndim)) - {axis})) - shape = tuple(cube.shape[a] for a in remaining_axes) + target_shape = tuple(cube.shape[a] for a in remaining_axes) if cube.has_lazy_data() and reference.has_lazy_data(): t_test = dask.array.stats.ttest_ind( @@ -752,7 +748,7 @@ def _calculate_t_test( **kwargs, ) # For some reason, t_test.pvalue is not a da.array, but a dask.Delayed - p_value_arr = da.from_delayed(t_test.pvalue, shape, cube.dtype) + p_value_arr = da.from_delayed(t_test.pvalue, target_shape, cube.dtype) p_value_arr = da.ma.masked_invalid(p_value_arr) else: # Avoid realizing cube.data From 9c7f297f3fd288cc36ff102d27bf5e8e5ff2ad68 Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Tue, 21 Jul 2026 11:04:11 +0200 Subject: [PATCH 22/33] Generalize t-test so it can handle multiple coordinates at once --- esmvalcore/preprocessor/_compare_with_refs.py | 105 ++++++++----- .../test_compare_with_refs.py | 146 ++++++++++++++---- 2 files changed, 185 insertions(+), 66 deletions(-) diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index 3583a97961..dc616afd94 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -3,6 +3,7 @@ from __future__ import annotations import logging +import string from functools import partial from typing import TYPE_CHECKING, Any, Literal @@ -15,7 +16,6 @@ from iris.common.metadata import CubeMetadata from iris.coords import AncillaryVariable, CellMethod from iris.cube import Cube, CubeList -from iris.exceptions import CoordinateMultiDimError from scipy.stats import ttest_ind, wasserstein_distance from esmvalcore.iris_helpers import ( @@ -35,7 +35,7 @@ if TYPE_CHECKING: from collections.abc import Iterable - from iris.coords import AuxCoord, Coord, DimCoord + from iris.coords import Coord from esmvalcore.preprocessor import PreprocessorFile @@ -620,7 +620,7 @@ def _get_emd( def t_test( products: set[PreprocessorFile] | Iterable[Cube], reference: Cube | None = None, - coordinate: str | AuxCoord | DimCoord = "time", + coords: Iterable[Coord] | Iterable[str] | None = None, **kwargs: Any, ) -> set[PreprocessorFile] | CubeList: """Perform t-test between dataset and reference dataset. @@ -631,6 +631,8 @@ def t_test( All input datasets need to have the same shape. To ensure this, the preprocessors :func:`esmvalcore.preprocessor.regrid` might be helpful. + Handles lazy data and masked data. + Notes ----- The reference dataset can be specified with the ``reference`` argument. If @@ -650,8 +652,9 @@ def t_test( :class:`~esmvalcore.preprocessor.PreprocessorFile` objects and exactly one dataset in ``products`` needs the facet ``reference_for_t_test: true``. Do not specify this argument in a recipe. - coordinate: - Coordinate axis of the input along which to compute the statistic. + coords: + Coordinates over which the t-test is calculated. If ``None``, calculate + the t-test over all coordinates, which results in a scalar cube. **kwargs: Additional keyword arguments passed to :func:`scipy.stats.ttest_ind` (not all input data are lazy) or :func:`dask.array.stats.ttest_ind` @@ -673,8 +676,6 @@ def t_test( ``reference_for_t_test: true`` if ``reference=None``; ``reference=None`` and the input products are given as iterable of :class:`~iris.cube.Cube` objects. - CoordinateMultiDimError - ``coordinate`` is not 1D. """ ref_product = None @@ -695,7 +696,7 @@ def t_test( # If input is an Iterable of Cube objects, calculate t-test for each element if all_cubes_given: cubes = [ - _calculate_t_test(c, reference, coordinate, **kwargs) + _calculate_t_test(c, reference, coords=coords, **kwargs) for c in products ] return CubeList(cubes) @@ -709,7 +710,7 @@ def t_test( cube = concatenate(product.cubes) # Calculate t-test - cube = _calculate_t_test(cube, reference, coordinate, **kwargs) + cube = _calculate_t_test(cube, reference, coords=coords, **kwargs) # Adapt metadata and provenance information if ref_product is not None: @@ -728,51 +729,60 @@ def t_test( def _calculate_t_test( cube: Cube, reference: Cube, - coordinate: str | AuxCoord | DimCoord, + *, + coords: Iterable[Coord] | Iterable[str] | None = None, **kwargs: Any, ) -> Cube: """Calculate the t-test and attach the p-value as ancillary variable to cube.""" cube = cube.copy() # do not modify input cube - if cube.coord(coordinate).ndim > 1: - raise CoordinateMultiDimError(cube.coord(coordinate)) - axis = cube.coord_dims(coordinate)[0] - remaining_axes = tuple(sorted(set(range(cube.ndim)) - {axis})) - target_shape = tuple(cube.shape[a] for a in remaining_axes) + # Ensure that data is not chunked along desired coords + coords = get_all_coords(cube, coords) + cube = rechunk_cube(cube, coords) + reference = rechunk_cube(reference, coords) + + axes = get_all_coord_dims(cube, coords) + n_axes = len(axes) if cube.has_lazy_data() and reference.has_lazy_data(): - t_test = dask.array.stats.ttest_ind( - da.ma.filled(cube.core_data(), np.nan), - da.ma.filled(reference.core_data(), np.nan), - axis=axis, + # da.apply_gufunc transposes the input array so that the axes given by + # the `axes` argument to da.apply_gufunc are the rightmost dimensions. + # Thus, we need to use `along_axes=(ndim-n_axes, ..., ndim-2, ndim-1)` + # for _get_pvalue_from_ttest_ind here. + axes_in_chunk = tuple(range(cube.ndim - n_axes, cube.ndim)) + + # The call signature depends also on the number of axes in `axes`, and + # will be (a,b,...)->(nbins) where a,b,... are the data dimensions that + # are collapsed, and nbins the number of bin centers + in_signature = f"({','.join(list(string.ascii_lowercase)[:n_axes])})" + p_value_arr = da.apply_gufunc( + _get_pvalue_from_ttest_ind, + f"{in_signature},{in_signature}->()", + cube.lazy_data(), + reference.lazy_data(), + axes=[axes, axes, ()], + output_dtypes=cube.dtype, + along_axes=axes_in_chunk, **kwargs, ) - # For some reason, t_test.pvalue is not a da.array, but a dask.Delayed - p_value_arr = da.from_delayed(t_test.pvalue, target_shape, cube.dtype) - p_value_arr = da.ma.masked_invalid(p_value_arr) else: # Avoid realizing cube.data - cube_data = np.ma.filled( - cube.lazy_data().compute() if cube.has_lazy_data() else cube.data, - np.nan, + cube_data = ( + cube.lazy_data().compute() if cube.has_lazy_data() else cube.data ) - ref_data = np.ma.filled( + ref_data = ( reference.lazy_data().compute() if reference.has_lazy_data() - else reference.data, - np.nan, + else reference.data + ) + p_value_arr = _get_pvalue_from_ttest_ind( + cube_data, + ref_data, + along_axes=axes, + **kwargs, ) - t_test = ttest_ind(cube_data, ref_data, axis=axis, **kwargs) - p_value_arr = np.ma.masked_invalid(t_test.pvalue) - - # Ensure that float dtype from parents is inherited - if ( - np.issubdtype(cube.dtype, np.floating) - and np.issubdtype(reference.dtype, np.floating) - and cube.dtype == reference.dtype - ): - p_value_arr = p_value_arr.astype(cube.dtype) + remaining_axes = tuple(sorted(set(range(cube.ndim)) - set(axes))) cube.add_ancillary_variable( AncillaryVariable( p_value_arr, @@ -784,3 +794,22 @@ def _calculate_t_test( ) return cube + + +@preserve_float_dtype +def _get_pvalue_from_ttest_ind( + arr: np.ndarray, + ref_arr: np.ndarray, + *, + along_axes: int | Iterable[int] | None, + **kwargs: Any, +) -> np.ndarray: + """Calculate :func:`scipy.stats.ttest_ind` and return p-value.""" + # To support masked arrays, first convert masked values to NaNs; + # afterwards, convert NaNs back to masked values + arr = np.ma.filled(arr, np.nan) + ref_arr = np.ma.filled(ref_arr, np.nan) + + t_test_result = ttest_ind(arr, ref_arr, axis=along_axes, **kwargs) + + return np.ma.masked_invalid(t_test_result.pvalue) diff --git a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py index 736058febc..c5ee399633 100644 --- a/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py +++ b/tests/unit/preprocessor/_compare_with_refs/test_compare_with_refs.py @@ -13,7 +13,7 @@ from cf_units import Unit from iris.coords import AncillaryVariable, AuxCoord, CellMeasure, CellMethod from iris.cube import Cube, CubeList -from iris.exceptions import CoordinateMultiDimError, CoordinateNotFoundError +from iris.exceptions import CoordinateNotFoundError from esmvalcore.preprocessor._compare_with_refs import ( bias, @@ -981,6 +981,7 @@ def assert_correct_p_value( cube: Cube, data: Any, # noqa: ANN401 *, + dims: tuple[int, ...], is_lazy: bool = False, ) -> None: assert len(cube.ancillary_variables()) == 1 @@ -992,14 +993,15 @@ def assert_correct_p_value( assert p_value.attributes == {} assert p_value.has_lazy_data() is is_lazy assert p_value.dtype == np.float32 + print("p-value:", p_value.data) if np.ma.is_masked(data): np.testing.assert_equal(p_value.data.mask, data.mask) np.testing.assert_allclose(p_value.data, data) + assert cube.ancillary_variable_dims(p_value) == dims @pytest.mark.parametrize("use_reference_product", [True, False]) def test_t_test_products(regular_cubes, ref_cubes, use_reference_product): # noqa: PLR0915 - """Test calculation of t_test with products.""" products = { PreprocessorFile(regular_cubes, "A", {"dataset": "a"}), PreprocessorFile(regular_cubes, "B", {"dataset": "b"}), @@ -1018,7 +1020,7 @@ def test_t_test_products(regular_cubes, ref_cubes, use_reference_product): # no expected_ancestors = set() reference = ref_cubes[0] - out_products = t_test(products, reference=reference) + out_products = t_test(products, reference=reference, coords=["time"]) assert isinstance(out_products, set) out_dict = products_set_to_dict(out_products) @@ -1039,6 +1041,7 @@ def test_t_test_products(regular_cubes, ref_cubes, use_reference_product): # no assert_correct_p_value( out_cube, [[1.0, 0.6666667], [0.42264974, 0.46547753]], + dims=(1, 2), ) assert product_a.wasderivedfrom.call_count == len(expected_ancestors) assert product_a.mock_ancestors == expected_ancestors @@ -1058,6 +1061,7 @@ def test_t_test_products(regular_cubes, ref_cubes, use_reference_product): # no assert_correct_p_value( out_cube, [[1.0, 0.6666667], [0.42264974, 0.46547753]], + dims=(1, 2), ) assert product_b.wasderivedfrom.call_count == len(expected_ancestors) assert product_b.mock_ancestors == expected_ancestors @@ -1080,16 +1084,17 @@ def test_t_test_products(regular_cubes, ref_cubes, use_reference_product): # no assert product_ref.mock_ancestors == set() +@pytest.mark.parametrize("coords_as_str", [True, False]) @pytest.mark.parametrize("lazy", [True, False]) -def test_t_test_cubes(regular_cubes, ref_cubes, lazy): - """Test calculation of t_test with cubes.""" +def test_t_test_cubes_1d(regular_cubes, ref_cubes, lazy, coords_as_str): cube = regular_cubes[0] ref_cube = ref_cubes[0] if lazy: - cube.data = cube.lazy_data() - ref_cube.data = ref_cube.lazy_data() + cube.data = cube.lazy_data().rechunk((1, 1, 1)) + ref_cube.data = ref_cube.lazy_data().rechunk((1, 1, 1)) + coords = ["time"] if coords_as_str else [cube.coord("time")] - out_cubes = t_test([cube], ref_cube) + out_cubes = t_test([cube], ref_cube, coords=coords) assert cube.has_lazy_data() is lazy assert ref_cube.has_lazy_data() is lazy @@ -1109,23 +1114,113 @@ def test_t_test_cubes(regular_cubes, ref_cubes, lazy): assert_correct_p_value( out_cube, [[1.0, 0.6666667], [0.42264974, 0.46547753]], + dims=(1, 2), is_lazy=lazy, ) +@pytest.mark.parametrize( + ("dims", "p_value", "out_dims"), + [ + ((0, 1), [0.4679941236972809, 0.3202061057090759], (2,)), + ((0, 2), [0.6890520453453064, 0.17230847477912903], (1,)), + ], +) +@pytest.mark.parametrize("lazy", [True, False]) +def test_t_test_cubes_2d( + regular_cubes, + ref_cubes, + lazy, + dims, + p_value, + out_dims, +): + cube = regular_cubes[0] + ref_cube = ref_cubes[0] + x_coord = AuxCoord([[0, 0], [0, 0]], var_name="x") + cube.add_aux_coord(x_coord, dims) + ref_cube.add_aux_coord(x_coord, dims) + if lazy: + cube.data = cube.lazy_data().rechunk((1, 1, 1)) + ref_cube.data = ref_cube.lazy_data().rechunk((1, 1, 1)) + + out_cubes = t_test([cube], reference=ref_cube, coords=["x"]) + + assert cube.has_lazy_data() is lazy + assert ref_cube.has_lazy_data() is lazy + + assert isinstance(out_cubes, CubeList) + assert len(out_cubes) == 1 + out_cube = out_cubes[0] + + assert out_cube.has_lazy_data() is lazy + assert out_cube.dtype == np.float32 + assert_allclose(out_cube.data, cube.data) + assert out_cube.var_name == "tas" + assert out_cube.standard_name == "air_temperature" + assert out_cube.units == "K" + assert out_cube.dim_coords == regular_cubes[0].dim_coords + assert out_cube.aux_coords == regular_cubes[0].aux_coords + assert_correct_p_value(out_cube, p_value, dims=out_dims, is_lazy=lazy) + + +@pytest.mark.parametrize( + ("coords", "p_value", "out_dims"), + [ + (("time", "latitude"), [0.4679941236972809, 0.3202061057090759], (2,)), + ( + ("time", "longitude"), + [0.6890520453453064, 0.17230847477912903], + (1,), + ), + ], +) @pytest.mark.parametrize("lazy", [True, False]) -def test_t_test_cubes_1d(regular_cubes, ref_cubes, lazy): - """Test calculation of t_test with 1D cubes.""" - cube = regular_cubes[0][0, :, 0] - ref_cube = ref_cubes[0][0, :, 0] +def test_t_test_cubes_multiple_coords( + regular_cubes, + ref_cubes, + lazy, + coords, + p_value, + out_dims, +): + cube = regular_cubes[0] + ref_cube = ref_cubes[0] if lazy: cube.data = cube.lazy_data() ref_cube.data = ref_cube.lazy_data() + out_cubes = t_test([cube], reference=ref_cube, coords=coords) + + assert cube.has_lazy_data() is lazy + assert ref_cube.has_lazy_data() is lazy + + assert isinstance(out_cubes, CubeList) + assert len(out_cubes) == 1 + out_cube = out_cubes[0] + + assert out_cube.has_lazy_data() is lazy + assert out_cube.dtype == np.float32 + assert_allclose(out_cube.data, cube.data) + assert out_cube.var_name == "tas" + assert out_cube.standard_name == "air_temperature" + assert out_cube.units == "K" + assert out_cube.dim_coords == regular_cubes[0].dim_coords + assert out_cube.aux_coords == regular_cubes[0].aux_coords + assert_correct_p_value(out_cube, p_value, dims=out_dims, is_lazy=lazy) + + +@pytest.mark.parametrize("lazy", [True, False]) +def test_t_test_cubes_scalar_result(regular_cubes, ref_cubes, lazy): + cube = regular_cubes[0] + ref_cube = ref_cubes[0] + if lazy: + cube.data = cube.lazy_data().rechunk((1, 1, 1)) + ref_cube.data = ref_cube.lazy_data().rechunk((1, 1, 1)) + out_cubes = t_test( [cube], ref_cube, - coordinate="latitude", equal_var=False, ) @@ -1144,12 +1239,16 @@ def test_t_test_cubes_1d(regular_cubes, ref_cubes, lazy): assert out_cube.units == "K" assert out_cube.dim_coords == cube.dim_coords assert out_cube.aux_coords == cube.aux_coords - assert_correct_p_value(out_cube, [0.4999999701976776], is_lazy=lazy) + assert_correct_p_value( + out_cube, + [0.20223067700862885], + dims=(), + is_lazy=lazy, + ) @pytest.mark.parametrize("ref_lazy", [True, False]) def test_t_test_cubes_partly_lazy(regular_cubes, ref_cubes, ref_lazy): - """Test calculation of t_test with cubes.""" cube = regular_cubes[0] ref_cube = ref_cubes[0] if ref_lazy: @@ -1157,7 +1256,7 @@ def test_t_test_cubes_partly_lazy(regular_cubes, ref_cubes, ref_lazy): else: cube.data = cube.lazy_data() - out_cubes = t_test([cube], ref_cube) + out_cubes = t_test([cube], ref_cube, coords=["time"]) assert cube.has_lazy_data() is not ref_lazy assert ref_cube.has_lazy_data() is ref_lazy @@ -1177,12 +1276,12 @@ def test_t_test_cubes_partly_lazy(regular_cubes, ref_cubes, ref_lazy): assert_correct_p_value( out_cube, [[1.0, 0.6666667], [0.42264974, 0.46547753]], + dims=(1, 2), ) @pytest.mark.parametrize("lazy", [True, False]) def test_t_test_cubes_masked(regular_cubes, ref_cubes, lazy): - """Test calculation of t_test with cubes (masked data).""" cube = regular_cubes[0] cube.data = np.ma.masked_inside(cube.data, 1.5, 3.5) cube.data = np.ma.masked_greater(cube.data, 6.5) @@ -1197,7 +1296,7 @@ def test_t_test_cubes_masked(regular_cubes, ref_cubes, lazy): cube.data = cube.lazy_data() ref_cube.data = ref_cube.lazy_data() - out_cubes = t_test([cube], ref_cube) + out_cubes = t_test([cube], ref_cube, coords=["time"]) assert cube.has_lazy_data() is lazy assert ref_cube.has_lazy_data() is lazy @@ -1217,20 +1316,11 @@ def test_t_test_cubes_masked(regular_cubes, ref_cubes, lazy): assert_correct_p_value( out_cube, np.ma.masked_invalid([[1.0, np.nan], [np.nan, np.nan]]), + dims=(1, 2), is_lazy=lazy, ) -def test_t_test_multidim_coord_fail(regular_cubes): - """Test t_test with multidimensional coordinate.""" - cube = regular_cubes[0] - x_coord = AuxCoord([[0, 0], [0, 0]], var_name="x") - cube.add_aux_coord(x_coord, (0, 1)) - msg = r"Multi-dimensional coordinate not supported: 'x'" - with pytest.raises(CoordinateMultiDimError, match=re.escape(msg)): - t_test([cube], reference=cube, coordinate=x_coord) - - def test_t_test_reference_none_cubes_fail(regular_cubes): """Test t_test with reference=None and cubes given.""" msg = ( From 32f57dd2e9ef0f209add73146d8822f49f27bd0b Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Tue, 21 Jul 2026 11:37:11 +0200 Subject: [PATCH 23/33] Add early check for presence of t-test reference in recipe --- esmvalcore/_recipe/check.py | 7 +++ esmvalcore/_recipe/recipe.py | 1 + tests/integration/recipe/test_recipe.py | 68 +++++++++++++++++++++++++ 3 files changed, 76 insertions(+) diff --git a/esmvalcore/_recipe/check.py b/esmvalcore/_recipe/check.py index bb935cbc0a..436c30b311 100644 --- a/esmvalcore/_recipe/check.py +++ b/esmvalcore/_recipe/check.py @@ -627,6 +627,13 @@ def _check_ref_attributes(products: set, *, step: str, attr_name: str) -> None: ) +reference_for_t_test_preproc = partial( + _check_ref_attributes, + step="t_test", + attr_name="reference_for_t_test", +) + + def statistics_preprocessors(settings: dict) -> None: """Check options of statistics preprocessors.""" mm_stats = ( diff --git a/esmvalcore/_recipe/recipe.py b/esmvalcore/_recipe/recipe.py index f15e61e92a..c974eda439 100644 --- a/esmvalcore/_recipe/recipe.py +++ b/esmvalcore/_recipe/recipe.py @@ -704,6 +704,7 @@ def _get_preprocessor_products( check.reference_for_bias_preproc(products) check.reference_for_distance_metric_preproc(products) + check.reference_for_t_test_preproc(products) _configure_multi_product_preprocessor( products=products, diff --git a/tests/integration/recipe/test_recipe.py b/tests/integration/recipe/test_recipe.py index 8dbe034e61..96dae3fe9a 100644 --- a/tests/integration/recipe/test_recipe.py +++ b/tests/integration/recipe/test_recipe.py @@ -3998,3 +3998,71 @@ def test_align_metadata_missing_arg(tmp_path, patched_datafinder, session): msg = "Missing required argument" with pytest.raises(ValueError, match=msg): get_recipe(tmp_path, content, session) + + +def test_t_test_no_ref(tmp_path, patched_datafinder, session): + content = dedent(""" + preprocessors: + test_t_test: + t_test: + + diagnostics: + diagnostic_name: + variables: + ta: + preprocessor: test_t_test + project: CMIP6 + mip: Amon + exp: historical + timerange: 20000101/20001231 + ensemble: r1i1p1f1 + grid: gn + additional_datasets: + - {dataset: CanESM5} + - {dataset: CESM2} + + scripts: null + """) + msg = ( + "Expected exactly 1 dataset with 'reference_for_t_test: true' in " + "products" + ) + with pytest.raises(RecipeError) as exc: + get_recipe(tmp_path, content, session) + assert str(exc.value) == INITIALIZATION_ERROR_MSG + assert msg in exc.value.failed_tasks[0].message + assert "found 0" in exc.value.failed_tasks[0].message + + +def test_t_test_two_refs(tmp_path, patched_datafinder, session): + content = dedent(""" + preprocessors: + test_t_test: + t_test: + + diagnostics: + diagnostic_name: + variables: + ta: + preprocessor: test_t_test + project: CMIP6 + mip: Amon + exp: historical + timerange: 20000101/20001231 + ensemble: r1i1p1f1 + grid: gn + additional_datasets: + - {dataset: CanESM5, reference_for_t_test: true} + - {dataset: CESM2, reference_for_t_test: true} + + scripts: null + """) + msg = ( + "Expected exactly 1 dataset with 'reference_for_t_test: true' in " + "products" + ) + with pytest.raises(RecipeError) as exc: + get_recipe(tmp_path, content, session) + assert str(exc.value) == INITIALIZATION_ERROR_MSG + assert msg in exc.value.failed_tasks[0].message + assert "found 2" in exc.value.failed_tasks[0].message From 980fdd6273dd36217d56177041cd73ba0af303ae Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Tue, 21 Jul 2026 12:07:22 +0200 Subject: [PATCH 24/33] Add docs for t-test calculation --- doc/recipe/preprocessor.rst | 91 +++++++++++++++++++++++++++++++++++-- 1 file changed, 88 insertions(+), 3 deletions(-) diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst index e64fad65f8..7acad0fb0a 100644 --- a/doc/recipe/preprocessor.rst +++ b/doc/recipe/preprocessor.rst @@ -2638,8 +2638,11 @@ This module contains the following preprocessor functions: * ``bias``: Calculate absolute or relative biases with respect to a reference dataset. -* ``distance_metric``: Calculate absolute or relative biases with respect to a - reference dataset. +* ``distance_metric``: Calculate distance metrics with respect to a reference + dataset. +* ``t_test``: Calculate t-test (null hypothesis: 2 independent samples have + identical mean) with respect to a reference dataset and attach the + corresponding p-value as ancillary variable to the data. ``bias`` -------- @@ -2937,8 +2940,90 @@ See also :func:`esmvalcore.preprocessor.distance_metric`. .. _Weighted Earth mover's distance: https://pythonot.github.io/ quickstart.html#computing-wasserstein-distance -.. _Other: +.. _t_test: + +``t_test`` +---------- + +This function performs a :func:`t-test ` (null +hypothesis: 2 independent samples have identical mean) with respect to a +reference dataset and attaches the corresponding p-value as ancillary variable +to the data. +For this, exactly one input dataset needs to be declared as +``reference_for_t_test: true`` in the recipe, e.g., + +.. code-block:: yaml + + datasets: + - {dataset: CanESM5, project: CMIP6, ensemble: r1i1p1f1, grid: gn} + - {dataset: CESM2, project: CMIP6, ensemble: r1i1p1f1, grid: gn} + - {dataset: MIROC6, project: CMIP6, ensemble: r1i1p1f1, grid: gn} + - {dataset: ERA-Interim, project: OBS6, tier: 3, type: reanaly, version: 1, + reference_for_t_test: true} + +In the example above, ERA-Interim is used as reference dataset for the t-test +calculation. + +It is also possible to use the output from the :ref:`multi-model statistics` or +:ref:`ensemble statistics` preprocessor as reference dataset. +In this case, make sure to use ``reference_for_metric: true`` for each dataset +that will be used to create the reference dataset and use the option +``keep_input_datasets: false`` for the multi-dataset preprocessor. +For example: + +.. code-block:: yaml + + datasets: + - {dataset: CanESM5, group: ref, reference_for_t_test: true} + - {dataset: CESM2, group: ref, reference_for_t_test: true} + - {dataset: MIROC6, group: notref} + preprocessors: + calculate_t_test: + custom_order: true + multi_model_statistics: + statistics: [mean] + span: overlap + groupby: [group] + keep_input_datasets: false + t_test: + +Here, the t-test for MIROC6 is calculated relative to the the multi-model mean +from the models CanESM5 and CESM2. + +All datasets need to have the same shape and coordinates. +To ensure this, the preprocessors :func:`esmvalcore.preprocessor.regrid` and/or +:func:`esmvalcore.preprocessor.regrid_time` might be helpful. + +The ``t_test`` preprocessor supports the following arguments in the recipe: + +* ``coords`` (:obj:`list` of :obj:`str`, default: ``None``): Coordinates over + which the t-test is calculated. + If ``None``, calculate the metric over all coordinates, which results in a + scalar cube. +* Other parameters are passed to :func:`scipy.stats.ttest_ind`. + +Example: + +.. code-block:: yaml + + preprocessors: + preproc_t_test: + t_test: + coords: [time] + equal_var: false + nan_policy: omit + +This will perform a `Welch's t-test +`__, which does not assume +equal variances between the samples, and omit all masked values for the t-test +calculation (the default behavior is to return masked values if at least one +input value is masked). + +See also :func:`esmvalcore.preprocessor.t_test`. + + +.. _Other: Other ===== From f8e65ec4ec40c57698ca96e6192959389ee53e21 Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Tue, 21 Jul 2026 12:19:04 +0200 Subject: [PATCH 25/33] Add example of p-value ancillary var --- doc/recipe/preprocessor.rst | 28 +++++++++++++++++++++++++++- 1 file changed, 27 insertions(+), 1 deletion(-) diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst index 7acad0fb0a..cc9e0f6f27 100644 --- a/doc/recipe/preprocessor.rst +++ b/doc/recipe/preprocessor.rst @@ -2949,7 +2949,33 @@ This function performs a :func:`t-test ` (null hypothesis: 2 independent samples have identical mean) with respect to a reference dataset and attaches the corresponding p-value as ancillary variable to the data. -For this, exactly one input dataset needs to be declared as + +With `Iris `__, the p-value can be accessed with the following code: + +.. code-block:: python + + import iris + cube = iris.load_cube("esmvaltool_output/recipe_t_test_20260721_095242/preproc/t-test/tas/CMIP6_BCC-ESM1_Amon_historical_r1i1p1f1_tas_gn_2000-2005.nc") + p_value = cube.ancillary_variable("p-value") + print(p_value) + +The returned p-value looks like this: + +.. code-block:: + + AncillaryVariable : p-value / (1) + data: [ + [0.21600565 , 0.18489039 , ..., 0.3673955 , 0.27788544 ], + [0.3407415 , 0.3905643 , ..., -- , 0.58700126 ], + ..., + [-- , -- , ..., -- , -- ], + [-- , -- , ..., -- , -- ]] + shape: (18, 36) + dtype: float32 + long_name: 'p-value' + var_name: 'pvalue' + +For this preprocessor, exactly one input dataset needs to be declared as ``reference_for_t_test: true`` in the recipe, e.g., .. code-block:: yaml From 69d2a9613d6d48e014e3633f6ddf6bc8ccce7dd9 Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Tue, 21 Jul 2026 12:19:58 +0200 Subject: [PATCH 26/33] metric -> t_test --- doc/recipe/preprocessor.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst index cc9e0f6f27..7ee6653461 100644 --- a/doc/recipe/preprocessor.rst +++ b/doc/recipe/preprocessor.rst @@ -2992,7 +2992,7 @@ calculation. It is also possible to use the output from the :ref:`multi-model statistics` or :ref:`ensemble statistics` preprocessor as reference dataset. -In this case, make sure to use ``reference_for_metric: true`` for each dataset +In this case, make sure to use ``reference_for_t_test: true`` for each dataset that will be used to create the reference dataset and use the option ``keep_input_datasets: false`` for the multi-dataset preprocessor. For example: From e6d1168497acbd56b01167f9aa416bbce51bed2d Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Tue, 21 Jul 2026 12:22:09 +0200 Subject: [PATCH 27/33] Improve docs --- doc/recipe/preprocessor.rst | 6 +++--- esmvalcore/preprocessor/_compare_with_refs.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst index 7ee6653461..bce797fec0 100644 --- a/doc/recipe/preprocessor.rst +++ b/doc/recipe/preprocessor.rst @@ -3017,9 +3017,9 @@ For example: Here, the t-test for MIROC6 is calculated relative to the the multi-model mean from the models CanESM5 and CESM2. -All datasets need to have the same shape and coordinates. -To ensure this, the preprocessors :func:`esmvalcore.preprocessor.regrid` and/or -:func:`esmvalcore.preprocessor.regrid_time` might be helpful. +All datasets need to have the same shape. +To ensure this, the preprocessors :func:`esmvalcore.preprocessor.regrid` might +be helpful. The ``t_test`` preprocessor supports the following arguments in the recipe: diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index dc616afd94..e1089f3254 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -629,7 +629,7 @@ def t_test( and reference dataset) have identical average (expected) values. All input datasets need to have the same shape. To ensure this, the - preprocessors :func:`esmvalcore.preprocessor.regrid` might be helpful. + preprocessors :func:`esmvalcore.preprocessor.regrid` might be helpful. Handles lazy data and masked data. From 1c73391a42665158af7f4810345d69152f92f472 Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Tue, 21 Jul 2026 12:22:54 +0200 Subject: [PATCH 28/33] Better docs --- doc/recipe/preprocessor.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst index bce797fec0..389de46db2 100644 --- a/doc/recipe/preprocessor.rst +++ b/doc/recipe/preprocessor.rst @@ -3034,7 +3034,7 @@ Example: .. code-block:: yaml preprocessors: - preproc_t_test: + welchs_t_test: t_test: coords: [time] equal_var: false From 33ace69b6879b165705efd44e4711cd00e55f0e8 Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Tue, 21 Jul 2026 12:23:58 +0200 Subject: [PATCH 29/33] Better example --- doc/recipe/preprocessor.rst | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst index 389de46db2..a083437bcd 100644 --- a/doc/recipe/preprocessor.rst +++ b/doc/recipe/preprocessor.rst @@ -3042,9 +3042,10 @@ Example: This will perform a `Welch's t-test `__, which does not assume -equal variances between the samples, and omit all masked values for the t-test -calculation (the default behavior is to return masked values if at least one -input value is masked). +equal variances between the samples +In addtion, all masked values will be omitted in the t-test calculation (the +default behavior is to return masked values if at least one input value is +masked). See also :func:`esmvalcore.preprocessor.t_test`. From 21862bcb4c878c12db63e9b774ec763ce83f4346 Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Tue, 21 Jul 2026 12:24:15 +0200 Subject: [PATCH 30/33] Typo --- doc/recipe/preprocessor.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst index a083437bcd..9604d8fb4c 100644 --- a/doc/recipe/preprocessor.rst +++ b/doc/recipe/preprocessor.rst @@ -3043,7 +3043,7 @@ Example: This will perform a `Welch's t-test `__, which does not assume equal variances between the samples -In addtion, all masked values will be omitted in the t-test calculation (the +In addition, all masked values will be omitted in the t-test calculation (the default behavior is to return masked values if at least one input value is masked). From 8fc32d036cf3ea03a2b862207059764c7d688a2d Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Tue, 21 Jul 2026 12:26:54 +0200 Subject: [PATCH 31/33] Better docstring --- esmvalcore/preprocessor/_compare_with_refs.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index e1089f3254..e78f7491b6 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -656,9 +656,7 @@ def t_test( Coordinates over which the t-test is calculated. If ``None``, calculate the t-test over all coordinates, which results in a scalar cube. **kwargs: - Additional keyword arguments passed to :func:`scipy.stats.ttest_ind` - (not all input data are lazy) or :func:`dask.array.stats.ttest_ind` - (all input data are lazy). + Additional keyword arguments passed to :func:`scipy.stats.ttest_ind`. Returns ------- @@ -667,7 +665,8 @@ def t_test( :class:`~esmvalcore.preprocessor.PreprocessorFile` objects if ``products`` is also one, a :class:`~iris.cube.CubeList` otherwise. The cubes contain an attached ancillary variable which is the p-value - w.r.t. the chosen t-test. + w.r.t. the chosen t-test which can be accessed with + ``cube.ancillary_variable("p-value")``. Raises ------ From dcbe055468039791fb2ced5d6f9422f10229fee2 Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Tue, 21 Jul 2026 12:31:56 +0200 Subject: [PATCH 32/33] Full stop missing --- doc/recipe/preprocessor.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/recipe/preprocessor.rst b/doc/recipe/preprocessor.rst index 9604d8fb4c..cb2747e520 100644 --- a/doc/recipe/preprocessor.rst +++ b/doc/recipe/preprocessor.rst @@ -3042,7 +3042,7 @@ Example: This will perform a `Welch's t-test `__, which does not assume -equal variances between the samples +equal variances between the samples. In addition, all masked values will be omitted in the t-test calculation (the default behavior is to return masked values if at least one input value is masked). From cff6911442e28215c3aeb23f37064bf18b55efbb Mon Sep 17 00:00:00 2001 From: Manuel Schlund Date: Tue, 21 Jul 2026 12:33:45 +0200 Subject: [PATCH 33/33] Better dosctring --- esmvalcore/preprocessor/_compare_with_refs.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/esmvalcore/preprocessor/_compare_with_refs.py b/esmvalcore/preprocessor/_compare_with_refs.py index e78f7491b6..b84350e81b 100644 --- a/esmvalcore/preprocessor/_compare_with_refs.py +++ b/esmvalcore/preprocessor/_compare_with_refs.py @@ -664,8 +664,8 @@ def t_test( Output datasets/cubes. Will be a :obj:`set` of :class:`~esmvalcore.preprocessor.PreprocessorFile` objects if ``products`` is also one, a :class:`~iris.cube.CubeList` otherwise. The - cubes contain an attached ancillary variable which is the p-value - w.r.t. the chosen t-test which can be accessed with + cubes contain an attached :class:`~iris.coords.AncillaryVariable` which + is the p-value w.r.t. the chosen t-test. It can be accessed with ``cube.ancillary_variable("p-value")``. Raises