You can run this notebook in a live session or view it on Github.
Implications of verify(dim)
#
This demo demonstrates how setting dim
in PerfectModelEnsemble.verify()
and HindcastEnsemble.verify()
alters the research question and prediction skill obtained, especially for spatial fields.
What’s the purpose of dim
?
dim
is designed the same way as in xarray.Dataset.mean()
. A given metric is calculated over dimension dim
, i.e. the result does not contain dim
anymore. dim
can be a string or list of strings.
3 ways to calculate aggregated prediction skill
Spatially aggregate first, then
verify
verify
on each grid point, then aggregate spatiallyverify
over spatial,member
andinit
ialization directly
# linting
%load_ext nb_black
%load_ext lab_black
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
import climpred
climpred.set_options(PerfectModel_persistence_from_initialized_lead_0=True)
/home/docs/checkouts/readthedocs.org/user_builds/climpred/conda/latest/lib/python3.9/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from .autonotebook import tqdm as notebook_tqdm
<climpred.options.set_options at 0x7effa4b92b80>
sample spatial initialized ensemble data and nomenclature#
We also have some sample output that contains gridded data on MPIOM grid:
initialized3d
: The initialized ensemble dataset withmember
s (1
, … ,4
) denoted byinit
s (initialization years:3014
,3061
,3175
,3237
) denoted bylead
years (1
, …,5
) denoted by.
x
andj
(think as longitude and latitude): spatial dimensions denotes by.
control3d
: The control dataset spanning overtime
(3000
, …,3049
) [not used here].
In PerfectModelEnsemble.verify()
, we calculate a on
containing spatial dimensions
and ensemble dimension
over dimensions
denoted by
.
# Sea surface temperature
initialized3d = climpred.tutorial.load_dataset("MPI-PM-DP-3D")
control3d = climpred.tutorial.load_dataset("MPI-control-3D")
v = "tos"
initialized3d["lead"].attrs = {"units": "years"}
/home/docs/checkouts/readthedocs.org/user_builds/climpred/conda/latest/lib/python3.9/site-packages/xarray/backends/plugins.py:71: RuntimeWarning: Engine 'cfgrib' loading failed:
Cannot find the ecCodes library
warnings.warn(f"Engine {name!r} loading failed:\n{ex}", RuntimeWarning)
# Create climpred PerfectModelEnsemble object.
pm = (
climpred.PerfectModelEnsemble(initialized3d).add_control(control3d)
# .generate_uninitialized() # needed for uninitialized in reference
)
/home/docs/checkouts/readthedocs.org/user_builds/climpred/checkouts/latest/climpred/utils.py:195: UserWarning: Assuming annual resolution starting Jan 1st due to numeric inits. Please change ``init`` to a datetime if it is another resolution. We recommend using xr.CFTimeIndex as ``init``, see https://climpred.readthedocs.io/en/stable/setting-up-data.html.
warnings.warn(
/home/docs/checkouts/readthedocs.org/user_builds/climpred/checkouts/latest/climpred/utils.py:195: UserWarning: Assuming annual resolution starting Jan 1st due to numeric inits. Please change ``init`` to a datetime if it is another resolution. We recommend using xr.CFTimeIndex as ``init``, see https://climpred.readthedocs.io/en/stable/setting-up-data.html.
warnings.warn(
subselecting North Atlantic#
pm = pm.isel(x=slice(80, 150), y=slice(40, 70))
pm.get_initialized().mean(["lead", "member", "init"])[v].plot(
yincrease=False, robust=True
)
<matplotlib.collections.QuadMesh at 0x7eff8f179a30>

verify_kwargs = dict(
metric="pearson_r", # 'rmse'
comparison="m2m",
reference=["persistence", "climatology"], # "uninitialized",
skipna=True,
)
dim = ["init", "member"]
spatial_dims = ["x", "y"]
Here, we use the anomaly correlation coefficient climpred.metrics._pearson_r()
. Also try root-mean-squared-error climpred.metrics._rmse()
.
Spatial aggregate first, then verify
#
Here, we first average the raw ensemble data over the spatial_dims
, i.e. creating North Atlantic averaged SST timeseries. On this timeseries, we calculate the prediction skill with PerfectModelEnsemble.verify()
.
where represents the average of
over dimension
This approach answers the question: “What’s the prediction skill of SSTs, which were averaged before over the North Atlantic? What is the skill at predicting average [metric] in [domain]?”
Used in Séférian et al. [2018], Spring and Ilyina [2020], Pohlmann et al. [2004].
Recommended by Goddard et al. [2013] section 2.1.4:
“We advocate verifying on at least two spatial scales: (1) the observational grid scales to which the model data is interpolated, and (2) smoothed or regional scales. The latter can be accomplished by suitable spatial-smoothing algorithms, such as simple averages or spectral filters.”
pm_aggregated = pm.mean(spatial_dims)
# see the NA SST timeseries
pm_aggregated.plot(show_members=True)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[7], line 3
1 pm_aggregated = pm.mean(spatial_dims)
2 # see the NA SST timeseries
----> 3 pm_aggregated.plot(show_members=True)
File ~/checkouts/readthedocs.org/user_builds/climpred/checkouts/latest/climpred/classes.py:486, in PredictionEnsemble.plot(self, variable, ax, show_members, cmap, x)
484 if cmap is None:
485 cmap = "tab10"
--> 486 return plot_ensemble_perfect_model(
487 self, variable=variable, ax=ax, show_members=show_members, cmap=cmap
488 )
File ~/checkouts/readthedocs.org/user_builds/climpred/checkouts/latest/climpred/graphics.py:354, in plot_ensemble_perfect_model(pm, variable, ax, show_members, cmap, x)
351 _cmap = mpl.cm.get_cmap(cmap, initialized.init.size)
353 for ii, i in enumerate(initialized.init.values):
--> 354 dsi = initialized.sel(init=i)
355 if uninitialized_present:
356 dsu = uninitialized.sel(init=i)
File ~/checkouts/readthedocs.org/user_builds/climpred/conda/latest/lib/python3.9/site-packages/xarray/core/dataarray.py:1527, in DataArray.sel(self, indexers, method, tolerance, drop, **indexers_kwargs)
1417 def sel(
1418 self: T_DataArray,
1419 indexers: Mapping[Any, Any] | None = None,
(...)
1423 **indexers_kwargs: Any,
1424 ) -> T_DataArray:
1425 """Return a new DataArray whose data is given by selecting index
1426 labels along the specified dimension(s).
1427
(...)
1525 Dimensions without coordinates: points
1526 """
-> 1527 ds = self._to_temp_dataset().sel(
1528 indexers=indexers,
1529 drop=drop,
1530 method=method,
1531 tolerance=tolerance,
1532 **indexers_kwargs,
1533 )
1534 return self._from_temp_dataset(ds)
File ~/checkouts/readthedocs.org/user_builds/climpred/conda/latest/lib/python3.9/site-packages/xarray/core/dataset.py:2565, in Dataset.sel(self, indexers, method, tolerance, drop, **indexers_kwargs)
2504 """Returns a new dataset with each array indexed by tick labels
2505 along the specified dimension(s).
2506
(...)
2562 DataArray.sel
2563 """
2564 indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "sel")
-> 2565 query_results = map_index_queries(
2566 self, indexers=indexers, method=method, tolerance=tolerance
2567 )
2569 if drop:
2570 no_scalar_variables = {}
File ~/checkouts/readthedocs.org/user_builds/climpred/conda/latest/lib/python3.9/site-packages/xarray/core/indexing.py:183, in map_index_queries(obj, indexers, method, tolerance, **indexers_kwargs)
181 results.append(IndexSelResult(labels))
182 else:
--> 183 results.append(index.sel(labels, **options))
185 merged = merge_sel_results(results)
187 # drop dimension coordinates found in dimension indexers
188 # (also drop multi-index if any)
189 # (.sel() already ensures alignment)
File ~/checkouts/readthedocs.org/user_builds/climpred/conda/latest/lib/python3.9/site-packages/xarray/core/indexes.py:473, in PandasIndex.sel(self, labels, method, tolerance)
471 else:
472 try:
--> 473 indexer = self.index.get_loc(label_value)
474 except KeyError as e:
475 raise KeyError(
476 f"not all values found in index {coord_name!r}. "
477 "Try setting the `method` keyword argument (example: method='nearest')."
478 ) from e
File ~/checkouts/readthedocs.org/user_builds/climpred/conda/latest/lib/python3.9/site-packages/xarray/coding/cftimeindex.py:468, in CFTimeIndex.get_loc(self, key, method, tolerance)
466 return self._get_string_slice(key)
467 else:
--> 468 return pd.Index.get_loc(self, key, method=method, tolerance=tolerance)
TypeError: get_loc() got an unexpected keyword argument 'method'

pm_aggregated_skill = pm_aggregated.verify(dim=dim, **verify_kwargs)
pm_aggregated_skill[v].plot(hue="skill")
/Users/aaron.spring/Coding/climpred/climpred/classes.py:1597: UserWarning: Calculate persistence from lead=1 instead of lead=0 (recommended).
warnings.warn(
[<matplotlib.lines.Line2D at 0x14d4628f0>,
<matplotlib.lines.Line2D at 0x14d4631c0>,
<matplotlib.lines.Line2D at 0x14b69cf10>]

verify
on each grid point, then aggregate#
Here, we first calculate prediction skill with PerfectModelEnsemble.verify()
and then average the SST skill over the North Atlantic region.
This approach answers the question: “What’s the prediction skill of North Altantic SST grid points averaged afterwards? What is the average skill of predicting [metric] in [domain]?”
Used in Frölicher et al. [2020].
grid_point_skill = pm.verify(dim=dim, **verify_kwargs)
grid_point_skill[v].plot(col="lead", row="skill", robust=True, yincrease=False)
/Users/aaron.spring/Coding/climpred/climpred/classes.py:1597: UserWarning: Calculate persistence from lead=1 instead of lead=0 (recommended).
warnings.warn(
<xarray.plot.facetgrid.FacetGrid at 0x14d8f1180>

grid_point_skill_aggregated = grid_point_skill.mean(spatial_dims)
grid_point_skill_aggregated[v].plot(hue="skill")
plt.show()

verify
over each grid point, member
and init
ialization directly#
Here, we directly calculate prediction skill with PerfectModelEnsemble.verify()
over spatial_dims
, member
and init
.
This approach answers the question: What’s the prediction skill of the North Atlantic SST pattern (calculated over all member
s and init
s)? What is the spatial skill of predicting [metric] over [domain]?
Used in Becker et al. [2014], Fransner et al. [2020].
skill_all_at_once = pm.verify(dim=dim + spatial_dims, **verify_kwargs)
skill_all_at_once[v].plot(hue="skill")
plt.show()
/Users/aaron.spring/Coding/climpred/climpred/classes.py:1597: UserWarning: Calculate persistence from lead=1 instead of lead=0 (recommended).
warnings.warn(

This approach yields very similar results to first calculating prediction skill over the spatial_dims
, i.e. calculating a pattern correlation and then doing a mean over ensemble dimensions member
and init
.
skill_spatial_then_ensemble_mean = pm.verify(dim=spatial_dims, **verify_kwargs).mean(
dim
)
skill_spatial_then_ensemble_mean[v].plot(hue="skill")
plt.show()
/Users/aaron.spring/Coding/climpred/climpred/classes.py:1597: UserWarning: Calculate persistence from lead=1 instead of lead=0 (recommended).
warnings.warn(

Summary#
The three approaches yield very different results based on the dim
keyword. Make you choose dim
according to the question you want to answer. And always compare your initialized
skill to the reference
skills, which also heavily depend on dim
. The small details matter a lot!
skills = xr.concat(
[pm_aggregated_skill, grid_point_skill_aggregated, skill_all_at_once], "method"
).assign_coords(method=["skill of aggregated", "grid skill aggregated", "all directly"])
skills[v].plot(hue="method", col="skill")
<xarray.plot.facetgrid.FacetGrid at 0x14d7f4460>

skills[v].plot(col="method", hue="skill")
<xarray.plot.facetgrid.FacetGrid at 0x14df5b970>

References#
Emily Becker, Huug van den Dool, and Qin Zhang. Predictability and Forecast Skill in NMME. Journal of Climate, 27(15):5891–5906, August 2014. doi:10.1175/JCLI-D-13-00597.1.
Filippa Fransner, François Counillon, Ingo Bethke, Jerry Tjiputra, Annette Samuelsen, Aleksi Nummelin, and Are Olsen. Ocean Biogeochemical Predictions—Initialization and Limits of Predictability. Frontiers in Marine Science, 2020. doi:10/gg22rr.
Thomas L. Frölicher, Luca Ramseyer, Christoph C. Raible, Keith B. Rodgers, and John Dunne. Potential predictability of marine ecosystem drivers. Biogeosciences, 17(7):2061–2083, April 2020. doi:10/ggr9mb.
L. Goddard, A. Kumar, A. Solomon, D. Smith, G. Boer, P. Gonzalez, V. Kharin, W. Merryfield, C. Deser, S. J. Mason, B. P. Kirtman, R. Msadek, R. Sutton, E. Hawkins, T. Fricker, G. Hegerl, C. a. T. Ferro, D. B. Stephenson, G. A. Meehl, T. Stockdale, R. Burgman, A. M. Greene, Y. Kushnir, M. Newman, J. Carton, I. Fukumori, and T. Delworth. A verification framework for interannual-to-decadal predictions experiments. Climate Dynamics, 40(1-2):245–272, January 2013. doi:10/f4jjvf.
Holger Pohlmann, Michael Botzet, Mojib Latif, Andreas Roesch, Martin Wild, and Peter Tschuck. Estimating the Decadal Predictability of a Coupled AOGCM. Journal of Climate, 17(22):4463–4472, November 2004. doi:10/d2qf62.
Aaron Spring and Tatiana Ilyina. Predictability Horizons in the Global Carbon Cycle Inferred From a Perfect-Model Framework. Geophysical Research Letters, 47(9):e2019GL085311, 2020. doi:10/ggtbv2.
Roland Séférian, Sarah Berthet, and Matthieu Chevallier. Assessing the Decadal Predictability of Land and Ocean Carbon Uptake. Geophysical Research Letters, March 2018. doi:10/gdb424.