Changeset 612

Timestamp:

Apr 4, 2012, 1:59:05 PM (13 years ago)

Author:

acolaitis

Message:

PYTHON
======

New options:

i) --mark lon,lat will superimpose a cross at lon,lat (floats) on a longitude latitude plot (mapmode 1)
ii) --lstyle let you customize linestyles for 1D plots. Works in the same way as --labels. Exemple: --lstyle "--r","-r","--b","-b"
iii) One can now ask for the variable -v slopexy or SLOPEXY to plot (or overplot with --var2 slopexy) the amplitude of the slope in mesoscale simulations
iv) One can now ask for the variable -v deltat or DELTAT to plot the temperature difference between surface temperature and an arbitrary model level. This automatically calls API with "tk,tsurf". One should use this in conjunction with -i 4 -l xxx so that the difference is made between 0m (surface) and xxx m.

Minor corrections:

Verification plot showing longitude latitude of a slice in a HGT map now also requires "display" to be true.
Winds can now be overplotted in --operation - maps.

Location:

trunk/UTIL/PYTHON

Files:

• myplot.py (modified) (4 diffs)
• myscript.py (modified) (2 diffs)
• planetoplot.py (modified) (13 diffs)
• pp.py (modified) (2 diffs)

trunk/UTIL/PYTHON/myplot.py

@@ -105 +105 @@
 
 ## Author: AC
+# Compute the norm of the winds
+# The corresponding variable to call is UV or uvmet (to use api)
 def windamplitude (nc):
     import numpy as np
@@ -123 +125 @@
        zvint=zv
     return np.sqrt(zuint**2 + zvint**2)
+
+## Author: AC
+# Compute the norm of the slope angles
+# The corresponding variable to call is SLOPEXY
+def slopeamplitude (nc):
+    import numpy as np
+    varinfile = nc.variables.keys()
+    if "slopex" in varinfile: zu=getfield(nc,'slopex')
+    elif "SLOPEX" in varinfile: zu=getfield(nc,'SLOPEX')
+    if "slopey" in varinfile: zv=getfield(nc,'slopey')
+    elif "SLOPEY" in varinfile: zv=getfield(nc,'SLOPEY')
+    znt,zny,znx = np.array(zu).shape
+    zuint = np.zeros([znt,zny,znx])
+    zvint = np.zeros([znt,zny,znx])
+    zuint=zu
+    zvint=zv
+    return np.sqrt(zuint**2 + zvint**2)
+
+## Author: AC
+# Compute the temperature difference between surface and first level.
+# API is automatically called to get TSURF and TK.
+# The corresponding variable to call is DELTAT 
+def deltat0t1 (nc):
+    import numpy as np
+    varinfile = nc.variables.keys()
+    if "tsurf" in varinfile: zu=getfield(nc,'tsurf')
+    elif "TSURF" in varinfile: zu=getfield(nc,'TSURF')
+    if "tk" in varinfile: zv=getfield(nc,'tk')
+    elif "TK" in varinfile: zv=getfield(nc,'TK')
+    znt,zny,znx = np.array(zu).shape
+    zuint = np.zeros([znt,zny,znx])
+    zuint=zu - zv[:,0,:,:]
+    return zuint
 
 ## Author: AS + TN + AC
@@ -814 +849 @@
              "TPOT":         "%.0f",\
              "TSURF":        "%.0f",\
+             "U_OUT1":       "%.0f",\
+             "T_OUT1":       "%.0f",\
              "def":          "%.1e",\
              "PTOT":         "%.0f",\
@@ -831 +868 @@
              "UM":           "%.0f",\
              "WIND":         "%.0f",\
+             "UVMET":         "%.0f",\
+             "UV":         "%.0f",\
              "ALBBARE":      "%.2f",\
              "TAU":          "%.1f",\

trunk/UTIL/PYTHON/myscript.py

@@ -39 +39 @@
     parser.add_option('--ylabel',       action='store',dest='ylab',       type="string",  default=None, help='customize the y-axis label')
     parser.add_option('--labels',       action='store',dest='labels',     type="string",  default=None, help='customize 1D curve labels. Str comma-separated. [None]')
+    parser.add_option('--lstyle',       action='store',dest='linestyle',  type="string",  default=None, help='customize 1D curve linestyles. String separated by commas. [None]')
 
     ### SPECIFIC FOR MAPPING [MAPMODE 1]
@@ -48 +49 @@
     parser.add_option('--blat',         action='store',dest='blat',      type="int",     default=None,  help='reference lat (or bounding lat for stere) [computed]')
     parser.add_option('--blon',         action='store',dest='blon',      type="int",     default=None,  help='reference lon [computed]')
+    parser.add_option('--mark',         action='append',dest='mark',  type="string",  default=None, help='superimpose a crossmark at given coords. lon,lat float separated by commas. [None]')
 
     ### SPECIFIC FOR SLICING [MAPMODE 0]

trunk/UTIL/PYTHON/planetoplot.py

@@ -64 +64 @@
            xlab=None,\
            ylab=None,\
-           lbls=None):
+           lbls=None,\
+           lstyle=None,\
+           cross=None):
 
     ####################################################################################################################
@@ -75 +77 @@
                        fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\
                        zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\
-                       getname,localtime,check_localtime,polarinterv,getsindex,define_axis,determineplot,readslices,bidimfind,getlschar,hole_bounds,windamplitude
+                       getname,localtime,check_localtime,polarinterv,getsindex,define_axis,determineplot,readslices,bidimfind,getlschar,hole_bounds,\
+                       windamplitude,slopeamplitude,deltat0t1
     from mymath import deg,max,min,mean,get_tsat,writeascii,fig2data,fig2img
     import matplotlib as mpl
@@ -148 +151 @@
     ### we get the names of variables to be read. in case only one available, we choose this one.
       varname=var[vvv] 
-      if ((varname not in nc.variables) and not ((typefile in ['meso']) and (varname in ['UV','uv','uvmet']))): 
+      if ((varname not in nc.variables) and not ((typefile in ['meso']) and (varname in ['UV','uv','uvmet','slopexy','SLOPEXY','deltat','DELTAT']))): 
           if len(varinfile) == 1:   varname = varinfile[0] 
           else:                     varname = False
@@ -205 +208 @@
           ###### principle: calculate correct indices then repopulate slon and slat
           if slon is not None or slat is not None:
-              if firstfile and save == 'png' and typefile == 'meso' and "HGT" in varinfile:   iwantawhereplot = nc     #show a topo map with a cross on the chosen point
+              if firstfile and save == 'png' and typefile == 'meso' and "HGT" in varinfile and display:   iwantawhereplot = nc     #show a topo map with a cross on the chosen point
               else:                                                                           iwantawhereplot = None   #do not show anything, just select indices
               numlon = 1 ; numlat = 1 
@@ -278 +281 @@
       all_namefile[k] = namefile
       all_time[k] = time
-      if var2: all_var2[k] = getfield(nc,var2)
+      if var2: 
+         if ((var2 in ['slopexy','SLOPEXY']) and (typefile in ['meso']) and (var2 not in nc.variables)): all_var2[k] = slopeamplitude(nc)
+         else: all_var2[k] = getfield(nc,var2)
       if winds: all_windu[k] = getfield(nc,uchar) ; all_windv[k] = getfield(nc,vchar)
       #####################
@@ -293 +298 @@
       elif ((varname in ['UV','uv','uvmet']) and (typefile in ['meso']) and (varname not in nc.variables)):
           all_var[k]=windamplitude(nc)
+      elif ((varname in ['slopexy','SLOPEXY']) and (typefile in ['meso']) and (varname not in nc.variables)):
+          all_var[k]=slopeamplitude(nc)
+      elif ((varname in ['DELTAT','deltat']) and (typefile in ['meso']) and (varname not in nc.variables)):
+          all_var[k]=deltat0t1(nc)
       else:
       ##### GENERAL STUFF HERE
@@ -310 +319 @@
                    if ((all_varname[k-1] in ['UV','uv','uvmet']) and (typefile in ['meso']) and (all_varname[k-1] not in ncref.variables)):
                       all_var[k] = windamplitude(ncref)
+                   elif ((all_varname[k-1] in ['slopexy','SLOPEXY']) and (typefile in ['meso']) and (all_varname[k-1] not in ncref.variables)):
+                      all_var[k] = slopeamplitude(ncref)
+                   elif ((all_varname[k-1] in ['DELTAT','deltat']) and (typefile in ['meso']) and (all_varname[k-1] not in ncref.variables)):
+                      all_var[k]=deltat0t1(ncref)
                    else:
                       all_var[k] = getfield(ncref,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] ; all_var2[k] = all_var2[k-1]
+                   if winds: all_windu[k] = getfield(ncref,uchar) ; all_windv[k] = getfield(ncref,vchar)
                 else: errormess("fileref is missing!") 
                 if ope == "-":     all_var[k+1]= all_var[k-1] - all_var[k]
@@ -321 +335 @@
                     all_var[k+1]= 100.*(all_var[k-1] - masked)/masked
                 all_varname[k+1] = all_varname[k] ; all_time[k+1] = all_time[k] ; all_namefile[k+1] = all_namefile[k] ; all_var2[k+1] = all_var2[k] ; numplot = numplot+2
+                if winds: all_windu[k+1] = all_windu[k-1]-all_windu[k] ; all_windv[k+1] = all_windv[k-1] - all_windv[k]
              elif ope in ["cat"]:
                 tabtime = all_time[0];tab = all_var[0];k = 1
@@ -349 +364 @@
     ### Time loop for plotting device
     nplot = 1 ; error = False 
-    linecycler = cycle(["-","--","-.",":"])
+    if lstyle is not None: linecycler = cycle(lstyle)
+    else: linecycler = cycle(["-","--","-.",":"])
     print "********************************************"
     while error is False:
@@ -492 +508 @@
 
                     if which == "unidim":
-                        lbl = ""
-                        if indexlat is not None:  lbl = lbl + " ix" + str(indexlat[0])
-                        if indexlon is not None:  lbl = lbl + " iy" + str(indexlon[0])
-                        if indexvert is not None: lbl = lbl + " iz" + str(indexvert[0])
-                        if indextime is not None: lbl = lbl + " it" + str(indextime[0])
-                        if lbl == "": lbl = namefiles[index_f]
-                        if lbls is not None: lbl=lbls[nplot-1]
+                        if lbls is not None: lbl=lbls[index_f]
+                        else:
+                           lbl = ""
+                           if indexlat is not None:  lbl = lbl + " ix" + str(indexlat[0])
+                           if indexlon is not None:  lbl = lbl + " iy" + str(indexlon[0])
+                           if indexvert is not None: lbl = lbl + " iz" + str(indexvert[0])
+                           if indextime is not None: lbl = lbl + " it" + str(indextime[0])
+                           if lbl == "": lbl = namefiles[index_f]
 
                         if mrate is not None: x = y  ## because swapaxes...
@@ -523 +540 @@
                             if mapmode == 1:       m.contourf( x, y, what_I_plot_frame, zelevels, cmap = palette, alpha=trans)
                             elif mapmode == 0:     contourf( x, y, what_I_plot_frame, zelevels, cmap = palette, alpha=trans)
+                            if cross is not None and mapmode == 1:
+                               idx,idy=m(cross[0][0],cross[0][1])
+                               mpl.pyplot.plot([idx],[idy],'+k',mew=0.5,ms=18)
                         else:
                             if mapmode == 1:       m.pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans)
@@ -553 +573 @@
                         zevminc, zevmaxc = calculate_bounds(what_I_plot_frame, vmin=min(what_I_plot_frame), vmax=max(what_I_plot_frame))
                         zelevels = np.linspace(zevminc,zevmaxc,ticks/2) #20)
-                        if var2 == 'HGT': zelevels = np.arange(-10000.,30000.,2000.)
+                        if var2 == 'HGT': zelevels = np.arange(-10000.,30000.,1000.)
                         if mapmode == 0:   
                             what_I_plot_frame, x, y = define_axis( lon,lat,vert,time,indexlon,indexlat,indexvert,\
                                                               itime,what_I_plot_frame, len(all_var2[index_f].shape),vertmode )
-                            cs = contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5)#, linestyles='solid')
-                        elif mapmode == 1:  cs = m.contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5)# , linestyles='solid')
+                            cs = contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33, alpha=0.5, linestyles='solid')
+                        elif mapmode == 1:  cs = m.contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33, alpha=0.5, linestyles='solid')
                     if which in ["regular","unidim"]:
 
@@ -645 +665 @@
         if zetitle[0] != "fill": zecb.ax.set_xlabel(zetitle[index_f]) ; zetitle[index_f]=""
 
-      
+
     ##########################################################################
     ### Save the figure in a file in the data folder or an user-defined folder

trunk/UTIL/PYTHON/pp.py

@@ -102 +102 @@
             ### winds or no winds
             if opt.winds            :  zefields = 'uvmet'
+            elif opt.var[j] in ['deltat','DELTAT'] : zefields = 'tk,TSURF'
             else                    :  zefields = ''
             ### var or no var
@@ -170 +171 @@
                 blat=opt.blat,blon=opt.blon,tsat=opt.tsat,flagnolow=opt.nolow,\
                 mrate=opt.rate,mquality=opt.quality,trans=opt.trans,zarea=opt.area,axtime=opt.axtime,\
-                redope=opt.redope,seevar=opt.seevar,xlab=opt.xlab,ylab=opt.ylab,lbls=separatenames(opt.labels))
+                redope=opt.redope,seevar=opt.seevar,xlab=opt.xlab,ylab=opt.ylab,lbls=separatenames(opt.labels),\
+                lstyle=separatenames(opt.linestyle),cross=readslices(opt.mark))
         print 'DONE: '+name
         system("rm -f to_be_erased")