Changeset 792 for trunk/UTIL

Timestamp:

Sep 21, 2012, 3:03:27 PM (13 years ago)

Author:

acolaitis

Message:

PYTHON GRAPHICS

• Added option --finddevil which finds the steepest pressure drop in the PSFC field of the file, and add a cross in it's position.

This new function is intended to provide lon and lat of a dust devil to planetoplot. This could be used to drive a slice plot at the corresponding location.

• Removed dumpbdy for LES files (no relaxation zone)
• Added --mark mode to tiled plot
• Added routines in mymath to convert back and forth between greyscale Image objects and 2D arrays

Location:

trunk/UTIL/PYTHON

Files:

• mymath.py (modified) (3 diffs)
• myplot.py (modified) (1 diff)
• myscript.py (modified) (1 diff)
• planetoplot.py (modified) (6 diffs)
• pp.py (modified) (1 diff)

trunk/UTIL/PYTHON/mymath.py

@@ -229 +229 @@
     return idx
 
+# Author: A.C.
 def fig2data ( fig ):
     import numpy
@@ -248 +249 @@
     return buf
 
+# Author: A.C.
 def fig2img ( fig ):
     import Image
@@ -260 +262 @@
     w, h, d = buf.shape
     return Image.fromstring( "RGBA", ( w ,h ), buf.tostring( ) )
+
+# Author: A.C.
+# Convert a single layer image object (greyscale) to an array
+def image2array(im):
+    import numpy as np
+    if im.mode not in ("L", "F"):
+        raise ValueError, ("can only convert single-layer images", im.mode)
+    if im.mode == "L":
+        a = np.fromstring(im.tostring(), np.uint8)
+    else:
+        a = np.fromstring(im.tostring(), np.float32)
+    a.shape = im.size[1], im.size[0]
+    return a
+
+# Author: A.C.
+# Convert a 2D array to a single layer image object (greyscale)
+def array2image(a):
+    import numpy as np
+    import Image
+    if a.dtype == np.uint8:
+        mode = "L"
+    elif a.dtype == np.float32:
+        mode = "F"
+    else:
+        raise ValueError, "unsupported image mode"
+    return Image.fromstring(mode, (a.shape[1], a.shape[0]), a.tostring())
+

trunk/UTIL/PYTHON/myplot.py

@@ -1541 +1541 @@
     zt=(zteta+220.)*(zptot/p0)**(r_cp)
     return zt
+
+# Author : A.C.
+# Find the lon and lat index of the dust devil with the largest pressure gradient
+# Steps :
+# 1/ convert the chosen PSFC frame to an image of the PSFC anomaly with respect to the mean
+# 2/ apply the Sobel operator
+# (The Sobel operator performs a 2-D spatial gradient measurement on an image and so emphasizes regions of high spatial frequency that correspond to edges.)
+# 3/ find the maximum of the resulting field 
+# 4/ find the points in a 5 pixel radius around the maximum for which the value of the Sobel transform is greater than half the maximum
+# 5/ define a slab of points encompassing the above selected points, including the potential points 'inside' them (if the above points are a hollow circle for example)
+# 6/ in this slab, find the point at which the surface pressure is minimum
+def find_devil(nc,indextime):
+    import numpy as np
+    from scipy import ndimage
+    from mymath import array2image,image2array
+
+    varinfile = nc.variables.keys()
+    if "PSFC" not in varinfile: errormess("You need PSFC in your file to find dust devils")
+    else: psfc_full=getfield(nc,'PSFC')
+    psfc,error=reducefield( psfc_full, d4=indextime)
+    psfcim=array2image(1000.*(psfc-psfc.mean()))
+    sx = ndimage.sobel(psfcim, axis=0, mode='constant') ; sy = ndimage.sobel(psfcim, axis=1, mode='constant')
+    sob = np.hypot(sx, sy)
+    zemax=np.max(sob)
+    goodvalues = sob[sob >= zemax/2]
+    ix = np.in1d(sob.ravel(), goodvalues).reshape(sob.shape)
+    idxs,idys=np.where(ix)
+    maxvalue = sob[sob == zemax]
+    ixmax = np.in1d(sob.ravel(), maxvalue[0]).reshape(sob.shape)
+    idxmax,idymax=np.where(ixmax)
+    valok=[]
+    for i in np.arange(len(idxs)):
+        a=np.sqrt((idxmax-idxs[i])**2 + (idymax-idys[i])**2)
+        if 0 < a <= 5.*np.sqrt(2.): valok.append(goodvalues[i])
+    ix = np.in1d(sob.ravel(), valok).reshape(sob.shape)
+    idxs,idys=np.where(ix)
+    hyperslab=psfc[np.min(idxs):np.max(idxs),np.min(idys):np.max(idys)]
+    idxsub,idysub=np.where(hyperslab==hyperslab.min())
+    idx=idxsub[0]+np.min(idxs) ; idy=idysub[0]+np.min(idys)
+    return np.int(idx),np.int(idy)

trunk/UTIL/PYTHON/myscript.py

@@ -52 +52 @@
     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 lon,lat [None]')
+    parser.add_option('--finddevil',       action='store_true',dest='mdevil',               default=False, help='superimpose a crossmark where the steepest dust devil is [False]')
 
     ### SPECIFIC FOR SLICING [MAPMODE 0]

trunk/UTIL/PYTHON/planetoplot.py

@@ -72 +72 @@
            lstyle=None,\
            cross=None,\
+           markdevil=False,\
            facwind=1.,\
            trycol=False,\
@@ -87 +88 @@
                        zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\
                        getname,localtime,check_localtime,polarinterv,getsindex,define_axis,determineplot,readslices,bidimfind,getlschar,hole_bounds,\
-                       getdimfromvar,select_getfield
+                       getdimfromvar,select_getfield,find_devil
     from mymath import deg,max,min,mean,writeascii,fig2data,fig2img
     import matplotlib as mpl
@@ -274 +275 @@
               lonp,latp = ( lon2d[indices[0,0,0],indices[0,0,1]] , lat2d[indices[0,0,0],indices[0,0,1]] )
           ### we get rid of boundary relaxation zone for plots. important to do that now and not before.
-          if typefile in ['meso'] and mapmode == 1: lon2d = dumpbdy(lon2d,6) ; lat2d = dumpbdy(lat2d,6)  
+          if typefile in ['meso'] and mapmode == 1:
+             if '9999' not in getattr(nc,'START_DATE'): lon2d = dumpbdy(lon2d,6) ; lat2d = dumpbdy(lat2d,6)  
           ### we read the keyword for vertical dimension. we take care for vertical staggering.
           if varname in ['PHTOT','W']:    vertdim='BOTTOM-TOP_PATCH_END_STAG'
@@ -589 +591 @@
                         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 ['meso'] and mapmode == 1:   what_I_plot_frame = dumpbdy(what_I_plot_frame,6,condition=True)
+                    if typefile in ['meso'] and mapmode == 1: 
+                       if '9999' not in getattr(nc,'START_DATE'): 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)
 
@@ -693 +696 @@
                             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)
                             elif mapmode == 0:     pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans)
-                        
+
+                        if (cross is not None or markdevil) and mapmode == 1:
+                            if cross is not None: idx,idy=m(cross[0][0],cross[0][1])
+                            elif markdevil:
+                                idx,idy=find_devil(nc,indextime)
+                                idx,idy=x[idx,idy],y[idx,idy]
+                            mpl.pyplot.plot([idx],[idy],'+k',mew=0.5,ms=18)
 
                         if colorb not in ['nobar','onebar']:        
@@ -713 +719 @@
                         if winds:
                             if typefile in ['meso']:
-                                [vecx_frame,vecy_frame] = [dumpbdy(vecx_frame,6,stag=uchar,condition=True), dumpbdy(vecy_frame,6,stag=vchar,condition=True)]
+                                if '9999' not in getattr(nc,'START_DATE') : [vecx_frame,vecy_frame] = [dumpbdy(vecx_frame,6,stag=uchar,condition=True), dumpbdy(vecy_frame,6,stag=vchar,condition=True)]
                                 key = True
                                 if fvar in ['UV','uv','uvmet']: key = False

trunk/UTIL/PYTHON/pp.py

@@ -176 +176 @@
                 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),\
-                lstyle=separatenames(opt.linestyle),cross=readslices(opt.mark),facwind=opt.facwind,\
+                lstyle=separatenames(opt.linestyle),cross=readslices(opt.mark),markdevil=opt.mdevil,facwind=opt.facwind,\
                 trycol=opt.trycol,streamflag=opt.stream,analysis=opt.analysis)
         print 'DONE: '+name