1 | #!/usr/bin/env python |
---|
2 | # -*- coding: utf-8 -*- |
---|
3 | import string |
---|
4 | import numpy as np |
---|
5 | import matplotlib.pyplot as plt |
---|
6 | import ffgrid2 |
---|
7 | from pylab import * |
---|
8 | from mpl_toolkits.basemap import Basemap |
---|
9 | from mpl_toolkits.basemap import shiftgrid, cm |
---|
10 | import netCDF4 |
---|
11 | import draw_map |
---|
12 | |
---|
13 | |
---|
14 | |
---|
15 | fichier=open('SSMIS_CH2_ANTARC_JANUARY2010.DAT','r') |
---|
16 | numlines = 0 |
---|
17 | for line in fichier: numlines += 1 |
---|
18 | |
---|
19 | fichier.close |
---|
20 | |
---|
21 | |
---|
22 | fichier=open('SSMIS_CH2_ANTARC_JANUARY2010.DAT','r') |
---|
23 | nbtotal=numlines |
---|
24 | |
---|
25 | iligne=0 |
---|
26 | lat=np.zeros([nbtotal],float) |
---|
27 | lon=np.zeros([nbtotal],float) |
---|
28 | jjr=np.zeros([nbtotal],float) |
---|
29 | zen=np.zeros([nbtotal],float) |
---|
30 | fov=np.zeros([nbtotal],float) |
---|
31 | ts=np.zeros([nbtotal],float) |
---|
32 | emis=np.zeros([nbtotal],float) |
---|
33 | tb=np.zeros([nbtotal],float) |
---|
34 | tup=np.zeros([nbtotal],float) |
---|
35 | tdn=np.zeros([nbtotal],float) |
---|
36 | trans=np.zeros([nbtotal],float) |
---|
37 | orog=np.zeros([nbtotal],float) |
---|
38 | |
---|
39 | while (iligne < nbtotal-1) : |
---|
40 | line=fichier.readline() |
---|
41 | # exemple : line = "0.22 2.3 5.0 6" |
---|
42 | liste = line.split() |
---|
43 | # exemple : listeCoord ['0.22', '2.3', '5.0', '6'] (liste de chaine de |
---|
44 | # caractÚres) |
---|
45 | lat[iligne] = float(liste[1]) |
---|
46 | lon[iligne] = float(liste[0]) |
---|
47 | jjr[iligne] = float(liste[4]) |
---|
48 | ts[iligne] = float(liste[10]) |
---|
49 | tb[iligne] = float(liste[13]) |
---|
50 | emis[iligne] = float(liste[16]) |
---|
51 | orog[iligne] = float(liste[13]) |
---|
52 | iligne=iligne+1 |
---|
53 | |
---|
54 | fichier.close |
---|
55 | |
---|
56 | |
---|
57 | |
---|
58 | dx=1.0 |
---|
59 | dy=1.0 |
---|
60 | x0, x1 = -180, 180 |
---|
61 | y0, y1 = -90, 90 |
---|
62 | |
---|
63 | monthly_outz=np.zeros([181,361],float) |
---|
64 | monthly_lon=np.zeros([361]) |
---|
65 | monthly_lat=np.zeros([181]) |
---|
66 | xx = lon |
---|
67 | yy = lat |
---|
68 | zz = tb |
---|
69 | zz0 = 100 |
---|
70 | zz1= 300 |
---|
71 | outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,zz0, zz1) |
---|
72 | monthly_outz=outz |
---|
73 | monthly_lon=outx |
---|
74 | monthly_lat=outy |
---|
75 | del outz, outx, outy, zz, xx, yy |
---|
76 | |
---|
77 | |
---|
78 | # ici je fais des cartes moyennes en melangeant les polars |
---|
79 | |
---|
80 | xx = lon_ssmis |
---|
81 | yy = lat_ssmis |
---|
82 | zz = 0.5*(emis_ssmis[1,:]+emis_ssmis[2,:]) |
---|
83 | outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,0.1, 1) |
---|
84 | monthly_outz_ssmis_polar[0,:,:]=outz |
---|
85 | del outz, outx, outy, zz |
---|
86 | |
---|
87 | zz = 0.5*(emis_ssmis[4,:]+emis_ssmis[5,:]) |
---|
88 | outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,0.1, 1) |
---|
89 | monthly_outz_ssmis_polar[1,:,:]=outz |
---|
90 | del outz, outx, outy, zz |
---|
91 | |
---|
92 | zz = 0.5*(emis_ssmis[6,:]+emis_ssmis[7,:]) |
---|
93 | outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,0.1, 1) |
---|
94 | monthly_outz_ssmis_polar[2,:,:]=outz |
---|
95 | del outz, outx, outy, zz, xx, yy |
---|
96 | |
---|
97 | # ici je fais des cartes moyennes des differences des polars |
---|
98 | xx = lon_ssmis |
---|
99 | yy = lat_ssmis |
---|
100 | zz = emis_ssmis[1,:]-emis_ssmis[2,:] |
---|
101 | outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,-0.05, 0.2) |
---|
102 | monthly_outz_ssmis_diff[0,:,:]=outz |
---|
103 | del outz, outx, outy, zz |
---|
104 | |
---|
105 | zz = emis_ssmis[4,:]-emis_ssmis[5,:] |
---|
106 | outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,-0.05, 0.2) |
---|
107 | monthly_outz_ssmis_diff[1,:,:]=outz |
---|
108 | del outz, outx, outy, zz |
---|
109 | |
---|
110 | zz = emis_ssmis[6,:]-emis_ssmis[7,:] |
---|
111 | outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,-0.05, 0.2) |
---|
112 | monthly_outz_ssmis_diff[2,:,:]=outz |
---|
113 | del outz, outx, outy, zz, xx, yy |
---|
114 | |
---|
115 | 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) |
---|
116 | 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) |
---|
117 | 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) |
---|
118 | |
---|
119 | 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) |
---|
120 | 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) |
---|
121 | 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) |
---|
122 | |
---|
123 | 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) |
---|
124 | 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) |
---|
125 | 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) |
---|
126 | 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) |
---|
127 | 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) |
---|
128 | 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) |
---|
129 | 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) |
---|
130 | 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) |
---|
131 | |
---|
132 | |
---|
133 | bins=arange(0.3,1,0.001) |
---|
134 | bb=(lat_ssmis >= 0) |
---|
135 | |
---|
136 | plt.hist(emis_ssmis[0,nonzero(bb)[0]], bins=bins,histtype='step', label='e50V',normed='True',color='#4BB5C1') |
---|
137 | plt.hist(emis_ssmis[1,nonzero(bb)[0]], bins=bins,histtype='step', label='e19V',normed='True',color='black') |
---|
138 | plt.hist(emis_ssmis[3,nonzero(bb)[0]], bins=bins,histtype='step', label='e22V',normed='True',color='#B9121B') |
---|
139 | plt.hist(emis_ssmis[4,nonzero(bb)[0]], bins=bins,histtype='step', label='e37V',normed='True',color='#9748D4') |
---|
140 | plt.hist(emis_ssmis[6,nonzero(bb)[0]], bins=bins,histtype='step', label='e91V',normed='True',color='#060DE5') |
---|
141 | plt.legend(loc='upper left') |
---|
142 | plt.show() |
---|
143 | plt.savefig('..\FIG\hist_ssmis_V_NH_'+le_mois+'.png') |
---|
144 | close() |
---|
145 | |
---|
146 | bins=arange(0.3,1,0.001) |
---|
147 | plt.hist(emis_ssmis[2,nonzero(bb)[0]], bins=bins,histtype='step', label='e19H',normed='True',color='black') |
---|
148 | plt.hist(emis_ssmis[5,nonzero(bb)[0]], bins=bins,histtype='step', label='e37H',normed='True',color='#9748D4') |
---|
149 | plt.hist(emis_ssmis[7,nonzero(bb)[0]], bins=bins,histtype='step', label='e91H',normed='True',color='#060DE5') |
---|
150 | plt.legend(loc='upper left') |
---|
151 | plt.show() |
---|
152 | plt.savefig('..\FIG\hist_ssmis_H_NH_'+le_mois+'.png') |
---|
153 | close() |
---|
154 | |
---|
155 | bins=arange(0.3,1,0.001) |
---|
156 | 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') |
---|
157 | 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') |
---|
158 | 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') |
---|
159 | plt.legend(loc='upper left') |
---|
160 | plt.show() |
---|
161 | plt.savefig('..\FIG\hist_ssmis_mpolar_NH_'+le_mois+'.png') |
---|
162 | close() |
---|
163 | |
---|
164 | |
---|
165 | # stats quotidienne autour de la station Thulé |
---|
166 | lat_stations=[76.32, 74.43, 78.13, 58.45, 68.6, 64.58] |
---|
167 | lon_stations=[-68.3, -94.59, 15.35, -78.08, 33.1, 40.5] |
---|
168 | nom_stations=['Thule', 'Resolute', 'Longyearbyen', 'Iqaluit', 'Murmansk', 'Arkhangelsk'] |
---|
169 | |
---|
170 | |
---|
171 | for sta in range(0,6): |
---|
172 | lat0=lat_stations[sta] |
---|
173 | lon0=lon_stations[sta] |
---|
174 | stat_jour=np.zeros([8,7,31],float) |
---|
175 | clear bb |
---|
176 | for canal in range(0,8): |
---|
177 | for jjr in range(0,31): |
---|
178 | jour_obs=jjr+1 |
---|
179 | bb=(jjr_ssmis == jour_obs) & (abs(lat_ssmis-lat0) < 2.) & (abs(lon_ssmis-lon0) < 2.) |
---|
180 | stat_jour[canal,0,jjr]=mean(emis_ssmis[canal,nonzero(bb)[0]]) |
---|
181 | stat_jour[canal,1,jjr]=std(emis_ssmis[canal,nonzero(bb)[0]]) |
---|
182 | stat_jour[canal,2,jjr]=size(nonzero(bb)) |
---|
183 | stat_jour[canal,3,jjr]=mean(ts_ssmis[nonzero(bb)[0]]) |
---|
184 | stat_jour[canal,4,jjr]=std(ts_ssmis[nonzero(bb)[0]]) |
---|
185 | stat_jour[canal,5,jjr]=mean(tb_ssmis[canal,nonzero(bb)[0]]) |
---|
186 | stat_jour[canal,6,jjr]=std(tb_ssmis[canal,nonzero(bb)[0]]) |
---|
187 | del bb |
---|
188 | np.save('STAT_SSMIS_'+nom_stations[sta]+'_'+le_mois+'.dat', stat_jour) |
---|
189 | del stat_jour |
---|
190 | |
---|
191 | mpolar_ssmis=np.zeros([3,nbtotal],float) |
---|
192 | mpolar_ssmis[0,:]=0.5*(emis_ssmis[1,:]+emis_ssmis[2,:]) |
---|
193 | mpolar_ssmis[1,:]=0.5*(emis_ssmis[4,:]+emis_ssmis[5,:]) |
---|
194 | mpolar_ssmis[2,:]=0.5*(emis_ssmis[6,:]+emis_ssmis[7,:]) |
---|
195 | |
---|
196 | mpolarTB_ssmis=np.zeros([3,nbtotal],float) |
---|
197 | mpolarTB_ssmis[0,:]=0.5*(tb_ssmis[1,:]+tb_ssmis[2,:]) |
---|
198 | mpolarTB_ssmis[1,:]=0.5*(tb_ssmis[4,:]+tb_ssmis[5,:]) |
---|
199 | mpolarTN_ssmis[2,:]=0.5*(tb_ssmis[6,:]+tb_ssmis[7,:]) |
---|
200 | |
---|
201 | for sta in range(0,6): |
---|
202 | lat0=lat_stations[sta] |
---|
203 | lon0=lon_stations[sta] |
---|
204 | stat2_jour=np.zeros([3,7,31],float) |
---|
205 | clear bb |
---|
206 | for canal in range(0,3): |
---|
207 | for jjr in range(0,31): |
---|
208 | jour_obs=jjr+1 |
---|
209 | bb=(jjr_ssmis == jour_obs) & (abs(lat_ssmis-lat0) < 2.) & (abs(lon_ssmis-lon0) < 2.) |
---|
210 | stat2_jour[canal,0,jjr]=mean(mpolar_ssmis[canal,nonzero(bb)[0]]) |
---|
211 | stat2_jour[canal,1,jjr]=std(mpolar_ssmis[canal,nonzero(bb)[0]]) |
---|
212 | stat2_jour[canal,2,jjr]=size(nonzero(bb)) |
---|
213 | stat2_jour[canal,3,jjr]=mean(ts_ssmis[nonzero(bb)[0]]) |
---|
214 | stat2_jour[canal,4,jjr]=std(ts_ssmis[nonzero(bb)[0]]) |
---|
215 | stat2_jour[canal,5,jjr]=mean(tb_ssmis[canal,nonzero(bb)[0]]) |
---|
216 | stat2_jour[canal,6,jjr]=std(tb_ssmis[canal,nonzero(bb)[0]]) |
---|
217 | del bb |
---|
218 | np.save('STAT_SSMIS-MPOLAR_'+nom_stations[sta]+'_'+le_mois+'.dat', stat_jour) |
---|
219 | del stat_jour |
---|
220 | |
---|
221 | # ecriture sous format nc |
---|
222 | from netCDF4 import Dataset |
---|
223 | rootgrp = Dataset('..\EMIS\EMIS_SSMIS_'+le_mois+'.nc', 'w', format='NETCDF4') |
---|
224 | |
---|
225 | rootgrp.createDimension('longitude', len(monthly_lon_ssmis)) |
---|
226 | rootgrp.createDimension('latitude', len(monthly_lat_ssmis)) |
---|
227 | rootgrp.createDimension('channels', 8) |
---|
228 | rootgrp.createDimension('bchannels', 3) |
---|
229 | |
---|
230 | # createVariable (nom de la variable, type, dimensions) |
---|
231 | # Si 1 dimension, ne pas oublier la virgule |
---|
232 | nclon = rootgrp.createVariable('longitude', 'f8', ('longitude',)) |
---|
233 | nclat = rootgrp.createVariable('latitude', 'f8', ('latitude',)) |
---|
234 | ncchan=rootgrp.createVariable('channels', 'f', ('channels',)) |
---|
235 | ncchan2=rootgrp.createVariable('bchannels', 'f', ('bchannels',)) |
---|
236 | nctemp = rootgrp.createVariable('emissivity', 'f8', ('channels','latitude', 'longitude')) |
---|
237 | nctemp2 = rootgrp.createVariable('emissivity melange polar', 'f8', ('bchannels','latitude', 'longitude')) |
---|
238 | |
---|
239 | nclon[:] = monthly_lon_ssmis |
---|
240 | nclat[:] = monthly_lat_ssmis |
---|
241 | ncchan[:]=[50,19.1,19.2,22,37.1,37.2,91.1,91.2] |
---|
242 | ncchan2[:]=[19,37,91] |
---|
243 | nctemp[:] = monthly_outz_ssmis |
---|
244 | nctemp2[:] = monthly_outz_ssmis_polar |
---|
245 | |
---|
246 | rootgrp.close() |
---|