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

Last change on this file since 1192 was 1070, checked in by aslmd, 12 years ago

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