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.

[1]:
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
[2]:
warnings.filterwarnings("ignore")

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

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

[5]:
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')
[5]:
[Text(0.5, 1.0, 'Example Data Coverage')]
../_images/examples_tropical-pacific-ssts_7_1.png

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.

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

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

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

[9]:
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 \alpha = 0.05.

[10]:
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)
[11]:
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())
../_images/examples_tropical-pacific-ssts_18_0.png

Root Mean Square Error of SSTs

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

[12]:
rmse = hindcast.compute_metric(metric='rmse')
[13]:
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=':')
../_images/examples_tropical-pacific-ssts_21_0.png
[ ]: