<jupyter_start><jupyter_text>S2D Tutorial: Reproduce Fig. 2 of Yeager et al. 2018**NOTE**: See exampels at https://climpred.readthedocs.io/en/stable/ for the latest updates on `climpred`. Many of the methods here might be outdated. For example, `climpred` has released a much more robust smoothing module that can be done with the `HindcastEnsemble` object instead of manually through `xesmf` before constructing the object.This `climpred` lesson analyzes CESM-DPLE SST forecasts to recreate a portion of Figure. 2 from Yeager et al. 2018![](https://i.imgur.com/y6RJE6E.png)<jupyter_code>import climpred

import fsspec
import numpy as np
import xarray as xr
import xesmf as xe

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

import warnings
# Generally not good practice, but we want to supress errors related
# to having NaNs on the grid.
warnings.filterwarnings("ignore")<jupyter_output><empty_output><jupyter_text>Load in CESM-DPLE and ERSSTv5---Note that I pre-processed the CESM-DPLE this by taking the ensemble mean globally to avoid the time it would take to compute live in a workshop. Using `dask_jobqueue` with a few nodes of workers, this only took a few seconds.<jupyter_code># filepath = '/glade/work/rbrady/workshops/climpred/CESM-DP-LE.SST.annmean.ensmean.nc'
# hind = xr.open_dataset(filepath)

# It not on Cheyenne, comment out above and uncomment below.
hind = xr.open_zarr(fsspec.get_mapper('gcs://climpred_workshop/CESM-DP-LE.SST.annmean.ensmean'))

print(hind)<jupyter_output><empty_output><jupyter_text>Rename to the dimension conventions from `climpred`. The package expects at the least a `lead` dimension and `init` dimension.<jupyter_code>hind = hind.rename({'L': 'lead', 'S': 'init'})<jupyter_output><empty_output><jupyter_text>Adjust the `init` dimension to accurately reflect the initialization year. The post-processed output available on Cheyenne uses the convention `S` = 1955 for a November 1954 initialization.<jupyter_code>hind['init'] = hind['init'] - 1<jupyter_output><empty_output><jupyter_text>Load in the ERSSTv5 monthly mean observations.<jupyter_code># filepath = '/glade/work/rbrady/workshops/climpred/ERSSTv5.mnmean.nc'
# obs = xr.open_dataset(filepath)

# It not on Cheyenne, comment out above and uncomment below.
obs = xr.open_zarr(fsspec.get_mapper('gcs://climpred_workshop/ERSSTv5.mnmean'))

obs = obs.rename({'sst': 'SST'})  # Matches variable name in CESM-DPLE
obs = obs.drop_vars(['time_bnds'])
obs.attrs = ''
print(obs)<jupyter_output><empty_output><jupyter_text>Convert observations to annual means.<jupyter_code>obs = obs.groupby('time.year').mean('time').rename({'year': 'time'})<jupyter_output><empty_output><jupyter_text>Generate anomalies by subtracting climatology period used in bias-correcting CESM-DPLE (1964-2014).<jupyter_code>obs = obs - obs.sel(time=slice(1964, 2014)).mean('time')<jupyter_output><empty_output><jupyter_text>Regrid Products to 5x5 Grid---This is done using an external package called `xesmf`, which wraps ESMF. `climpred` has a `smoothing` module that is still in development that will automate this in the future.<jupyter_code># Rename dimensions for what `xesmf` expects.
hind = hind.rename({'nlat': 'y',
                    'nlon': 'x',
                    'TLAT': 'lat',
                    'TLONG': 'lon'})

obs = obs.rename_dims({'lat': 'y', 'lon': 'x'})
hindcast = climpred.HindcastEnsemble(hind)
hindcast = hindcast.add_observations(obs)<jupyter_output><empty_output><jupyter_text>Generate a target grid that is a rectilinear 5x5 degree mesh.<jupyter_code>target_grid = xe.util.grid_global(5, 5)<jupyter_output><empty_output><jupyter_text>Regrid CESM-DPLE following `xesmf` procedure.<jupyter_code>regridder = xe.Regridder(hind, target_grid, 'bilinear', periodic=True)
hind_regridded = regridder(hind)<jupyter_output><empty_output><jupyter_text>Do the same for observations.<jupyter_code>regridder = xe.Regridder(obs, target_grid, 'bilinear', periodic=True)
obs_regridded = regridder(obs)<jupyter_output><empty_output><jupyter_text>We can do a quick check with `xarray`'s native plotting methods to see that the regridded mesh looks as we expect.<jupyter_code>obs_regridded.SST.sel(time=1998).plot()<jupyter_output><empty_output><jupyter_text>Unsmoothed Forecasts---We'll start out by replicating the first two rows of Fig. 2 from the Yeager BAMS paper without applying smoothing. `climpred` is still relatively early in development and does not handle datetime alignment for temporally smoothed persistence forecasts yet.<jupyter_code># Create `HindcastEnsemble` object
hindcast = climpred.HindcastEnsemble(hind_regridded)
hindcast = hindcast.add_observations(obs_regridded)
print(hindcast)<jupyter_output><empty_output><jupyter_text>We can verify the forecasts against the ERSST observations by using the `.verify()` method. Strings for different metrics can be found here: https://climpred.readthedocs.io/en/latest/metrics.html. We use the `e2o` comparison to compare the ensemble mean to observations over the dim `init`. We use the alignment `maximize` to select as many inits as possible that verify at each lead.<jupyter_code>acc_init = hindcast.verify(metric='acc', comparison='e2o', dim='init', alignment='maximize')
pval_init = hindcast.verify(metric='pval', comparison='e2o', dim='init', alignment='maximize')<jupyter_output><empty_output><jupyter_text>ACC plots for individual years.<jupyter_code># Set up 3 subplots as Robinson projections.
f, axs = plt.subplots(ncols=3,
                      subplot_kw=dict(projection=ccrs.Robinson()),
                      figsize=(18, 10))

# Loop through the three axes and leads 3, 5, and 7.
for ax, lead in zip(axs.flatten(), [3, 5, 7]):

    p = ax.pcolormesh(acc_init.lon,
                     acc_init.lat,
                     acc_init['SST'].sel(lead=lead),
                     vmin=-1, vmax=1,  # Colorbar bounds.
                     cmap='RdYlBu_r',
                     transform=ccrs.PlateCarree()) # Transform needed to transform to Robinson projection.

    # Only keep lon/lat with alpha of 0.10 for stippling.
    lons = pval_init['lon'].where(pval_init['SST'] <= 0.10)
    lats = pval_init['lat'].where(pval_init['SST'] <= 0.10)
    ax.scatter(lons.sel(lead=lead),
               lats.sel(lead=lead),
               marker='.', s=1, color='k',
               transform=ccrs.PlateCarree())
    
    ax.add_feature(cfeature.LAND, color='k', zorder=4) # Add land and put on top layer
    ax.set(title=f'Lead Year {lead}')

f.tight_layout()
# Makes one colorbar for the multiple subplots.
f.colorbar(p, ax=axs.ravel().tolist(),
           orientation='horizontal',
           fraction=0.03, pad=0.05,
           label='Anomaly Correlation Coefficient')<jupyter_output><empty_output><jupyter_text>$\Delta$ACC Plots for the same years---p values not included here since we didn't compute z-score significance.<jupyter_code># Compute reference ACC from persistence forecast.
acc = hindcast.verify(metric='acc', comparison='e2o', dim='init', alignment='maximize', reference='persistence')
pval = hindcast.verify(metric='pval', comparison='e2o', dim='init', alignment='maximize', reference='persistence')

# Find delta ACC relative to persistence.
deltaACC = acc.sel(skill='initialized') - acc.sel(skill='persistence')
f, axs = plt.subplots(ncols=3,
                      subplot_kw=dict(projection=ccrs.Robinson()),
                      figsize=(18,10))

for ax, lead in zip(axs.flatten(), [3, 5, 7]):
    p = ax.contourf(deltaACC.lon,
                   deltaACC.lat,
                   deltaACC.SST.sel(lead=lead),
                   levels=np.arange(-0.5, 0.51, 0.05),
                   cmap='RdYlBu_r', extend='both',
                   transform=ccrs.PlateCarree())
    ax.add_feature(cfeature.LAND, color='k', zorder=4)
    ax.set(title=f'Lead Year {lead}')

f.tight_layout()
f.colorbar(p, ax=axs.ravel().tolist(),
           orientation='horizontal',
           fraction=0.03, pad=0.05,
           label='$\Delta$ACC')<jupyter_output><empty_output><jupyter_text>Smoothed Forecasts---Here, we compute the ACC and MSSS for temporally smoothed forecasts. We don't compute the persistence reference for the smoothed forecast, as that functionality is not yet supported. **Note**: We have to make direct calls to the temporal smoothing module, although this will be abstracted away into a method of the `HindcastEnsemble` class in a future version release.<jupyter_code># Smooth in 5-year leads (LY1-5, 2-6, and so on.)
hind_ts = climpred.smoothing.temporal_smoothing(hind_regridded, smooth_kws={'lead': 5}, rename_dim=False)
obs_ts = climpred.smoothing.temporal_smoothing(obs_regridded, smooth_kws={'time': 5}, rename_dim=False)<jupyter_output><empty_output><jupyter_text>Re-instantiate the `HindcastEnsemble` object with the smoothed versions.<jupyter_code>hindcast = climpred.HindcastEnsemble(hind_ts)
hindcast = hindcast.add_observations(obs_ts)
print(hindcast)
acc = hindcast.verify(metric='acc', comparison='e2o', dim='init', alignment='maximize')
pval = hindcast.verify(metric='pval', comparison='e2o', dim='init', alignment='maximize')<jupyter_output><empty_output><jupyter_text>Call utility function to rename lead dimension to readable temporal averages, e.g. "1-5", "2-6". Again, this will be abstracted away in a future version release.<jupyter_code>acc = climpred.smoothing._reset_temporal_axis(acc, smooth_kws={'lead': 5})
pval = climpred.smoothing._reset_temporal_axis(pval, smooth_kws={'lead': 5})<jupyter_output><empty_output><jupyter_text>Plot as in Fig. 2 a-c of Yeager et al. Note that we did not do the FDR test for field significance.<jupyter_code># Set up 3 subplots as Plate Carree projections.
f, axs = plt.subplots(ncols=3,
                      subplot_kw=dict(projection=ccrs.PlateCarree()),
                      figsize=(18,10))

# Loop through the three axes and leads 1-5, 3-7, and 5-9.
for ax, lead in zip(axs.flatten(), ['1-5', '3-7', '5-9']):
    p = ax.pcolormesh(acc.lon,
                   acc.lat,
                   acc.SST.sel(lead=lead),
                   vmin=-1, vmax=1,  # Colorbar bounds.
                   cmap='RdYlBu_r')
    
    # Only keep lon/lat with alpha of 0.10 for stippling.
    lons = pval['lon'].where(pval['SST'] <= 0.10)
    lats = pval['lat'].where(pval['SST'] <= 0.10)
    ax.scatter(lons.sel(lead=lead),
               lats.sel(lead=lead),
               marker='.', s=1, color='k',
               transform=ccrs.PlateCarree())
    
    ax.add_feature(cfeature.LAND, color='k', zorder=4) # Add land and put on top layer.
    ax.set(title=f'Lead Year {lead}')

f.tight_layout()
# Makes one colorbar for the multiple subplots.
f.colorbar(p, ax=axs.ravel().tolist(),
           orientation='horizontal',
           fraction=0.03, pad=0.05,
           label='Anomaly Correlation Coefficient')<jupyter_output><empty_output><jupyter_text>Bonus: FDR Test---We don't yet have field significance tests built into `climpred`, but here is a way you might implement it manually.<jupyter_code>from statsmodels.stats.multitest import multipletests<jupyter_output><empty_output><jupyter_text>We "stack" the three dimensions into one big dimension, since `multipletests` just expects a singular dimension.<jupyter_code>pval = pval.stack(grid=['x', 'y', 'lead']).SST<jupyter_output><empty_output><jupyter_text>The function returns a number of variables, but we just care about the corrected pvalues. Run `help(multipletests)` to see the documentation.<jupyter_code>_, pval_corrected, _, _ = multipletests(pval, alpha=0.1)<jupyter_output><empty_output><jupyter_text>Create `DataArray` matching the dimensions and coordinates of original `pval` DataArray, since `multipletests` returns a `numpy` array.<jupyter_code>new_p = xr.DataArray(pval_corrected, dims=pval.dims, coords=pval.coords)
new_p = new_p.unstack('grid') # Unstack back to 3 dimensions.
# Only keep lon/lat with alpha of 0.10 for stippling.
lons = new_p.lon.where(new_p <= 0.1)
lats = new_p.lat.where(new_p <= 0.1)
# Set up 3 subplots as Plate Carree projections.
f, axs = plt.subplots(ncols=3,
                      subplot_kw=dict(projection=ccrs.PlateCarree()),
                      figsize=(18,10))

# Loop through the three axes and leads 1-5, 3-7, and 5-9.
for ax, lead in zip(axs.flatten(), ['1-5', '3-7', '5-9']):
    p = ax.pcolormesh(acc.lon,
                   acc.lat,
                   acc.SST.sel(lead=lead),
                   vmin=-1, vmax=1,  # Colorbar bounds.
                   cmap='RdYlBu_r')
    ax.scatter(lons.sel(lead=lead),
               lats.sel(lead=lead),
               marker='.', s=1, color='k',
               transform=ccrs.PlateCarree())
    
    ax.add_feature(cfeature.LAND, color='k', zorder=4) # Add land and put on top layer.
    ax.set(title=f'Lead Year {lead}')

f.tight_layout()
# Makes one colorbar for the multiple subplots.
f.colorbar(p, ax=axs.ravel().tolist(),
           orientation='horizontal',
           fraction=0.03, pad=0.05,
           label='Anomaly Correlation Coefficient')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12308279282716521, 'fraction_numerical': 0.015490270888973675, 'mean_word_length': 4.167981072555205, 'pattern_counts': {'":': 0, '<': 103, '<?xml version=': 0, '>': 97, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '325912', 'license': 'MIT', 'n_tokens_mistral': 4606, 'n_tokens_neox': 4319, 'n_words': 1108, 'path': 'notebooks/S2D/Yeager_BAMS.ipynb', 'repo_name': 'bradyrx/climpred_workshop', 'size': '13105'}