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.

[1]:
%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:

[2]:
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.

[3]:
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.

[4]:
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.

[5]:
# 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'.

[6]:
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.

[7]:
# 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.

[8]:
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-11-15 21:17:10
[9]:
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()
_images/quick-start_16_0.png

We can also check error in our forecasts.

[10]:
init = hindcast.compute_metric(metric='rmse')
persistence = hindcast.compute_persistence(metric='rmse')
[11]:
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()
_images/quick-start_19_0.png