Calculating ENSO with Xarray

Overview

In this notebook, we will:

Load SST data from the CESM2 model
Mask data using .where()
Compute climatologies and anomalies using .groupby()
Use .rolling() to compute moving average
Compute, normalize, and plot the Niño 3.4 Index

Prerequisites

Concepts	Importance	Notes
Introduction to Xarray	Necessary
Computation and Masking	Necessary

Time to learn: 20 minutes

Imports

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import xarray as xr
from pythia_datasets import DATASETS

Compute the anomalies

We first group by month and subtract the mean SST at each point in the Niño 3.4 region. We then compute the weighted average over the region to obtain the anomalies:

gb = tos_nino34.tos.groupby('time.month')
tos_nino34_anom = gb - gb.mean(dim='time')
index_nino34 = tos_nino34_anom.weighted(tos_nino34.areacello).mean(dim=['lat', 'lon'])

Now, smooth the anomalies with a 5-month running mean:

index_nino34_rolling_mean = index_nino34.rolling(time=5, center=True).mean()

index_nino34.plot(size=8)
index_nino34_rolling_mean.plot()
plt.legend(['anomaly', '5-month running mean anomaly'])
plt.title('SST anomaly over the Niño 3.4 region');

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In [9], line 1
----> 1 index_nino34.plot(size=8)
      2 index_nino34_rolling_mean.plot()
      3 plt.legend(['anomaly', '5-month running mean anomaly'])

File /usr/share/miniconda/envs/pythia-book-dev/lib/python3.8/site-packages/xarray/plot/accessor.py:46, in DataArrayPlotAccessor.__call__(self, **kwargs)
     44 @functools.wraps(dataarray_plot.plot, assigned=("__doc__", "__annotations__"))
     45 def __call__(self, **kwargs) -> Any:
---> 46     return dataarray_plot.plot(self._da, **kwargs)

File /usr/share/miniconda/envs/pythia-book-dev/lib/python3.8/site-packages/xarray/plot/dataarray_plot.py:312, in plot(darray, row, col, col_wrap, ax, hue, subplot_kws, **kwargs)
    308     plotfunc = hist
    310 kwargs["ax"] = ax
--> 312 return plotfunc(darray, **kwargs)

File /usr/share/miniconda/envs/pythia-book-dev/lib/python3.8/site-packages/xarray/plot/dataarray_plot.py:517, in line(darray, row, col, figsize, aspect, size, ax, hue, x, y, xincrease, yincrease, xscale, yscale, xticks, yticks, xlim, ylim, add_legend, _labels, *args, **kwargs)
    513 ylabel = label_from_attrs(yplt, extra=y_suffix)
    515 _ensure_plottable(xplt_val, yplt_val)
--> 517 primitive = ax.plot(xplt_val, yplt_val, *args, **kwargs)
    519 if _labels:
    520     if xlabel is not None:

File /usr/share/miniconda/envs/pythia-book-dev/lib/python3.8/site-packages/matplotlib/axes/_axes.py:1664, in Axes.plot(self, scalex, scaley, data, *args, **kwargs)
   1662 lines = [*self._get_lines(*args, data=data, **kwargs)]
   1663 for line in lines:
-> 1664     self.add_line(line)
   1665 if scalex:
   1666     self._request_autoscale_view("x")

File /usr/share/miniconda/envs/pythia-book-dev/lib/python3.8/site-packages/matplotlib/axes/_base.py:2340, in _AxesBase.add_line(self, line)
   2337 if line.get_clip_path() is None:
   2338     line.set_clip_path(self.patch)
-> 2340 self._update_line_limits(line)
   2341 if not line.get_label():
   2342     line.set_label(f'_child{len(self._children)}')

File /usr/share/miniconda/envs/pythia-book-dev/lib/python3.8/site-packages/matplotlib/axes/_base.py:2363, in _AxesBase._update_line_limits(self, line)
   2359 def _update_line_limits(self, line):
   2360     """
   2361     Figures out the data limit of the given line, updating self.dataLim.
   2362     """
-> 2363     path = line.get_path()
   2364     if path.vertices.size == 0:
   2365         return

File /usr/share/miniconda/envs/pythia-book-dev/lib/python3.8/site-packages/matplotlib/lines.py:1031, in Line2D.get_path(self)
   1029 """Return the `~matplotlib.path.Path` associated with this line."""
   1030 if self._invalidy or self._invalidx:
-> 1031     self.recache()
   1032 return self._path

File /usr/share/miniconda/envs/pythia-book-dev/lib/python3.8/site-packages/matplotlib/lines.py:659, in Line2D.recache(self, always)
    657 if always or self._invalidx:
    658     xconv = self.convert_xunits(self._xorig)
--> 659     x = _to_unmasked_float_array(xconv).ravel()
    660 else:
    661     x = self._x

File /usr/share/miniconda/envs/pythia-book-dev/lib/python3.8/site-packages/matplotlib/cbook/__init__.py:1369, in _to_unmasked_float_array(x)
   1367     return np.ma.asarray(x, float).filled(np.nan)
   1368 else:
-> 1369     return np.asarray(x, float)

TypeError: float() argument must be a string or a number, not 'cftime._cftime.DatetimeNoLeap'

Compute the standard deviation of the SST in the Nino 3.4 region, over the entire time period of the data array:

std_dev = tos_nino34.tos.std()
std_dev

<xarray.DataArray 'tos' ()>
array(1.8436452, dtype=float32)

Then we’ll normalize the values by dividing the rolling mean by the standard deviation of the SST in the Niño 3.4 region:

normalized_index_nino34_rolling_mean = index_nino34_rolling_mean / std_dev

Visualize the computed Niño 3.4 index

We will highlight values in excess of \(\pm\)0.5, roughly corresponding to El Niño (warm) and La Niña (cold) events.

fig = plt.figure(figsize=(12, 6))

plt.fill_between(
    normalized_index_nino34_rolling_mean.time.data,
    normalized_index_nino34_rolling_mean.where(
        normalized_index_nino34_rolling_mean >= 0.4
    ).data,
    0.4,
    color='red',
    alpha=0.9,
)
plt.fill_between(
    normalized_index_nino34_rolling_mean.time.data,
    normalized_index_nino34_rolling_mean.where(
        normalized_index_nino34_rolling_mean <= -0.4
    ).data,
    -0.4,
    color='blue',
    alpha=0.9,
)

normalized_index_nino34_rolling_mean.plot(color='black')
plt.axhline(0, color='black', lw=0.5)
plt.axhline(0.4, color='black', linewidth=0.5, linestyle='dotted')
plt.axhline(-0.4, color='black', linewidth=0.5, linestyle='dotted')
plt.title('Niño 3.4 Index');

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In [12], line 3
fig = plt.figure(figsize=(12, 6))
----> 3 plt.fill_between(
   normalized_index_nino34_rolling_mean.time.data,
   normalized_index_nino34_rolling_mean.where(
       normalized_index_nino34_rolling_mean >= 0.4
   ).data,
   0.4,
   color='red',
   alpha=0.9,
)
plt.fill_between(
   normalized_index_nino34_rolling_mean.time.data,
   normalized_index_nino34_rolling_mean.where(
   (...)
   alpha=0.9,
)
normalized_index_nino34_rolling_mean.plot(color='black')

File /usr/share/miniconda/envs/pythia-book-dev/lib/python3.8/site-packages/matplotlib/pyplot.py:2526, in fill_between(x, y1, y2, where, interpolate, step, data, **kwargs)
@_copy_docstring_and_deprecators(Axes.fill_between)
def fill_between(
       x, y1, y2=0, where=None, interpolate=False, step=None, *,
       data=None, **kwargs):
-> 2526     return gca().fill_between(
       x, y1, y2=y2, where=where, interpolate=interpolate, step=step,
       **({"data": data} if data is not None else {}), **kwargs)

File /usr/share/miniconda/envs/pythia-book-dev/lib/python3.8/site-packages/matplotlib/__init__.py:1423, in _preprocess_data.<locals>.inner(ax, data, *args, **kwargs)
@functools.wraps(func)
def inner(ax, *args, data=None, **kwargs):
   if data is None:
-> 1423         return func(ax, *map(sanitize_sequence, args), **kwargs)
   bound = new_sig.bind(ax, *args, **kwargs)
   auto_label = (bound.arguments.get(label_namer)
                 or bound.kwargs.get(label_namer))

File /usr/share/miniconda/envs/pythia-book-dev/lib/python3.8/site-packages/matplotlib/axes/_axes.py:5367, in Axes.fill_between(self, x, y1, y2, where, interpolate, step, **kwargs)
def fill_between(self, x, y1, y2=0, where=None, interpolate=False,
                step=None, **kwargs):
-> 5367     return self._fill_between_x_or_y(
       "x", x, y1, y2,
       where=where, interpolate=interpolate, step=step, **kwargs)

File /usr/share/miniconda/envs/pythia-book-dev/lib/python3.8/site-packages/matplotlib/axes/_axes.py:5272, in Axes._fill_between_x_or_y(self, ind_dir, ind, dep1, dep2, where, interpolate, step, **kwargs)
       kwargs["facecolor"] = \
           self._get_patches_for_fill.get_next_color()
# Handle united data, such as dates
-> 5272 ind, dep1, dep2 = map(
   ma.masked_invalid, self._process_unit_info(
       [(ind_dir, ind), (dep_dir, dep1), (dep_dir, dep2)], kwargs))
for name, array in [
       (ind_dir, ind), (f"{dep_dir}1", dep1), (f"{dep_dir}2", dep2)]:
   if array.ndim > 1:

File /usr/share/miniconda/envs/pythia-book-dev/lib/python3.8/site-packages/numpy/ma/core.py:2366, in masked_invalid(a, copy)
   cls = type(a)
else:
-> 2366     condition = ~(np.isfinite(a))
   cls = MaskedArray
result = a.view(cls)

TypeError: ufunc 'isfinite' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

Summary

We have applied a variety of Xarray’s selection, grouping, and statistical functions to compute and visualize an important climate index.

Resources and References

Niño 3.4 index
Matplotlib’s fill_between method
Matplotlib’s axhline method (see also its analogous axvline method)

Calculating ENSO with Xarray

Overview

Prerequisites

Imports

The Niño 3.4 Index

Select the Niño 3.4 region

Compute the anomalies

Visualize the computed Niño 3.4 index

Summary

Resources and References

About

Give Feedback

Contribute