Warming of the oceans¶
Roberto De Almeida [email protected]
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:
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:
- calculate the (monthly) seasonal cycle
- remove it from the original time series
- average over the whole globe
- 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
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
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
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]: