source: trunk/UTIL/PYTHON/powerlaw/finddd.py @ 1062

Last change on this file since 1062 was 1053, checked in by aslmd, 12 years ago

UTIL PYTHON. tools to study dust devils in LES.
  • Property svn:executable set to *
File size: 15.6 KB
Line 
1#! /usr/bin/env python
2from ppclass import pp
3from ppplot import plot2d,save,writeascii
4import numpy as np
5from scipy.ndimage.measurements import minimum_position
6import matplotlib.pyplot as mpl
7import ppcompute
8
9###############################################################################
10# FIND DUST DEVILS
11# filefile --> file
12# dx --> horizontal resolution
13# dt_out --> frequency of outputs (s)
14# lt_start --> starting local time
15# halolim --> limit for halos
16# method --> method used to detect vortices
17# plotplot --> plot or not
18# save --> save results or not
19###############################################################################
20def finddd(filefile,\
21           timelist,\
22           dx=50.,\
23           dt_out=50.,\
24           lt_start=8.,\
25           halolim=30.,\
26           method=1,\
27           plotplot=False,\
28           save=True):
29
30    ###############################################################################
31    ########################## FOR METHOD 1 FOR METHOD 1 ##########################
32    ## FACLIST : how many sigmas below mean we start to consider this could be a vortex
33    ## ... < 3.0 dubious low-intensity minima are caught
34    ## ... 3 ~ 3.1 is a bit too low but helps? because size of pressure drop could be an underestimate of actual dust devil
35    ## ... 3.2 ~ 3.3 is probably right, this is the one we choose for exploration [3.25]
36    ## ... NB: values 3.2 to 3.6 yields the same number of devils but measured sizes could vary a bit
37    ## ... > 3.6 makes strong minima disappear... especially >3.8
38    ## ... EVENTUALLY 3.5-4 is better for problematic cases
39    ###############################################################################
40    faclist = [3.75]
41
42    ## (see below) NEIGHBOR_FAC is the multiple of std used to evaluate size
43    ## --> 1: limit not discriminative enough. plus does not separate neighbouring vortices.
44    ##        ... but interesting: gives an exponential law (because vortices are artificially merged?)
45    ## --> 2.7: very good for method 1. corresponds usually to ~0.3
46    ## --> 2: so-so. do not know what to think.
47    neighbor_fac = 2.7
48
49    ###############################################################################
50    if save:
51        myfile1 = open(filefile+'m'+str(method)+'_'+'1.txt', 'w')
52        myfile2 = open(filefile+'m'+str(method)+'_'+'2.txt', 'w')
53    ###############################################################################
54   
55    ## mean and std calculations
56    ## -- std is used in both methods for limits
57    ## -- mean is only used in method 1
58    print "calculate mean and std, please wait."
59    ## -- get time series of 2D surface pressure
60    psfc = pp(file=filefile,var="PSFC").getf()
61    ## -- calculate mean and standard deviation
62    ## -- ... calculating std at all time is not right!
63    ## -- ... for mean value though, similar results with both methods
64    mean = np.mean(psfc,dtype=np.float64)
65    std = np.std(psfc,dtype=np.float64)
66    ## some information about inferred limits
67    print "**************************************************************"
68    print "MEAN",mean
69    print "STD",std
70    print "LIMIT FOR PRESSURE MINIMUM:",-np.array(faclist)*std
71    print "LIMIT FOR EVALUATING SIZE OF A GIVEN LOCAL MINIMUM",-neighbor_fac*std
72    print "**************************************************************"
73   
74    ## LOOP ON TIME
75    for time in timelist:
76   
77     ## get 2D surface pressure at a given time
78     ## (this is actually so quick we don't use psfc above)
79     psfc2d = pp(file=filefile,var="PSFC",t=time).getf()
80     #ustm = pp(file=filefile,var="USTM",t=time).getf()
81
82     ## MAIN ANALYSIS. LOOP ON FAC.
83     for fac in faclist:
84     #fac = 3.75
85     #for method in [2,1]:
86 
87      ## initialize arrays
88      tabij = [] ; tabsize = [] ; tabdrop = []
89      tabijcenter = [] ; tabijvortex = [] ; tabijnotconv = [] ; tabdim = []
90   
91      ################ FIND RELEVANT POINTS TO BE ANALYZED
92      ## lab is 1 for points to be treated by minimum_position routine
93      ## we set elements at 1 where pressure is under mean-fac*std
94      ## otherwise we set to 0 because this means we are close enough to mean pressure (background)
95      lab = np.zeros(psfc2d.shape)
96      if method == 1:
97          # method 1: standard deviation
98          lab[np.where(psfc2d < mean-fac*std)] = 1
99      else:
100          # method 2: polynomial fit
101          # ... tried smooth but too difficult, not accurate and too expensive
102          deg = 5 #plutot bien (deg 10 ~pareil) mais loupe les gros (~conv cell)
103          #deg = 2 #pas mal mais false positive (pareil 3-4 mm si un peu mieux)
104          #        #(OK now with fixing the false positive bug)
105          nx = psfc2d.shape[1] ; ny = psfc2d.shape[0]
106          xxx = np.array(range(nx)) ; yyy = np.array(range(ny))
107          anopsfc2d = psfc2d*0. ; polypsfc2d = psfc2d*0.
108          for iii in range(0,nx,1):
109             poly = np.poly1d(np.polyfit(yyy,psfc2d[iii,:],deg))
110             polypsfc2d[iii,:] = poly(yyy)
111          for jjj in range(0,ny,1):
112             poly = np.poly1d(np.polyfit(xxx,psfc2d[:,jjj],deg))
113             polypsfc2d[:,jjj] = 0.5*polypsfc2d[:,jjj] + 0.5*poly(xxx)
114          ## smooth a little to avoid 'crosses' (plus, this removes PBC problems)
115          polypsfc2d = ppcompute.smooth2diter(polypsfc2d,n=deg)
116          # compute anomaly and find points to be explored
117          anopsfc2d = psfc2d - polypsfc2d
118          limlim = fac*std ## same as method 1
119          lab[np.where(anopsfc2d < -limlim)] = 1
120   
121      ## while there are still points to be analyzed...
122      while 1 in lab:
123        ## ... get the point with the minimum field values
124        if method == 1:
125            ij = minimum_position(psfc2d,labels=lab)
126        else:
127            ij = minimum_position(anopsfc2d,labels=lab)
128        ## ... store the indexes of the point in tabij
129        tabij.append(ij)
130        ## ... remove the point from labels to be further explored by minimum_position
131        lab[ij] = 0
132   
133      ################ GET SIZES BASED ON THOSE FOUND POINTS
134      ## reslab is the same as lab
135      ## except for scanning purpose we keep the information
136      ## about how went the detection
137      ## --> and we set a lower fac
138      ## --> above a high fac is good to catch only strong vortices
139      ## --> but here a casual fac=3 is better to get accurate sizes
140      ## --> or even lower as shown by plotting reslab
141      reslab = np.zeros(psfc2d.shape)
142      if method == 1:
143          reslab[np.where(psfc2d < mean-neighbor_fac*std)] = 1
144      else:
145          reslab[np.where(anopsfc2d < -neighbor_fac*std)] = 1
146     
147      ## initialize halomax and while loop
148      ## HALOMAX : maximum halo defined around a minima to evaluate the size
149      ## ... halomax must be large enough to encompass a vortex
150      ## ... but not too large otherwise neighboring vortex are caught
151      halomax = 3 ; notconv = 9999 ; yorgl = 9999
152   
153      ## WHILE LOOP on HALOMAX exploration
154      while ( notconv > 0 and halomax < halolim and yorgl != 0 ):
155       # now browse through all points caught in previous loop
156       for ij in tabij:
157        ## ... OK. take each indexes found before with minimum_position
158        i,j = ij[0],ij[1]
159        ## EITHER
160        ## ... if reslab is already 0, we do not have to do anything
161        ## ... because this means point is already part of detected vortex
162        ## OR
163        ## ... if the ij couple is already in a previously detected vortex
164        ## ... we don't actually need to consider it again
165        ## ... this is necessary otherwise (sometimes a lot) of false positives
166        if reslab[i,j] <= 0 or ij in tabijvortex:
167          pass
168        else:
169          ## ... then define a growing halo around this point
170          ## ... we make the halo grow until convergence of
171          ## ... the number of under-limit points (i.e. with lab=1) within halo
172          ## ... which means the whole vortex is encompassed
173          ## ... we start with a halo of 1 point around minimum point
174          ## ... and we end when we reached halomax which likely means dubious case
175          halo = 1 ; nmesh = -9999. ; prevmesh = 9999. ; notconverged = False
176          while nmesh != prevmesh and halo <= halomax:
177              ## ... save the number of vortex points calculated at previous iteration
178              prevmesh = nmesh
179              ## ... define a halo around the minimum point
180              minx,maxx,miny,maxy = i-halo,i+halo+1,j-halo,j+halo+1
181              ## ... treat the boundary case (TBD: periodic boundary conditions)
182              if minx < 0: minx = 0
183              if miny < 0: miny = 0
184              if maxx > psfc2d.shape[0]: maxx = psfc2d.shape[0]
185              if maxy > psfc2d.shape[1]: maxy = psfc2d.shape[1]
186              ## ... define the patch, made of the halo of points
187              patch = reslab[minx:maxx,miny:maxy]       
188              ## ... count how many 1 are inside this patch
189              ## ... these are points for which value is >0
190              ## ... because close to mean is 0 and already caught is <0
191              ## ... and not converged are >1 so still to be scanned
192              nmesh = len(np.where(patch >= 1)[0])
193              ## ... in case widening halo adds 1-2 points only
194              ## ... we consider vortex is encompassed
195              ## ... this helps to separate neighboring vortices
196              ## ... and only yields marginal inaccuracies
197              if halo > 3 and abs(nmesh-prevmesh)<=2: prevmesh = nmesh
198              ## ... if with halo=1 we caught a 1-point vortex, end here because spurious
199              if halo == 1 and nmesh == 1: prevmesh = nmesh
200              ## ... with the last halo=halomax test, we can detect not-converged cases
201              if halo == halomax and nmesh != prevmesh: notconverged = True
202              ## ... increment halo size before the next iteration
203              halo = halo + 1
204          ## ... now we got the grid points encompassed by the vortex in nmesh
205          ## ... then to continue scanning we store results in reslab
206          if notconverged:
207            ## multiply reslab by 2.
208            ## --> if it is >1, point could be part of another vortex found later
209            ## --> if it is <0, point is already within a caught vortex, so reslab should remain <0
210            ## --> if it is =0, this should remain 0 because close to average
211            reslab[i,j] = reslab[i,j]*2
212            tabijnotconv.append(ij)
213          else:
214            ## OK. this is most likely an actual vortex. we get the drop.
215            ## we multiply by mesh area, then square to get approx. size of vortex
216            if method == 1:
217                drop = -psfc2d[i,j]+mean
218            else:
219                drop = -anopsfc2d[i,j]
220            size = int(np.sqrt(nmesh*dx*dx))
221            ## if vortex is too small (i.e. too close to mesh grid resolution) we don't store information
222            ## ... we just remove the patch from labels to be further explored
223            facdx = 2.
224            if size < facdx*dx:
225                reslab[minx:maxx,miny:maxy] = 0
226            ## otherwise it is a VORTEX! we store info in arrays
227            else:
228                ## we put the vortex points into tabijvortex so that
229                ## a (i,j) couple into a vortex is not considered in further explorations
230                ## -- also we evaluate the x-size and y-size of the vortex (max distance to center)
231                ## -- this can be useful to detect not-so-round fake vortices (convective gusts?)
232                maxw = -9999 ; maxh = -9999 ; maxu = -9999
233                for iii in range(minx,maxx):
234                 for jjj in range(miny,maxy):
235                  if reslab[iii,jjj] != 0:
236                     tabijvortex.append((iii,jjj))
237                     width = np.abs(iii-ij[0])
238                     if width > maxw: maxw = width
239                     height = np.abs(jjj-ij[1])
240                     if height > maxh: maxh = height
241                     #ustar = ustm[iii,jjj]
242                     #if ustar > maxu: maxu = ustar
243                #print size,drop,ustar
244                ## store info in dedicated arrays
245                tabdim.append((maxw*dx,maxh*dx))
246                tabsize.append(size)
247                tabdrop.append(drop)
248                tabijcenter.append(ij)
249                #print "... VORTEX!!!! size %.0f drop %.1f coord %.0f %.0f" % (size,drop,i,j)
250                ## we remove the patch from labels to be further explored
251                ## we multiply reslab by -1 to plot detection maps
252                reslab[minx:maxx,miny:maxy] = patch*-1
253     
254       ## count how many points are not converged and left to be analyzed
255       notconv = len(np.where(reslab > 1)[0])
256       yorgl = len(np.where(reslab == 1)[0])
257   
258       ## increment halomax
259       ## to speed-up the increment is slightly increasing with considered halomax
260       halomax = halomax + halomax / 2
261   
262      ## just for simpler plots.
263      reslab[reslab > 2] = 4
264      reslab[reslab < -2] = -4
265
266      ## give some info to the user
267      if len(tabsize) > 0:
268        nvortex = len(tabsize)
269        maxsize = np.max(tabsize)
270        maxdrop = np.max(tabdrop)
271      else:
272        nvortex = 0
273        maxsize = 0
274        maxdrop = 0.
275      notconv = len(np.where(reslab > 1)[0])
276      print "t=%3.0f / n=%2.0f / s_max=%4.0f / d_max=%4.1f / halo_out=%3.0f / notconvp=%3.1f" \
277            % (time,nvortex,maxsize,maxdrop,halomax,100.*notconv/float(reslab.size))       
278   
279      ## save results in a text file
280      if save:
281          # convert t in local time
282          ttt = lt_start + time*dt_out/3700.     
283          # write files
284          myfile2.write( "%5.2f ; %5.0f ; %5.0f ; %7.2f\n" % (ttt,nvortex,maxsize,maxdrop) )
285          for iii in range(len(tabsize)):
286              myfile1.write( "%5.2f ; %5.0f ; %7.2f ; %5.0f ; %5.0f\n" \
287              % (ttt,tabsize[iii],tabdrop[iii],tabdim[iii][0],tabdim[iii][1]) )
288   
289      #### PLOT PLOT PLOT PLOT
290      if nvortex>0 and plotplot:
291       mpl.figure(figsize=(12,8))
292       myplot = plot2d()
293       myplot.absc = np.array(range(psfc2d.shape[1]))*dx/1000.
294       myplot.ordi = np.array(range(psfc2d.shape[0]))*dx/1000.
295       myplot.title = str(nvortex)+" vortices found (indicated diameter / pressure drop)"
296       myplot.xlabel = "x distance (km)"
297       myplot.ylabel = "y distance (km)"
298       if method > 0:
299       #if method == 1:
300           #myplot.field = ustm
301           myplot.field = psfc2d
302           #myplot.field = polypsfc2d
303           myplot.vmin = -2.*std
304           myplot.vmax = +2.*std
305           myplot.vmin = mean - 6.*std
306           myplot.vmax = mean + 6.*std
307       else:
308           myplot.field = anopsfc2d
309           myplot.vmin = -1.5
310           myplot.vmax = 0.5
311       myplot.fmt = "%.1f"
312       myplot.div = 20
313       myplot.colorb = "spectral"
314       myplot.make()
315     
316       ### ANNOTATIONS
317       for iii in range(len(tabsize)):
318        ij = tabijcenter[iii]
319        coord1 = ij[1]*dx/1000.
320        coord2 = ij[0]*dx/1000.
321        txt = "%.0f/%.1f" % (tabsize[iii],tabdrop[iii])
322        mpl.annotate(txt,xy=(coord1,coord2),
323             xytext=(-10,-30),textcoords='offset points',ha='center',va='bottom',\
324             bbox=dict(boxstyle='round,pad=0.2', fc='yellow', alpha=0.3),\
325             arrowprops=dict(arrowstyle='->', connectionstyle='arc3,rad=-0.3',color='red'),\
326             size='small')
327   
328       ###show detection
329       #lev = [-4,-2,-1,0,1,2,4]
330       #lev = [-1,0]
331       #mpl.contourf(myplot.absc,myplot.ordi,reslab,alpha=0.9,cmap=mpl.cm.get_cmap("binary_r"),levels=lev)
332   
333       ### SHOW OR SAVE IN FILE
334       mpl.show()
335       #save(mode="png",filename="detectm"+"_"+str(time)+"_"+str(method),folder="detect/",includedate=False)
336   
337    ## close data files
338    if save:
339        myfile1.close()
340        myfile2.close()
Note: See TracBrowser for help on using the repository browser.