source: trunk/src/scripts_Laura/read_SSMIS_test.py @ 20

Last change on this file since 20 was 20, checked in by lahlod, 10 years ago

modif_script_SSMIS

File size: 14.5 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3import string
4import numpy as np
5import matplotlib.pyplot as plt
6from pylab import *
7from mpl_toolkits.basemap import Basemap
8from mpl_toolkits.basemap import shiftgrid, cm
9import netCDF4
10import ffgrid2
11
12
13################ fichiers par mois - CH1 ###################################################
14
15f1 = '/net/dedale/usr/dedale/surf/lelod/ANTARC/SSMIS_CH'
16f3 = '_ANTARC_JUNE2010.DAT'
17#date=np.array(['JANUARY', 'FEBRUARY', 'MARCH', 'APRIL', 'MAY', 'JUNE', 'JULY'])
18channel = np.array([12, 13, 15, 16, 17, 18])
19numlines = np.zeros([len(channel)],int)
20
21for ich in range (0, len(channel)):
22     print channel[ich] 
23     f = f1 + str(channel[ich]) + f3
24     fichier = open(f, 'r')
25     numlines[ich] = 0
26     for line in fichier: numlines[ich] += 1
27
28
29     fichier.close()
30
31ich = 2 # 37GHz, H polar
32fichier = open(f1 + str(channel[ich]) + f3, 'r')
33ssmis = np.zeros([18, numlines[ich]], float)
34for iligne in range (0,numlines[ich]):
35    line = fichier.readline()
36    liste = line.split()
37    for j in range(0,18):
38        ssmis[j,iligne] = float(liste[j])
39
40
41    fichier.close
42
43
44ssch15_JUN=ssmis
45lon15_JUN=ssch15_JUN[0,:]
46lat15_JUN=ssch15_JUN[1,:]
47jjr15_JUN=ssch15_JUN[4,:]
48ts15_JUN=ssch15_JUN[8,:]
49emis15_JUN=ssch15_JUN[14,:]
50tb15_JUN=ssch15_JUN[13,:]
51tup15_JUN=ssch15_JUN[16,:]
52tdn15_JUN=ssch15_JUN[15,:]
53trans15_JUN=ssch15_JUN[17,:]
54orog15_JUN=ssch15_JUN[11,:]
55
56
57ich = 3 # 37GHz, V polar
58fichier = open(f1 + str(channel[ich]) + f3, 'r')
59ssmis = np.zeros([18, numlines[ich]], float)
60for iligne in range (0,numlines[ich]):
61    line = fichier.readline()
62    liste = line.split()
63    for j in range(0,18):
64        ssmis[j,iligne] = float(liste[j])
65
66
67    fichier.close
68
69
70ssch16_JUN=ssmis
71lon16_JUN=ssch16_JUN[0,:]
72lat16_JUN=ssch16_JUN[1,:]
73jjr16_JUN=ssch16_JUN[4,:]
74ts16_JUN=ssch16_JUN[8,:]
75emis16_JUN=ssch16_JUN[14,:]
76tb16_JUN=ssch16_JUN[13,:]
77tup16_JUN=ssch16_JUN[16,:]
78tdn16_JUN=ssch16_JUN[15,:]
79trans16_JUN=ssch16_JUN[17,:]
80orog16_JUN=ssch16_JUN[11,:]
81
82
83ich = 0 # 19.35GHz, H polar
84fichier = open(f1 + str(channel[ich]) + f3, 'r')
85ssmis = np.zeros([18, numlines[ich]], float)
86for iligne in range (0,numlines[ich]-1):
87    line = fichier.readline()
88    liste = line.split()
89    for j in range(0,18):
90        ssmis[j,iligne] = float(liste[j])
91
92
93    fichier.close
94
95
96ssch12_JUN=ssmis
97lon12_JUN=ssch12_JUN[0,:]
98lat12_JUN=ssch12_JUN[1,:]
99jjr12_JUN=ssch12_JUN[4,:]
100ts12_JUN=ssch12_JUN[8,:]
101emis12_JUN=ssch12_JUN[14,:]
102tb12_JUN=ssch12_JUN[13,:]
103tup12_JUN=ssch12_JUN[16,:]
104tdn12_JUN=ssch12_JUN[15,:]
105trans12_JUN=ssch12_JUN[17,:]
106orog12_JUN=ssch12_JUN[11,:]
107
108
109ich = 1 # 19.35GHz, V polar
110fichier = open(f1 + str(channel[ich]) + f3, 'r')
111ssmis = np.zeros([18, numlines[ich]], float)
112for iligne in range (0,numlines[ich]-1):
113    line = fichier.readline()
114    liste = line.split()
115    for j in range(0,18):
116        ssmis[j,iligne] = float(liste[j])
117
118
119    fichier.close
120
121
122ssch13_JUN=ssmis
123lon13_JUN=ssch13_JUN[0,:]
124lat13_JUN=ssch13_JUN[1,:]
125jjr13_JUN=ssch13_JUN[4,:]
126ts13_JUN=ssch13_JUN[8,:]
127emis13_JUN=ssch13_JUN[14,:]
128tb13_JUN=ssch13_JUN[13,:]
129tup13_JUN=ssch13_JUN[16,:]
130tdn13_JUN=ssch13_JUN[15,:]
131trans13_JUN=ssch13_JUN[17,:]
132orog13_JUN=ssch13_JUN[11,:]
133
134
135ich = 4 # 91.66GHz, V polar
136fichier = open(f1 + str(channel[ich]) + f3, 'r')
137ssmis = np.zeros([18, numlines[ich]], float)
138for iligne in range (0,numlines[ich]-1):
139    line = fichier.readline()
140    liste = line.split()
141    for j in range(0,18):
142        ssmis[j,iligne] = float(liste[j])
143
144
145    fichier.close
146
147
148ssch17_JUN=ssmis
149lon17_JUN=ssch17_JUN[0,:]
150lat17_JUN=ssch17_JUN[1,:]
151jjr17_JUN=ssch17_JUN[4,:]
152ts17_JUN=ssch17_JUN[8,:]
153emis17_JUN=ssch17_JUN[14,:]
154tb17_JUN=ssch17_JUN[13,:]
155tup17_JUN=ssch17_JUN[16,:]
156tdn17_JUN=ssch17_JUN[15,:]
157trans17_JUN=ssch17_JUN[17,:]
158orog17_JUN=ssch17_JUN[11,:]
159
160
161ich = 5 # 91.66GHz, V polar
162fichier = open(f1 + str(channel[ich]) + f3, 'r')
163ssmis = np.zeros([18, numlines[ich]], float)
164for iligne in range (0,numlines[ich]-1):
165    line = fichier.readline()
166    liste = line.split()
167    for j in range(0,18):
168        ssmis[j,iligne] = float(liste[j])
169
170
171    fichier.close
172
173
174ssch18_JUN=ssmis
175lon18_JUN=ssch18_JUN[0,:]
176lat18_JUN=ssch18_JUN[1,:]
177jjr18_JUN=ssch18_JUN[4,:]
178ts18_JUN=ssch18_JUN[8,:]
179emis18_JUN=ssch18_JUN[14,:]
180tb18_JUN=ssch18_JUN[13,:]
181tup18_JUN=ssch18_JUN[16,:]
182tdn18_JUN=ssch18_JUN[15,:]
183trans18_JUN=ssch18_JUN[17,:]
184orog18_JUN=ssch18_JUN[11,:]
185
186########################
187# EVOLUTION TEMPORELLE #
188########################
189## Autour de "Dome C" (lon=123.23;lat=-75.06) - mois: JUNE ##
190## Autre zone de glace de mer / de continent (lon entre -40 et -60 // lat entre -75 et -85)##
191## ch17 ##
192#bbzone_CH17_JUN = nonzero((lon17_JUN>=120.23)&(lon17_JUN<=126.23)&(lat17_JUN>=-78.06)&(lat17_JUN<=-72.06))
193bbzone_CH17_JUN = nonzero((lon17_JUN >= -60.) & (lon17_JUN <= -40.) & (lat17_JUN >= -85.) & (lat17_JUN <= -75.))
194emis17_JUN_moy = np.zeros([30],float)
195tb17_JUN_moy = np.zeros([30],float)
196ts17_JUN_moy = np.zeros([30],float)
197
198## ch18 ##
199#bbzone_CH18_JUN = nonzero((lon18_JUN >= 120.23) & (lon18_JUN <= 126.23) & (lat18_JUN >= -78.06) & (lat18_JUN <= -72.06))
200bbzone_CH18_JUN = nonzero((lon18_JUN >= -60.) & (lon18_JUN <= -40.) & (lat18_JUN >= -85.) & (lat18_JUN <= -75.))
201emis18_JUN_moy = np.zeros([30],float)
202tb18_JUN_moy = np.zeros([30],float)
203ts18_JUN_moy = np.zeros([30],float)
204
205## ch 16 ##
206#bbzone_CH16_JUN = nonzero((lon16_JUN >= 120.23) & (lon16_JUN <= 126.23) & (lat16_JUN >= -78.06) & (lat16_JUN <= -72.06))
207#emis16_JUN_moy = np.zeros([30],float)
208#tb16_JUN_moy = np.zeros([30],float)
209#ts16_JUN_moy = np.zeros([30],float)
210
211for ijr in range (0,30):
212     print 'jour=', ijr+1
213     ## ch17 ##
214     ind_jr17 = np.where(jjr17_JUN[bbzone_CH17_JUN]==ijr+1)[0]
215     a = emis17_JUN[bbzone_CH17_JUN][ind_jr17]
216     b = tb17_JUN[bbzone_CH17_JUN][ind_jr17]
217     c = ts17_JUN[bbzone_CH17_JUN][ind_jr17]
218     emis17_JUN_moy[ijr] = mean(a[nonzero((a!=-500.)&(a<=1.))])
219     tb17_JUN_moy[ijr] = mean(b[nonzero((b!=-500.)&(b!=0.))])
220     ts17_JUN_moy[ijr] = mean(c[nonzero((c!=-500.)&(c!=0.))])
221     ## ch18 ##
222     ind_jr18 = np.where(jjr18_JUN[bbzone_CH18_JUN]==ijr+1)[0]
223     d = emis18_JUN[bbzone_CH18_JUN][ind_jr18]
224     e = tb18_JUN[bbzone_CH18_JUN][ind_jr18]
225     f = ts18_JUN[bbzone_CH18_JUN][ind_jr18]
226     emis18_JUN_moy[ijr] = mean(d[nonzero((d != -500.) & (d <= 1.))])
227     tb18_JUN_moy[ijr] = mean(e[nonzero((e != -500.) & (e != 0.))])
228     ts18_JUN_moy[ijr] = mean(f[nonzero((f != -500.) & (f != 0.))])
229     ## ch16 ##
230#     ind_jr16 = np.where(jjr16_JUN[bbzone_CH16_JUN]==ijr+1)[0]
231#     g = emis16_JUN[bbzone_CH16_JUN][ind_jr16]
232#     h = tb16_JUN[bbzone_CH16_JUN][ind_jr16]
233#     l = ts16_JUN[bbzone_CH16_JUN][ind_jr16]
234#     emis16_JUN_moy[ijr] = mean(g[nonzero((g != -500.) & (g <= 1.))])
235#     tb16_JUN_moy[ijr] = mean(h[nonzero((h != -500.) & (h != 0.))])
236#     ts16_JUN_moy[ijr] = mean(l[nonzero((l != -500.) & (l != 0.))])
237
238
239## ch17 - ch18 ##
240fig=plt.figure()
241plt.subplot(3,1,1)
242plt.plot(arange(0,30,1),emis17_JUN_moy,label='91.66GHz-V',c='r')
243plt.plot(arange(0,30,1),emis18_JUN_moy,label='91.66GHz-H',c='b')
244plt.xticks(arange(1,31,1))
245grid(True)
246ylabel('emissivity')
247legend(loc=7)
248plt.title('SSMIS - JUNE 2010')
249plt.subplot(3,1,2)
250plt.plot(arange(0,30,1),tb17_JUN_moy,label='91.66GHz-V',c='r')
251plt.plot(arange(0,30,1),tb18_JUN_moy,label='91.66GHz-H',c='b')
252plt.xticks(arange(1,31,1))
253grid(True)
254ylabel('brightness temperature')
255legend(loc=7)
256plt.subplot(3,1,3)
257plt.plot(arange(0,30,1),ts17_JUN_moy,label='91.66GHz-V',c='r')
258plt.plot(arange(0,30,1),ts18_JUN_moy,label='91.66GHz-H',c='b')
259plt.xticks(arange(1,31,1))
260grid(True)
261ylabel('surface temperature')
262legend(loc=7)
263fig.show()
264
265## ch16 ##
266fig=plt.figure()
267plt.subplot(3,1,1)
268plt.plot(arange(0,30,1),emis16_JUN_moy,label='37GHz-V',c='r')
269plt.plot(arange(0,30,1),emis17_JUN_moy,label='91.66GHz-V',c='g')
270plt.xticks(arange(1,31,1))
271grid(True)
272ylabel('emissivity')
273legend(loc=7)
274plt.title('SSMIS - JUNE 2010')
275plt.subplot(3,1,2)
276plt.plot(arange(0,30,1),tb16_JUN_moy,label='37GHz-V',c='r')
277plt.plot(arange(0,30,1),tb17_JUN_moy,label='91.66GHz-V',c='g')
278plt.xticks(arange(1,31,1))
279grid(True)
280ylabel('brightness temperature')
281legend(loc=7)
282plt.subplot(3,1,3)
283plt.plot(arange(0,30,1),ts16_JUN_moy,label='37GHz-V',c='r')
284plt.plot(arange(0,30,1),ts17_JUN_moy,label='91.66GHz-V',c='g')
285plt.xticks(arange(1,31,1))
286grid(True)
287ylabel('surface temperature')
288legend(loc=7)
289fig.show()
290
291
292
293tb17_JUN_anom = np.zeros([30],float)
294tb18_JUN_anom = np.zeros([30],float)
295for ijr in range (0,30):
296     print 'jour=', ijr + 1
297     tb17_JUN_anom[ijr]=tb17_JUN_moy[ijr]-tb17_JUN_moy.mean()
298     tb18_JUN_anom[ijr]=tb18_JUN_moy[ijr]-tb18_JUN_moy.mean()
299
300
301fig=plt.figure()
302plt.plot(arange(0,30,1),tb17_JUN_anom,label='91.66Ghz-V',c='r')
303plt.plot(arange(0,30,1),tb18_JUN_anom,label='91.66GHz-H',c='b')
304plt.plot(arange(0,31,1),np.zeros([31]),'--',c='k')
305ylabel('Tb anomaly')
306legend()
307twinx().plot(arange(0,30,1),ts17_JUN_moy,label='surface temperature',c='g')
308ylabel('surface temperature')
309plt.xticks(arange(0, 30, 1), arange(1, 31, 1))
310grid(True)
311plt.title('SSMIS - JUNE 2010')
312fig.show()
313
314
315
316
317################
318# CARTOGRAPHIE #
319################
320## Etude sur l'Antarctique ##
321dx=5.
322dy=5.
323x0, x1 = -180, 180
324y0, y1 = -90, -30
325
326## ch17 ##
327bbemis_ch17_JUN=nonzero((emis17_JUN!=-500.)&(emis17_JUN<1.)&(emis17_JUN>0.))
328bbtb_ch17_JUN=nonzero((tb17_JUN!=-500.)&(tb17_JUN!=0.))
329OUTZCH17_JUN = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx)),30], float)
330outzch17_JUN = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx))], float)
331lonch17_JUN=np.zeros([len(np.arange(x0, x1+1, dx))], float)
332latch17_JUN=np.zeros([len(np.arange(y0, y1+1, dy))], float)
333
334## ch18 ##
335bbemis_ch18_JUN=nonzero((emis18_JUN!=-500.)&(emis18_JUN<1.)&(emis18_JUN>0.))
336bbtb_ch18_JUN=nonzero((tb18_JUN!=-500.)&(tb18_JUN!=0.))
337OUTZCH18_JUN = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx)),30], float)
338outzch18_JUN = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx))], float)
339lonch18_JUN=np.zeros([len(np.arange(x0, x1+1, dx))], float)
340latch18_JUN=np.zeros([len(np.arange(y0, y1+1, dy))], float)
341
342for ijr in range(0,30):
343    print 'jour=', ijr+1
344    ## ch17 ##
345    ind_jr17_JUN = np.where(jjr17_JUN[bbtb_ch17_JUN] == ijr+1)[0]
346    xx = lon17_JUN[bbtb_ch17_JUN][ind_jr17_JUN]
347    yy = lat17_JUN[bbtb_ch17_JUN][ind_jr17_JUN]
348    zz = tb17_JUN[bbtb_ch17_JUN][ind_jr17_JUN]
349    zz0 = min(zz)
350    zz1 = max(zz)
351    outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0, x1, y0, y1, zz0, zz1)
352    outzch17_JUN = outz
353    lonch17_JUN = outx
354    latch17_JUN = outy
355    OUTZCH17_JUN[:,:,ijr] = outzch17_JUN[:,:]
356    ## ch18 ##
357    ind_jr18_JUN = np.where(jjr18_JUN[bbtb_ch18_JUN] == ijr+1)[0]
358    xx = lon18_JUN[bbtb_ch18_JUN][ind_jr18_JUN]
359    yy = lat18_JUN[bbtb_ch18_JUN][ind_jr18_JUN]
360    zz = tb18_JUN[bbtb_ch18_JUN][ind_jr18_JUN]
361    zz0 = min(zz)
362    zz1 = max(zz)
363    outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0, x1, y0, y1, zz0, zz1)
364    outzch18_JUN = outz
365    lonch18_JUN = outx
366    latch18_JUN = outy
367    OUTZCH18_JUN[:,:,ijr] = outzch18_JUN[:,:]
368
369   
370## ch17 ##
371mean_outzch17_JUN = np.zeros([len(latch17_JUN), len(lonch17_JUN)], float)
372for ilon in range (0,len(lonch17_JUN)):
373    for ilat in range (0,len(latch17_JUN)):
374        mean_outzch17_JUN[ilat,ilon] = OUTZCH17_JUN[ilat,ilon,:].mean()
375
376
377## ch18 ##
378mean_outzch18_JUN = np.zeros([len(latch18_JUN), len(lonch18_JUN)], float)
379for ilon in range (0,len(lonch18_JUN)):
380    for ilat in range (0,len(latch18_JUN)):
381        mean_outzch18_JUN[ilat,ilon] = OUTZCH18_JUN[ilat,ilon,:].mean()
382
383
384## ch17 ##
385tbch17_anom_JUN = np.zeros([len(latch17_JUN),len(lonch17_JUN),30], float)
386for ijr in range (0,30):
387    tbch17_anom_JUN[:,:,ijr] = OUTZCH17_JUN[:,:,ijr] - mean_outzch17_JUN[:,:]
388
389
390## ch18 ##
391tbch18_anom_JUN = np.zeros([len(latch18_JUN),len(lonch18_JUN),30], float)
392for ijr in range (0,30):
393    tbch18_anom_JUN[:,:,ijr] = OUTZCH18_JUN[:,:,ijr] - mean_outzch18_JUN[:,:]
394
395
396for ijr in range (23, 30):
397     figure()
398     plt.ion()
399     m = Basemap(llcrnrlon=-60, urcrnrlon=-40, llcrnrlat=-85, urcrnrlat=-75, projection='cyl', resolution='c', fix_aspect=True)
400     m.drawcoastlines(linewidth = 1)
401     m.drawparallels(np.arange(-85., 75., 2))
402     m.drawmeridians(np.arange(-60., -40., 2))
403     #m.fillcontinents()
404     clevs = arange(-25., 5., 0.1)
405     xii,yii = m(*np.meshgrid(lonch17_JUN, latch17_JUN))
406     cs = m.contourf(xii, yii, tbch17_anom_JUN[:,:,ijr], clevs, cmap=cm.s3pcpn_l_r)
407     cbar = colorbar(cs)
408     cbar.set_label('Tb anomaly SSMIS CH17 - JUNE')
409     xticks(np.arange(-60., -40., 2))
410     yticks(np.arange(-85., 75., 2))
411     plt.savefig('/usr/home/lahlod/twice_d/figures_output_ANTARC/SSMIS/mean_tb_anomaly_map_'+str(ijr+1)+'JUN_ch17_SSMIS_Zoom_zone2')
412
413
414figure()
415plt.ion()
416m = Basemap(llcrnrlon=-180, urcrnrlon=180, llcrnrlat=-90, urcrnrlat=-30, projection='cyl', resolution='c', fix_aspect=True)
417m.drawcoastlines(linewidth = 1)
418m.drawparallels(np.arange(-90., 90., 20))
419m.drawmeridians(np.arange(-180., 180., 20))
420#m.fillcontinents()
421xii,yii = m(*np.meshgrid(lonch17_JUN, latch17_JUN))
422
423## BIAIS ch17-ch18 ##
424biais_JUN = mean_outzch17_JUN - mean_outzch18_JUN
425clevs = arange(0., 45., 1.)
426cs = m.contourf(xii, yii, biais_JUN, clevs, cmap=cm.s3pcpn_l_r)
427cbar = colorbar(cs)
428cbar.set_label('Bias [CH17-CH18] of Mean Tb JUNE - SSMIS')
429
430## VARIANCE ch17 - ch18 ##
431std1 = np.zeros([len(latch17_JUN), len(lonch17_JUN)], float)
432for ilat in range (0, len(latch17_JUN)):
433    for ilon in range (0, len(lonch17_JUN)):
434        std1[ilat, ilon] = (mean_outzch17_JUN[ilat, ilon]**2)-(mean_outzch18_JUN[ilat, ilon]**2)
435
436
437N = len(lonch17_JUN)*len(latch17_JUN)
438std_JUN = std1/N
439clevs = arange(0., 25., 0.1)
440cs = m.contourf(xii, yii, std_JUN, clevs, cmap=cm.s3pcpn_l_r)
441cbar = colorbar(cs)
442cbar.set_label('Stantard deviation [CH17-CH18] of Mean Tb JUNE - SSMIS')
443
444
445
446##########################
447# DIAGRAMME DE HOVMOLLER #
448##########################
449
450# shape(tbch17_anom_JUN) = [ilat, ilon, ijr]
451
452## tranche de latitude étudiée ##
453#bbtranche17_JUN = nonzero((latch17_JUN == -75.))
454bbtranche17_JUN = nonzero((latch17_JUN >= -85.) & (latch17_JUN <= -75))
455#bbtranche17_JUN = nonzero((latch17_JUN >= -90.) & (latch17_JUN <= -85))
456mean_tbch17_anom_JUN = np.zeros([len(lonch17_JUN), 30], float)
457for ilon in range (0, len(lonch17_JUN)):
458    for ijr in range (0,30):
459        mean_tbch17_anom_JUN[ilon,ijr] = tbch17_anom_JUN[bbtranche17_JUN][:,ilon,ijr].mean()
460
461
462y_time, x_space = m(*np.meshgrid(arange(0,30,1), lonch17_JUN))
463fig = plt.figure()
464plt.pcolor(x_space, y_time, mean_tbch17_anom_JUN, cmap=cm.s3pcpn_l_r, vmin = -30., vmax = 25.)
465plt.axis([-180., 180., 0, 29])
466cb = plt.colorbar()
467cb.set_label('Tb anomaly - SSMIS CH17')
468plt.xticks(arange(-180.,220.,40))
469plt.yticks(arange(0, 30, 1), arange(1,31,1))
470plt.xlabel('longitude')
471plt.ylabel('JUNE 2010')
472
Note: See TracBrowser for help on using the repository browser.