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 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.single.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 0x7f22ae17d150>

../_images/b4cd7e2b0ac3a0c6f92029eb8eea108ee16ca44886dcc36b993a1a28d7eb8813.png

ep_pcs.plot(col='mode')

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

../_images/c45c6fcd8cc92ea009ed4ad0cbd0ad1b7fd191836be54b2deb187a9e9ae69006.png

ep_expvar

<xarray.DataArray 'explained_variance_ratio' (mode: 4)> Size: 32B
array([0.32945427, 0.08714154, 0.05324459, 0.04714413])
Coordinates:
  * mode     (mode) int64 32B 1 2 3 4
Attributes: (12/15)
    model:          EOF analysis
    software:       xeofs
    version:        3.0.3
    date:           2024-10-17 22:59:23
    n_modes:        4
    center:         True
    ...             ...
    sample_name:    sample
    feature_name:   feature
    random_state:   None
    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/525e4a6f58ebe66ffcddcbd47ece2b62c3c31bac49cc6cc100541b42a0723281.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/07d53283eb78f448e9c5c15e1061c8880c125ea45f03cf9f5931f9308ab3cfc2.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.single.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 0x7f22aba61710>

../_images/b1ba1c34a1dee106dea3a4a6445b2b8bfcf37d534ca1d74ef3eb55be36c5674e.png

np_pcs.plot(col='mode')

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

../_images/48a30a87ce3d0b582f6558030f3f9434feb35edf1016e651a111dd55928cccb9.png

np_expvar

<xarray.DataArray 'explained_variance_ratio' (mode: 4)> Size: 32B
array([0.17530566, 0.11049583, 0.0722191 , 0.06694374])
Coordinates:
  * mode     (mode) int64 32B 1 2 3 4
Attributes: (12/15)
    model:          EOF analysis
    software:       xeofs
    version:        3.0.3
    date:           2024-10-17 22:59:24
    n_modes:        4
    center:         True
    ...             ...
    sample_name:    sample
    feature_name:   feature
    random_state:   None
    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/89bc18eb0ad3b6f20a0d18688e42ff1bcd45632f6a7ab948b8d2048d695cf8fd.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