source: trunk/MESOSCALE/PLOT/MINIMAL/map_latlon.pro @ 134

Last change on this file since 134 was 131, checked in by aslmd, 14 years ago

mars

M 130 mars/libf/phymars/callradite.F
M 130 mars/README
no more need to modify callradite.F prior to compilation
[but still dimradmars.h must be modified]

--> in callradite.F we now have

-- DEFAULT name_iaer(1) is "dust_conrath"
-- IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq"
-- IF (water.AND.activice) name_iaer(2) = "h2o_ice"

M 130 000-USERS
small meeting about SVN for "planeto" team

M 130 mars/libf/phymars/meso_physiq.F
purely cosmetic changes

M 130 mesoscale/PLOT/MINIMAL/map_latlon.pro
M 130 mesoscale/PLOT/SPEC/GW/gravitwave2.pro
M 130 mesoscale/PLOT/SPEC/GW/gravitwaveprof.pro
M 130 mesoscale/PLOT/SPEC/GW/gravitwave_inc.pro
A 0 mesoscale/PLOT/SPEC/MAP/defs/THARSIS_newphys_bugnoRAC_RACgcm_TAUTES.map_uvt_inc.pro
M 130 mesoscale/PLOT/SPEC/MAP/map_uvt.pro
graphic stuff for GW studies and H2O cycle
File size: 17.2 KB
Line 
1pro map_latlon, $
2        what_I_plot, $                          ; 2D field
3        lon, $                                  ; 1D latitude OR 2D
4        lat, $                                  ; 1D longitude OR 2D
5        minfield=minfield_init, $               ; minimum value of plotted field (=0: calculate)
6        maxfield=maxfield_init, $               ; maximum value of plotted field (=0: calculate)
7        overcontour=overcontour, $              ; another 2D field to overplot with contour lines (=0: no)
8        overvector_x=overvector_x, $            ; wind vector - x component (=0: no)
9        overvector_y=overvector_y, $            ; wind vector - y component (=0: no)
10        ct=pal, $                               ; color table (33-rainbow is default)
11        colors=colors, $                        ; number of colors/levels (32 is default)
12        title=title_user, $                     ; title of the plot ('' is default)
13        format=format                           ; format of colorbar annotations ('(F6.2)' is default)
14
15;--------------------------------------------
16;
17; A general routine to plot lat/lon maps
18;
19; USE:  - map_latlon, what_I_plot       
20;       - map_latlon, what_I_plot, lon, lat
21;       - ... and see options above
22;
23;--------------------------------------------
24; A. Spiga, August 2007
25;--------------------------------------------
26
27
28
29
30;---------------------------
31; a few permanent settings
32;---------------------------
33;
34missing_value=1.e30
35flag_cb='true'
36;-- optional:
37;-- user defined parameters
38; do not forget to type '.compile map_latlon' to update changes
39SPAWN, 'touch param_plot.idl'
40@param_plot.idl
41;
42;------------------
43
44
45;------------------
46; user settings
47;------------------
48if (n_elements(lon) eq 0) then begin
49        lon=findgen(n_elements(what_I_plot(*,0)))
50endif
51if (n_elements(lat) eq 0) then lat=findgen(n_elements(what_I_plot(0,*)))
52if (n_elements(minfield_init) eq 0) then minfield_init=0.
53if (n_elements(maxfield_init) eq 0) then maxfield_init=0.
54if (n_elements(overcontour) eq 0) then overcontour=0.
55if ((n_elements(overvector_x) eq 0) or (n_elements(overvector_y) eq 0)) then begin
56        overvector_x=0.
57        overvector_y=0.
58endif
59if (n_elements(colors) eq 0) then colors=32
60if (n_elements(title_user) eq 0) then title_user=''
61if (n_elements(subtitle_user) eq 0) then subtitle_user=''
62if (n_elements(pal) eq 0) then pal=33
63if (n_elements(format) eq 0) then format='(F6.2)'
64if (format eq '') then format='(F6.2)'
65;if (n_elements(xtitleset) eq 0) then xtitleset='longitude'
66;if (n_elements(ytitleset) eq 0) then ytitleset='latitude'
67if (n_elements(title_axis) ne 2) then begin
68        xtitleset='longitude'
69        ytitleset='latitude'
70endif else begin
71        xtitleset=title_axis(0)
72        ytitleset=title_axis(1)
73endelse
74if (n_elements(isotropic) eq 0) then isotropic='true'
75if (n_elements(coord2d) eq 0) then coord2d='false'
76
77;------------------
78; limits
79;------------------
80;if (minfield_init*maxfield_init eq 0) then begin   
81if (minfield_init eq maxfield_init) then begin
82        ;---different min/max for each layer
83        w=where(abs(what_I_plot) lt missing_value)
84        if (w(0) ne -1) then begin
85                minfield=min(what_I_plot[w])
86                maxfield=max(what_I_plot[w])
87;goodscale, $
88;     what_I_plot, $
89;     minfield, $
90;     maxfield
91        endif else begin
92                what_I_plot=0.
93                print, 'no valid values in this field'
94        endelse
95endif else begin
96        ;---user-defined limits
97        minfield=minfield_init
98        maxfield=maxfield_init
99        ;;w=where(what_I_plot lt minfield) & if (w(0) ne -1) then what_I_plot[w]=minfield
100        ;;w=where(what_I_plot gt maxfield) & if (w(0) ne -1) then what_I_plot[w]=maxfield
101endelse
102
103if ((minfield+maxfield eq 0) and (abs(minfield) ne abs(maxfield))) then begin
104        print, 'nothing to plot ... skipping this plot ...'
105endif else begin
106
107
108;-----------------------------------------
109; fix the possible longitude ugliness :)
110; ...get rid of the -180/180 limit
111;-----------------------------------------
112
113;;;;; PB WITH 2D COORDINATES
114;w=where(lon eq min(lon))
115;if (w(0) ne 0) then begin
116;       lon[0:w(0)-1]=lon[0:w(0)-1]-360
117;endif
118
119
120;---------------
121; adapt ticks
122;---------------
123if (n_elements(windowx) eq 0) then windowx=[min(lon),max(lon)]
124if (n_elements(windowy) eq 0) then windowy=[min(lat),max(lat)]
125if (n_elements(intervalx) eq 0) then intervalx=round((windowx(1)-windowx(0))/8.)
126if (n_elements(intervaly) eq 0) then intervaly=round((windowy(1)-windowy(0))/5.)
127
128
129;------------------
130; plot window
131;------------------
132
133        ;lons=lon & lats=lat
134        ;;--------------------------------------------
135        ;;
136        ;; ANOMALY
137        ;;
138        ;w=where( (lon ge windowx(0)) and (lon le windowx(1)) ) & w2=where( (lat ge windowy(0)) and (lat le windowy(1)) )
139        ;yeah=fltarr(n_elements(w),n_elements(w2)) & for i=0,n_elements(w)-1 do for j=0,n_elements(w2)-1 do yeah(i,j)=what_I_plot(w(i),w2(j))
140        ;what_I_plot = yeah - mean(yeah)
141        ;lon=lon[w] & lat=lat[w2]
142        ;;--------------------------------------------
143
144
145        ; to ensure the right limits
146        dumb_what_I_plot=what_I_plot
147        dumb_what_I_plot(*,*)=minfield & print, minfield
148        dumb_what_I_plot(0,0)=maxfield & print, maxfield
149
150;;
151;;
152;subtitle_hunhun = subtitle_user
153;subtitle_user = ''
154;;
155;;
156
157
158pospos = [0.10, 0.12, 0.90, 0.92] ;; iso
159;pospos = [0.15, 0.35, 0.95, 0.95] ;; non-iso
160
161if (isotropic eq 'true') then begin
162        loadct, 0, /silent
163        if ( coord2d ne 'polar' ) then begin
164                ;;trace non polaire
165                contour, /nodata, $
166                         /isotropic, $
167                         position=pospos, $
168                         dumb_what_I_plot,      lon,                            lat, $
169                         title=title_user,      xtitle=xtitleset,               ytitle=ytitleset, $
170                         xrange=windowx,                yrange=windowy, $
171                         xtickinterval=intervalx,       ytickinterval=intervaly, $
172                         subtitle=subtitle_user, color=0
173        endif else begin
174                ;;trace avec MAP_SET -- MAP_SET -- MAP_SET
175                print,'------------------- Polar plot'
176                contour, /nodata, $
177                         /isotropic, $
178                         /overplot, $
179                         dumb_what_I_plot,      lon,                            lat, $
180                         title=title_user,      xtitle=xtitleset,               ytitle=ytitleset, $
181                         xrange=windowx,                yrange=windowy, $
182                         xtickinterval=intervalx,       ytickinterval=intervaly, $
183                         subtitle=subtitle_user, color=0
184        endelse
185        subtitle_user = ''
186endif else begin
187        pospos = [0.15, 0.35, 0.95, 0.95] ;; non-iso
188        loadct, 0, /silent
189        if ( coord2d ne 'polar' ) then begin
190                ;;trace non polaire
191                contour, /nodata, $
192                         position=pospos, $
193                         dumb_what_I_plot,      lon,                            lat, $
194                         title=title_user,      xtitle=xtitleset,               ytitle=ytitleset, $
195                         xrange=windowx,                yrange=windowy, $
196                         xtickinterval=intervalx,       ytickinterval=intervaly, $
197                         color=0,subtitle=subtitle_user
198        endif else begin
199                ;;trace avec MAP_SET -- MAP_SET -- MAP_SET
200                print,'------------------- Polar plot'
201                contour, /nodata, $
202                         /overplot, $
203                         dumb_what_I_plot,      lon,                            lat, $
204                         title=title_user,      xtitle=xtitleset,               ytitle=ytitleset, $
205                         xrange=windowx,                yrange=windowy, $
206                         xtickinterval=intervalx,       ytickinterval=intervaly, $
207                         color=0,subtitle=subtitle_user
208        endelse
209endelse
210
211;------------------
212; plot field
213;------------------
214
215; to avoid spurious blank zones
216nlevels=colors
217levp=minfield+(findgen(nlevels)/float(nlevels-1))*(maxfield-minfield)
218
219loadct,pal,/silent
220;reverse_ct
221if (pal eq 0) then reverse_ct
222contour, what_I_plot, $
223        lon,lat, $
224;       nlevels=colors, $
225        levels=levp, $
226        /cell_fill, $
227        max_value=maxfield, $
228        min_value=minfield, $   
229        /overplot
230
231
232;nlat=128.
233;oplot, lat*0.+nlat*100./1000., linestyle=1
234
235;;; ANOMALY
236;lon=lons
237;lat=lats
238
239;------------------
240; colorbar
241;------------------
242
243;if (flag_cb eq 'true') then begin
244;        format_colorbar=format
245;        colorbar, $
246;                position=[0.15, 0.85, 0.95, 0.90], $
247;                divisions=8, $
248;                range=[float(minfield),float(maxfield)], $
249;                format=format_colorbar
250;endif
251;goto, skip_new
252
253if (flag_cb eq 'true') then begin
254format_colorbar=format
255
256;;; ndiv peut etre regle
257if (n_elements(ndiv) eq 0) then ndiv=10
258
259;if (n_elements(poscb) eq 0) then poscb=0.95
260;if (n_elements(poscb) eq 0) then poscb=0.80
261if (n_elements(poscb) eq 0) then poscb=0.70
262
263     if (isotropic eq 'true') then begin
264       ;poscb=0.85 &
265       pos=[poscb, 0.12, poscb+0.03, 0.92]
266       ;poscb=0.95 & pos=[poscb, 0.12, poscb+0.03, 0.92]
267;       poscb=0.96 & pos=[poscb, 0.12, poscb+0.03, 0.92]
268;       poscb=0.05 & pos=[poscb, 0.12, poscb+0.03, 0.92]
269               ;colorbar = Obj_New("COLORBAR",$
270               ; Range=[float(minfield),float(maxfield)],$
271               ; /vertical, $
272               ; Ncolors=255,$
273               ; charsize=0.85,$
274               ; format=format_colorbar,$
275               ; major=ndiv,$
276               ; ticklen=-0.15,$
277               ; position=pos)
278               ;colorbar->Draw
279               ;Obj_Destroy, colorbar
280        colorbar, $
281                /vertical, $
282                position=pos, $
283                range=[float(minfield),float(maxfield)], $
284                divisions=ndiv, $
285                ncolors=255,$
286;               charsize=0.85,$
287                ticklen=-0.15,$
288                format=format_colorbar
289
290xyouts, pos(2), 0.01, subtitle_user, charsize=0.8, alignment=1., /NORMAL
291
292        if (n_elements(nam2) ne 0) then begin
293                thedate=STRSPLIT(nam2, /EXTRACT)
294                thedate=STRSPLIT(thedate(2),"'",/EXTRACT)
295                thedate=STRSPLIT(thedate(0),"_",/EXTRACT)
296                thedate=STRSPLIT(thedate(0),"-",/EXTRACT)
297                xyouts, pos(2), 0.01, 'LMD Martian Mesoscale Model - d'+thedate(2)+'/m'+thedate(1), charsize=0.7, alignment=1., /NORMAL
298        endif else begin
299                thedate=''
300        endelse
301;;
302;;
303;;
304;xyouts, pos(2), 0.01, 'LMD Mars Mesoscale Model / '+subtitle_hunhun, charsize=0.95, alignment=1., /NORMAL
305;;
306;;
307;;
308
309     endif else begin
310
311        pos=[0.15, 0.17, 0.95, 0.20]
312        if (subtitle_user eq '') then pos=[0.15, 0.20, 0.95, 0.23]
313
314               ;colorbar = Obj_New("COLORBAR",$
315               ; Range=[float(minfield),float(maxfield)],$
316               ; Ncolors=255,$
317               ; charsize=0.85,$
318               ; format=format_colorbar,$
319               ; major=ndiv,$
320               ; ticklen=-0.15,$
321               ; position=pos)
322               ;colorbar->Draw
323               ;Obj_Destroy, colorbar
324        colorbar, $
325                position=pos, $
326                range=[float(minfield),float(maxfield)], $
327                divisions=ndiv, $
328                ncolors=255,$
329;               charsize=0.85,$
330                ticklen=-0.15,$
331                format=format_colorbar
332
333        if (n_elements(nam2) ne 0) then begin
334                thedate=STRSPLIT(nam2, /EXTRACT)
335                thedate=STRSPLIT(thedate(2),"'",/EXTRACT)
336                thedate=STRSPLIT(thedate(0),"_",/EXTRACT)
337                thedate=STRSPLIT(thedate(0),"-",/EXTRACT)
338                xyouts, 0.95, pos(1)-0.08, 'LMD Martian Mesoscale Model - d'+thedate(2)+'/m'+thedate(1), charsize=0.7, alignment=1., /NORMAL
339        endif else begin
340                thedate=''
341        endelse
342
343      endelse   
344endif
345
346skip_new:
347
348;;;;;;;;;;;;;;;;;;;;;;;;
349;;;;;;;;;;;;;;;;;;;;;;;;
350;;; ajouter un point
351;loadct, 0
352;xyouts, 353.87-360, -1.88, '+', ALIGNMENT=0.5, CHARSIZE=2, CHARTHICK=2
353;;;;;;;;;;;;;;;;;;;;;;;;
354;;;;;;;;;;;;;;;;;;;;;;;;
355
356;goto, nocontour
357;print, 'contour'
358;--------------------
359; overplot contour
360;--------------------
361if (n_elements(overcontour) ne 1) then begin
362
363        w=where(abs(overcontour) lt missing_value)
364        if (w(0) ne -1) then begin
365                min_contour=min(overcontour[w])
366                max_contour=max(overcontour[w])
367        endif
368
369
370; to avoid spurious blank zones
371nlevels=colors/2
372if (n_elements(lev) eq 0) then lev=min_contour+(findgen(nlevels)/float(nlevels-1))*(max_contour-min_contour)
373
374;lev=round(lev)
375;yeah=UNIQ(lev) & lev=lev[yeah]
376
377colcont = 75
378if (pal eq 16) then colcont=0 ;noir
379if (pal eq 4) then colcont=255 ;blanc
380
381if (isotropic eq 'true') then begin
382loadct,0,/silent
383;reverse_ct
384contour, overcontour, $
385        lon,lat, $
386/overplot, $
387;       nlevels=colors/2, $
388        levels=lev, $
389        C_LINESTYLE = (lev LT 0.0), $   ;; dotted lines if negative
390        max_value=max_contour, $
391        min_value=min_contour, $       
392        c_labels=findgen(n_elements(lev))*1., $
393;        c_charsize=!P.charsize/2, $
394        color=colcont;75;0;50;75;100tropclair;150;0;, $
395;       /noerase, $
396;       xtitle=xtitleset, $
397;       xrange=windowx, $
398;       ytitle=ytitleset, $
399;       yrange=windowy, $
400;       position=pospos, $
401;       /isotropic, $
402;       xtickinterval=intervalx, $
403;       ytickinterval=intervaly
404endif else begin
405loadct,0,/silent
406contour, overcontour, $
407        lon,lat, $
408/overplot, $
409;        nlevels=colors/2, $
410        levels=lev, $
411        C_LINESTYLE = (lev LT 0.0), $   ;; dotted lines if negative
412        max_value=max_contour, $
413        min_value=min_contour, $
414        c_labels=findgen(n_elements(lev))*0., $
415        color=colcont;, $
416;        /noerase, $
417;        xtitle=xtitleset, $
418;        xrange=windowx, $
419;        ytitle=ytitleset, $
420;        yrange=windowy, $
421;;       position=[0.15, 0.15, 0.95, 0.75], $
422;        xtickinterval=intervalx, $
423;        ytickinterval=intervaly
424endelse
425
426
427endif
428nocontour:
429
430;--------------------
431; overplot wind
432;--------------------
433if (n_elements(overvector_x)*n_elements(overvector_y) ne 1) then begin
434
435        w=where(abs(overvector_x) ge missing_value)  ; u and v have the same missing values ...
436        if (w(0) ne -1) then begin
437                overvector_x[w]=0.
438                overvector_y[w]=0.
439        endif
440               
441
442len=2.0
443if (n_elements(windex) eq 0) then windex=20.
444ref_mag=round(windex)
445mytext=STRTRIM(STRING(ref_mag),2)+'ms!E-1!N'
446
447if (n_elements(stride) eq 0) then begin
448;stride=3               ; STRIDE: Pick every nth vector for plotting. (1)
449;stride=2
450stride=5
451endif else begin
452stride=round(stride)
453endelse
454       
455ref_pos=[0.10, 0.90]   
456ref_pos=[0.15, 0.10]
457ref_pos=[0.20, 0.05]
458;ref_pos=[0.20, 0.02]
459
460s=size(lon)
461if (s[0] eq 1) then begin
462 zeu=overvector_x
463 zev=overvector_y
464 longs=lon
465 lats=lat
466endif else begin
467 minlat=min(lat) & maxlat=max(lat) & minlon=min(lon) & maxlon=max(lon)
468 npoints=2*n_elements(lon(*,0))  ;; trop de points, mais au moins on ne perd rien (2 et 3 donnent results similaires)
469 TRIANGULATE, lon, lat, tr
470 zeu = GRIDDATA( lon, lat, overvector_x, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN )
471 zev = GRIDDATA( lon, lat, overvector_y, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN )
472 longs =  minlon + (maxlon - minlon)*findgen(npoints)/float(npoints-1)
473 lats =  minlat + (maxlat - minlat)*findgen(npoints)/float(npoints-1)
474endelse
475
476wlon=where( (longs gt windowx(0)) and (longs lt windowx(1)) )
477wlat=where( (lats gt windowy(0)) and (lats lt windowy(1)) )
478
479if ((wlon(0) eq -1) or (wlat(0) eq -1)) then begin
480        print, 'wrong window settings'
481        print, windowx(0), windowx(1)
482        print, longs
483        print, windowy(0), windowy(1)
484        print, lats
485        stop
486endif
487
488myu=dblarr(n_elements(longs),n_elements(lats))
489myv=dblarr(n_elements(longs),n_elements(lats))
490
491myu(*,*)=!Values.F_NAN
492myv(*,*)=!Values.F_NAN
493
494for i=0,n_elements(wlon)-1 do begin
495for j=0,n_elements(wlat)-1 do begin     
496        myu(wlon(i),wlat(j))=zeu(wlon(i),wlat(j))
497        myv(wlon(i),wlat(j))=zev(wlon(i),wlat(j))
498endfor 
499endfor 
500
501        ;; eliminate vector for colorbar clarity ?
502        ;nxs=1
503        ;nxe=n_elements(longs(*))-2
504        ;nys=1
505        ;nye=n_elements(lats(*))-2
506        ;myu=overvector_x(nxs:nxe,nys:nye)
507        ;myv=overvector_y(nxs:nxe,nys:nye)
508        ;longs=lon(nxs:nxe)
509        ;lats=lat(nys:nye)
510
511loadct,0,/silent
512VECTOR, myu, myv, longs, lats, LENGTH=len, REF_MAG=ref_mag,$
513        STRIDE=stride, TYPE=2, HEAD_LEN=0.3, ANGLE=30, $
514        REF_POS=ref_pos, REF_TEXT=mytext,$
515        ALIGN=0.5, COLOR=1, /OVERPLOT;, /NOERASE
516
517yeah:
518endif
519
520
521endelse
522
523
524
525goto, why_contour_end
526;;;;;;;;;;;;;;;;;;;;;;
527;;;;;;;;;;;;;;;;;;;;;;
528;;;;;;;;;;;;;;;;;;;;;;
529
530if (n_elements(overcontour) ne 1) then begin
531
532        w=where(abs(overcontour) lt missing_value)
533        if (w(0) ne -1) then begin
534                min_contour=min(overcontour[w])
535                max_contour=max(overcontour[w])
536        endif
537
538
539; to avoid spurious blank zones
540nlevels=colors/2
541if (n_elements(lev) eq 0) then lev=min_contour+(findgen(nlevels)/float(nlevels-1))*(max_contour-min_contour)
542;lev=round(lev)
543;yeah=UNIQ(lev) & lev=lev[yeah]
544
545
546if (isotropic eq 'true') then begin
547loadct,0,/silent
548contour, overcontour, $
549        lon,lat, $
550;       nlevels=colors/2, $
551        levels=lev, $
552        C_LINESTYLE = (lev LT 0.0), $   ;; dotted lines if negative
553        max_value=max_contour, $
554        min_value=min_contour, $
555        c_labels=findgen(nlevels), $ ;*0., $
556        color=0, $
557        /noerase, $
558        xtitle=xtitleset, $
559        xrange=windowx, $
560        ytitle=ytitleset, $
561        yrange=windowy, $
562;        position=[0.15, 0.15, 0.95, 0.75], $
563        /isotropic, $
564        xtickinterval=intervalx, $
565        ytickinterval=intervaly
566endif else begin
567loadct,0,/silent
568contour, overcontour, $
569        lon,lat, $
570;       nlevels=colors/2, $
571        levels=lev, $
572        C_LINESTYLE = (lev LT 0.0), $   ;; dotted lines if negative
573        max_value=max_contour, $
574        min_value=min_contour, $
575        c_labels=findgen(nlevels), $
576        color=0, $
577        /noerase, $
578        xtitle=xtitleset, $
579        xrange=windowx, $
580        ytitle=ytitleset, $
581        yrange=windowy, $
582        position=[0.15, 0.15, 0.95, 0.75], $
583        xtickinterval=intervalx, $
584        ytickinterval=intervaly
585endelse
586
587
588endif
589
590why_contour_end:
591
592end
Note: See TracBrowser for help on using the repository browser.