Changeset 418

Timestamp:

Nov 23, 2011, 4:39:01 PM (13 years ago)

Author:

acolaitis

Message:

PYTHON. Moved the nolow mode to an option (--nolow) and re-arranged the way -H worked in planetoplot, by introducing a new function (hole_bounds). Bounds is called in normal cases and hole_bounds in -H cases. Minor corrections in mcs.py for missing values comming from zrecast.

Location:

trunk/UTIL/PYTHON

Files:

• mcs.py (modified) (7 diffs)
• myplot.py (modified) (1 diff)
• myscript.py (modified) (1 diff)
• planetoplot.py (modified) (4 diffs)
• pp.py (modified) (1 diff)

trunk/UTIL/PYTHON/mcs.py

@@ -60 +60 @@
                     fields  = varznames, \
                     predifined = 'mcs')[0]
-
+   else:filename=opt.file
    # Files
 
@@ -122 +122 @@
    for zn in varznames:
        load=getfield(nc,zn)
-       if opt.zrecast:
-          load=np.ma.masked_where(load < -1.e-30,load)
-          load.set_fill_value([np.NaN])
-          load=load.filled()
+       load=np.ma.masked_where(load < -1.e-20,load)
+       load.set_fill_value([np.NaN])
+       load=load.filled()
+       load=np.ma.masked_invalid(load)
+       load.set_fill_value([np.NaN])
+       load=load.filled()
        varz.append(load)
+       load=0.
        print "Found: "+zn+" with dimensions: "
        print np.array(varz[n]).shape
@@ -193 +196 @@
               varchk[vv,:,:,:,:]=variable[istart:istop+1,:,:,:]
            vv=vv+1
+       varchk=np.ma.masked_invalid(varchk)
+       varchk.set_fill_value([np.NaN])
        varchktime=time[istart:istop+1]
        ndays=np.floor(varchktime[len(varchktime)-1])-np.floor(varchktime[0])
@@ -199 +204 @@
        ntmichk=ntimemin[i,:,:]
        ntmachk=ntimemax[i,:,:]
-       dtmichk.set_fill_value(np.NaN)
-       dtmachk.set_fill_value(np.NaN)
-       ntmichk.set_fill_value(np.NaN)
-       ntmachk.set_fill_value(np.NaN)
+       dtmichk.set_fill_value([np.NaN])
+       dtmachk.set_fill_value([np.NaN])
+       ntmichk.set_fill_value([np.NaN])
+       ntmachk.set_fill_value([np.NaN])
        dtmichk=dtmichk.filled()
        dtmachk=dtmachk.filled()
@@ -215 +220 @@
        vartmpnight=np.zeros([nzvar,ndays,nz,ny,nx])
        vartmpday=np.zeros([nzvar,ndays,nz,ny,nx])
+       vartmpnight=np.ma.masked_invalid(vartmpnight)
+       vartmpday=np.ma.masked_invalid(vartmpday)
+       vartmpnight.set_fill_value([np.NaN])
+       vartmpday.set_fill_value([np.NaN])
+
        nd=0
        print "     .time indices MCS grid -> diagfi grid."
@@ -239 +249 @@
                        vartmpday[vv,nd,:,iy,ix]=np.NaN
                    elif nitdtmichk > nitdtmachk:
-                       vartmpday[vv,nd,:,iy,ix]=(np.mean(varchk[vv,daystart+nitdtmichk:daystop+1,:,iy,ix],axis=0)+np.mean(varchk[vv,daystart:daystart+nitdtmachk+1,:,iy,ix],axis=0))/2.
+                       vartmpday[vv,nd,:,iy,ix]=(np.ma.mean(varchk[vv,daystart+nitdtmichk:daystop+1,:,iy,ix],axis=0)+np.ma.mean(varchk[vv,daystart:daystart+nitdtmachk+1,:,iy,ix],axis=0))/2.
                    else:
-                       vartmpday[vv,nd,:,iy,ix]=np.mean(varchk[vv,daystart+nitdtmichk:daystart+nitdtmachk+1,:,iy,ix],axis=0)
+                       vartmpday[vv,nd,:,iy,ix]=np.ma.mean(varchk[vv,daystart+nitdtmichk:daystart+nitdtmachk+1,:,iy,ix],axis=0)
                    if ((nitntmichk is np.NaN) or (nitntmachk is np.NaN)):
                        vartmpnight[vv,nd,:,iy,ix]=np.NaN
                    elif nitntmichk > nitntmachk:
-                       vartmpnight[vv,nd,:,iy,ix]=(np.mean(varchk[vv,daystart+nitntmichk:daystop+1,:,iy,ix],axis=0)+np.mean(varchk[vv,daystart:daystart+nitntmachk+1,:,iy,ix],axis=0))/2.
+                       vartmpnight[vv,nd,:,iy,ix]=(np.ma.mean(varchk[vv,daystart+nitntmichk:daystop+1,:,iy,ix],axis=0)+np.ma.mean(varchk[vv,daystart:daystart+nitntmachk+1,:,iy,ix],axis=0))/2.
                    else:                                                            
-                       vartmpnight[vv,nd,:,iy,ix]=np.mean(varchk[vv,daystart+nitntmichk:daystart+nitntmachk+1,:,iy,ix],axis=0)
+                       vartmpnight[vv,nd,:,iy,ix]=np.ma.mean(varchk[vv,daystart+nitntmichk:daystart+nitntmachk+1,:,iy,ix],axis=0)
                 iy=iy+1
              ix=ix+1
@@ -256 +266 @@
        vartmpdaymasked=np.ma.masked_invalid(vartmpday)
        vartmpnightmasked=np.ma.masked_invalid(vartmpnight)
+       vartmpdaymasked.set_fill_value([np.NaN])
+       vartmpnightmasked.set_fill_value([np.NaN])
        for vv in np.arange(nzvar):
           vardayout[vv,i,:,:,:]=np.ma.mean(vartmpdaymasked[vv,:,:,:,:],axis=0)

trunk/UTIL/PYTHON/myplot.py

@@ -651 +651 @@
     return what_I_plot
 
+
+## Author : AC
+def hole_bounds(what_I_plot,zevmin,zevmax):
+    import numpy as np
+    zi=0
+    for i in what_I_plot:
+        zj=0
+        for j in i:
+            if ((j < zevmin) or (j > zevmax)):what_I_plot[zi,zj]=np.NaN
+            zj=zj+1
+        zi=zi+1
+
+    return what_I_plot
+
 ## Author: AS
 def zoomset (wlon,wlat,zoom):

trunk/UTIL/PYTHON/myscript.py

@@ -16 +16 @@
     parser.add_option('-M', '--max',    action='append',dest='vmax',     type="float",   default=None,  help='bounding maximum value (append)') 
     parser.add_option('-H', '--hole',   action='store_true',dest='hole',                 default=False, help='holes above max and below min')
+    parser.add_option('--nolow',        action='store_true',dest='nolow',                default=False, help='plot only minima and maxima [false]')
 
     ### VERTICAL INTERPOLATION

trunk/UTIL/PYTHON/planetoplot.py

@@ -52 +52 @@
            yintegral=False,\
            blat=None,\
-           tsat=False):
+           tsat=False,\
+           nolow=False):
 
 
@@ -64 +65 @@
                        fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\
                        zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\
-                       getname,localtime,polarinterv,getsindex,define_axis,determineplot,readslices,bidimfind,getlschar
+                       getname,localtime,polarinterv,getsindex,define_axis,determineplot,readslices,bidimfind,getlschar,hole_bounds
     from mymath import deg,max,min,mean,get_tsat,writeascii
     import matplotlib as mpl
@@ -323 +324 @@
                ##### 2D field
                if len(what_I_plot.shape) is 2:
+                 if hole: what_I_plot = hole_bounds(what_I_plot,zevmin,zevmax)
+                 else: what_I_plot = bounds(what_I_plot,zevmin,zevmax)
+                 if nolow:  what_I_plot = nolow(what_I_plot)
                  if not tile:
-                     if not hole: what_I_plot = bounds(what_I_plot,zevmin,zevmax)
                      #zelevels = np.linspace(zevmin*(1. + 1.e-7),zevmax*(1. - 1.e-7)) #,num=20)
                      zelevels = np.linspace(zevmin,zevmax,num=ticks)
@@ -330 +333 @@
                      elif mapmode == 0:     contourf( x, y, what_I_plot, zelevels, cmap = palette)
                  else:
-                     if hole:  what_I_plot = nolow(what_I_plot)
                      if mapmode == 1:       m.pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )
                      elif mapmode == 0:     pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax )

trunk/UTIL/PYTHON/pp.py

@@ -135 +135 @@
                 ope=opt.operat,fileref=reffile,minop=opt.vminope,maxop=opt.vmaxope,titleref=opt.titref,\
                 invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy,yintegral=opt.column,\
-                blat=opt.blat,tsat=opt.tsat)
+                blat=opt.blat,tsat=opt.tsat,nolow=opt.nolow)
         print 'DONE: '+name
         system("rm -f to_be_erased")