Roberto De Almeida rob@marinexplore.com
This notebook was used to generate the graph for the article Warming Ocean Threatens Sea Life, using data from the ICOADS program. Data used for the plots was downloaded from the Earth System Research Laboratory at NOAA using the wget
Linux utility:
%%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:
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:
%%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
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[Gt=month_reg@mod] 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:
%%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
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,x=@ave,y=@ave,t=1-jan-1900:31-dec-2012@sbx:13] yes? LET std_dev = sst[d=2,x=@ave,y=@ave,t=1-jan-1900:31-dec-2012@sbx: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.
%%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
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[x=@ave,y=@ave,t=1-jan-1900:31-dec-2012@sbx:13] yes? LET atl_sd = if sst_basin eq 0 then sst[d=2] yes? LET atl_std_dev = atl_sd[x=@ave,y=@ave,t=1-jan-1900:31-dec-2012@sbx:13] yes? yes? LET pac_ssta = if sst_basin eq 1 then ssta[d=1] yes? LET pac_mean = pac_ssta[x=@ave,y=@ave,t=1-jan-1900:31-dec-2012@sbx:13] yes? LET pac_sd = if sst_basin eq 1 then sst[d=2] yes? LET pac_std_dev = pac_sd[x=@ave,y=@ave,t=1-jan-1900:31-dec-2012@sbx: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
from pupynere import netcdf_file
import coards
import numpy as np
Using pupynere we load the data from the NetCDF files created using Ferret:
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:
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:
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')