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
EOFs with NumPy	Helpful

Time to learn: 30 minutes

Imports¶

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import xeofs as xe

Accessing and preparing the data¶

We will use the NOAA Extended Reconstructed Sea Surface Temperature version 5 (ERSSTv5) monthly gridded dataset Huang et al., 2017, which is accessible using OPeNDAP. More information on using OPeNDAP to access NOAA data can be found here.

data_url = 'https://psl.noaa.gov/thredds/dodsC/Datasets/noaa.ersst.v5/sst.mnmean.nc'

sst = xr.open_dataset(data_url).sst
sst

Check that the data looks as expected:

sst.isel(time=0).plot()

Before we modify the data, let’s do an EOF analysis on the whole dataset:

s_model = xe.single.EOF(n_modes=4, use_coslat=True)
s_model.fit(sst, dim='time')
s_eofs = s_model.components()
s_pcs = s_model.scores()
s_expvar = s_model.explained_variance_ratio()

s_eofs.plot(col='mode')

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

s_pcs.plot(col='mode')

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

s_expvar

EOF1 explains 83% of the variance, and the map shows interhemispheric asymmetry. The corresponding PC has a period of one year, which we can see more clearly by only plotting a few years:

s_pcs.sel(mode=1, time=slice('1900', '1903')).plot(figsize=(8, 3))

This mode is showing the seasonal cycle. This is interesting, but it obfuscates other modes. If we want to study the other ways Earth’s climate varies, we should remove the seasonal cycle from our data. Here we compute this (calling it the SST anomaly) by subtracting out the average of each month using Xarray’s .groupby() method:

sst_clim = sst.groupby('time.month')
ssta = sst_clim - sst_clim.mean(dim='time')

The remaining 3 EOFs show a combination of the long-term warming trend, the seasonal cycle (EOF analyses do not cleanly separate physical modes), and other internal variability. The warming trend is also interesting (see the CMIP6 Cookbook), but here we want to pull out some modes of internal/natural variability. We can detrend the data by removing the global average SST anomaly.

def global_average(data):
    weights = np.cos(np.deg2rad(data.lat))
    data_weighted = data.weighted(weights)
    return data_weighted.mean(dim=['lat', 'lon'], skipna=True)

ssta_dt = (ssta - global_average(ssta)).squeeze()

Let’s find the global EOFs again but with the deseasonalized, detrended data:

ds_model = xe.single.EOF(n_modes=4, use_coslat=True)
ds_model.fit(ssta_dt, dim='time')
ds_eofs = ds_model.components()
ds_pcs = ds_model.scores()
ds_expvar = ds_model.explained_variance_ratio()

ds_eofs.plot(col='mode')

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

ds_pcs.plot(col='mode')

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

ds_expvar

Now we can see some modes of variability! EOF1 looks like ENSO or IPO, and EOF2 is probably picking up a pattern of the recent temperature trend where the Southern Ocean and southeastern Pacific are slightly cooling. EOF3 and EOF4 appear to be showing some decadal modes of variability (PDO and maybe AMO), among other things. There is a lot going on in each of these maps, so to get a clearer index of some modes, we can restrict our domain.

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 0x7f09bc06cb00>

ep_pcs.plot(col='mode')

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

ep_expvar

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)')

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)')

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 0x7f09b12c8290>

np_pcs.plot(col='mode')

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

np_expvar

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)')

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

References¶

Huang, B., Thorne, P. W., Banzon, V. F., Boyer, T., Chepurin, G., Lawrimore, J. H., Menne, M. J., Smith, T. M., Vose, R. S., & Zhang, H.-M. (2017). Extended Reconstructed Sea Surface Temperature, Version 5 (ERSSTv5): Upgrades, Validations, and Intercomparisons. Journal of Climate, 30(20), 8179–8205. 10.1175/JCLI-D-16-0836.1
Rieger, N., & Levang, S. J. (2024). xeofs: Comprehensive EOF analysis in Python with xarray. Journal of Open Source Software, 9(93), 6060. 10.21105/joss.06060

Foundations

EOFs with NumPy