source: trunk/src/python_script/read_amsua_ch2_gc.py @ 6

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


fichier=open('/net/dedale/usr/dedale/surf/lelod/ANTARC/AMSUA_CH2_ANTARC_JUNE2010.DAT','r')
numlines = 0
for line in fichier: numlines += 1

fichier.close


fichier=open('/net/dedale/usr/dedale/surf/lelod/ANTARC/AMSUA_CH2_ANTARC_JUNE2010.DAT','r')
nbtotal=numlines

iligne=0
lat_a2=np.zeros([nbtotal],float)
lon_a2=np.zeros([nbtotal],float)
jjr_a2=np.zeros([nbtotal],float)
zen_a2=np.zeros([nbtotal],float)
fov_a2=np.zeros([nbtotal],float)
ts_a2=np.zeros([nbtotal],float)
emis_a2=np.zeros([nbtotal],float)
tb_a2=np.zeros([nbtotal],float)
tup_a2=np.zeros([nbtotal],float)
tdn_a2=np.zeros([nbtotal],float)
trans_a2=np.zeros([nbtotal],float)
orog_a2=np.zeros([nbtotal],float)
mask_a2=np.zeros([nbtotal],float)
pos_a2=np.zeros([nbtotal],float)

while (iligne <= nbtotal-1) :
         line=fichier.readline()
         # exemple : line = "0.22 2.3 5.0 6"
         liste = line.split()
         # exemple : listeCoord ['0.22', '2.3', '5.0', '6'] (liste de chaine de caract?es)
         lon_a2[iligne] = float(liste[0])
         lat_a2[iligne] = float(liste[1])
         jjr_a2[iligne] = float(liste[4])
         pos_a2[iligne] = float(liste[7])
         zen_a2[iligne] = float(liste[8])
         ts_a2[iligne] = float(liste[10])
         orog_a2[iligne] = float(liste[13])
         mask_a2[iligne] = float(liste[14])
         tb_a2[iligne] = float(liste[15])
         emis_a2[iligne] = float(liste[16])
         tdn_a2[iligne] = float(liste[17])
         tup_a2[iligne] = float(liste[18])
         trans_a2[iligne] = float(liste[19])
         iligne=iligne+1
         

fichier.close


vnx=np.array([ 200.,  100.,   60.,   50.,   40.,   30.,   25.,   25.,   20.,20.,   18.,   16.,   15.,   12.,   12.,   12.,   12.,   12.,
12.,   12.,   12.,   12.,   12.,   12.,   12.,   12.,   12.,
12.,   12.,   12.,   12.,   12.,   12.,   12.,   12.,   12.,
12.,   12.,   12.,   12.])

#correr vnx
###################""
#matrices altitude
x=lon_a2
y=lat_a2
z=orog_a2
z0=0
z1=5000

agrid_t, apgrid_t, angrid_t, anngrid_t, asigma_grid_t, xvec, yvec, a_t = newgrid_z.newgrid(x, y, z, vnx,z0, z1)

#matrices ts mask
x=lon_a2
y=lat_a2
z=ts_a2
m=mask_a2
z0=100
z1=400
y11 = -60

import newgrid_zml
tsgrid_ta260, tspgrid_ta260, tsngrid_ta260, tsnngrid_ta260, tssigma_grid_ta260, xvec, yvec, ts_ta260, mm_ta260 = newgrid_zml.newgrid(x, y, z, m, vnx,z0, z1, y11)

#matrices ts mask
x=lon_a2
y=lat_a2
z=ts_a2
m=mask_a2
z0=100
z1=400
y11 = -50

import newgrid_zml
tsgrid_ta2, tspgrid_ta2, tsngrid_ta2, tsnngrid_ta2, tssigma_grid_ta2, xvec, yvec, ts_ta2, mm_ta2 = newgrid_zml.newgrid(x, y, z, m, vnx,z0, z1, y11)


#################################################"
#matrices tb mask
x=lon_a2
y=lat_a2
z=tb_a2
m=mask_a2
z0=100
z1=300
y11 = -60

import newgrid_zml
tbgrid_ta260, tbpgrid_ta260, tbngrid_ta260, tbnngrid_ta260, tbsigma_grid_ta260, xvec60, yvec60, tb_ta260, mm_ta260 = newgrid_zml.newgrid(x, y, z, m, vnx,z0, z1, y11)

#################################################"
#matrices tb mask
x=lon_a2
y=lat_a2
z=tb_a2
m=mask_a2
z0=100
z1=300
y11 = -50

import newgrid_zml
tbgrid_t, tbpgrid_t, tbngrid_t, tbnngrid_t, tbsigma_grid_t, xvec, yvec, tb_t, mm_t = newgrid_zml.newgrid(x, y, z, m, vnx,z0, z1, y11)

nscan = 1
import newgridns_zm
tbgrid_1n, tbpgrid_1n, tbngrid_1n, tbnngrid_1n, tbsigma_grid_1n, xvec, yvec, tb_1n, mm_1n = newgridns_z.newgridns(x, y, z, m, vnx, nscan, pos, z0, z1)

nscan = 2
tbgrid_2n, tbpgrid_2n, tbngrid_2n, tbnngrid_2n, tbsigma_grid_2n, xvec, yvec, tb_2n = newgridns_z.newgridns(x, y, z, vnx, nscan, pos, z0, z1)

nscan = 3
tbgrid_3n, tbpgrid_3n, tbngrid_3n, tbnngrid_3n, tbsigma_grid_3n, xvec, yvec,tb_3n = newgridns_z.newgridns(x, y, z, vnx, nscan, pos, z0, z1)


######################################################################
x=lon_a2
y=lat_a2
z=tup_a2
m=mask_a2
z0=0
z1=70
y11 = -50

import newgrid_zml
tupgrid_ta2, tuppgrid_ta2, tupngrid_ta2, tupnngrid_ta2, tupsigma_grid_ta2, xvec, yvec, tup_ta2, mm_ta2 = newgrid_zml.newgrid(x, y, z, m, vnx,z0, z1, y11)

##################################################################
x=lon_a2
y=lat_a2
z=tup_a2
m=mask_a2
z0=0
z1=70
y11 = -60

import newgrid_zml
tupgrid_ta260, tuppgrid_ta260, tupngrid_ta260, tupnngrid_ta260, tupsigma_grid_ta260, xvec60, yvec60, tup_ta260, mm_ta260 = newgrid_zml.newgrid(x, y, z, m, vnx,z0, z1, y11)

#matrix mask
x=lon_a2
y=lat_a2
z=mask_a2
m=mask_a2
z0=0
z1=1
y11 = -60

import newgrid_zml
mgrid_t60, mpgrid_t60, mngrid_t60, mnngrid_t60, msigma_grid_t60, xvec60, yvec60, m_t60, mm_t60 = newgrid_zml.newgrid(x, y, z, m, vnx,z0, z1, y11)

nscan = 1
import newgridns_zm
tbgrid_1n, tbpgrid_1n, tbngrid_1n, tbnngrid_1n, tbsigma_grid_1n, xvec, yvec, tb_1n, mm_1n = newgridns_z.newgridns(x, y, z, m, vnx, nscan, pos, z0, z1)

nscan = 2
tbgrid_2n, tbpgrid_2n, tbngrid_2n, tbnngrid_2n, tbsigma_grid_2n, xvec, yvec, tb_2n = newgridns_z.newgridns(x, y, z, vnx, nscan, pos, z0, z1)

nscan = 3
tbgrid_3n, tbpgrid_3n, tbngrid_3n, tbnngrid_3n, tbsigma_grid_3n, xvec, yvec,tb_3n = newgridns_z.newgridns(x, y, z, vnx, nscan, pos, z0, z1)

##################################################

#matrices tr
x=lon
y=lat
z=trans
z0=0
z1=1


import newgrid_z
trgrid_t, trpgrid_t, trngrid_t, trnngrid_t, trsigma_grid_t, xvec, yvec, tr_t = newgrid_z.newgrid(x, y, z, vnx,z0, z1)

nscan = 1
import newgridns_z
trgrid_1n, trpgrid_1n, trngrid_1n, trnngrid_1n, trsigma_grid_1n, xvec, yvec, tr_1n = newgridns_z.newgridns(x, y, z, vnx, nscan, pos, z0, z1)

nscan = 2
trgrid_2n, trpgrid_2n, trngrid_2n, trnngrid_2n, trsigma_grid_2n, xvec, yvec, tr_2n = newgridns_z.newgridns(x, y, z, vnx, nscan, pos, z0, z1)

nscan = 3
trgrid_3n, trpgrid_3n, trngrid_3n, trnngrid_3n, trsigma_grid_3n, xvec, yvec,tr_3n = newgridns_z.newgridns(x, y, z, vnx, nscan, pos, z0, z1)

########
#matrices z=emis

z=emis
egrid_t, eepgrid_t, engrid_t, enngrid_t, esigma_grid_t, xvec, yvec = newgrid.newgrid(x, y, z, vnx)

nscan = 1
egrid_1n, eepgrid_1n, engrid_1n, enngrid_1n, esigma_grid_1n, xvec, yvec = newgridns.newgridns(x, y, z, vnx, nscan, pos)

nscan = 2
egrid_2n, eepgrid_2n, engrid_2n, enngrid_2n, esigma_grid_2n, xvec, yvec = newgridns.newgridns(x, y, z, vnx, nscan, pos)

nscan = 3
egrid_3n, eepgrid_3n, engrid_3n, enngrid_3n, esigma_grid_3n, xvec, yvec = newgridns.newgridns(x, y, z, vnx, nscan, pos)

#######################################################
# tb terre

z0=0
z1=1
z=mask

import newgrid_z
mzgrid_t, mzzpgrid_t, mngrid_t, mnngrid_t, msigma_grid_t, xvec, yvec,mzz_t = newgrid_z.newgrid(x, y, z, vnx, z0, z1)

mpgrid=np.round(mzzpgrid_t)
mer=mpgrid
mer[mer==0]=2
mer=mer-1

mertb=tbpgrid_t*mer #tb mer et glace de mer
mertb_1n=tbpgrid_1n*mer #tb mer et glace de mer pos1
mertb_2n=tbpgrid_2n*mer #tb mer et glace de mer pos2
mertb_3n=tbpgrid_3n*mer #tb mer et glace de mer pos3



tbpgrid_tc=tbpgrid_t*mpgrid #tb continent
tbpgrid_1nc=tbpgrid_1n*mpgrid #tb continent pos1
tbpgrid_2nc=tbpgrid_2n*mpgrid #tb continent pos2
tbpgrid_3nc=tbpgrid_3n*mpgrid #tb continent pos3


#################################################"
#vec 
import newvec_zml
x=lon
y=lat
z=tb
m=mask
z0=100
z1=300
y11 = -50
import newvec_x
tb_t50, alt_t50, mm_t50, lo_t50, la_t50 =newvec_x.newgrid(x, y, z, mask, orog, z0, z1, y11)


y11 = -60
import newvecns_zml
nscan = 1
tb_1na, mm_tl1n = newvecns_zml.newgridns(x, y, z, m, vnx, nscan, pos, z0, z1, y11)

nscan = 2
tb_2na, mm_tl2n = newvecns_zml.newgridns(x, y, z, m, vnx, nscan, pos, z0, z1, y11)

nscan = 3
tb_3na, mm_tl3n = newvecns_zml.newgridns(x, y, z, m, vnx, nscan, pos, z0, z1, y11)


##################################
#matrices tb mask
x=lon
y=lat
z=ts
m=mask
z0=100
z1=305
y11 = -50

import newgrid_zml
tsgrid_t, tspgrid_t, tsngrid_t, tsnngrid_t, tssigma_grid_t, xvec, yvec, ts_t, mmts_t = newgrid_zml.newgrid(x, y, z, m, vnx,z0, z1, y11)

nscan = 1
import newgridns_zm
tbgrid_1n, tbpgrid_1n, tbngrid_1n, tbnngrid_1n, tbsigma_grid_1n, xvec, yvec, tb_1n, mm_1n = newgridns_z.newgridns(x, y, z, m, vnx, nscan, pos, z0, z1)

nscan = 2
tbgrid_2n, tbpgrid_2n, tbngrid_2n, tbnngrid_2n, tbsigma_grid_2n, xvec, yvec, tb_2n = newgridns_z.newgridns(x, y, z, vnx, nscan, pos, z0, z1)

nscan = 3
tbgrid_3n, tbpgrid_3n, tbngrid_3n, tbnngrid_3n, tbsigma_grid_3n, xvec, yvec,tb_3n = newgridns_z.newgridns(x, y, z, vnx, nscan, pos, z0, z1)


#####################

dx=0.1
dy=1.0
x0, x1 = -180, 180
y0, y1 = -90, 90

monthly_outz=np.zeros([40,3600],float)
monthly_lon=np.zeros([3600])
monthly_lat=np.zeros([40])
xx = lon
yy = lat
zz = tb
zz0 = 100
zz1= 300
outz, outx, outy, Nil = ffgrid3.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,zz0, zz1)
x=outx
y=outy
z = np.transpose(outz)
z[Nil]=600
z=np.transpose(z)

del outz, outx, outy, zz, xx, yy


# ici je fais des cartes moyennes en melangeant les polars

xx = lon_ssmis
yy = lat_ssmis
zz = 0.5*(emis_ssmis[1,:]+emis_ssmis[2,:])
outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,0.1, 1)
monthly_outz_ssmis_polar[0,:,:]=outz    
del outz, outx, outy, zz

zz = 0.5*(emis_ssmis[4,:]+emis_ssmis[5,:])
outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,0.1, 1)
monthly_outz_ssmis_polar[1,:,:]=outz    
del outz, outx, outy, zz

zz = 0.5*(emis_ssmis[6,:]+emis_ssmis[7,:])
outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,0.1, 1)
monthly_outz_ssmis_polar[2,:,:]=outz    
del outz, outx, outy, zz, xx, yy

# ici je fais des cartes moyennes des differences des polars
xx = lon_ssmis
yy = lat_ssmis
zz = emis_ssmis[1,:]-emis_ssmis[2,:]
outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,-0.05, 0.2)
monthly_outz_ssmis_diff[0,:,:]=outz    
del outz, outx, outy, zz

zz = emis_ssmis[4,:]-emis_ssmis[5,:]
outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,-0.05, 0.2)
monthly_outz_ssmis_diff[1,:,:]=outz    
del outz, outx, outy, zz

zz = emis_ssmis[6,:]-emis_ssmis[7,:]
outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,-0.05, 0.2)
monthly_outz_ssmis_diff[2,:,:]=outz    
del outz, outx, outy, zz, xx, yy

draw_map.draw(monthly_outz_ssmis_polar[0,:], monthly_lon_ssmis, monthly_lat_ssmis, '..\FIG\mean_ssmis_19GHzmpolar_HN_'+le_mois+'.png', '',0.6,1.01,0.002,cm.s3pcpn_l_r)
draw_map.draw(monthly_outz_ssmis_polar[1,:], monthly_lon_ssmis, monthly_lat_ssmis, '..\FIG\mean_ssmis_37GHzmpolar_HN_'+le_mois+'.png', '',0.6,1.01,0.002,cm.s3pcpn_l_r)
draw_map.draw(monthly_outz_ssmis_polar[2,:], monthly_lon_ssmis, monthly_lat_ssmis, '..\FIG\mean_ssmis_85GHzmpolar_HN_'+le_mois+'.png', '',0.6,1.01,0.002,cm.s3pcpn_l_r)

draw_map.draw(monthly_outz_ssmis_diff[0,:], monthly_lon_ssmis, monthly_lat_ssmis, '..\FIG\mean_ssmis_19V-H_HN_'+le_mois+'.png', '',0,0.2,0.002,cm.s3pcpn_l_r)
draw_map.draw(monthly_outz_ssmis_diff[1,:], monthly_lon_ssmis, monthly_lat_ssmis, '..\FIG\mean_ssmis_37V-H_HN_'+le_mois+'.png', '',0,0.2,0.002,cm.s3pcpn_l_r)
draw_map.draw(monthly_outz_ssmis_diff[2,:], monthly_lon_ssmis, monthly_lat_ssmis, '..\FIG\mean_ssmis_85V-H_HN_'+le_mois+'.png', '',0,0.2,0.002,cm.s3pcpn_l_r)

draw_map.draw(monthly_outz_ssmis[0,:], monthly_lon_ssmis, monthly_lat_ssmis, '..\FIG\mean_ssmis_50V_HN_'+le_mois+'.png', '',0.6,1.01,0.002,cm.s3pcpn_l_r)
draw_map.draw(monthly_outz_ssmis[1,:], monthly_lon_ssmis, monthly_lat_ssmis, '..\FIG\mean_ssmis_19V_HN_'+le_mois+'.png', '',0.6,1.01,0.002,cm.s3pcpn_l_r)
draw_map.draw(monthly_outz_ssmis[2,:], monthly_lon_ssmis, monthly_lat_ssmis, '..\FIG\mean_ssmis_19H_HN_'+le_mois+'.png', '',0.6,1.01,0.002,cm.s3pcpn_l_r)
draw_map.draw(monthly_outz_ssmis[3,:], monthly_lon_ssmis, monthly_lat_ssmis, '..\FIG\mean_ssmis_22V_HN_'+le_mois+'.png', '',0.6,1.01,0.002,cm.s3pcpn_l_r)
draw_map.draw(monthly_outz_ssmis[4,:], monthly_lon_ssmis, monthly_lat_ssmis, '..\FIG\mean_ssmis_37V_HN_'+le_mois+'.png', '',0.6,1.01,0.002,cm.s3pcpn_l_r)
draw_map.draw(monthly_outz_ssmis[5,:], monthly_lon_ssmis, monthly_lat_ssmis, '..\FIG\mean_ssmis_37H_HN_'+le_mois+'.png', '',0.6,1.01,0.002,cm.s3pcpn_l_r)
draw_map.draw(monthly_outz_ssmis[6,:], monthly_lon_ssmis, monthly_lat_ssmis, '..\FIG\mean_ssmis_85V_HN_'+le_mois+'.png', '',0.6,1.01,0.002,cm.s3pcpn_l_r)
draw_map.draw(monthly_outz_ssmis[7,:], monthly_lon_ssmis, monthly_lat_ssmis, '..\FIG\mean_ssmis_85H_HN_'+le_mois+'.png', '',0.6,1.01,0.002,cm.s3pcpn_l_r)


bins=arange(0.3,1,0.001)
bb=(lat_ssmis >= 0) 
 
plt.hist(emis_ssmis[0,nonzero(bb)[0]], bins=bins,histtype='step', label='e50V',normed='True',color='#4BB5C1')
plt.hist(emis_ssmis[1,nonzero(bb)[0]], bins=bins,histtype='step', label='e19V',normed='True',color='black')
plt.hist(emis_ssmis[3,nonzero(bb)[0]], bins=bins,histtype='step', label='e22V',normed='True',color='#B9121B')
plt.hist(emis_ssmis[4,nonzero(bb)[0]], bins=bins,histtype='step', label='e37V',normed='True',color='#9748D4')
plt.hist(emis_ssmis[6,nonzero(bb)[0]], bins=bins,histtype='step', label='e91V',normed='True',color='#060DE5')
plt.legend(loc='upper left')
plt.show()
plt.savefig('..\FIG\hist_ssmis_V_NH_'+le_mois+'.png')
close()

bins=arange(0.3,1,0.001)
plt.hist(emis_ssmis[2,nonzero(bb)[0]], bins=bins,histtype='step', label='e19H',normed='True',color='black')
plt.hist(emis_ssmis[5,nonzero(bb)[0]], bins=bins,histtype='step', label='e37H',normed='True',color='#9748D4')
plt.hist(emis_ssmis[7,nonzero(bb)[0]], bins=bins,histtype='step', label='e91H',normed='True',color='#060DE5')
plt.legend(loc='upper left')
plt.show()
plt.savefig('..\FIG\hist_ssmis_H_NH_'+le_mois+'.png')
close()

bins=arange(0.3,1,0.001)
plt.hist(0.5*(emis_ssmis[1,nonzero(bb)[0]]+emis_ssmis[2,nonzero(bb)[0]]), bins=bins,histtype='step', label='e19',normed='True',color='black')
plt.hist(0.5*(emis_ssmis[4,nonzero(bb)[0]]+emis_ssmis[5,nonzero(bb)[0]]), bins=bins,histtype='step', label='e37',normed='True',color='#9748D4')
plt.hist(0.5*(emis_ssmis[6,nonzero(bb)[0]]+emis_ssmis[7,nonzero(bb)[0]]), bins=bins,histtype='step', label='e91',normed='True',color='#060DE5')
plt.legend(loc='upper left')
plt.show()
plt.savefig('..\FIG\hist_ssmis_mpolar_NH_'+le_mois+'.png')
close()


# stats quotidienne autour de la station Thulé
lat_stations=[76.32, 74.43, 78.13, 58.45, 68.6, 64.58]
lon_stations=[-68.3, -94.59, 15.35, -78.08, 33.1, 40.5]
nom_stations=['Thule', 'Resolute', 'Longyearbyen', 'Iqaluit', 'Murmansk', 'Arkhangelsk']


for sta in range(0,6):
    lat0=lat_stations[sta]
    lon0=lon_stations[sta]
    stat_jour=np.zeros([8,7,31],float)
    clear bb
    for canal in range(0,8):
        for jjr in range(0,31):
            jour_obs=jjr+1
            bb=(jjr_ssmis == jour_obs) & (abs(lat_ssmis-lat0) < 2.) & (abs(lon_ssmis-lon0) < 2.)
            stat_jour[canal,0,jjr]=mean(emis_ssmis[canal,nonzero(bb)[0]])
            stat_jour[canal,1,jjr]=std(emis_ssmis[canal,nonzero(bb)[0]])            stat_jour[canal,2,jjr]=size(nonzero(bb))
            stat_jour[canal,3,jjr]=mean(ts_ssmis[nonzero(bb)[0]])
            stat_jour[canal,4,jjr]=std(ts_ssmis[nonzero(bb)[0]])
            stat_jour[canal,5,jjr]=mean(tb_ssmis[canal,nonzero(bb)[0]])
            stat_jour[canal,6,jjr]=std(tb_ssmis[canal,nonzero(bb)[0]])
            del bb
    np.save('STAT_SSMIS_'+nom_stations[sta]+'_'+le_mois+'.dat', stat_jour)
    del stat_jour

mpolar_ssmis=np.zeros([3,nbtotal],float)
mpolar_ssmis[0,:]=0.5*(emis_ssmis[1,:]+emis_ssmis[2,:])
mpolar_ssmis[1,:]=0.5*(emis_ssmis[4,:]+emis_ssmis[5,:])
mpolar_ssmis[2,:]=0.5*(emis_ssmis[6,:]+emis_ssmis[7,:])

mpolarTB_ssmis=np.zeros([3,nbtotal],float)
mpolarTB_ssmis[0,:]=0.5*(tb_ssmis[1,:]+tb_ssmis[2,:])
mpolarTB_ssmis[1,:]=0.5*(tb_ssmis[4,:]+tb_ssmis[5,:])
mpolarTN_ssmis[2,:]=0.5*(tb_ssmis[6,:]+tb_ssmis[7,:])

for sta in range(0,6):
    lat0=lat_stations[sta]
    lon0=lon_stations[sta]
    stat2_jour=np.zeros([3,7,31],float)
    clear bb
    for canal in range(0,3):
        for jjr in range(0,31):
            jour_obs=jjr+1
            bb=(jjr_ssmis == jour_obs) & (abs(lat_ssmis-lat0) < 2.) & (abs(lon_ssmis-lon0) < 2.)
            stat2_jour[canal,0,jjr]=mean(mpolar_ssmis[canal,nonzero(bb)[0]])
            stat2_jour[canal,1,jjr]=std(mpolar_ssmis[canal,nonzero(bb)[0]])
            stat2_jour[canal,2,jjr]=size(nonzero(bb))
            stat2_jour[canal,3,jjr]=mean(ts_ssmis[nonzero(bb)[0]])
            stat2_jour[canal,4,jjr]=std(ts_ssmis[nonzero(bb)[0]])
            stat2_jour[canal,5,jjr]=mean(tb_ssmis[canal,nonzero(bb)[0]])            stat2_jour[canal,6,jjr]=std(tb_ssmis[canal,nonzero(bb)[0]])
            del bb
    np.save('STAT_SSMIS-MPOLAR_'+nom_stations[sta]+'_'+le_mois+'.dat', stat_jour)
    del stat_jour

# ecriture sous format nc
from netCDF4 import Dataset
rootgrp = Dataset('..\EMIS\EMIS_SSMIS_'+le_mois+'.nc', 'w', format='NETCDF4')
 
rootgrp.createDimension('longitude', len(monthly_lon_ssmis))
rootgrp.createDimension('latitude', len(monthly_lat_ssmis))
rootgrp.createDimension('channels', 8)
rootgrp.createDimension('bchannels', 3)

# createVariable (nom de la variable, type, dimensions)
# Si 1 dimension, ne pas oublier la virgule
nclon = rootgrp.createVariable('longitude', 'f8', ('longitude',))
nclat = rootgrp.createVariable('latitude', 'f8', ('latitude',))
ncchan=rootgrp.createVariable('channels', 'f', ('channels',))
ncchan2=rootgrp.createVariable('bchannels', 'f', ('bchannels',))
nctemp = rootgrp.createVariable('emissivity', 'f8', ('channels','latitude', 'longitude'))
nctemp2 = rootgrp.createVariable('emissivity melange polar', 'f8', ('bchannels','latitude', 'longitude'))

nclon[:] = monthly_lon_ssmis
nclat[:] = monthly_lat_ssmis
ncchan[:]=[50,19.1,19.2,22,37.1,37.2,91.1,91.2]
ncchan2[:]=[19,37,91]
nctemp[:] = monthly_outz_ssmis
nctemp2[:] = monthly_outz_ssmis_polar

rootgrp.close()