This notebook tries to reproduce some of the graphs from [1] using their original data.
[1] Rohde et al., Geoinfor Geostat: An Overview 2013, 1:1, http://dx.doi.org/10.4172/gigs.1000101
In [1]:
%pylab inline
from IPython.display import Image
Let's grab raw data for the Berkeley Averaging method:
In [2]:
# Only execute this if you want to regenerate the downloaded file
# open("data/Full_TAVG_complete.txt", "w").write(urllib2.urlopen("http://berkeleyearth.lbl.gov/auto/Global/Full_TAVG_complete.txt").read())
In [3]:
D = loadtxt("data/Full_TAVG_complete.txt", comments="%")
# Read from file:
T0 = 8.79
# Construct floating point years using years + months
years_int = D[:, 0]
months = D[:, 1]
years = years_int + (months-1) / 12
temp = D[:, 2] + T0
unc = D[:, 3]
temp1 = D[:, 4] + T0
unc1 = D[:, 5]
temp5 = D[:, 6] + T0
unc5 = D[:, 7]
temp10 = D[:, 8] + T0
unc10 = D[:, 9]
temp20 = D[:, 10] + T0
unc20 = D[:, 11]
Monthly averages
In [4]:
fill_between(years, temp+unc, temp-unc, color="gray")
plot(years, temp, "k-")
xlabel("Year")
ylabel("Monthly temperature [$^\circ C$]")
title("Berkeley Averaging method")
xlim([1750, 2013])
ylim([0, 15])
grid();
Annual temperature (this should be equivalent to the Figure 5. top, copied as the second graph below for comparison):
In [5]:
fill_between(years, temp1+unc1, temp1-unc1, color="gray")
plot(years, temp1, "k-")
xlabel("Year")
ylabel("Annual temperature [$^\circ C$]")
title("Berkeley Averaging method")
xlim([1750, 2013])
ylim([6, 11])
grid()
show()
Image(filename="data/annual.png")
Out[5]:
Decadal temperature (this should be equivalent to the Figure 5. bottom, copied as the second graph below for comparison):
In [6]:
fill_between(years, temp10+unc10, temp10-unc10, color="gray")
plot(years, temp10, "k-")
xlabel("Year")
ylabel("Decadal temperature [$^\circ C$]")
title("Berkeley Averaging method")
xlim([1750, 2013])
ylim([7, 10])
grid()
show()
Image(filename="data/decadal.png")
Out[6]:
The same plot, but also plot years of volcanic eruptions:
In [7]:
fill_between(years, temp10+unc10, temp10-unc10, color="gray")
plot(years, temp10, "k-")
erupations = [1783, 1815, 1835]
for e in erupations:
plot([e, e], [0, 100], "r-")
xlabel("Year")
ylabel("Decadal temperature [$^\circ C$]")
title("Berkeley Averaging method")
xlim([1750, 2013])
ylim([7, 10])
grid();
Try to estimate the average warming between 1980 and 2010:
In [8]:
from scipy.optimize import curve_fit
idx1980 = sum(years < 1980)
idx2010 = sum(years < 2010)
Y0 = 1980
def f(x, a, b):
return a*(x-Y0)+b
(a, b), pcov = curve_fit(f, years[idx1980:idx2010], temp[idx1980:idx2010], p0=(1, 1), sigma=unc[idx1980:idx2010])
adev = sqrt(pcov[0, 0])
bdev = sqrt(pcov[1, 1])
print "par dev"
print a, adev
print b, bdev
In [9]:
plot(years, temp, "k-")
xlabel("Year")
ylabel("Monthly temperature [$^\circ C$]")
title("Berkeley Averaging method")
plot(years, a*(years-Y0)+b)
fill_between(years, (a+adev)*(years-Y0)+b+bdev, (a-adev)*(years-Y0)+b-bdev, color="gray")
xlim([1980, 2010])
ylim([8, 10.5]);
Temperature difference between 1980 and 2010:
In [10]:
dT = a * 30 # "a" is C/year, so we multiply by 30 years
dT
Out[10]:
In [12]:
dTdt = a*10 # C / decade
dTdt
Out[12]:
In [ ]: