See the RSS notebook for introduction.
#!wget http://www.remss.com/data/msu/data/netcdf/uat4_tb_v03r03_avrg_chTLT_197812_201308.nc3.nc
#!mv uat4_tb_v03r03_avrg_chTLT_197812_201308.nc3.nc data/
#!wget http://www.remss.com/data/msu/data/netcdf/uat4_tb_v03r03_avrg_chTMT_197812_201308.nc3.nc
#!mv uat4_tb_v03r03_avrg_chTMT_197812_201308.nc3.nc data/
#!wget http://www.remss.com/data/msu/data/netcdf/uat4_tb_v03r03_anom_chTLT_197812_201308.nc3.nc
#!mv uat4_tb_v03r03_anom_chTLT_197812_201308.nc3.nc data/
%pylab inline
import urllib2
import os
from IPython.display import Image
def download(url, dir):
"""Saves file 'url' into 'dir', unless it already exists."""
filename = os.path.basename(url)
fullpath = os.path.join(dir, filename)
if os.path.exists(fullpath):
print "Already downloaded:", filename
else:
print "Downloading:", filename
open(fullpath, "w").write(urllib2.urlopen(url).read())
Populating the interactive namespace from numpy and matplotlib
download("http://www.remss.com/data/msu/weighting_functions/std_atmosphere_wt_function_chan_TTS.txt", "data")
download("http://www.remss.com/data/msu/weighting_functions/std_atmosphere_wt_function_chan_TLS.txt", "data")
download("http://www.remss.com/data/msu/weighting_functions/std_atmosphere_wt_function_chan_tlt_land.txt", "data")
download("http://www.remss.com/data/msu/weighting_functions/std_atmosphere_wt_function_chan_tlt_ocean.txt", "data")
download("http://www.remss.com/data/msu/weighting_functions/std_atmosphere_wt_function_chan_tmt_land.txt", "data")
download("http://www.remss.com/data/msu/weighting_functions/std_atmosphere_wt_function_chan_tmt_ocean.txt", "data")
Already downloaded: std_atmosphere_wt_function_chan_TTS.txt Already downloaded: std_atmosphere_wt_function_chan_TLS.txt Already downloaded: std_atmosphere_wt_function_chan_tlt_land.txt Already downloaded: std_atmosphere_wt_function_chan_tlt_ocean.txt Already downloaded: std_atmosphere_wt_function_chan_tmt_land.txt Already downloaded: std_atmosphere_wt_function_chan_tmt_ocean.txt
D = loadtxt("data/std_atmosphere_wt_function_chan_TTS.txt", skiprows=6)
h = D[:, 1]
wTTS = D[:, 5]
D = loadtxt("data/std_atmosphere_wt_function_chan_TLS.txt", skiprows=6)
assert max(abs(h-D[:, 1])) < 1e-12
wTLS = D[:, 5]
D = loadtxt("data/std_atmosphere_wt_function_chan_tlt_land.txt", skiprows=7)
assert max(abs(h-D[:, 1])) < 1e-12
wTLT_land = D[:, 5]
D = loadtxt("data/std_atmosphere_wt_function_chan_tlt_ocean.txt", skiprows=7)
assert max(abs(h-D[:, 1])) < 1e-12
wTLT_ocean = D[:, 5]
D = loadtxt("data/std_atmosphere_wt_function_chan_tmt_land.txt", skiprows=7)
assert max(abs(h-D[:, 1])) < 1e-12
wTMT_land = D[:, 5]
D = loadtxt("data/std_atmosphere_wt_function_chan_tmt_ocean.txt", skiprows=7)
assert max(abs(h-D[:, 1])) < 1e-12
wTMT_ocean = D[:, 5]
figure(figsize=(3, 8))
plot(wTLS, h/1000, label="TLS")
plot(wTTS, h/1000, label="TTS")
plot(wTMT_ocean, h/1000, label="TMT ocean")
plot(wTMT_land, h/1000, label="TMT land")
plot(wTLT_ocean, h/1000, label="TLT ocean")
plot(wTLT_land, h/1000, label="TLT land")
xlim([0, 0.2])
ylim([0, 50])
legend()
ylabel("Height [km]")
show()
Image(url="http://www.ssmi.com/msu/img/wt_func_plot_for_web_2012.all_channels2.png", embed=True)
<IPython.core.display.Image at 0x26b6e10>
from netCDF4 import Dataset
from numpy.ma import average
rootgrp = Dataset('data/uat4_tb_v03r03_avrg_chTLT_197812_201308.nc3.nc')
list(rootgrp.variables)
[u'longitude', u'latitude', u'time', u'climatology_time', u'longitude_bounds', u'latitude_bounds', u'time_bounds', u'climatology_time_bounds', u'brightness_temperature', u'brightness_temperature_climatology', u'offset_values', u'target_factor_values', u'tb_factor_values', u'msu_amsu_offsets', u'satellites_used']
# 144 values, interval [-180, 180]
longitude = rootgrp.variables["longitude"][:]
# 72 values, interval [-90, 90]
latitude = rootgrp.variables["latitude"][:]
# 144 rows of [min, max]
longitude_bounds = rootgrp.variables["longitude_bounds"][:]
# 72 rows of [min, max]
latitude_bounds = rootgrp.variables["latitude_bounds"][:]
# time in days, 1978 - today
time = rootgrp.variables["time"][:]
# time in years
years = time / 365.242 + 1978
# 12 values: time in days for 12 months in a year
time_climatology = rootgrp.variables["climatology_time"][:]
# (time, latitude, longitude)
brightness_temperature = rootgrp.variables["brightness_temperature"][:]
# (time_climatology, latitude, longitude)
brightness_temperature_climatology = rootgrp.variables["brightness_temperature_climatology"][:]
The weighting function (based on area) is given by: $$ w_\theta = \sin{\pi\over144}\cos\theta $$ $$ \sum_\theta w_\theta = 1 $$
w_theta = sin(pi/144) * cos(latitude*pi/180)
sum(w_theta)
0.99999988
Tavg = average(brightness_temperature, axis=2)
Tavg = average(Tavg, axis=1, weights=w_theta)
plot(years, Tavg-273.15)
xlabel("Year")
ylabel("T [C]")
title("TLT (Temperature Lower Troposphere)")
show()
The temperature oscillates each year. To calculate the "anomaly", we subtract from each month its average temperature:
Tanom = empty(Tavg.shape)
for i in range(12):
Tanom[i::12] = Tavg[i::12] - average(Tavg[i::12])
We calculate linear fit
from scipy.optimize import curve_fit
# Skip the first year, start from 1979, that's why you see the "12" here and below:
Y0 = years[12]
def f(x, a, b):
return a*(x-Y0)+b
(a, b), pcov = curve_fit(f, years[12:], Tanom[12:])
adev = sqrt(pcov[0, 0])
bdev = sqrt(pcov[1, 1])
print "par dev"
print a, adev
print b, bdev
par dev 0.0128077506785 0.000847402784383 -0.220550651996 0.0169189854372
And compare against official graph + trend. As can be seen, the agreement is perfect:
from matplotlib.ticker import MultipleLocator
figure(figsize=(6.6, 3.5))
plot(years, Tanom, "b-", lw=0.7)
plot(years, a*(years-Y0)+b, "b-", lw=0.7, label="Trend = $%.3f \pm %.3f$ K/decade" % (a*10, adev*10))
xlim([1979, 2014])
ylim([-1.2, 1.2])
gca().xaxis.set_minor_locator(MultipleLocator(1))
legend()
xlabel("Year")
ylabel("Temperature Anomaly [K]")
title("TLT (Temperature Lower Troposphere)")
show()
Image(url="http://www.remss.com/data/msu/graphics/TLT/plots/RSS_TS_channel_TLT_Global_Land_And_Sea_v03_3.png", embed=True)
<IPython.core.display.Image at 0x30f3dd0>
Let's show some plots.
from mpl_toolkits.basemap import Basemap
X, Y = meshgrid(longitude, latitude)
idx = 12*30
#X, Y = meshgrid(concatenate([longitude, [longitude[0]]]), concatenate([latitude, [latitude[0]]]))
N = 12
print 'Temperatures (%s)' % int(years[idx])
Tall=(brightness_temperature[idx:idx+N, :, :]-273.15)
vmin = Tall.min()
vmax = Tall.max()
from mpl_toolkits.basemap import Basemap
for i in range(N):
T = brightness_temperature[idx+i, :, :]-273.15
figure(figsize=(18, 6))
subplot(1, 3, 1)
m = Basemap(projection='ortho',lat_0=45, lon_0=-100, resolution='l')
m.drawcoastlines()
m.pcolor(X, Y, T, latlon=True, vmin=vmin, vmax=vmax)
subplot(1, 3, 2)
m = Basemap(projection='hammer',lon_0=0)
m.drawcoastlines()
p = m.pcolor(X, Y, T, latlon=True, vmin=vmin, vmax=vmax)
cbar = m.colorbar(p, location='bottom', pad="5%")
cbar.set_label('$^\circ C$')
subplot(1, 3, 3)
m = Basemap(llcrnrlon=0,llcrnrlat=-90,urcrnrlon=360,urcrnrlat=90,projection='mill')
m.drawcoastlines()
m.drawparallels(np.arange(-80,81,20),labels=[1,1,0,0])
m.drawmeridians(np.arange(0,360,60),labels=[0,0,0,1])
m.pcolor(X, Y, T, latlon=True, vmin=vmin, vmax=vmax)
suptitle("Month: %d" % (i+1))
show()
Temperatures (2008)
Remove seasonality (=calculate anomaly):
T = brightness_temperature
T2 = empty(T.shape)
for i in range(12):
T2[i::12, :, :] = T[i::12, :, :] - average(T[i::12, :, :], axis=0)
m = Basemap(projection='hammer',lon_0=0)
m.drawcoastlines()
p = m.pcolor(X, Y, T2[12*20, :, :], latlon=True)
cbar = m.colorbar(p, location='bottom', pad="5%")
cbar.set_label('$^\circ C$')
Calculate trends:
import numpy.ma as ma
def trend(years, series):
Y0 = years[12]
def f(x, a, b):
return a*(x-Y0)+b
(a, b), pcov = curve_fit(f, years[12:], series[12:])
if isinstance(pcov, float):
return ma.masked
adev = sqrt(pcov[0, 0])
bdev = sqrt(pcov[1, 1])
return a*10
T3 = ma.empty(T2.shape[1:])
print T3.shape
for i in range(size(T2, 1)):
print i
for j in range(size(T2, 2)):
T3[i, j] = trend(years, T2[:, i, j])
(72, 144) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
m = Basemap(projection='hammer',lon_0=0)
m.drawcoastlines()
p = m.pcolor(X, Y, T3, latlon=True, vmin=-0.6, vmax=0.6)
cbar = m.colorbar(p, location='bottom', pad="5%")
cbar.set_label('Trend [$^\circ$C/decade]')
Average over longitude:
T4_TLT = average(T3, axis=1)
plot(latitude, T4_TLT)
[<matplotlib.lines.Line2D at 0x52a6450>]
rootgrp = Dataset('data/uat4_tb_v03r03_avrg_chTMT_197812_201308.nc3.nc')
brightness_temperature = rootgrp.variables["brightness_temperature"][:]
T = brightness_temperature
T2 = empty(T.shape)
for i in range(12):
T2[i::12, :, :] = T[i::12, :, :] - average(T[i::12, :, :], axis=0)
T3 = ma.empty(T2.shape[1:])
print T3.shape
for i in range(size(T2, 1)):
print i
for j in range(size(T2, 2)):
T3[i, j] = trend(years, T2[:, i, j])
T4_TMT = average(T3, axis=1)
(72, 144) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
rootgrp = Dataset('data/uat4_tb_v03r03_avrg_chTLS_197812_201308.nc3.nc')
brightness_temperature = rootgrp.variables["brightness_temperature"][:]
T = brightness_temperature
T2 = empty(T.shape)
for i in range(12):
T2[i::12, :, :] = T[i::12, :, :] - average(T[i::12, :, :], axis=0)
T3 = ma.empty(T2.shape[1:])
print T3.shape
for i in range(size(T2, 1)):
print i
for j in range(size(T2, 2)):
T3[i, j] = trend(years, T2[:, i, j])
T4_TLS = average(T3, axis=1)
(72, 144) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
profile= wTLT_ocean/max(wTLT_ocean) * T4_TLT[:, None]+wTMT_ocean/max(wTMT_ocean) * T4_TMT[:, None]+wTLS/max(wTLS) * T4_TLS[:, None]
This should be compared with Fig. 2H in:
Santer, B. D., Painter, J. F., Bonfils, C., Mears, C. a., Solomon, S., Wigley, T. M. L., … Wentz, F. J. (2013). Human and natural influences on the changing thermal structure of the atmosphere. Proceedings of the National Academy of Sciences, 6–11. doi:10.1073/pnas.1305332110
imshow(profile[:, :50], vmin=-0.3, vmax=0.3)
colorbar()
<matplotlib.colorbar.Colorbar instance at 0x3406e60>