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

Last change on this file since 432 was 432, checked in by acolaitis, 13 years ago

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