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".
%pylab inline
Populating the interactive namespace from numpy and matplotlib
import urllib2
from pandas import read_table, DataFrame
from IPython.display import Image
open("data/N_06_area.txt", "w").write(urllib2.urlopen("ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/Jun/N_06_area.txt").read())
data = read_table("data/N_06_area.txt", skipfooter=16, sep="\s+")
data
year | mo | data_type | region | extent | area | |
---|---|---|---|---|---|---|
0 | 1979 | 6 | Goddard | N | 12.59 | 9.28 |
1 | 1980 | 6 | Goddard | N | 12.31 | 8.93 |
2 | 1981 | 6 | Goddard | N | 12.57 | 8.96 |
3 | 1982 | 6 | Goddard | N | 12.69 | 9.40 |
4 | 1983 | 6 | Goddard | N | 12.36 | 9.20 |
5 | 1984 | 6 | Goddard | N | 12.20 | 9.02 |
6 | 1985 | 6 | Goddard | N | 12.40 | 8.88 |
7 | 1986 | 6 | Goddard | N | 12.10 | 8.93 |
8 | 1987 | 6 | Goddard | N | 12.57 | 9.32 |
9 | 1988 | 6 | Goddard | N | 12.02 | 9.62 |
10 | 1989 | 6 | Goddard | N | 12.31 | 9.90 |
11 | 1990 | 6 | Goddard | N | 11.68 | 9.12 |
12 | 1991 | 6 | Goddard | N | 12.23 | 9.60 |
13 | 1992 | 6 | Goddard | N | 12.13 | 9.89 |
14 | 1993 | 6 | Goddard | N | 11.99 | 9.20 |
15 | 1994 | 6 | Goddard | N | 12.10 | 9.62 |
16 | 1995 | 6 | Goddard | N | 11.55 | 8.86 |
17 | 1996 | 6 | Goddard | N | 12.10 | 9.78 |
18 | 1997 | 6 | Goddard | N | 11.91 | 9.14 |
19 | 1998 | 6 | Goddard | N | 11.85 | 9.11 |
20 | 1999 | 6 | Goddard | N | 12.10 | 9.18 |
21 | 2000 | 6 | Goddard | N | 11.71 | 8.99 |
22 | 2001 | 6 | Goddard | N | 11.69 | 9.01 |
23 | 2002 | 6 | Goddard | N | 11.69 | 9.13 |
24 | 2003 | 6 | Goddard | N | 11.77 | 9.05 |
25 | 2004 | 6 | Goddard | N | 11.51 | 9.18 |
26 | 2005 | 6 | Goddard | N | 11.29 | 8.74 |
27 | 2006 | 6 | Goddard | N | 11.06 | 8.34 |
28 | 2007 | 6 | Goddard | N | 11.49 | 8.15 |
29 | 2008 | 6 | Goddard | N | 11.36 | 8.52 |
30 | 2009 | 6 | Goddard | N | 11.46 | 8.92 |
31 | 2010 | 6 | Goddard | N | 10.82 | 8.02 |
32 | 2011 | 6 | Goddard | N | 10.99 | 8.20 |
33 | 2012 | 6 | Goddard | N | 10.92 | 7.77 |
34 | 2013 | 6 | NRTSI-G | N | 11.58 | 8.57 |
avg = data[2:32].extent.mean()
avg
11.890333333333333
anomaly = (data.extent - avg) / avg * 100
Let's calculate linear fit:
d = DataFrame({"anomaly": anomaly.values}, index=data.year)
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
par dev -0.365124884301 0.0319512189999 5.95201384078 0.631794456031
And plot the results. The first plot is ours, the second is the original:
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)
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:
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
par dev -0.0434145658263 0.00379910644314 12.598047619 0.0751224668034
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:
open("data/N_07_area.txt", "w").write(urllib2.urlopen("ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/Jul/N_07_area.txt").read())
data = read_table("data/N_07_area.txt", skipfooter=16, sep="\s+")
data
year | mo | data_type | region | extent | area | |
---|---|---|---|---|---|---|
0 | 1979 | 7 | Goddard | N | 10.47 | 6.59 |
1 | 1980 | 7 | Goddard | N | 10.39 | 6.45 |
2 | 1981 | 7 | Goddard | N | 10.62 | 6.23 |
3 | 1982 | 7 | Goddard | N | 10.75 | 6.79 |
4 | 1983 | 7 | Goddard | N | 10.91 | 6.79 |
5 | 1984 | 7 | Goddard | N | 10.15 | 6.19 |
6 | 1985 | 7 | Goddard | N | 10.09 | 5.89 |
7 | 1986 | 7 | Goddard | N | 10.47 | 6.53 |
8 | 1987 | 7 | Goddard | N | 9.98 | 6.84 |
9 | 1988 | 7 | Goddard | N | 10.04 | 6.93 |
10 | 1989 | 7 | Goddard | N | 10.38 | 7.24 |
11 | 1990 | 7 | Goddard | N | 9.62 | 6.44 |
12 | 1991 | 7 | Goddard | N | 9.68 | 6.66 |
13 | 1992 | 7 | Goddard | N | 10.61 | 7.14 |
14 | 1993 | 7 | Goddard | N | 9.66 | 6.17 |
15 | 1994 | 7 | Goddard | N | 10.22 | 6.85 |
16 | 1995 | 7 | Goddard | N | 9.15 | 6.05 |
17 | 1996 | 7 | Goddard | N | 10.36 | 7.36 |
18 | 1997 | 7 | Goddard | N | 9.59 | 6.41 |
19 | 1998 | 7 | Goddard | N | 9.62 | 6.38 |
20 | 1999 | 7 | Goddard | N | 9.59 | 6.49 |
21 | 2000 | 7 | Goddard | N | 9.75 | 6.31 |
22 | 2001 | 7 | Goddard | N | 9.22 | 6.22 |
23 | 2002 | 7 | Goddard | N | 9.49 | 6.34 |
24 | 2003 | 7 | Goddard | N | 9.46 | 6.06 |
25 | 2004 | 7 | Goddard | N | 9.60 | 6.43 |
26 | 2005 | 7 | Goddard | N | 8.93 | 5.81 |
27 | 2006 | 7 | Goddard | N | 8.67 | 5.71 |
28 | 2007 | 7 | Goddard | N | 8.13 | 5.03 |
29 | 2008 | 7 | Goddard | N | 8.99 | 5.74 |
30 | 2009 | 7 | Goddard | N | 8.80 | 5.77 |
31 | 2010 | 7 | Goddard | N | 8.36 | 5.27 |
32 | 2011 | 7 | Goddard | N | 7.91 | 5.04 |
33 | 2012 | 7 | Goddard | N | 7.93 | 4.76 |
34 | 2013 | 7 | NRTSI-G | N | 8.45 | 5.41 |
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
par dev -0.0716022408964 0.00637574388393 10.8183809524 0.126072173927
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)
This seems to agree exactly.
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
year | mo | data_type | region | extent | area | |
---|---|---|---|---|---|---|
0 | 1979 | 6 | Goddard | S | 14.74 | 11.13 |
1 | 1980 | 6 | Goddard | S | 13.24 | 9.48 |
2 | 1981 | 6 | Goddard | S | 13.92 | 10.07 |
3 | 1982 | 6 | Goddard | S | 13.99 | 10.41 |
4 | 1983 | 6 | Goddard | S | 13.35 | 9.67 |
5 | 1984 | 6 | Goddard | S | 13.94 | 10.55 |
6 | 1985 | 6 | Goddard | S | 13.99 | 10.42 |
7 | 1986 | 6 | Goddard | S | 13.31 | 9.62 |
8 | 1987 | 6 | Goddard | S | 13.31 | 9.80 |
9 | 1988 | 6 | Goddard | S | 13.60 | 10.19 |
10 | 1989 | 6 | Goddard | S | 14.17 | 10.73 |
11 | 1990 | 6 | Goddard | S | 13.75 | 10.26 |
12 | 1991 | 6 | Goddard | S | 13.52 | 9.87 |
13 | 1992 | 6 | Goddard | S | 13.34 | 9.90 |
14 | 1993 | 6 | Goddard | S | 13.61 | 10.28 |
15 | 1994 | 6 | Goddard | S | 14.13 | 10.63 |
16 | 1995 | 6 | Goddard | S | 13.84 | 10.37 |
17 | 1996 | 6 | Goddard | S | 14.38 | 11.00 |
18 | 1997 | 6 | Goddard | S | 13.67 | 10.12 |
19 | 1998 | 6 | Goddard | S | 13.98 | 10.25 |
20 | 1999 | 6 | Goddard | S | 14.22 | 10.74 |
21 | 2000 | 6 | Goddard | S | 14.41 | 10.95 |
22 | 2001 | 6 | Goddard | S | 13.91 | 10.47 |
23 | 2002 | 6 | Goddard | S | 12.94 | 9.60 |
24 | 2003 | 6 | Goddard | S | 14.48 | 11.07 |
25 | 2004 | 6 | Goddard | S | 14.36 | 10.80 |
26 | 2005 | 6 | Goddard | S | 13.94 | 10.20 |
27 | 2006 | 6 | Goddard | S | 13.98 | 10.49 |
28 | 2007 | 6 | Goddard | S | 13.87 | 10.39 |
29 | 2008 | 6 | Goddard | S | 14.52 | 11.31 |
30 | 2009 | 6 | Goddard | S | 14.41 | 11.01 |
31 | 2010 | 6 | Goddard | S | 15.00 | 11.56 |
32 | 2011 | 6 | Goddard | S | 13.78 | 10.46 |
33 | 2012 | 6 | Goddard | S | 14.20 | 10.70 |
34 | 2013 | 6 | NRTSI-G | S | 14.65 | 11.12 |
avg = data[2:32].extent.mean()
anomaly = (data.extent - avg) / avg * 100
d = DataFrame({"anomaly": anomaly.values}, index=data.year)
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
par dev 0.138044879066 0.0519420666416 -2.14778042176 1.02708787853
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)