%%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

%%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[Gt=month_reg@mod]

! 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

%%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,x=@ave,y=@ave,t=1-jan-1900:31-dec-2012@sbx:13]
LET std_dev = sst[d=2,x=@ave,y=@ave,t=1-jan-1900:31-dec-2012@sbx:13]

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

quit
EOF

%%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[x=@ave,y=@ave,t=1-jan-1900:31-dec-2012@sbx:13]
LET atl_sd = if sst_basin eq 0 then sst[d=2]
LET atl_std_dev = atl_sd[x=@ave,y=@ave,t=1-jan-1900:31-dec-2012@sbx:13]

LET pac_ssta = if sst_basin eq 1 then ssta[d=1]
LET pac_mean = pac_ssta[x=@ave,y=@ave,t=1-jan-1900:31-dec-2012@sbx:13]
LET pac_sd = if sst_basin eq 1 then sst[d=2]
LET pac_std_dev = pac_sd[x=@ave,y=@ave,t=1-jan-1900:31-dec-2012@sbx:13]

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

quit
EOF

from pupynere import netcdf_file
import coards
import numpy as np

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'][:]

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} 

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