Changeset 1279 for trunk/UTIL


Ignore:
Timestamp:
Jun 2, 2014, 11:03:13 PM (10 years ago)
Author:
aslmd
Message:

PLANETOPLOT. added vorticity and divergence computations. corrected perturbation calculations which was not working in previous commit. added examples of dynamical analysis. simplified handling of save filename and extension. simplified handling of projections (plus changed a few settings). added a makesave method to ppplot objects.

Location:
trunk/UTIL/PYTHON/planetoplot_v2
Files:
5 added
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/UTIL/PYTHON/planetoplot_v2/pp.py

    r1084 r1279  
    5151parser.add_option('-m','--mult',action='store',dest='mult',type="float",default=None,help="multiplicative factor on field")
    5252parser.add_option('-a','--add',action='store',dest='add',type="float",default=None,help="additive factor on field")
    53 parser.add_option('-o','--output',action='store',dest='filename',type="string",default="myplot",help="name of output files")
     53parser.add_option('-o','--output',action='store',dest='filename',type="string",default=None,help="name of output files")
    5454parser.add_option('-d','--directory',action='store',dest='folder',type="string",default="./",help="directory of output files")
    5555parser.add_option('-s','--changetime',action='store',dest='changetime',type="string",default=None,\
    … …  
    6868parser.add_option('-H','--trans',action='store',dest='trans',type="float",default=1.0,help="float for transparency (0=transp,1=opaque)")
    6969parser.add_option('-Z','--logy',action='store_true',dest='logy',default=False,help="set log for vertical axis")
    70 parser.add_option('-O','--save',action='store',dest='out',type="string",default="gui",help="save mode: 'gui' 'png' 'pdf' 'eps' 'svg' 'ps'")
     70parser.add_option('-O','--save',action='store',dest='out',type="string",default=None,help="save mode: 'gui' 'png' 'pdf' 'eps' 'svg' 'ps'")
    7171parser.add_option('-V','--void',action='store_true',dest='void',default=False,help="no colorbar, no title, no labels")
    7272parser.add_option('-U','--units',action='append',dest='units',type="string",default=None,help="units for the field")
    … …  
    8989# -- 2D plot
    9090parser.add_option('-C','--colorbar',action='append',dest='colorbar',type="string",default=None,help="[2D] colormap: http://micropore.files.wordpress.com/2010/06/colormaps.png")
    91 parser.add_option('-P','--proj',action='append',dest='proj',type="string",default=None,help="[2D] map projection: 'cyl' 'npstere' 'spstere' 'ortho' 'moll' 'robin' 'lcc' 'laea' 'merc' 'noproj'")
     91parser.add_option('-P','--proj',action='store',dest='proj',type="string",default=None,help="[2D] map projection: 'cyl' 'npstere' 'spstere' 'ortho' 'moll' 'robin' 'lcc' 'laea' 'merc'")
    9292parser.add_option('-B','--back',action='append',dest='back',type="string",default=None,help='[2D] predefined map background (cf. set_back.txt)')
    9393parser.add_option('-A','--area',action='append',dest='area',type="string",default=None,help='[2D] predefined region of mapping (cf. set_area.txt)')
    … …  
    154154user.folder = opt.folder
    155155user.out = opt.out
    156 # if noproj is given for proj, no map mode
    157 if opt.proj is not None:
    158     if 'noproj' in opt.proj:
    159         user.noproj = True
    160 # if user wants to give a name, we drop the indication of date
    161 if opt.filename != "myplot":
    162     user.includedate = False
     156user.proj = opt.proj
    163157# define plot
    164158user.defineplot()
    … …  
    174168command = ""
    175169for arg in sys.argv: command = command + arg + ' '
    176 try:
     170if opt.filename is not None:
     171  try:
    177172    f = open(opt.folder+'/'+opt.filename+'.sh', 'w')
    178173    f.write(command)   
    179 except IOError:
     174  except IOError:
    180175    print "!! WARNING !! pp.py command not saved. Probably do not have permission to write here."

  • trunk/UTIL/PYTHON/planetoplot_v2/ppclass.py

    r1219 r1279  
    168168                      verbose=False,\
    169169                      quiet=False,\
    170                       noproj=False,\
    171170                      superpose=False,\
    172171                      plotin=None,\
    173172                      forcedimplot=-1,\
    174                       out="gui",\
    175                       filename="myplot",\
     173                      out=None,\
     174                      filename=None,\
    176175                      folder="./",\
    177                       includedate=True,\
     176                      includedate=False,\
    178177                      res=150.,\
    179178                      xlabel=None,ylabel=None,\
    … …  
    233232        self.verbose = verbose
    234233        self.quiet = quiet
    235         self.noproj = noproj
    236234        self.plotin = plotin
    237235        self.superpose = superpose
    … …  
    272270    def printstatus(self):
    273271      if not self.quiet:
    274         if self.filename == "THIS_IS_A_CLONE":
    275             pass
    276         else:
    277             print "**** PPCLASS. Done step: " + self.status
     272        print "**** PPCLASS. Done step: " + self.status
    278273
    279274    # print attributes
    … …  
    303298            self.kind3d = other.kind3d
    304299            self.verbose = other.verbose
    305             self.noproj = other.noproj
     300            self.quiet = other.quiet
    306301            self.plotin = other.plotin
    307302            self.superpose = other.superpose
    … …  
    376371               setattr(the_clone,k,v)
    377372        the_clone.verbose = False
    378         the_clone.filename = "THIS_IS_A_CLONE" # trick to avoid additional outputs
     373        the_clone.quiet = True # trick to avoid additional outputs
    379374        the_clone.define()
    380375        for i in range(self.nfin):
    … …  
    390385                setattr(obj,k,v)
    391386        the_clone.status = "retrieved"
    392         the_clone.filename = self.filename
     387        the_clone.quiet = self.quiet
    393388        return the_clone
    394389
    … …  
    886881                    # specific 2d plot stuff
    887882                    if dp == 2:
     883                        plobj.proj = self.proj
     884                        plobj.svx = self.svx
     885                        plobj.svy = self.svy
    888886                        # -- light grey background for missing values
    889887                        if type(plobj.f).__name__ in 'MaskedArray': plobj.axisbg = '0.75'
    890888                        # -- define if it is a map or a plot
    891889                        plobj.mapmode = ( obj.method_x+obj.method_y == "freefree" \
    892                                        and "grid points" not in obj.name_x \
    893                                        and not self.noproj )
     890                                       and "grid points" not in obj.name_x )
    894891                    # possible user-defined plot settings shared by all plots
    895892                    if self.div is not None: plobj.div = self.div
    … …  
    919916                    # -- 2D specific
    920917                    elif dp == 2:
    921                         if self.proj is not None and not self.noproj: plobj.proj = self.proj
    922918                        if self.vmin is not None: plobj.vmin = self.vmin
    923919                        if self.vmax is not None: plobj.vmax = self.vmax
    924920                        if self.trans is not None: plobj.trans = self.trans
    925921                        if self.back is not None: plobj.back = self.back
    926                         plobj.svx = self.svx
    927                         plobj.svy = self.svy
    928922                    # finally append plot object
    929923                    self.p.append(plobj)
    … …  
    10141008            self.customplot = self.p[0].f.ndim == 2 \
    10151009                        and self.p[0].mapmode == True \
    1016                         and self.p[0].proj not in ["ortho"]
     1010                        and self.p[0].proj not in ["ortho","robin"]
    10171011            if self.customplot:
    10181012                margin = 0.07
    … …  
    10261020            self.howmanyplots = self.plotin.howmanyplots
    10271021            self.customplot = self.plotin.customplot
     1022            self.filename = self.plotin.filename
    10281023        # LOOP on all subplots
    10291024        # NB: cannot use 'for pl in self.p' if self.plotin not None
    … …  
    10871082            mpl.close()
    10881083        # SAVE A PICKLE FILE WITH THE self.p ARRAY OF OBJECTS
    1089         if self.verbose: print "**** Saving session in "+self.filename + ".ppobj"
    1090         savfile = self.folder + "/" + self.filename + ".ppobj"
    1091         try:
    1092             filehandler = open(savfile, 'w')
    1093             pickle.dump(self.p, filehandler)
    1094         except IOError:
    1095             if self.verbose: print "!! WARNING !! Saved object file not written. Probably do not have permission to write here."
     1084        if self.filename is not None:
     1085          if self.verbose: print "**** Saving session in "+self.filename + ".ppobj"
     1086          savfile = self.folder + "/" + self.filename + ".ppobj"
     1087          try:
     1088              filehandler = open(savfile, 'w')
     1089              pickle.dump(self.p, filehandler)
     1090          except IOError:
     1091              if self.verbose: print "!! WARNING !! Saved object file not written. Probably do not have permission to write here."
    10961092        return self
    10971093
    … …  
    12351231                try: self.p[iii].legend = opt.legend[0]
    12361232                except: pass
    1237             ###
    1238             try: self.p[iii].proj = opt.proj[iii]
    1239             except:
    1240                 try: self.p[iii].proj = opt.proj[0]
    1241                 except: pass
    1242             ###
     1233            ####
    12431234            try: self.p[iii].back = opt.back[iii]
    12441235            except:

  • trunk/UTIL/PYTHON/planetoplot_v2/ppcompute.py

    r1219 r1279  
    123123    mm = mean(field,axis=axis)
    124124    # include back the reduced dim
    125     if field.ndim == 4: mm = np.tile(mm,(field.shape[axis],1,1,1))
    126     elif field.ndim == 3: mm = np.tile(mm,(field.shape[axis],1,1))
    127     elif field.ndim == 2: mm = np.tile(mm,(field.shape[axis],1))
    128     # array has right shape but not in good order: fix this.
    129     mm = np.reshape(mm,field.shape)
     125    if field.ndim == 4:
     126      if axis == 0:   mm = mm[np.newaxis,:,:,:]
     127      elif axis == 1: mm = mm[:,np.newaxis,:,:]
     128      elif axis == 2: mm = mm[:,:,np.newaxis,:]
     129      elif axis == 3: mm = mm[:,:,:,np.newaxis]
     130    else:
     131      print "not supported. yet!"
     132      exit()
     133    # repeat the calculated mean on this dimension
     134    mm = np.repeat(mm,field.shape[axis],axis=axis)
    130135    # compute perturbations
    131136    field = field - mm
    … …  
    314319    return time
    315320
     321#########################
     322#### FLOW DIAGNOSTIC ####
     323#########################
     324
     325def divort(u,v,lon,lat,rad):
     326  ####################
     327  # compute divergence and vorticity
     328  # div,vorti = divort(u,v,lon,lat,rad)
     329  ####################
     330  # u: zonal wind
     331  # v: meridional wind
     332  # lon: longitude
     333  # lat: latitude
     334  # rad: radius of the planet
     335  ####################
     336
     337  # compute gradients -- with normalized coordinates
     338  du_y,du_x = np.gradient(u)
     339  dv_y,dv_x = np.gradient(v)
     340
     341  # convert to radian
     342  la = lon*np.pi/180.
     343  ph = lat*np.pi/180.
     344
     345  # ensure 2D arrays for lon,lat
     346  if lon.ndim == 1 and lat.ndim == 1:
     347    [lar,phr] = np.meshgrid(la,ph)
     348  elif lon.ndim == 2 and lat.ndim == 2:
     349    lar,phr = la,ph
     350  else:
     351    print "ppcompute. there is a problem with lat/lon rank. stop."
     352    exit()
     353
     354  # compute normalized gradients for lat/lon grid
     355  dump,dla_x = np.gradient(lar)
     356  dph_y,dump = np.gradient(phr)
     357   
     358  # compute cartesian differential coordinates
     359  dx = dla_x*rad*np.cos(phr)
     360  dy = dph_y*rad
     361   
     362  # eventually compute cartesian derivatives
     363  du_dx = du_x / dx
     364  du_dy = du_y / dy
     365  dv_dx = dv_x / dx
     366  dv_dy = dv_y / dy
     367   
     368  # vorticity
     369  vorti = dv_dx - du_dy
     370  div = du_dx + dv_dy
     371
     372  return div,vorti
     373
     374

  • trunk/UTIL/PYTHON/planetoplot_v2/ppplot.py

    r1211 r1279  
    299299    return
    300300
     301# a necessary addition for people used to matplotlib
     302def show():
     303    mpl.show()
     304
    301305# a generic function to show (GUI) or save a plot (PNG,EPS,PDF,...)
    302306# -------------------------------
    303 def save(mode="gui",filename="plot",folder="./",includedate=True,res=150,custom=False):
    304     if mode != "nothing":
     307def save(mode=None,filename=None,folder="./",includedate=False,res=150,custom=False):
     308    # no filename or no mode set --> graphical user interface
     309    if filename is None or mode is None:
     310      show()
     311    # otherwise --> an image is saved
     312    else:
    305313      # a few settings
    306314      possiblesave = ['eps','ps','svg','png','jpg','pdf']
    307315      # now the main instructions
    308       if mode == "gui":
    309           mpl.show()
    310       elif mode in possiblesave:
     316      if mode in possiblesave:
    311317          ## name of plot
    312318          name = folder+'/'+filename
    … …  
    327333          print "!! ERROR !! File format not supported. Supported: ",possiblesave
    328334
    329 # a necessary addition for people used to matplotlib
    330 def show():
    331     mpl.show()
     335
    332336
    333337##################################
    … …  
    554558            ax.xaxis.set_major_locator(MaxNLocator(self.nxticks))
    555559        else:
    556             print "!! WARNING. in logx mode, ticks are set automatically."
     560            # in logx mode, ticks are set automatically
     561            pass
    557562        if not self.logy:
    558563            ax.yaxis.set_major_locator(MaxNLocator(self.nyticks))
    559564        else:
    560             print "!! WARNING. in logy mode, ticks are set automatically."
     565            # in logy mode, ticks are set automatically
     566            pass
    561567        ## specific modulo labels
    562568        if self.modx is not None:
    … …  
    568574        self.make()
    569575        mpl.show()
     576    # makesave = make + save
     577    # -------------------------------
     578    def makesave(self,mode="png",filename="plot",includedate=False,res=150):
     579        self.make()
     580        save(mode=mode,filename=filename,includedate=includedate,res=res)
     581        close()
     582
    570583
    571584################################
    … …  
    648661        ### PRE-SETTINGS
    649662        ############################################################################################
     663        # if no projection is set, set mapmode to False
     664        if self.proj is None:
     665            self.mapmode = False
    650666        # set dummy xy axis if not defined
    651667        if self.x is None:
    … …  
    727743                ax.xaxis.set_major_locator(MaxNLocator(self.nxticks))
    728744            else:
    729                 print "!! WARNING. in logx mode, ticks are set automatically."
     745                pass
     746                #print "!! WARNING. in logx mode, ticks are set automatically."
    730747            if not self.logy:
    731748                ax.yaxis.set_major_locator(MaxNLocator(self.nyticks))
    732749            else:
    733                 print "!! WARNING. in logy mode, ticks are set automatically."
     750                pass
     751                #print "!! WARNING. in logy mode, ticks are set automatically."
    734752            ## specific modulo labels
    735753            if self.modx is not None:
    … …  
    774792            if self.proj == "cyl":
    775793                format = '%.0f'
     794                partab = np.r_[-90.,-60.,-30.,0.,30.,60.,90.]
    776795            # ... global projections
    777796            elif self.proj in ["ortho","moll","robin"]:
    … …  
    782801                if self.proj in ["robin"]: steplon = 90.
    783802                mertab = np.r_[-360.:360.:steplon]
    784                 partab = np.r_[-90.:90.+steplat:steplat]
     803                #partab = np.r_[-90.:90.+steplat:steplat]
     804                partab = np.r_[-90.,-30.,0.,30.,90.]
    785805                if self.proj == "ortho":
    786806                    merlab = [0,0,0,0] ; parlab = [0,0,0,0]
    … …  
    847867            # contour field. first line contour then shaded contour.
    848868            if self.c is not None:
     869                #zelevelsc = np.arange(900.,1100.,5.)
    849870                objC2 = m.contour(x, y, what_I_contour, \
    850871                            zelevelsc, colors = ccol, linewidths = cline)
    851872                #mpl.clabel(objC2, inline=1, fontsize=10,manual=True,fmt='-%2.0f$^{\circ}$C',colors='r')
     873                #mpl.clabel(objC2, inline=0, fontsize=8, fmt='%.0f',colors='r', inline_spacing=0)
    852874            m.contourf(x, y, what_I_plot, zelevels, cmap = palette, alpha = self.trans, antialiased=True)
    853875        ############################################################################################
    … …  
    856878        if self.trans > 0. and self.showcb:
    857879            ## draw colorbar. settings are different with projections. or if not mapmode.
    858             if not self.mapmode: orientation=zeorientation ; frac = 0.075 ; pad = 0.03 ; lu = 0.5
     880            #if not self.mapmode: orientation=zeorientation ; frac = 0.075 ; pad = 0.03 ; lu = 0.5
     881            if not self.mapmode: orientation=zeorientation ; frac = 0.15 ; pad = 0.04 ; lu = 0.5
    859882            elif self.proj in ['moll']: orientation="horizontal" ; frac = 0.08 ; pad = 0.03 ; lu = 1.0
     883            elif self.proj in ['robin']: orientation="horizontal" ; frac = 0.07 ; pad = 0.1 ; lu = 1.0
    860884            elif self.proj in ['cyl']: orientation="vertical" ; frac = 0.023 ; pad = 0.03 ; lu = 0.5
    861885            else: orientation = zeorientation ; frac = zefrac ; pad = 0.03 ; lu = 0.5
    … …  
    932956        self.make()
    933957        mpl.show()
     958
     959    # makesave = make + save
     960    # -------------------------------
     961    def makesave(self,mode="png",filename="plot",includedate=False,res=150):
     962        self.make()
     963        save(mode=mode,filename=filename,includedate=includedate,res=res)
     964        close()
     965
Note: See TracChangeset for help on using the changeset viewer.