source: trunk/src/scripts_Laura/test_AMSUB.py @ 45

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import string
import numpy as np
import matplotlib.pyplot as plt
from pylab import *
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import shiftgrid, cm
import netCDF4


########################
# EVOLUTION TEMPORELLE #
########################

## Autour de "Dome C" (lon=123.23;lat=-75.06) - mois: JUNE ##

nb_jr = np.array([31, 28, 31])
## ch1 = 89 GHz - V ##
bbzone_CH1_JAN = nonzero((lon1_JAN >= 120.23) & (lon1_JAN <= 126.23) & (lat1_JAN >= -78.06) & (lat1_JAN <= -72.06))
bbzone_CH1_FEB = nonzero((lon1_FEB >= 120.23) & (lon1_FEB <= 126.23) & (lat1_FEB >= -78.06) & (lat1_FEB <= -72.06))
bbzone_CH1_MAR = nonzero((lon1_MAR >= 120.23) & (lon1_MAR <= 126.23) & (lat1_MAR >= -78.06) & (lat1_MAR <= -72.06))
tb1_JAN_moy = np.zeros([nb_jr[0]],float)
tb1_FEB_moy = np.zeros([nb_jr[1]],float)
tb1_MAR_moy = np.zeros([nb_jr[2]],float)
## ch2 = 150 GHz - V ##
bbzone_CH2_JAN = nonzero((lon2_JAN >= 120.23) & (lon2_JAN <= 126.23) & (lat2_JAN >= -78.06) & (lat2_JAN <= -72.06))
bbzone_CH2_FEB = nonzero((lon2_FEB >= 120.23) & (lon2_FEB <= 126.23) & (lat2_FEB >= -78.06) & (lat2_FEB <= -72.06))
bbzone_CH2_MAR = nonzero((lon2_MAR >= 120.23) & (lon2_MAR <= 126.23) & (lat2_MAR >= -78.06) & (lat2_MAR <= -72.06))
tb2_JAN_moy = np.zeros([nb_jr[0]],float)
tb2_FEB_moy = np.zeros([nb_jr[1]],float)
tb2_MAR_moy = np.zeros([nb_jr[2]],float)

imo = 0 # JANUARY
for ijr in range (0,nb_jr[imo]):
     print 'jour=', ijr+1
     ## ch1 ##
     ind_jr1 = np.where(jjr1_JAN[bbzone_CH1_JAN]==ijr+1)[0]
     b = tb1_JAN[bbzone_CH1_JAN][ind_jr1]
     tb1_JAN_moy[ijr] = mean(b[nonzero((b!=-500.)&(b!=0.))])
     ## ch2 ##
     ind_jr2 = np.where(jjr2_JAN[bbzone_CH2_JAN]==ijr+1)[0]
     e = tb2_JAN[bbzone_CH2_JAN][ind_jr2]
     tb2_JAN_moy[ijr] = mean(e[nonzero((e != -500.) & (e != 0.))])


imo = 1 # FEBUARY
for ijr in range (0,nb_jr[imo]):
     print 'jour=', ijr+1
     ## ch1 ##
     ind_jr1 = np.where(jjr1_FEB[bbzone_CH1_FEB]==ijr+1)[0]
     b = tb1_FEB[bbzone_CH1_FEB][ind_jr1]
     tb1_FEB_moy[ijr] = mean(b[nonzero((b!=-500.)&(b!=0.))])
     ## ch2 ##
     ind_jr2 = np.where(jjr2_FEB[bbzone_CH2_FEB]==ijr+1)[0]
     e = tb2_FEB[bbzone_CH2_FEB][ind_jr2]
     tb2_FEB_moy[ijr] = mean(e[nonzero((e != -500.) & (e != 0.))])


imo = 2 # MARCH
for ijr in range (0,nb_jr[imo]):
     print 'jour=', ijr+1
     ## ch1 ##
     ind_jr1 = np.where(jjr1_MAR[bbzone_CH1_MAR]==ijr+1)[0]
     b = tb1_MAR[bbzone_CH1_MAR][ind_jr1]
     tb1_MAR_moy[ijr] = mean(b[nonzero((b!=-500.)&(b!=0.))])
     ## ch2 ##
     ind_jr2 = np.where(jjr2_MAR[bbzone_CH2_MAR]==ijr+1)[0]
     e = tb2_MAR[bbzone_CH2_MAR][ind_jr2]
     tb2_MAR_moy[ijr] = mean(e[nonzero((e != -500.) & (e != 0.))])


## plot tb en fonction du temps - ch1 et ch2 ##
tb1 = np.append(tb1_JAN_moy, tb1_FEB_moy)
tb1_month = np.append(tb1, tb1_MAR_moy)
tb2 = np.append(tb2_JAN_moy, tb2_FEB_moy)
tb2_month = np.append(tb2, tb2_MAR_moy)
vec_jours = np.arange(0,90,1)
fig = plt.figure()
plt.subplot(2,1,1)
plt.plot(vec_jours, tb1_month, label = '89GHz-V', c = 'r')
plt.plot(vec_jours, tb2_month, label = '150GHz-V', c = 'b')
plt.xticks(np.array([0, 31, 59, 90]), np.array(['JAN', 'FEB', 'MAR']))
grid(True)
ylabel('brightness temperature')
legend()
plt.title('AMSUB - 2010')
## plot grain-index en fonction du temps ##
imo = 0 # JANUARY
grain_ind = 1 - (tb2_month / tb1_month)
plt.subplot(2,1,2)
plt.plot(vec_jours, grain_ind, c = 'k')
plt.xticks(np.array([0, 31, 59, 90]), np.array(['JAN', 'FEB', 'MAR']))
grid(True)
ylabel('Grain-index') 
fig.show()


## difference entre CH1 et CH2 ##
std_1_2 = np.zeros([len(vec_jours)], float)
for ijr in range (0, len(vec_jours)):
     std_1_2[ijr] = sqrt((tb1_month[ijr]**2 - tb2_month[ijr]**2) / len(vec_jours))


fig = plt.figure()
plt.plot(vec_jours, std_1_2, c = 'k')
fig.show()