source: trunk/mesoscale/LMD_MM_MARS/SRC/ARWpost/idl/map_latlon.pro @ 43

Last change on this file since 43 was 11, checked in by aslmd, 15 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien

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