1 | #!/usr/bin/env python |
---|
2 | # -*- coding: utf-8 -*- |
---|
3 | import string |
---|
4 | import numpy as np |
---|
5 | import matplotlib.pyplot as plt |
---|
6 | from pylab import * |
---|
7 | from mpl_toolkits.basemap import Basemap |
---|
8 | from mpl_toolkits.basemap import shiftgrid, cm |
---|
9 | from netCDF4 import Dataset |
---|
10 | import ffgrid2 |
---|
11 | |
---|
12 | |
---|
13 | |
---|
14 | |
---|
15 | ## ATTENTION: lecture des fichiers NETCDF AMSUA - definition des variables dans "read_netcdf_AMSUA_test.py" ## |
---|
16 | ## lon |
---|
17 | ## lat |
---|
18 | ## jrs_MONTH |
---|
19 | ## chan |
---|
20 | len_month = np.array([31, 28, 31, 30, 31, 30, 31]) |
---|
21 | ######################### |
---|
22 | ## moyennes mensuelles ## |
---|
23 | ######################### |
---|
24 | monthly_emis_JAN = np.zeros([len(chan), len(lat), len(lon)], float) |
---|
25 | monthly_emis_FEB = np.zeros([len(chan), len(lat), len(lon)], float) |
---|
26 | monthly_emis_MAR = np.zeros([len(chan), len(lat), len(lon)], float) |
---|
27 | monthly_emis_APR = np.zeros([len(chan), len(lat), len(lon)], float) |
---|
28 | monthly_emis_MAY = np.zeros([len(chan), len(lat), len(lon)], float) |
---|
29 | monthly_emis_JUN = np.zeros([len(chan), len(lat), len(lon)], float) |
---|
30 | monthly_emis_JUL = np.zeros([len(chan), len(lat), len(lon)], float) |
---|
31 | for ich in range (0, len(chan)): |
---|
32 | for ilon in range (0, len(lon)): |
---|
33 | for ilat in range (0, len(lat)): |
---|
34 | monthly_emis_JAN[ich, ilat, ilon] = mean(emis_JAN[ich, ilat, ilon, :][nonzero(isnan(emis_JAN[ich, ilat, ilon, :]) == False)]) |
---|
35 | monthly_emis_FEB[ich, ilat, ilon] = mean(emis_FEB[ich, ilat, ilon, :][nonzero(isnan(emis_FEB[ich, ilat, ilon, :]) == False)]) |
---|
36 | monthly_emis_MAR[ich, ilat, ilon] = mean(emis_MAR[ich, ilat, ilon, :][nonzero(isnan(emis_MAR[ich, ilat, ilon, :]) == False)]) |
---|
37 | monthly_emis_APR[ich, ilat, ilon] = mean(emis_APR[ich, ilat, ilon, :][nonzero(isnan(emis_APR[ich, ilat, ilon, :]) == False)]) |
---|
38 | monthly_emis_MAY[ich, ilat, ilon] = mean(emis_MAY[ich, ilat, ilon, :][nonzero(isnan(emis_MAY[ich, ilat, ilon, :]) == False)]) |
---|
39 | monthly_emis_JUN[ich, ilat, ilon] = mean(emis_JUN[ich, ilat, ilon, :][nonzero(isnan(emis_JUN[ich, ilat, ilon, :]) == False)]) |
---|
40 | monthly_emis_JUL[ich, ilat, ilon] = mean(emis_JUL[ich, ilat, ilon, :][nonzero(isnan(emis_JUL[ich, ilat, ilon, :]) == False)]) |
---|
41 | |
---|
42 | monthly_emis = [monthly_emis_JAN, monthly_emis_FEB, monthly_emis_MAR, monthly_emis_APR, monthly_emis_MAY, monthly_emis_JUN, monthly_emis_JUL] |
---|
43 | |
---|
44 | ############################## |
---|
45 | ## plot moyennes mensuelles ## |
---|
46 | ############################## |
---|
47 | ## evolution de l'emissivite moyenne mensuelle dans une tranche de latitude ## |
---|
48 | plt.ion() |
---|
49 | for ilat in range (2, 7): |
---|
50 | plt.figure() |
---|
51 | plt.plot(lon, monthly_emis_JAN[0, ilat, :], 'r--', label = 'JAN') |
---|
52 | plt.plot(lon, monthly_emis_FEB[0, ilat, :], 'g--', label = 'FEB') |
---|
53 | plt.plot(lon, monthly_emis_MAR[0, ilat, :], 'b--', label = 'MAR') |
---|
54 | plt.plot(lon, monthly_emis_APR[0, ilat, :], 'm', label = 'APR') |
---|
55 | plt.plot(lon, monthly_emis_MAY[0, ilat, :], 'y', label = 'MAY') |
---|
56 | plt.plot(lon, monthly_emis_JUN[0, ilat, :], 'k', label = 'JUN') |
---|
57 | plt.plot(lon, monthly_emis_JUL[0, ilat, :], 'c', label = 'JUL') |
---|
58 | plt.xticks(np.arange(-180, 200, 20), rotation = 45) |
---|
59 | legend(bbox_to_anchor = (0., 1.02, 1., 0.01), loc = 3, ncol = 7) |
---|
60 | plt.ylabel('emissivity') |
---|
61 | plt.xlabel('longitude') |
---|
62 | plt.title('latitude: ' + str(lat[ilat]) + ' - CH1 - AMSUA') |
---|
63 | grid(True) |
---|
64 | plt.savefig('/usr/home/lahlod/twice_d/figures_output_ANTARC/AMSUA/emis_lon_lat' + str(lat[ilat]) + '_plt_FEB-APR-JUL_ch1_AMSUA.png') |
---|
65 | |
---|
66 | |
---|
67 | ## evolution biais dans une tranche de latitude ## |
---|
68 | plt.figure() |
---|
69 | lat_color = np.array(['r', 'm', 'b', 'g', 'c']) |
---|
70 | for ilat in range (2, 7): |
---|
71 | plt.plot(lon, (monthly_emis_JAN[0, ilat, :] / monthly_emis_JUL[0, ilat, :]), lat_color[ilat-2], label = str(lat[ilat])) |
---|
72 | |
---|
73 | legend(bbox_to_anchor = (0., 1.04), loc = 3, ncol = 5) |
---|
74 | plt.xticks(np.arange(-180, 200, 20), rotation = 45) |
---|
75 | plt.ylabel('emissivity bias') |
---|
76 | plt.xlabel('longitude') |
---|
77 | grid(True) |
---|
78 | plt.title('emiss(JAN)/emiss(JUL)') |
---|
79 | plt.plot(lon, np.ones([len(lon)]), 'k--') |
---|
80 | |
---|
81 | |
---|
82 | ################ |
---|
83 | ## emissivite ## |
---|
84 | ################ |
---|
85 | Njr = np.sum(len_month) |
---|
86 | emis_glob = np.zeros([len(chan), len(lat), len(lon), Njr], float) |
---|
87 | for ich in range (0, len(chan)): |
---|
88 | for ilon in range (0, len(lon)): |
---|
89 | for ilat in range (0, len(lat)): |
---|
90 | e1 = np.append(emis_JAN[ich, ilat, ilon, :], emis_FEB[ich, ilat, ilon, :]) |
---|
91 | e2 = np.append(e1, emis_MAR[ich, ilat, ilon, :]) |
---|
92 | e3 = np.append(e2, emis_APR[ich, ilat, ilon, :]) |
---|
93 | e4 = np.append(e3, emis_MAY[ich, ilat, ilon, :]) |
---|
94 | e5 = np.append(e4, emis_JUN[ich, ilat, ilon, :]) |
---|
95 | e6 = np.append(e5, emis_JUL[ich, ilat, ilon, :]) |
---|
96 | emis_glob[ich, ilat, ilon, :] = e6 |
---|
97 | |
---|
98 | y_time, x_space = np.meshgrid(arange(0, Njr, 1), lon) |
---|
99 | for ilat in range (2, 7): |
---|
100 | plt.figure() |
---|
101 | plt.title('emissivity - AMSUA - latitude ' + str(lat[ilat])) |
---|
102 | for ich in range (0, len(chan)): |
---|
103 | plt.subplot(2, 1, ich+1) |
---|
104 | plt.pcolor(x_space, y_time, emis_glob[ich, ilat, :, :], cmap=cm.s3pcpn_l_r, vmin = 0.5, vmax = 0.95) |
---|
105 | plt.axis([-180., 180., 0, Njr]) |
---|
106 | cb = plt.colorbar() |
---|
107 | cb.set_label('CH ' + str(chan[ich]) + 'GHz') |
---|
108 | plt.yticks(np.array([0, 31, 59, 90, 120, 151, 181]), np.array(['JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', 'JUL'])) |
---|
109 | plt.xticks(np.arange(-180, 200, 20), rotation = 30) |
---|
110 | plt.grid(True) |
---|
111 | plt.savefig('/usr/home/lahlod/twice_d/figures_output_ANTARC/AMSUA/hov_emis_lat' + str(lat[ilat]) + '_ch1-ch15_AMSUA.png') |
---|
112 | |
---|
113 | |
---|
114 | |
---|
115 | ######################### |
---|
116 | ## anomalie emissivite ## |
---|
117 | ######################### |
---|
118 | emis_anom_glob = np.zeros([len(chan), len(lat), len(lon), Njr], float) |
---|
119 | for ich in range (0, len(chan)): |
---|
120 | for ilon in range (0, len(lon)): |
---|
121 | for ilat in range (0, len(lat)): |
---|
122 | for indjr in range (15, Njr - 15): # moyenne glissante centree |
---|
123 | emis_anom_glob[ich, ilat, ilon, indjr] = emis_glob[ich, ilat, ilon, indjr][nonzero(emis_glob[ich, ilat, ilon, indjr] != 0.)] - mean(emis_glob[ich, ilat, ilon, indjr - 15 : indjr + 15][nonzero(emis_glob[ich, ilat, ilon, indjr - 15 : indjr + 15] != 0.)]) |
---|
124 | |
---|
125 | |
---|
126 | |
---|
127 | ################################################################### |
---|
128 | ## variation saisonniÚre emissivite en fonction de type de glace ## |
---|
129 | ################################################################### |
---|
130 | |
---|
131 | #### definition des zones d etaude #### |
---|
132 | |
---|
133 | # zone1 et 2 = lon ([30E - 70E] et [90E - 160E]) tranche de lat [70S] |
---|
134 | lat1 = np.where(lat == -70.)[0][0] |
---|
135 | lon1 = np.where(lon == 30.)[0][0] |
---|
136 | lon2 = np.where(lon == 75.)[0][0] |
---|
137 | lon3 = np.where(lon == 90.)[0][0] |
---|
138 | lon4 = np.where(lon == 165.)[0][0] |
---|
139 | # zone3 = lon ([100W - 80W] tranche de lat [80S] |
---|
140 | lat2 = np.where(lat == -80.)[0][0] |
---|
141 | lon5 = np.where(lon == -100.)[0][0] |
---|
142 | lon6 = np.where(lon == -75.)[0][0] |
---|
143 | # zone4 = lon ([20W - 20E] tranche de lat [75S] |
---|
144 | lat3 = np.where(lat == -75.)[0][0] |
---|
145 | lon7 = np.where(lon == -20.)[0][0] |
---|
146 | lon8 = np.where(lon == 25.)[0][0] |
---|
147 | # zone5 = lon ([90E - 160E] tranche de lat [80S] |
---|
148 | lon3 = np.where(lon == 90.)[0][0] |
---|
149 | lon4 = np.where(lon == 165.)[0][0] |
---|
150 | # zone6 = lon ([180W - 160W] tranche de lat [75S] |
---|
151 | lon9 = np.where(lon == -180.)[0][0] |
---|
152 | lon10 = np.where(lon == -155.)[0][0] |
---|
153 | # zone7 = lon ([180E - 160E] tranche de lat [80S] |
---|
154 | lon11 = np.where(lon == 160.)[0][0] |
---|
155 | lon12 = np.where(lon == 180.)[0][0] |
---|
156 | # zone8 = lon ([180E - 160E] tranche de lat [75S] |
---|
157 | lon11 = np.where(lon == 160.)[0][0] |
---|
158 | lon12 = np.where(lon == 180.)[0][0] |
---|
159 | # zone9 = lon ([40E - 90E] tranche de lat [80S] |
---|
160 | lon13 = np.where(lon == 40.)[0][0] |
---|
161 | lon14 = np.where(lon == 95.)[0][0] |
---|
162 | # zone10 = lon ([120W - 140W] tranche de lat [80S] |
---|
163 | lon15 = np.where(lon == -140.)[0][0] |
---|
164 | lon16 = np.where(lon == -115.)[0][0] |
---|
165 | |
---|
166 | #### emissivite #### |
---|
167 | |
---|
168 | emis_glob_jr = np.zeros([len(chan), Njr], float) |
---|
169 | for ich in range (0, len(chan)): |
---|
170 | for indjr in range (0, Njr): |
---|
171 | emis_glob_jr[ich, indjr] = mean(emis_glob[ich, lat2, lon3 : lon4, indjr]) |
---|
172 | |
---|
173 | |
---|
174 | plt.ion() |
---|
175 | plt.figure() |
---|
176 | color_chan = np.array(['b', 'r']) |
---|
177 | for ich in range (0, len(chan)): |
---|
178 | plt.plot(emis_glob_jr[ich, :][nonzero(emis_glob_jr[ich, :] != 0.)], color_chan[ich], label = str(chan[ich]) + 'GHz') |
---|
179 | |
---|
180 | legend(loc = 4) |
---|
181 | plt.xlim(0, 212) |
---|
182 | plt.xticks(np.array([0, 31, 59, 90, 120, 151, 181]), np.array(['JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', 'JUL']), rotation = 30) |
---|
183 | plt.ylabel('emissivitty') |
---|
184 | grid(True) |
---|
185 | plt.title('Seasonal variation of emissivity - zone[120W-140W;80S] - AMSUA (nadir)') |
---|
186 | |
---|
187 | #### anomalie d emissivite (sur une moyenne glissante) #### |
---|
188 | |
---|
189 | emis_anom_glob_jr = np.zeros([len(chan), Njr], float) |
---|
190 | for ich in range (0, len(chan)): |
---|
191 | for indjr in range (0, Njr): |
---|
192 | emis_anom_glob_jr[ich, indjr] = mean(emis_anom_glob[ich, lat3, lon11 : lon12, indjr][nonzero(isnan(emis_anom_glob[ich, lat3, lon11 : lon12, indjr]) == False)]) |
---|
193 | |
---|
194 | |
---|
195 | plt.ion() |
---|
196 | plt.figure() |
---|
197 | color_chan = np.array(['b', 'r']) |
---|
198 | for ich in range (0, len(chan)): |
---|
199 | plt.plot(emis_anom_glob_jr[ich, :], color_chan[ich], label = str(chan[ich]) + 'GHz') |
---|
200 | |
---|
201 | legend(loc = 4) |
---|
202 | plt.xlim(0, 212) |
---|
203 | plt.xticks(np.array([0, 31, 59, 90, 120, 151, 181]), np.array(['JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', 'JUL']), rotation = 30) |
---|
204 | plt.ylabel('emissivity anomaly') |
---|
205 | grid(True) |
---|
206 | plt.title('Seasonal variation of emissivity anomaly - zone[30E-70E;70S] - AMSUA (nadir)') |
---|
207 | |
---|
208 | #### gradient de frequence #### |
---|
209 | |
---|
210 | monthly_emis = np.array([monthly_emis_JAN, monthly_emis_FEB, monthly_emis_MAR, monthly_emis_APR, monthly_emis_MAY, monthly_emis_JUN, monthly_emis_JUL]) |
---|
211 | monthly_emis_mo = np.zeros([len(month), len(chan)], float) |
---|
212 | for imo in range (0, len(month)): |
---|
213 | for ich in range(0, len(chan)): |
---|
214 | monthly_emis_mo[imo, ich] = mean(monthly_emis[imo, ich, lat2, lon15 : lon16]) |
---|
215 | |
---|
216 | |
---|
217 | monthly_emis_diff = monthly_emis_mo[:, 0] / monthly_emis_mo[:, 1] |
---|
218 | |
---|
219 | |
---|
220 | plt.ion() |
---|
221 | plt.figure() |
---|
222 | color_chan = np.array(['b', 'r']) |
---|
223 | subplot(2,1,1) |
---|
224 | for ich in range (0, len(chan)): |
---|
225 | plt.plot(monthly_emis_mo[:, ich], color_chan[ich], label = str(chan[ich]) + 'GHz') |
---|
226 | |
---|
227 | plt.ylabel('emissivity') |
---|
228 | legend(loc = 4) |
---|
229 | grid(True) |
---|
230 | plt.title('Seasonal emissivity bias - zone[140W-120W;80S] - AMSUA (nadir)') |
---|
231 | plt.xticks(np.arange(0, 7, 1), np.array(['JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', 'JUL']), rotation = 30) |
---|
232 | subplot(2,1,2) |
---|
233 | plot(monthly_emis_diff, 'k', label = 'bias') |
---|
234 | plot(np.ones([len(month)]), 'g--') |
---|
235 | plt.ylabel('bias (CH1/CH15)') |
---|
236 | plt.xticks(np.arange(0, 7, 1), np.array(['JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN', 'JUL']), rotation = 30) |
---|
237 | grid(True) |
---|
238 | |
---|
239 | |
---|
240 | |
---|
241 | |
---|
242 | ############################### |
---|
243 | ## biais emissivite par mois ## |
---|
244 | ############################### |
---|
245 | biais_month = np.zeros([len(month), len(lat), len(lon)], float) |
---|
246 | for imo in range (0, len(month)): |
---|
247 | biais_month[imo, :, :] = monthly_emis[imo][0, :, :] - monthly_emis[imo][1, :, :] |
---|
248 | |
---|
249 | for imo in range (0, len(month)): |
---|
250 | plt.figure() |
---|
251 | plt.title('AMSUA [CH1-CH2] - '+str(month[imo])+' - (lat -80)') |
---|
252 | vlines(lon, [0], biais_month[imo, 2, :]) |
---|
253 | xlabel('longitude') |
---|
254 | ylabel('bias(emissivity)') |
---|
255 | xlim(-180, 180) |
---|
256 | xticks(np.arange(-180., 200., 20), rotation = 30) |
---|
257 | grid(True) |
---|
258 | plt.savefig('/usr/home/lahlod/twice_d/figures_output_ANTARC/AMSUA/vlines_bias_emis_ch1-ch15_lat-80_'+str(month[imo])+'_AMSUA.png') |
---|
259 | |
---|
260 | |
---|
261 | ############################# |
---|
262 | ## var emissivite par mois ## |
---|
263 | ############################# |
---|
264 | var = np.zeros([len(lat), len(lon)], float) |
---|
265 | for ilon in range (0, len(lon)): |
---|
266 | for ilat in range (0, len(lat)): |
---|
267 | var[ilat, ilon] = np.sum((biais_month[:, ilat, ilon] - mean(biais_month[:, ilat, ilon]))**2) |
---|
268 | |
---|
269 | var = var/7 |
---|
270 | for ilat in range (2, 7): |
---|
271 | plt.figure() |
---|
272 | vlines(lon, [0], var[ilat, :]) |
---|
273 | xlabel('longitude') |
---|
274 | ylabel('var(emiss_bias)') |
---|
275 | plt.title('AMSUA [CH1-CH2] - lat '+str(lat[ilat]) ) |
---|
276 | xlim(-180, 180) |
---|
277 | xticks(np.arange(-180., 200., 20), rotation = 30) |
---|
278 | grid(True) |
---|
279 | plt.savefig('/usr/home/lahlod/twice_d/figures_output_ANTARC/AMSUA/vlines_var_bias-emis_ch1-ch15_lat-'+str(lat[ilat])+'_AMSUA.png') |
---|
280 | |
---|
281 | |
---|
282 | |
---|
283 | |
---|
284 | |
---|
285 | |
---|
286 | |
---|
287 | |
---|
288 | |
---|
289 | |
---|
290 | |
---|
291 | |
---|
292 | |
---|
293 | |
---|
294 | |
---|
295 | |
---|
296 | |
---|
297 | |
---|
298 | plt.figure() |
---|
299 | plt.ion() |
---|
300 | #vlines(biais_month[0, 2]) |
---|
301 | #plt.subplot(2, 1, 1) |
---|
302 | m = Basemap(llcrnrlon=-180, urcrnrlon=180, llcrnrlat=-90, urcrnrlat=-30, projection='cyl', resolution='c', fix_aspect=True) |
---|
303 | m.drawcoastlines(linewidth = 1) |
---|
304 | m.drawparallels(np.arange(-90., -30., 20)) |
---|
305 | m.drawmeridians(np.arange(-180., 180., 20)) |
---|
306 | xii,yii = m(*np.meshgrid(lon, lat)) |
---|
307 | clevs = (arange(0., 0.05, 0.0001)) |
---|
308 | cs = m.contourf(xii, yii, var, clevs, cmap=cm.s3pcpn_l_r) |
---|
309 | cbar = colorbar(cs) |
---|
310 | cbar.set_label('CH1 - 23.8GHz') |
---|
311 | plt.title('var(emis) JUL - AMSUA') |
---|
312 | xticks(arange(-180, 200, 20), rotation = 45) |
---|
313 | yticks(arange(-90, -10, 20)) |
---|
314 | plt.subplot(2, 1, 2) |
---|
315 | m = Basemap(llcrnrlon=-180, urcrnrlon=180, llcrnrlat=-90, urcrnrlat=-30, projection='cyl', resolution='c', fix_aspect=True) |
---|
316 | m.drawcoastlines(linewidth = 1) |
---|
317 | m.drawparallels(np.arange(-90., -30., 20)) |
---|
318 | m.drawmeridians(np.arange(-180., 180., 20)) |
---|
319 | xii,yii = m(*np.meshgrid(lon, lat)) |
---|
320 | clevs = (arange(0., 0.0035, 0.00001)) |
---|
321 | cs = m.contourf(xii, yii, var_JUL[1, :, :], clevs, cmap=cm.s3pcpn_l_r) |
---|
322 | cbar = colorbar(cs) |
---|
323 | cbar.set_label('CH15 - 89GHz') |
---|
324 | xticks(arange(-180, 200, 20), rotation = 45) |
---|
325 | yticks(arange(-90, -10, 20)) |
---|
326 | |
---|
327 | |
---|
328 | |
---|
329 | |
---|
330 | |
---|