climpred: analysis of ensemble forecast models for climate prediction

Version 1 Release

We currently only support annual forecasts, but our focus is to support sub-annual (e.g., seasonal, monthly, weekly, daily) in our next major release (v2.0.0). We provide a host of deterministic metrics, as well as some probabilistic metrics, although the latter have not been tested rigorously. We support both perfect-model and hindcast prediction ensembles, and provide PerfectModelEnsemble and HindcastEnsemble classes to make analysis easier.

See quick start and our examples to get started.

Installation

You can install the latest release of climpred using pip or conda:

pip install climpred

conda install -c conda-forge climpred

You can also install the bleeding edge (pre-release versions) by cloning this repository and running pip install . --upgrade in the main directory

Getting Started

• Overview: Why climpred?
• Scope of climpred
• Quick Start
• Examples

Overview: Why climpred?

There are many packages out there related to computing metrics on initialized geoscience predictions. However, we didn’t find any one package that unified all our needs.

Output from decadal climate prediction experiments is difficult to work with. A typical output file could contain the dimensions initialization, lead time, ensemble member, latitude, longitude, depth. climpred leverages the labeled dimensions of xarray to handle the headache of bookkeeping for you. We offer HindcastEnsemble and PerfectModelEnsemble objects that carry references (e.g., control runs, reconstructions, uninitialized ensembles) along with your decadal prediction output.

When computing lead-dependent skill scores, climpred handles all of the lag-correlating for you. We offer a suite of vectorized deterministic and probabilistic metrics that can be applied to time series and grids. It’s as easy as adding your decadal prediction output to an object and running compute: HindcastEnsemble.compute_metric(metric='rmse').

Scope of climpred

climpred aims to be the primary package used to analyze output from initialized dynamical forecast models, ranging from short-term weather forecasts to decadal climate forecasts. The code base will be driven entirely by the geoscientific prediction community through open source development. It leverages xarray to keep track of core prediction ensemble dimensions (e.g., ensemble member, initialization date, and lead time) and dask to perform out-of-memory computations on large datasets.

The primary goal of climpred is to offer a comprehensive set of analysis tools for assessing the forecasts relative to references (e.g., observations, reanalysis products, control runs, baseline forecasts). This will range from simple deterministic and probabilistic verification metrics—such as mean absolute error and various skill scores—to more advanced analysis methods, such as relative entropy and mutual information. climpred expects users to handle their domain-specific post-processing of model output, so that the package can focus on the actual analysis of forecasts.

Finally, the climpred documentation will serve as a repository of unified analysis methods through jupyter notebook examples, and will also collect relevant references and literature.

Quick Start

The easiest way to get up and running is to load in one of our example datasets (or load in some data of your own) and to convert them to either a HindcastEnsemble or PerfectModelEnsemble object.

climpred provides example datasets from the MPI-ESM-LR decadal prediction ensemble and the CESM decadal prediction ensemble. See our examples to see some analysis cases.

%matplotlib inline
import matplotlib.pyplot as plt
import xarray as xr

from climpred import HindcastEnsemble
import climpred

You can view the datasets available to be loaded with the load_datasets() command without passing any arguments:

climpred.tutorial.load_dataset()

'MPI-control-1D': area averages for the MPI control run of SST/SSS.
'MPI-control-3D': lat/lon/time for the MPI control run of SST/SSS.
'MPI-PM-DP-1D': perfect model decadal prediction ensemble area averages of SST/SSS/AMO.
'MPI-PM-DP-3D': perfect model decadal prediction ensemble lat/lon/time of SST/SSS/AMO.
'CESM-DP-SST': hindcast decadal prediction ensemble of global mean SSTs.
'CESM-DP-SSS': hindcast decadal prediction ensemble of global mean SSS.
'CESM-DP-SST-3D': hindcast decadal prediction ensemble of eastern Pacific SSTs.
'CESM-LE': uninitialized ensemble of global mean SSTs.
'MPIESM_miklip_baseline1-hind-SST-global': hindcast initialized ensemble of global mean SSTs
'MPIESM_miklip_baseline1-hist-SST-global': uninitialized ensemble of global mean SSTs
'MPIESM_miklip_baseline1-assim-SST-global': assimilation in MPI-ESM of global mean SSTs
'ERSST': observations of global mean SSTs.
'FOSI-SST': reconstruction of global mean SSTs.
'FOSI-SSS': reconstruction of global mean SSS.
'FOSI-SST-3D': reconstruction of eastern Pacific SSTs

From here, loading a dataset is easy. Note that you need to be connected to the internet for this to work – the datasets are being pulled from the climpred-data repository. Once loaded, it is cached on your computer so you can reload extremely quickly. These datasets are very small (< 1MB each) so they won’t take up much space.

hind = climpred.tutorial.load_dataset('CESM-DP-SST')
obs = climpred.tutorial.load_dataset('ERSST')

Make sure your prediction ensemble’s dimension labeling conforms to climpred’s standards. In other words, you need an init, lead, and (optional) member dimension. Make sure that your init and lead dimensions align. E.g., a November 1st, 1954 initialization should be labeled as init=1954 so that the lead=1 forecast is 1955.

print(hind)

<xarray.Dataset>
Dimensions:  (init: 64, lead: 10, member: 10)
Coordinates:
  * lead     (lead) int32 1 2 3 4 5 6 7 8 9 10
  * member   (member) int32 1 2 3 4 5 6 7 8 9 10
  * init     (init) float32 1954.0 1955.0 1956.0 1957.0 ... 2015.0 2016.0 2017.0
Data variables:
    SST      (init, lead, member) float64 ...

We’ll quickly process the data to create anomalies. CESM-DPLE’s drift-correction occurs over 1964-2014, so we’ll remove that from the observations.

# subtract climatology
obs = obs - obs.sel(time=slice(1964, 2014)).mean()

We can now create a HindcastEnsemble object and add our reference and name it 'Obs'.

hindcast = HindcastEnsemble(hind)
hindcast = hindcast.add_reference(obs, 'Obs')
print(hindcast)

<climpred.HindcastEnsemble>
Initialized Ensemble:
    SST      (init, lead, member) float64 ...
Obs:
    SST      (time) float32 -0.40146065 -0.35238647 ... 0.34601402 0.45021248
Uninitialized:
    None

We’ll remove a linear trend so that it doesn’t artificially boost our predictability. Note that climpred objects (HindcastEnsemble and PerfectModelEnsemble) can have any arbitrary xarray function applied to them. Here, we use the xarray .apply() function to apply our climpred trend removal function.

# Apply the `rm_trend` function twice to detrend our obs over time and
# detrend our initialized forecasts over init. The objects ignore an xarray
# operation if the dimension doesn't exist for the given dataset.
hindcast = hindcast.apply(climpred.stats.rm_trend, dim='time')
hindcast = hindcast.apply(climpred.stats.rm_trend, dim='init')
print(hindcast)

<climpred.HindcastEnsemble>
Initialized Ensemble:
    SST      (init, lead, member) float64 0.005165 0.03014 ... 0.1842 0.1812
Obs:
    SST      (time) float32 -0.061960407 -0.023283795 ... 0.072058104 0.165859
Uninitialized:
    None

Now we’ll quickly calculate skill and persistence. We have a variety of possible metrics to use.

init = hindcast.compute_metric(metric='acc')
persistence = hindcast.compute_persistence(metric='acc')
print(init)

<xarray.Dataset>
Dimensions:  (lead: 10)
Coordinates:
  * lead     (lead) int64 1 2 3 4 5 6 7 8 9 10
Data variables:
    SST      (lead) float64 0.6778 0.5476 0.4527 ... 0.1393 -0.03366 -0.1084
Attributes:
    prediction_skill:              calculated by climpred https://climpred.re...
    skill_calculated_by_function:  compute_hindcast
    number_of_initializations:     64
    number_of_members:             10
    metric:                        pearson_r
    comparison:                    e2r
    units:                         None
    created:                       2019-12-17 21:08:21

plt.style.use('fivethirtyeight')
f, ax = plt.subplots(figsize=(8, 3))
init.SST.plot(marker='o', markersize=10, label='skill')
persistence.SST.plot(marker='o', markersize=10, label='persistence',
                     color='#a9a9a9')
plt.legend()
ax.set(title='Global Mean SST Predictability',
       ylabel='Anomaly \n Correlation Coefficient',
       xlabel='Lead Year')
plt.show()

We can also check error in our forecasts.

init = hindcast.compute_metric(metric='rmse')
persistence = hindcast.compute_persistence(metric='rmse')

plt.style.use('fivethirtyeight')
f, ax = plt.subplots(figsize=(8, 3))
init.SST.plot(marker='o', markersize=10, label='initialized forecast')
persistence.SST.plot(marker='o', markersize=10, label='persistence',
                     color='#a9a9a9')
plt.legend()
ax.set(title='Global Mean SST Forecast Error',
       ylabel='RMSE',
       xlabel='Lead Year')
plt.show()

Examples

Demo of Perfect Model Predictability Functions

This demo demonstrates climpred’s capabilities for a perfect-model framework ensemble simulation.

What’s a perfect-model framework simulation?

A perfect-model framework uses a set of ensemble simulations that are based on a General Circulation Model (GCM) or Earth System Model (ESM) alone. There is no use of any reanalysis, reconstruction, or data product to initialize the decadal prediction ensemble. An arbitrary number of members are initialized from perturbed initial conditions (the “ensemble”), and the control simulation can be viewed as just another member.

How to compare predictability skill score: As no observational data interferes with the random climate evolution of the model, we cannot use an observation-based reference for computing skill scores. Therefore, we can compare the members with one another (m2m), against the ensemble mean (m2e), or against the control (m2c). We can also compare the ensemble mean to the control (e2c). See the comparisons page for more information.

When to use perfect-model frameworks:

• You don’t have a sufficiently long observational record to use as a reference.
• You want to avoid biases between model climatology and reanalysis climatology.
• You want to avoid sensitive reactions of biogeochemical cycles to disruptive changes in ocean physics due to assimilation.
• You want to delve into process understanding of predictability in a model without outside artifacts.

import warnings

import cartopy.crs as ccrs
import cartopy.feature as cfeature
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

import climpred

warnings.filterwarnings("ignore")

Load sample data

Here we use a subset of ensembles and members from the MPI-ESM-LR (CMIP6 version) esmControl simulation of an early state. This corresponds to vga0214 from year 3000 to 3300.

1-dimensional output

Our 1D sample output contains datasets of time series of certain spatially averaged area (‘global’, ‘North_Atlantic’) and temporally averaged period (‘ym’, ‘DJF’, …) for some lead years (1, …, 20).

ds: The ensemble dataset of all members (1, …, 10), inits (initialization years: 3014, 3023, …, 3257), areas, periods, and lead years.

control: The control dataset with the same areas and periods, as well as the years 3000 to 3299.

ds = climpred.tutorial.load_dataset('MPI-PM-DP-1D')
control = climpred.tutorial.load_dataset('MPI-control-1D')

# Add to climpred PerfectModelEnsemble object.
pm = climpred.PerfectModelEnsemble(ds)
pm = pm.add_control(control)
print(pm)

<climpred.PerfectModelEnsemble>
Initialized Ensemble:
    tos      (period, lead, area, init, member) float32 ...
    sos      (period, lead, area, init, member) float32 ...
    AMO      (period, lead, area, init, member) float32 ...
Control:
    tos      (period, time, area) float32 ...
    sos      (period, time, area) float32 ...
    AMO      (period, time, area) float32 ...
Uninitialized:
    None

We’ll sub-select annual means (‘ym’) of sea surface temperature (‘tos’) in the North Atlantic.

# Currently cannot sub-select variables. Easiest way is to just use drop, or if there's lots
# of variables, select them before creating the object.
pm = pm.sel(area='North_Atlantic', period='ym').drop(['sos', 'AMO']).reset_coords(drop=True)

Bootstrapping with Replacement

Here, we bootstrap the ensemble with replacement [Goddard et al. 2013] to compare the initialized ensemble to an “uninitialized” counterpart and a persistence forecast. The visualization is based on those used in [Li et al. 2016]. The p-value demonstrates the probability that the uninitialized or persistence beats the initialized forecast based on N=100 bootstrapping with replacement.

for metric in ['acc', 'rmse']:
    bootstrapped = pm.bootstrap(metric=metric, comparison='m2e', bootstrap=100, sig=95)
    # Hacky fix that needs to be dealt with in a PR.
    # climpred objects return a dataset. graphics module wants a DataArray but looks
    # for the attributes that are attached to the Dataset.
    bs = bootstrapped['tos']
    bs.attrs = bootstrapped.attrs
    climpred.graphics.plot_bootstrapped_skill_over_leadyear(bs, sig=95)
    plt.title(' '.join(['SST', 'North Atlantic', 'Annual:', metric]),fontsize=18)
    plt.ylabel(metric)
    plt.show()

bootstrapping iteration: 100%|██████████| 100/100 [00:39<00:00,  2.55it/s]

bootstrapping iteration: 100%|██████████| 100/100 [00:37<00:00,  2.69it/s]

Computing Skill with Different Comparison Methods

Here, we use compute_perfect_model to compute the Anomaly Correlation Coefficient (ACC) with different comparison methods. This generates different ACC values by design. See the comparisons page for a description of the various ways to compute skill scores for a perfect-model framework.

for c in ['e2c','m2c','m2e','m2m']:
    pm.compute_metric(metric='acc', comparison=c)['tos'].plot(label=c)
# Persistence computation for a baseline.
pm.compute_persistence(metric='acc')['tos'].plot(label='persistence', ls=':')
plt.ylabel('ACC')
plt.xticks(np.arange(1,21))
plt.legend()
plt.title('Different forecast-reference comparisons for pearson_r \n lead to systematically different magnitude of skill score')
plt.show()

3-dimensional output (maps)

We also have some sample output that contains gridded time series on the curvilinear MPI grid. Our compute functions (compute_perfect_model, compute_persistence) are indifferent to any dimensions that exist in addition to init, member, and lead. In other words, the functions are set up to make these computations on a grid, if one includes lat, lon, lev, depth, etc.

ds3d: The ensemble dataset of members (1, 2, 3, 4), inits (initialization years: 3014, 3061, 3175, 3237), and lead years (1, 2, 3, 4, 5).

control3d: The control dataset spanning (3000, …, 3049).

Note: These are very small subsets of the actual MPI simulations so that we could host the sample output maps on Github.

# Sea surface temperature
ds3d = climpred.tutorial.load_dataset('MPI-PM-DP-3D') \
                        .sel(init=3014) \
                        .expand_dims('init')['tos']
control3d = climpred.tutorial.load_dataset('MPI-control-3D')['tos']

# Create climpred PerfectModelEnsemble object.
pm = climpred.PerfectModelEnsemble(ds3d)
pm = pm.add_control(control3d)
print(pm)

<climpred.PerfectModelEnsemble>
Initialized Ensemble:
    tos      (init, lead, member, y, x) float32 nan nan nan nan ... nan nan nan
Control:
    tos      (time, y, x) float32 ...
Uninitialized:
    None

Maps of Skill by Lead Year

pm.compute_metric(metric='rmse', comparison='m2e')['tos'].T.plot(col='lead', robust=True, yincrease=False)

<xarray.plot.facetgrid.FacetGrid at 0x11abb3cf8>

References

1. Bushuk, Mitchell, Rym Msadek, Michael Winton, Gabriel Vecchi, Xiaosong Yang, Anthony Rosati, and Rich Gudgel. “Regional Arctic Sea–Ice Prediction: Potential versus Operational Seasonal Forecast Skill.” Climate Dynamics, June 9, 2018. https://doi.org/10/gd7hfq.
2. Collins, Matthew, and Sinha Bablu. “Predictability of Decadal Variations in the Thermohaline Circulation and Climate.” Geophysical Research Letters 30, no. 6 (March 22, 2003). https://doi.org/10/cts3cr.
3. Goddard, Lisa, et al. “A verification framework for interannual-to-decadal predictions experiments.” Climate Dynamics 40.1-2 (2013): 245-272.
4. Griffies, S. M., and K. Bryan. “A Predictability Study of Simulated North Atlantic Multidecadal Variability.” Climate Dynamics 13, no. 7–8 (August 1, 1997): 459–87. https://doi.org/10/ch4kc4.
5. Hawkins, Ed, Steffen Tietsche, Jonathan J. Day, Nathanael Melia, Keith Haines, and Sarah Keeley. “Aspects of Designing and Evaluating Seasonal-to-Interannual Arctic Sea-Ice Prediction Systems.” Quarterly Journal of the Royal Meteorological Society 142, no. 695 (January 1, 2016): 672–83. https://doi.org/10/gfb3pn.
6. Li, Hongmei, Tatiana Ilyina, Wolfgang A. Müller, and Frank Sienz. “Decadal Predictions of the North Atlantic CO2 Uptake.” Nature Communications 7 (March 30, 2016): 11076. https://doi.org/10/f8wkrs.
7. Pohlmann, Holger, Michael Botzet, Mojib Latif, Andreas Roesch, Martin Wild, and Peter Tschuck. “Estimating the Decadal Predictability of a Coupled AOGCM.” Journal of Climate 17, no. 22 (November 1, 2004): 4463–72. https://doi.org/10/d2qf62.

Hindcast Predictions of Equatorial Pacific SSTs

In this example, we evaluate hindcasts (retrospective forecasts) of sea surface temperatures in the eastern equatorial Pacific from CESM-DPLE. These hindcasts are evaluated against a forced ocean–sea ice simulation that initializes the model.

See the quick start for an analysis of time series (rather than maps) from a hindcast prediction ensemble.

import warnings

import cartopy.crs as ccrs
import cartopy.feature as cfeature
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

import climpred
from climpred import HindcastEnsemble

warnings.filterwarnings("ignore")

We’ll load in a small region of the eastern equatorial Pacific for this analysis example.

climpred.tutorial.load_dataset()

'MPI-control-1D': area averages for the MPI control run of SST/SSS.
'MPI-control-3D': lat/lon/time for the MPI control run of SST/SSS.
'MPI-PM-DP-1D': perfect model decadal prediction ensemble area averages of SST/SSS/AMO.
'MPI-PM-DP-3D': perfect model decadal prediction ensemble lat/lon/time of SST/SSS/AMO.
'CESM-DP-SST': hindcast decadal prediction ensemble of global mean SSTs.
'CESM-DP-SSS': hindcast decadal prediction ensemble of global mean SSS.
'CESM-DP-SST-3D': hindcast decadal prediction ensemble of eastern Pacific SSTs.
'CESM-LE': uninitialized ensemble of global mean SSTs.
'MPIESM_miklip_baseline1-hind-SST-global': hindcast initialized ensemble of global mean SSTs
'MPIESM_miklip_baseline1-hist-SST-global': uninitialized ensemble of global mean SSTs
'MPIESM_miklip_baseline1-assim-SST-global': assimilation in MPI-ESM of global mean SSTs
'ERSST': observations of global mean SSTs.
'FOSI-SST': reconstruction of global mean SSTs.
'FOSI-SSS': reconstruction of global mean SSS.
'FOSI-SST-3D': reconstruction of eastern Pacific SSTs

hind = climpred.tutorial.load_dataset('CESM-DP-SST-3D')['SST']
recon = climpred.tutorial.load_dataset('FOSI-SST-3D')['SST']
print(hind)

<xarray.DataArray 'SST' (init: 64, lead: 10, nlat: 37, nlon: 26)>
[615680 values with dtype=float32]
Coordinates:
    TLAT     (nlat, nlon) float64 ...
    TLONG    (nlat, nlon) float64 ...
  * init     (init) float32 1954.0 1955.0 1956.0 1957.0 ... 2015.0 2016.0 2017.0
  * lead     (lead) int32 1 2 3 4 5 6 7 8 9 10
    TAREA    (nlat, nlon) float64 ...
Dimensions without coordinates: nlat, nlon

These two example products cover a small portion of the eastern equatorial Pacific.

ax = plt.axes(projection=ccrs.Orthographic(-80, 0))
p = ax.pcolormesh(recon.TLONG, recon.TLAT, recon.mean('time'),
              transform=ccrs.PlateCarree(), cmap='twilight')
ax.add_feature(cfeature.LAND, color='#d3d3d3')
ax.set_global()
plt.colorbar(p, label='Sea Surface Temperature [degC]')
ax.set(title='Example Data Coverage')

[Text(0.5, 1.0, 'Example Data Coverage')]

We first need to remove the same climatology that was used to drift-correct the CESM-DPLE. Then we’ll create a detrended version of our two products to assess detrended predictability.

# Remove 1964-2014 climatology.
recon = recon - recon.sel(time=slice(1964, 2014)).mean('time')

Although functions can be called directly in climpred, we suggest that you use our classes (HindcastEnsemble and PerfectModelEnsemble) to make analysis code cleaner.

hindcast = HindcastEnsemble(hind)
hindcast = hindcast.add_reference(recon, 'Reconstruction')
print(hindcast)

<climpred.HindcastEnsemble>
Initialized Ensemble:
    SST      (init, lead, nlat, nlon) float32 ...
Reconstruction:
    SST      (time, nlat, nlon) float32 0.0029411316 0.0013866425 ... 1.4646168
Uninitialized:
    None

I’ll also detrend the reconstruction over its time dimension and initialized forecast ensemble over init.

# Apply the `rm_trend` function twice to detrend our obs over time and
# detrend our initialized forecasts over init. The objects ignore an xarray
# operation if the dimension doesn't exist for the given dataset.
hindcast = hindcast.apply(climpred.stats.rm_trend, dim='init')
hindcast = hindcast.apply(climpred.stats.rm_trend, dim='time')

Anomaly Correlation Coefficient of SSTs

We can now compute the ACC over all leads and all grid cells.

predictability = hindcast.compute_metric(metric='acc')
print(predictability)

<xarray.Dataset>
Dimensions:  (lead: 10, nlat: 37, nlon: 26)
Coordinates:
    TLONG    (lead, nlat, nlon) float64 250.8 251.9 253.1 ... 276.7 277.8 278.9
    TAREA    (lead, nlat, nlon) float64 3.661e+13 3.661e+13 ... 3.714e+13
  * nlat     (nlat) int64 0 1 2 3 4 5 6 7 8 9 ... 27 28 29 30 31 32 33 34 35 36
  * nlon     (nlon) int64 0 1 2 3 4 5 6 7 8 9 ... 16 17 18 19 20 21 22 23 24 25
  * lead     (lead) int64 1 2 3 4 5 6 7 8 9 10
    TLAT     (lead, nlat, nlon) float64 -9.75 -9.75 -9.75 ... -0.1336 -0.1336
Data variables:
    SST      (lead, nlat, nlon) float32 0.54588705 0.53977644 ... 0.088374
Attributes:
    prediction_skill:              calculated by climpred https://climpred.re...
    skill_calculated_by_function:  compute_hindcast
    number_of_initializations:     64
    metric:                        pearson_r
    comparison:                    e2r
    units:                         None
    created:                       2019-11-21 15:52:28

We use the pval keyword to get associated p-values for our ACCs. We can then mask our final maps based on .

significance = hindcast.compute_metric(metric='p_pval')

# Mask latitude and longitude by significance for stippling.
siglat = significance.TLAT.where(significance.SST <= 0.05)
siglon = significance.TLONG.where(significance.SST <= 0.05)

p = predictability.SST.plot.pcolormesh(x='TLONG', y='TLAT',
                                       transform=ccrs.PlateCarree(),
                                       col='lead', col_wrap=5,
                                       subplot_kws={'projection': ccrs.PlateCarree(),
                                                    'aspect': 3},
                                       cbar_kwargs={'label': 'Anomaly Correlation Coefficient'},
                                       vmin=-0.7, vmax=0.7,
                                       cmap='RdYlBu_r')
for i, ax in enumerate(p.axes.flat):
    ax.add_feature(cfeature.LAND, color='#d3d3d3', zorder=4)
    ax.gridlines(alpha=0.3, color='k', linestyle=':')
    # Add significance stippling
    ax.scatter(siglon.isel(lead=i),
               siglat.isel(lead=i),
               color='k',
               marker='.',
               s=1.5,
               transform=ccrs.PlateCarree())

Root Mean Square Error of SSTs

We can also check error in our forecasts, just by changing the metric keyword.

rmse = hindcast.compute_metric(metric='rmse')

p = rmse.SST.plot.pcolormesh(x='TLONG', y='TLAT',
                            transform=ccrs.PlateCarree(),
                            col='lead', col_wrap=5,
                            subplot_kws={'projection': ccrs.PlateCarree(),
                                        'aspect': 3},
                            cbar_kwargs={'label': 'Root Mean Square Error (degC)'},
                            cmap='Purples')
for ax in p.axes.flat:
    ax.add_feature(cfeature.LAND, color='#d3d3d3', zorder=4)
    ax.gridlines(alpha=0.3, color='k', linestyle=':')

Diagnosing Potential Predictability

This demo demonstrates climpred’s capabilities to diagnose areas containing potentially predictable variations from a control or reference alone without requiring multi-member, multi-initialization simulations. This notebook identifies the slow components of internal variability that indicate potential predictability. Here, we showcase a set of methods to show regions indicating probabilities for decadal predictability.

import warnings
%matplotlib inline
import climpred
warnings.filterwarnings("ignore")

# Sea surface temperature
varname='tos'
control3d = climpred.tutorial.load_dataset('MPI-control-3D')[varname].load()

Diagnostic Potential Predictability (DPP)

We can first use the [Resplandy 2015] and [Seferian 2018] method for computing the unbiased DPP by not chunking the time dimension.

# calculate DPP with m=10
DPP10 = climpred.stats.dpp(control3d, m=10, chunk=False)
# calculate a threshold by random shuffling (based on bootstrapping with replacement at 95% significance level)
threshold = climpred.bootstrap.dpp_threshold(control3d,
                                             m=10,
                                             chunk=False,
                                             bootstrap=10,
                                             sig=95)
# plot grid cells where DPP above threshold
DPP10.where(DPP10 > threshold).plot(yincrease=False, vmin=-0.1, vmax=0.6, cmap='viridis')

<matplotlib.collections.QuadMesh at 0x7fa3315bb908>

Now, we can turn on chunking (the default for this function) to use the [Boer 2004] method.

# chunk = True signals the Boer 2004 method
DPP10 = climpred.stats.dpp(control3d, m=10, chunk=True)
threshold = climpred.bootstrap.dpp_threshold(control3d,
                                             m=10,
                                             chunk=True,
                                             bootstrap=50,
                                             sig=95)
DPP10.where(DPP10>0).plot(yincrease=False, vmin=-0.1, vmax=0.6, cmap='viridis')

<matplotlib.collections.QuadMesh at 0x7fa331c7e0b8>

Variance-Weighted Mean Period

A periodogram is computed based on a control simulation to extract the mean period of variations, which are weighted by the respective variance. Regions with a high mean period value indicate low-frequency variations with are potentially predictable [Branstator 2010].

vwmp = climpred.stats.varweighted_mean_period(control3d, dim='time')
threshold = climpred.bootstrap.varweighted_mean_period_threshold(control3d,
                                                             bootstrap=10)
vwmp.where(vwmp > threshold).plot(yincrease=False, robust=True)

<matplotlib.collections.QuadMesh at 0x7fa3131b7240>

Lag-1 Autocorrelation

The lag-1 autocorrelation also indicates where slower modes of variability occur by identifying regions with high temporal correlation [vonStorch 1999].

# use climpred.bootstrap._bootstrap_func to wrap any stats function
threshold = climpred.bootstrap._bootstrap_func(climpred.stats.autocorr,control3d,'time',bootstrap=100)
corr_ef = climpred.stats.autocorr(control3d, dim='time')
corr_ef.where(corr_ef>threshold).plot(yincrease=False, robust=False)

<matplotlib.collections.QuadMesh at 0x7fa3213b8cc0>

Decorrelation time

Taking the lagged correlation further over all lags, the decorrelation time shows the time after which the autocorrelation fell beyond its e-folding [vonStorch 1999]

threshold = climpred.bootstrap._bootstrap_func(climpred.stats.decorrelation_time,control3d,'time',bootstrap=100)
decorr_time = climpred.stats.decorrelation_time(control3d)
decorr_time.where(decorr_time>threshold).plot(yincrease=False, robust=False)

<matplotlib.collections.QuadMesh at 0x7fa3312afa58>

Verify diagnostic potential predictability in predictability simulations

Do we find predictability in the areas highlighted above also in perfect-model experiments?

ds3d = climpred.tutorial.load_dataset('MPI-PM-DP-3D')[varname].load()

bootstrap_skill = climpred.bootstrap.bootstrap_perfect_model(ds3d,
                                          control3d,
                                          metric='rmse',
                                          comparison='m2e',
                                          bootstrap=20)

init_skill = bootstrap_skill.sel(results='skill',kind='init')
# p value: probability that random uninitialized forecasts perform better than initialized
p = bootstrap_skill.sel(results='p',kind='uninit')

init_skill.where(p<=.05).plot(col='lead', robust=True, yincrease=False)

<xarray.plot.facetgrid.FacetGrid at 0x7fa311a2c898>

The metric rmse is negatively oriented, e.g. higher values show large disprepancy between members and hence less skill.

As suggested by DPP, the variance-weighted mean period and autocorrelation, also in slight perturbed initial values ensembles there is predictability in the North Atlantic, North Pacific and Southern Ocean in sea-surface temperatures.

References

1. Boer, Georges J. “Long time-scale potential predictability in an ensemble of coupled climate models.” Climate dynamics 23.1 (2004): 29-44.
2. Resplandy, Laure, R. Séférian, and L. Bopp. “Natural variability of CO2 and O2 fluxes: What can we learn from centuries‐long climate models simulations?.” Journal of Geophysical Research: Oceans 120.1 (2015): 384-404.
3. Séférian, Roland, Sarah Berthet, and Matthieu Chevallier. “Assessing the Decadal Predictability of Land and Ocean Carbon Uptake.” Geophysical Research Letters, March 15, 2018. https://doi.org/10/gdb424.
4. Branstator, Grant, and Haiyan Teng. “Two Limits of Initial-Value Decadal Predictability in a CGCM.” Journal of Climate 23, no. 23 (August 27, 2010): 6292–6311. https://doi.org/10/bwq92h.
5. Storch, H. v, and Francis W. Zwiers. Statistical Analysis in Climate Research. Cambridge ; New York: Cambridge University Press, 1999.

Temporal and spatial smoothing

This demo demonstrates climpred’s capabilities to postprocess decadal prediction output before skill verification. Here, we showcase a set of methods to smooth out noise in the spatial and temporal domain.

import warnings
%matplotlib inline
import climpred
warnings.filterwarnings("ignore")

# Sea surface temperature
varname='tos'
ds3d = climpred.tutorial.load_dataset('MPI-PM-DP-3D')[varname]
control3d = climpred.tutorial.load_dataset('MPI-control-3D')[varname]

Temporal smoothing

In order to reduce temporal noise, decadal predictions are recommended to take multi-year averages [Goddard2013].

ds3d_ts = climpred.smoothing.temporal_smoothing(ds3d,smooth_kws={'lead':4})
control3d_ts = climpred.smoothing.temporal_smoothing(control3d, smooth_kws={'time':4})

climpred.prediction.compute_perfect_model(ds3d_ts,
                                          control3d_ts,
                                          metric='rmse',
                                          comparison='m2e') \
                   .plot(col='lead', robust=True, yincrease=False)

<xarray.plot.facetgrid.FacetGrid at 0x1088d1518>

Compare to without smoothing:

climpred.prediction.compute_perfect_model(ds3d,
                                          control3d,
                                          metric='rmse',
                                          comparison='m2e') \
                   .plot(col='lead', vmax=.69, yincrease=False)

<xarray.plot.facetgrid.FacetGrid at 0x1249b5278>

Note: When using temporal_smoothing on compute_hindcast, set rename_dim=False and after calculating the skill _reset_temporal_axis to get proper labeling of the lead dimension.

hind = climpred.tutorial.load_dataset('CESM-DP-SST-3D').load()['SST']
reconstruction = climpred.tutorial.load_dataset('FOSI-SST-3D').load()['SST']
# get anomaly reconstruction
reconstruction = reconstruction - reconstruction.mean('time')

hind_ts = climpred.smoothing.temporal_smoothing(hind,smooth_kws={'lead':4},rename_dim=False)
reconstruction_ts = climpred.smoothing.temporal_smoothing(reconstruction, smooth_kws={'time':4},rename_dim=False)

s = climpred.prediction.compute_hindcast(hind_ts,
                                          reconstruction_ts,
                                          metric='rmse',
                                          comparison='e2r')
s = climpred.smoothing._reset_temporal_axis(s,smooth_kws={'lead':4})
s.plot(col='lead', robust=True)

<xarray.plot.facetgrid.FacetGrid at 0x1206d6ef0>

Spatial smoothing

In order to reduce spatial noise, global decadal predictions are recommended to get regridded to a 5 degree longitude x 5 degree latitude grid as recommended [Goddard2013].

ds3d_ss = climpred.smoothing.spatial_smoothing_xesmf(ds3d,d_lon_lat_kws={'lon':5, 'lat':5})
control3d_ss = climpred.smoothing.spatial_smoothing_xesmf(control3d, d_lon_lat_kws={'lon':5,'lat':5})

Create weight file: bilinear_220x256_36x73.nc
Reuse existing file: bilinear_220x256_36x73.nc

climpred.prediction.compute_perfect_model(ds3d_ss,
                                          control3d_ss,
                                          metric='rmse',
                                          comparison='m2e') \
                   .plot(col='lead', robust=True, yincrease=True)

<xarray.plot.facetgrid.FacetGrid at 0x1220bf588>

Alternatively, also climpred.smoothing.spatial_smoothing_xrcoarsen aggregates gridcells like xr_coarsen.

smooth_goddard2013 creates 4-year means and 5x5 degree regridding as suggested in [Goddard2013].

climpred.smoothing.smooth_goddard_2013(ds3d).coords

Reuse existing file: bilinear_220x256_36x73.nc

Coordinates:
  * lead     (lead) <U3 '1-4' '2-5'
  * init     (init) int64 3014 3061 3175 3237
  * member   (member) int64 1 2 3 4
  * lon      (lon) float64 -180.0 -175.0 -170.0 -165.0 ... 170.0 175.0 180.0
  * lat      (lat) float64 -83.97 -78.97 -73.97 -68.97 ... 81.03 86.03 91.03

References

1. Goddard, L., A. Kumar, A. Solomon, D. Smith, G. Boer, P. Gonzalez, V. Kharin, et al. “A Verification Framework for Interannual-to-Decadal Predictions Experiments.” Climate Dynamics 40, no. 1–2 (January 1, 2013): 245–72. https://doi.org/10/f4jjvf.

User Guide

• Setting Up Your Dataset
• PredictionEnsemble Objects
• Comparisons
• Metrics
• Prediction Terminology
• Baseline Forecasts

Setting Up Your Dataset

climpred relies on a consistent naming system for xarray dimensions. This allows things to run more easily under-the-hood.

Prediction ensembles are expected at the minimum to contain dimensions init and lead. init is the initialization dimension, that relays the time steps at which the ensemble was initialized. lead is the lead time of the forecasts from initialization. Another crucial dimension is member, which holds the various ensemble members. Any additional dimensions will be passed through climpred without issue: these could be things like lat, lon, depth, etc.

Control runs, references, and observational products are expected to contain the time dimension at the minimum. For best use of climpred, their time dimension should cover the full length of init from the accompanying prediction ensemble, if possible. These products can also include additional dimensions, such as lat, lon, depth, etc.

See the below table for a summary of dimensions used in climpred, and data types that climpred supports for them.

short_name	types	long_name
lead	int	lead timestep after initialization [init]
init	int	initialization: start date of experiment
member	int, str	ensemble member

PredictionEnsemble Objects

One of the major features of climpred is our objects that are based upon the PredictionEnsemble class. We supply users with a HindcastEnsemble and PerfectModelEnsemble object. We encourage users to take advantage of these high-level objects, which wrap all of our core functions. These objects don’t comprehensively cover all functions yet, but eventually we’ll deprecate direct access to the function calls in favor of the lightweight objects.

Briefly, we consider a HindcastEnsemble to be one that is initialized from some observational-like product (e.g., assimilated data, reanalysis products, or a model reconstruction). Thus, this object is built around comparing the initialized ensemble to various observational products. In contrast, a PerfectModelEnsemble is one that is initialized off of a model control simulation. These forecasting systems are not meant to be compared directly to real-world observations. Instead, they provide a contained model environment with which to theoretically study the limits of predictability. You can read more about the terminology used in climpred here.

Let’s create a demo object to explore some of the functionality and why they are much smoother to use than direct function calls.

%matplotlib inline
import matplotlib.pyplot as plt
import xarray as xr

from climpred import HindcastEnsemble
import climpred

We can pull in some sample data that is packaged with climpred.

climpred.tutorial.load_dataset()

'MPI-control-1D': area averages for the MPI control run of SST/SSS.
'MPI-control-3D': lat/lon/time for the MPI control run of SST/SSS.
'MPI-PM-DP-1D': perfect model decadal prediction ensemble area averages of SST/SSS/AMO.
'MPI-PM-DP-3D': perfect model decadal prediction ensemble lat/lon/time of SST/SSS/AMO.
'CESM-DP-SST': hindcast decadal prediction ensemble of global mean SSTs.
'CESM-DP-SSS': hindcast decadal prediction ensemble of global mean SSS.
'CESM-DP-SST-3D': hindcast decadal prediction ensemble of eastern Pacific SSTs.
'CESM-LE': uninitialized ensemble of global mean SSTs.
'MPIESM_miklip_baseline1-hind-SST-global': hindcast initialized ensemble of global mean SSTs
'MPIESM_miklip_baseline1-hist-SST-global': uninitialized ensemble of global mean SSTs
'MPIESM_miklip_baseline1-assim-SST-global': assimilation in MPI-ESM of global mean SSTs
'ERSST': observations of global mean SSTs.
'FOSI-SST': reconstruction of global mean SSTs.
'FOSI-SSS': reconstruction of global mean SSS.
'FOSI-SST-3D': reconstruction of eastern Pacific SSTs

HindcastEnsemble

We’ll start out with a HindcastEnsemble demo, followed by a PerfectModelEnsemble case.

hind = climpred.tutorial.load_dataset('CESM-DP-SST') # CESM-DPLE hindcast ensemble output.
obs = climpred.tutorial.load_dataset('ERSST') # ERSST observations.
recon = climpred.tutorial.load_dataset('FOSI-SST') # Reconstruction simulation that initialized CESM-DPLE.

CESM-DPLE was drift-corrected prior to uploading the output, so we just need to subtract the climatology over the same period for our other products before building the object.

obs = obs - obs.sel(time=slice(1964, 2014)).mean('time')
recon = recon - recon.sel(time=slice(1964, 2014)).mean('time')

Now we instantiate the HindcastEnsemble object and append all of our products to it.

hindcast = HindcastEnsemble(hind) # Instantiate object by passing in our initialized ensemble.
print(hindcast)

<climpred.HindcastEnsemble>
Initialized Ensemble:
    SST      (init, lead, member) float64 ...
References:
    None
Uninitialized:
    None

Now we just use the add_ methods to attach other objects. See the API here. Note that we strive to make our conventions follow those of ``xarray``’s. For example, we don’t allow inplace operations. One has to run hindcast = hindcast.add_reference(...) to modify the object upon later calls rather than just hindcast.add_reference(...).

hindcast = hindcast.add_reference(recon, 'reconstruction')
hindcast = hindcast.add_reference(obs, 'ERSST')

print(hindcast)

<climpred.HindcastEnsemble>
Initialized Ensemble:
    SST      (init, lead, member) float64 ...
reconstruction:
    SST      (time) float64 -0.05064 -0.0868 -0.1396 ... 0.3023 0.3718 0.292
ERSST:
    SST      (time) float32 -0.40146065 -0.35238647 ... 0.34601402 0.45021248
Uninitialized:
    None

You can apply most standard xarray functions directly to our objects! climpred will loop through the objects and apply the function to all applicable xarray.Datasets within the object. If you reference a dimension that doesn’t exist for the given xarray.Dataset, it will ignore it. This is useful, since the initialized ensemble is expected to have dimension init, while other products have dimension time (see more here).

Let’s start by taking the ensemble mean of the initialized ensemble so our metric computations don’t have to take the extra time on that later. I’m just going to use deterministic metrics here, so we don’t need the individual ensemble members. Note that above our initialized ensemble had a member dimension, and now it is reduced.

hindcast = hindcast.mean('member')
print(hindcast)

<climpred.HindcastEnsemble>
Initialized Ensemble:
    SST      (init, lead) float64 -0.2121 -0.1637 -0.1206 ... 0.7286 0.7532
reconstruction:
    SST      (time) float64 -0.05064 -0.0868 -0.1396 ... 0.3023 0.3718 0.292
ERSST:
    SST      (time) float32 -0.40146065 -0.35238647 ... 0.34601402 0.45021248
Uninitialized:
    None

We still have a trend in all of our products, so we could also detrend them as well.

hindcast.get_reference('reconstruction').SST.plot()

[<matplotlib.lines.Line2D at 0x110386710>]

from scipy.signal import detrend

I’m going to transpose this first since my initialized ensemble has dimensions ordered (init, lead) and scipy.signal.detrend is applied over the last axis. I’d like to detrend over the init dimension rather than lead dimension.

hindcast = hindcast.transpose().apply(detrend)

And it looks like everything got detrended by a linear fit! That wasn’t too hard.

hindcast.get_reference('reconstruction').SST.plot()

[<matplotlib.lines.Line2D at 0x129a77c18>]

hindcast.get_initialized().isel(lead=0).SST.plot()

[<matplotlib.lines.Line2D at 0x129af50f0>]

Now that we’ve done our pre-processing, let’s quickly compute some metrics. Check the metrics page here for all the keywords you can use. The API is currently pretty simple for the HindcastEnsemble. You can essentially compute standard skill metrics and a reference persistence forecast.

If you just pass a metric, it’ll compute the skill metric against all references and return a dictionary with keys of the names the user entered when adding them.

hindcast.compute_metric(metric='mse')

{'reconstruction': <xarray.Dataset>
 Dimensions:  (lead: 10)
 Coordinates:
   * lead     (lead) int32 1 2 3 4 5 6 7 8 9 10
 Data variables:
     SST      (lead) float64 0.005091 0.009096 0.008964 ... 0.01103 0.01261
 Attributes:
     prediction_skill:              calculated by climpred https://climpred.re...
     skill_calculated_by_function:  compute_hindcast
     number_of_initializations:     64
     metric:                        mse
     comparison:                    e2r
     created:                       2019-12-15 21:22:07,
 'ERSST': <xarray.Dataset>
 Dimensions:  (lead: 10)
 Coordinates:
   * lead     (lead) int32 1 2 3 4 5 6 7 8 9 10
 Data variables:
     SST      (lead) float64 0.003606 0.005651 0.006373 ... 0.007823 0.009009
 Attributes:
     prediction_skill:              calculated by climpred https://climpred.re...
     skill_calculated_by_function:  compute_hindcast
     number_of_initializations:     64
     metric:                        mse
     comparison:                    e2r
     created:                       2019-12-15 21:22:07}

One can also directly call individual references to compare to. Here we leverage xarray’s plotting method to compute Mean Absolute Error and the Anomaly Correlation Coefficient for both our reference products, as well as the equivalent metrics computed for persistence forecasts for each of those metrics.

import numpy as np

plt.style.use('ggplot')
plt.style.use('seaborn-talk')

RECON_COLOR = '#1b9e77'
OBS_COLOR = '#7570b3'

f, axs = plt.subplots(nrows=2, figsize=(8, 8), sharex=True)

for ax, metric in zip(axs.ravel(), ['mae', 'acc']):
    handles = []
    for product, color in zip(['reconstruction', 'ERSST'], [RECON_COLOR, OBS_COLOR]):
        p1, = hindcast.compute_metric(product, metric=metric).SST.plot(ax=ax,
                                                                 marker='o',
                                                                 color=color,
                                                                 label=product,
                                                                 linewidth=2)
        p2, = hindcast.compute_persistence(product, metric=metric).SST.plot(ax=ax,
                                                                      color=color,
                                                                      linestyle='--',
                                                                      label=product + ' persistence')
        handles.append(p1)
        handles.append(p2)
    ax.set_title(metric.upper())


axs[0].set_ylabel('Mean Error [degC]')
axs[1].set_ylabel('Correlation Coefficient')
axs[0].set_xlabel('')
axs[1].set_xlabel('Lead Year')
axs[1].set_xticks(np.arange(10)+1)

# matplotlib/xarray returning weirdness for the legend handles.
handles = [i.get_label() for i in handles]

# a little trick to put the legend on the outside.
plt.legend(handles, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.0)

plt.suptitle('CESM Decadal Prediction Large Ensemble Global SSTs', fontsize=16)
plt.show()

PerfectModelEnsemble

We’ll now play around a bit with the PerfectModelEnsemble object, using sample data from the MPI perfect model configuration.

from climpred import PerfectModelEnsemble

ds = climpred.tutorial.load_dataset('MPI-PM-DP-1D') # initialized ensemble from MPI
control = climpred.tutorial.load_dataset('MPI-control-1D') # base control run that initialized it

print(ds)

<xarray.Dataset>
Dimensions:  (area: 3, init: 12, lead: 20, member: 10, period: 5)
Coordinates:
  * lead     (lead) int64 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
  * period   (period) object 'DJF' 'JJA' 'MAM' 'SON' 'ym'
  * area     (area) object 'global' 'North_Atlantic' 'North_Atlantic_SPG'
  * init     (init) int64 3014 3023 3045 3061 3124 ... 3175 3178 3228 3237 3257
  * member   (member) int64 0 1 2 3 4 5 6 7 8 9
Data variables:
    tos      (period, lead, area, init, member) float32 ...
    sos      (period, lead, area, init, member) float32 ...
    AMO      (period, lead, area, init, member) float32 ...

pm = climpred.PerfectModelEnsemble(ds)
pm = pm.add_control(control)
print(pm)

<climpred.PerfectModelEnsemble>
Initialized Ensemble:
    tos      (period, lead, area, init, member) float32 ...
    sos      (period, lead, area, init, member) float32 ...
    AMO      (period, lead, area, init, member) float32 ...
Control:
    tos      (period, time, area) float32 ...
    sos      (period, time, area) float32 ...
    AMO      (period, time, area) float32 ...
Uninitialized:
    None

Our objects are carrying sea surface temperature (tos), sea surface salinity (sos), and the Atlantic Multidecadal Oscillation index (AMO). Say we just want to look at skill metrics for temperature and salinity over the North Atlantic in JJA. We can just call a few easy xarray commands to filter down our object.

pm = pm.drop('AMO').sel(area='North_Atlantic', period='JJA')

Now we can easily compute for a host of metrics. Here I just show a number of deterministic skill metrics comparing all individual members to the initialized ensemble mean. See comparisons for more information on the comparison keyword.

METRICS = ['mse', 'rmse', 'mae', 'acc',
           'nmse', 'nrmse', 'nmae', 'msss']

result = []
for metric in METRICS:
    result.append(pm.compute_metric(metric, comparison='m2e'))

result = xr.concat(result, 'metric')
result['metric'] = METRICS

# Leverage the `xarray` plotting wrapper to plot all results at once.
result.to_array().plot(col='metric',
                       hue='variable',
                       col_wrap=4,
                       sharey=False,
                       sharex=True)

<xarray.plot.facetgrid.FacetGrid at 0x124b53e10>

It is useful to compare the initialized ensemble to an uninitialized run. See terminology for a description on “uninitialized” simulations. This gives us information about how initializations lead to enhanced predictability over knowledge of external forcing, whereas a comparison to persistence just tells us how well a dynamical forecast simulation does in comparison to a naive method. We can use the generate_uninitialized() method to bootstrap the control run and create a pseudo-ensemble that approximates what an uninitialized ensemble would look like.

pm = pm.generate_uninitialized()
print(pm)

<climpred.PerfectModelEnsemble>
Initialized Ensemble:
    tos      (lead, init, member) float32 ...
    sos      (lead, init, member) float32 ...
Control:
    tos      (time) float32 ...
    sos      (time) float32 ...
Uninitialized:
    tos      (init, member, lead) float32 13.856491 13.286671 ... 13.289608
    sos      (init, member, lead) float32 33.210403 33.172962 ... 33.176933

pm = pm.drop('tos') # Just assess for salinity.

Here I plot the ACC for the initialized, uninitialized, and persistence forecasts for North Atlantic sea surface salinity in JJA. I add circles to the lines if the correlations are statistically significant for .

# ACC for initialized ensemble
acc = pm.compute_metric('acc')
acc.sos.plot(color='red')
acc.where(pm.compute_metric('p_pval') <= 0.05).sos.plot(marker='o', linestyle='None', color='red', label='initialized')

# ACC for 'uninitialized' ensemble
acc = pm.compute_uninitialized('acc')
acc.sos.plot(color='gray')
acc.where(pm.compute_uninitialized('p_pval') <= 0.05).sos.plot(marker='o', linestyle='None', color='gray', label='uninitialized')

# ACC for persistence forecast
acc = pm.compute_persistence('acc')
acc.sos.plot(color='k', linestyle='--')
acc.where(pm.compute_persistence('p_pval') <= 0.05).sos.plot(marker='o', linestyle='None', color='k', label='persistence')

plt.legend()

<matplotlib.legend.Legend at 0x124cee588>

Metrics

All high-level functions have an optional metric argument that can be called to determine which metric is used in computing predictability (potential predictability or prediction skill).

Note

We use the phrase ‘observations’ o here to refer to the ‘truth’ data to which we compare the forecast f. These metrics can also be applied in reference to a control simulation, reconstruction, observations, etc. This would just change the resulting score from referencing skill to referencing potential predictability.

Internally, all metric functions require forecast and reference as inputs. The dimension dim is set by compute_hindcast() or compute_perfect_model() to specify over which dimensions the metric is applied. See Comparisons.

Deterministic

Deterministic metrics quantify the level to which the forecast predicts the observations. These metrics are just a special case of probabilistic metrics where a value of 100% is assigned to the forecasted value [Jolliffe2011].

Core Metrics

Pearson Anomaly Correlation Coefficient (ACC)

keyword: 'pearson_r', 'pr', 'pacc', 'acc'

A measure of the linear association between the forecast and observations that is independent of the mean and variance of the individual distributions [Jolliffe2011]. climpred uses the Pearson correlation coefficient.

climpred.metrics._pearson_r(forecast, reference, dim=None, **metric_kwargs)

Pearson’s Anomaly Correlation Coefficient (ACC).

Note

Use metric pearson_r_p_value to get the corresponding pvalue.

Parameters:

• forecast (xarray object) – forecast
• reference (xarray object) – reference
• dim (str) – dimension(s) to perform metric over. Automatically set by compute_.
• metric_kwargs – (optional) weights, skipna, see xskillscore.pearson_r

Range:

• perfect: 1
• min: -1

See also

• xskillscore.pearson_r
• xskillscore.pearson_r_p_value

Spearman Anomaly Correlation Coefficient (SACC)

keyword: 'spearman_r', 'sacc', 'sr'

A measure of how well the relationship between two variables can be described using a monotonic function.

climpred.metrics._spearman_r(forecast, reference, dim=None, **metric_kwargs)

Spearman’s Anomaly Correlation Coefficient (SACC).

Note

Use metric spearman_r_p_value to get the corresponding pvalue.

Parameters:

• forecast (xarray object) – forecast
• reference (xarray object) – reference
• dim (str) – dimension(s) to perform metric over. Automatically set by compute_.
• metric_kwargs – (optional) weights, skipna, see xskillscore.spearman_r

Range:

• perfect: 1
• min: -1

See also

• xskillscore.spearman_r
• xskillscore.spearman_r_p_value

Mean Squared Error (MSE)

keyword: 'mse'

The average of the squared difference between forecasts and observations. This incorporates both the variance and bias of the estimator.

climpred.metrics._mse(forecast, reference, dim=None, **metric_kwargs)

Mean Sqaure Error (MSE).

Parameters:

• forecast (xarray object) – forecast
• reference (xarray object) – reference
• dim (str) – dimension(s) to perform metric over. Automatically set by compute_.
• metric_kwargs – (optional) weights, skipna, see xskillscore.mse

Range:

• perfect: 0
• min: 0
• max: ∞

See also

• xskillscore.mse

Root Mean Square Error (RMSE)

keyword: 'rmse'

The square root of the average of the squared differences between forecasts and observations [Jolliffe2011]. It puts a greater influence on large errors than small errors, which makes this a good choice if large errors are undesirable or one wants to be a more conservative forecaster.

climpred.metrics._rmse(forecast, reference, dim=None, **metric_kwargs)

Root Mean Sqaure Error (RMSE).

Parameters:

• forecast (xarray object) – forecast
• reference (xarray object) – reference
• dim (str) – dimension(s) to perform metric over. Automatically set by compute_.
• metric_kwargs – (optional) weights, skipna, see xskillscore.rmse

Range:

• perfect: 0
• min: 0
• max: ∞

See also

• xskillscore.rmse

Mean Absolute Error (MAE)

keyword: 'mae'

The average of the absolute differences between forecasts and observations [Jolliffe2011]. A more robust measure of forecast accuracy than root mean square error or mean square error which is sensitive to large outlier forecast errors [EOS].

climpred.metrics._mae(forecast, reference, dim=None, **metric_kwargs)

Mean Absolute Error (MAE).

Parameters:

• forecast (xarray object) – forecast
• reference (xarray object) – reference
• dim (str) – dimension(s) to perform metric over. Automatically set by compute_.
• metric_kwargs – (optional) weights, skipna, see xskillscore.mae

Range:

• perfect: 0
• min: 0
• max: ∞

See also

• xskillscore.mae

Median Absolute Deviation (MAD)

keyword: 'mad'

The median of the absolute differences between forecasts and observations.

climpred.metrics._mad(forecast, reference, dim=None, **metric_kwargs)

Median Absolute Deviation (MAD).

Parameters:

• forecast (xarray object) – forecast
• reference (xarray object) – reference
• dim (str) – dimension(s) to perform metric over. Automatically set by compute_.
• metric_kwargs – (optional) weights, skipna, see xskillscore.mad

Range:

• perfect: 0
• min: 0
• max: ∞

See also

• xskillscore.mad

Derived Metrics

Distance-based metrics like mse can be normalized to 1. The normalization factor depends on the comparison type choosen, eg. the distance between an ensemble member and the ensemble mean is half the distance of an ensemble member with other ensemble members. (see climpred.metrics._get_norm_factor()).

Normalized Mean Square Error (NMSE)

keyword: 'nmse','nev'

climpred.metrics._nmse(forecast, reference, dim=None, **metric_kwargs)

Normalized MSE (NMSE) also known as Normalized Ensemble Variance (NEV).

Parameters:

• forecast (*) –
• reference (*) –
• dim (*) – dimension to apply metric to
• comparison (*) – name comparison needed for normalization factor fac, see climpred.metrics._get_norm_factor() (internally required to be added via **metric_kwargs)
• metric_kwargs (*) – (optional) weights, skipna, see xskillscore.mse

Range:

• 0: perfect forecast: 0
• 0 - 1: better than climatology forecast
• > 1: worse than climatology forecast

References

• Griffies, S. M., and K. Bryan. “A Predictability Study of Simulated North Atlantic Multidecadal Variability.” Climate Dynamics 13, no. 7–8 (August 1, 1997): 459–87. https://doi.org/10/ch4kc4.

Normalized Mean Absolute Error (NMAE)

keyword: 'nmae'

climpred.metrics._nmae(forecast, reference, dim=None, **metric_kwargs)

Normalized Ensemble Mean Absolute Error metric.

Parameters:

• forecast (*) –
• reference (*) –
• dim (*) – dimension to apply metric to
• comparison (*) – name comparison needed for normalization factor fac, see climpred.metrics._get_norm_factor() (internally required to be added via **metric_kwargs)
• metric_kwargs (*) – (optional) weights, skipna, see xskillscore.mae

Range:

• 0: perfect forecast: 0
• 0 - 1: better than climatology forecast
• > 1: worse than climatology forecast

References

• Griffies, S. M., and K. Bryan. “A Predictability Study of Simulated North Atlantic Multidecadal Variability.” Climate Dynamics 13, no. 7–8 (August 1, 1997): 459–87. https://doi.org/10/ch4kc4.

Normalized Root Mean Square Error (NRMSE)

keyword: 'nrmse'

climpred.metrics._nrmse(forecast, reference, dim=None, **metric_kwargs)

Normalized Root Mean Square Error (NRMSE) metric.

Parameters:

• forecast (*) –
• reference (*) –
• dim (*) – dimension to apply metric to
• comparison (*) – name comparison needed for normalization factor fac, see climpred.metrics._get_norm_factor() (internally required to be added via **metric_kwargs)
• metric_kwargs (*) – (optional) weights, skipna, see xskillscore.rmse

Range:

• 0: perfect forecast
• 0 - 1: better than climatology forecast
• > 1: worse than climatology forecast

References

• Bushuk, Mitchell, Rym Msadek, Michael Winton, Gabriel Vecchi, Xiaosong Yang, Anthony Rosati, and Rich Gudgel. “Regional Arctic Sea–Ice Prediction: Potential versus Operational Seasonal Forecast Skill.” Climate Dynamics, June 9, 2018. https://doi.org/10/gd7hfq.
• Hawkins, Ed, Steffen Tietsche, Jonathan J. Day, Nathanael Melia, Keith Haines, and Sarah Keeley. “Aspects of Designing and Evaluating Seasonal-to-Interannual Arctic Sea-Ice Prediction Systems.” Quarterly Journal of the Royal Meteorological Society 142, no. 695 (January 1, 2016): 672–83. https://doi.org/10/gfb3pn.

Mean Square Skill Score (MSSS)

keyword: 'msss','ppp'

climpred.metrics._ppp(forecast, reference, dim=None, **metric_kwargs)

Prognostic Potential Predictability (PPP) metric.

Parameters:

• forecast (xarray object) – forecast
• reference (xarray object) – reference
• dim (str) – dimension(s) to perform metric over. Automatically set by compute_.
• comparison (str) – name comparison needed for normalization factor fac, climpred.metrics._get_norm_factor() (internally required to be added via **metric_kwargs)
• metric_kwargs – (optional) weights, skipna, see xskillscore.mse

Range:

• 1: perfect forecast
• positive: better than climatology forecast
• negative: worse than climatology forecast

References

• Griffies, S. M., and K. Bryan. “A Predictability Study of Simulated North Atlantic Multidecadal Variability.” Climate Dynamics 13, no. 7–8 (August 1, 1997): 459–87. https://doi.org/10/ch4kc4.
• Pohlmann, Holger, Michael Botzet, Mojib Latif, Andreas Roesch, Martin Wild, and Peter Tschuck. “Estimating the Decadal Predictability of a Coupled AOGCM.” Journal of Climate 17, no. 22 (November 1, 2004): 4463–72. https://doi.org/10/d2qf62.
• Bushuk, Mitchell, Rym Msadek, Michael Winton, Gabriel Vecchi, Xiaosong Yang, Anthony Rosati, and Rich Gudgel. “Regional Arctic Sea–Ice Prediction: Potential versus Operational Seasonal Forecast Skill. Climate Dynamics, June 9, 2018. https://doi.org/10/gd7hfq.

Mean Absolute Percentage Error (MAPE)

keyword: 'mape'

The mean of the absolute differences between forecasts and observations normalized by observations.

climpred.metrics._mape(forecast, reference, dim=None, **metric_kwargs)

Mean Absolute Percentage Error (MAPE).

rac{|f-o|}{|o|}

Args:

forecast (xarray object): forecast reference (xarray object): reference dim (str): dimension(s) to perform metric over.

Automatically set by compute_.

metric_kwargs: (optional) weights, skipna, see xskillscore.mape

Range:

• perfect: 0
• min: 0
• max: ∞

See also:

• xskillscore.mape

Symmetric Mean Absolute Percentage Error (sMAPE)

keyword: 'smape'

The mean of the absolute differences between forecasts and observations normalized by their sum.

climpred.metrics._smape(forecast, reference, dim=None, **metric_kwargs)

symmetric Mean Absolute Percentage Error (sMAPE).

rac{|f-o|}{|f|+|o|}

Args:

forecast (xarray object): forecast reference (xarray object): reference dim (str): dimension(s) to perform metric over.

Automatically set by compute_.

metric_kwargs: (optional) weights, skipna, see xskillscore.smape

Range:

• perfect: 0
• min: 0
• max: 1

See also:

• xskillscore.smape

Unbiased ACC

keyword: 'uacc'

climpred.metrics._uacc(forecast, reference, dim=None, **metric_kwargs)

Bushuk’s unbiased ACC (uACC).

Parameters:

• forecast (*) –
• reference (*) –
• dim (*) – dimension to apply metric to
• comparison (*) – name comparison needed for normalization factor fac, see climpred.metrics._get_norm_factor() (internally required to be added via **metric_kwargs)
• metric_kwargs (*) – (optional) weights, skipna, see xskillscore.mse

Range:

• 1: perfect
• 0 - 1: better than climatology

References

• Bushuk, Mitchell, Rym Msadek, Michael Winton, Gabriel Vecchi, Xiaosong Yang, Anthony Rosati, and Rich Gudgel. “Regional Arctic Sea–Ice Prediction: Potential versus Operational Seasonal Forecast Skill. Climate Dynamics, June 9, 2018. https://doi.org/10/gd7hfq.

Murphy decomposition metrics

[Murphy1988] relates the MSSS with ACC and unconditional bias.

Standard Ratio

keyword: 'std_ratio'

climpred.metrics._std_ratio(forecast, reference, dim=None, **metric_kwargs)

Ratio of standard deviations of reference over forecast.

Parameters:

• forecast (xarray object) – forecast
• reference (xarray object) – reference
• dim (str) – dimension(s) to perform metric over. Automatically set by compute_.

References

• https://www-miklip.dkrz.de/about/murcss/

Unconditional Bias

keyword: 'bias', 'unconditional_bias', 'u_b'

climpred.metrics._bias(forecast, reference, dim=None, **metric_kwargs)

Unconditional bias.

Parameters:

• forecast (xarray object) – forecast
• reference (xarray object) – reference
• dim (str) – dimension(s) to perform metric over. Automatically set by compute_.

Range:

• pos: positive bias
• neg: negative bias
• perfect: 0

References

• https://www.cawcr.gov.au/projects/verification/
• https://www-miklip.dkrz.de/about/murcss/

Bias Slope

keyword: 'bias_slope'

climpred.metrics._bias_slope(forecast, reference, dim=None, **metric_kwargs)

Bias slope between reference and forecast standard deviations.

Parameters:

• forecast (xarray object) – forecast
• reference (xarray object) – reference
• dim (str) – dimension(s) to perform metric over. Automatically set by compute_.

References

• https://www-miklip.dkrz.de/about/murcss/

Conditional Bias

keyword: 'conditional_bias', c_b'

climpred.metrics._conditional_bias(forecast, reference, dim=None, **metric_kwargs)

Conditional bias between forecast and reference.

Parameters:

• forecast (xarray object) – forecast
• reference (xarray object) – reference
• dim (str) – dimension(s) to perform metric over. Automatically set by compute_.

References

• https://www-miklip.dkrz.de/about/murcss/

Murphy’s Mean Square Skill Score

keyword: 'msss_murphy'

climpred.metrics._msss_murphy(forecast, reference, dim=None, **metric_kwargs)

Murphy’s Mean Square Skill Score (MSSS).

Parameters:

• forecast (xarray object) – forecast
• reference (xarray object) – reference
• dim (str) – dimension(s) to perform metric over. Automatically set by compute_.
• metric_kwargs – (optional) weights, skipna

References

• https://www-miklip.dkrz.de/about/murcss/
• Murphy, Allan H. “Skill Scores Based on the Mean Square Error and Their Relationships to the Correlation Coefficient.” Monthly Weather Review 116, no. 12 (December 1, 1988): 2417–24. https://doi.org/10/fc7mxd.

Probabilistic

keyword: 'crps'

climpred.metrics._crps(forecast, reference, **metric_kwargs)

Continuous Ranked Probability Score (CRPS) is the probabilistic MSE.

Parameters:

• forecast (*) – forecast with member dim
• reference (*) – references without member dim
• metric_kwargs (*) – weights, see properscoring.crps_ensemble

Range:

• perfect: 0
• min: 0
• max: ∞

References

• Matheson, James E., and Robert L. Winkler. “Scoring Rules for Continuous Probability Distributions.” Management Science 22, no. 10 (June 1, 1976): 1087–96. https://doi.org/10/cwwt4g.

See also

• properscoring.crps_ensemble
• xskillscore.crps_ensemble

keyword: 'crpss'

climpred.metrics._crpss(forecast, reference, **metric_kwargs)

Continuous Ranked Probability Skill Score

Note

When assuming a gaussian distribution of forecasts, use default gaussian=True. If not gaussian, you may specify the distribution type, xmin/xmax/tolerance for integration (see xskillscore.crps_quadrature).

Parameters:

• forecast (*) – forecast with member dim
• reference (*) – references without member dim
• gaussian (*) – Assuming gaussian distribution for baseline skill. Default: True (optional)
• cdf_or_dist (*) – distribution to assume if not gaussian. default: scipy.stats.norm
• xmin, xmax, tol (*) – only relevant if not gaussian (see xskillscore.crps_quadrature)

Range:

• perfect: 1
• pos: better than climatology forecast
• neg: worse than climatology forecast

References

• Matheson, James E., and Robert L. Winkler. “Scoring Rules for Continuous Probability Distributions.” Management Science 22, no. 10 (June 1, 1976): 1087–96. https://doi.org/10/cwwt4g.
• Gneiting, Tilmann, and Adrian E Raftery. “Strictly Proper Scoring Rules, Prediction, and Estimation.” Journal of the American Statistical Association 102, no. 477 (March 1, 2007): 359–78. https://doi.org/10/c6758w.

Example

>>> compute_perfect_model(ds, control, metric='crpss')
>>> compute_perfect_model(ds, control, metric='crpss', gaussian=False,
                          cdf_or_dist=scipy.stats.norm, xminimum=-10,
                          xmaximum=10, tol=1e-6)

See also

• properscoring.crps_ensemble
• xskillscore.crps_ensemble

keyword: 'crpss_es'

climpred.metrics._crpss_es(forecast, reference, **metric_kwargs)

CRPSS Ensemble Spread.

Parameters:

• forecast (*) – forecast with member dim
• reference (*) – references without member dim
• metric_kwargs (*) – weights, skipna used for mse

References

• Kadow, Christopher, Sebastian Illing, Oliver Kunst, Henning W. Rust, Holger Pohlmann, Wolfgang A. Müller, and Ulrich Cubasch. “Evaluation of Forecasts by Accuracy and Spread in the MiKlip Decadal Climate Prediction System.” Meteorologische Zeitschrift, December 21, 2016, 631–43. https://doi.org/10/f9jrhw.

Range:

• perfect: 0
• else: negative

keyword: 'brier_score', 'brier', 'bs'

climpred.metrics._brier_score(forecast, reference, **metric_kwargs)

Brier score for forecasts on binary reference.

..math:

BS(f, o) = (f - o)^2

Parameters:

• forecast (*) – forecast with member dim
• reference (*) – references without member dim
• func (*) – function to be applied to reference and forecasts and then mean(‘member’) to get forecasts and reference in interval [0,1]. (required to be added via **metric_kwargs)

Reference:

• Brier, Glenn W. “VERIFICATION OF FORECASTS EXPRESSED IN TERMS OF

PROBABILITY.” Monthly Weather Review 78, no. 1 (1950). https://doi.org/10.1175/1520-0493(1950)078<0001:VOFEIT>2.0.CO;2.

Example

>>> def pos(x): return x > 0
>>> compute_perfect_model(ds, control, metric='brier_score', func=pos)

See also

• properscoring.brier_score
• xskillscore.brier_score

keyword: 'threshold_brier_score', 'tbs'

climpred.metrics._threshold_brier_score(forecast, reference, **metric_kwargs)

Brier scores of an ensemble for exceeding given thresholds. Provide threshold via metric_kwargs.

Range:

• perfect: 0
• min: 0
• max: 1

Parameters:

• forecast (*) – forecast with member dim
• reference (*) – references without member dim
• threshold (*) – Threshold to check exceedance, see properscoring.threshold_brier_score (required to be added via **metric_kwargs)

References

• Brier, Glenn W. “VERIFICATION OF FORECASTS EXPRESSED IN TERMS OF

PROBABILITY.” Monthly Weather Review 78, no. 1 (1950). https://doi.org/10.1175/1520-0493(1950)078<0001:VOFEIT>2.0.CO;2.

Example

>>> compute_perfect_model(ds, control,
                          metric='threshold_brier_score', threshold=.5)

See also

• properscoring.threshold_brier_score
• xskillscore.threshold_brier_score

User-defined metrics

You can also construct your own metrics via the climpred.metrics.Metric class.

Metric(name, function, positive, …[, …])

Master class for all metrics.

First, write your own metric function, similar to the existing ones with required arguments forecast, reference, dim=None, and **metric_kwargs:

from climpred.metrics import Metric

def _my_msle(forecast, reference, dim=None, **metric_kwargs):
    """Mean squared logarithmic error (MSLE).
    https://peltarion.com/knowledge-center/documentation/modeling-view/build-an-ai-model/loss-functions/mean-squared-logarithmic-error."""
    return ( (np.log(forecast + 1) + np.log(reference + 1) ) ** 2).mean(dim)

Then initialize this metric function with climpred.metrics.Metric:

_my_msle = Metric(
    name='my_msle',
    function=_my_msle,
    probabilistic=False,
    positive=False,
    unit_power=0,
    )

Finally, compute skill based on your own metric:

skill = compute_perfect_model(ds, control, metric='rmse', comparison=_my_msle)

Once you come up with an useful metric for your problem, consider contributing this metric to climpred, so all users can benefit from your metric, see contributing.

References

[EOS]

https://eos.org/opinions/climate-and-other-models-may-be-more-accurate-than-reported

[Jolliffe2011]

(1, 2, 3, 4) Ian T. Jolliffe and David B. Stephenson. Forecast Verification: A Practitioner’s Guide in Atmospheric Science. John Wiley & Sons, Ltd, Chichester, UK, December 2011. ISBN 978-1-119-96000-3 978-0-470-66071-3. URL: http://doi.wiley.com/10.1002/9781119960003.

[Murphy1988]

Allan H. Murphy. Skill Scores Based on the Mean Square Error and Their Relationships to the Correlation Coefficient. Monthly Weather Review, 116(12):2417–2424, December 1988. https://doi.org/10/fc7mxd.

Comparisons

Forecast skill is always evaluated against a reference for verification. In ESM-based predictions, it is common to compare the ensemble mean forecast against the reference.

In hindcast ensembles compute_hindcast(), this ensemble mean forecast (comparison='e2r') is expected to perform better than individual ensemble members (comparison='m2r') as the chaotic component of forecasts is expected to be suppressed by this averaging, while the memory of the system sustains. [Boer2016] HindcastEnsemble skill is computed by default as the ensemble mean forecast against the reference (comparison='e2r').

In perfect-model frameworks compute_perfect_model(), there are even more ways of comparisons. [Seferian2018] shows comparison of the ensemble members against the control run (comparison='m2c') and ensemble members against all other ensemble members (comparison='m2m'). Furthermore, using the ensemble mean forecast can be also verified against one control member (comparison='e2c') or all members (comparison='m2e') as done in [Griffies1997]. Perfect-model framework comparison defaults to the ensemble mean forecast verified against each member in turns (comparison='m2e').

These different comparisons demand for a normalization factor to arrive at a normalized skill of 1, when skill saturation is reached (ref: metrics).

While HindcastEnsemble skill is computed over all initializations init of the hindcast, the resulting skill is a mean forecast skill over all initializations. PerfectModelEnsemble skill is computed over a supervector comprised of all initializations and members, which allows the computation of the ACC-based skill [Bushuk2018], but also returns a mean forecast skill over all initializations. The supervector approach shown in [Bushuk2018] and just calculating a distance-based metric like rmse over the member dimension as in [Griffies1997] yield very similar results.

Compute over dimension

The optional argument dim defines over which dimension a metric is computed. We can apply a metric over dim from ['init', 'member', ['member', 'init']] in compute_perfect_model() and ['init', 'member'] in compute_hindcast(). The resulting skill is then reduced by this dim. Therefore, applying a metric over dim='member' creates a skill for all initializations individually. This can show the initial conditions dependence of skill. Likewise when computing skill over 'init', we get skill for each member. This dim argument is different from the comparison argument which just specifies how forecast and reference are defined. However, this above logic applies to deterministic metrics. Probabilistic metrics need to be applied to the member dimension and comparison from ['m2c', 'm2m'] in compute_perfect_model() and 'm2r' comparison in compute_hindcast(). Using a probabilistic metric automatically switches internally to using dim='member'.

HindcastEnsemble

keyword: 'e2r'

_e2r(ds, reference[, stack_dims])

Compare the ensemble mean forecast to a reference in HindcastEnsemble.

keyword: 'm2r'

_m2r(ds, reference[, stack_dims])

Compares each member individually to a reference in HindcastEnsemble.

PerfectModelEnsemble

keyword: 'm2e'

_m2e(ds[, stack_dims])

Compare all members to ensemble mean while leaving out the reference in

keyword: 'm2c'

_m2c(ds[, control_member, stack_dims])

Compare all other members forecasts to control member verification.

keyword: 'm2m'

_m2m(ds[, stack_dims])

Compare all members to all others in turn while leaving out the verification member.

keyword: 'e2c'

_e2c(ds[, control_member, stack_dims])

Compare ensemble mean forecast to control member verification.

User-defined comparisons

You can also construct your own comparisons via the Comparison class.

Comparison(name, function, hindcast, …[, …])

Master class for all comparisons.

First, write your own comparison function, similar to the existing ones. If a comparison should also be used for probabilistic metrics, use stack_dims to return forecast with member dimension and reference without. For deterministic metric, return forecast and reference with identical dimensions:

from climpred.comparisons import Comparison, _drop_members

def _my_m2median_comparison(ds, stack_dims=True):
    """Identical to m2e but median."""
    reference_list = []
    forecast_list = []
    supervector_dim = 'member'
    for m in ds.member.values:
        forecast = _drop_members(ds, rmd_member=[m]).median('member')
        reference = ds.sel(member=m).squeeze()
        forecast_list.append(forecast)
        reference_list.append(reference)
    reference = xr.concat(reference_list, supervector_dim)
    forecast = xr.concat(forecast_list, supervector_dim)
    forecast[supervector_dim] = np.arange(forecast[supervector_dim].size)
    reference[supervector_dim] = np.arange(reference[supervector_dim].size)
    return forecast, reference

Then initialize this comparison function with Comparison:

__my_m2median_comparison = Comparison(
    name='m2me',
    function=_my_m2median_comparison,
    probabilistic=False,
    hindcast=False)

Finally, compute skill based on your own comparison:

skill = compute_perfect_model(ds, control, metric='rmse', comparison=__my_m2median_comparison)

Once you come up with an useful comparison for your problem, consider contributing this comparison to climpred, so all users can benefit from your comparison, see contributing.

References

[Boer2016]

Boer, G. J., D. M. Smith, C. Cassou, F. Doblas-Reyes, G. Danabasoglu, B. Kirtman, Y. Kushnir, et al. “The Decadal Climate Prediction Project (DCPP) Contribution to CMIP6.” Geosci. Model Dev. 9, no. 10 (October 25, 2016): 3751–77. https://doi.org/10/f89qdf.

[Bushuk2018]

(1, 2) Mitchell Bushuk, Rym Msadek, Michael Winton, Gabriel Vecchi, Xiaosong Yang, Anthony Rosati, and Rich Gudgel. Regional Arctic sea–ice prediction: potential versus operational seasonal forecast skill. Climate Dynamics, June 2018. https://doi.org/10/gd7hfq.

[Griffies1997]

(1, 2) Griffies and K. Bryan. A predictability study of simulated North Atlantic multidecadal variability. Climate Dynamics, 13(7-8):459–487, August 1997. https://doi.org/10/ch4kc4.

[Seferian2018]

Roland Séférian, Sarah Berthet, and Matthieu Chevallier. Assessing the Decadal Predictability of Land and Ocean Carbon Uptake. Geophysical Research Letters, March 2018. https://doi.org/10/gdb424.

Prediction Terminology

Terminology is often confusing and highly variable amongst those that make predictions in the geoscience community. Here we define some common terms in climate prediction and how we use them in climpred.

Simulation Design

Initialized Ensemble

Perfect Model Experiment: m ensemble members are initialized from a control simulation at n randomly chosen initialization dates and integrated for l lead years [Griffies1997] (PerfectModelEnsemble).

Hindcast Ensemble: m ensemble members are initialized from a reference simulation (generally a reconstruction from reanalysis) at n initialization dates and integrated for l lead years [Boer2016] (HindcastEnsemble).

Uninitialized Ensemble

In this framework, an uninitialized ensemble is one that is generated by perturbing initial conditions only at one point in the historical run. These are generated via micro (round-off error perturbations) or macro (starting from completely different restart files) methods. Uninitialized ensembles are used to approximate the magnitude of internal climate variability and to confidently extract the forced response (ensemble mean) in the climate system.

In climpred, we use uninitialized ensembles as a baseline for how important (reoccurring) initializations are for lending predictability to the system. Some modeling centers (such as NCAR) provide a dynamical uninitialized ensemble (the CESM Large Ensemble) along with their initialized prediction system (the CESM Decadal Prediction Large Ensemble). If this isn’t available, one can approximate the unintiailized response by bootstrapping a control simulation.

Reconstruction:

Reconstruction/Assimilation: A “reconstruction” is a model solution that uses observations in some capacity to approximate historical conditions. This could be done via a forced simulation, such as an OMIP run that uses a dynamical ocean/sea ice core with reanalysis forcing from atmospheric winds. This could also be a fully data assimilative model, which assimilates observations into the model solution.

Predictability vs. Prediction skill

(Potential) Predictability: This characterizes the “ability to be predicted” rather than the current “ability to predict.” One acquires this by computing a metric (like the anomaly correlation coefficient (ACC)) between the prediction ensemble and a verification member (in a perfect-model setup) or the reconstruction that initialized it (in a hindcast setup) [Meehl2013].

(Prediction) Skill: This characterizes the current ability of the ensemble forecasting system to predict the real world. This is derived by computing the ACC between the prediction ensemble and observations of the real world [Meehl2013].

Forecasting

Hindcast: Retrospective forecasts of the past initialized from a reconstruction integrated under external forcing [Boer2016].

Prediction: Forecasts initialized from a reconstruction integrated into the future with external forcing [Boer2016].

Projection An estimate of the future climate that is dependent on the externally forced climate response, such as anthropogenic greenhouse gases, aerosols, and volcanic eruptions [Meehl2013].

References

[Griffies1997]

Griffies, S. M., and K. Bryan. “A Predictability Study of Simulated North Atlantic Multidecadal Variability.” Climate Dynamics 13, no. 7–8 (August 1, 1997): 459–87. https://doi.org/10/ch4kc4

[Boer2016]

(1, 2, 3) Boer, G. J., Smith, D. M., Cassou, C., Doblas-Reyes, F., Danabasoglu, G., Kirtman, B., Kushnir, Y., Kimoto, M., Meehl, G. A., Msadek, R., Mueller, W. A., Taylor, K. E., Zwiers, F., Rixen, M., Ruprich-Robert, Y., and Eade, R.: The Decadal Climate Prediction Project (DCPP) contribution to CMIP6, Geosci. Model Dev., 9, 3751-3777, https://doi.org/10.5194/gmd-9-3751-2016, 2016.

[Meehl2013]

(1, 2, 3) Meehl, G. A., Goddard, L., Boer, G., Burgman, R., Branstator, G., Cassou, C., … & Karspeck, A. (2014). Decadal climate prediction: an update from the trenches. Bulletin of the American Meteorological Society, 95(2), 243-267. https://doi.org/10.1175/BAMS-D-12-00241.1.

Baseline Forecasts

To quantify the quality of an initialized forecast, it is useful to judge it against some simple baseline forecast. climpred currently supports a persistence forecast, but future releases will allow computation of other baseline forecasts. Consider opening a Pull Request to get it implemented more quickly.

Persistence Forecast: Whatever is observed at the time of initialization is forecasted to persist into the forecast period [Jolliffe2012]. You can compute this directly via compute_persistence() or as a method of HindcastEnsemble and PerfectModelEnsemble.

Damped Persistence Forecast: (Not Implemented) The amplitudes of the anomalies reduce in time exponentially at a time scale of the local autocorrelation [Yuan2016].

Climatology: (Not Implemented) The average values at the temporal forecast resolution (e.g., annual, monthly) over some long period, which is usually 30 years [Jolliffe2012].

Random Mechanism: (Not Implemented) A probability distribution is assigned to the possible range of the variable being forecasted, and a sequence of forecasts is produced by taking a sequence of independent values from that distribution [Jolliffe2012]. This would be similar to computing an uninitialized forecast, using climpred’s compute_uninitialized() function.

References

[Jolliffe2012]

(1, 2, 3) Jolliffe, Ian T., and David B. Stephenson, eds. Forecast verification: a practitioner’s guide in atmospheric science. John Wiley & Sons, 2012.

[Yuan2016]

Yuan, Xiaojun, et al. “Arctic sea ice seasonal prediction by a linear Markov model.” Journal of Climate 29.22 (2016): 8151-8173.

Help & Reference

• API Reference
• What’s New
• Helpful Links
• Publications Using climpred
• Contribution Guide
• Release Procedure
• Contributors

API Reference

This page provides an auto-generated summary of climpred’s API. For more details and examples, refer to the relevant chapters in the main part of the documentation.

High-Level Classes

A primary feature of climpred is our prediction ensemble objects, HindcastEnsemble and PerfectModelEnsemble. Users can append their initialized ensemble to these classes, as well as an arbitrary number of references (assimilations, reconstructions, observations), control runs, and uninitialized ensembles.

HindcastEnsemble

A HindcastEnsemble is a prediction ensemble that is initialized off of some form of observations (an assimilation, renanalysis, etc.). Thus, it is anticipated that forecasts are verified against observation-like references. Read more about the terminology here.

HindcastEnsemble(xobj)

An object for climate prediction ensembles initialized by a data-like product.

Add and Retrieve Datasets

HindcastEnsemble.__init__(xobj)	Create a HindcastEnsemble object by inputting output from a prediction ensemble in xarray format.
HindcastEnsemble.add_reference(xobj, name)	Add a reference product for comparison to the initialized ensemble.
HindcastEnsemble.add_uninitialized(xobj)	Add a companion uninitialized ensemble for comparison to references.
HindcastEnsemble.get_initialized()	Returns the xarray dataset for the initialized ensemble.
HindcastEnsemble.get_reference([name])	Returns the given reference(s).
HindcastEnsemble.get_uninitialized()	Returns the xarray dataset for the uninitialized ensemble.

Analysis Functions

HindcastEnsemble.compute_metric([refname, …])	Compares the initialized ensemble to a given reference.
HindcastEnsemble.compute_persistence([…])	Compute a simple persistence forecast for a reference.
HindcastEnsemble.compute_uninitialized([…])	Compares the uninitialized ensemble to a given reference.

Pre-Processing

HindcastEnsemble.smooth([smooth_kws])

Smooth all entries of PredictionEnsemble in the same manner to be able to still calculate prediction skill afterwards.

PerfectModelEnsemble

A PerfectModelEnsemble is a prediction ensemble that is initialized off of a control simulation for a number of randomly chosen initialization dates. Thus, forecasts cannot be verified against real-world observations. Instead, they are compared to one another and to the original control run. Read more about the terminology here.

PerfectModelEnsemble(xobj)

An object for “perfect model” climate prediction ensembles.

Add and Retrieve Datasets

PerfectModelEnsemble.__init__(xobj)	Create a PerfectModelEnsemble object by inputting output from the control run in xarray format.
PerfectModelEnsemble.add_control(xobj)	Add the control run that initialized the climate prediction ensemble.
PerfectModelEnsemble.get_initialized()	Returns the xarray dataset for the initialized ensemble.
PerfectModelEnsemble.get_control()	Returns the control as an xarray dataset.
PerfectModelEnsemble.get_uninitialized()	Returns the xarray dataset for the uninitialized ensemble.

Analysis Functions

PerfectModelEnsemble.bootstrap([metric, …])	Bootstrap ensemble simulations with replacement.
PerfectModelEnsemble.compute_metric([…])	Compares the initialized ensemble to the control run.
PerfectModelEnsemble.compute_persistence([…])	Compute a simple persistence forecast for the control run.
PerfectModelEnsemble.compute_uninitialized([…])	Compares the bootstrapped uninitialized run to the control run.

Generate Data

PerfectModelEnsemble.generate_uninitialized()

Generate an uninitialized ensemble by bootstrapping the initialized prediction ensemble.

Direct Function Calls

A user can directly call functions in climpred. This requires entering more arguments, e.g. the initialized ensemble Dataset/xarray.core.dataarray.DataArray directly as well as a reference product or control run. Our objects HindcastEnsemble and PerfectModelEnsemble wrap most of these functions, making the analysis process much simpler. Once we have wrapped all of the functions in their entirety, we will likely depricate the ability to call them directly.

Bootstrap

bootstrap_compute(hind, reference[, hist, …])	Bootstrap compute with replacement.
bootstrap_hindcast(hind, hist, reference[, …])	Bootstrap compute with replacement. Wrapper of
bootstrap_perfect_model(ds, control[, …])	Bootstrap compute with replacement. Wrapper of
bootstrap_uninit_pm_ensemble_from_control(ds, …)	Create a pseudo-ensemble from control run.
bootstrap_uninitialized_ensemble(hind, hist)	Resample uninitialized hindcast from historical members.
dpp_threshold(control[, sig, bootstrap, dim])	Calc DPP significance levels from re-sampled dataset.
varweighted_mean_period_threshold(control[, …])	Calc the variance-weighted mean period significance levels from re-sampled dataset.

Prediction

compute_hindcast(hind, reference[, metric, …])	Compute a predictability skill score against a reference
compute_perfect_model(ds, control[, metric, …])	Compute a predictability skill score for a perfect-model framework simulation dataset.
compute_persistence(hind, reference[, …])	Computes the skill of a persistence forecast from a simulation.
compute_uninitialized(uninit, reference[, …])	Compute a predictability score between an uninitialized ensemble and a reference.

Metrics

Metric(name, function, positive, …[, …])	Master class for all metrics.
_get_norm_factor(comparison)	Get normalization factor with respect to the type of comparison used for

Comparisons

Comparison(name, function, hindcast, …[, …])

Master class for all comparisons.

Statistics

autocorr(ds[, lag, dim, return_p])	Calculate the lagged correlation of time series.
corr(x, y[, dim, lag, return_p])	Computes the Pearson product-moment coefficient of linear correlation.
decorrelation_time(da[, r, dim])	Calculate the decorrelaton time of a time series.
dpp(ds[, dim, m, chunk])	Calculates the Diagnostic Potential Predictability (dpp)
rm_poly(ds, order[, dim])	Returns xarray object with nth-order fit removed.
rm_trend(da[, dim])	Remove linear trend from time series.
varweighted_mean_period(da[, dim])	Calculate the variance weighted mean period of time series based on xrft.power_spectrum.

Tutorial

load_dataset([name, cache, cache_dir, …])

Load example data or a mask from an online repository.

What’s New

climpred v1.2.0 (2019-12-17)

Depreciated

• Abbreviation pval depreciated. Use p_pval for pearson_r_p_value instead. (GH#264) Aaron Spring.

New Features

• Users can now pass a custom metric or comparison to compute functions. (GH#268) Aaron Spring.

  • See user-defined-metrics and user-defined-comparisons.

• New deterministic metrics (see metrics). (GH#264) Aaron Spring.

  • Spearman ranked correlation (spearman_r)
  • Spearman ranked correlation p-value (spearman_r_p_value)
  • Mean Absolute Deviation (mad)
  • Mean Absolute Percent Error (mape)
  • Symmetric Mean Absolute Percent Error (smape)

• Users can now apply arbitrary xarray methods to HindcastEnsemble and PerfectModelEnsemble. (GH#243) Riley X. Brady.

  • See the Prediction Ensemble objects demo page.

• Add “getter” methods to HindcastEnsemble and PerfectModelEnsemble to retrieve xarray datasets from the objects. (GH#243) Riley X. Brady.

  >>> hind = climpred.tutorial.load_dataset('CESM-DP-SST')
  >>> ref = climpred.tutorial.load_dataset('ERSST')
  >>> hindcast = climpred.HindcastEnsemble(hind)
  >>> hindcast = hindcast.add_reference(ref, 'ERSST')
  >>> print(hindcast)
  <climpred.HindcastEnsemble>
  Initialized Ensemble:
      SST      (init, lead, member) float64 ...
  ERSST:
      SST      (time) float32 ...
  Uninitialized:
      None
  >>> print(hindcast.get_initialized())
  <xarray.Dataset>
  Dimensions:  (init: 64, lead: 10, member: 10)
  Coordinates:
  * lead     (lead) int32 1 2 3 4 5 6 7 8 9 10
  * member   (member) int32 1 2 3 4 5 6 7 8 9 10
  * init     (init) float32 1954.0 1955.0 1956.0 1957.0 ... 2015.0 2016.0 2017.0
  Data variables:
      SST      (init, lead, member) float64 ...
  >>> print(hindcast.get_reference('ERSST'))
  <xarray.Dataset>
  Dimensions:  (time: 61)
  Coordinates:
  * time     (time) int64 1955 1956 1957 1958 1959 ... 2011 2012 2013 2014 2015
  Data variables:
      SST      (time) float32 ...
  

• metric_kwargs can be passed to Metric. (GH#264) Aaron Spring.

  • See metric_kwargs under metrics.

Bug Fixes

• compute_metric() doesn’t drop coordinates from the initialized hindcast ensemble anymore. (GH#258) Aaron Spring.
• Metric uacc does not crash when ppp negative anymore. (GH#264) Aaron Spring.
• Update xskillscore to version 0.0.9 to fix all-NaN issue with pearson_r and pearson_r_p_value when there’s missing data. (GH#269) Riley X. Brady.

Internals/Minor Fixes

• Rewrote varweighted_mean_period() based on xrft. Changed time_dim to dim. Function no longer drops coordinates. (GH#258) Aaron Spring

• Add dim='time' in dpp(). (GH#258) Aaron Spring

• Comparisons m2m, m2e rewritten to not stack dims into supervector because this is now done in xskillscore. (GH#264) Aaron Spring

• Add tqdm progress bar to bootstrap_compute(). (GH#244) Aaron Spring

• Remove inplace behavior for HindcastEnsemble and PerfectModelEnsemble. (GH#243) Riley X. Brady

  • See demo page on prediction ensemble objects

• Added tests for chunking with dask. (GH#258) Aaron Spring

• Fix test issues with esmpy 8.0 by forcing esmpy 7.1 (GH#269). Riley X. Brady

• Rewrote metrics and comparisons as classes to accomodate custom metrics and comparisons. (GH#268) Aaron Spring

  • See user-defined-metrics and user-defined-comparisons.

Documentation

• Add examples notebook for temporal and spatial smoothing. (GH#244) Aaron Spring
• Add documentation for computing a metric over a specified dimension. (GH#244) Aaron Spring
• Update API to be more organized with individual function/class pages. (GH#243) Riley X. Brady.
• Add page describing the HindcastEnsemble and PerfectModelEnsemble objects more clearly. (GH#243) Riley X. Brady
• Add page for publications and helpful links. (GH#270) Riley X. Brady.

climpred v1.1.0 (2019-09-23)

Features

• Write information about skill computation to netcdf attributes(GH#213) Aaron Spring
• Temporal and spatial smoothing module (GH#224) Aaron Spring
• Add metrics brier_score, threshold_brier_score and crpss_es (GH#232) Aaron Spring
• Allow compute_hindcast and compute_perfect_model to specify which dimension dim to calculate metric over (GH#232) Aaron Spring

Bug Fixes

• Correct implementation of probabilistic metrics from xskillscore in compute_perfect_model, bootstrap_perfect_model, compute_hindcast and bootstrap_hindcast, now requires xskillscore>=0.05 (GH#232) Aaron Spring

Internals/Minor Fixes

• Rename .stats.DPP to dpp (GH#232) Aaron Spring
• Add matplotlib as a main dependency so that a direct pip installation works (GH#211) Riley X. Brady.
• climpred is now installable from conda-forge (GH#212) Riley X. Brady.
• Fix erroneous descriptions of sample datasets (GH#226) Riley X. Brady.
• Benchmarking time and peak memory of compute functions with asv (GH#231) Aaron Spring

Documentation

• Add scope of package to docs for clarity for users and developers. (GH#235) Riley X. Brady.

climpred v1.0.1 (2019-07-04)

Bug Fixes

• Accomodate for lead-zero within the lead dimension (GH#196) Riley X. Brady.
• Fix issue with adding uninitialized ensemble to HindcastEnsemble object (GH#199) Riley X. Brady.
• Allow max_dof keyword to be passed to compute_metric and compute_persistence for HindcastEnsemble (GH#199) Riley X. Brady.

Internals/Minor Fixes

• Force xskillscore version 0.0.4 or higher to avoid ImportError (GH#204) Riley X. Brady.
• Change max_dfs keyword to max_dof (GH#199) Riley X. Brady.
• Add testing for HindcastEnsemble and PerfectModelEnsemble (GH#199) Riley X. Brady

climpred v1.0.0 (2019-07-03)

climpred v1.0.0 represents the first stable release of the package. It includes HindcastEnsemble and PerfectModelEnsemble objects to perform analysis with. It offers a suite of deterministic and probabilistic metrics that are optimized to be run on single time series or grids of data (e.g., lat, lon, and depth). Currently, climpred only supports annual forecasts.

Features

• Bootstrap prediction skill based on resampling with replacement consistently in ReferenceEnsemble and PerfectModelEnsemble. (GH#128) Aaron Spring
• Consistent bootstrap function for climpred.stats functions via bootstrap_func wrapper. (GH#167) Aaron Spring
• many more metrics: _msss_murphy, _less and probabilistic _crps, _crpss (GH#128) Aaron Spring

Bug Fixes

• compute_uninitialized now trims input data to the same time window. (GH#193) Riley X. Brady
• rm_poly now properly interpolates/fills NaNs. (GH#192) Riley X. Brady

Internals/Minor Fixes

• The climpred version can be printed. (GH#195) Riley X. Brady
• Constants are made elegant and pushed to a separate module. (GH#184) Andrew Huang
• Checks are consolidated to their own module. (GH#173) Andrew Huang

Documentation

• Documentation built extensively in multiple PRs.

climpred v0.3 (2019-04-27)

climpred v0.3 really represents the entire development phase leading up to the version 1 release. This was done in collaboration between Riley X. Brady, Aaron Spring, and Andrew Huang. Future releases will have less additions.

Features

• Introduces object-oriented system to climpred, with classes ReferenceEnsemble and PerfectModelEnsemble. (GH#86) Riley X. Brady

• Expands bootstrapping module for perfect-module configurations. (GH#78, GH#87) Aaron Spring

• Adds functions for computing Relative Entropy (GH#73) Aaron Spring

• Sets more intelligible dimension expectations for climpred (GH#98, GH#105) Riley X. Brady and Aaron Spring:

  • init: initialization dates for the prediction ensemble
  • lead: retrospective forecasts from prediction ensemble; returned dimension for prediction calculations
  • time: time dimension for control runs, references, etc.
  • member: ensemble member dimension.

• Updates open_dataset to display available dataset names when no argument is passed. (GH#123) Riley X. Brady

• Change ReferenceEnsemble to HindcastEnsemble. (GH#124) Riley X. Brady

• Add probabilistic metrics to climpred. (GH#128) Aaron Spring

• Consolidate separate perfect-model and hindcast functions into singular functions. (GH#128) Aaron Spring

• Add option to pass proxy through to open_dataset for firewalled networks. (GH#138) Riley X. Brady

Bug Fixes

• xr_rm_poly can now operate on Datasets and with multiple variables. It also interpolates across NaNs in time series. (GH#94) Andrew Huang
• Travis CI, treon, and pytest all run for automated testing of new features. (GH#98, GH#105, GH#106) Riley X. Brady and Aaron Spring
• Clean up check_xarray decorators and make sure that they work. (GH#142) Andrew Huang
• Ensures that help() returns proper docstring even with decorators. (GH#149) Andrew Huang
• Fixes bootstrap so p values are correct. (GH#170) Aaron Spring

Internals/Minor Fixes

• Adds unit testing for all perfect-model comparisons. (GH#107) Aaron Spring
• Updates CESM-LE uninitialized ensemble sample data to have 34 members. (GH#113) Riley X. Brady
• Adds MPI-ESM hindcast, historical, and assimilation sample data. (GH#119) Aaron Spring
• Replaces check_xarray with a decorator for checking that input arguments are xarray objects. (GH#120) Andrew Huang
• Add custom exceptions for clearer error reporting. (GH#139) Riley X. Brady
• Remove “xr” prefix from stats module. (GH#144) Riley X. Brady
• Add codecoverage for testing. (GH#152) Riley X. Brady
• Update exception messages for more pretty error reporting. (GH#156) Andrew Huang
• Add pre-commit and flake8/black check in CI. (GH#163) Riley X. Brady
• Change loadutils module to tutorial and open_dataset to load_dataset. (GH#164) Riley X. Brady
• Remove predictability horizon function to revisit for v2. (GH#165) Riley X. Brady
• Increase code coverage through more testing. (GH#167) Aaron Spring
• Consolidates checks and constants into modules. (GH#173) Andrew Huang

climpred v0.2 (2019-01-11)

Name changed to climpred, developed enough for basic decadal prediction tasks on a perfect-model ensemble and reference-based ensemble.

climpred v0.1 (2018-12-20)

Collaboration between Riley Brady and Aaron Spring begins.

Helpful Links

We hope to curate in the climpred documentation a comprehensive report of terminology, best practices, analysis methods, etc. in the prediction community. Here we suggest other resources for initialized prediction of the Earth system to round out the information provided in our documentation.

Forecast Verification

• CAWCR Forecast Verification Overview: A nice overview of forecast verification, including a suite of metrics and their derivation.

Publications Using climpred

Below is a list of publications that have made use of climpred in their analysis. You can nod to climpred, e.g., in your acknowledgements section to help build the community. The main developers of the package intend to release a manuscript documenting climpred in 2020 with a citable DOI, so this can be referenced in the future.

Feel free to open a Pull Request to add your publication to the list!

2019

• Brady, R. X., Lovenduski, N. S., Yeager, S. G., Long, M. C., & Lindsay, K. (2019, October 10). Skillful multiyear predictions of ocean acidification in the California Current System. https://doi.org/10.31223/osf.io/3m2h7

Contribution Guide

Contributions are highly welcomed and appreciated. Every little help counts, so do not hesitate! You can make a high impact on climpred just by using it and reporting issues.

The following sections cover some general guidelines regarding development in climpred for maintainers and contributors. Nothing here is set in stone and can’t be changed. Feel free to suggest improvements or changes in the workflow.

Contribution links

• Contribution Guide
  • Feature requests and feedback
  • Report bugs
  • Fix bugs
  • Write documentation
  • Preparing Pull Requests

Feature requests and feedback

We are eager to hear about your requests for new features and any suggestions about the API, infrastructure, and so on. Feel free to submit these as issues with the label “feature request.”

Please make sure to explain in detail how the feature should work and keep the scope as narrow as possible. This will make it easier to implement in small PRs.

Report bugs

Report bugs for climpred in the issue tracker with the label “bug”.

If you are reporting a bug, please include:

• Your operating system name and version.
• Any details about your local setup that might be helpful in troubleshooting, specifically the Python interpreter version, installed libraries, and climpred version.
• Detailed steps to reproduce the bug.

If you can write a demonstration test that currently fails but should passm that is a very useful commit to make as well, even if you cannot fix the bug itself.

Fix bugs

Look through the GitHub issues for bugs.

Talk to developers to find out how you can fix specific bugs.

Write documentation

climpred could always use more documentation. What exactly is needed?

• More complementary documentation. Have you perhaps found something unclear?
• Docstrings. There can never be too many of them.
• Example notebooks with different Earth System Models, lead times, etc. – they’re all very appreciated.

You can also edit documentation files directly in the GitHub web interface, without using a local copy. This can be convenient for small fixes.

Our documentation is written in reStructuredText. You can follow our conventions in already written documents. Some helpful guides are located here and here.

Note

Build the documentation locally with the following command:

$ conda env update -f ci/environment-dev-3.6.yml
$ cd docs
$ make html

The built documentation should be available in the docs/build/.

If you need to add new functions to the API, run sphinx-autogen -o api api.rst from the docs/source directory and add the functions to api.rst.

Preparing Pull Requests

1. Fork the climpred GitHub repository. It’s fine to use climpred as your fork repository name because it will live under your user.

2. Clone your fork locally using git, connect your repository to the upstream (main project), and create a branch:

   $ git clone git@github.com:YOUR_GITHUB_USERNAME/climpred.git
   $ cd climpred
   $ git remote add upstream git@github.com:bradyrx/climpred.git
   
   # now, to fix a bug or add feature create your own branch off "master":
   
   $ git checkout -b your-bugfix-feature-branch-name master
   

   If you need some help with Git, follow this quick start guide: https://git.wiki.kernel.org/index.php/QuickStart

3. Install dependencies into a new conda environment:

   $ conda env update -f ci/environment-dev-3.7.yml
   $ conda activate climpred-dev
   

4. Make an editable install of climpred by running:

   $ pip install -e .
   

5. Install pre-commit and its hook on the climpred repo:

   $ pip install --user pre-commit
   $ pre-commit install
   

   Afterwards pre-commit will run whenever you commit.

   https://pre-commit.com/ is a framework for managing and maintaining multi-language pre-commit hooks to ensure code-style and code formatting is consistent.

   Now you have an environment called climpred-dev that you can work in. You’ll need to make sure to activate that environment next time you want to use it after closing the terminal or your system.

   You can now edit your local working copy and run/add tests as necessary. Please follow PEP-8 for naming. When committing, pre-commit will modify the files as needed, or will generally be quite clear about what you need to do to pass the commit test.

6. Break your edits up into reasonably sized commits.

   $ git commit -a -m “<commit message>” $ git push -u

7. Run all the tests

   Now running tests is as simple as issuing this command:

   $ coverage run --source climpred -m py.test
   

   This command will run tests via the “pytest” tool against Python 3.6.

8. Create a new changelog entry in CHANGELOG.rst:

   • The entry should be entered as:

   <description> (:pr:`#<pull request number>`) `<author's names>`_

   where <description> is the description of the PR related to the change and <pull request number> is the pull request number and <author's names> are your first and last names.

   • Add yourself to list of authors at the end of CHANGELOG.rst file if not there yet, in alphabetical order.

1. Add yourself to the contributors <https://climpred.readthedocs.io/en/latest/contributors.html>_ list via docs/source/contributors.rst.

1. Finally, submit a pull request through the GitHub website using this data:

   head-fork: YOUR_GITHUB_USERNAME/climpred
   compare: your-branch-name
   
   base-fork: bradyrx/climpred
   base: master
   

Note that you can create the Pull Request while you’re working on this. The PR will update as you add more commits. climpred developers and contributors can then review your code and offer suggestions.

Release Procedure

We follow semantic versioning, e.g., v1.0.0. A major version causes incompatible API changes, a minor version adds functionality, and a patch covers bug fixes.

1. Create a new branch release-vX.x.x with the version for the release.

• Update CHANGELOG.rst
• Make sure all new changes, features are reflected in the documentation.

1. Open a new pull request for this branch targeting master

2. After all tests pass and the PR has been approved, merge the PR into master

3. Tag a release and push to github:

   $ git tag -a v1.0.0 -m "Version 1.0.0"
   $ git push origin master --tags
   

4. Build and publish release on PyPI:

   $ git clean -xfd  # remove any files not checked into git
   $ python setup.py sdist bdist_wheel --universal  # build package
   $ twine upload dist/*  # register and push to pypi
   

5. Update the stable branch (used by ReadTheDocs):

   $ git checkout stable
   $ git rebase master
   $ git push -f origin stable
   $ git checkout master
   

6. Update climpred conda-forge feedstock

• Fork climpred-feedstock repository

• Clone this fork and edit recipe:

  $ git clone git@github.com:username/climpred-feedstock.git
  $ cd climpred-feedstock
  $ cd recipe
  $ # edit meta.yaml
  

• Update version
• Get sha256 from pypi.org for climpred
• Fill in the rest of information as described here

• Commit and submit a PR

Contributors

Core Developers

• Riley X. Brady (github)
• Aaron Spring (github)

Contributors

• Andrew Huang (github)

For a list of all the contributions, see the github contribution graph.