Original data with detailed description is available at: http://nsidc.org/data/docs/noaa/g02135_seaice_index/ We are following the section 3., subsection "Monthly Sea Ice Extent Anomaly Graphs".
In [1]:
%pylab inline
In [2]:
import urllib2
from pandas import read_table, DataFrame
from IPython.display import Image
In [3]:
open("data/N_06_area.txt", "w").write(urllib2.urlopen("ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/Jun/N_06_area.txt").read())
In [4]:
data = read_table("data/N_06_area.txt", skipfooter=16, sep="\s+")
data
Out[4]:
In [5]:
avg = data[2:32].extent.mean()
avg
Out[5]:
In [6]:
anomaly = (data.extent - avg) / avg * 100
Let's calculate linear fit:
In [7]:
d = DataFrame({"anomaly": anomaly.values}, index=data.year)
In [8]:
years = array(d.anomaly.index)
yvals = d.anomaly.values
Y0 = years[0]
def f(x, a, b):
return a*(x-Y0)+b
from scipy.optimize import curve_fit
(a, b), pcov = curve_fit(f, years, yvals)
adev = sqrt(pcov[0, 0])
bdev = sqrt(pcov[1, 1])
print "par dev"
print a, adev
print b, bdev
And plot the results. The first plot is ours, the second is the original:
In [9]:
figure(figsize=(6, 3))
plot(years, a*(years-Y0)+b, "k--", label="x")
plot(years, yvals, "kx-")
xlim([1970, 2020])
ylim([-25, 25])
title("Northern Hemisphere Extent Anomalies Jun 2013")
text(1971, -23, "1981-2010 mean = %.1f million sq km" % avg)
ylabel("%")
xlabel("slope = %.1f(+/-%.1f) %% per decade" % (a*10, adev*10))
show()
Image(url="ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/Jun/N_06_plot.png", embed=True)
Out[9]:
Looks like the fit agrees exactly, only when it is rounded to one decimal digit, we obtained a=3.7, not 3.6, but that's probably just a typo in label for the original graph. Finally, for the standard deviation we got 0.3, not 0.6. That might be caused by using different algorithm for the linear regression.
Let's plot the extent directly for June:
In [10]:
years = data.year.values
yvals = data.extent.values
Y0 = years[0]
def f(x, a, b):
return a*(x-Y0)+b
from scipy.optimize import curve_fit
(a, b), pcov = curve_fit(f, years, yvals)
adev = sqrt(pcov[0, 0])
bdev = sqrt(pcov[1, 1])
print "par dev"
print a, adev
print b, bdev
In [11]:
figure(figsize=(6, 4))
plot(years, a*(years-Y0)+b, "b-", lw=2)
plot(years, yvals, "k.-", lw=2)
xlim([1978, 2013])
title("Average Monthly Arctic Sea Ice Extent\nJune 1979-2013")
ylabel("Extent (million square kilometers)")
xlabel("Year")
show()
And July, which can be compared with another original graph:
In [12]:
open("data/N_07_area.txt", "w").write(urllib2.urlopen("ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/Jul/N_07_area.txt").read())
In [13]:
data = read_table("data/N_07_area.txt", skipfooter=16, sep="\s+")
data
Out[13]:
In [14]:
years = data.year.values
yvals = data.extent.values
Y0 = years[0]
def f(x, a, b):
return a*(x-Y0)+b
from scipy.optimize import curve_fit
(a, b), pcov = curve_fit(f, years, yvals)
adev = sqrt(pcov[0, 0])
bdev = sqrt(pcov[1, 1])
print "par dev"
print a, adev
print b, bdev
In [15]:
figure(figsize=(6, 4))
plot(years, a*(years-Y0)+b, "b-", lw=2)
plot(years, yvals, "k.-", lw=2)
xlim([1978, 2013])
ylim([7, 12])
title("Average Monthly Arctic Sea Ice Extent\nJuly 1979-2013")
ylabel("Extent (million square kilometers)")
xlabel("Year")
show()
Image(url="http://nsidc.org/arcticseaicenews/files/2013/08/Figure3.png", embed=True, width=400)
Out[15]:
This seems to agree exactly.
In [16]:
open("data/S_06_area.txt", "w").write(urllib2.urlopen("ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/Jun/S_06_area.txt").read())
data = read_table("data/S_06_area.txt", skipfooter=7, sep="\s+")
data
Out[16]:
In [17]:
avg = data[2:32].extent.mean()
anomaly = (data.extent - avg) / avg * 100
d = DataFrame({"anomaly": anomaly.values}, index=data.year)
In [18]:
years = array(d.anomaly.index)
yvals = d.anomaly.values
Y0 = years[0]
def f(x, a, b):
return a*(x-Y0)+b
from scipy.optimize import curve_fit
(a, b), pcov = curve_fit(f, years, yvals)
adev = sqrt(pcov[0, 0])
bdev = sqrt(pcov[1, 1])
print "par dev"
print a, adev
print b, bdev
In [19]:
figure(figsize=(6, 3))
plot(years, a*(years-Y0)+b, "k--", label="x")
plot(years, yvals, "kx-")
xlim([1970, 2020])
ylim([-25, 25])
title("Southern Hemisphere Extent Anomalies Jun 2013")
text(1971, -23, "1981-2010 mean = %.1f million sq km" % avg)
ylabel("%")
xlabel("slope = %.1f(+/-%.1f) %% per decade" % (a*10, adev*10))
show()
Image(url="ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/Jun/S_06_plot.png", embed=True)
Out[19]: