Satellite Temperature Data

We use raw data that we downloaded from:

http://www.drroyspencer.com/latest-global-temperatures/



In [1]:

    
%pylab inline
import urllib2
from IPython.display import Image









    



Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.



In [2]:

    
# We don't use this dataset:
# open("data/uahncdc.lt", "w").write(urllib2.urlopen("http://vortex.nsstc.uah.edu/data/msu/t2lt/uahncdc.lt").read())



In [3]:

    
# Only execute this if you want to regenerate the downloaded file
# open("data/tltglhmam_5.5", "w").write(urllib2.urlopen("http://vortex.nsstc.uah.edu/data/msu/t2lt/tltglhmam_5.5").read())

Load the data:



In [4]:

    
D = loadtxt("data/tltglhmam_5.5", skiprows=5, comments="DECADAL")



In [5]:

    
years_int = D[:, 0]
months = D[:, 1]
years = years_int + (months-1) / 12
temp = D[:, 2]

Let's first reproduce the temperature graph from the website:



In [6]:

    
def running_average(temp):
    avg = temp.copy()
    n = 6
    for i in range(n, size(temp)-n-1):
        avg[i] = average(temp[i-n:i+n+1])
    avg[:n] = NaN
    avg[-n:] = NaN
    return avg

figure(figsize=(16, 8))
plot(years, temp, "bo-", label="monthly temperatures")
plot(years, running_average(temp), "r-", lw=3, label="running, centered 13-month average")
xlabel("Year")
ylabel("Monthly temperature [$^\circ C$]")
title("NOAA")
xlim([1979, 2014])
ylim([-0.7, 0.8])
grid()
legend(loc="upper left")
show()
Image(filename="data/UAH_LT_1979_thru_June_2013_v5.6.png")









    












    Out[6]:

The graph is the same, so we have read and interpreted the data correctly, including the running average.

Now let's fit a linear trend using the monthly data:



In [7]:

    
from scipy.optimize import curve_fit
Y0 = years[0]
def f(x, a, b):
    return a*(x-Y0)+b
(a, b), pcov = curve_fit(f, years, temp)
adev = sqrt(pcov[0, 0])
bdev = sqrt(pcov[1, 1])
print "par              dev"
print a, adev
print b, bdev









    



par              dev
0.0137645887899 0.000866761831126
-0.232619879518 0.0172750909687



In [8]:

    
figure(figsize=(16, 8))
plot(years, temp, "bo-", label="monthly temperatures")
plot(years, running_average(temp), "r-", lw=3, label="running, centered 13-month average")
plot(years, a*(years-Y0)+b, "k-", lw=2, label="decadal average: $%.3f \pm %.3f$ $^\circ$C / decade" % (a*10, adev*10))
fill_between(years, (a+adev)*(years-Y0)+b+bdev, (a-adev)*(years-Y0)+b-bdev, color="gray")
xlabel("Year")
ylabel("Monthly temperature [$^\circ C$]")
title("NOAA")
xlim([1979, 2014])
ylim([-0.7, 0.8])
grid()
legend(loc="upper left");

The decadal trend is:



In [9]:

    
dTdt = a*10 # C / decade
print "Decadal trend =", dTdt









    



Decadal trend = 0.137645887899