source: trunk/SRC/ToBeReviewed/PLOTS/DESSINE/plt.pro @ 231

Last change on this file since 231 was 231, checked in by pinsard, 17 years ago

improvements/corrections of some *.pro headers

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 27.4 KB
RevLine 
[2]1;+
2;
[231]3; @file_comments
[152]4; Draw horizontal graph (map) with CONTOUR procedure
[2]5;
[231]6; @categories
[157]7; Graphics
[2]8;
[231]9; @param TAB1 {in}{required}
[152]10; The field whose we want to make the horizontal map can be:
[67]11;         1) an array. if needed, its mean along the z and t direction
12;         will be automatically performed.
[142]13;         2) a structure respecting all criterions specified by
[231]14;         <pro>litchamp</pro> cf. IDL> xhelp,'litchamp'
15;
16; @param PARAM2 {in}{optional}{default=min of tab1 (on ocean points)}
[142]17; Min value we want to consider in the contour's drawing.
18;         Note: Could also be the type of plot that can be only 'xy' for plt
[2]19;
[231]20; @param PARAM3 {in}{optional}{default=max of tab1 (on ocean points)}
[142]21; Max value we want to consider in the contour's drawing.
22;         Note: if param2 is defined as 'xy' then param3 is used to define the min (see param2)
[2]23;
[231]24; @param PARAM4 {in}{optional}{default='xy'}
[142]25; Type of plot (can be only 'xy' for plt).
26;         Note: if param2 is defined as 'xy' then param4 is used to define the max (see param3)
[2]27;
[163]28; @keyword BOXZOOM {type=vector}
[231]29; Vector indicating the geographic zone on which we want to cut the map.
[142]30;  If BOXZOOM has :
31; 1 element : The extraction is made on [lon1, lon2, lat1, lat2, 0.,boxzoom[0]]
32; 2 elements: The extraction is made on [lon1, lon2, lat1, lat2, boxzoom[0],boxzoom[1]]
33; 4 elements: The extraction is made on [Boxzoom, 0, max([gdept, gdepw])]
34; 5 elements: The extraction is made on [Boxzoom[0:3], 0, Boxzoom[4]]
35; 6 elements: The extraction is made on Boxzoom
[2]36;
[142]37; Where lon1, lon2,lat1,lat2 are global variables defined at the last domdef!
[2]38;
[231]39; @keyword REALCONT
[142]40; Allow to draw continents defined in IDL. REALCONT can have 2 form:
[163]41;           /REALCONT: we draw continents in place of the mask
[231]42;           REALCONT=2 we draw continent's contour over the mask (this allows
[142]43;                      to see if the mask correspond at real continents).
[2]44;
[231]45; @keyword CONTINTERVALLE
46; When CONTOUR is activated, it is the value between 2 isolines
47; traced by a trait. So it can be different from the one specified by INTERVALLE which,
48; in this case, does not control colored isolines in color anymore. If there is noone
49; specified min, we choose a contour min which goes well with the specified interval!
50; If this keyword is not specified, we trace 20 isolines from the min to the max.
[2]51;
[163]52; @keyword CONTLABEL {type=integer}
53; It is an integer n. When CONTOUR is activated, if n
[231]54; is different of 0, choose the label type corresponding to n cases for
55; the traced by a traitisolines. To specify the type of label of the
[142]56; colored contour, use LABEL
[2]57;
[163]58; @keyword CONTMAX {default=The max of the array passed in the keyword CONTOUR (on ocean points)}
[231]59; When CONTOUR is activated, max value we want to consider in the isoline
[163]60; traced by a trait's line
[2]61;
[163]62; @keyword CONTMIN {default=The min of the array passed in the keyword CONTOUR (on ocean points)}
[231]63; When CONTOUR is activated, min value we want to consider in the isoline
[163]64; traced by a trait's line.
[2]65;
[163]66; @keyword CONTNLEVEL {default=20}
67; When  CONTOUR is activated, it is the number of contours
[142]68; traced by a trait for drawing (active if
[163]69; CONTLABEL=0)
[67]70;
[231]71; @keyword CONTOUR
72; If we want to trace contours of a different field than the one
73; whose we have the colored drawing (by example E-P in color and QSR in contours).
[163]74; It must be a field respecting same characteristics than the argument number one of plt.
[67]75;
[163]76; @keyword GRIDTYPE 'U','T','V','W' or 'F' to specify possibly the grid on which the field is joined.
[142]77; Comment: In this case, we should use a structure like field.
[2]78;
[231]79; @keyword INTERVALLE
80; Value of an internal between to isoline. If there is none specified min,
81; we choose a min contour which goes well with the specified interval!. If this keyword is not
82; specified, we trace 20 isoline from the min to the max. Comment: When CONTOUR is activated,
83; INTERVALLE only specify the interval between 2 colored isoline. To specify the interval
[142]84; between 2 isolines traced by a trait, use CONTINTERVALLE.
[67]85;
[231]86; @keyword INV
87; Invert the color vector used to color the graph
[142]88;               (without use the black, the white and the used palette)
[231]89;
[163]90; @keyword LABEL {type=integer}
[231]91; It is an integer n. If n different of 0, it choose the label's type
92; corresponding to cases n. cf <pro>label</pro>
93; Comment: When CONTOUR is activated, it only specify the label's type for colored isolines.
[163]94; For these one traced by a trait, use CONTLABEL.
[2]95;       
[163]96; @keyword MAP {default=Cylindrical projection}
[142]97; We use it when we want to do a projection.
98; This keyword can be of two types:
[163]99;     MAP=[P0lat,P0lon,Rot]. For the description of these 3 values (see the online help of MAP_SET).
[231]100;     /MAP: In this case, map is automatically calculated have the value:
[2]101;          map = [0, (lon1+lon2)/2., 0]
[142]102; Comment: A good way to choose the type of the projection we want to do is to have a look at IDL demo:
[2]103; IDL> demo
[142]104; Then choose earth sciences and  mapping.
[231]105; Comment2: By default it is a cylindrical projection which is effectuated (with or without the keyword map).
[142]106; If we want an other projection, MAP must be activated and we have to add the keyword: /nom_projection.
107; For example, for a polar projection centered on the south pole:
[2]108; IDL> domdef,-180,180,-90,-45
109; IDL> plt, tab, /stereo,map=[-90,0,0]
110;
[163]111; @keyword LABMAP {default=labmap=1}
112; Corresponds to label keywords of map_set.
[67]113;
[231]114; @keyword MAXIN
[152]115; to specify the max value we want to plot with a keyword instead of with the
116; input parameter max. If max is defined by both, parameter and keyword, the
117; keyword is retained.
[2]118;
[231]119; @keyword MININ
[152]120; to specify the min value we want to plot with a keyword instead of with the
121; input parameter min. If min is defined by both, parameter and keyword, the
122; keyword is retained.
[2]123;
[163]124; @keyword NLEVEL {default=20}
[142]125; Number of contour to draw by default =20. active if
126; LABEL=0 or is not specified.
[2]127;
[231]128; @keyword NOTRI
129; To force not to use the triangulation. Beware, in this case, the
130; drawing only works if the grid is undeformed (It means that each point of a longitude
131; give one latitude and each point of a latitude give one longitude) except if we use
132; the keyword CELL_FILL=2.
[142]133; Comment: if the field contain points !values.f_nan, then we even do a triangulation.
[2]134;
[231]135; @keyword OVERPLOT
136; To make a plot over an other one.
137; Comment: Contrarily to the use of CONTOUR or VECTEUR, the use of this keyword
[142]138; does not the caption and/or the color bar.
[2]139;
[231]140; @keyword STRICTFILL
141; Activate this keyword to that the filling of contours be
142; precisely done between the min and the max specified letting values inferior at the
[142]143; specified min and values superior at the specified max in white.
[2]144;
[163]145; @keyword STYLE {default=style=0}
146; Contour's style to adopt to draw isolines.
[2]147;
[163]148; @keyword VECTEUR {type=struct}
[142]149; It is a structure composed by 2 elements containing the 2 fields U
150; and V of values of the zonal and meridian component of the vector of the fields
151; to draw. These fields can be an array or a structure.
152;       For example: vecteur={matriceu:lec('unsurface'),matricev:lec('vnsurface')}
[163]153;       Comment: name of elements of vectors are inconsequential.
[142]154;       vecteur={u:lec('unsurface'),v:lec('vnsurface')} is also appropriated.
[2]155;
[142]156; @keyword _EXTRA
[231]157; Used to pass keywords
158;
159; @keyword CONT_NOFILL
160; Activate it not to fill the point mask to let them transparent!
[142]161; Comment: Nevertheless, we trace mask's contour.
[231]162;
[142]163; @keyword USETRI
164; To force using triangulation.
[231]165;
166; @keyword MASKFILL
167;
[142]168; @keyword DUPLICATE
[2]169;
[142]170; @keyword DECIMATETRI
[2]171;
[142]172; @uses
173; common.pro
[2]174;
[231]175; @history
[157]176; Sebastien Masson (smasson\@lodyc.jussieu.fr)
[142]177; 7/1999
178; Sebastien Masson 08/02/2000 checkfield and
179; notri keyword (or triangule = -1) .
[231]180;
[142]181; @version
182; $Id$
[231]183;
184; @todo seb
[142]185; Changer param "tab1",keyword "REMPLI", "UNLABSUR", "UNSUR2",
186; "UNVECTSUR"+ quelques trucs dans la routine.
[2]187;-
[231]188;
189PRO plt, tab1, param2, param3, param4, REALCONT = realcont, CONTOUR = contour $
[67]190         , INTERVALLE = intervalle, INV = inv, GRIDTYPE = gridtype, BOXZOOM = boxzoom $
191         , CONTINTERVALLE = contintervalle, LABEL = label, CONTLABEL = contlabel $
192         , STYLE = style, CONTMAX = contmax, CONTMIN = contmin, NLEVEL = nlevel, CONTNLEVEL = contnlevel $
193         , VECTEUR = vecteur, MAP = map, MININ = minin, MAXIN = maxin $
194         , CONT_NOFILL = cont_nofill $
195         , USETRI = usetri, NOTRI = notri, MASKFILL = maskfill $
196         , DUPLICATE = duplicate, STRICTFILL = strictfill, OVERPLOT = overplot $
197         , DECIMATETRI = decimatetri,  LABMAP = labmap, _extra = ex
[114]198;
199  compile_opt idl2, strictarrsubs
200;
[67]201@cm_4mesh
202@cm_4data
203  IF NOT keyword_set(key_forgetold) THEN BEGIN
204@updatenew
205@updatekwd
206  ENDIF
[2]207;--------------------------------------------------------------
[142]208  tempsun = systime(1)          ; To key_performance
[2]209;--------------------------------------------------------------
210;--------------------------------------------------------------
[142]211; I) preparation of the graphic environment and small verifications.
[2]212;--------------------------------------------------------------
213;--------------------------------------------------------------
[142]214; I1) verification of the grid's type associated to tab1.
[2]215;--------------------------------------------------------------
[67]216  if keyword_set(gridtype) then vargrid = gridtype
217  if keyword_set(vecteur) AND (NOT keyword_set(gridtype)) then BEGIN
218    vargrid = litchamp(tab1, /grid)
219    if vargrid eq '' then BEGIN
220      vargrid = xquestion('What is the grid associated to the data to contour?', 'T', /chkwidget)
221      vargrid = strupcase(vargrid)
222    endif
223  ENDIF
[2]224;--------------------------------------------------------------
[142]225; I2) Reading of the field and checkup.
[2]226;--------------------------------------------------------------
[67]227  if keyword_set(boxzoom) AND n_elements(contour) ne 4  then BEGIN
228    savedbox = 1b
229    saveboxparam, 'boxparam4plt.dat'
230  END
[142]231  if n_elements(param2) NE 0 then type = param2
232  if n_elements(param3) NE 0 then min = param3
233  if n_elements(param4) NE 0 then max = param4
[231]234  if n_elements(minin) NE 0 then min = minin
[67]235  if n_elements(maxin) NE 0 then max = maxin
236  checktypeminmax, 'plt', TYPE = type, MIN = min, MAX = max, _extra = ex
[231]237  z2d = checkfield(tab1, 'plt', TYPE = type, BOXZOOM = boxzoom, DIREC = direc, VECTEUR = vecteur, _extra = ex)
238  if z2d[0] EQ -1 then BEGIN
[67]239    IF keyword_set(savedbox) THEN restoreboxparam, 'boxparam4plt.dat'
240    return
241  ENDIF
242  IF n_elements(usetri) EQ 0 THEN BEGIN
243; do we have holes in the triangulation?
[231]244    holeintri = n_elements(triangles_list)/3 LT (jpi-1+keyword_set(key_periodic))*(jpj-1)*2 ; -> the triangulation must be used to draw the field
[67]245; do we have a triangulation?
246    wehavetri = triangles_list[0] NE -1 ; -> the triangulation must be used to draw the continents
247; if we make a map, are we periodic and nx=jpi ?
248    CASE strupcase(vargrid) OF
249      'T':nx = nxt
250      'W':nx = nxt
251      'U':nx = nxu
252      'V':nx = nxv
253      'F':nx = nxf
254    ENDCASE
[231]255    mapperio = keyword_set(map)*keyword_set(key_periodic)*(nx eq jpi)
[67]256    usetri = (wehavetri*(wehavetri+holeintri+mapperio+keyword_set(key_irregular))) < 2
257    IF n_elements(notri) NE 0 THEN usetri = 0 > (2-notri)
258  ENDIF
[2]259;--------------------------------------------------------------
[142]260; I3) reinitialization. !p.,!x.,!y.
261; Comment: We do not reinitializate when we call bck plt in loop to use contour.
[2]262;--------------------------------------------------------------
[67]263  if n_elements(contour) ne 4 AND NOT keyword_set(overplot) then reinitplt, /z, /invert
[2]264;---------------------------------------------------------------
[142]265; I4) attribution of the mask and of longitude and latitude arrays.
[2]266;---------------------------------------------------------------
[67]267  IF  strupcase(vargrid) EQ 'W' THEN profond = firstzw NE 0 $
[231]268  ELSE profond = firstzt NE 0
[67]269; do we need to extract now the triangulation that will be use for
[142]270; contouring the field?
[67]271  if keyword_set(profond) $
272    OR (usetri EQ 0 AND (vargrid EQ 'T' OR vargrid EQ 'W')) $
273    OR (usetri NE 2 AND (vargrid NE 'T' AND vargrid NE 'W')) THEN BEGIN
274    grille, mask, glam, gphi, gdep, nx, ny, nz $
275      , firstx, firsty, firstz, lastx $
276      , lasty, lastz, /forplt, _extra = ex
277  ENDIF ELSE BEGIN
278    grille, mask, glam, gphi, gdep, nx, ny, nz $
279      , firstx, firsty, firstz, lastx, lasty $
280      , lastz, TRI = trifield, /forplt, _extra = ex
281  ENDELSE
[2]282;----------------------------------------------------------------------------
[231]283; I5 determination of the mi:min and of the ma:max of z2d in the same way
[142]284; as max: max and min: min for the drawing.
[2]285;-----------------------------------------------------------------------------
[67]286  masknan = finite(z2d)
[231]287  nan = total(masknan) NE n_elements(z2d)
[142]288; Do we need to do an autoscale ???
[67]289  autoscale = testvar(var = min) EQ testvar(var = max) AND NOT keyword_set(intervalle)
290  determineminmax, z2d, mask, mi, ma, glam, gphi, MININ = min, MAXIN = max $
291    , nan = nan, INTERVALLE = intervalle, usetri = usetri, _extra = ex
292  if z2d[0] EQ -1 THEN GOTO, sortie
[142]293; We do an autoscale if needed.
[67]294  if autoscale then autoscale, min, max, intervalle
[2]295;--------------------------------------------------------------
296;--------------------------------------------------------------
[231]297; II) We put the drawing in its place on the window or the page
[142]298; and possible opening of the window or of the page.
[2]299;--------------------------------------------------------------
[67]300  if n_elements(contour) NE 4 AND NOT keyword_set(overplot) THEN $
301    placedessin, 'plt', posfenetre, posbar $
[2]302    , CONTOUR = contour, VECTEUR = vecteur, MAP = map, DIREC = direc, _extra = ex
303;--------------------------------------------------------------
304;--------------------------------------------------------------
[142]305; III) Covering of the drawing (labels,style,axis)
[2]306;------------------------------------------------------------
307;--------------------------------------------------------------
308;
309;--------------------------------------------------------------
[142]310; III1) Choice of labels
[2]311;-----------------------------------------------------------
[67]312  if keyword_set(intervalle) AND NOT keyword_set(label) then label = 1
313  if keyword_set(label) eq 0 then cas = 0 else cas = label
314  label, cas, min, max, ncontour, level_z2d, colnumb, NLEVEL = nlevel $
315    , INTERVALLE = intervalle, strictfill = strictfill
[2]316;--------------------------------------------------------------
[142]317; III2) Choice of style
[2]318;-----------------------------------------------------------
[67]319  if not keyword_set(style) then style = 0
320  style, style, level_z2d, linestyle, thick
321  if keyword_set(inv) then colnumb = reverse(colnumb)
[2]322;-----------------------------------------------------------
[142]323; III3) Definition of axis
[2]324;----------------------------------------------------------
[67]325  if NOT keyword_set(overplot) THEN axe, 'xy', _EXTRA = ex
[2]326;--------------------------------------------------------------
327;--------------------------------------------------------------
328;------------------------------------------------------------
[142]329; IV) Drawing
[2]330;--------------------------------------------------------------
331;--------------------------------------------------------------
[67]332;--------------------------------------------------------------
[142]333; extrapolation of field on lands and setup of min, max values
[67]334;--------------------------------------------------------------
335  if keyword_set(nan) then begin
336    z2d[where(masknan EQ 0)] = max
337  ENDIF ELSE masknan = 1
338; filling the mask values
339; we fill only masknan or we fill mask*masknan?
340  IF keyword_set(nan) AND keyword_set(cont_nofill) THEN $
341    z2d = remplit(z2d, nite = 1+(vargrid NE 'T' AND vargrid NE 'W') $
342                  , mask = masknan, _extra = ex) $
343   ELSE z2d = remplit(z2d, nite = (1 + (vargrid NE 'T' AND vargrid NE 'W') $
344                                   + keyword_set(nan)) $
345                      *(1-keyword_set(cont_nofill)) $
346                      *(1-(n_elements(maskfill) NE 0)) $
347                , mask = mask*masknan, _extra = ex)
348  IF keyword_set(strictfill) EQ 0 AND n_elements(maskfill) EQ 0 then $
349    z2d = min > z2d <  max
350  if n_elements(maskfill) NE 0 then BEGIN
351    z2d = temporary(z2d)*mask*masknan
352    if maskfill NE 0 then z2d = temporary(z2d) + maskfill*(1-mask*masknan)
353  ENDIF
354;----------------------------------------------------------
[171]355; check mask and triangulation according to the grid type and NaN
[231]356;----------------------------------------------------------
357; If we make a drawing at depth, we must redefine a triangulation
[171]358; on the zoom because the land/sea mask at depth may differ from
359; the one at the surface.
360; As this triangulation will be used to draw the land/sea mask with
361; the good shape, we use tmask to compute it even if the data are not
362; located on the T grid
[67]363  if (keyword_set(profond) OR keyword_set(cont_nofill)) $
364    AND ((usetri GE 1 AND (vargrid EQ 'T' OR vargrid EQ 'W')) $
365         OR (usetri EQ 2 AND (vargrid NE 'T' AND vargrid NE 'W'))) then BEGIN
366    trifield = triangule(tmask[firstx:lastx, firsty:lasty, firstz] $
367                         , coinmonte = coinmontemask $
368                         , coindescend = coindescendmask $
369                         , keep_cont = cont_nofill, _extra = ex)
370    indicezoommask = (lindgen(jpi, jpj))[firstx:lastx, firsty:lasty]
371  ENDIF
372; triangulation for nan mask
373  if keyword_set(nan) then BEGIN
[171]374    trinan = triangule(masknan, /keep_cont, coinmonte = coinmontenan $
[67]375                       , coindescend = coindescendnan)
376    indicezoomnan = (lindgen(jpi, jpj))[firstx:lastx, firsty:lasty]
377  ENDIF
378;  IF n_elements(twin_corners_up) EQ 0 THEN coinmontemask = -1 ELSE coinmontemask = twin_corners_up
379;  IF n_elements(twin_corners_dn) EQ 0 THEN coindescendmask =  -1 ELSE coindescendmask = twin_corners_dn
[231]380  if vargrid EQ 'T' OR vargrid EQ 'W' then BEGIN
[67]381    glammsk = glam
382    gphimsk = gphi
[171]383  ENDIF ELSE begin
[231]384; decoupe terre: To draw the coast in a clean way, we try to take additionally
[171]385; points to draw the land. Like that, we do not see gap between T and U/V/F grid.
386; It is what decoupeterre do. We also redefine trimsk.
387    maskorg = mask
[67]388    decoupeterre, mask, glammsk, gphimsk, type = 'xy'$
389      , TRI = trimsk, usetri = usetri, indicezoom = indicezoommask $
390      , coinmonte = coinmontemask, coindescend = coindescendmask $
391      , _EXTRA = ex
392  ENDELSE
[2]393;------------------------------------------------------------
[142]394; IV1) Choice of type of drawing.
[2]395;------------------------------------------------------------
[67]396  typetrace = 'classique'
[231]397  if keyword_set(map) AND key_onearth then BEGIN
[142]398; Call of  mapset when we want to do projections.
[67]399    IF n_elements(map) NE 3 THEN map = [0, ((lon1+lon2)/2.) MOD 360, 0]
400    typetrace = 'projection'
401    map_lat = map[0]
402    map_lon = map[1]
403    map_rot = map[2]
404    if chkstru(ex, 'TITLE') then begin
[231]405      maptitre = ex.title
[67]406      ex.title = ''
407    endif
408    map_set, map_lat, map_lon, map_rot, _extra = ex, position = posfenetre, /iso $
409      , limit = [lat1, lon1, lat2, lon2], /noborder
410    if n_elements(maptitre) ne 0 then ex.title = maptitre
411    if n_elements(trifield) GE 2 then trifield = ciseauxtri(trifield, glam, gphi, _EXTRA = ex)
412    if n_elements(trimsk) GE 2 then trimsk = ciseauxtri(trimsk, glammsk, gphimsk, _EXTRA = ex)
[231]413    if n_elements(trinan) GE 2 then BEGIN
[67]414      trinan = ciseauxtri(trinan, glam, gphi, _EXTRA = ex)
415      if trinan[0] EQ -1 then undefine, trinan
416    endif
[231]417  ENDIF ELSE BEGIN
[142]418;To axes of coordinates be considerated.
[67]419    if !x.type EQ 0 AND n_elements(contour) LE 4 then $
[231]420      plot, [0], [0], /nodata, xstyle = 5, ystyle = 5, title = '', subtitle = '', /noerase
[67]421    if keyword_set(key_periodic) then BEGIN
[231]422; In this case, triangulation is closed in x and cover all the sphere.
[142]423; We have to cut it at the level where we cut the sphere to make the drawing.
[67]424      if n_elements(trifield) GE 2 then trifield = ciseauxtri(trifield, glam, gphi, _EXTRA = ex)
425      if n_elements(trimsk) GE 2 then trimsk = ciseauxtri(trimsk, glammsk, gphimsk, _EXTRA = ex)
426      if n_elements(trinan) GE 2 then trinan = ciseauxtri(trinan, glam, gphi, _EXTRA = ex)
427    ENDIF
428  endelse
[2]429;------------------------------------------------------------
[142]430; IV2) Contours and coloring
[2]431;------------------------------------------------------------
[67]432  if keyword_set(duplicate)  then BEGIN
[2]433; pour marina uniquement ATTENTION C'EST TRES MAL CODE
[67]434    lon = glam[*, 0]
435    decalage = max(lon)-min(lon)+(lon-shift(lon, 1))[n_elements(lon)-1]
436    !x.range[1] = !x.range[1]+(duplicate-1)*decalage
437    for i = 1, duplicate-1 do BEGIN
438      z2d = [z2d, z2d]
439      gphi = [gphi, gphi]
440      mask = [mask, mask]
441      gphimsk = [gphimsk, gphimsk]
442      glam = [glam, glam+i*decalage]
443      glammsk = [glammsk, glammsk]
444    ENDFOR
445  endif
446; save, glam, gphi, trifield, file = 'tri.dat'
447  if keyword_set(decimatetri) then BEGIN
[142]448    tempsdeux = systime(1)      ; To key_performance
[67]449    IF n_elements(trimsk) EQ 0 THEN trimsk = trifield
450    Verts = transpose([[(temporary(glam))[*]], [(temporary(gphi))[*]] $
451                       , [(temporary(z2d))[*]]])
452    Conn = [replicate(3, 1, n_elements(trifield)/3), trifield]
453    Result = mesh_decimate(temporary(verts), (temporary(Conn))[*] $
454                           , Connout, vertices = Vertsout $
455                           , percent_vertices = decimatetri)
[231]456    connout = reform(connout, 4,  n_elements(connout)/4, /over)
[67]457    trifield = (temporary(connout))[1:3, *]
458    glam = reform(Vertsout[0, *])
459    gphi = reform(Vertsout[1, *])
460    z2d = reform(Vertsout[2, *])
461    undefine, Vertsout
[231]462    print, 'temps decimatetri', systime(1)-tempsdeux
[67]463  ENDIF
464  pltbase, z2d, glam, gphi $
465    , mask, glammsk, gphimsk, trichamp = trifield, trimsk = trimsk, /forplt $
466    , level_z2d, colnumb, contour = contour, usetri = usetri, realcont = realcont $
[2]467    , overplot = keyword_set(overplot)+keyword_set(map) $
[67]468    , c_linestyle = linestyle, c_labels = 1-(indgen(n_elements(level_z2d)) MOD 2) $
469    , c_thick = thick, cont_nofill = cont_nofill, nan = nan $
[2]470    , coinmontemask = coinmontemask, coindescendmask = coindescendmask $
471    , coinmontenan = coinmontenan, coindescendnan = coindescendnan $
472    , indicezoommask = indicezoommask, indicezoomnan = indicezoomnan $
[171]473    , maskorg = maskorg, masknan = masknan, trinan = trinan, _extra = ex
[2]474;------------------------------------------------------------
[142]475; IV3) Recall of plt in loop when contour is activated.
[2]476;------------------------------------------------------------
[67]477  if n_elements(contour) eq 4 then BEGIN ; c''est la 2eme fois que je passe ds pltt
478    contour = {mietma:[mi, ma], unit:varunit, inter:intervalle} ; je renvoie le min, le max et l''unite
479    return
480  endif
[231]481  if keyword_set(contour) THEN BEGIN
[67]482    pourlegende = [1, 1, 1, 1]
483    oldattributs = saveatt()
484    oldcolnumb = colnumb
485    plt, contour, contmin, contmax, CONTOUR = pourlegende, /NOERASE $
486      , USETRI = usetri $
487      , INTERVALLE = contintervalle, LABEL = contlabel, STYLE = style $
488      , NLEVEL = contnlevel, DUPLICATE = duplicate, STRICTFILL = strictfill $
489      , MASKFILL = maskfill, _extra = ex
490    restoreatt, oldattributs
491    colnumb = oldcolnumb
492  ENDIF
[2]493;------------------------------------------------------------
494;------------------------------------------------------------
[142]495; V) Small functions
[2]496;------------------------------------------------------------
497;------------------------------------------------------------
498;
499;------------------------------------------------------------
[142]500; V1) Possible add of vectors in double exposure.
[2]501;------------------------------------------------------------
[231]502  if keyword_set(vecteur) then BEGIN
[67]503    oldattributs = saveatt()
504    ajoutvect, vecteur, vectlegende, _extra = ex
505    restoreatt, oldattributs
506  ENDIF
[2]507;
[67]508  if keyword_set(overplot) then GOTO, fini
[2]509;------------------------------------------------------------
[142]510; V2) Trace the line of change of date, the equator and the greenwich meridian, ...
[2]511;------------------------------------------------------------
[67]512  if NOT keyword_set(map) AND key_onearth then meridienparallele, 'xy'
[2]513;------------------------------------------------------------
[142]514; V3) To trace IDL's continents
[2]515;------------------------------------------------------------
[67]516  if keyword_set(realcont) then BEGIN
[2]517; si noease est passe de _extra, on s''assure qu''il est a 1
[67]518    if chkstru(ex, 'NOERASE') then begin
519      oldnoerase = ex.noerase
520      ex.noerase = 1
521    ENDIF
522    if chkstru(ex, 'coast_thick') then mlinethick = ex.coast_thick ELSE mlinethick = 1
523    if chkstru(ex, 'coast_color') then mcolor = ex.coast_color ELSE mcolor = 0
524    IF NOT keyword_set(map) THEN $
525      map_set, 0, ((lon1+lon2)/2.) MOD 360, 0, position = posfenetre $
526      ,  limit = [lat1, lon1, lat2, lon2], /NOERASE, /noborder, color = 0 $
527      , _extra = ex
528    if realcont NE 2 AND NOT keyword_set(cont_nofill) then BEGIN
529      if chkstru(ex, 'cont_color') then cntcol = ex.coast_color $
530      ELSE cntcol = (!d.n_colors - 1) <  255
531      map_continents, /fill_continents, color = cntcol, _extra = ex, /noerase
532    ENDIF
533    map_continents, /continents, color = mcolor, MLINETHICK = mlinethick, /noerase, _extra = ex
534    if chkstru(ex, 'NOERASE') THEN  ex.noerase = oldnoerase
535  ENDIF
[2]536;------------------------------------------------------------
[142]537; V4) caption + display of these.
[2]538;------------------------------------------------------------
[67]539  legende, mi, ma, 'xy', CONTOUR = pourlegende, VECTLEGENDE = vectlegende, INTERVALLE = intervalle, DIREC = direc, _EXTRA = ex
[231]540  if n_elements(ex) NE 0 then BEGIN
[142]541; To keep frame's axes in black
[67]542    if (where(tag_names(ex) EQ 'COLOR'))[0] NE -1 then ex.COLOR = coast_color
543  endif
[2]544;
[67]545  case typetrace of
546    'classique': $
547      plot, [0], [0], /nodata, /noerase, color = 0, xstyle = 1, ystyle = 1, _extra = ex
[231]548    'projection': BEGIN
[67]549      if chkstru(ex, 'NOERASE') then begin
550        oldnoerase = ex.noerase
551        ex.noerase = 1
552      endif
553      if chkstru(ex, 'SUBTITLE') then !p.subtitle = ex.SUBTITLE
554      if n_elements(maptitre) ne 0 then ex.title = maptitre
555      map_set, map_lat, map_lon, map_rot, _extra = ex, /iso, limit = [lat1, lon1, lat2, lon2]$
556        , /NOERASE, /noborder, title = !p.title, color = 0
[2]557;         map_proj_info, numproj, /current
558;         map_proj_info, numproj, name = nomproj
559;         if nomproj EQ 'Mercator' OR nomproj EQ 'Cylindrical' OR nomproj EQ 'LambertConic'$
560;          OR nomproj EQ 'Gnomic' OR nomproj EQ 'AlbersEqualAreaConic' $
561;          OR nomproj EQ 'TransverseMercator' OR nomproj EQ 'MillerCylindrical' $
562;          OR nomproj EQ 'LambertConicEllipsoid' then map_grid, box_axes=1,latdel=10,londel=10 $
[67]563;         ELSE map_grid, charsize = 0.75, /label,latalign = 1,lonalign
564;         = 1,latdel = 10,londel = 30
565      IF n_elements(labmap) EQ 0 THEN labmap = 1
566      map_grid, charsize = 0.75, label = labmap, latalign = 1, lonalign = 1, latdel = 10, londel = 30, color = 0, _extra = ex
567      if chkstru(ex, 'NOERASE') THEN  ex.noerase = oldnoerase
568    end
569  endcase
[2]570;------------------------------------------------------------
[142]571; V5) Colorbar.
[2]572;------------------------------------------------------------
[67]573  colnumb = colnumb[0:ncontour-1-keyword_set(strictfill)]
574  barrecouleur, colnumb, min,  max, (ncontour-keyword_set(strictfill))/2 $
575                , position = posbar, _extra = ex
[2]576;------------------------------------------------------------
577;------------------------------------------------------------
[142]578; VI) Possible print.
[2]579;------------------------------------------------------------
580;------------------------------------------------------------
581fini:
[67]582  terminedessin, _extra = ex
[2]583;------------------------------------------------------------
584;------------------------------------------------------------
585sortie:
[67]586  if keyword_set(savedbox) THEN restoreboxparam, 'boxparam4plt.dat'
[231]587  if keyword_set(key_performance) NE 0 THEN print, 'temps plt', systime(1)-tempsun
[2]588;------------------------------------------------------------
589;------------------------------------------------------------
[67]590  return
[2]591end
592
593
Note: See TracBrowser for help on using the repository browser.