One of the standard global SST products is readily available from NOAA: http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.html. We will use the monthly-averaged version. We need two files: the actual SST, and the land mask.

In [18]:
%matplotlib inline

import os
import time
import urllib

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

from mpl_toolkits.basemap import Basemap

import netCDF4 as nc

from pycurrents.codas import to_day
from pycurrents.num import rangeslice
from pycurrents.num import eof
In [19]:
fnames = ['sst.mnmean.nc', 'lsmask.nc']
urlbase = 'ftp://ftp.cdc.noaa.gov/Datasets/noaa.oisst.v2/'

for fname in fnames:
    if not os.path.exists(fname):
        urllib.request.urlretrieve(urlbase + fname, fname)
        
lsmasknc = nc.Dataset(fnames[1])

# strip the time dimension off the mask
seamask = lsmasknc.variables['mask'][0].astype(bool)
lsmasknc.close()
landmask = np.logical_not(seamask)
print(landmask.shape)

sstnc = nc.Dataset(fnames[0])

# Read all the 1-D dimensional arrays.
time = sstnc.variables['time'][:]
lat = sstnc.variables['lat'][:]                       
lon = sstnc.variables['lon'][:]

# Read the 3-D SST array, and mask it.
sst = np.ma.array(sstnc.variables['sst'][:])
sst[:, landmask] = np.ma.masked
sstnc.close()

(180, 360)
In [20]:
dday2000 = time - to_day(1800, 2000, 1, 1, 0, 0, 0)
yeartime = 2000 + dday2000 / 365.25
print(yeartime[0], yeartime[-1])

1981.9164955509925 2019.1622176591375
In [21]:
# Quick plot of on month against indices (sanity check).
plt.pcolormesh(sst[10, ::-1])
Out[21]:
<matplotlib.collections.QuadMesh at 0x139373f98>

Take out the time mean, and fill the land with zeros.

In [22]:
sstmean = sst.mean(axis=0)
sstdm = sst - sstmean

sstdmz = sstdm.filled(0)

Exercise for the reader: instead of taking out just the mean, take out the mean annual cycle as well.

The grid is uniform in latitude and longitude, so high-latitude cells represent smaller areas than low-latitude cells. We would like the EOF spatial functions to be approximately orthogonal in an area-integrated sense. Because this orthogonality condition involves the product of two EOFs, we weight the data by $\sqrt{\cos\theta}$ so that the product is weighted by $\cos\theta$, where $\theta$ is latitude.

In [23]:
latweights = np.sqrt(np.cos(np.deg2rad(lat)))
latweights = latweights[np.newaxis, :, np.newaxis]
sstdmz *= latweights
In [24]:
ssteof = eof.EOF(sstdmz)
In [25]:
fig, ax = plt.subplots()
ax.plot(ssteof.s, 'o')
ax.margins(0.05)
In [26]:
fig, ax = plt.subplots()
ax.plot(ssteof.percent_var()[:10], 'o')
print(ssteof.percent_var()[:10])

[84.07269878  3.75713812  2.48118546  1.3238733   0.80777787  0.51227754
  0.42653752  0.30876349  0.26307074  0.24013626]
In [27]:
fig, ax = plt.subplots()
ax.plot(yeartime, ssteof.u[:, 5:7])
Out[27]:
[<matplotlib.lines.Line2D at 0x11f015128>,
 <matplotlib.lines.Line2D at 0x11efb8ac8>]

We now take the reshaped spatial patterns, mask them with the landmask, undo the latitude weighting, and multiply them by the singular values.

In [28]:
pats = np.ma.array(ssteof.v_reshaped)
pats[:, landmask] = np.ma.masked
pats /= latweights
pats = pats *  ssteof.s[:, np.newaxis, np.newaxis]
In [29]:
cmap = plt.get_cmap('RdBu_r')

Now make a quick contour plot of the spatial pattern, without bothering to use a map projection.

In [30]:
# quick and dirty: using global pats and ssteof
def contour_pattern(i):
    # normalize so that std of time function is unity, and pattern
    # is in degrees.
    std_t = ssteof.u[:, i].std()
    pat = pats[i] * std_t
    mag = pat.std()
    cticker = mpl.ticker.MaxNLocator(nbins=30, symmetric=True)
    cticker.create_dummy_axis()
    cticker.set_bounds(-2*mag, 2*mag)
    clevs = cticker()

    fig, ax = plt.subplots()
    cs = ax.contourf(lon, lat, pat, levels=clevs,
                 cmap=cmap,
                 extend='both')
    ax.patch.set_facecolor('gray')
    fig.colorbar(cs, ax=ax, shrink=0.9)
In [31]:
contour_pattern(0)

Last, let's make a nicer plot for each of the first few EOFs, with a map for the spatial pattern.

In [32]:
gmap = Basemap(lon_0=180, projection="kav7", resolution='c')
X, Y = gmap(*np.meshgrid(lon, lat))
def draw_map(gmap, ax):
    gmap.drawmapboundary(ax=ax)
    gmap.drawcoastlines(ax=ax)
    gmap.fillcontinents(ax=ax)
    gmap.drawmeridians(np.arange(0, 360.01, 30), ax=ax)
    gmap.drawparallels(np.arange(-60, 60.01, 30), ax=ax)
In [33]:
def map_pattern(i, pats, ssteof, yeartime, gmap):
    # normalize so that std of time function is unity, and pattern
    # is in degrees.
    std_t = ssteof.u[:, i].std()
    pat = pats[i] * std_t
    mag = pat.std()
    cticker = mpl.ticker.MaxNLocator(nbins=30, symmetric=True)
    cticker.create_dummy_axis()
    cticker.set_bounds(-2*mag, 2*mag)
    clevs = cticker()

    #fig, ax = plt.subplots()
    #cs = ax.contourf(lon, lat, pat, levels=clevs,
    #             cmap=cmap,
    #             extend='both')
    #fig.colorbar(cs, ax=ax, shrink=0.9)

    fig, axs = plt.subplots(nrows=2, figsize=(7, 10))
    
    draw_map(gmap, ax=axs[1])
    cs = gmap.contourf(X, Y, pat, levels=clevs, 
                   cmap=cmap,
                   extend='both', ax=axs[1])
    cbar = fig.colorbar(cs, ax=axs[1], orientation='vertical',
                        shrink=0.8,
                        ticks=plt.MaxNLocator(nbins=6, symmetric=True))
    cbar.set_label('$^\circ$C')
    ax = axs[0]
    ax.plot(yeartime, ssteof.u[:, i] / std_t)
    ax.set_title("EOF %d: %.1f%% of variance" % 
                 (i, ssteof.percent_var()[i]))
In [34]:
for i in range(5):
    map_pattern(i, pats, ssteof, yeartime, gmap)