source: trunk/MESOSCALE_DEV/PLOT/MINIMAL/map_latlon.pro @ 937

Last change on this file since 937 was 167, checked in by aslmd, 14 years ago

MESOSCALE: correction makemeso (LES). ajout a README-svn. tests tempetes de poussiere avec JF.

File size: 17.3 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
222if (pal eq 21) then reverse_ct
223contour, what_I_plot, $
224        lon,lat, $
225;       nlevels=colors, $
226        levels=levp, $
227        /cell_fill, $
228        max_value=maxfield, $
229        min_value=minfield, $
230        /overplot
231
232
233;nlat=128.
234;oplot, lat*0.+nlat*100./1000., linestyle=1
235
236;;; ANOMALY
237;lon=lons
238;lat=lats
239
240;------------------
241; colorbar
242;------------------
243
244;if (flag_cb eq 'true') then begin
245;        format_colorbar=format
246;        colorbar, $
247;                position=[0.15, 0.85, 0.95, 0.90], $
248;                divisions=8, $
249;                range=[float(minfield),float(maxfield)], $
250;                format=format_colorbar
251;endif
252;goto, skip_new
253
254if (flag_cb eq 'true') then begin
255format_colorbar=format
256
257;;; ndiv peut etre regle
258if (n_elements(ndiv) eq 0) then ndiv=10
259
260;if (n_elements(poscb) eq 0) then poscb=0.95
261;if (n_elements(poscb) eq 0) then poscb=0.80
262if (n_elements(poscb) eq 0) then poscb=0.70
263
264     if (isotropic eq 'true') then begin
265       ;poscb=0.85 &
266       pos=[poscb, 0.12, poscb+0.03, 0.92]
267       ;poscb=0.95 & pos=[poscb, 0.12, poscb+0.03, 0.92]
268;       poscb=0.96 & pos=[poscb, 0.12, poscb+0.03, 0.92]
269;       poscb=0.05 & pos=[poscb, 0.12, poscb+0.03, 0.92]
270               ;colorbar = Obj_New("COLORBAR",$
271               ; Range=[float(minfield),float(maxfield)],$
272               ; /vertical, $
273               ; Ncolors=255,$
274               ; charsize=0.85,$
275               ; format=format_colorbar,$
276               ; major=ndiv,$
277               ; ticklen=-0.15,$
278               ; position=pos)
279               ;colorbar->Draw
280               ;Obj_Destroy, colorbar
281        colorbar, $
282                /vertical, $
283                position=pos, $
284                range=[float(minfield),float(maxfield)], $
285                divisions=ndiv, $
286                ncolors=255,$
287;               charsize=0.85,$
288                ticklen=-0.15,$
289                format=format_colorbar
290
291xyouts, pos(2), 0.01, subtitle_user, charsize=0.8, alignment=1., /NORMAL
292
293        if (n_elements(nam2) ne 0) then begin
294                thedate=STRSPLIT(nam2, /EXTRACT)
295                thedate=STRSPLIT(thedate(2),"'",/EXTRACT)
296                thedate=STRSPLIT(thedate(0),"_",/EXTRACT)
297                thedate=STRSPLIT(thedate(0),"-",/EXTRACT)
298                xyouts, pos(2), 0.01, 'LMD Martian Mesoscale Model - d'+thedate(2)+'/m'+thedate(1), charsize=0.7, alignment=1., /NORMAL
299        endif else begin
300                thedate=''
301        endelse
302;;
303;;
304;;
305;xyouts, pos(2), 0.01, 'LMD Mars Mesoscale Model / '+subtitle_hunhun, charsize=0.95, alignment=1., /NORMAL
306;;
307;;
308;;
309
310     endif else begin
311
312        pos=[0.15, 0.17, 0.95, 0.20]
313        if (subtitle_user eq '') then pos=[0.15, 0.20, 0.95, 0.23]
314
315               ;colorbar = Obj_New("COLORBAR",$
316               ; Range=[float(minfield),float(maxfield)],$
317               ; Ncolors=255,$
318               ; charsize=0.85,$
319               ; format=format_colorbar,$
320               ; major=ndiv,$
321               ; ticklen=-0.15,$
322               ; position=pos)
323               ;colorbar->Draw
324               ;Obj_Destroy, colorbar
325        colorbar, $
326                position=pos, $
327                range=[float(minfield),float(maxfield)], $
328                divisions=ndiv, $
329                ncolors=255,$
330;               charsize=0.85,$
331                ticklen=-0.15,$
332                format=format_colorbar
333
334        if (n_elements(nam2) ne 0) then begin
335                thedate=STRSPLIT(nam2, /EXTRACT)
336                thedate=STRSPLIT(thedate(2),"'",/EXTRACT)
337                thedate=STRSPLIT(thedate(0),"_",/EXTRACT)
338                thedate=STRSPLIT(thedate(0),"-",/EXTRACT)
339                xyouts, 0.95, pos(1)-0.08, 'LMD Martian Mesoscale Model - d'+thedate(2)+'/m'+thedate(1), charsize=0.7, alignment=1., /NORMAL
340        endif else begin
341                thedate=''
342        endelse
343
344      endelse   
345endif
346
347skip_new:
348
349;;;;;;;;;;;;;;;;;;;;;;;;
350;;;;;;;;;;;;;;;;;;;;;;;;
351;;; ajouter un point
352;loadct, 0
353;xyouts, 353.87-360, -1.88, '+', ALIGNMENT=0.5, CHARSIZE=2, CHARTHICK=2
354;;;;;;;;;;;;;;;;;;;;;;;;
355;;;;;;;;;;;;;;;;;;;;;;;;
356
357;goto, nocontour
358;print, 'contour'
359;--------------------
360; overplot contour
361;--------------------
362if (n_elements(overcontour) ne 1) then begin
363
364        w=where(abs(overcontour) lt missing_value)
365        if (w(0) ne -1) then begin
366                min_contour=min(overcontour[w])
367                max_contour=max(overcontour[w])
368        endif
369
370
371; to avoid spurious blank zones
372nlevels=colors/2
373if (n_elements(lev) eq 0) then lev=min_contour+(findgen(nlevels)/float(nlevels-1))*(max_contour-min_contour)
374
375;lev=round(lev)
376;yeah=UNIQ(lev) & lev=lev[yeah]
377
378colcont = 150 ;75
379if (pal eq 16) then colcont=0 ;noir
380if (pal eq 4) then colcont=255 ;blanc
381
382if (isotropic eq 'true') then begin
383loadct,0,/silent
384;reverse_ct
385contour, overcontour, $
386        lon,lat, $
387/overplot, $
388;       nlevels=colors/2, $
389        levels=lev, $
390        C_LINESTYLE = (lev LT 0.0), $   ;; dotted lines if negative
391        max_value=max_contour, $
392        min_value=min_contour, $       
393        c_labels=findgen(n_elements(lev))*0., $
394;        c_charsize=!P.charsize/2, $
395        color=colcont;75;0;50;75;100tropclair;150;0;, $
396;       /noerase, $
397;       xtitle=xtitleset, $
398;       xrange=windowx, $
399;       ytitle=ytitleset, $
400;       yrange=windowy, $
401;       position=pospos, $
402;       /isotropic, $
403;       xtickinterval=intervalx, $
404;       ytickinterval=intervaly
405endif else begin
406loadct,0,/silent
407contour, overcontour, $
408        lon,lat, $
409/overplot, $
410;        nlevels=colors/2, $
411        levels=lev, $
412        C_LINESTYLE = (lev LT 0.0), $   ;; dotted lines if negative
413        max_value=max_contour, $
414        min_value=min_contour, $
415        c_labels=findgen(n_elements(lev))*0., $
416        color=colcont;, $
417;        /noerase, $
418;        xtitle=xtitleset, $
419;        xrange=windowx, $
420;        ytitle=ytitleset, $
421;        yrange=windowy, $
422;;       position=[0.15, 0.15, 0.95, 0.75], $
423;        xtickinterval=intervalx, $
424;        ytickinterval=intervaly
425endelse
426
427
428endif
429nocontour:
430
431;--------------------
432; overplot wind
433;--------------------
434if (n_elements(overvector_x)*n_elements(overvector_y) ne 1) then begin
435
436        w=where(abs(overvector_x) ge missing_value)  ; u and v have the same missing values ...
437        if (w(0) ne -1) then begin
438                overvector_x[w]=0.
439                overvector_y[w]=0.
440        endif
441               
442
443len=2.0
444if (n_elements(windex) eq 0) then windex=20.
445ref_mag=round(windex)
446mytext=STRTRIM(STRING(ref_mag),2)+'ms!E-1!N'
447
448if (n_elements(stride) eq 0) then begin
449;stride=3               ; STRIDE: Pick every nth vector for plotting. (1)
450;stride=2
451stride=5
452endif else begin
453stride=round(stride)
454endelse
455       
456ref_pos=[0.10, 0.90]   
457ref_pos=[0.15, 0.10]
458ref_pos=[0.20, 0.05]
459;ref_pos=[0.20, 0.02]
460
461s=size(lon)
462if (s[0] eq 1) then begin
463 zeu=overvector_x
464 zev=overvector_y
465 longs=lon
466 lats=lat
467endif else begin
468 minlat=min(lat) & maxlat=max(lat) & minlon=min(lon) & maxlon=max(lon)
469 npoints=2*n_elements(lon(*,0))  ;; trop de points, mais au moins on ne perd rien (2 et 3 donnent results similaires)
470 TRIANGULATE, lon, lat, tr
471 zeu = GRIDDATA( lon, lat, overvector_x, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN )
472 zev = GRIDDATA( lon, lat, overvector_y, /LINEAR, triangles=tr, dimension=npoints, MISSING=!VALUES.F_NAN )
473 longs =  minlon + (maxlon - minlon)*findgen(npoints)/float(npoints-1)
474 lats =  minlat + (maxlat - minlat)*findgen(npoints)/float(npoints-1)
475endelse
476
477wlon=where( (longs gt windowx(0)) and (longs lt windowx(1)) )
478wlat=where( (lats gt windowy(0)) and (lats lt windowy(1)) )
479
480if ((wlon(0) eq -1) or (wlat(0) eq -1)) then begin
481        print, 'wrong window settings'
482        print, windowx(0), windowx(1)
483        print, longs
484        print, windowy(0), windowy(1)
485        print, lats
486        stop
487endif
488
489myu=dblarr(n_elements(longs),n_elements(lats))
490myv=dblarr(n_elements(longs),n_elements(lats))
491
492myu(*,*)=!Values.F_NAN
493myv(*,*)=!Values.F_NAN
494
495for i=0,n_elements(wlon)-1 do begin
496for j=0,n_elements(wlat)-1 do begin     
497        myu(wlon(i),wlat(j))=zeu(wlon(i),wlat(j))
498        myv(wlon(i),wlat(j))=zev(wlon(i),wlat(j))
499endfor 
500endfor 
501
502        ;; eliminate vector for colorbar clarity ?
503        ;nxs=1
504        ;nxe=n_elements(longs(*))-2
505        ;nys=1
506        ;nye=n_elements(lats(*))-2
507        ;myu=overvector_x(nxs:nxe,nys:nye)
508        ;myv=overvector_y(nxs:nxe,nys:nye)
509        ;longs=lon(nxs:nxe)
510        ;lats=lat(nys:nye)
511
512loadct,0,/silent
513VECTOR, myu, myv, longs, lats, LENGTH=len, REF_MAG=ref_mag,$
514        STRIDE=stride, TYPE=2, HEAD_LEN=0.3, ANGLE=30, $
515        REF_POS=ref_pos, REF_TEXT=mytext,$
516        ALIGN=0.5, COLOR=1, /OVERPLOT;, /NOERASE
517
518yeah:
519endif
520
521
522endelse
523
524
525
526goto, why_contour_end
527;;;;;;;;;;;;;;;;;;;;;;
528;;;;;;;;;;;;;;;;;;;;;;
529;;;;;;;;;;;;;;;;;;;;;;
530
531if (n_elements(overcontour) ne 1) then begin
532
533        w=where(abs(overcontour) lt missing_value)
534        if (w(0) ne -1) then begin
535                min_contour=min(overcontour[w])
536                max_contour=max(overcontour[w])
537        endif
538
539
540; to avoid spurious blank zones
541nlevels=colors/2
542if (n_elements(lev) eq 0) then lev=min_contour+(findgen(nlevels)/float(nlevels-1))*(max_contour-min_contour)
543;lev=round(lev)
544;yeah=UNIQ(lev) & lev=lev[yeah]
545
546
547if (isotropic eq 'true') then begin
548loadct,0,/silent
549contour, overcontour, $
550        lon,lat, $
551;       nlevels=colors/2, $
552        levels=lev, $
553        C_LINESTYLE = (lev LT 0.0), $   ;; dotted lines if negative
554        max_value=max_contour, $
555        min_value=min_contour, $
556        c_labels=findgen(nlevels), $ ;*0., $
557        color=0, $
558        /noerase, $
559        xtitle=xtitleset, $
560        xrange=windowx, $
561        ytitle=ytitleset, $
562        yrange=windowy, $
563;        position=[0.15, 0.15, 0.95, 0.75], $
564        /isotropic, $
565        xtickinterval=intervalx, $
566        ytickinterval=intervaly
567endif else begin
568loadct,0,/silent
569contour, overcontour, $
570        lon,lat, $
571;       nlevels=colors/2, $
572        levels=lev, $
573        C_LINESTYLE = (lev LT 0.0), $   ;; dotted lines if negative
574        max_value=max_contour, $
575        min_value=min_contour, $
576        c_labels=findgen(nlevels), $
577        color=0, $
578        /noerase, $
579        xtitle=xtitleset, $
580        xrange=windowx, $
581        ytitle=ytitleset, $
582        yrange=windowy, $
583        position=[0.15, 0.15, 0.95, 0.75], $
584        xtickinterval=intervalx, $
585        ytickinterval=intervaly
586endelse
587
588
589endif
590
591why_contour_end:
592
593end
Note: See TracBrowser for help on using the repository browser.