source: trunk/UTIL/PYTHON/planetoplot.py @ 435

Last change on this file since 435 was 435, checked in by aslmd, 13 years ago

GRAPHICS: minor fixes and comments for movies.

  • Property svn:executable set to *
File size: 25.7 KB
Line 
1#######################
2##### PLANETOPLOT #####
3#######################
4
5### A. Spiga     -- LMD -- 06~09/2011 -- General building and mapping capabilities
6### T. Navarro   -- LMD -- 10~11/2011 -- Improved use for GCM and added sections + 1Dplot capabilities
7### A. Colaitis  -- LMD --    11/2011 -- Mostly minor improvements and inter-plot operation capabilities + zrecast interpolation for gcm
8### A. Spiga     -- LMD --    11/2011 -- Extended multivar subplot capabilities + cosmetic changes + general cleaning and tests
9### A. Colaitis  -- LMD --    12/2011 -- Added movie capability [mencoder must be installed]
10
11def planetoplot (namefiles,\
12           level=0,\
13           vertmode=0,\
14           proj=None,\
15           back=None,\
16           target=None,
17           stride=3,\
18           var=None,\
19           colorb="def",\
20           winds=False,\
21           addchar=None,\
22           interv=[0,1],\
23           vmin=None,\
24           vmax=None,\
25           tile=False,\
26           zoom=None,\
27           display=True,\
28           hole=False,\
29           save="gui",\
30           anomaly=False,\
31           var2=None,\
32           ndiv=10,\
33           mult=1.,\
34           zetitle="fill",\
35           slon=None,\
36           slat=None,\
37           svert=None,\
38           stime=None,\
39           outputname=None,\
40           resolution=200,\
41           ope=None,\
42           fileref=None,\
43           minop=0.,\
44           maxop=0.,\
45           titleref="fill",\
46           invert_y=False,\
47           xaxis=[None,None],\
48           yaxis=[None,None],\
49           ylog=False,\
50           yintegral=False,\
51           blat=None,\
52           tsat=False,\
53           flagnolow=False,\
54           mrate=None):
55
56
57    ####################################################################################################################
58    ### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
59
60    #################################
61    ### Load librairies and functions
62    from netCDF4 import Dataset
63    from myplot import getcoord2d,define_proj,makeplotres,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\
64                       fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\
65                       zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\
66                       getname,localtime,polarinterv,getsindex,define_axis,determineplot,readslices,bidimfind,getlschar,hole_bounds
67    from mymath import deg,max,min,mean,get_tsat,writeascii,fig2data,fig2img
68    import matplotlib as mpl
69    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel, title
70    from matplotlib.cm import get_cmap
71    import numpy as np
72    from numpy.core.defchararray import find
73    from videosink import VideoSink
74
75    ################################
76    ### Preliminary stuff
77    ################################
78    print "********************************************"
79    print "********** WELCOME TO PLANETOPLOT **********"
80    print "********************************************"
81    if not isinstance(namefiles, np.ndarray): namefiles = [namefiles]
82    if not isinstance(var, np.ndarray):       var = [var]
83
84    ################################
85    ### Which plot needs to be done?
86    ################################
87    nlon, nlat, nvert, ntime, mapmode, nslices = determineplot(slon, slat, svert, stime)
88    vlon = None ; vlat = None
89    if slon is not None: vlon = slon[0][0]
90    if slat is not None: vlat = slat[0][0]
91    if mrate is not None: mapmode = 0 ### THIS IS DIRTY. MAPMPODE=1 SHOULD WORK WITH MOVIES. WILL BE FIXED.
92    if mapmode == 0:       winds=False
93    elif mapmode == 1:     
94        if svert is None:  svert = readslices(str(level)) ; nvert=1
95    zelen = len(namefiles)*len(var)
96    numplot = zelen*nslices
97    print "********** FILES, SLICES, VARS, TOTAL PLOTS: ", len(namefiles), nslices, len(var), numplot
98    if ope is not None:
99        if fileref is not None:       zelen = zelen + 2
100        elif "var" in ope:            zelen = zelen + 1
101    all_var  = [[]]*zelen ; all_var2  = [[]]*zelen ; all_title = [[]]*zelen ; all_varname = [[]]*zelen ; all_namefile = [[]]*zelen ; all_time = [[]]*zelen
102 
103    #################################################################################################
104    ### Loop over the files + vars initially separated by commas to be plotted on the same figure ###
105    #################################################################################################
106    k = 0 ; firstfile = True
107    for nnn in range(len(namefiles)):
108     for vvv in range(len(var)): 
109
110      print "********** LOOP..... THIS IS SUBPLOT NUMBER.....",k
111
112      ######################
113      ### Load NETCDF object
114      namefile = namefiles[nnn] ; print "********** THE NAMEFILE IS....", namefile
115      nc  = Dataset(namefile)
116
117      ##################################
118      ### Initial checks and definitions
119      ### ... TYPEFILE
120      typefile = whatkindfile(nc)                                 
121      if typefile in ['mesoideal']:   mapmode=0;winds=False
122      if firstfile: print "********** MAPMODE: ", mapmode
123      if firstfile:                 typefile0 = typefile
124      elif typefile != typefile0:   errormess("Not the same kind of files !", [typefile0, typefile])
125      ### ... VAR
126      varname=var[vvv]
127      print "********** THE VAR IS....",varname, var2
128      if varname not in nc.variables: varname = False
129      ### ... WINDS
130      if winds:                                                   
131         [uchar,vchar,metwind] = getwinddef(nc)             
132         if uchar == 'not found': winds = False
133      if not varname and not winds: errormess("please set at least winds or var",printvar=nc.variables)
134      ### ... COORDINATES, could be moved below
135      [lon2d,lat2d] = getcoorddef(nc)
136      ### ... PROJECTION
137      if ((proj == None) and (typefile not in ['mesoideal'])):   proj = getproj(nc)                 
138
139##########################################################
140      if typefile == "gcm":
141          lat = nc.variables["latitude"][:] ; lon = nc.variables["longitude"][:]
142          if "Time" in nc.variables:      time = nc.variables["Time"][:]
143          elif "time" in nc.variables:    time = nc.variables["time"][:]
144          else:                           errormess("no time axis found.")
145          vert = nc.variables["altitude"][:]
146      elif typefile in ['meso','mesoapi','geo','mesoideal']:
147          if vlon is not None or vlat is not None:   indices = bidimfind(lon2d,lat2d,vlon,vlat) ; print '********** INDICES: ', indices
148          if slon is not None: slon[0][0] = indices[0] ; slon[0][1] = indices[0]
149          if slat is not None: slat[0][0] = indices[1] ; slat[0][1] = indices[1]
150          if varname in ['PHTOT','W']:    vertdim='BOTTOM-TOP_PATCH_END_STAG'
151          else:                           vertdim='BOTTOM-TOP_PATCH_END_UNSTAG'
152          if varname in ['V']:  latdim='SOUTH-NORTH_PATCH_END_STAG'
153          else:                 latdim='SOUTH-NORTH_PATCH_END_UNSTAG'
154          if varname in ['U']:  londim='WEST-EAST_PATCH_END_STAG'
155          else:                 londim='WEST-EAST_PATCH_END_UNSTAG'
156          lon = np.arange(0,getattr(nc,londim),1) ; lat = np.arange(0,getattr(nc,latdim),1)
157          if "Times" in nc.variables:time = np.arange(0,len(nc.variables["Times"]),1)
158          elif "Time" in nc.variables:time = np.arange(0,len(nc.variables["Time"]),1)
159          if typefile in ['geo']:   vert = [0.] ; stime = readslices(str(0))
160          else:
161              if vertmode is None:  vertmode=0
162              if vertmode == 0:     vert = np.arange(0,getattr(nc,vertdim),1)
163              else:                 vert = nc.variables["vert"][:]
164       #if firstfile:
165       #   lat0 = lat
166       #elif len(lat0) != len(lat):
167       #   errormess("Not the same latitude lengths !", [len(lat0), len(lat)])
168       #elif sum((lat == lat0) == False) != 0:
169       #   errormess("Not the same latitudes !", [lat,lat0])
170       ## Faire d'autre checks sur les compatibilites entre fichiers!!
171##########################################################
172
173      if firstfile:
174         ##########################
175         ### Define plot boundaries
176         ### todo: possible areas in latinterv in argument (ex: "Far_South_Pole")
177         if proj in ["npstere","spstere"]: [wlon,wlat] = polarinterv(lon2d,lat2d)
178         elif proj in ["lcc","laea"]:      [wlon,wlat] = wrfinterv(lon2d,lat2d)
179         else:                             [wlon,wlat] = simplinterv(lon2d,lat2d)
180         if zoom:                          [wlon,wlat] = zoomset(wlon,wlat,zoom) 
181
182      all_varname[k] = varname
183      all_namefile[k] = namefile
184      all_time[k] = time
185      if var2: all_var2[k] = getfield(nc,var2)
186      ##### SPECIFIC
187      if varname in ["temp","t","T_nadir_nit","T_nadir_day","temp_day","temp_night"] and tsat:
188          tt=getfield(nc,varname) ; print "computing Tsat-T, I ASSUME Z-AXIS IS PRESSURE"
189          if type(tt).__name__=='MaskedArray':  tt.set_fill_value([np.NaN]) ; tinput=tt.filled()
190          else:                                 tinput=tt
191          all_var[k]=get_tsat(vert,tinput,zlon=lon,zlat=lat,zalt=vert,ztime=time)
192      else:
193      ##### GENERAL STUFF HERE
194          all_var[k] = getfield(nc,varname)
195      print "********** all_var[k].shape", all_var[k].shape
196      k += 1
197      firstfile = False
198      #### End of for namefile in namefiles
199
200    ##################################
201    ### Operation on files
202    if ope is not None:
203        print ope
204        if "var" not in ope:
205             if len(var) > 1: errormess("for this operation... please set only one var !")
206             if ope in ["-","+"]:
207                if fileref is not None:   all_var[k] = getfield(Dataset(fileref),all_varname[k-1]) ; all_varname[k] = all_varname[k-1] ; all_time[k] = all_time[k-1] ; all_namefile[k] = all_namefile[k-1]
208                else:                     errormess("fileref is missing!") 
209                if ope == "-":     all_var[k+1]= all_var[k-1] - all_var[k]
210                elif ope == "+":   all_var[k+1]= all_var[k-1] + all_var[k]
211                all_varname[k+1] = all_varname[k] ; all_time[k+1] = all_time[k] ; all_namefile[k+1] = all_namefile[k] ; numplot = numplot+2
212             elif ope in ["cat"]:
213                tab = all_var[0];k = 1
214                while k != len(namefiles)-1: 
215                    tab = np.append(tab,all_var[k],axis=0) ; k += 1
216                all_time[0] = np.arange(0,len(tab),1) ### AS: time reference is too simplistic, should be better
217                all_var[0] = np.array(tab) ; numplot = 1
218             else: errormess(ope+" : non-implemented operation. Check pp.py --help")
219        else:
220             if len(namefiles) > 1: errormess("for this operation... please set only one file !") 
221             if len(var) > 2:       errormess("not sure this works for more than 2 vars... please check.")
222             if   ope in ["div_var"]: all_var[k] = all_var[k-2] / all_var[k-1] ; insert = '_div_'
223             elif ope in ["mul_var"]: all_var[k] = all_var[k-2] * all_var[k-1] ; insert = '_mul_'
224             elif ope in ["add_var"]: all_var[k] = all_var[k-2] + all_var[k-1] ; insert = '_add_'
225             elif ope in ["sub_var"]: all_var[k] = all_var[k-2] - all_var[k-1] ; insert = '_sub_'
226             else:                    errormess(ope+" : non-implemented operation. Check pp.py --help")
227             numplot = numplot + 1 ; all_time[k] = all_time[k-1] ; all_namefile[k] = all_namefile[k-1]
228             all_varname[k] = all_varname[k-2] + insert + all_varname[k-1] 
229
230    ##################################
231    ### Open a figure and set subplots
232    fig = figure()
233    subv,subh = definesubplot( numplot, fig ) 
234 
235    #################################
236    ### Time loop for plotting device
237    nplot = 1;error = False
238    print "********************************************"
239    while error is False:
240       print "********** NPLOT", nplot
241     
242       ### General plot settings
243       if nplot > numplot: break
244       if numplot > 1:  subplot(subv,subh,nplot)
245
246       ### Map projection                   
247       if mapmode == 1:
248           m = define_proj(proj,wlon,wlat,back=back,blat=blat)
249           x, y = m(lon2d, lat2d)
250
251       ####################################################################
252       ## get all indexes to be taken into account for this subplot and then reduce field
253       ## We plot 1) all lon slices 2) all lat slices 3) all vert slices 4) all time slices and then go to the next slice
254       indexlon  = getsindex(slon,(nplot-1)%nlon,lon)
255       indexlat  = getsindex(slat,((nplot-1)//nlon)%nlat,lat)
256       indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert) 
257       if ope is not None:
258           if fileref is not None:      index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%(len(namefiles)+2)  ## OK only 1 var,  see test in the beginning
259           elif "var" in ope:           index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%(len(var)+1)        ## OK only 1 file, see test in the beginning
260       else:                            yeah = len(namefiles)*len(var) ; index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%yeah
261       time = all_time[index_f]
262       if stime is not None:
263           if stime[0][0] < 0:
264               if typefile in ['mesoapi','meso']:
265                   for i in range(len(time)):  time[i] = localtime ( interv[0]+time[i]*interv[1], 0.5*(wlon[0]+wlon[1]) )
266                   print "OK... WORKING WITH LOCAL TIMES"
267               else: errormess("local times not supported for GCM files. not too hard to modify the code though.")
268       if mrate is not None:                 indextime = None 
269       elif mapmode == 1 and stime is None:  indextime = 1
270       else:                                 indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time)
271       ltst = None 
272       if typefile in ['mesoapi','meso'] and indextime is not None:  ltst = localtime ( interv[0]+indextime*interv[1], 0.5*(wlon[0]+wlon[1]) ) 
273       print "********** index lon, lat, vert, time ",indexlon,indexlat,indexvert,indextime
274       ####################################################################
275
276       ticks = ndiv + 1
277
278       #### Contour plot -- to be activated with movies
279       if var2 and mrate is None:
280           what_I_plot, error = reducefield(all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert, yint=yintegral, alt=vert)
281           #what_I_plot = what_I_plot*mult
282           if not error:
283              if mapmode == 1:
284                  if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
285              zevmin, zevmax = calculate_bounds(what_I_plot)
286              zelevels = np.linspace(zevmin,zevmax,ticks) #20)
287              if var2 == 'HGT':  zelevels = np.arange(-10000.,30000.,2000.)
288              if mapmode == 0:
289                  if typefile in ['mesoideal']:    what_I_plot = dumpbdy(what_I_plot,0,stag='W')
290                  itime=indextime
291                  if len(what_I_plot.shape) is 3:itime=[0]
292                  what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
293                        itime,what_I_plot, len(all_var2[index_f].shape),vertmode)
294              ### If we plot a 2-D field
295              if len(what_I_plot.shape) is 2:
296                  #zelevels=[-10.,-8.,-6.,-4.,-2.,0.,2.,4.,6.,8.,10.]
297                  cs = contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5)
298                  #clabel(cs,zelevels,inline=3,fmt='%1.1f',fontsize=7)
299              ### If we plot a 1-D field
300              elif len(what_I_plot.shape) is 1:
301                  plot(what_I_plot,x) 
302           else:
303              errormess("There is an error in reducing field !")
304
305       #### Shaded plot
306       varname = all_varname[index_f]
307       if varname:
308           what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert , yint=yintegral, alt=vert, anomaly=anomaly)
309           what_I_plot = what_I_plot*mult
310           if not error: 
311               fvar = varname
312               if anomaly: fvar = 'anomaly'
313               ##### MAPMODE-SPECIFIC SETTINGS #####
314               if mapmode == 1:
315                   if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
316               elif mapmode == 0:
317                   itime=indextime
318                   if len(what_I_plot.shape) is 3:itime=[0]
319                   what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
320                         itime,what_I_plot, len(all_var[index_f].shape),vertmode)
321                   zxmin, zxmax = xaxis ; zymin, zymax = yaxis
322                   if zxmin is not None: mpl.pyplot.xlim(xmin=zxmin)
323                   if zxmax is not None: mpl.pyplot.xlim(xmax=zxmax)
324                   if zymin is not None: mpl.pyplot.ylim(ymin=zymin) 
325                   if zymax is not None: mpl.pyplot.ylim(ymax=zymax)
326                   if invert_y:     lima,limb = mpl.pyplot.ylim() ; mpl.pyplot.ylim(limb,lima)
327                   if ylog:         mpl.pyplot.semilogy()
328               if (fileref is not None) and (index_f is numplot-1):    zevmin, zevmax = calculate_bounds(what_I_plot,vmin=minop,vmax=maxop)
329               else:                                                   zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
330               if colorb in ["def","nobar"]:   palette = get_cmap(name=defcolorb(fvar.upper()))
331               elif (fileref is not None) and (index_f is numplot-1): palette = get_cmap(name="RdBu_r")
332               else:                           palette = get_cmap(name=colorb)
333               ##### simple 2D field and movies of 2D fields
334               if len(what_I_plot.shape) >= 2:
335                 if (len(what_I_plot.shape) is 3 and mrate is None): 
336                     errormess("3D field in input but not rate specified for movie (use --rate RATE or specify --time TIME). Exit.")
337                 istart=0 
338                 if mrate is not None: iend=len(time)-1
339                 else:                 iend=istart
340                 imov=istart
341                 while imov <= iend:
342                    what_I_plot_frame = what_I_plot
343                    if mrate is not None:      what_I_plot_frame = what_I_plot[imov,:,:] ; print "-> frame ",imov+1
344                    if hole:                   what_I_plot_frame = hole_bounds(what_I_plot_frame,zevmin,zevmax)
345                    else:                      what_I_plot_frame = bounds(what_I_plot_frame,zevmin,zevmax)
346                    if flagnolow:              what_I_plot_frame = nolow(what_I_plot_frame)
347                    if not tile:
348                        #zelevels = np.linspace(zevmin*(1. + 1.e-7),zevmax*(1. - 1.e-7)) #,num=20)
349                        zelevels = np.linspace(zevmin,zevmax,num=ticks)
350                        if imov is 0:          print np.array(x).shape, np.array(y).shape, np.array(what_I_plot_frame).shape
351                        if mapmode == 1:       m.contourf( x, y, what_I_plot_frame, zelevels, cmap = palette)
352                        elif mapmode == 0:     contourf( x, y, what_I_plot_frame, zelevels, cmap = palette)
353                    else:
354                        if mapmode == 1:       m.pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax )
355                        elif mapmode == 0:     pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax )
356                    if colorb != 'nobar':       
357                        if (fileref is not None) and (index_f is numplot-1):   daformat = "%.3f"
358                        else:                                                  daformat = fmtvar(fvar.upper())
359                        colorbar( fraction=0.05,pad=0.03,format=daformat,\
360                                  ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,20])),extend='neither',spacing='proportional' ) 
361                        #mframe=mpl.pyplot.gci().to_rgba(mpl.pyplot.gci().get_array())
362                        #mframe=mpl.pyplot.imshow()
363                        #mframe=fig2data(mpl.pyplot.gcf())
364                        mframe=fig2img(mpl.pyplot.gcf())
365                    if ((mrate is not None) and (imov is 0)):
366                        W,H = mpl.pyplot.gcf().canvas.get_width_height()
367                        video = VideoSink((H,W), "test", rate=mrate, byteorder="rgba")
368                    if mrate is not None: 
369                        video.run(mframe)
370                        mpl.pyplot.close()
371                    imov=imov+1
372                 if mrate is not None: video.close
373               ##### 1D field
374               elif len(what_I_plot.shape) is 1:
375                 plot(x,what_I_plot)
376                 if save == 'txt':  writeascii(np.transpose(what_I_plot),'profile'+str(nplot)+'.txt')
377
378               #### Other cases: (maybe plot 3-D field one day or movie ??)
379               else:
380                 print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!! dimensions: ",what_I_plot.shape
381                 errormess("Are you sure you did not forget to prescribe a dimension ?")
382           else:
383               errormess("There is an error in reducing field !")
384
385       ### Vector plot --- a adapter
386       if winds:
387           vecx, error = reducefield( getfield(nc,uchar), d4=indextime, d3=indexvert , yint=yintegral , alt=vert)
388           vecy, error = reducefield( getfield(nc,vchar), d4=indextime, d3=indexvert , yint=yintegral , alt=vert)
389           #what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert )
390           if not error:
391               if typefile in ['mesoapi','meso']:   
392                   [vecx,vecy] = [dumpbdy(vecx,6,stag=uchar), dumpbdy(vecy,6,stag=vchar)]
393                   key = True
394               elif typefile in ['gcm']:           
395                   key = False
396               if metwind:  [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
397               if varname == False:       colorvec = definecolorvec(back)
398               else:                  colorvec = definecolorvec(colorb)
399               vectorfield(vecx, vecy,\
400                          x, y, stride=stride, csmooth=2,\
401                          #scale=15., factor=300., color=colorvec, key=key)
402                          scale=20., factor=250., color=colorvec, key=key)
403                                            #200.         ## or csmooth=stride
404   
405       ### Next subplot
406       basename = getname(var=varname,winds=winds,anomaly=anomaly)
407       basename = basename + getstralt(nc,level) 
408       if typefile in ['mesoapi','meso']:
409            if slon is not None: basename = basename + "_lon_" + str(int(lon2d[indices[1],indices[0]]))
410            if slat is not None: basename = basename + "_lat_" + str(int(lat2d[indices[1],indices[0]]))
411            plottitle = basename
412            if addchar: 
413                [addchar,gogol,gogol2] = getlschar ( all_namefile[index_f] )
414                plottitle = plottitle + addchar + "_LT"
415            else:        plottitle = plottitle + "_LT"
416            if ltst is not None: plottitle = plottitle + str(ltst)
417       else:
418            if fileref is not None:
419                if index_f is numplot-1:     plottitle = basename+' '+"fig(1) "+ope+" fig(2)"
420                elif index_f is numplot-2:   plottitle = basename+' '+fileref
421                else:                        plottitle = basename+' '+namefiles[0]#index_f]
422            else:                            plottitle = basename+' '+namefiles[0]#index_f]
423       if mult != 1:                         plottitle = str(mult) + "*" + plottitle
424       if zetitle != "fill":                 
425          plottitle = zetitle
426          if titleref is "fill":             titleref=zetitle
427          if fileref is not None:
428             if index_f is numplot-2:        plottitle = titleref
429             if index_f is numplot-1:        plottitle = "fig(1) "+ope+" fig(2)"
430#       if indexlon is not None:      plottitle = plottitle + " lon: " + str(min(lon[indexlon])) +" "+ str(max(lon[indexlon]))
431#       if indexlat is not None:      plottitle = plottitle + " lat: " + str(min(lat[indexlat])) +" "+ str(max(lat[indexlat]))
432#       if indexvert is not None:     plottitle = plottitle + " vert: " + str(min(vert[indexvert])) +" "+ str(max(vert[indexvert]))
433#       if indextime is not None:     plottitle = plottitle + " time: " + str(min(time[indextime])) +" "+ str(max(time[indextime]))
434       title( plottitle )
435       if nplot >= numplot: error = True
436       nplot += 1
437
438 
439
440
441
442
443
444     
445    ##########################################################################
446    ### Save the figure in a file in the data folder or an user-defined folder
447    if outputname is None:
448       if typefile in ['meso','mesoapi']:   prefix = getprefix(nc)
449       elif typefile in ['gcm']:            prefix = 'LMD_GCM_'
450       else:                                prefix = ''
451    ###
452       zeplot = prefix + basename
453       if addchar:         zeplot = zeplot + addchar
454       if numplot <= 0:    zeplot = zeplot + "_LT"+str(abs(numplot))
455    ###
456       if not target:      zeplot = namefile[0:find(namefile,'wrfout')] + zeplot
457       else:               zeplot = target + "/" + zeplot 
458    ###
459    else:
460       zeplot=outputname
461
462    pad_inches_value = 0.35
463    print "********** SAVE ", save
464    if save == 'png': 
465        if display: makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True)  ## a miniature
466        makeplotres(zeplot,res=resolution,pad_inches_value=pad_inches_value,disp=False)
467    elif save in ['eps','svg','pdf']:     makeplotres(zeplot,pad_inches_value=pad_inches_value,disp=False,ext=save)
468    elif save == 'gui':                   show()
469    elif save == 'txt':                   print "Saved results in txt file." 
470    else: 
471         print "INFO: save mode not supported. using gui instead."
472         show()
473
474    ###############
475    ### Now the end
476    return zeplot
Note: See TracBrowser for help on using the repository browser.