source: trunk/MESOSCALE_DEV/PLOT/MINIMAL/section.pro @ 357

Last change on this file since 357 was 85, checked in by aslmd, 14 years ago

LMD_MM_MARS et LMD_LES_MARS: ajout des routines IDL pour tracer les sorties --> voir mesoscale/PLOT

File size: 11.0 KB
Line 
1pro section, $
2        what_I_plot, $                          ; 2D field
3        space, $                                ; horizontal coordinate
4        altitude, $                             ; altitude coordinate
5        minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
6        maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
7        minspace=minspace, $                    ; minimum value of space window (=0: calculate)
8        maxspace=maxspace, $                    ; maximum value of space window (=0: calculate)
9        overcontour=overcontour, $              ; another 2D field to overplot with contour lines (=0: no)
10        overvector_x=overvector_x, $            ; wind vector - x component (=0: no)
11        overvector_y=overvector_y, $            ; wind vector - y component (=0: no)
12        colors=colors, $                        ; number of colors/levels (32 is default)
13        title_plot=title_user, $                ; title of the plot ('Profile' is default)
14        title_axis=title_axis, $                ; title of the [x,y] axis (['Field','Altitude'] is default)
15        ct=pal, $                               ; color table (33-rainbow is default)
16topo=topography, $
17        format=format                           ; format of colorbar annotations ('(F6.2)' is default)
18
19       
20;--------------------------------------------
21;
22; A general routine to plot H/V section maps
23;
24; USE:  - section, what_I_plot 
25;       - section, what_I_plot, space, altitude
26;       - ... and see options above
27;
28;--------------------------------------------
29; A. Spiga, September 2007
30;--------------------------------------------
31
32
33;---------------------------
34; a few permanent settings
35;---------------------------
36;
37missing_value=1.e30
38flag_cb='true'
39SPAWN, 'touch param_plot.idl'
40space_save=space
41@param_plot.idl
42;
43;------------------
44
45
46
47;------------------
48; user settings
49;------------------
50if (n_elements(title_user) eq 0) then title='Section' else title=title_user
51if (n_elements(subtitle_user) eq 0) then subtitle_user=''
52if (n_elements(title_axis) ne 2) then begin
53        xtitle='Horizontal coordinate'
54        ytitle='Altitude'
55endif else begin
56        xtitle=title_axis(0)
57        ytitle=title_axis(1)
58endelse
59
60if (n_elements(space) eq 0) then begin
61        space=findgen(n_elements(what_I_plot(*,0)))
62endif
63if (n_elements(altitude) eq 0) then begin
64        altitude=findgen(n_elements(what_I_plot(0,*)))
65endif
66
67if (n_elements(minfield_init) eq 0) then minfield_init=0.
68if (n_elements(maxfield_init) eq 0) then maxfield_init=0.
69if (n_elements(minspace) eq 0) then minspace=0.
70if (n_elements(maxspace) eq 0) then maxspace=0.
71if (n_elements(minalt) eq 0) then minalt=min(altitude)
72if (n_elements(maxalt) eq 0) then maxalt=max(altitude)
73if (n_elements(colors) eq 0) then colors=32
74if (n_elements(title) eq 0) then title=''
75if (n_elements(pal) eq 0) then pal=33
76if (n_elements(format) eq 0) then format='(F6.2)'
77if (format eq '') then format='(F6.2)'
78
79if ((n_elements(overvector_x) eq 0) or (n_elements(overvector_y) eq 0)) then begin
80        overvector_x=0.
81        overvector_y=0.
82endif
83if (n_elements(overcontour) eq 0) then overcontour=0.
84
85
86;
87; dilatation factor (for winds mainly)
88;
89;if (n_elements(overvector_x)*n_elements(overvector_y) ne 1) then begin
90
91if (n_elements(spacekm) eq 0) then spacekm='false'
92
93if (n_elements(factor) ne 0.) then begin
94        if (factor eq 0.) then stop
95        space=space/factor
96        if (spacekm ne 'false') then xtitle=xtitle+' (km)'
97        xtitle=xtitle+' - divided by '+string(factor,'(I0)')
98endif else begin
99        factor=1.
100endelse
101;endif
102       
103;------------------
104; limits
105;------------------
106;if (minfield_init*maxfield_init eq 0) then begin   
107if (minfield_init eq maxfield_init) then begin
108        ;---different min/max for each layer
109        w=where(abs(what_I_plot) lt missing_value)
110        if (w(0) ne -1) then begin
111                minfield=min(what_I_plot[w])
112                maxfield=max(what_I_plot[w])
113        endif else begin
114                what_I_plot=0.
115                print, 'no valid values in this field'
116                stop
117        endelse
118endif else begin
119        ;---user-defined limits
120        minfield=minfield_init
121        maxfield=maxfield_init
122        ;;w=where(what_I_plot lt minfield) & if (w(0) ne -1) then what_I_plot[w]=minfield
123        ;;w=where(what_I_plot gt maxfield) & if (w(0) ne -1) then what_I_plot[w]=maxfield
124endelse
125
126if ((minfield+maxfield eq 0) and (abs(minfield) ne abs(maxfield))) then begin
127        print, 'nothing to plot ... skipping this plot ...'
128endif else begin
129
130
131
132;-----------------------------------------
133; fix the possible longitude ugliness :)
134; ...get rid of the -180/180 limit
135;-----------------------------------------
136
137if (spacekm ne 'true') then begin
138        w=where(space eq min(space))
139        if (w(0) ne 0) then begin
140                space[0:w(0)-1]=space[0:w(0)-1]-360
141        endif
142endif   
143
144
145;---------------
146; adapt ticks
147;---------------
148if ((minspace eq 0) and (maxspace eq 0)) then begin
149        minspace=min(space)
150        maxspace=max(space)
151endif
152if (n_elements(intervalx) eq 0) then intervalx=round((maxspace-minspace)/8.)
153if (n_elements(intervaly) eq 0) then intervaly=round((maxalt-minalt)/5.)
154
155
156;------------------
157; plot window
158;------------------
159
160        ; to ensure the right limits
161        dumb_what_I_plot=what_I_plot
162        dumb_what_I_plot(*,*)=minfield & print, minfield
163        dumb_what_I_plot(0,0)=maxfield & print, maxfield
164
165
166;xticks=sindgen(21)
167;xticks(*)=' '
168;;xticks(0)='-120'
169;;xticks(2)='-119'
170;;xticks(4)='-118'
171;;xticks(6)='-117'
172;;xticks(8)='-116'
173;;xticks(10)='-115'
174;;xticks(12)='-114'
175;;xticks(14)='-113'
176;;xticks(16)='-112'
177;;xticks(18)='-111'
178;;xticks(20)='-110'
179;xticks(0)='-110'
180;xticks(2)='-109'
181;xticks(4)='-108'
182;xticks(6)='-107'
183;xticks(8)='-106'
184;xticks(10)='-105'
185;xticks(12)='-104'
186;xticks(14)='-103'
187;xticks(16)='-102'
188;xticks(18)='-101'
189;xticks(20)='-100'
190;print, xticks
191
192loadct,0,/silent 
193contour, $
194        /nodata, $
195        dumb_what_I_plot, $
196        space, $
197        altitude, $
198        title=title, $
199        xtitle=xtitle, $
200        xrange=[minspace,maxspace], $
201        ytitle=ytitle, $
202        yrange=[minalt,maxalt], $
203        xtickinterval=intervalx, $
204        ytickinterval=intervaly, $
205;       position=[0.15, 0.15, 0.95, 0.75], $ ;;OLD OLD
206        position=[0.15, 0.35, 0.95, 0.95], $
207        subtitle=subtitle_user, $
208;/isotropic, $
209xtickname=xticks, $
210        color=0
211
212;if (n_elements(topography) gt 1) then begin
213;oplot, space, topography/1000.
214;print, max(topography)
215;endif
216
217;------------------
218; plot field
219;------------------
220
221; to avoid spurious blank zones
222nlevels=colors
223levp=minfield+(findgen(nlevels)/float(nlevels-1))*(maxfield-minfield)
224
225loadct,pal,/silent
226;reverse_ct
227contour, what_I_plot, $
228        space, $
229        altitude, $
230;       nlevels=colors, $
231        levels=levp, $
232        /cell_fill, $
233        max_value=maxfield, $
234        min_value=minfield, $   
235        /overplot
236
237
238;------------------
239; colorbar
240;------------------
241
242;if (flag_cb eq 'true') then begin
243;        format_colorbar=format
244;        colorbar, $
245;                position=[0.15, 0.85, 0.95, 0.90], $
246;                divisions=8, $
247;                range=[float(minfield),float(maxfield)], $
248;                format=format_colorbar
249;endif
250;goto, skip_new
251
252
253if (flag_cb eq 'true') then begin
254
255format_colorbar=format
256ndiv=10
257;ndiv=8
258pos=[0.15, 0.17, 0.95, 0.20]
259if (subtitle_user eq '') then pos=[0.15, 0.20, 0.95, 0.23]
260
261;colorbar = Obj_New("COLORBAR",$
262;        Range=[minfield,maxfield],$
263;        Ncolors=255,$
264;        charsize=0.85,$
265;        format=format_colorbar,$
266;        major=ndiv,$
267;        ticklen=-0.15,$
268;       position=pos)
269;colorbar->Draw
270;Obj_Destroy, colorbar
271
272        colorbar, $
273;                /vertical, $
274                position=pos, $
275                range=[float(minfield),float(maxfield)], $
276                divisions=ndiv, $
277                ncolors=255,$
278;               charsize=0.85,$
279                ticklen=-0.15,$
280                format=format_colorbar
281
282endif
283
284if (n_elements(nam2) ne 0) then thedate=STRSPLIT(nam2, /EXTRACT) & thedate=STRSPLIT(thedate(2),"'",/EXTRACT) & thedate=STRSPLIT(thedate(0),"_",/EXTRACT) & thedate=STRSPLIT(thedate(0),"-",/EXTRACT) & xyouts, 0.95, pos(1)-0.08, 'LMD Martian Mesoscale Model - d'+thedate(2)+'/m'+thedate(1), charsize=0.7, alignment=1., /NORMAL
285
286skip_new:
287
288;--------------------
289; overplot contour
290;--------------------
291if (n_elements(overcontour) ne 1) then begin
292
293        w=where(abs(overcontour) lt missing_value)
294        if (w(0) ne -1) then begin
295                min_contour=min(overcontour[w])
296                max_contour=max(overcontour[w])
297        endif
298
299
300; to avoid spurious blank zones
301nlevels=colors/2
302if (n_elements(lev) eq 0) then lev=min_contour+(findgen(nlevels)/float(nlevels-1))*(max_contour-min_contour)
303nlevels=n_elements(lev)
304
305;lev=round(lev)
306;yeah=UNIQ(lev) & lev=lev[yeah]
307
308
309loadct,0,/silent
310contour, $
311        overcontour, $
312        space,altitude, $
313;       nlevels=colors/2, $
314        levels=lev, $
315        C_LINESTYLE = (lev LT 0.0), $   ;; dotted lines if negative
316        max_value=max_contour, $
317        min_value=min_contour, $       
318        c_labels=findgen(nlevels), $
319        color=0, $ ;255, $ ;0, $
320        /noerase, $
321        /overplot
322;        xtickinterval=intervalx, $
323;        ytickinterval=intervaly, $
324;;/isotropic, $
325;;xtickname=xticks, $
326;       xtitle=xtitle, $
327;       xrange=[minspace,maxspace], $
328;       ytitle=ytitle, $
329;       yrange=[minalt,maxalt], $
330;position=[0.15, 0.15, 0.95, 0.75]
331
332
333
334endif
335
336;--------------------
337; overplot wind
338;--------------------
339if (n_elements(overvector_x)*n_elements(overvector_y) ne 1) then begin
340
341        w=where(abs(overvector_x) ge missing_value)  ; u and v have the same missing values ...
342        if (w(0) ne -1) then begin
343                overvector_x[w]=0.
344                overvector_y[w]=0.
345        endif
346               
347
348len=2.0
349if (n_elements(windex) eq 0) then begin
350        ref_mag=round(max(sqrt(overvector_x^2+overvector_y^2))/2)
351        ;ref_mag=round(median(sqrt(overvector_x^2+overvector_y^2)))
352ref_mag=15.     ;;ce qui est ci dessus deconne !!!
353endif else begin
354        ref_mag=round(windex)
355endelse
356mytext=STRTRIM(STRING(ref_mag),2)+'ms!E-1!N'
357
358mytext=''
359
360if (n_elements(stride) eq 0) then begin
361stride=3                ; STRIDE: Pick every nth vector for plotting. (1)
362stride=2
363stride=5
364endif else begin
365stride=round(stride)
366endelse
367
368
369ref_pos=[0.10, 0.90]   
370ref_pos=[0.20, 0.05]
371
372myu=overvector_x
373myv=overvector_y
374
375longs=space
376lats=altitude
377
378wlon=where( (space gt minspace) and (space lt maxspace) )
379wlat=where( (altitude gt minalt) and (altitude lt maxalt) )
380
381if ((wlon(0) eq -1) or (wlat(0) eq -1)) then begin
382        print, 'wrong window settings'
383        print, space
384        print, 'lim ',minspace,maxspace
385        print, altitude
386        print, 'lim ',minalt,maxalt
387        stop
388endif
389
390myu=dblarr(n_elements(space),n_elements(altitude))
391myv=dblarr(n_elements(space),n_elements(altitude))
392
393myu(*,*)=!Values.F_NAN
394myv(*,*)=!Values.F_NAN
395
396for i=0,n_elements(wlon)-1 do begin
397for j=0,n_elements(wlat)-1 do begin
398        myu(wlon(i),wlat(j))=overvector_x(wlon(i),wlat(j))
399        myv(wlon(i),wlat(j))=overvector_y(wlon(i),wlat(j))
400endfor
401endfor
402
403
404        ;; eliminate vector for colorbar clarity ?
405        ;nxs=1
406        ;nxe=n_elements(longs(*))-2
407        ;nys=1
408        ;nye=n_elements(lats(*))-2
409        ;myu=overvector_x(nxs:nxe,nys:nye)
410        ;myv=overvector_y(nxs:nxe,nys:nye)
411        ;longs=space(nxs:nxe)
412        ;lats=altitude(nys:nye)
413
414
415;;---------------------------
416;;
417;factor=float(mean(overvector_x^2))/float(mean(overvector_y^2))
418;factor=sqrt(factor)
419;factor=factor/2
420;
421;factor=6.
422;factor=12.
423;factor=15.
424;
425;print, 'X/Z correction factor', factor
426;
427;;rapport ordre de grandeur entre mouvements verticaux et horizontaux
428;;
429;;---------------------------
430
431
432myv=myv*factor
433
434
435loadct,0,/silent
436VECTOR, myu, myv, longs, lats, LENGTH=len, REF_MAG=ref_mag,$
437        STRIDE=stride, TYPE=2, HEAD_LEN=0.3, ANGLE=30, $
438;       REF_POS=ref_pos, $
439;       REF_TEXT=mytext,$
440        ALIGN=0.5, COLOR=1, /OVERPLOT;, /NOERASE
441
442endif
443
444
445endelse
446
447space=space_save        ;; careful with sequential call to section.pro
448                        ;; combined with modified space (e.g. with factor)
449                        ;; -- this save fix the issue
450
451
452end
Note: See TracBrowser for help on using the repository browser.