In [1]:
%%bash

if [ ! -f sst.mean.nc ]; then
    wget ftp://ftp.cdc.noaa.gov/Datasets/coads/2degree/enh/sst.mean.nc 2>/dev/null
fi

if [ ! -f sst.stddev.nc ]; then
    wget ftp://ftp.cdc.noaa.gov/Datasets/coads/2degree/enh/sst.stddev.nc 2>/dev/null
fi

We now use the Ferret application to:

  1. calculate the (monthly) seasonal cycle
  2. remove it from the original time series
  3. average over the whole globe
  4. smooth with a 13-point box moving average

One important detail about the Ferret averaging is that it is automatically based by area. This is important because the grid cells in the Tropics are larger than the grid cells close to the poles. Simply averaging by number would result in values that are below the true average.

Here we perform steps 1-2, calculating the anomalies from the seasonal cycle:

In [2]:
%%bash

ferret << EOF

USE "sst.mean.nc"

! here we limit the period to 1900-01-01 to 2012-12-31, same as the one 
! used for the Scientific American article
LET sst_interest = sst[d=1,t=1-jan-1900:31-dec-2012]

! create seasonal cycle
USE climatological_axes
CANCEL DATA climatological_axes
LET sstc = sst_interest[[email protected]]

! calculate anomalies from the climatology
LET ssta = sst[d=1] - sstc[gt=sst[d=1]@asn]

set memory/size=1000
save/file=ssta.mean.nc/clobber ssta

quit
EOF
yes? 
yes? USE "sst.mean.nc"
yes? 
yes? ! here we limit the period to 1900-01-01 to 2012-12-31, same as the one 
yes? ! used for the Scientific American article
yes? LET sst_interest = sst[d=1,t=1-jan-1900:31-dec-2012]
yes? 
yes? ! create seasonal cycle
yes? USE climatological_axes
yes? CANCEL DATA climatological_axes
yes? LET sstc = sst_interest[[email protected]]
yes? 
yes? ! calculate anomalies from the climatology
yes? LET ssta = sst[d=1] - sstc[gt=sst[d=1]@asn]
yes? 
yes? set memory/size=1000
yes? save/file=ssta.mean.nc/clobber ssta
yes? 
yes? quit
 	NOAA/PMEL TMAP
 	FERRET v6.82  
 	Linux 2.6.32-279.1.1.el6.x86_64 64-bit - 08/03/12
 	24-Mar-13 20:57     

 *** NOTE: regarding /usr/local/ferret/go/climatological_axes.cdf ...
 *** NOTE: Climatological axes SEASONAL_REG, MONTH_REG, MONTH_IRREG, MONTH_GREGORIAN, MONTH_NOLEAP, MONTH_360_DAY, MONTH_ALL_LEAP  defined
 Cached data cleared from memory
 LISTing to file ssta.mean.nc

And now steps 3-4, averaging over the globe and smoothing the resulting time series with the 13-point box:

In [3]:
%%bash

ferret << EOF

use "ssta.mean.nc"
use "sst.stddev.nc"

! average over x/y, and apply 13-point box moving average
LET mean = ssta[d=1,[email protected],[email protected],t=1-jan-1900:[email protected]:13]
LET std_dev = sst[d=2,[email protected],[email protected],t=1-jan-1900:[email protected]:13]

save/file=global.nc/clobber mean, std_dev

quit
EOF
yes? 
yes? use "ssta.mean.nc"
yes? use "sst.stddev.nc"
yes? 
yes? ! average over x/y, and apply 13-point box moving average
yes? LET mean = ssta[d=1,[email protected],[email protected],t=1-jan-1900:[email protected]:13]
yes? LET std_dev = sst[d=2,[email protected],[email protected],t=1-jan-1900:[email protected]:13]
yes? 
yes? save/file=global.nc/clobber mean, std_dev
yes? 
yes? quit
 	NOAA/PMEL TMAP
 	FERRET v6.82  
 	Linux 2.6.32-279.1.1.el6.x86_64 64-bit - 08/03/12
 	24-Mar-13 20:57     

 LISTing to file global.nc

We now repeat the process for the Atlantic and Pacific basins separately. The Ferret mailing list has an example of how to do that. We need a file which identifies the ocean basins.

In [4]:
%%bash

ferret << EOF

use "ssta.mean.nc"
use "sst.stddev.nc"

! file contaning the different ocean basins; we need to regrid it to
! the same grid as the SST
USE "http://dods.ipsl.jussieu.fr/cgi-bin/nph-dods/ocmip/grids/LEV_grid.nc"
USE "http://dods.ipsl.jussieu.fr/cgi-bin/nph-dods/ocmip/specials/LEV_Basins.nc"
LET lev_basin = reshape(basin[d=LEV_Basins.nc],area[d=LEV_grid.nc])
LET sst_basin = lev_basin[Gxy=ssta[d=1]]

LET atl_ssta = if sst_basin eq 0 then ssta[d=1]
LET atl_mean = atl_ssta[[email protected],[email protected],t=1-jan-1900:[email protected]:13]
LET atl_sd = if sst_basin eq 0 then sst[d=2]
LET atl_std_dev = atl_sd[[email protected],[email protected],t=1-jan-1900:[email protected]:13]

LET pac_ssta = if sst_basin eq 1 then ssta[d=1]
LET pac_mean = pac_ssta[[email protected],[email protected],t=1-jan-1900:[email protected]:13]
LET pac_sd = if sst_basin eq 1 then sst[d=2]
LET pac_std_dev = pac_sd[[email protected],[email protected],t=1-jan-1900:[email protected]:13]

save/file=atlantic.nc/clobber atl_mean, atl_std_dev
save/file=pacific.nc/clobber pac_mean, pac_std_dev

quit
EOF
yes? 
yes? use "ssta.mean.nc"
yes? use "sst.stddev.nc"
yes? 
yes? ! file contaning the different ocean basins; we need to regrid it to
yes? ! the same grid as the SST
yes? USE "http://dods.ipsl.jussieu.fr/cgi-bin/nph-dods/ocmip/grids/LEV_grid.nc"
yes? USE "http://dods.ipsl.jussieu.fr/cgi-bin/nph-dods/ocmip/specials/LEV_Basins.nc"
yes? LET lev_basin = reshape(basin[d=LEV_Basins.nc],area[d=LEV_grid.nc])
yes? LET sst_basin = lev_basin[Gxy=ssta[d=1]]
yes? 
yes? LET atl_ssta = if sst_basin eq 0 then ssta[d=1]
yes? LET atl_mean = atl_ssta[[email protected],[email protected],t=1-jan-1900:[email protected]:13]
yes? LET atl_sd = if sst_basin eq 0 then sst[d=2]
yes? LET atl_std_dev = atl_sd[[email protected],[email protected],t=1-jan-1900:[email protected]:13]
yes? 
yes? LET pac_ssta = if sst_basin eq 1 then ssta[d=1]
yes? LET pac_mean = pac_ssta[[email protected],[email protected],t=1-jan-1900:[email protected]:13]
yes? LET pac_sd = if sst_basin eq 1 then sst[d=2]
yes? LET pac_std_dev = pac_sd[[email protected],[email protected],t=1-jan-1900:[email protected]:13]
yes? 
yes? save/file=atlantic.nc/clobber atl_mean, atl_std_dev
yes? save/file=pacific.nc/clobber pac_mean, pac_std_dev
yes? 
yes? quit
 	NOAA/PMEL TMAP
 	FERRET v6.82  
 	Linux 2.6.32-279.1.1.el6.x86_64 64-bit - 08/03/12
 	24-Mar-13 20:57     

 *** NOTE: unsupported ordering of axes in variable BOUNDS_LON
 *** NOTE: The default ordering will be used
 *** NOTE: unsupported ordering of axes in variable BOUNDS_LAT
 *** NOTE: The default ordering will be used
 LISTing to file atlantic.nc
 LISTing to file pacific.nc

For making the graph, Python was used. In addition to the standard library, the script depends on Numpy (a numerical library for Python), pupynere (a library for reading/writing NetCDF files in Python), and coards (for parsing dates from the NetCDF file).

In [5]:
from pupynere import netcdf_file
import coards
import numpy as np

Using pupynere we load the data from the NetCDF files created using Ferret:

In [6]:
f1 = netcdf_file("global.nc", maskandscale=True)
f2 = netcdf_file("atlantic.nc", maskandscale=True)
f3 = netcdf_file("pacific.nc", maskandscale=True)

# convert time to pylab internal units
time = date2num([coards.parse(v, f1.variables['TIME'].units) for v in f1.variables['TIME']])

mean = f1.variables['MEAN'][:]
std = f1.variables['STD_DEV'][:]
atl_mean = f2.variables['ATL_MEAN'][:]
atl_std = f2.variables['ATL_STD_DEV'][:]
pac_mean = f3.variables['PAC_MEAN'][:]
pac_std = f3.variables['PAC_STD_DEV'][:]

This is a list of El Niño and La Niña events, as well as the date where human population reached 2 to 7 billion:

In [7]:
elnino = [2009, 2006, 2004, 2002, 1997, 1994, 1993, 1992, 1991, 1987, 1982, 1977, 1972, 1969, 1965, 1963, 1957, 1953, 1951, 1946, 1941, 1940, 1932, 1925, 1923, 1919, 1918, 1914, 1913, 1911, 1905, 1902, 1896][:-1]
lanina = [2008, 2007, 2000, 1998, 1996, 1988, 1975, 1974, 1973, 1971, 1970, 1964, 1956, 1955, 1950, 1947, 1938, 1924, 1921, 1917, 1916, 1910, 1909, 1908, 1906, 1893, 1892][:-2]
population = {2: 1927, 3: 1960, 4: 1974, 5: 1987, 6: 1999, 7:2011} 

And this is the plot of the data:

In [8]:
fig = figure(figsize=(15, 5))

ax = fig.gca()
title('Global, Atlantic and Pacific sea surface temperature from 1900 to present')

ax.fill_between(time, mean-std, mean+std, color='k', alpha=0.05)
ax.plot_date(time, atl_mean, 'b-', linewidth=1)
ax.plot_date(time, pac_mean, 'r-', linewidth=1)
ax.plot_date(time, mean, 'w-', linewidth=4)
ax.plot_date(time, mean, 'k-', linewidth=2)
ax.plot_date(time, mean+std, 'w-', linewidth=1)
ax.plot_date(time, mean-std, 'w-', linewidth=1)
ax.set_ylabel(u'°C')

handles = [
        Rectangle((0, 0), 1, 1, fc="k", ec='k', linewidth=2),
        Rectangle((0, 0), 1, 1, fc="b", ec='b', linewidth=1),
        Rectangle((0, 0), 1, 1, fc="r", ec='r', linewidth=1)]
labels = ["Global", "Atlantic", "Pacific"]
legend1 = legend(handles, labels, loc=2, prop={'size':9})

# WW1
x = date2num([datetime.datetime(1914, 7, 28), datetime.datetime(1918, 11, 11)])
y = [-1.5, 1.5]
ax.fill([x[0], x[1], x[1], x[0]], [y[0], y[0], y[1], y[1]], color='k', alpha=0.1, linewidth=0)
ax.text(0.5*(x[0]+x[1]), 1.3, 'WWI', horizontalalignment='center', verticalalignment='center')

# WW2
x = date2num([datetime.datetime(1939, 9, 1), datetime.datetime(1945, 5, 8)])
y = [-1.5, 1.5]
ax.fill([x[0], x[1], x[1], x[0]], [y[0], y[0], y[1], y[1]], color='k', alpha=0.1, linewidth=0)
ax.text(0.5*(x[0]+x[1]), 1.3, 'WWII', horizontalalignment='center', verticalalignment='center')

handles = []
labels = []

# El Ninos
for year in elnino:
    x = date2num(datetime.datetime(year, 12, 31))
    ax.fill([x, x], y, color='r', alpha=0.1, linewidth=2)
handles.append(Rectangle((0, 0), 1, 1, fc="r", alpha=0.1, linewidth=0))
labels.append('El Nino event')

# La Ninas
for year in lanina:
    x = date2num(datetime.datetime(year, 12, 31))
    ax.fill([x, x], y, color='b', alpha=0.1, linewidth=2)
handles.append(Rectangle((0, 0), 1, 1, fc="b", alpha=0.1, linewidth=0))
labels.append('La Nina event')

# population
for billions, year in population.items():
    x = date2num(datetime.datetime(year, 6, 1))
    ax.fill([x, x], y, color='#00aa00', alpha=0.1, linewidth=2)
    ax.text(x, 1.3, '%dB' % billions, horizontalalignment='center', verticalalignment='center', color='#00aa00', alpha=0.9)
handles.append(Rectangle((0, 0), 1, 1, fc="#00aa00", alpha=0.1, linewidth=0))
labels.append('Population')

# Pinatubo 15-06-91
x = date2num(datetime.datetime(1991, 6, 15))
ax.fill([x, x], y, color='k', alpha=0.5, linewidth=2)
handles.append(Rectangle((0, 0), 1, 1, fc="k", alpha=0.5, linewidth=0))
labels.append('Volcanic eruption')

# Katmai-Novarupta 07-06-1912
x = date2num(datetime.datetime(1912, 6, 7))
ax.fill([x, x], y, color='k', alpha=0.5, linewidth=2)

legend(handles, labels, loc=4, prop={'size':9})
ax.add_artist(legend1)

axis((693598.0, 734596.0, -1.5, 1.5))

fig.savefig('fig1a.png')
fig.savefig('fig1a.pdf')
In [8]: