Remote Sensing Systems (RSS, http://www.ssmi.com/) provide machine readable curated datasets of satellite measurements, and the website also explains how they were obtained, processed etc. The temperature data is called MSU (Microwave Sounding Units), that operated between 1978-2005, and AMSU (Advanced Microwave Sounding Units) from 1998. They provide 4 main datasets:
The AMSU also provides channels 10-14 (datasets available from RSS), which measure temperatures higher in the stratosphere than the highest MSU channel (4).
In [1]:
#!wget http://data.remss.com/msu/data/netcdf/uat4_tb_v03r03_avrg_chtlt_197812_201504.nc3.nc
#!mv uat4_tb_v03r03_avrg_chtlt_197812_201504.nc3.nc data/
In [2]:
#!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/
In [3]:
%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())
In [4]:
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")
In [5]:
D = loadtxt("data/std_atmosphere_wt_function_chan_TTS.txt", skiprows=6)
h = D[:, 1]
wTTS = D[:, 5]
In [6]:
D = loadtxt("data/std_atmosphere_wt_function_chan_TLS.txt", skiprows=6)
assert max(abs(h-D[:, 1])) < 1e-12
wTLS = D[:, 5]
In [7]:
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]
In [8]:
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]
In [9]:
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]
In [10]:
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]
In [11]:
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)
Out[11]:
In [12]:
from netCDF4 import Dataset
from numpy.ma import average
In [13]:
rootgrp = Dataset('data/uat4_tb_v03r03_avrg_chtlt_197812_201504.nc3.nc')
list(rootgrp.variables)
Out[13]:
In [14]:
# 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"][:]
We need to calculate the element area (on a unit sphere) as follows: $$ S_{\theta\phi} = \int_{\theta_{min}}^{\theta_{max}} \int_{\phi_{min}}^{\phi_{max}} \sin\theta\, d \theta d \phi = (\cos\theta_{min} - \cos\theta_{max})(\phi_{max} - \phi_{min}) $$ Note that $-180 \le \phi \le 180$ is longitude and $0 \le \theta \le 180$ is something like latitude. Introducing $\Delta\theta = \theta_{max} - \theta_{min}$, $\Delta\phi = \phi_{max} - \phi_{min}$ and $\theta = {\theta_{max} + \theta_{min} \over 2}$ we can write: $$ S_{\theta\phi} = (\cos(\theta-{\Delta\theta\over2}) - \cos(\theta+{\Delta\theta\over 2})) \Delta \phi = 2 \Delta\phi \sin\theta\, \sin{\Delta\theta\over 2} $$ For $\Delta\theta = \Delta\phi = 2.5 {\pi\over 180} = {\pi\over 72}$ we finally obtain: $$ S_\theta = 2 {\pi\over 72} \sin {\pi\over 2\cdot72} \, \sin\theta = {\pi\over 36} \sin {\pi\over 144} \, \sin\theta $$ Finally, we would like to use $\theta$ for latitude, so we need to substitute $\theta \to \theta + {\pi\over 2}$: $$ S_\theta = {\pi\over 36} \sin {\pi\over 144} \, \sin(\theta+{\theta\over 2}) = {\pi\over 36} \sin {\pi\over 144} \, \cos\theta $$
As a check, we calculate the surface of the unit sphere (equal to $4\pi$): $$ \sum_{\theta}144S_\theta = 4\pi $$
In [15]:
S_theta = pi / 36 * sin(pi/144) * cos(latitude*pi/180)
sum(144 * S_theta)-4*pi
Out[15]:
Let's create averaging weights that are normalized to 1 as follows: $$ w_\theta = S_\theta {144\over 4\pi} = \sin{\pi\over144}\cos\theta $$ $$ \sum_\theta w_\theta = 1 $$
In [16]:
w_theta = sin(pi/144) * cos(latitude*pi/180)
sum(w_theta)
Out[16]:
In [17]:
Tavg = average(brightness_temperature, axis=2)
Tavg = average(Tavg, axis=1, weights=w_theta)
In [18]:
#Tavg_ = average(brightness_temperature, axis=2)
#Tavg = ma.empty(Tavg_.shape[0])
#for i in range(size(Tavg_, 0)):
# Tavg[i] = ma.sum(Tavg_[i, :] * w_theta) / sum(w_theta[~Tavg_[i, :].mask])
#print sum(w_theta), sum(w_theta[~Tavg_[i, :].mask]), Tavg_[i, :].mask
In [19]:
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:
In [20]:
Tanom = empty(Tavg.shape)
for i in range(12):
Tanom[i::12] = Tavg[i::12] - average(Tavg[i::12])
We calculate linear fit
In [21]:
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:], p0=(0., 0.))
if pcov == inf:
adev = bdev = 0
else:
adev = sqrt(pcov[0, 0])
bdev = sqrt(pcov[1, 1])
print "par dev"
print a, adev
print b, bdev
And compare against official graph + trend. As can be seen, the agreement is perfect:
In [22]:
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, 2016])
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)
Out[22]:
In [53]:
def maverage(T):
return sum(T)/size(T)
#return average(T)
def mwaverage(T, weights):
return sum(T*weights)
#return average(T, weights=weights)
def maverage(T):
return (sum(T**4)/size(T))**(1./4)
def mwaverage(T, weights):
return sum(T**4 * weights)**(1./4)
def my_average(T):
Tavg = ma.empty(T.shape[:2])
for i in range(size(Tavg, 0)):
for j in range(size(Tavg, 1)):
Tavg[i, j] = maverage(T[i, j, :])
if isnan(Tavg[i, j]):
Tavg[i, j] = ma.masked
#Tavg = average(T, axis=2)
#Tavg = average(Tavg, axis=1, weights=w_theta)
Tavg2 = ma.empty(Tavg.shape[0])
for i in range(size(Tavg2)):
Tavg2[i] = mwaverage(Tavg[i, :], weights=w_theta)
if isnan(Tavg2[i]):
Tavg2[i] = ma.masked
return Tavg2
In [54]:
from utils import arit, arit2, rmean
def doit(longitude, latitude, mafield):
field = asarray(mafield, dtype="double")
mask = asarray(~mafield.mask, dtype="int32")
# arit() gives almost the same results, just a 0.04 off (or more accurate...)
return rmean(asarray(longitude, dtype="double"), asarray(latitude, dtype="double"),
field, mask, -100.)
Tavg2 = ma.empty(brightness_temperature.shape[0])
for i in range(size(brightness_temperature, 0)):
Tavg2[i] = doit(longitude, latitude, brightness_temperature[i, :, :])
if Tavg2[i] == 0 or isnan(Tavg2[i]):
Tavg2[i] = ma.masked
print "abs(Tavg-Tavg2) =", abs(Tavg-Tavg2).max()
In [25]:
plot(years, Tavg-273.15, label="orig")
plot(years, Tavg2-273.15, label="new")
xlabel("Year")
ylabel("T [C]")
legend()
title("TLT (Temperature Lower Troposphere)")
show()
In [26]:
#Tavg2 = my_average(brightness_temperature)
Tanom2 = ma.empty(Tavg2.shape)
for i in range(12):
Tanom2[i::12] = Tavg2[i::12] - average(Tavg2[i::12])
In [27]:
from matplotlib.ticker import MultipleLocator
figure(figsize=(6.6, 3.5))
plot(years, Tanom, "b-", lw=0.7, label="arit. average")
plot(years, Tanom2, "r-", lw=0.7, label="T^{-100} average")
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([-4.2, 4.2])
gca().xaxis.set_minor_locator(MultipleLocator(1))
legend(loc="lower right")
xlabel("Year")
ylabel("Temperature Anomaly [K]")
title("TLT (Temperature Lower Troposphere)")
show()
The anomaly that we did above can also be downloaded directly, we plot it to show that it agrees, but we won't be using it, as we prefer to work with the original data and do all the analysis ourselves for clarity.
In [28]:
rootgrp = Dataset('data/uat4_tb_v03r03_anom_chTLT_197812_201308.nc3.nc')
list(rootgrp.variables)
Out[28]:
In [29]:
# 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"][:]
# 12 values: time in days for 12 months in a year
time_climatology = rootgrp.variables["climatology_time"][:]
# (time, latitude, longitude)
brightness_temperature_anomaly = rootgrp.variables["brightness_temperature_anomaly"][:]
# (time_climatology, latitude, longitude)
brightness_temperature_climatology = rootgrp.variables["brightness_temperature_climatology"][:]
In [30]:
Tavg = average(brightness_temperature_anomaly, axis=2)
Tavg = average(Tavg, axis=1, weights=w_theta)
figure(figsize=(6.6, 3.5))
plot(time / 365.242 + 1978, Tavg, lw=0.5)
xlim([1979, 2014])
ylim([-1.2, 1.2])
xlabel("Year")
ylabel("Temperature Anomaly [K]")
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)
Out[30]:
In [30]: