One can look at the sea ice extent (first part of this notebook) or the sea ice volume (second part of the notebook).
The least amount of ice during the year happens in September, so we plot the sea ice extent (in square km) between 1979 - 2014 in September
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
from scipy.signal import argrelextrema
from scipy.optimize import curve_fit
def running_average(temp, n=6):
avg = empty(size(temp), dtype=temp.dtype)
for i in range(n, size(temp)-n):
avg[i] = average(temp[i-n:i+n+1])
avg[:n] = NaN
avg[-n:] = NaN
return avg
In [3]:
open("data/N_09_area.txt", "w").write(urllib2.urlopen("ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/Sep/N_09_area.txt").read())
In [4]:
data = read_table("data/N_09_area.txt", skipfooter=17, sep="\s+")
data
Out[4]:
Let's plot the extent directly for September together with a linear fit from 1996-2014. I chose the year 1996 because there was a local maximum in the ice extent that year. That way we get as large slope as possible, i.e. ice free summer as soon as possible, assuming linear relationship (which might not hold).
In [5]:
years = data.year.values
yvals = data.extent.values
Y0 = years[0]
def f(x, a, b):
return a*(x-Y0)+b
(a, b), pcov = curve_fit(f, years[17:], yvals[17:])
adev = sqrt(pcov[0, 0])
bdev = sqrt(pcov[1, 1])
print "par dev"
print a, adev
print b, bdev
In [6]:
figure(figsize=(9, 6))
years2 = linspace(1978, 2050, 2)
plot(years2, a*(years2-Y0)+b-1.2, "r-", lw=1)
plot(years, yvals, "kd-", lw=2, label="ice extent")
plot(years2, a*(years2-Y0)+b, "b-", lw=2, label="average")
xlim([1978, 2050])
ylim([0, 9])
grid()
title("Average Monthly Arctic Sea Ice Extent\nSeptember 1979-2014")
ylabel("Extent (million square kilometers)")
xlabel("Year")
legend()
show()
Years at which the lines cross the $x$-axis:
In [7]:
print "blue line:", -b/a + Y0
print "red line:", -(b-1.2)/a + Y0
The average (blue line) points to 2043. This is only the ice extent, which is different from the volume (that depends on the thickness), see below for a volume graph, which seems to have a faster slope. The relationship might be non linear, the fit is not a great one (the volume graph below seems to provide a much better model for the fit). Also, the first ice free summer will happen when there is a local minimum, so one should really fit the local minima instead (red line), which points to 2035.
Data were taken from:
http://psc.apl.uw.edu/research/projects/arctic-sea-ice-volume-anomaly/data/
In [8]:
data = read_table("data/PIOMAS.vol.daily.1979.2015.Current.v2.1.dat", sep="\s+")
data
Out[8]:
In [9]:
years = data.Year.values
days = data["#day"].values
yvals = data.Vol.values
time = years + (days - 1) / 365.
In [10]:
# Find local minima (smooth the function first, because it is noisy)
idx = argrelextrema(running_average(yvals, n=14), less)
In [11]:
Y0 = years[0]
def f(x, a, b):
return a*(x-Y0)+b
(a, b), pcov = curve_fit(f, time[idx], yvals[idx])
adev = sqrt(pcov[0, 0])
bdev = sqrt(pcov[1, 1])
print "par dev"
print a, adev
print b, bdev
In [12]:
figure(figsize=(9, 6))
years2 = linspace(1978, 2050, 2)
plot(years2, a*(years2-Y0)+b-4, "g-", lw=1)
plot(years2, a*(years2-Y0)+b-2.7, "r-", lw=1)
plot(years2, a*(years2-Y0)+b+2.5, "g-", lw=1)
plot(time, yvals, "k-", lw=1, label="ice volume")
plot(time, running_average(yvals, n=365*2), "c-", lw=2, label="ice volume (5 year central average)")
plot(time[idx], yvals[idx], "bo", label="yearly minima")
plot(years2, a*(years2-Y0)+b, "b-", lw=2, label="average of minima")
xlim([1978, 2050])
ylim([0, 35])
grid()
title("Average Daily Arctic Sea Ice Volume\n 1979-2015")
ylabel("Volume [$10^3 \mathrm{km}^3$]")
xlabel("Year")
legend()
show()
Years at which the lines cross the $x$-axis:
In [13]:
print "blue line:", -b/a + Y0
print "red line:", -(b-2.7)/a + Y0
print "green lines:", -(b-4)/a + Y0, -(b+2.5)/a + Y0
The blue line is the average, green lines are the envelope, and the red line is chosen so that it captures the few recent minima (i.e. all the minima that are the biggest anomalies, years 1982, 2010, 2011, 2012, except the 1981 minimum).
The envelope (green lines) points to years 2020 - 2040. The average (blue line) points to 2032. The red line (i.e. if there are some anomalies like in the years 2010-2012) points to 2024.