Changeset 792 for trunk/UTIL/PYTHON/myplot.py

Timestamp:

Sep 21, 2012, 3:03:27 PM (12 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

File:

• trunk/UTIL/PYTHON/myplot.py (modified) (1 diff)

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)