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

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

modifs Laura

File size: 23.2 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
9from netCDF4 import Dataset
10import ffgrid2
11
12
13
14
15########################
16# EVOLUTION TEMPORELLE #
17########################
18
19## ZONE1 = Autour de "Dome C" (lon=123.23;lat=-75.06) - mois: JUNE ##
20## ch17 ##
21#bbzone_CH17_JUN = nonzero((lon17_JUN >= 120.23) & (lon17_JUN <= 126.23) & (lat17_JUN >= -78.06) & (lat17_JUN <= -72.06))
22
23## ch18 ##
24#bbzone_CH18_JUN = nonzero((lon18_JUN >= 120.23) & (lon18_JUN <= 126.23) & (lat18_JUN >= -78.06) & (lat18_JUN <= -72.06))
25
26
27## ZONE2 = Autre zone de glace de mer / de continent (lon entre -40 et -60 // lat entre -75 et -85) ##
28## ch17 ##
29bbzone_CH17_FEB = nonzero((lon17_FEB >= -60.) & (lon17_FEB <= -40.) & (lat17_FEB >= -85.) & (lat17_FEB <= -75.))
30bbzone_CH17_APR = nonzero((lon17_APR >= -60.) & (lon17_APR <= -40.) & (lat17_APR >= -85.) & (lat17_APR <= -75.))
31bbzone_CH17_MAY = nonzero((lon17_MAY >= -60.) & (lon17_MAY <= -40.) & (lat17_MAY >= -85.) & (lat17_MAY <= -75.))
32bbzone_CH17_JUN = nonzero((lon17_JUN >= -60.) & (lon17_JUN <= -40.) & (lat17_JUN >= -85.) & (lat17_JUN <= -75.))
33bbzone_CH17_JUL = nonzero((lon17_JUL >= -60.) & (lon17_JUL <= -40.) & (lat17_JUL >= -85.) & (lat17_JUL <= -75.))
34
35## ch18 ##
36#bbzone_CH18_JUN = nonzero((lon18_JUN >= -60.) & (lon18_JUN <= -40.) & (lat18_JUN >= -85.) & (lat18_JUN <= -75.))
37
38## ch 16 ##
39#bbzone_CH16_JUN = nonzero((lon16_JUN >= 120.23) & (lon16_JUN <= 126.23) & (lat16_JUN >= -78.06) & (lat16_JUN <= -72.06))
40
41
42## FEBUARY ##
43emis17_FEB_moy = np.zeros([28],float)
44tb17_FEB_moy = np.zeros([28],float)
45ts17_FEB_moy = np.zeros([28],float)
46orog17_FEB_moy = np.zeros([28],float)
47for ijr in range (0,28):
48     print 'jour=', ijr+1
49     ## ch17 ##
50     ind_jr17_FEB = np.where(jjr17_FEB[bbzone_CH17_FEB]==ijr+1)[0]
51     a = emis17_FEB[bbzone_CH17_FEB][ind_jr17_FEB]
52     b = tb17_FEB[bbzone_CH17_FEB][ind_jr17_FEB]
53     c = ts17_FEB[bbzone_CH17_FEB][ind_jr17_FEB]
54     d = orog17_FEB[bbzone_CH17_FEB][ind_jr17_FEB]
55     emis17_FEB_moy[ijr] = mean(a[nonzero((a!=-500.)&(a<=1.))])
56     tb17_FEB_moy[ijr] = mean(b[nonzero((b!=-500.)&(b!=0.))])
57     ts17_FEB_moy[ijr] = mean(c[nonzero((c!=-500.)&(c!=0.))])
58     orog17_FEB_moy[ijr] = mean(d[nonzero((d!=-500.))])
59
60
61## APRIL ##
62emis17_APR_moy = np.zeros([30],float)
63tb17_APR_moy = np.zeros([30],float)
64ts17_APR_moy = np.zeros([30],float)
65orog17_APR_moy = np.zeros([30],float)
66for ijr in range (0,30):
67     print 'jour=', ijr+1
68     ind_jr17_APR = np.where(jjr17_APR[bbzone_CH17_APR]==ijr+1)[0]
69     a = emis17_APR[bbzone_CH17_APR][ind_jr17_APR]
70     b = tb17_APR[bbzone_CH17_APR][ind_jr17_APR]
71     c = ts17_APR[bbzone_CH17_APR][ind_jr17_APR]
72     d = orog17_APR[bbzone_CH17_APR][ind_jr17_APR]
73     emis17_APR_moy[ijr] = mean(a[nonzero((a!=-500.)&(a<=1.))])
74     tb17_APR_moy[ijr] = mean(b[nonzero((b!=-500.)&(b!=0.))])
75     ts17_APR_moy[ijr] = mean(c[nonzero((c!=-500.)&(c!=0.))])
76     orog17_APR_moy[ijr] = mean(d[nonzero((d!=-500.))])
77
78
79## MAY ##
80emis17_MAY_moy = np.zeros([31],float)
81tb17_MAY_moy = np.zeros([31],float)
82ts17_MAY_moy = np.zeros([31],float)
83orog17_MAY_moy = np.zeros([31],float)
84for ijr in range (0,31):
85     print 'jour=', ijr+1
86     ind_jr17_MAY = np.where(jjr17_MAY[bbzone_CH17_MAY]==ijr+1)[0]
87     a = emis17_MAY[bbzone_CH17_MAY][ind_jr17_MAY]
88     b = tb17_MAY[bbzone_CH17_MAY][ind_jr17_MAY]
89     c = ts17_MAY[bbzone_CH17_MAY][ind_jr17_MAY]
90     d = orog17_MAY[bbzone_CH17_MAY][ind_jr17_MAY]
91     emis17_MAY_moy[ijr] = mean(a[nonzero((a!=-500.)&(a<=1.))])
92     tb17_MAY_moy[ijr] = mean(b[nonzero((b!=-500.)&(b!=0.))])
93     ts17_MAY_moy[ijr] = mean(c[nonzero((c!=-500.)&(c!=0.))])
94     orog17_MAY_moy[ijr] = mean(d[nonzero((d!=-500.))])
95
96
97## JUNE ##
98emis17_JUN_moy = np.zeros([30],float)
99tb17_JUN_moy = np.zeros([30],float)
100ts17_JUN_moy = np.zeros([30],float)
101#orog17_JUN_moy = np.zeros([30],float)
102#emis18_JUN_moy = np.zeros([30],float)
103#tb18_JUN_moy = np.zeros([30],float)
104#ts18_JUN_moy = np.zeros([30],float)
105#emis16_JUN_moy = np.zeros([30],float)
106#tb16_JUN_moy = np.zeros([30],float)
107#ts16_JUN_moy = np.zeros([30],float)
108for ijr in range (0,30):
109     print 'jour=', ijr+1
110     ind_jr17_JUN = np.where(jjr17_JUN[bbzone_CH17_JUN]==ijr+1)[0]
111     a = emis17_JUN[bbzone_CH17_JUN][ind_jr17_JUN]
112     b = tb17_JUN[bbzone_CH17_JUN][ind_jr17_JUN]
113     c = ts17_JUN[bbzone_CH17_JUN][ind_jr17_JUN]
114     d = orog17_JUN[bbzone_CH17_JUN][ind_jr17_JUN]
115     emis17_JUN_moy[ijr] = mean(a[nonzero((a!=-500.)&(a<=1.))])
116     tb17_JUN_moy[ijr] = mean(b[nonzero((b!=-500.)&(b!=0.))])
117     ts17_JUN_moy[ijr] = mean(c[nonzero((c!=-500.)&(c!=0.))])
118#     orog17_JUN_moy[ijr] = mean(d[nonzero((d!=-500.))])
119     ## ch18 ##
120#     ind_jr18 = np.where(jjr18_JUN[bbzone_CH18_JUN]==ijr+1)[0]
121#     d = emis18_JUN[bbzone_CH18_JUN][ind_jr18]
122#     e = tb18_JUN[bbzone_CH18_JUN][ind_jr18]
123#     f = ts18_JUN[bbzone_CH18_JUN][ind_jr18]
124#     emis18_JUN_moy[ijr] = mean(d[nonzero((d != -500.) & (d <= 1.))])
125#     tb18_JUN_moy[ijr] = mean(e[nonzero((e != -500.) & (e != 0.))])
126#     ts18_JUN_moy[ijr] = mean(f[nonzero((f != -500.) & (f != 0.))])
127     ## ch16 ##
128#     ind_jr16 = np.where(jjr16_JUN[bbzone_CH16_JUN]==ijr+1)[0]
129#     g = emis16_JUN[bbzone_CH16_JUN][ind_jr16]
130#     h = tb16_JUN[bbzone_CH16_JUN][ind_jr16]
131#     l = ts16_JUN[bbzone_CH16_JUN][ind_jr16]
132#     emis16_JUN_moy[ijr] = mean(g[nonzero((g != -500.) & (g <= 1.))])
133#     tb16_JUN_moy[ijr] = mean(h[nonzero((h != -500.) & (h != 0.))])
134#     ts16_JUN_moy[ijr] = mean(l[nonzero((l != -500.) & (l != 0.))])
135
136
137## JULY ##
138emis17_JUL_moy = np.zeros([31],float)
139tb17_JUL_moy = np.zeros([31],float)
140ts17_JUL_moy = np.zeros([31],float)
141orog17_JUL_moy = np.zeros([31],float)
142for ijr in range (0,31):
143     print 'jour=', ijr+1
144     ind_jr17_JUL = np.where(jjr17_JUL[bbzone_CH17_JUL]==ijr+1)[0]
145     a = emis17_JUL[bbzone_CH17_JUL][ind_jr17_JUL]
146     b = tb17_JUL[bbzone_CH17_JUL][ind_jr17_JUL]
147     c = ts17_JUL[bbzone_CH17_JUL][ind_jr17_JUL]
148     d = orog17_JUL[bbzone_CH17_JUL][ind_jr17_JUL]
149     emis17_JUL_moy[ijr] = mean(a[nonzero((a!=-500.)&(a<=1.))])
150     tb17_JUL_moy[ijr] = mean(b[nonzero((b!=-500.)&(b!=0.))])
151     ts17_JUL_moy[ijr] = mean(c[nonzero((c!=-500.)&(c!=0.))])
152     orog17_JUL_moy[ijr] = mean(d[nonzero((d!=-500.))])
153
154
155## plot evolution temporelle ##
156## FEBUARY ##
157fig = plt.figure()
158plt.subplot(3, 1, 1)
159plt.title('FEBUARY')
160plt.plot(arange(0, 28, 1), emis17_FEB_moy, c = 'b')
161plt.ylabel('Emissivity')
162xticks(arange(0, 28, 1), arange(1, 29, 1))
163xlim (0, 28)
164grid(True)
165plt.subplot(3, 1, 2)
166plt.plot(arange(0, 28, 1), tb17_FEB_moy, c = 'g', label = 'Tb')
167ylabel('Brightness temperature')
168xticks(arange(0, 28, 1), arange(1, 29, 1))
169xlim (0, 28)
170grid(True)
171plt.subplot(3, 1, 3)
172plt.plot(arange(0, 28, 1), ts17_FEB_moy, c = 'r', label = 'Ts')
173plt.ylabel('Surface temperature')
174xticks(arange(0, 28, 1), arange(1, 29, 1))
175xlim (0, 28)
176grid(True)
177fig.show()
178
179## from APRIL to JULY ##
180emis1 = np.append(emis17_APR_moy, emis17_MAY_moy)
181emis2 = np.append(emis1, emis17_JUN_moy)
182emis17_moy = np.append(emis2, emis17_JUL_moy)
183tb1 = np.append(tb17_APR_moy, tb17_MAY_moy)
184tb2 = np.append(tb1, tb17_JUN_moy)
185tb17_moy = np.append(tb2, tb17_JUL_moy)
186ts1 = np.append(ts17_APR_moy, ts17_MAY_moy)
187ts2 = np.append(ts1, ts17_JUN_moy)
188ts17_moy = np.append(ts2, ts17_JUL_moy)
189fig = plt.figure()
190plt.subplot(3, 1, 1)
191plt.plot(arange(0, 122, 1), emis17_moy, c = 'b')
192plt.ylabel('Emissivity')
193xticks(array([0, 30, 61, 91, 122]), date[1:])
194xlim (0, 122)
195grid(True)
196plt.subplot(3, 1, 2)
197plt.plot(arange(0, 122, 1), tb17_moy, c = 'g', label = 'Tb')
198ylabel('Brightness temperature')
199xticks(array([0, 30, 61, 91, 122]), date[1:])
200xlim (0, 122)
201grid(True)
202plt.subplot(3, 1, 3)
203plt.plot(arange(0, 122, 1), ts17_moy, c = 'r', label = 'Ts')
204plt.ylabel('Surface temperature')
205xticks(array([0, 30, 61, 91, 122]), date[1:])
206xlim (0, 122)
207grid(True)
208fig.show()
209
210## JUNE ##
211## ch17 - ch18 ##
212fig=plt.figure()
213plt.subplot(3,1,1)
214plt.plot(arange(0,30,1),emis17_JUN_moy,label='91.66GHz-V',c='r')
215plt.plot(arange(0,30,1),emis18_JUN_moy,label='91.66GHz-H',c='b')
216plt.xticks(arange(1,31,1))
217grid(True)
218ylabel('emissivity')
219legend(loc=7)
220plt.title('SSMIS - JUNE 2010')
221plt.subplot(3,1,2)
222plt.plot(arange(0,30,1),tb17_JUN_moy,label='91.66GHz-V',c='r')
223plt.plot(arange(0,30,1),tb18_JUN_moy,label='91.66GHz-H',c='b')
224plt.xticks(arange(1,31,1))
225grid(True)
226ylabel('brightness temperature')
227legend(loc=7)
228plt.subplot(3,1,3)
229plt.plot(arange(0,30,1),ts17_JUN_moy,label='91.66GHz-V',c='r')
230plt.plot(arange(0,30,1),ts18_JUN_moy,label='91.66GHz-H',c='b')
231plt.xticks(arange(1,31,1))
232grid(True)
233ylabel('surface temperature')
234legend(loc=7)
235fig.show()
236## ch16 - ch17 ##
237fig=plt.figure()
238plt.subplot(3,1,1)
239plt.plot(arange(0,30,1),emis16_JUN_moy,label='37GHz-V',c='r')
240plt.plot(arange(0,30,1),emis17_JUN_moy,label='91.66GHz-V',c='g')
241plt.xticks(arange(1,31,1))
242grid(True)
243ylabel('emissivity')
244legend(loc=7)
245plt.title('SSMIS - JUNE 2010')
246plt.subplot(3,1,2)
247plt.plot(arange(0,30,1),tb16_JUN_moy,label='37GHz-V',c='r')
248plt.plot(arange(0,30,1),tb17_JUN_moy,label='91.66GHz-V',c='g')
249plt.xticks(arange(1,31,1))
250grid(True)
251ylabel('brightness temperature')
252legend(loc=7)
253plt.subplot(3,1,3)
254plt.plot(arange(0,30,1),ts16_JUN_moy,label='37GHz-V',c='r')
255plt.plot(arange(0,30,1),ts17_JUN_moy,label='91.66GHz-V',c='g')
256plt.xticks(arange(1,31,1))
257grid(True)
258ylabel('surface temperature')
259legend(loc=7)
260fig.show()
261
262
263## calcul anomalie de Tb
264## FEBUARY ##
265tb17_FEB_anom = np.zeros([28], float)
266for ijr in range (0,28):
267     print 'jour=', ijr + 1
268     tb17_FEB_anom[ijr]=tb17_FEB_moy[ijr]-mean(tb17_FEB_moy)
269
270## APRIL ##
271tb17_APR_anom = np.zeros([30], float)
272for ijr in range (0,30):
273     print 'jour=', ijr + 1
274     tb17_APR_anom[ijr]=tb17_APR_moy[ijr]-mean(tb17_APR_moy)
275
276## MAY ##
277bbnan = nonzero(isnan(tb17_MAY_moy) == False)
278tb17_MAY_anom = np.zeros([31], float)
279for ijr in range (0,31):
280     print 'jour=', ijr + 1
281     tb17_MAY_anom[ijr]=tb17_MAY_moy[ijr]-mean(tb17_MAY_moy[bbnan])
282
283## JUNE ##
284tb17_JUN_anom = np.zeros([30],float)
285#tb18_JUN_anom = np.zeros([30],float)
286for ijr in range (0,30):
287     print 'jour=', ijr + 1
288     tb17_JUN_anom[ijr]=tb17_JUN_moy[ijr]-mean(tb17_JUN_moy)
289#     tb18_JUN_anom[ijr]=tb18_JUN_moy[ijr]-mean(tb18_JUN_moy)
290
291## JULY ##
292tb17_JUL_anom = np.zeros([31],float)
293for ijr in range (0,30):
294     print 'jour=', ijr + 1
295     tb17_JUL_anom[ijr]=tb17_JUL_moy[ijr]-mean(tb17_JUL_moy)
296
297
298
299a1 = np.append(tb17_APR_anom, tb17_MAY_anom)
300a2 = np.append(a1, tb17_JUN_anom)
301tb17_month_anom = np.append(a2, tb17_JUL_anom)
302
303
304## for APRIL to JULY ##
305fig = plt.figure()
306plt.plot(arange(0, 122, 1), tb17_month_anom, c = 'k')
307plt.plot(arange(0, 122, 1), np.zeros([122]), '--', c = 'r')
308ylabel('Tb anomaly (K)')
309xlim(0, 122)
310xticks(array([0, 30, 61, 91 122]), date[1:])
311grid(True)
312plt.title('SSMIS - zone2')
313fig.show()
314
315## FEBRUARY ##
316fig = plt.figure()
317plt.plot(arange(0, 28, 1), tb17_FEB_anom, c = 'k')
318plt.plot(arange(0, 28, 1), np.zeros([28]), '--', c = 'r')
319ylabel('Tb anomaly (K)')
320xlim(0, 27)
321xticks(arange(0, 28, 1), arange(1, 29, 1))
322grid(True)
323plt.xlabel('FABRUARY')
324plt.title('SSMIS - zone2')
325fig.show()
326
327## JUNE ##
328fig=plt.figure()
329plt.plot(arange(0,30,1),tb17_JUN_anom,label='91.66Ghz-V',c='r')
330plt.plot(arange(0,30,1),tb18_JUN_anom,label='91.66GHz-H',c='b')
331plt.plot(arange(0,31,1),np.zeros([31]),'--',c='k')
332ylabel('Tb anomaly')
333legend()
334twinx().plot(arange(0,30,1),ts17_JUN_moy,label='surface temperature',c='g')
335ylabel('surface temperature')
336plt.xticks(arange(0, 30, 1), arange(1, 31, 1))
337grid(True)
338plt.title('SSMIS - JUNE 2010')
339fig.show()
340
341
342
343
344################
345# CARTOGRAPHIE #
346################
347## Etude sur l'Antarctique ##
348dx=5.
349dy=5.
350x0, x1 = -180, 180
351y0, y1 = -90, -30
352
353
354## FEBRUARY ##
355bbtb_ch17_FEB = nonzero((tb17_FEB != -500.) & (tb17_FEB != 0.))
356OUTZCH17_FEB = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx)),28], float)
357outzch17_FEB = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx))], float)
358lonch17_FEB = np.zeros([len(np.arange(x0, x1+1, dx))], float)
359latch17_FEB = np.zeros([len(np.arange(y0, y1+1, dy))], float)
360for ijr in range(0,28):
361    print 'jour=', ijr+1
362    ## ch17 ##
363    ind_jr17_FEB = np.where(jjr17_FEB[bbtb_ch17_FEB] == ijr+1)[0]
364    xx = lon17_FEB[bbtb_ch17_FEB][ind_jr17_FEB]
365    yy = lat17_FEB[bbtb_ch17_FEB][ind_jr17_FEB]
366    zz = tb17_FEB[bbtb_ch17_FEB][ind_jr17_FEB]
367    zz0 = min(zz)
368    zz1 = max(zz)
369    outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0, x1, y0, y1, zz0, zz1)
370    outzch17_FEB = outz
371    lonch17_FEB = outx
372    latch17_FEB = outy
373    OUTZCH17_FEB[:,:,ijr] = outzch17_FEB[:,:]
374
375## APRIL ##
376bbtb_ch17_APR = nonzero((tb17_APR != -500.) & (tb17_APR != 0.))
377OUTZCH17_APR = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx)),30], float)
378outzch17_APR = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx))], float)
379lonch17_APR = np.zeros([len(np.arange(x0, x1+1, dx))], float)
380latch17_APR = np.zeros([len(np.arange(y0, y1+1, dy))], float)
381for ijr in range(0,30):
382    print 'jour=', ijr+1
383    ## ch17 ##
384    ind_jr17_APR = np.where(jjr17_APR[bbtb_ch17_APR] == ijr+1)[0]
385    xx = lon17_APR[bbtb_ch17_APR][ind_jr17_APR]
386    yy = lat17_APR[bbtb_ch17_APR][ind_jr17_APR]
387    zz = tb17_APR[bbtb_ch17_APR][ind_jr17_APR]
388    zz0 = min(zz)
389    zz1 = max(zz)
390    outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0, x1, y0, y1, zz0, zz1)
391    outzch17_APR = outz
392    lonch17_APR = outx
393    latch17_APR = outy
394    OUTZCH17_APR[:,:,ijr] = outzch17_APR[:,:]
395
396## MAY ##
397bbtb_ch17_MAY = nonzero((tb17_MAY != -500.) & (tb17_MAY != 0.))
398OUTZCH17_MAY = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx)),31], float)
399outzch17_MAY = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx))], float)
400lonch17_MAY = np.zeros([len(np.arange(x0, x1+1, dx))], float)
401latch17_MAY = np.zeros([len(np.arange(y0, y1+1, dy))], float)
402for ijr in range(0,31):
403    print 'jour=', ijr+1
404    ## ch17 ##
405    ind_jr17_MAY = np.where(jjr17_MAY[bbtb_ch17_MAY] == ijr+1)[0]
406    xx = lon17_MAY[bbtb_ch17_MAY][ind_jr17_MAY]
407    yy = lat17_MAY[bbtb_ch17_MAY][ind_jr17_MAY]
408    zz = tb17_MAY[bbtb_ch17_MAY][ind_jr17_MAY]
409    zz0 = min(zz)
410    zz1 = max(zz)
411    outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0, x1, y0, y1, zz0, zz1)
412    outzch17_MAY = outz
413    lonch17_MAY = outx
414    latch17_MAY = outy
415    OUTZCH17_MAY[:,:,ijr] = outzch17_MAY[:,:]
416
417## JUNE ##
418## ch17 ##
419bbtb_ch17_JUN = nonzero((tb17_JUN != -500.) & (tb17_JUN != 0.))
420OUTZCH17_JUN = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx)),30], float)
421outzch17_JUN = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx))], float)
422lonch17_JUN = np.zeros([len(np.arange(x0, x1+1, dx))], float)
423latch17_JUN = np.zeros([len(np.arange(y0, y1+1, dy))], float)
424## ch18 ##
425bbemis_ch18_JUN=nonzero((emis18_JUN!=-500.)&(emis18_JUN<1.)&(emis18_JUN>0.))
426bbtb_ch18_JUN=nonzero((tb18_JUN!=-500.)&(tb18_JUN!=0.))
427OUTZCH18_JUN = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx)),30], float)
428outzch18_JUN = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx))], float)
429lonch18_JUN=np.zeros([len(np.arange(x0, x1+1, dx))], float)
430latch18_JUN=np.zeros([len(np.arange(y0, y1+1, dy))], float)
431for ijr in range(0,30):
432    print 'jour=', ijr+1
433    ## ch17 ##
434    ind_jr17_JUN = np.where(jjr17_JUN[bbtb_ch17_JUN] == ijr+1)[0]
435    xx = lon17_JUN[bbtb_ch17_JUN][ind_jr17_JUN]
436    yy = lat17_JUN[bbtb_ch17_JUN][ind_jr17_JUN]
437    zz = tb17_JUN[bbtb_ch17_JUN][ind_jr17_JUN]
438    zz0 = min(zz)
439    zz1 = max(zz)
440    outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0, x1, y0, y1, zz0, zz1)
441    outzch17_JUN = outz
442    lonch17_JUN = outx
443    latch17_JUN = outy
444    OUTZCH17_JUN[:,:,ijr] = outzch17_JUN[:,:]
445    ## ch18 ##
446    ind_jr18_JUN = np.where(jjr18_JUN[bbtb_ch18_JUN] == ijr+1)[0]
447    xx = lon18_JUN[bbtb_ch18_JUN][ind_jr18_JUN]
448    yy = lat18_JUN[bbtb_ch18_JUN][ind_jr18_JUN]
449    zz = tb18_JUN[bbtb_ch18_JUN][ind_jr18_JUN]
450    zz0 = min(zz)
451    zz1 = max(zz)
452    outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0, x1, y0, y1, zz0, zz1)
453    outzch18_JUN = outz
454    lonch18_JUN = outx
455    latch18_JUN = outy
456    OUTZCH18_JUN[:,:,ijr] = outzch18_JUN[:,:]
457
458## JULY ##
459bbtb_ch17_JUL = nonzero((tb17_JUL != -500.) & (tb17_JUL != 0.))
460OUTZCH17_JUL = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx)),31], float)
461outzch17_JUL = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx))], float)
462lonch17_JUL = np.zeros([len(np.arange(x0, x1+1, dx))], float)
463latch17_JUL = np.zeros([len(np.arange(y0, y1+1, dy))], float)
464for ijr in range(0,31):
465    print 'jour=', ijr+1
466    ## ch17 ##
467    ind_jr17_JUL = np.where(jjr17_JUL[bbtb_ch17_JUL] == ijr+1)[0]
468    xx = lon17_JUL[bbtb_ch17_JUL][ind_jr17_JUL]
469    yy = lat17_JUL[bbtb_ch17_JUL][ind_jr17_JUL]
470    zz = tb17_JUL[bbtb_ch17_JUL][ind_jr17_JUL]
471    zz0 = min(zz)
472    zz1 = max(zz)
473    outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0, x1, y0, y1, zz0, zz1)
474    outzch17_JUL = outz
475    lonch17_JUL = outx
476    latch17_JUL = outy
477    OUTZCH17_JUL[:,:,ijr] = outzch17_JUL[:,:]
478
479
480## calcul de la climatologie moyenne sur le mois et anomalie en chaque lon/lat ##
481## FEBRUARY ##
482mean_outzch17_FEB = np.zeros([len(latch17_FEB), len(lonch17_FEB)], float)
483for ilon in range (0,len(lonch17_FEB)):
484    for ilat in range (0,len(latch17_FEB)):
485        mean_outzch17_FEB[ilat,ilon] = mean(OUTZCH17_FEB[ilat,ilon,:])
486
487tbch17_anom_FEB = np.zeros([len(latch17_FEB),len(lonch17_FEB),28], float)
488for ijr in range (0,28):
489    tbch17_anom_FEB[:,:,ijr] = OUTZCH17_FEB[:,:,ijr] - mean_outzch17_FEB[:,:]
490
491## APRIL ##
492mean_outzch17_APR = np.zeros([len(latch17_APR), len(lonch17_APR)], float)
493for ilon in range (0,len(lonch17_APR)):
494    for ilat in range (0,len(latch17_APR)):
495        mean_outzch17_APR[ilat,ilon] = mean(OUTZCH17_APR[ilat,ilon,:])
496
497tbch17_anom_APR = np.zeros([len(latch17_APR),len(lonch17_APR),30], float)
498for ijr in range (0,30):
499    tbch17_anom_APR[:,:,ijr] = OUTZCH17_APR[:,:,ijr] - mean_outzch17_APR[:,:]
500
501## MAY ##
502mean_outzch17_MAY = np.zeros([len(latch17_MAY), len(lonch17_MAY)], float)
503for ilon in range (0,len(lonch17_MAY)):
504    for ilat in range (0,len(latch17_MAY)):
505        mean_outzch17_MAY[ilat,ilon] = mean(OUTZCH17_MAY[ilat,ilon,:])
506
507tbch17_anom_MAY = np.zeros([len(latch17_MAY),len(lonch17_MAY),31], float)
508for ijr in range (0,31):
509    tbch17_anom_MAY[:,:,ijr] = OUTZCH17_MAY[:,:,ijr] - mean_outzch17_MAY[:,:]
510
511
512## JUNE ##
513## ch17 ##
514mean_outzch17_JUN = np.zeros([len(latch17_JUN), len(lonch17_JUN)], float)
515for ilon in range (0,len(lonch17_JUN)):
516    for ilat in range (0,len(latch17_JUN)):
517        mean_outzch17_JUN[ilat,ilon] = OUTZCH17_JUN[ilat,ilon,:].mean()
518
519tbch17_anom_JUN = np.zeros([len(latch17_JUN),len(lonch17_JUN),30], float)
520for ijr in range (0,30):
521    tbch17_anom_JUN[:,:,ijr] = OUTZCH17_JUN[:,:,ijr] - mean_outzch17_JUN[:,:]
522
523## ch18 ##
524mean_outzch18_JUN = np.zeros([len(latch18_JUN), len(lonch18_JUN)], float)
525for ilon in range (0,len(lonch18_JUN)):
526    for ilat in range (0,len(latch18_JUN)):
527        mean_outzch18_JUN[ilat,ilon] = OUTZCH18_JUN[ilat,ilon,:].mean()
528
529tbch18_anom_JUN = np.zeros([len(latch18_JUN),len(lonch18_JUN),30], float)
530for ijr in range (0,30):
531    tbch18_anom_JUN[:,:,ijr] = OUTZCH18_JUN[:,:,ijr] - mean_outzch18_JUN[:,:]
532
533
534
535for ijr in range (23, 30):
536     figure()
537     plt.ion()
538     m = Basemap(llcrnrlon=-60, urcrnrlon=-40, llcrnrlat=-85, urcrnrlat=-75, projection='cyl', resolution='c', fix_aspect=True)
539     m.drawcoastlines(linewidth = 1)
540     m.drawparallels(np.arange(-85., 75., 2))
541     m.drawmeridians(np.arange(-60., -40., 2))
542     #m.fillcontinents()
543     clevs = arange(-25., 5., 0.1)
544     xii,yii = m(*np.meshgrid(lonch17_JUN, latch17_JUN))
545     cs = m.contourf(xii, yii, tbch17_anom_JUN[:,:,ijr], clevs, cmap=cm.s3pcpn_l_r)
546     cbar = colorbar(cs)
547     cbar.set_label('Tb anomaly SSMIS CH17 - JUNE')
548     xticks(np.arange(-60., -40., 2))
549     yticks(np.arange(-85., 75., 2))
550     plt.savefig('/usr/home/lahlod/twice_d/figures_output_ANTARC/SSMIS/mean_tb_anomaly_map_'+str(ijr+1)+'JUN_ch17_SSMIS_Zoom_zone2')
551
552
553figure()
554plt.ion()
555m = Basemap(llcrnrlon=-180, urcrnrlon=180, llcrnrlat=-90, urcrnrlat=-30, projection='cyl', resolution='c', fix_aspect=True)
556m.drawcoastlines(linewidth = 1)
557m.drawparallels(np.arange(-90., -30., 20))
558m.drawmeridians(np.arange(-180., 180., 20))
559#m.fillcontinents()
560xii,yii = m(*np.meshgrid(lonch17_FEB, latch17_FEB))
561plt.xticks(arange(-180, 200, 20))
562plt.yticks(arange(-90, -30, 20))
563## ch17 ##
564clevs = arange(150., 270., 1.)
565cs = m.contourf(xii, yii, mean_outzch17_FEB, clevs, cmap=cm.s3pcpn_l_r)
566cbar = colorbar(cs)
567cbar.set_label('Mean Tb FEB [CH17] - SSMIS')
568## BIAIS ch17-ch18 ##
569biais_JUN = mean_outzch17_JUN - mean_outzch18_JUN
570clevs = arange(0., 45., 1.)
571cs = m.contourf(xii, yii, biais_JUN, clevs, cmap=cm.s3pcpn_l_r)
572cbar = colorbar(cs)
573cbar.set_label('Bias [CH17-CH18] of Mean Tb JUNE - SSMIS')
574
575## VARIANCE ch17 - ch18 ##
576std1 = np.zeros([len(latch17_JUN), len(lonch17_JUN)], float)
577for ilat in range (0, len(latch17_JUN)):
578    for ilon in range (0, len(lonch17_JUN)):
579        std1[ilat, ilon] = (mean_outzch17_JUN[ilat, ilon]**2)-(mean_outzch18_JUN[ilat, ilon]**2)
580
581
582N = len(lonch17_JUN)*len(latch17_JUN)
583std_JUN = std1/N
584clevs = arange(0., 25., 0.1)
585cs = m.contourf(xii, yii, std_JUN, clevs, cmap=cm.s3pcpn_l_r)
586cbar = colorbar(cs)
587cbar.set_label('Stantard deviation [CH17-CH18] of Mean Tb JUNE - SSMIS')
588
589
590
591##########################
592# DIAGRAMME DE HOVMOLLER #
593##########################
594# shape(tbch17_anom_JUN) = [ilat, ilon, ijr]
595## FEBRUARY ##
596bbtranche17_FEB = nonzero((latch17_FEB >= -85.) & (latch17_FEB <= -75))
597mean_tbch17_anom_FEB = np.zeros([len(lonch17_FEB), 28], float)
598for ilon in range (0, len(lonch17_FEB)):
599    for ijr in range (0,28):
600        mean_tbch17_anom_FEB[ilon,ijr] = mean(tbch17_anom_FEB[bbtranche17_FEB][:,ilon,ijr])
601
602y_time, x_space = np.meshgrid(arange(0,28,1), lonch17_FEB)
603fig = plt.figure()
604plt.pcolor(x_space, y_time, mean_tbch17_anom_FEB, cmap=cm.s3pcpn_l_r, vmin = -10., vmax = 15.)
605plt.axis([-180., 180., 0, 28])
606cb = plt.colorbar()
607cb.set_label('Tb anomaly - SSMIS CH17')
608plt.xticks(arange(-180.,200.,40))
609plt.yticks(arange(0, 28, 1), arange (1, 29, 1))
610plt.xlabel('longitude')
611plt.ylabel('FEBRUARY 2010')
612
613## MAY ##
614bbtranche17_MAY = nonzero((latch17_MAY >= -85.) & (latch17_MAY <= -75))
615mean_tbch17_anom_MAY = np.zeros([len(lonch17_MAY), 31], float)
616for ilon in range (0, len(lonch17_MAY)):
617    for ijr in range (0,31):
618        mean_tbch17_anom_MAY[ilon,ijr] = mean(tbch17_anom_MAY[bbtranche17_MAY][:,ilon,ijr])
619
620y_time, x_space = np.meshgrid(arange(0,31,1), lonch17_MAY)
621fig = plt.figure()
622plt.pcolor(x_space, y_time, mean_tbch17_anom_MAY, cmap=cm.s3pcpn_l_r, vmin = -12., vmax = 25.)
623plt.axis([-180., 180., 0, 30])
624cb = plt.colorbar()
625cb.set_label('Tb anomaly - SSMIS CH17')
626plt.xticks(arange(-180.,200.,40))
627plt.yticks(arange(0, 30, 1), arange(1,31,1))
628plt.xlabel('longitude')
629plt.ylabel('MAY 2010')
630
631## JUNE ##
632bbtranche17_JUN = nonzero((latch17_JUN >= -85.) & (latch17_JUN <= -75))
633mean_tbch17_anom_JUN = np.zeros([len(lonch17_JUN), 30], float)
634for ilon in range (0, len(lonch17_JUN)):
635    for ijr in range (0,30):
636        mean_tbch17_anom_JUN[ilon,ijr] = tbch17_anom_JUN[bbtranche17_JUN][:,ilon,ijr].mean()
637
638
639y_time, x_space = np.meshgrid(arange(0,30,1), lonch17_APR)
640fig = plt.figure()
641plt.pcolor(x_space, y_time, mean_tbch17_anom_APR, cmap=cm.s3pcpn_l_r, vmin = -30., vmax = 25.)
642plt.axis([-180., 180., 0, 29])
643cb = plt.colorbar()
644cb.set_label('Tb anomaly - SSMIS CH17')
645plt.xticks(arange(-180.,220.,40))
646plt.yticks(arange(0, 30, 1), arange(1,31,1))
647plt.xlabel('longitude')
648plt.ylabel('JUNE 2010')
649
Note: See TracBrowser for help on using the repository browser.