1 | PRO plt_map, cmd, iplot, win, iover, landscape |
---|
2 | |
---|
3 | @common |
---|
4 | @com_eg |
---|
5 | ;; |
---|
6 | ;;------------------------------------- |
---|
7 | ;; |
---|
8 | ;; make window plot |
---|
9 | ;; |
---|
10 | ;; EG 19-2-99 |
---|
11 | ;;------------------------------------- |
---|
12 | ; |
---|
13 | ; |
---|
14 | ; cmd = command line of window |
---|
15 | ; win = position of window (petit[,,,]) |
---|
16 | ; iover = overlay index |
---|
17 | IF debug_w THEN print, ' ' |
---|
18 | IF debug_w THEN print, ' ENTER plt_map...' |
---|
19 | ; |
---|
20 | ; empty window case for on=2 |
---|
21 | ; |
---|
22 | ioverchk = iover |
---|
23 | |
---|
24 | cmd_wrk = cmd |
---|
25 | |
---|
26 | IF cmd.on EQ 2 THEN return |
---|
27 | ; |
---|
28 | ; incompatible options of plt_def and post_it |
---|
29 | ; |
---|
30 | trend_typp = cmd.trend |
---|
31 | ; |
---|
32 | ; ============= |
---|
33 | ; 1.read data |
---|
34 | ; ============= |
---|
35 | ; |
---|
36 | ; type of plot |
---|
37 | |
---|
38 | ; time series of rms deviation on a sigma surface? |
---|
39 | st_rms = 0 |
---|
40 | IF strpos(cmd.plt, 'st') NE -1 AND strpos(cmd.var, '@s') GT 4 THEN BEGIN |
---|
41 | st_rms = 1 |
---|
42 | sig_surface = strmid(cmd.var, strpos(cmd.var, '@s')+2, strlen(cmd.var)) |
---|
43 | cmd.var = strmid(cmd.var, 0, strpos(cmd.var, '@s')) |
---|
44 | legendsurf = ' on sigma='+sig_surface |
---|
45 | ENDIF |
---|
46 | |
---|
47 | CASE cmd.plt OF |
---|
48 | 'x': cmd.plt = cmd.plt+'_global' |
---|
49 | 'y': cmd.plt = cmd.plt+'_global' |
---|
50 | 'z': cmd.plt = cmd.plt+'_global' |
---|
51 | 's': cmd.plt = cmd.plt+'_global' |
---|
52 | 'yz': cmd.plt = cmd.plt+'_global' |
---|
53 | ELSE: |
---|
54 | ENDCASE |
---|
55 | |
---|
56 | ; Definition of plttyp, hotyp, dimplot |
---|
57 | datyp = def_dptyp(cmd) |
---|
58 | |
---|
59 | plttyp = datyp.plttyp |
---|
60 | hotyp = datyp.hotyp |
---|
61 | dimplot = datyp.dimplot |
---|
62 | pltztyp = datyp.pltztyp |
---|
63 | plt1dtyp = datyp.plt1dtyp |
---|
64 | splot = datyp.splot |
---|
65 | |
---|
66 | hotypchk = hotyp |
---|
67 | |
---|
68 | ; stat computations: |
---|
69 | |
---|
70 | CASE cmd.plt OF |
---|
71 | 'xyt': BEGIN & plttyp = 'plt' & hotyp = 'xyt' & dimplot = 2 & END |
---|
72 | ELSE: |
---|
73 | ENDCASE |
---|
74 | |
---|
75 | ; xy plot on density surface (cmd.var = <var>@s<sig_surf> |
---|
76 | legsurf = '' |
---|
77 | sig_surf = 0. |
---|
78 | IF strpos(cmd.var, '@s') GT 4 THEN BEGIN |
---|
79 | splot = 1 |
---|
80 | sig_surf = strmid(cmd.var, strpos(cmd.var, '@s')+2, strlen(cmd.var)) |
---|
81 | cmd.var = strmid(cmd.var, 0, strpos(cmd.var, '@s')) |
---|
82 | legsurf = ' on sigma='+sig_surf |
---|
83 | ENDIF |
---|
84 | |
---|
85 | ; contour of coasts |
---|
86 | |
---|
87 | CASE splot OF |
---|
88 | 1: c_cont = 0 |
---|
89 | 0: c_cont = (!d.n_colors-1) < 255 |
---|
90 | ENDCASE |
---|
91 | |
---|
92 | ; scatter plot y=f(x) |
---|
93 | |
---|
94 | IF strpos(cmd.var, '=f(') GE 1 THEN BEGIN |
---|
95 | CASE cmd.plt OF |
---|
96 | 'xy': cmd.plt = cmd.plt+'_global' |
---|
97 | ELSE: |
---|
98 | ENDCASE |
---|
99 | plttyp = 'yfx' |
---|
100 | dimplot = 1 |
---|
101 | ENDIF |
---|
102 | |
---|
103 | ; binning plot y=bin(x) |
---|
104 | |
---|
105 | IF strpos(cmd.var, '=bin(') GE 1 THEN BEGIN |
---|
106 | CASE cmd.plt OF |
---|
107 | 'xy': cmd.plt = cmd.plt+'_global' |
---|
108 | ELSE: |
---|
109 | ENDCASE |
---|
110 | plttyp = 'ybinx' |
---|
111 | dimplot = 1 |
---|
112 | ENDIF |
---|
113 | |
---|
114 | ; vector plot |
---|
115 | |
---|
116 | IF strpos(cmd.var, '[') EQ 0 THEN BEGIN |
---|
117 | vecplot = 1 |
---|
118 | ; sampling of vector |
---|
119 | IF vector_sample GE 2 THEN BEGIN |
---|
120 | vect_sample = ',UNVECTSUR=strtrim(string(vector_sample),2)' |
---|
121 | ENDIF ELSE vect_sample = '' |
---|
122 | ENDIF ELSE BEGIN |
---|
123 | vecplot = 0 |
---|
124 | ENDELSE |
---|
125 | |
---|
126 | ; spectrum |
---|
127 | |
---|
128 | spectplt = 0 |
---|
129 | IF strpos(cmd.plt, '@s') GE 1 THEN BEGIN |
---|
130 | IF hotyp EQ 't' THEN spectplt = 1 ELSE BEGIN |
---|
131 | print, ' *** spectrum [@s] applies to 1D time serie only (t_<box>)' |
---|
132 | stop |
---|
133 | ENDELSE |
---|
134 | cmd.plt = strmid(cmd.plt, 0, strpos(cmd.plt, '@s')) |
---|
135 | ENDIF |
---|
136 | |
---|
137 | ; wavelet |
---|
138 | |
---|
139 | waveplt = 0 |
---|
140 | IF strpos(cmd.plt, '@w') GE 1 THEN BEGIN |
---|
141 | IF hotyp EQ 't' THEN waveplt = 1 ELSE BEGIN |
---|
142 | print, ' *** wavelet [@w] applies to 1D time serie only (t_<box>)' |
---|
143 | stop |
---|
144 | ENDELSE |
---|
145 | cmd.plt = strmid(cmd.plt, 0, strpos(cmd.plt, '@w')) |
---|
146 | ENDIF |
---|
147 | |
---|
148 | ; running std dev |
---|
149 | |
---|
150 | run_stddev = 0 |
---|
151 | IF strpos(cmd.plt, '@r') GE 1 THEN BEGIN |
---|
152 | IF hotyp NE 't' THEN BEGIN |
---|
153 | print, ' *** running [@r<win>] std dev applies to 1D time serie only (t_<box> & 1m@t412)' |
---|
154 | stop |
---|
155 | ENDIF |
---|
156 | run_stddev = strmid(cmd.plt, strpos(cmd.plt, '@r')+2, strlen(cmd.plt)) |
---|
157 | cmd.plt = strmid(cmd.plt, 0, strpos(cmd.plt, '@r')) |
---|
158 | ENDIF |
---|
159 | |
---|
160 | |
---|
161 | ; define grids |
---|
162 | ; read if new type + triangule |
---|
163 | ; |
---|
164 | cmd1_back = cmd |
---|
165 | sw_diffg = 0 |
---|
166 | def_grid, cmd |
---|
167 | IF debug_w THEN print, ' after def_grid in plt_map: cmd.grid/varexp/vargrid ',cmd.grid,' ', varexp, ' ', vargrid |
---|
168 | |
---|
169 | ; vertical average if 3D field |
---|
170 | vert_switch = 0 |
---|
171 | IF plttyp EQ 'plt' AND vert_type ne '0' THEN vert_switch = 1 |
---|
172 | IF plttyp EQ 'pltt' AND vert_type ne '0' THEN vert_switch = 1 |
---|
173 | IF plttyp EQ 'plt1d' AND vert_type ne '0' THEN vert_switch = 1 |
---|
174 | IF plt1dtyp EQ 'z' THEN vert_switch = 0 |
---|
175 | |
---|
176 | ; read all data if netcdf output |
---|
177 | |
---|
178 | IF debug_w THEN print, ' force_all_data_read in pltmap: ', force_all_data_read |
---|
179 | |
---|
180 | IF cmd.out NE 'cdf' THEN BEGIN |
---|
181 | all_data_str = '' |
---|
182 | ENDIF ELSE BEGIN |
---|
183 | all_data_str = ',/ALL_DATA' |
---|
184 | END |
---|
185 | |
---|
186 | ; =========== |
---|
187 | ; Read data |
---|
188 | ; =========== |
---|
189 | |
---|
190 | really_1m_st = 1 |
---|
191 | IF cmd.timave EQ '1y' AND strmid(cmd.plt, 0, 2) EQ 'st' THEN BEGIN |
---|
192 | @densit_pltmap_read |
---|
193 | ENDIF ELSE BEGIN |
---|
194 | pltcmd = 'field = data_read(cmd,'''+hotyp+''','''+plttyp+''','+string(dimplot)+','+string(iover)+all_data_str+', ZMTYP = '''+cmd.plt+''')' |
---|
195 | printf, nulhis, strcompress(pltcmd) |
---|
196 | printf, nulhis, ' ' |
---|
197 | IF execute(pltcmd) EQ 0 THEN stop |
---|
198 | ENDELSE |
---|
199 | |
---|
200 | ; data specific treatments |
---|
201 | IF cmd.var EQ 'sla' THEN field.origin = 'diff' |
---|
202 | IF cmd.var EQ 'sosstsst' AND cmd.exp EQ 'HadISST1' THEN BEGIN |
---|
203 | ; limit data to -1.8 C |
---|
204 | idx = where(field.data EQ valmask) |
---|
205 | field.data = field.data > (-1.8) |
---|
206 | IF idx[0] NE -1 THEN field.data(idx) = valmask |
---|
207 | ENDIF |
---|
208 | |
---|
209 | IF n_elements(field.data) EQ 1 THEN return |
---|
210 | |
---|
211 | ; change cmd.var if macro |
---|
212 | CASE STRMID(cmd.var, 0, 2) OF |
---|
213 | '@@': BEGIN |
---|
214 | cmd.var = field.name |
---|
215 | cmd.var = STRMID(cmd.var, 2, 15) |
---|
216 | IF debug_w THEN print, 'cmd.var : ', cmd.var |
---|
217 | END |
---|
218 | ELSE: |
---|
219 | ENDCASE |
---|
220 | IF vecplot EQ 1 THEN BEGIN ; vectors case |
---|
221 | CASE STRMID(cmd.var2, 0, 2) OF |
---|
222 | '@@': begin |
---|
223 | cmd.var2 = STRMID(cmd.var2, 2, 15) |
---|
224 | END |
---|
225 | ELSE: |
---|
226 | ENDCASE |
---|
227 | ENDIF |
---|
228 | |
---|
229 | ; spectum plot |
---|
230 | IF spectplt EQ 1 THEN plttyp = 'spec' |
---|
231 | ; wavelet plot |
---|
232 | IF waveplt EQ 1 THEN plttyp = 'wavelet' |
---|
233 | |
---|
234 | xlogax = '' |
---|
235 | |
---|
236 | ; density pojection plot (splot=1) |
---|
237 | IF splot EQ 1 THEN BEGIN |
---|
238 | ; keep grid attributes |
---|
239 | izminmesh_b = izminmesh |
---|
240 | izmaxmesh_b = izmaxmesh |
---|
241 | jpk_b = jpk |
---|
242 | jpkglo_b = jpkglo |
---|
243 | gdept_b = gdept |
---|
244 | gdepw_b = gdepw |
---|
245 | e3t_b = e3t |
---|
246 | e3w_b = e3w |
---|
247 | tmask_b = tmask |
---|
248 | ; IF cmd.var EQ 'vosigvol' THEN xlogax = ',xlog = 1' |
---|
249 | ; modify vertical grid -> sigma |
---|
250 | n_sig = (sig_max - sig_min)/sig_del + 1 |
---|
251 | z_sig = sig_min+findgen(n_sig)*sig_del |
---|
252 | izminmesh = 0 |
---|
253 | izmaxmesh = n_sig-1 |
---|
254 | jpk = long(izmaxmesh-izminmesh+1) |
---|
255 | jpkglo = jpk |
---|
256 | gdept = z_sig |
---|
257 | gdepw = gdept |
---|
258 | e3t = shift(gdept, 1)-gdept |
---|
259 | e3t[0] = e3t[1] |
---|
260 | e3w = e3t |
---|
261 | prof1 = sig_min |
---|
262 | prof2 = sig_max |
---|
263 | ; grille,mask,glam,gphi,gdep,nx,ny,nz,premierx,premiery,premierz,dernierx,derniery,dernierz ; added PDW 12/5/04 |
---|
264 | ; tmask = mask ; added PDW 12/5/04 |
---|
265 | tmask = reform(reform(tmask[*, *, 0], jpi*jpj)#replicate(1, jpk), jpi, jpj, jpk) |
---|
266 | domdef |
---|
267 | END |
---|
268 | |
---|
269 | IF debug_w THEN print, ' cmd after data_read in plt_map = ', cmd |
---|
270 | ; |
---|
271 | ; ====================== |
---|
272 | ; 2. window attributes |
---|
273 | ; ====================== |
---|
274 | ; |
---|
275 | |
---|
276 | ; find field attributes |
---|
277 | ; --------------------- |
---|
278 | |
---|
279 | fldatt = {name:'', assos:'', min: 0.0, max: 0.0, int: 0.0, mult: 0.0, add: 0.0, unit: '', mid: 0.0, homin:0.0, homax:0.0, min1d:0.0, max1d:0.0, dmax:0.0} |
---|
280 | |
---|
281 | ; which variable to look for |
---|
282 | |
---|
283 | CASE plttyp OF |
---|
284 | 'ybinx': BEGIN |
---|
285 | ; if var3 exists then var = var1/var3 |
---|
286 | IF var3_ybinx NE "" THEN var_att = cmd.var+"/"+var3_ybinx ELSE var_att = cmd.var |
---|
287 | END |
---|
288 | ELSE: var_att = cmd.var |
---|
289 | ENDCASE |
---|
290 | |
---|
291 | IF debug_w THEN print, ' Looking for min/max/iso attributes for variable :', var_att |
---|
292 | |
---|
293 | ; extremum |
---|
294 | fldext = fld_pltext(var_att, plttyp, dimplot, hotyp) |
---|
295 | |
---|
296 | fldatt.name = fldext.name |
---|
297 | fldatt.min = fldext.min |
---|
298 | fldatt.max = fldext.max |
---|
299 | fldatt.homin = fldext.homin |
---|
300 | fldatt.homax = fldext.homax |
---|
301 | fldatt.min1d = fldext.min1d |
---|
302 | fldatt.max1d = fldext.max1d |
---|
303 | fldatt.dmax = fldext.dmax |
---|
304 | fldatt.assos = fldext.assos |
---|
305 | |
---|
306 | ; interval + change of unit |
---|
307 | fldint = fld_pltint(var_att, plttyp, dimplot, hotyp) |
---|
308 | |
---|
309 | ; help, fldint, /struct |
---|
310 | |
---|
311 | fldatt.int = fldint.int |
---|
312 | fldatt.mult = fldint.mult |
---|
313 | fldatt.add = fldint.add |
---|
314 | fldatt.unit = fldint.unit |
---|
315 | fldatt.mid = fldint.mid |
---|
316 | |
---|
317 | IF debug_w THEN print, 'plt_map 1 :', fldatt.mult, fldatt.add, fldatt.min, fldatt.max |
---|
318 | |
---|
319 | IF fldatt.mult NE 1. OR fldatt.add NE 0. THEN field.units = fldatt.unit |
---|
320 | fld = field.data*fldatt.mult + fldatt.add |
---|
321 | IF (where(field.data EQ valmask))[0] NE -1 THEN $ |
---|
322 | fld[where(field.data EQ valmask)] = valmask |
---|
323 | |
---|
324 | ; field 2 if needed |
---|
325 | |
---|
326 | IF plttyp EQ 'yfx' OR plttyp EQ 'ybinx' OR vecplot EQ 1 THEN BEGIN |
---|
327 | fldatt2 = fldatt |
---|
328 | |
---|
329 | fldext2 = fld_pltext(cmd.var2, plttyp, dimplot, hotyp) |
---|
330 | |
---|
331 | fldatt2.name = fldext2.name |
---|
332 | fldatt2.min = fldext2.min |
---|
333 | fldatt2.max = fldext2.max |
---|
334 | fldatt2.homin = fldext2.homin |
---|
335 | fldatt2.homax = fldext2.homax |
---|
336 | fldatt2.min1d = fldext2.min1d |
---|
337 | fldatt2.max1d = fldext2.max1d |
---|
338 | fldatt2.dmax = fldext2.dmax |
---|
339 | |
---|
340 | ; interval + change of unit |
---|
341 | fldint2 = fld_pltint(cmd.var2, plttyp, dimplot, hotyp) |
---|
342 | |
---|
343 | fldatt2.int = fldint2.int |
---|
344 | fldatt2.mult = fldint2.mult |
---|
345 | fldatt2.add = fldint2.add |
---|
346 | fldatt2.unit = fldint2.unit |
---|
347 | fldatt2.mid = fldint2.mid |
---|
348 | |
---|
349 | IF fldatt2.mult NE 1. OR fldatt2.add NE 0. THEN field.units2 = fldatt2.unit |
---|
350 | fld2 = field.data2*fldatt2.mult + fldatt2.add |
---|
351 | IF (where(field.data2 EQ valmask))[0] NE -1 THEN $ |
---|
352 | fld2[where(field.data2 EQ valmask)] = valmask |
---|
353 | |
---|
354 | ENDIF |
---|
355 | |
---|
356 | ; field 3 if needed |
---|
357 | |
---|
358 | IF plttyp EQ 'ybinx' AND var3_ybinx NE "" THEN BEGIN |
---|
359 | fldatt3 = fldatt |
---|
360 | |
---|
361 | fldext3 = fld_pltext(var3_ybinx, plttyp, dimplot, hotyp) |
---|
362 | |
---|
363 | fldatt3.name = fldext3.name |
---|
364 | fldatt3.min = fldext3.min |
---|
365 | fldatt3.max = fldext3.max |
---|
366 | fldatt3.homin = fldext3.homin |
---|
367 | fldatt3.homax = fldext3.homax |
---|
368 | fldatt3.min1d = fldext3.min1d |
---|
369 | fldatt3.max1d = fldext3.max1d |
---|
370 | fldatt3.dmax = fldext3.dmax |
---|
371 | |
---|
372 | ; interval + change of unit |
---|
373 | fldint3 = fld_pltint(var3_ybinx, plttyp, dimplot, hotyp) |
---|
374 | |
---|
375 | fldatt3.int = fldint3.int |
---|
376 | fldatt3.mult = fldint3.mult |
---|
377 | fldatt3.add = fldint3.add |
---|
378 | fldatt3.unit = fldint3.unit |
---|
379 | fldatt3.mid = fldint3.mid |
---|
380 | |
---|
381 | IF fldatt3.mult NE 1. OR fldatt3.add NE 0. THEN field.units3 = fldatt3.unit |
---|
382 | fld3 = field.data3*fldatt3.mult + fldatt3.add |
---|
383 | IF (where(field.data3 EQ valmask))[0] NE -1 THEN $ |
---|
384 | fld3[where(field.data3 EQ valmask)] = valmask |
---|
385 | |
---|
386 | ENDIF |
---|
387 | |
---|
388 | ; If running std dev change min/max |
---|
389 | IF run_stddev GT 0 THEN BEGIN |
---|
390 | fldatt.homin = 0. |
---|
391 | fldatt.homax = fldext.dmax |
---|
392 | ENDIF |
---|
393 | |
---|
394 | ; legend text |
---|
395 | ; ----------- |
---|
396 | |
---|
397 | ; IF plttyp EQ "pltt" THEN BEGIN |
---|
398 | CASE strmid(cmd.trend, 0, 1) OF |
---|
399 | '1': trd_txt = ' trend' |
---|
400 | '2': trd_txt = ' drift' |
---|
401 | '3': trd_txt = ' inverse trend' |
---|
402 | '4': trd_txt = ' anomaly' |
---|
403 | '5': trd_txt = ' filter' |
---|
404 | '7': trd_txt = ' month sel' |
---|
405 | ELSE: trd_txt = '' |
---|
406 | ENDCASE |
---|
407 | ; ENDIF ELSE trd_txt = '' |
---|
408 | IF field_int EQ 1 AND dimplot EQ 1 THEN BEGIN |
---|
409 | trd_txt = trd_txt+' integral' |
---|
410 | ENDIF |
---|
411 | |
---|
412 | varlegend = field.legend+trd_txt+' ('+field.units+')' |
---|
413 | |
---|
414 | CASE plttyp OF |
---|
415 | 'yfx': varlegend2 = field.legend2+trd_txt+' ('+field.units2+')' |
---|
416 | 'ybinx': BEGIN |
---|
417 | IF strmid(cmd2.trend, 0, 1) EQ '4' THEN varlegend = field.legend2+trd_txt+' ('+field.units2+')' $ |
---|
418 | ELSE varlegend = field.legend2+' ('+field.units2+')' |
---|
419 | varlegend2 = field.legend+trd_txt+' ('+field.units+')' |
---|
420 | END |
---|
421 | ELSE: BEGIN |
---|
422 | IF run_stddev GT 0 THEN trd_txt = trd_txt+' running Std Dev ['+string(run_stddev, format = '(I3)')+']' |
---|
423 | END |
---|
424 | ENDCASE |
---|
425 | |
---|
426 | ; date text |
---|
427 | |
---|
428 | CASE strmid(cmd.timave, 0, 2) OF |
---|
429 | '1m': BEGIN |
---|
430 | mn = def_month(cmd.timave, cmd.date1) |
---|
431 | IF strmid(cmd.timave, 0, 3) EQ '1mm' THEN $ |
---|
432 | date_txt = mn+' ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')' $ |
---|
433 | ELSE date_txt = mn+' '+strmid(cmd.date1, 0, strlen(cmd.date1)-2) |
---|
434 | END |
---|
435 | '3m': BEGIN |
---|
436 | mn = def_month(cmd.timave, cmd.date1) |
---|
437 | IF strmid(cmd.timave, 0, 3) EQ '3mm' THEN $ |
---|
438 | date_txt = mn+' ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')' $ |
---|
439 | ELSE date_txt = mn+' '+strmid(cmd.date1, 0, strlen(cmd.date1)-2) |
---|
440 | END |
---|
441 | '6m': BEGIN |
---|
442 | mn = def_month(cmd.timave, cmd.date1) |
---|
443 | IF strmid(cmd.timave, 0, 3) EQ '6mm' THEN $ |
---|
444 | date_txt = mn+' ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')' $ |
---|
445 | ELSE date_txt = mn+' '+strmid(cmd.date1, 0, strlen(cmd.date1)-2) |
---|
446 | END |
---|
447 | ELSE:date_txt = cmd.timave+' '+cmd.date1 |
---|
448 | ENDCASE |
---|
449 | |
---|
450 | CASE plttyp OF |
---|
451 | 'pltt':BEGIN |
---|
452 | date_txt = cmd.timave |
---|
453 | IF strmid(cmd.timave, 1, 2) EQ 'mm' THEN $ |
---|
454 | date_txt = cmd.timave+' ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')' |
---|
455 | END |
---|
456 | 'yfx': BEGIN |
---|
457 | IF hotyp NE '-' THEN date_txt = cmd.timave+' '+cmd.date1+' '+cmd.spec |
---|
458 | END |
---|
459 | 'ybinx': BEGIN |
---|
460 | IF hotyp NE '-' THEN date_txt = cmd.timave+' '+cmd.date1+' '+cmd.spec |
---|
461 | END |
---|
462 | 'wavelet': BEGIN |
---|
463 | date_txt = cmd.timave+' Wavelet' |
---|
464 | IF strmid(cmd.timave, 1, 2) EQ 'mm' THEN $ |
---|
465 | date_txt = cmd.timave+' Wavelet ('+strmid(cmd.date1, 3, strlen(cmd.date1)-3)+')' |
---|
466 | END |
---|
467 | ELSE: |
---|
468 | ENDCASE |
---|
469 | |
---|
470 | ; min/max/int |
---|
471 | |
---|
472 | ; defined if fraction of x/y.range is added to plot domain |
---|
473 | |
---|
474 | ; IF free_1d_minmax EQ 'no' THEN fraction = 0. ELSE fraction = 1.0 |
---|
475 | fraction = 0. |
---|
476 | |
---|
477 | CASE dimplot OF |
---|
478 | 1: BEGIN |
---|
479 | CASE plttyp OF |
---|
480 | 'pltt': BEGIN & minc = fldatt.homin & maxc = fldatt.homax & END |
---|
481 | 'wavelet': BEGIN & minc = fldatt.homin & maxc = fldatt.homax & END |
---|
482 | 'plt1d': BEGIN & minc = fldatt.min1d & maxc = fldatt.max1d & END |
---|
483 | 'yfx': BEGIN ; y=f(x) scatter plot |
---|
484 | IF long(strmid(cmd.trend, 0, 1)) GT 0 THEN BEGIN |
---|
485 | IF run_stddev EQ 0 THEN BEGIN |
---|
486 | minc = -fldatt.dmax & maxc = fldatt.dmax |
---|
487 | minc2 = -fldatt2.dmax & maxc2 = fldatt2.dmax |
---|
488 | ENDIF ELSE BEGIN |
---|
489 | minc = 0. & maxc = fldatt.dmax |
---|
490 | minc2 = 0. & maxc2 = fldatt2.dmax |
---|
491 | ENDELSE |
---|
492 | ENDIF ELSE BEGIN |
---|
493 | minc = fldatt.min1d & maxc = fldatt.max1d |
---|
494 | minc2 = fldatt2.min1d & maxc2 = fldatt2.max1d |
---|
495 | ENDELSE |
---|
496 | END |
---|
497 | 'ybinx': BEGIN ; y=bin(x) binning plot |
---|
498 | IF long(strmid(cmd.trend, 0, 1)) GT 0 THEN BEGIN |
---|
499 | minc = -fldatt.dmax & maxc = fldatt.dmax |
---|
500 | minc2 = fldatt2.min1d & maxc2 = fldatt2.max1d |
---|
501 | ENDIF ELSE BEGIN |
---|
502 | minc = fldatt.min1d & maxc = fldatt.max1d |
---|
503 | minc2 = fldatt2.min1d & maxc2 = fldatt2.max1d |
---|
504 | ENDELSE |
---|
505 | IF debug_w THEN print, ' minc/maxc for 1d fld1/fld2', minc, maxc, minc2, maxc2 |
---|
506 | END |
---|
507 | ELSE: BEGIN & minc = 0. & maxc = 0. & END |
---|
508 | ENDCASE |
---|
509 | fldatt.int = !VALUES.F_NAN |
---|
510 | END |
---|
511 | 2: BEGIN |
---|
512 | CASE plttyp OF |
---|
513 | 'pltt': BEGIN & minc = fldatt.homin & maxc = fldatt.homax & END |
---|
514 | ELSE: BEGIN & minc = fldatt.min & maxc = fldatt.max & END |
---|
515 | ENDCASE |
---|
516 | END |
---|
517 | ENDCASE |
---|
518 | |
---|
519 | IF cmd.var NE fld_prev OR cmd.trend NE '0' THEN BEGIN |
---|
520 | ; field.origin EQ 'diff' OR |
---|
521 | print, '' |
---|
522 | CASE dimplot OF |
---|
523 | 2: BEGIN |
---|
524 | print, ' '+cmd.var, ' plot attributes ('+plttyp+') : min/max/int=', $ |
---|
525 | minc, maxc, fldatt.int |
---|
526 | IF ( fldatt.mult NE 1.0 OR fldatt.add NE 0.0 ) THEN $ |
---|
527 | print, ' --> Modify field using : ', 'mult/add=' , $ |
---|
528 | fldatt.mult, fldatt.add, ' to obtain : ', fldatt.unit |
---|
529 | END |
---|
530 | 1: BEGIN |
---|
531 | print, ' '+cmd.var, ' plot attributes ('+plttyp+') : min/max=',minc, maxc |
---|
532 | IF ( fldatt.mult NE 1.0 OR fldatt.add NE 0.0 ) THEN $ |
---|
533 | print, ' --> Modify field using : ', 'mult/add=' , $ |
---|
534 | fldatt.mult, fldatt.add, ' to obtain : ', fldatt.unit |
---|
535 | END |
---|
536 | 0: BEGIN |
---|
537 | print, ' '+cmd.var, ' plot attributes ('+plttyp+') : min/max=',minc, maxc |
---|
538 | IF ( fldatt.mult NE 1.0 OR fldatt.add NE 0.0 ) THEN $ |
---|
539 | print, ' --> Modify field using : ', 'mult/add=' , $ |
---|
540 | fldatt.mult, fldatt.add, ' to obtain : ', fldatt.unit |
---|
541 | print, ' '+cmd.var2, ' plot attributes ('+plttyp+') : min/max=',minc2, maxc2 |
---|
542 | IF ( fldatt2.mult NE 1.0 OR fldatt2.add NE 0.0 ) THEN $ |
---|
543 | print, ' --> Modify field using : ', 'mult/add=' , $ |
---|
544 | fldatt2.mult, fldatt2.add, ' to obtain : ', fldatt2.unit |
---|
545 | END |
---|
546 | ELSE: |
---|
547 | ENDCASE |
---|
548 | print, '' |
---|
549 | |
---|
550 | ENDIF |
---|
551 | IF ( fldatt.mult EQ -1.0) THEN BEGIN |
---|
552 | temp_m = minc |
---|
553 | minc = -maxc |
---|
554 | maxc = -temp_m |
---|
555 | ENDIF |
---|
556 | IF finite(fldatt.int) EQ 0 THEN BEGIN |
---|
557 | intcmd = '' |
---|
558 | colbarfor = '' |
---|
559 | fmt = '(f6.1)' |
---|
560 | ENDIF ELSE BEGIN |
---|
561 | intcmd = ',int = '+string(fldatt.int) |
---|
562 | IF fldatt.int LT 0.1 THEN BEGIN |
---|
563 | fmt = '(f6.2)' |
---|
564 | ENDIF ELSE IF fldatt.int LT 1 THEN BEGIN |
---|
565 | fmt = '(f5.1)' |
---|
566 | ENDIF ELSE BEGIN |
---|
567 | fmt = '(f5.0)' |
---|
568 | ENDELSE |
---|
569 | IF long(fldatt.int) NE 0 THEN BEGIN |
---|
570 | IF fldatt.int/long(fldatt.int) NE 1 THEN fmt = '(f5.1)' |
---|
571 | ENDIF |
---|
572 | colbarfor = ', format = '''+fmt+'''' |
---|
573 | ENDELSE |
---|
574 | |
---|
575 | ; if window > 1 or overlay > 1 noerase |
---|
576 | |
---|
577 | IF win[2] NE 1 OR iover GT 1 THEN BEGIN |
---|
578 | noerase = 1 |
---|
579 | ENDIF ELSE BEGIN |
---|
580 | noerase = 0 |
---|
581 | ENDELSE |
---|
582 | |
---|
583 | ; choose window number |
---|
584 | |
---|
585 | CASE multi_win OF |
---|
586 | 1: window_number = ', window='+string(iplot) |
---|
587 | ELSE: window_number = '' |
---|
588 | ENDCASE |
---|
589 | |
---|
590 | ; color label control |
---|
591 | |
---|
592 | labstr = '' |
---|
593 | IF iover EQ 1 THEN BEGIN |
---|
594 | nofill = 1-shading |
---|
595 | ENDIF ELSE BEGIN |
---|
596 | nofill = 1 |
---|
597 | ENDELSE |
---|
598 | |
---|
599 | nocoltxt = '' |
---|
600 | IF nofill EQ 1 THEN nocoltxt = ',NOCOULEUR=1' |
---|
601 | |
---|
602 | IF dimplot NE 2 THEN nofill = 1 |
---|
603 | |
---|
604 | nocolorbar = 0 |
---|
605 | IF col_palette EQ 'no' THEN nocolorbar = 1 |
---|
606 | IF pal_type EQ '2dom' THEN nocolorbar = 1 |
---|
607 | IF iover GT 1 THEN nocolorbar = 1 |
---|
608 | IF dimplot NE 2 THEN nocolorbar = 1 |
---|
609 | |
---|
610 | IF cmd.out NE 'cdf' THEN BEGIN |
---|
611 | readpal = 0 |
---|
612 | IF nofill EQ 0 AND dimplot GT 1 THEN readpal = 1 |
---|
613 | IF waveplt EQ 1 THEN readpal = 2 |
---|
614 | IF field.origin EQ 'div' AND plttyp NE 'pltt' THEN readpal = 3 |
---|
615 | ;; Necessary for correlation plots (overlay of 2D fields) |
---|
616 | IF dimplot GT 1 AND nover GT 1 AND iover GT 1 THEN readpal = 1 |
---|
617 | IF readpal GE 1 THEN BEGIN |
---|
618 | lec_pal_gmt, cmd.var, c_anot_str, fmt, found, readpal |
---|
619 | colbarfor = ', format = '''+fmt+'''' |
---|
620 | IF found EQ 1 THEN BEGIN |
---|
621 | labstr = ',label=3' |
---|
622 | ; if mincmd/naxcmd not defined then use one from palette |
---|
623 | ; if define take min of 2 mins and max of 2 maxs |
---|
624 | min_palette = min(levels_gmt) |
---|
625 | max_palette = max(levels_gmt) |
---|
626 | IF finite(minc) NE 0 THEN minc = min([minc, min_palette]) ELSE minc = min_palette |
---|
627 | IF finite(maxc) NE 0 THEN maxc = max([maxc, max_palette]) ELSE maxc = max_palette |
---|
628 | print, ' '+cmd.var, ' plot attributes ('+plttyp+') modified via color palette: min/max=', $ |
---|
629 | minc, maxc |
---|
630 | ENDIF |
---|
631 | ENDIF |
---|
632 | ENDIF |
---|
633 | colbar = colbarfor |
---|
634 | |
---|
635 | mincmd = '' |
---|
636 | maxcmd = '' |
---|
637 | |
---|
638 | IF finite(minc) NE 0 THEN BEGIN |
---|
639 | mincmd = ','+string(minc) |
---|
640 | maxcmd = ','+string(maxc) |
---|
641 | ENDIF |
---|
642 | |
---|
643 | ; contour annotation + other color bar controls |
---|
644 | |
---|
645 | IF n_elements(c_anot_str) EQ 0 THEN BEGIN |
---|
646 | c_anot = '' |
---|
647 | ENDIF ELSE BEGIN |
---|
648 | c_anot = ',c_annotation=c_anot_str' |
---|
649 | ; sampling of colorbar annotation |
---|
650 | IF col_palette EQ 'yes' THEN BEGIN |
---|
651 | sampling = win[0]*win[1] |
---|
652 | ; idx_cb = lindgen((ncont_gmt-2)/sampling)*sampling+1 |
---|
653 | idx_cb0 = lindgen(ncont_gmt-1)+1 |
---|
654 | idx_cb = lindgen(ncont_gmt-2)+2 |
---|
655 | idx_cb1 = lindgen(ncont_gmt-2)+1 |
---|
656 | colbar = colbar+',discret=coul_gmt[idx_cb1],div=n_elements(idx_cb0)-1,cb_label=levels_gmt[idx_cb0]' |
---|
657 | ; c_anot = ',c_annotation=levels_gmt[idx_cb0]' |
---|
658 | ENDIF |
---|
659 | ENDELSE |
---|
660 | ; |
---|
661 | ; vertical grid type |
---|
662 | ; |
---|
663 | IF debug_w THEN print, ' splot=', splot |
---|
664 | CASE mesh_type OF |
---|
665 | 'oce': BEGIN |
---|
666 | IF splot EQ 1 THEN BEGIN |
---|
667 | type_yz = ',type_yz =''sigma''' |
---|
668 | zoom_txt = ', zoom='+string(sig_max) |
---|
669 | ENDIF ELSE BEGIN |
---|
670 | type_yz = ',type_yz=''m''' |
---|
671 | zoom_txt = ', zoom='+string(zoom_z) |
---|
672 | ENDELSE |
---|
673 | END |
---|
674 | 'atm': BEGIN & type_yz = ',type_yz=''Pa''' & zoom_txt = ', zoom='+string(pres_max) & END |
---|
675 | ELSE:type_yz = '' |
---|
676 | ENDCASE |
---|
677 | |
---|
678 | ; continents color, thickness |
---|
679 | ; c_cont = 100 |
---|
680 | ; continent fill |
---|
681 | ; real continent drawn |
---|
682 | strcont ='' |
---|
683 | IF mesh_type NE 'oce' THEN BEGIN |
---|
684 | IF cont_fill EQ 0 THEN strcont = ',/CONT_NOFILL, nite=0' |
---|
685 | IF cont_real GE 1 THEN strcont = strcont+', REALCONT = '+string(strtrim(cont_real, 2))+', COAST_THICK = 2' |
---|
686 | ENDIF |
---|
687 | |
---|
688 | ; calendar type |
---|
689 | |
---|
690 | cal_typ = '' |
---|
691 | IF calendar_type GT 1 THEN cal_typ = ',ndayspm=calendar_type' |
---|
692 | |
---|
693 | ; titles options |
---|
694 | ; By default, SAXO puts a title and a subtitle for the first plot |
---|
695 | ; Then, if you overlay another field, title and subtitle are set to '' |
---|
696 | CASE title_type OF |
---|
697 | "T": tit_str = ',subtitle='''','+marge_option |
---|
698 | "S": tit_str = ',title='''','+marge_option |
---|
699 | "TS": tit_str = ','+marge_option |
---|
700 | "off":tit_str = ',title='''',subtitle='''','+marge_option |
---|
701 | ENDCASE |
---|
702 | |
---|
703 | ; fill_space option |
---|
704 | |
---|
705 | filltxt = '' |
---|
706 | IF fill_space EQ 1 THEN filltxt = ',/rempli' |
---|
707 | |
---|
708 | ; axis charsize option (x/ychartxt read from plt_def) |
---|
709 | |
---|
710 | ; common command string to plots |
---|
711 | |
---|
712 | com_strplt = ',petit = ['+string(win[0])+','+string(win[1])+','+string(win[2])+']'+filltxt+nocoltxt+', LANDSCAPE = '+string(landscape)+', NOCOLORBAR = '+string(nocolorbar)+', NOERASE = '+string(noerase)+look+labstr+c_anot+colbar+cal_typ+',cont_thick=2'+type_yz+window_number+tit_str+xlogax+', xcharsize='+string(xchartxt)+', ycharsize='+string(ychartxt) |
---|
713 | |
---|
714 | ; add contour_options not overlay |
---|
715 | |
---|
716 | IF iover EQ 1 THEN com_strplt = com_strplt+contour_options |
---|
717 | |
---|
718 | IF debug_w THEN print, ' Cmd before making output in plt_map = ', cmd |
---|
719 | |
---|
720 | IF debug_w THEN print, ' We are in plt_map: ' |
---|
721 | @zoom_indices_debug |
---|
722 | |
---|
723 | ; ============================ |
---|
724 | ; 3. make output (data/plot) |
---|
725 | ; ============================ |
---|
726 | ; |
---|
727 | |
---|
728 | CASE cmd.out OF |
---|
729 | 'data': write_data = long((iplot-1)*10+win[2]) ; write 1D data ascii in trends.pro |
---|
730 | 'tcdf': write_data = -long((iplot-1)*10+win[2]) ; write 1D data netcdf in trends.pro |
---|
731 | 'cdf': write_data = 0 |
---|
732 | ELSE: write_data = 0 |
---|
733 | ENDCASE |
---|
734 | |
---|
735 | ; make output |
---|
736 | |
---|
737 | CASE cmd.out OF |
---|
738 | 'tv': BEGIN ; $$$$$$$$$$$$$$$$$$$$$$ tvnplot $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ |
---|
739 | erase |
---|
740 | |
---|
741 | CASE (size(fld))[0] OF |
---|
742 | 2: fldtv = fld |
---|
743 | 3: BEGIN |
---|
744 | niveau = xquestion('3d data : which level ?', '1', /chkwidget) |
---|
745 | fldtv = fld[*, *, niveau] |
---|
746 | END |
---|
747 | 4: BEGIN |
---|
748 | niveau = xquestion('3d data : which level ?', '1', /chkwidget) |
---|
749 | timeplot = xquestion('4d data : which date ?', '1', /chkwidget) |
---|
750 | fldtv = fld[*, *, niveau, timeplot] |
---|
751 | END |
---|
752 | ELSE: |
---|
753 | ENDCASE |
---|
754 | |
---|
755 | CASE cmd.grid OF |
---|
756 | ; 'T': tvnplot, fldtv*tmask |
---|
757 | 'T': tvplus, fldtv, min = fldatt.min, max = fldatt.max |
---|
758 | 'U': tvplus, fldtv, min = fldatt.min, max = fldatt.max |
---|
759 | 'V': tvplus, fldtv, min = fldatt.min, max = fldatt.max |
---|
760 | 'W': tvplus, fldtv, min = fldatt.min, max = fldatt.max |
---|
761 | ELSE: tvplus, fldtv, min = fldatt.min, max = fldatt.max |
---|
762 | ENDCASE |
---|
763 | END |
---|
764 | |
---|
765 | 'cdf': BEGIN ; $$$$$$$$$$$$$$$$$$$$$$ netCDF output $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ |
---|
766 | |
---|
767 | ; build structures containing grids |
---|
768 | cdf_description = nc_build(cmd, field, field.direc, vargrid) |
---|
769 | ; build netCDF file (here if whole data) |
---|
770 | IF hotyp EQ '-' OR hotyp EQ 'xyt' THEN BEGIN |
---|
771 | fldcdf = {data:field.data, units:field.units, short_name:cmd.var, long_name:field.legend, missing_value:valmask, direc:field.direc} |
---|
772 | file_out = cmd.var+'_'+strtrim(string(FORMAT = '(I2.2)', iplot), 2)+'.nc' |
---|
773 | pltcmd ='nc_put, fldcdf, file_out, ncdf_db = hom_idl'+cdf_description |
---|
774 | printf, nulhis, strcompress(pltcmd) |
---|
775 | res = execute(pltcmd) |
---|
776 | ENDIF |
---|
777 | END |
---|
778 | 'tcdf': BEGIN ; $$$$$$$$$$$$$$$$$$$$$$ 1D netCDF output $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ |
---|
779 | mask_z, fld, cmd, boite_plt1d, dimplot |
---|
780 | fld = checkfield(fld, 'pltt', type = datyp.hotyp, boite = boite_plt1d) |
---|
781 | IF debug_w THEN print, 'size(fld)', size(fld) |
---|
782 | IF debug_w THEN print, 'boite_plt1d', boite_plt1d |
---|
783 | IF cmd.trend GT 0 THEN BEGIN |
---|
784 | fld = trends(fld, cmd.trend, datyp.hotyp) |
---|
785 | add_txt = trd_txt |
---|
786 | ENDIF ELSE print, 'You need to have a trend to make 1D netcdf output' |
---|
787 | END |
---|
788 | ELSE: BEGIN ; $$$$$$$$$$$$$$$$$$$$$$$$$$ make plot $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ |
---|
789 | |
---|
790 | CASE plttyp OF |
---|
791 | 'plt': BEGIN |
---|
792 | ; |
---|
793 | ; Map plot |
---|
794 | ; -------- |
---|
795 | ; |
---|
796 | print, ' '+cmd.var+' data to plot min and max : ', $ |
---|
797 | min(fld(where (fld NE valmask)), /nan), max(fld(where (fld NE valmask)), /nan) |
---|
798 | print, '' |
---|
799 | |
---|
800 | map_cmd = '' |
---|
801 | |
---|
802 | niveau = 1 |
---|
803 | |
---|
804 | IF vecplot EQ 0 THEN BEGIN |
---|
805 | ; eddy field ? (if so, remove zonal mean) |
---|
806 | IF strpos(cmd.plt, '@z') GT 1 THEN BEGIN |
---|
807 | fldzm = moyenne(fld, 'x', boite=[20,380,-90,90], /NAN) |
---|
808 | fldzm = replicate(1, (size(fld))(1))#fldzm |
---|
809 | fld = fld-fldzm |
---|
810 | filleg = ' Eddy ' |
---|
811 | ENDIF ELSE BEGIN |
---|
812 | filleg = '' |
---|
813 | ENDELSE |
---|
814 | |
---|
815 | IF hotyp EQ 'xyt' THEN date_txt = 'monthly' |
---|
816 | |
---|
817 | pltcmd = 'plt,{a:fld, d:'''+date_txt+''', n:'''+filleg+varlegend+' '+legbox+legsurf+''', u:'''+field.units+''', g:'''+vargrid+'''}'+mincmd+maxcmd+intcmd+com_strplt+strcont+map_cmd |
---|
818 | ;; For 2D plots with overlays (correlation and proba for instance) |
---|
819 | ;; Scratch the titles and define contours style for the second plot |
---|
820 | |
---|
821 | IF nover GT 1 AND iover GT 1 THEN BEGIN |
---|
822 | contour_style = ',c_thick=2,c_linestyle=2,c_charthick=2' |
---|
823 | CASE title_type OF |
---|
824 | 'T': pltcmd = pltcmd+',title='''''+contour_style |
---|
825 | 'S': pltcmd = pltcmd+',subtitle='''''+contour_style |
---|
826 | 'TS': pltcmd = pltcmd+',title='''',subtitle='''''+contour_style |
---|
827 | 'off': pltcmd = pltcmd+contour_style |
---|
828 | END |
---|
829 | ENDIF |
---|
830 | ENDIF ELSE BEGIN |
---|
831 | |
---|
832 | fld = {a:fld, g:strmid(cmd.grid, 0, 1)} |
---|
833 | fld2 = {a:fld2, g:strmid(cmd.grid, 1, 1)} |
---|
834 | fldv = {data1: fld, data2: fld2} |
---|
835 | vargrid = strmid(cmd.grid, 2, 1) |
---|
836 | |
---|
837 | pltcmd = 'ajoutvect,fldv'+vect_sample+com_strplt+strcont+map_cmd |
---|
838 | |
---|
839 | ENDELSE |
---|
840 | |
---|
841 | IF debug_w THEN print, strcompress(pltcmd) |
---|
842 | printf, nulhis, strcompress(pltcmd) |
---|
843 | res = execute(pltcmd(0)) |
---|
844 | |
---|
845 | ; draw nino boxes |
---|
846 | IF nino_plot EQ 1 THEN BEGIN |
---|
847 | xn4 = [160, 210, 210, 160, 160] |
---|
848 | yn4 = [5, 5, -5, -5, 5] |
---|
849 | plot, xn4, yn4, /noerase |
---|
850 | xn3 = [210, 270, 270, 210, 210] |
---|
851 | yn3 = [5, 5, -5, -5, 5] |
---|
852 | plot, xn3, yn3, /noerase |
---|
853 | xn1 = [270, 280, 280, 270, 270] |
---|
854 | yn1 = [-5, -5, -10, -10, -5] |
---|
855 | plot, xn1, yn1, /noerase |
---|
856 | xn2 = [270, 280, 280, 270, 270] |
---|
857 | yn2 = [0, 0, -5, -5, 0] |
---|
858 | plot, xn2, yn2, /noerase |
---|
859 | ENDIF |
---|
860 | |
---|
861 | END |
---|
862 | 'pltz': BEGIN |
---|
863 | ; |
---|
864 | ; Vertical section/mean plot |
---|
865 | ; -------------------------- |
---|
866 | ; |
---|
867 | g = gphit |
---|
868 | t = tmask |
---|
869 | |
---|
870 | ; if already 1D, reform fld |
---|
871 | sz = size(fld) |
---|
872 | IF sz[0] EQ 2 THEN BEGIN |
---|
873 | IF sz[1] EQ 1 OR sz[2] EQ 1 THEN fld = reform(fld, sz[1]*sz[2]) |
---|
874 | ENDIF |
---|
875 | IF sz[0] EQ 3 THEN BEGIN |
---|
876 | IF sz[1] EQ 1 THEN fld = reform(fld, sz[2], sz[3]) |
---|
877 | IF sz[2] EQ 1 THEN fld = reform(fld, sz[1], sz[3]) |
---|
878 | IF sz[3] EQ 1 THEN fld = reform(fld, sz[1], sz[2]) |
---|
879 | ENDIF |
---|
880 | |
---|
881 | ; mask fld |
---|
882 | mask_z, fld, cmd, boite_pltz, dimplot, legz |
---|
883 | printf, nulhis, ' boite_pltz = ', boite_pltz |
---|
884 | |
---|
885 | print, ' '+cmd.var+' data to plot min and max : ', $ |
---|
886 | min(fld(where (fld NE valmask)), /NAN), max(fld(where (fld NE valmask)), /NAN) |
---|
887 | print, '' |
---|
888 | IF vecplot EQ 0 THEN BEGIN |
---|
889 | |
---|
890 | IF cmd.var EQ 'vozonbsf' THEN xindstr = ', /xindex' ELSE xindstr = '' |
---|
891 | |
---|
892 | pltcmd = 'pltz,{a:fld, d:'''+date_txt+''', n:'''+varlegend+' '+legz+''', u:'''+field.units+''', g:'''+vargrid+'''}'+mincmd+maxcmd+intcmd+',/'+pltztyp+', boite=boite_pltz'+zoom_txt+com_strplt+xindstr |
---|
893 | |
---|
894 | |
---|
895 | ENDIF ELSE BEGIN |
---|
896 | |
---|
897 | fld = {a:fld, g:strmid(cmd.grid, 0, 1)} |
---|
898 | fld2 = {a:fld2, g:strmid(cmd.grid, 1, 1)} |
---|
899 | fldv = {data1: fld, data2: fld2} |
---|
900 | vargrid = strmid(cmd.grid, 2, 1) |
---|
901 | |
---|
902 | pltcmd = 'ajoutvectz,fldv'+com_strplt+strcont+',type='''+strmid(cmd.plt, 0, 2)+''', boite=boite_pltz' |
---|
903 | |
---|
904 | ENDELSE |
---|
905 | |
---|
906 | IF debug_w THEN print, pltcmd |
---|
907 | printf, nulhis, strcompress(pltcmd) |
---|
908 | res = execute(pltcmd(0)) |
---|
909 | |
---|
910 | |
---|
911 | ; overlay bowl if sig_bowl=1 |
---|
912 | IF sig_bowl EQ 1 THEN begin |
---|
913 | ; define line color, thickness and type |
---|
914 | iover = 2 |
---|
915 | overc = overlay_type(iover, dimplot) |
---|
916 | ; type of latitude axis |
---|
917 | IF strmid(cmd.plt, 0, 1) EQ 'y' AND lat_axis EQ 'sin' THEN $ |
---|
918 | sin = ',/sin' ELSE sin = '' |
---|
919 | plt1dtyp = strmid(pltztyp, 0, 1) |
---|
920 | noerase = 1 |
---|
921 | com_strplt = ', NOERASE = '+string(noerase) |
---|
922 | pltcmd = 'plt1d,field.bowl,/'+plt1dtyp+', boite=boite_pltz'+overc+sin+com_strplt |
---|
923 | IF debug_w THEN print, pltcmd |
---|
924 | |
---|
925 | res = execute(pltcmd(0)) |
---|
926 | ENDIF |
---|
927 | |
---|
928 | gphit = g |
---|
929 | tmask = t |
---|
930 | |
---|
931 | END |
---|
932 | 'plt1d': BEGIN |
---|
933 | ; |
---|
934 | ; 1D, x, y, z plot |
---|
935 | ; ---------------- |
---|
936 | ; |
---|
937 | g = gphit |
---|
938 | t = tmask |
---|
939 | ; if already 1D, reform fld |
---|
940 | sz = size(fld) |
---|
941 | IF sz[0] EQ 2 THEN BEGIN |
---|
942 | IF sz[1] EQ 1 OR sz[2] EQ 1 THEN fld = reform(fld, sz[1]*sz[2]) |
---|
943 | ENDIF |
---|
944 | IF sz[0] EQ 3 THEN BEGIN |
---|
945 | IF sz[1] EQ 1 THEN fld = reform(fld, sz[2], sz[3]) |
---|
946 | IF sz[2] EQ 1 THEN fld = reform(fld, sz[1], sz[3]) |
---|
947 | IF sz[3] EQ 1 THEN fld = reform(fld, sz[1], sz[2]) |
---|
948 | ENDIF |
---|
949 | ; mask fld |
---|
950 | mask_z, fld, cmd, boite_plt1d, dimplot, legz |
---|
951 | |
---|
952 | IF debug_w THEN print, boite_plt1d |
---|
953 | niveau = 1 |
---|
954 | ; define line color, thickness and type |
---|
955 | overc = overlay_type(iover, dimplot) |
---|
956 | ; type of latitude axis |
---|
957 | IF strmid(cmd.plt, 0, 1) EQ 'y' AND lat_axis EQ 'sin' THEN $ |
---|
958 | sin = ',/sin' ELSE sin = '' |
---|
959 | |
---|
960 | print, ' '+cmd.var+' data to plot min and max : ', $ |
---|
961 | min(fld(where (fld NE valmask))), max(fld(where (fld NE valmask))) |
---|
962 | print, '' |
---|
963 | take_log = ',take_log=0' |
---|
964 | IF cmd.var EQ 'vosigvol' AND n_elements(strsplit(cmd.exp,'-', /EXTRACT)) EQ 1 THEN BEGIN ; don't take log if difference plot |
---|
965 | take_log = ',take_log=1' |
---|
966 | ENDIF |
---|
967 | pltcmd = 'plt1d,{a:fld, d:'''+date_txt+''', n:'''+varlegend+' '+legz+''', u:'''+field.units+''', g:'''+vargrid+'''},'''+plt1dtyp+''''+mincmd+maxcmd+intcmd+', boite=boite_plt1d'+overc+sin+com_strplt+take_log |
---|
968 | printf, nulhis, 'boite_plt1d=',boite_plt1d |
---|
969 | printf, nulhis, strcompress(pltcmd) |
---|
970 | IF debug_w THEN print, pltcmd |
---|
971 | res = execute(pltcmd(0)) |
---|
972 | |
---|
973 | @legend_overlay |
---|
974 | |
---|
975 | gphit = g |
---|
976 | tmask = t |
---|
977 | |
---|
978 | END |
---|
979 | 'pltt': BEGIN |
---|
980 | ; |
---|
981 | ; hovmoeller |
---|
982 | ; ----------- |
---|
983 | ; |
---|
984 | g = gphit |
---|
985 | t = tmask |
---|
986 | |
---|
987 | ; mask fld |
---|
988 | mask_z, fld, cmd, boite_pltt, dimplot, legz |
---|
989 | ; define line color, thickness and type |
---|
990 | overc = overlay_type(iover, dimplot) |
---|
991 | ; define additional text for pltt |
---|
992 | @add_txt_pltt |
---|
993 | ; repeat for seasonal mean |
---|
994 | CASE cmd.timave OF |
---|
995 | '1mm': rep_txt = ',repeat_c=nb_cycles,DATE_FORMAT="%M"' |
---|
996 | ELSE: rep_txt = '' |
---|
997 | ENDCASE |
---|
998 | |
---|
999 | print, ' '+cmd.var+' data to plot min and max : ', $ |
---|
1000 | min(fld(where (fld NE valmask))), max(fld(where (fld NE valmask))) |
---|
1001 | print, '' |
---|
1002 | vardate = 'toto' |
---|
1003 | grille, mask, glam, gphi, gdep, nx, ny,nz |
---|
1004 | IF st_rms EQ 1 THEN BEGIN ; time series of rms deviation on a sigma surface |
---|
1005 | @densit_pltmap_plot |
---|
1006 | ENDIF ELSE BEGIN |
---|
1007 | pltcmd = 'pltt,{a:fld, n:'''+date_txt+' '+varlegend+' '+legbox+''', u:'''+field.units+filleg+''', g:'''+vargrid+'''}, '''+hotyp+''''+mincmd+maxcmd+intcmd+',boite=boite_pltt'+com_strplt+overc+filtxt+',TREND_TYPE='+cmd.trend+rep_txt ; +',/integration' |
---|
1008 | ENDELSE |
---|
1009 | |
---|
1010 | printf, nulhis, 'boite_pltt=',boite_pltt |
---|
1011 | printf, nulhis, strcompress(pltcmd) |
---|
1012 | IF debug_w THEN print, pltcmd |
---|
1013 | res = execute(pltcmd(0)) |
---|
1014 | IF iover EQ 4 THEN print, 'There might be a pb with the legends !' |
---|
1015 | ; legend if plot=t |
---|
1016 | IF hotyp EQ 't' THEN BEGIN |
---|
1017 | ; positions of the legend depend on nb_cycles |
---|
1018 | ; it is different from 1 only with 1mm case |
---|
1019 | IF cmd.timave NE '1mm' THEN nb_cycles = 1 |
---|
1020 | @legend_overlay |
---|
1021 | ENDIF |
---|
1022 | gphit = g |
---|
1023 | tmask = t |
---|
1024 | END |
---|
1025 | |
---|
1026 | 'yfx': BEGIN |
---|
1027 | @yfx |
---|
1028 | END |
---|
1029 | 'ybinx': BEGIN |
---|
1030 | @ybinx |
---|
1031 | END |
---|
1032 | 'spec': BEGIN |
---|
1033 | ; |
---|
1034 | ; Spectrum plot y=fft(x) |
---|
1035 | ; --------------------- |
---|
1036 | ; |
---|
1037 | t = tmask |
---|
1038 | ; plot domain |
---|
1039 | mask_z, fld, cmd, boite_pltspec, dimplot, legspec |
---|
1040 | ; define line color and type |
---|
1041 | overc = overlay_type(iover, dimplot) |
---|
1042 | vardate = date_txt |
---|
1043 | varname = varlegend |
---|
1044 | varunit = cmd.exp+' '+cmd.timave+' '+cmd.date1+'-'+cmd.spec+' '+varlegend+' '+legspec |
---|
1045 | |
---|
1046 | ; make mean in box |
---|
1047 | fld = checkfield(fld, 'pltt', type = 't', boite = boite_pltspec, _extra = ex) |
---|
1048 | ; apply trends |
---|
1049 | IF long(cmd.trend) GT 0 THEN fld = trends(fld, long(cmd.trend), 't') |
---|
1050 | |
---|
1051 | ; make spectrum |
---|
1052 | time_interval = time[1]-time[0] |
---|
1053 | |
---|
1054 | print, ' Max plot spectrum in plt_def : days/month/year ', max_spec, max_spec/30, max_spec/360 |
---|
1055 | print, ' lenght of time serie (years) = ', long(time(jpt-1)-time(0)+time_interval)/360 |
---|
1056 | print, ' spectrum window (years)/chunks = ', spec_win/360, long(time(jpt-1)-time(0)+time_interval)/spec_win |
---|
1057 | print, ' ' |
---|
1058 | |
---|
1059 | ; number of time windowsš |
---|
1060 | nt_win = long((time(jpt-1)-time(0)+time_interval)/spec_win) |
---|
1061 | idx_win_size = long(spec_win/time_interval) |
---|
1062 | IF idx_win_size*nt_win NE jpt THEN BEGIN |
---|
1063 | print, ' *** Warning : spec_win must divide lenght of time serie ', idx_win_size, jpt |
---|
1064 | return |
---|
1065 | ENDIF |
---|
1066 | ; mean of idx_win_size chunks |
---|
1067 | spect = spectrum(fld[0:idx_win_size-1], time_interval) |
---|
1068 | FOR it = 2, nt_win DO BEGIN |
---|
1069 | idx1 = idx_win_size*(it-1) |
---|
1070 | idx2 = idx1 + idx_win_size - 1 |
---|
1071 | specc = spectrum(fld[idx1:idx2], time_interval) |
---|
1072 | spect = spect + specc |
---|
1073 | ENDFOR |
---|
1074 | spect = spect/nt_win |
---|
1075 | ; build new time array |
---|
1076 | idx = where (spect[0, *] LE max_spec) |
---|
1077 | jpt = n_elements(idx) |
---|
1078 | fld = reverse(reform(spect[1, idx], jpt)) |
---|
1079 | time = findgen(jpt) |
---|
1080 | time = reverse(reform(spect[0, idx], jpt)) |
---|
1081 | IF max(time) GT 20*360 THEN BEGIN |
---|
1082 | time = time/360 |
---|
1083 | xtitle = 'Years' |
---|
1084 | ENDIF ELSE IF max(time) GT 5*360 THEN BEGIN |
---|
1085 | time = time/30 |
---|
1086 | xtitle = 'Months' |
---|
1087 | ENDIF ELSE xtitle = 'Days' |
---|
1088 | ; time/space filter ? |
---|
1089 | IF strpos(cmd.plt, '@f') GT 1 THEN BEGIN |
---|
1090 | filter = long(strmid(cmd.plt, strpos(cmd.plt, '@f')+3, strlen(cmd.plt)-strpos(cmd.plt, '@f')-3)) |
---|
1091 | fildirec = strmid(cmd.plt, strpos(cmd.plt, '@f')+2, 1) |
---|
1092 | IF fildirec EQ 't' THEN BEGIN |
---|
1093 | xtitle = xtitle+' (filter '+strtrim(string(filter), 2)+' warning : discrete filter)' |
---|
1094 | fld = smooth(fld, filter) |
---|
1095 | fld[0:filter/2-1] = 0. |
---|
1096 | fld[(size(fld))[1]-filter/2-1:(size(fld))[1]-1] = 0. |
---|
1097 | ENDIF |
---|
1098 | ENDIF |
---|
1099 | ; plot min/max |
---|
1100 | min1 = min(fld) |
---|
1101 | max1 = max(fld) |
---|
1102 | min2 = 0 |
---|
1103 | max2 = max(time) |
---|
1104 | !x.range = [min2-abs(max2-min2)/5.,max2+abs(max2-min2)/5.] |
---|
1105 | !y.range = [min1-abs(max1-min1)/5.,max1+abs(max1-min1)/5.] |
---|
1106 | |
---|
1107 | ; plot |
---|
1108 | ytitle = 'Power spectrum (window='+strtrim(string(spec_win/360), 2)+'y)' |
---|
1109 | pltcmd = 'splot,time,fld,xstyle=1,ystyle=1,title=varunit,xtitle=xtitle,ytitle=ytitle'+overc+com_strplt |
---|
1110 | printf, nulhis, 'boite_pltspec=',boite_pltspec |
---|
1111 | printf, nulhis, strcompress(pltcmd) |
---|
1112 | IF debug_w THEN print, pltcmd |
---|
1113 | res = execute(pltcmd(0)) |
---|
1114 | |
---|
1115 | tmask = t |
---|
1116 | |
---|
1117 | END |
---|
1118 | 'wavelet': BEGIN |
---|
1119 | ; |
---|
1120 | ; Wavelet plot |
---|
1121 | ; ------------ |
---|
1122 | ; |
---|
1123 | t = tmask |
---|
1124 | ; plot domain |
---|
1125 | mask_z, fld, cmd, boite_pltspec, dimplot, legspec |
---|
1126 | ; define line color and type |
---|
1127 | vardate = date_txt |
---|
1128 | varname = varlegend + ' '+date_txt |
---|
1129 | varunit = cmd.exp+' '+cmd.timave+' '+cmd.date1+'-'+cmd.spec+' '+varlegend+' '+legspec |
---|
1130 | ; make mean in box |
---|
1131 | fld = checkfield(fld, 'pltt', type = 't', boite = boite_pltspec, _extra = ex) |
---|
1132 | ; apply trends |
---|
1133 | IF long(cmd.trend) GT 0 THEN fld = trends(fld, long(cmd.trend), 't') |
---|
1134 | |
---|
1135 | ; make spectrum |
---|
1136 | time_interval = time[1]-time[0] |
---|
1137 | |
---|
1138 | print, ' Max plot spectrum in plt_def : days/month/year ', max_spec, max_spec/30, max_spec/360 |
---|
1139 | print, ' lenght of time serie (years) = ', long(time(jpt-1)-time(0)+time_interval)/360 |
---|
1140 | print, ' ' |
---|
1141 | ; make wavelet |
---|
1142 | wave = wavelet(fld,time_interval,period=period,coi=coi,/pad,signif=signif) |
---|
1143 | power = abs(wave)^2 |
---|
1144 | nscale = n_elements(period) |
---|
1145 | period = period/365 ; to have period in years |
---|
1146 | ; compute significance |
---|
1147 | signif = rebin(transpose(signif),jpt,nscale) |
---|
1148 | ; plot wavelet |
---|
1149 | mincmd = ','+string(min(power)) |
---|
1150 | maxcmd = ','+string(max(power)) |
---|
1151 | boite_pltspec = [0, 0, min(period), max(period)] |
---|
1152 | ; domdef, boite_pltspec |
---|
1153 | lat1r = lat1 |
---|
1154 | lat2r = lat2 |
---|
1155 | lat1 = 0 |
---|
1156 | lat2 = max_spec/360 |
---|
1157 | key_onearth = 0 |
---|
1158 | pltcmd = 'plttg,power,period'+mincmd+maxcmd+intcmd+',boite=boite_pltspec,reverse_y=1,nocontour=1'+com_strplt |
---|
1159 | printf, nulhis, 'boite_pltspec=',boite_pltspec |
---|
1160 | printf, nulhis, strcompress(pltcmd) |
---|
1161 | IF debug_w THEN print, pltcmd |
---|
1162 | res = execute(pltcmd(0)) |
---|
1163 | contour,abs(wave)^2/signif,time,period, /overplot,level=1.0,c_annot='95%' |
---|
1164 | plot,time,coi/365, noclip = 0, /noerase |
---|
1165 | tmask = t |
---|
1166 | lat1 = lat1r |
---|
1167 | lat2 = lat2r |
---|
1168 | key_onearth = 1 |
---|
1169 | END |
---|
1170 | ELSE: BEGIN |
---|
1171 | print, ' unknown projection plot ', cmd.plt, ' / plot type = ', plttyp |
---|
1172 | |
---|
1173 | END |
---|
1174 | ENDCASE |
---|
1175 | |
---|
1176 | END |
---|
1177 | |
---|
1178 | ENDCASE |
---|
1179 | |
---|
1180 | fld_prev = cmd.var |
---|
1181 | ; |
---|
1182 | ; reset incompatible options of plt_def |
---|
1183 | |
---|
1184 | cmd.trend = trend_typp |
---|
1185 | |
---|
1186 | ; |
---|
1187 | ; reset vertical grid after density bining |
---|
1188 | |
---|
1189 | IF splot EQ 1 THEN BEGIN |
---|
1190 | |
---|
1191 | izminmesh = izminmesh_b |
---|
1192 | izmaxmesh = izmaxmesh_b |
---|
1193 | jpk = jpk_b |
---|
1194 | jpkglo = jpkglo_b |
---|
1195 | gdept = gdept_b |
---|
1196 | gdepw = gdepw_b |
---|
1197 | e3t = e3t_b |
---|
1198 | e3w = e3w_b |
---|
1199 | tmask = tmask_b |
---|
1200 | |
---|
1201 | ENDIF |
---|
1202 | IF debug_w THEN print, ' ...EXIT plt_map' |
---|
1203 | |
---|
1204 | END |
---|