A Demonstration of Python’s Power

Here’s just a quick demonstration of how to accomplish a pretty interesting task, computing and plotting Oceanic Niño Index, in Python. We won’t explain much of what’s going on here, but just want to show how much you can accomplish in Python.

First we just bring in some tools from a variety of Python libraries, using the import command. Python is really powerful at assembling various tools together like lego bricks.

A whole bunch of imports

%matplotlib inline
import warnings
warnings.simplefilter("ignore") # Silence warnings
import xarray as xr
import numpy as np
from matplotlib import pyplot as plt

Read sea surface temperature data

ds = xr.open_dataset("./data/NOAA_NCDC_ERSST_v3b_SST.nc")
ds

<xarray.Dataset>
Dimensions:  (lat: 89, lon: 180, time: 684)
Coordinates:
  * lat      (lat) float32 -88.0 -86.0 -84.0 -82.0 -80.0 ... 82.0 84.0 86.0 88.0
  * lon      (lon) float32 0.0 2.0 4.0 6.0 8.0 ... 350.0 352.0 354.0 356.0 358.0
  * time     (time) datetime64[ns] 1960-01-15 1960-02-15 ... 2016-12-15
Data variables:
    sst      (time, lat, lon) float32 ...
Attributes:
    Conventions:  IRIDL
    source:       https://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NCDC/.ERSST/...
    history:      extracted and cleaned by Ryan Abernathey for Research Compu...

sst = ds.sst

Plot mean sst along time dimension

sst.mean(dim='time').plot(vmin=-2, vmax=30)

<matplotlib.collections.QuadMesh at 0x12636e9b0>

png

Plot mean sst along time and longitude dimensions

sst.mean(dim=('time', 'lon')).plot()

[<matplotlib.lines.Line2D at 0x12643d080>]

png

Compute Zonal Anomaly

sst_zonal_time_mean = sst.mean(dim=('time', 'lon'))

(sst.mean(dim='lon') - sst_zonal_time_mean).T.plot()

<matplotlib.collections.QuadMesh at 0x10e264908>

png

Plot Data at a specific grid point

sst.sel(lon=230, lat=0, method='nearest').plot()
sst.sel(lon=230, lat=45, method='nearest').plot()

[<matplotlib.lines.Line2D at 0x12671af28>]

png

Compute Climatologies

sst_clim = sst.groupby('time.month').mean(dim='time')
sst_clim.sel(lon=230, lat=45, method='nearest').plot()


[<matplotlib.lines.Line2D at 0x120c5e048>]

png

sst_clim.mean(dim='lon').T.plot.contourf(levels=np.arange(-2,30))

<matplotlib.contour.QuadContourSet at 0x120da6e10>

png

sst_anom = sst.groupby('time.month') - sst_clim
sst_anom.sel(lon=230, lat=45, method='nearest').plot()

[<matplotlib.lines.Line2D at 0x1204b50b8>]

png

sst_anom.std(dim='time').plot()

<matplotlib.collections.QuadMesh at 0x1207ac2b0>

png

Compute El Niño (La Niña) Index

https://www.ncdc.noaa.gov/teleconnections/enso/indicators/sst.php

El Niño (La Niña) is a phenomenon in the equatorial Pacific Ocean characterized by a five consecutive 3-month running mean of sea surface temperature (SST) anomalies in the Niño 3.4 region that is above (below) the threshold of +0.5°C (-0.5°C). This standard of measure is known as the Oceanic Niño Index (ONI).

nino-regions

sst_anom_nino34 = sst_anom.sel(lat=slice(-5, 5), lon=slice(190, 240))
sst_anom_nino34[0].plot()

<matplotlib.collections.QuadMesh at 0x12098db00>

png

sst_anom_nino34_mean = sst_anom_nino34.mean(dim=('lon', 'lat'))
oni = sst_anom_nino34_mean.rolling(time=3).mean(dim='time')

fig, ax = plt.subplots()
sst_anom_nino34_mean.plot(ax=ax, label='raw')
oni.plot(ax=ax, label='smoothed')
ax.grid()

png

# create a categorical  dataarray
nino34 = xr.full_like(oni, 'none', dtype='U4')
nino34[oni >= 0.5] = 'nino'
nino34[oni <= -0.5] = 'nina'
nino34


<xarray.DataArray (time: 684)>
array(['none', 'none', 'none', ..., 'nina', 'nina', 'nina'], dtype='<U4')
Coordinates:
  * time     (time) datetime64[ns] 1960-01-15 1960-02-15 ... 2016-12-15
    month    (time) int64 1 2 3 4 5 6 7 8 9 10 11 ... 2 3 4 5 6 7 8 9 10 11 12

sst_nino_composite = sst_anom.groupby(nino34.rename('nino34')).mean(dim='time')

sst_nino_composite.sel(nino34='nino').plot()

<matplotlib.collections.QuadMesh at 0x120e547f0>

png

sst_nino_composite.sel(nino34='nina').plot()

<matplotlib.collections.QuadMesh at 0x122e98f98>

png

nino_ds = xr.Dataset({'nino34': nino34, 'oni': oni}).drop('month')
nino_ds.to_netcdf('./data/nino34_index.nc')