source: trunk/src/scripts_Laura/test_laura.py @ 55

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



###################################################
# definir variables du script 'read_data_test.py' #
###################################################


## sur un zone de l'ANTARC (autour de DomeC) ##
bbzone_CH1_JAN=nonzero((lon1_JAN>=100.)&(lon1_JAN<=150.)&(lat1_JAN>=-85.)&(lat1_JAN<=-65.))
bbzone_CH2_JAN=nonzero((lon2_JAN>=100.)&(lon2_JAN<=150.)&(lat2_JAN>=-85.)&(lat2_JAN<=-65.))
bbzone_CH3_JAN=nonzero((lon3_JAN>=100.)&(lon3_JAN<=150.)&(lat3_JAN>=-85.)&(lat3_JAN<=-65.))
bbzone_CH4_JAN=nonzero((lon4_JAN>=100.)&(lon4_JAN<=150.)&(lat4_JAN>=-85.)&(lat4_JAN<=-65.))
bbzone_CH1_JUL=nonzero((lon1_JUL>=100.)&(lon1_JUL<=150.)&(lat1_JUL>=-85.)&(lat1_JUL<=-65.))
bbzone_CH2_JUL=nonzero((lon2_JUL>=100.)&(lon2_JUL<=150.)&(lat2_JUL>=-85.)&(lat2_JUL<=-65.))
bbzone_CH3_JUL=nonzero((lon3_JUL>=100.)&(lon3_JUL<=150.)&(lat3_JUL>=-85.)&(lat3_JUL<=-65.))
bbzone_CH4_JUL=nonzero((lon4_JUL>=100.)&(lon4_JUL<=150.)&(lat4_JUL>=-85.)&(lat4_JUL<=-65.))
## sur Dome C (lon=123.23;lat=-75.06) - Pour JANUARY ##
ind_lon=np.where(lon_JAN==123.23)[0] 
d=abs(abs(lat_JAN[ind_lon])-75.06)
ind_d=np.where(d==min(d))[0]
ind_lat=lat_JAN[ind_lon][ind_d][0] # coord lat du Dome C
lonDC=lon_JAN[ind_lon][ind_d][0]
latDC=ind_lat
emis_JAN[ind_lon][ind_d]
## sur une zone de l'ANTARC (en altitude) ## 
bbalt_CH1_JAN=nonzero((orog1_JAN>=2000.)&(orog1_JAN<=4000.))
bbalt_CH2_JAN=nonzero((orog2_JAN>=2000.)&(orog2_JAN<=4000.))
bbalt_CH3_JAN=nonzero((orog3_JAN>=2000.)&(orog3_JAN<=4000.))
bbalt_CH4_JAN=nonzero((orog4_JAN>=2000.)&(orog4_JAN<=4000.))
## sur une zone de l'ANTARC (en altitude) ## 
bbalt_CH1_JAN=nonzero((orog1_JAN>=2000.)&(orog1_JAN<=4000.))
bbalt_CH2_JAN=nonzero((orog2_JAN>=2000.)&(orog2_JAN<=4000.))
bbalt_CH3_JAN=nonzero((orog3_JAN>=2000.)&(orog3_JAN<=4000.))
bbalt_CH4_JAN=nonzero((orog4_JAN>=2000.)&(orog4_JAN<=4000.))


nbjr=np.array(31,28,31,30,31,30,31)
i=0 # JANUARY
emis1_JAN_moy=np.zeros([31],float)
emis2_JAN_moy=np.zeros([31],float)
emis3_JAN_moy=np.zeros([31],float)
emis4_JAN_moy=np.zeros([31],float)
#emis1_JUL_moy=np.zeros([31],float)
#emis2_JUL_moy=np.zeros([31],float)
#emis3_JUL_moy=np.zeros([31],float)
#emis4_JUL_moy=np.zeros([31],float)
ts1_JAN_moy=np.zeros([31],float)
ts2_JAN_moy=np.zeros([31],float)
ts3_JAN_moy=np.zeros([31],float)
ts4_JAN_moy=np.zeros([31],float)
#ts1_JUL_moy=np.zeros([31],float)
#ts2_JUL_moy=np.zeros([31],float)
#ts3_JUL_moy=np.zeros([31],float)
#ts4_JUL_moy=np.zeros([31],float)
tb1_JAN_moy=np.zeros([31],float)
tb2_JAN_moy=np.zeros([31],float)
tb3_JAN_moy=np.zeros([31],float)
tb4_JAN_moy=np.zeros([31],float)
#tb1_JUL_moy=np.zeros([31],float)
#tb2_JUL_moy=np.zeros([31],float)
#tb3_JUL_moy=np.zeros([31],float)
#tb4_JUL_moy=np.zeros([31],float)
orog1_JAN_moy=np.zeros([31],float)
orog2_JAN_moy=np.zeros([31],float)
orog3_JAN_moy=np.zeros([31],float)
orog4_JAN_moy=np.zeros([31],float)
#orog1_JUL_moy=np.zeros([31],float)
#orog2_JUL_moy=np.zeros([31],float)
#orog3_JUL_moy=np.zeros([31],float)
#orog4_JUL_moy=np.zeros([31],float)

for ijr in range(0,nbjr[i]):
     print 'jour=', ijr+1
     ind_jr1=np.where(jjr1_JAN[bbalt_CH1_JAN]==ijr+1)[0]
     ind_jr2=np.where(jjr2_JAN[bbalt_CH2_JAN]==ijr+1)[0]
     ind_jr3=np.where(jjr3_JAN[bbalt_CH3_JAN]==ijr+1)[0]
     ind_jr4=np.where(jjr4_JAN[bbalt_CH4_JAN]==ijr+1)[0]
     a=emis1_JAN[bbalt_CH1_JAN][ind_jr1]
     b=emis2_JAN[bbalt_CH2_JAN][ind_jr2]
     c=emis3_JAN[bbalt_CH3_JAN][ind_jr3]
     d=emis4_JAN[bbalt_CH4_JAN][ind_jr4]
     e=ts1_JAN[bbalt_CH1_JAN][ind_jr1]
     f=ts2_JAN[bbalt_CH2_JAN][ind_jr2]
     g=ts3_JAN[bbalt_CH3_JAN][ind_jr3]
     h=ts4_JAN[bbalt_CH4_JAN][ind_jr4]
     l=tb1_JAN[bbalt_CH1_JAN][ind_jr1]
     m=tb2_JAN[bbalt_CH2_JAN][ind_jr2]
     n=tb3_JAN[bbalt_CH3_JAN][ind_jr3]
     o=tb4_JAN[bbalt_CH4_JAN][ind_jr4]
     p=orog1_JAN[bbalt_CH1_JAN][ind_jr1]
     q=orog2_JAN[bbalt_CH2_JAN][ind_jr2]
     r=orog3_JAN[bbalt_CH3_JAN][ind_jr3]
     s=orog4_JAN[bbalt_CH4_JAN][ind_jr4]
     emis1_JAN_moy[ijr]=mean(a[nonzero((a!=-500.)&(a<=1.))])
     emis2_JAN_moy[ijr]=mean(b[nonzero((b!=-500.)&(b<=1.))])
     emis3_JAN_moy[ijr]=mean(c[nonzero((c!=-500.)&(c<=1.))])
     emis4_JAN_moy[ijr]=mean(d[nonzero((d!=-500.)&(d<=1.))])
     ts1_JAN_moy[ijr]=mean(e[nonzero((e!=-500.))])
     ts2_JAN_moy[ijr]=mean(f[nonzero((f!=-500.))])
     ts3_JAN_moy[ijr]=mean(g[nonzero((g!=-500.))])
     ts4_JAN_moy[ijr]=mean(h[nonzero((h!=-500.))])   
     tb1_JAN_moy[ijr]=mean(l[nonzero((l!=-500.))])
     tb2_JAN_moy[ijr]=mean(m[nonzero((m!=-500.))])
     tb3_JAN_moy[ijr]=mean(n[nonzero((n!=-500.))])
     tb4_JAN_moy[ijr]=mean(o[nonzero((o!=-500.))]) 
     orog1_JAN_moy[ijr]=mean(p[nonzero((p!=-500.))])
     orog2_JAN_moy[ijr]=mean(q[nonzero((q!=-500.))])
     orog3_JAN_moy[ijr]=mean(r[nonzero((r!=-500.))])
     orog4_JAN_moy[ijr]=mean(s[nonzero((s!=-500.))])  



i=1 # FEBRUARY
emis_FEB_moy=np.zeros([nbjr[i]],float)
tb_FEB_moy=np.zeros([nbjr[i]],float)
ts_FEB_moy=np.zeros([nbjr[i]],float)
for ijr in range(0,nbjr[i]):
     #print 'jour=', ijr+1
     ind_jr1=np.where(jjr_FEB[bbplat_FEB]==ijr+1)[0] 
#     ind_jr2=np.where(jjr_FEB[bb_tb_FEB]==ijr+1)[0] 
  #   ind_jr3=np.where(jjr_FEB[bb_ts_FEB]==ijr+1)[0] 
     emis_FEB_moy[ijr]=emis_FEB[ind_jr1][nonzero(emis_FEB[ind_jr1]<=1.)].mean()
     tb_FEB_moy[ijr]=tb_FEB[ind_jr1][nonzero(tb_FEB[ind_jr1]!=-500.)].mean()
     ts_FEB_moy[ijr]=ts_FEB[ind_jr1][nonzero(ts_FEB[ind_jr1]!=-500.)].mean()


i=2 # MARCH
emis_MAR_moy=np.zeros([nbjr[i]],float)
tb_MAR_moy=np.zeros([nbjr[i]],float)
ts_MAR_moy=np.zeros([nbjr[i]],float)
for ijr in range(0,nbjr[i]):
     #print 'jour=', ijr+1
     ind_jr1=np.where(jjr_MAR[bb_emis_MAR]==ijr+1)[0] 
#     ind_jr2=np.where(jjr_MAR[bb_tb_MAR]==ijr+1)[0] 
#     ind_jr3=np.where(jjr_MAR[bb_ts_MAR]==ijr+1)[0] 
     emis_MAR_moy[ijr]=emis_MAR[ind_jr1][nonzero(emis_MAR[ind_jr1]<=1.)].mean()
     tb_MAR_moy[ijr]=tb_MAR[ind_jr1][nonzero(tb_MAR[ind_jr1]!=500.)].mean()
     ts_MAR_moy[ijr]=ts_MAR[ind_jr1][nonzero(ts_MAR[ind_jr1]!=-500.)].mean()



fig=plt.figure()
plt.subplot(3,1,1)
plt.plot(arange(1,32,1),emis1_JAN_moy,label='ch1',c='r')
plt.plot(arange(1,32,1),emis2_JAN_moy,label='ch2',c='y')
plt.plot(arange(1,32,1),emis3_JAN_moy,label='ch3',c='g')
plt.plot(arange(1,32,1),emis4_JAN_moy,label='ch4',c='b')
plt.xticks(arange(1,32,5))
grid(True)
ylabel('emissivity')
legend(loc=7)
plt.title('JANUARY 2010')
plt.subplot(3,1,2)
plt.plot(arange(1,32,1),ts1_JAN_moy,label='ch1',c='r')
plt.plot(arange(1,32,1),ts2_JAN_moy,label='ch2',c='y')
plt.plot(arange(1,32,1),ts3_JAN_moy,label='ch3',c='g')
plt.plot(arange(1,32,1),ts4_JAN_moy,label='ch4',c='b')
plt.xticks(arange(1,32,5))
grid(True)
ylabel('surface temp')
legend(loc=7)
plt.subplot(3,1,3)
plt.plot(arange(1,32,1),tb1_JAN_moy,label='ch1',c='r')
plt.plot(arange(1,32,1),tb2_JAN_moy,label='ch2',c='y')
plt.plot(arange(1,32,1),tb3_JAN_moy,label='ch3',c='g')
plt.plot(arange(1,32,1),tb4_JAN_moy,label='ch4',c='b')
plt.xticks(arange(1,32,5))
grid(True)
ylabel('brightness temp')
legend(loc=7)
#plt.subplot(4,1,4)
plt.plot(arange(1,32,1),orog1_JAN_moy,label='ch1',c='r')
plt.plot(arange(1,32,1),orog2_JAN_moy,label='ch2',c='y')
plt.plot(arange(1,32,1),orog3_JAN_moy,label='ch3',c='g')
plt.plot(arange(1,32,1),orog4_JAN_moy,label='ch4',c='b')
plt.xticks(arange(1,32,5))
grid(True)
ylabel('orography')
legend(loc=7)
fig.show()



i=0
plt.subplot(3,1,1)
plt.plot(arange(1,nbjr[i]+1,1),emis_JAN_moy,'o-',c='r',label='emis')
legend()
plt.subplot(3,1,2)
plt.plot(arange(1,nbjr[i]+1,1),tb_JAN_moy,'+-',c='b', label='Tb')
legend()
plt.subplot(3,1,3)
plt.plot(arange(1,nbjr[i]+1,1),ts_JAN_moy,'.-',c='y', label='Ts')
plt.legend()
plt.title('JANUARY')
plt.hold()
plt.show()

i=1
fig=plt.figure()
plt.subplot(3,1,1)
plt.plot(arange(1,nbjr[i]+1,1),emis_FEB_moy,'o-',c='r',label='emis')
legend()
plt.subplot(3,1,2)
plt.plot(arange(1,nbjr[i]+1,1),tb_FEB_moy,'+-',c='b', label='Tb')
legend()
plt.subplot(3,1,3)
plt.plot(arange(1,nbjr[i]+1,1),ts_FEB_moy,'.-',c='y', label='Ts')
plt.legend()
plt.title('FEBRUARY')
fig.show()

e1=np.append(emis_JAN_moy,emis_FEB_moy)
emis_month=np.append(emis_month1,emis_MAR_moy)
ts1=np.append(ts_JAN_moy,ts_FEB_moy)
ts_month=np.append(ts1,ts_MAR_moy)
tb1=np.append(tb_JAN_moy,tb_FEB_moy)
tb_month=np.append(tb1,tb_MAR_moy)

fig=plt.figure()
plt.subplot(3,1,1)
plt.plot(arange(0,len(emis_month)),emis_month,'o-',c='r',label='emissivity')
plt.legend()
plt.grid()
plt.subplot(3,1,2)
plt.plot(arange(0,len(ts_month)),ts_month,'+-',c='g',label='Ts')
plt.legend()
plt.grid()
plt.subplot(3,1,3)
plt.plot(arange(0,len(tb_month)),tb_month,'.-',c='y',label='Tb')
plt.legend()
plt.grid()
fig.show()