#!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/

#!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/

"""
Explanation: 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:

TLT (Temperature Lower Troposphere): MSU channel 2 by subtracting measurements made at different angles from each other
TMT (Temperature Middle Troposphere): MSU channel 2
TTS (Temperature Troposphere Stratosphere): MSU channel 3
TLS (Temperature Lower Stratosphere): MSU channel 4

The AMSU also provides channels 10-14 (datasets available from RSS), which measure temperatures higher in the stratosphere than the highest MSU channel (4).
End of explanation
"""


%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())

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")

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)

"""
Explanation: Weight functions
End of explanation
"""


from netCDF4 import Dataset
from numpy.ma import average

rootgrp = Dataset('data/uat4_tb_v03r03_avrg_chtlt_197812_201504.nc3.nc')
list(rootgrp.variables)

# 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"][:]

"""
Explanation: Netcdf data
End of explanation
"""


S_theta = pi / 36 * sin(pi/144) * cos(latitude*pi/180)
sum(144 * S_theta)-4*pi

"""
Explanation: 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
$$
End of explanation
"""


w_theta = sin(pi/144) * cos(latitude*pi/180)
sum(w_theta)

Tavg = average(brightness_temperature, axis=2)
Tavg = average(Tavg, axis=1, weights=w_theta)

#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

plot(years, Tavg-273.15)
xlabel("Year")
ylabel("T [C]")
title("TLT (Temperature Lower Troposphere)")
show()

"""
Explanation: 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
$$
End of explanation
"""


Tanom = empty(Tavg.shape)
for i in range(12):
    Tanom[i::12] = Tavg[i::12] - average(Tavg[i::12])

"""
Explanation: The temperature oscillates each year. To calculate the "anomaly", we subtract from each month its average temperature:
End of explanation
"""


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

"""
Explanation: We calculate linear fit
End of explanation
"""


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)

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

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()

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()

#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])

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()

"""
Explanation: And compare against official graph + trend. As can be seen, the agreement is perfect:
End of explanation
"""


rootgrp = Dataset('data/uat4_tb_v03r03_anom_chTLT_197812_201308.nc3.nc')
list(rootgrp.variables)

# 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"][:]

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)

"""
Explanation: Appendix
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.
End of explanation
"""


