Finding Climate Modes with EOFs

Overview

In this notebook, we will identify and plot a few different modes of climate variability with the help of an EOF package that interfaces with Xarray called xeofs.

Prerequisites

Concepts	Importance	Notes
Intro to Xarray	Necessary
Intro to EOFs	Helpful

Time to learn: 30 minutes

Imports

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.path as mpath
from matplotlib.colors import CenteredNorm
from scipy import stats, signal
from cartopy import crs as ccrs, feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import xeofs as xe

El Niño Southern Oscillation (ENSO)

Here we restrict our domain to the equatorial Pacific. Note that ENSO is commonly defined using an index of SST anomaly over a region of the equatorial Pacific (e.g., the Oceanic Niño Index (ONI)) instead of an EOF. You can read more about ENSO here.

ep_ssta_dt = ssta_dt.where((ssta_dt.lat < 30) & (ssta_dt.lat > -30) & (ssta_dt.lon > 120) & (ssta_dt.lon < 290), drop=True)

ep_model = xe.models.EOF(n_modes=4, use_coslat=True)
ep_model.fit(ep_ssta_dt, dim='time')
ep_eofs = ep_model.components()
ep_pcs = ep_model.scores()
ep_expvar = ep_model.explained_variance_ratio()

ep_eofs.plot(col='mode')

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

../_images/1f15188c4fcdd7d81973ffc6764012c5a8e52a6b0302aa3ec544368cde64c5f0.png

ep_pcs.plot(col='mode')

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

../_images/7b22a4ee2afe18512db35758a435953458c7033615349944af40aca614f28618.png

ep_expvar

<xarray.DataArray 'explained_variance_ratio' (mode: 4)> Size: 32B
array([0.32955431, 0.08712222, 0.05315633, 0.04717546])
Coordinates:
  * mode     (mode) int64 32B 1 2 3 4
Attributes: (12/16)
    model:          EOF analysis
    software:       xeofs
    version:        2.3.2
    date:           2024-06-07 22:12:03
    n_modes:        4
    center:         True
    ...             ...
    feature_name:   feature
    random_state:   None
    verbose:        False
    compute:        True
    solver:         auto
    solver_kwargs:  {}

fig, ax = plt.subplots(1, 1, figsize=(10, 2), dpi=130)
plt.fill_between(ep_pcs.time, ep_pcs.isel(mode=0).where(ep_pcs.isel(mode=0) > 0), color='r')
plt.fill_between(ep_pcs.time, ep_pcs.isel(mode=0).where(ep_pcs.isel(mode=0) < 0), color='b')
plt.ylabel('PC')
plt.xlabel('Year')
plt.xlim(ep_pcs.time.min(), ep_pcs.time.max())
plt.grid(linestyle=':')
plt.title('ENSO Index (detrended equatorial Pacific SSTA EOF1)')

Text(0.5, 1.0, 'ENSO Index (detrended equatorial Pacific SSTA EOF1)')

../_images/b088de0e03bb2595986c47a36f06c5260d856945d1de2ecd7d24202ff7978b95.png

Compare to the ONI:

fig, ax = plt.subplots(1, 1, figsize=(10, 2), dpi=130)
plt.fill_between(ep_pcs.time, ep_pcs.isel(mode=0).where(ep_pcs.isel(mode=0) > 0), color='r')
plt.fill_between(ep_pcs.time, ep_pcs.isel(mode=0).where(ep_pcs.isel(mode=0) < 0), color='b')
plt.ylabel('PC')
plt.xlabel('Year')
plt.xlim(ep_pcs.time.sel(time='1950-01').squeeze(), ep_pcs.time.max())
plt.grid(linestyle=':')
plt.title('ENSO Index (detrended equatorial Pacific SSTA EOF1)')

Text(0.5, 1.0, 'ENSO Index (detrended equatorial Pacific SSTA EOF1)')

../_images/b59484f6718f8ad9a784fe7e791d7e45229a8ca6ba7632cd6b8a6c25c1598a7d.png

ONI

Pacific Decadal Oscillation (PDO)

Here we restrict our domain to the North Pacific. You can read more about PDO here.

np_ssta_dt = ssta_dt.where((ssta_dt.lat < 70) & (ssta_dt.lat > 20) & (ssta_dt.lon > 120) & (ssta_dt.lon < 260), drop=True)

np_model = xe.models.EOF(n_modes=4, use_coslat=True)
np_model.fit(np_ssta_dt, dim='time')
np_eofs = np_model.components()
np_pcs = np_model.scores()
np_expvar = np_model.explained_variance_ratio()

np_eofs.plot(col='mode')

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

../_images/094f56da0c6b4163da9c2635bd0f75a0f6f20be74b4fbe5b2bdff2c63ce5b5f4.png

np_pcs.plot(col='mode')

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

../_images/99e946f70e919fa6a7d6b8aca6ab72b23eb94128a56c8835ac689f945b301050.png

np_expvar

<xarray.DataArray 'explained_variance_ratio' (mode: 4)> Size: 32B
array([0.17464065, 0.11074887, 0.072165  , 0.06702556])
Coordinates:
  * mode     (mode) int64 32B 1 2 3 4
Attributes: (12/16)
    model:          EOF analysis
    software:       xeofs
    version:        2.3.2
    date:           2024-06-07 22:12:05
    n_modes:        4
    center:         True
    ...             ...
    feature_name:   feature
    random_state:   None
    verbose:        False
    compute:        True
    solver:         auto
    solver_kwargs:  {}

fig, ax = plt.subplots(1, 1, figsize=(10, 2), dpi=130)
plt.fill_between(np_pcs.time, np_pcs.isel(mode=0).where(np_pcs.isel(mode=0) > 0), color='r')
plt.fill_between(np_pcs.time, np_pcs.isel(mode=0).where(np_pcs.isel(mode=0) < 0), color='b')
plt.plot(np_pcs.time, np_pcs.isel(mode=0).rolling(time=48, center=True).mean(), color='k', linewidth=2)
plt.ylabel('PC')
plt.xlabel('Year')
plt.xlim(np_pcs.time.min(), np_pcs.time.max())
plt.grid(linestyle=':')
plt.title('PDO Index (detrended North Pacific SSTA EOF1)')

Text(0.5, 1.0, 'PDO Index (detrended North Pacific SSTA EOF1)')

../_images/36f370cfc4505edbe82096e01e138f2bc0573fb6faf0637e61d05b1ea438f6fd.png

Summary

In this notebook, we demonstrated a basic workflow for performing an EOF analysis on gridded SST data using the xeofs package. We plotted the PCs associated with ENSO and PDO using deseasonalized, detrended SSTs.

What’s next?

In the future, additional notebooks may use EOFs to recreate published figures, give an overview of other EOF packages, or explore variations of the EOF method.

Resources and references

Scientific description of the ERSSTv5 data set: Huang et al. (2017), doi:10.1175/jcli-d-16-0836.1
xeofs documentation
Paper describing the xeofs software: Rieger et al., (2024) doi:10.21105/joss.06060