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

Last change on this file since 457 was 457, checked in by acolaitis, 14 years ago

PYTHON.

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