Index: trunk/UTIL/PYTHON/planetoplot.py
===================================================================
--- trunk/UTIL/PYTHON/planetoplot.py	(revision 448)
+++ trunk/UTIL/PYTHON/planetoplot.py	(revision 451)
@@ -8,4 +8,5 @@
 ### A. Spiga     -- LMD -- 11~12/2011 -- Extended multivar subplot capabilities + cosmetic changes + general cleaning and tests
 ### A. Colaitis  -- LMD --    12/2011 -- Added movie capability [mencoder must be installed]
+### A. Spiga     -- LMD --    12/2011 -- Added HTML animated page capability + general tests of consistency [winds, etc...]
 
 def planetoplot (namefiles,\
@@ -50,4 +51,5 @@
            yintegral=False,\
            blat=None,\
+           blon=None,\
            tsat=False,\
            flagnolow=False,\
@@ -68,5 +70,5 @@
     from mymath import deg,max,min,mean,get_tsat,writeascii,fig2data,fig2img
     import matplotlib as mpl
-    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel, title
+    from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel, title, close
     from matplotlib.cm import get_cmap
     import numpy as np
@@ -97,4 +99,5 @@
            stime = readslices(str(0)) ; ntime=1 ## this is a default choice
            print "WELL... nothing about time axis. I took default: first time reference stored in file."
+    if mrate is not None and len(var) > 1: errormess("multivar not allowed in movies. should be fixed soon!")
     zelen = len(namefiles)*len(var)
     numplot = zelen*nslices
@@ -103,5 +106,5 @@
         if fileref is not None:       zelen = zelen + 2
         elif "var" in ope:            zelen = zelen + 1
-    all_var  = [[]]*zelen ; all_var2  = [[]]*zelen ; all_title = [[]]*zelen ; all_varname = [[]]*zelen ; all_namefile = [[]]*zelen ; all_time = [[]]*zelen
+    all_var  = [[]]*zelen ; all_var2  = [[]]*zelen ; all_title = [[]]*zelen ; all_varname = [[]]*zelen ; all_namefile = [[]]*zelen ; all_time = [[]]*zelen ; all_windu = [[]]*zelen ; all_windv = [[]]*zelen
  
     #################################################################################################
@@ -191,4 +194,5 @@
       all_time[k] = time
       if var2: all_var2[k] = getfield(nc,var2)
+      if winds: all_windu[k] = getfield(nc,uchar) ; all_windv[k] = getfield(nc,vchar)
       ##### SPECIFIC
       if varname in ["temp","t","T_nadir_nit","T_nadir_day","temp_day","temp_night"] and tsat:
@@ -219,5 +223,5 @@
              elif ope in ["cat"]:
                 tab = all_var[0];k = 1
-                while k != len(namefiles)-1: 
+                while k != len(namefiles):
                     tab = np.append(tab,all_var[k],axis=0) ; k += 1
                 all_time[0] = np.arange(0,len(tab),1) ### AS: time reference is too simplistic, should be better
@@ -252,5 +256,5 @@
 
        ### Map projection                    
-       if mapmode == 1:     m = define_proj(proj,wlon,wlat,back=back,blat=blat) ; x, y = m(lon2d, lat2d)
+       if mapmode == 1:     m = define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon) ; x, y = m(lon2d, lat2d)
        elif mapmode ==0:    m = None ; x = None ; y = None
 
@@ -264,4 +268,5 @@
            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
            elif "var" in ope:           index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%(len(var)+1)        ## OK only 1 file, see test in the beginning
+           elif "cat" in ope:           index_f = 0
        else:                            yeah = len(namefiles)*len(var) ; index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%yeah
        time = all_time[index_f]
@@ -284,7 +289,10 @@
                                              yint=yintegral, alt=vert, anomaly=anomaly )
            what_I_plot = what_I_plot*mult
-       if var2:      ### what is contoured
+       if var2:      ### what is contoured.
            what_I_plot_contour, error = reducefield( all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert, \
                                              yint=yintegral, alt=vert )
+       if winds:     ### what is plot as vectors.
+           vecx, error = reducefield( all_windu[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert)
+           vecy, error = reducefield( all_windv[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert)
        ####################################################################
 
@@ -325,4 +333,7 @@
                         if mrate is None:                                   what_I_plot_frame = what_I_plot
                         else:                                               what_I_plot_frame = what_I_plot[imov,:,:]
+                        if winds:
+                            if mrate is None:                                   vecx_frame = vecx ; vecy_frame = vecy
+                            else:                                               vecx_frame = vecx[imov,:,:] ; vecy_frame = vecy[imov,:,:]
                     elif which == "contour":  
                         if mrate is None or what_I_plot_contour.ndim < 3:   what_I_plot_frame = what_I_plot_contour
@@ -330,6 +341,6 @@
                     if mrate is not None:     
                         if mapmode == 1: 
-                            m = define_proj(proj,wlon,wlat,back=back,blat=blat)  ## this is dirty, defined above but out of imov loop
-                            x, y = m(lon2d, lat2d)                               ## this is dirty, defined above but out of imov loop
+                            m = define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon)  ## this is dirty, defined above but out of imov loop
+                            x, y = m(lon2d, lat2d)                                         ## this is dirty, defined above but out of imov loop
                     if typefile in ['mesoapi','meso'] and mapmode == 1:   what_I_plot_frame = dumpbdy(what_I_plot_frame,6,condition=True)
                     if typefile in ['mesoideal']:    what_I_plot_frame = dumpbdy(what_I_plot_frame,0,stag='W',condition=dumped_vert_stag)
@@ -362,4 +373,17 @@
                             colorbar( fraction=0.05,pad=0.03,format=daformat,\
                                       ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,20])),extend='neither',spacing='proportional' ) 
+                        if winds:
+                            if typefile in ['mesoapi','meso']:
+                                [vecx_frame,vecy_frame] = [dumpbdy(vecx_frame,6,stag=uchar,condition=True), dumpbdy(vecy_frame,6,stag=vchar,condition=True)]
+                                key = True
+                            elif typefile in ['gcm']:
+                                key = False
+                            if metwind:  [vecx_frame,vecy_frame] = m.rotate_vector(vecx_frame, vecy_frame, lon2d, lat2d)
+                            if var:       colorvec = definecolorvec(back)
+                            else:         colorvec = definecolorvec(colorb)
+                            vectorfield(vecx_frame, vecy_frame, x, y, stride=stride, csmooth=2,\
+                                             #scale=15., factor=300., color=colorvec, key=key)
+                                             scale=20., factor=250., color=colorvec, key=key)
+                                                              #200.         ## or csmooth=stride
                     elif which == "contour":
                         zevminc, zevmaxc = calculate_bounds(what_I_plot_frame)
@@ -375,12 +399,25 @@
                     if which == "regular":
                         if mrate is not None:
-                           figframe=mpl.pyplot.gcf()
-                           if mquality:   figframe.set_dpi(600.)
-                           else:          figframe.set_dpi(200.)
-                           mframe=fig2img(figframe)
-                           if ((mrate is not None) and (imov is 0)):
-                              moviename='test' ;W,H = figframe.canvas.get_width_height()
-                              video = VideoSink((H,W), moviename, rate=mrate, byteorder="rgba")
-                           video.run(mframe) ; mpl.pyplot.close()
+                           ### THIS IS A MENCODER MOVIE
+                           if mrate > 0:
+                             figframe=mpl.pyplot.gcf()
+                             if mquality:   figframe.set_dpi(600.)
+                             else:          figframe.set_dpi(200.)
+                             mframe=fig2img(figframe)
+                             if imov == 0:
+                                moviename='movie' ;W,H = figframe.canvas.get_width_height()
+                                video = VideoSink((H,W), moviename, rate=mrate, byteorder="rgba")
+                             video.run(mframe) ; close()
+                             if imov == iend: video.close()                            
+                           ### THIS IS A WEBPAGE MOVIE 
+                           else:
+                             nameframe = "image"+str(1000+imov)
+                             makeplotres(nameframe,res=100.,disp=False) ; close()
+                             if imov == 0: myfile = open("zepics", 'w')
+                             myfile.write("modImages["+str(imov)+"] = '"+nameframe+"_100.png';"+ '\n')
+                             if imov == iend:
+                                 myfile.write("first_image = 0;"+ '\n')
+                                 myfile.write("last_image = "+str(iend)+";"+ '\n')
+                                 myfile.close()
                         if var2:  which = "contour"
                         imov = imov+1
@@ -388,6 +425,4 @@
                         which = "regular"
 
-                 if mrate is not None: video.close()
-
                ##### 1D field
                elif len(what_I_plot.shape) is 1:
@@ -395,5 +430,5 @@
                  if save == 'txt':  writeascii(np.transpose(what_I_plot),'profile'+str(nplot)+'.txt')
 
-               #### Other cases: (maybe plot 3-D field one day or movie ??)
+               #### Other cases: (maybe plot 3-D field one day ??)
                else:
                  print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!! dimensions: ",what_I_plot.shape
@@ -402,27 +437,8 @@
                errormess("There is an error in reducing field !")
 
-       ### Vector plot --- a adapter
-       if winds:
-           vecx, error = reducefield( getfield(nc,uchar), d4=indextime, d3=indexvert , yint=yintegral , alt=vert)
-           vecy, error = reducefield( getfield(nc,vchar), d4=indextime, d3=indexvert , yint=yintegral , alt=vert)
-           #what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert )
-           if not error:
-               if typefile in ['mesoapi','meso']:    
-                   [vecx,vecy] = [dumpbdy(vecx,6,stag=uchar,condition=True), dumpbdy(vecy,6,stag=vchar,condition=True)]
-                   key = True
-               elif typefile in ['gcm']:            
-                   key = False
-               if metwind:  [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)
-               if varname == False:       colorvec = definecolorvec(back)
-               else:                  colorvec = definecolorvec(colorb)
-               vectorfield(vecx, vecy,\
-                          x, y, stride=stride, csmooth=2,\
-                          #scale=15., factor=300., color=colorvec, key=key)
-                          scale=20., factor=250., color=colorvec, key=key)
-                                            #200.         ## or csmooth=stride
-    
        ### Next subplot
        basename = getname(var=varname,winds=winds,anomaly=anomaly)
-       basename = basename + getstralt(nc,level)  
+       basename = basename + getstralt(nc,level) 
+       if mrate is not None: basename = "movie_" + basename 
        if typefile in ['mesoapi','meso']:
             if slon is not None: basename = basename + "_lon_" + str(int(lon2d[indices[1],indices[0]]))
