source: trunk/UTIL/PYTHON/powerlaw/find_devils.py @ 963

Last change on this file since 963 was 943, checked in by aslmd, 12 years ago

clean and organized UTIL/PYTHON folder
  • Property svn:executable set to *
File size: 10.7 KB
Line 
1#! /usr/bin/env python
2
3def detsize( xx, res=1, thres=3, loga=False ):
4    import numpy as np
5    import math
6    size = []
7    sizecalc = 1
8    diff = np.asarray( np.roll(xx,-1) - xx )
9    for i in diff:
10        if abs(i) > 1:
11            if sizecalc >= thres: 
12                if loga: addthis = math.log(sizecalc*res)
13                else:    addthis = sizecalc*res
14                size.append(addthis)
15            sizecalc = 1
16        else:
17            sizecalc += 1
18    return size
19
20def getsize(filename):
21
22    import numpy as np
23    from scipy.ndimage.measurements import minimum_position
24    from netCDF4 import Dataset
25    import matplotlib.pyplot as plt
26    import myplot as myp
27   
28    ### LOAD NETCDF DATA
29    nc = Dataset(filename) 
30    psfc = nc.variables["PSFC"]
31   
32    ### LOOP on TIME
33    ### NB: a same event could be counted several times...
34    shape = np.array(psfc).shape
35    allsizesx = []
36    allsizesy = []
37    depression = []
38    stride = 1 #5
39    stride = 20
40    #stride = 50
41    #stride = 100
42    start = 0
43    start = stride
44    for i in range(start,shape[0],stride):
45
46        psfc2d = np.array ( psfc [ i, : , : ] )
47   
48        ############### CRITERION
49        ave = np.mean(psfc2d,dtype=np.float64)  ## dtype otherwise inaccuracy
50
51        #limdp = -0.2 ## on en loupe pas mal
52        #where = np.where(psfc2d - ave < limdp)  ## comme le papier Phoenix
53   
54        std = np.std(psfc2d,dtype=np.float64)   ## dtype otherwise inaccuracy
55        fac = 4.   ## how many sigmas. not too low, otherwise vortices are not caught. 4 good choice.
56                   ## 2.5 clearly too low, 3.5 not too bad, 4 probably good
57        fac = 3.5
58        fac = 3.2
59        ##fac = 2.5
60        #fac = 3. ## final choice
61        #fac = 2.5
62        lim = ave - fac*std
63        where = np.where(psfc2d < lim)
64        ############### END CRITERION
65
66        depression = np.append(depression,np.ravel(psfc2d[where])-ave)
67   
68        lab = np.zeros(np.array(psfc2d).shape) ## points to be treated by the minimum_position routine
69        lab[where] = 1.  ## do not treat points close to 'mean' (background) pressure
70
71        draw = False
72        #draw = True
73        if draw:
74        ##################################################################################
75            vmin = -0.3
76            vmax =  0.0
77            ndiv = 3
78            palette = plt.get_cmap(name="YlGnBu") 
79            what_I_plot = psfc2d-ave
80            zevmin, zevmax = myp.calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)
81            what_I_plot = myp.bounds(what_I_plot,zevmin,zevmax)
82            zelevels = np.linspace(zevmin,zevmax)
83            fig = plt.figure(figsize=(16,8)) 
84            subv,subh = myp.definesubplot(2,fig)
85            plt.subplot(subv,subh,1) 
86            plt.contourf(what_I_plot,zelevels,cmap=palette) 
87            plt.colorbar(fraction=0.05,pad=0.03,format="%.1f",\
88                           ticks=np.linspace(zevmin,zevmax,ndiv+1),\
89                           extend='both',spacing='proportional')
90            plt.subplot(subv,subh,2)
91            palette = plt.get_cmap(name="hot") #"binary")
92            #lab[np.where(lab > 0.99)] = np.NaN
93            #plt.pcolor(lab,cmap=palette)
94            plt.contourf(lab,2,cmap=palette)
95            plt.show() 
96        ##################################################################################
97 
98        xx = []
99        yy = []
100        while 1 in lab:
101            p = minimum_position(psfc2d,labels=lab)
102            lab[p] = 0 ## once a minimum has been found in a grid point, do not search here again.
103            if p[0] not in xx: xx.append(p[0]) ## if x coordinate not yet in the list add it
104            if p[1] not in yy: yy.append(p[1]) ## if y coordinate not yet in the list add it
105        xx.sort()
106        yy.sort()
107        ### now xx and yy are sorted arrays containing grid points with pressure minimum
108       
109        ######## DETERMINE SIZE OF STRUCTURES
110        ######## this is rather brute-force...
111        sizex = detsize( xx, res = 10, loga=False, thres=2 )
112        sizey = detsize( yy, res = 10, loga=False, thres=2 )
113        #sizex = detsize( xx, res = 10, loga=False, thres=3 )
114        #sizey = detsize( yy, res = 10, loga=False, thres=3 )
115        sizex = detsize( xx, res = 15, loga=False, thres=2 )
116        sizey = detsize( yy, res = 15, loga=False, thres=2 )
117        ###
118        print sizex, sizey
119        #if ( mym.max(sizex) > mym.max(sizey) ): sizey = sizex  ### un peu limite dans certains cas
120        if (len(sizex) > len(sizey))     : sizey = sizex        ### plus fidele mais petit souci lorsque PBC
121        elif (len(sizex) == len(sizey))  : 
122          if ( mym.max(sizex) > mym.max(sizey) ): sizey = sizex
123          else                                  : sizex = sizey
124        else                            : sizex = sizey
125        allsizesx = np.append(allsizesx,sizex)
126        allsizesy = np.append(allsizesy,sizey)
127        print i, ' on ', shape[0], ' caught ', len(sizex), ' vortices ', sizex
128        ########
129 
130    allsizesx.sort()
131    allsizesy.sort()
132   
133    return allsizesx, allsizesy, depression
134
135#########################################################################
136#########################################################################
137
138import matplotlib.pyplot as plt
139import pickle
140import numpy as np
141import matplotlib.mlab as mlab
142import mymath as mym
143import myplot as myp
144
145import plfit
146import plplot
147import randht
148import plpva
149
150save = True
151#save = False
152pression = False
153#pression = True
154
155filename = "/home/aymeric/Big_Data/psfc_f18.nc"
156
157if save:
158    allsizesx, allsizesy, depression = getsize(filename)
159    ### sauvegarde texte pour inspection
160    mym.writeascii(allsizesx,'allsizex.txt')
161    mym.writeascii(allsizesy,'allsizey.txt')
162    mym.writeascii(depression,'alldepression.txt')
163    ### sauvegarde binaire pour utilisation python
164    myfile = open('allsizex.bin', 'wb')
165    pickle.dump(allsizesx, myfile)
166    myfile.close()
167    myfile = open('allsizey.bin', 'wb')
168    pickle.dump(allsizesy, myfile)
169    myfile.close()
170    myfile = open('alldepression.bin', 'wb')
171    pickle.dump(depression, myfile)
172    myfile.close()
173
174### load files
175myfile = open('allsizex.bin', 'r')
176allsizesx = pickle.load(myfile)
177myfile = open('allsizey.bin', 'r')
178allsizesy = pickle.load(myfile)
179myfile = open('alldepression.bin', 'r')
180depression = pickle.load(myfile)
181depression = np.array(abs(depression))#*1000.
182
183### sizes
184#plothist = np.append(allsizesx,allsizesy)
185plothist = allsizesx
186if pression: plothist = depression
187plothist.sort()
188print 'mean ', np.mean(plothist,dtype=np.float64) 
189print 'std ', np.std(plothist,dtype=np.float64)
190print 'max ', np.max(plothist)
191print 'min ', np.min(plothist)
192print 'len ', len(plothist)
193
194
195### MAKE BINS
196nbins = 100
197zebins = [2.0]
198#nbins = 8
199#zebins = [19.0]
200#nbins = 15
201#zebins = [11.] ##12 non mais donne un peu la meme chose
202#zebins = [20.] ##20 non car trop pres du premier
203nbins = 100
204zebins = [2./np.sqrt(2.)]  ## ne pas tomber sur une dizaine ronde
205nbins = 200
206
207if pression: zebins = [0.3]
208
209for i in range(0,nbins):  zebins.append(zebins[i]*np.sqrt(2))
210zebins = np.array(zebins)
211#### select reasonable bins for DD
212if not pression:
213    zebins = zebins [ zebins > 15. ] 
214    #zebins = zebins [ zebins > 20. ]
215    zebins = zebins [ zebins > 25. ]
216    zebins = zebins [ zebins < 1000. ]
217else:
218    zebins = zebins [ zebins < 10. ]
219print 'corrected bins ',zebins
220
221#### HISTOGRAM
222plt.figure(1)
223plt.hist(    plothist,\
224             log=True,\
225             bins=zebins,\
226#             cumulative=-1,\
227             normed=True,\
228        )
229plt.xscale('log')
230if pression:    plt.xlabel('Pressure (Pa)')
231else:           plt.xlabel('Diameter (m)')
232plt.ylabel('Population (normalized)')
233if pression: prefix="p"
234else:        prefix=""
235myp.makeplotres(prefix+"histogram",res=200,disp=False)
236plt.close(1)
237
238### COMPARED HISTOGRAMS
239### --- FIT WITH POWER LAW
240if pression: [alpha, xmin, L] = plfit.plfit(plothist,'xmin',0.3)
241else:        [alpha, xmin, L] = plfit.plfit(plothist,'limit',20.)
242print alpha,xmin
243
244#a = plpva.plpva(plothist,0.75,'xmin',0.75)
245#print a
246
247#### DEUXIEME ROUTINE
248####IL FAUT UTILISER LE DISCRET POUR LA TAILLE !!!
249#if pression:   myplfit = plfit.plfit(plothist,verbose=True,xmin=0.75)
250#else:          myplfit = plfit.plfit(plothist,verbose=True,xmin=20.)
251#myplfit.plotppf()
252#plt.show()
253#exit()
254
255
256#plt.figure(1)
257#h = plplot.plplot(plothist,xmin,alpha)
258#myp.makeplotres(prefix+"fit",res=200,disp=False)
259plt.figure(2)
260### --- POWER LAW (factor does not really matter)
261power = (xmin/2.2)*np.array(randht.randht(10000,'powerlaw',alpha))
262#power = (xmin/2.2)*np.array(randht.randht(10000,'cutoff',alpha,10.)) ##marche pas si trop grand
263print 'mean ', np.mean(power,dtype=np.float64)
264### --- EXPONENTIAL LAW
265expo = randht.randht(10000,'exponential',1./(np.mean(power,dtype=np.float64)*1.00))
266print 'mean ', np.mean(expo,dtype=np.float64)
267### --- PLOT
268plt.hist(    [plothist,power,expo],\
269             label=['LES vortices','Power law '+'{:.1f}'.format(alpha),'Exponential law'],\
270             log=True,\
271             bins=zebins,\
272#            cumulative=-1,\
273             normed=True,\
274        )
275plt.legend()
276plt.xscale('log')
277if pression:    plt.xlabel('Pressure (Pa)')
278else:           plt.xlabel('Diameter (m)')
279plt.ylabel('Population (normalized)')
280myp.makeplotres(prefix+"comparison",res=200,disp=False)
281plt.close(2)
282
283########################
284########################
285zebins = [30.,42.,60.,84.,120.,170.,240.,340.]
286plothist = []
287plothist = np.append(plothist,30 *np.ones(306))
288plothist = np.append(plothist,42 *np.ones(58) )
289plothist = np.append(plothist,60 *np.ones(66) )
290plothist = np.append(plothist,84 *np.ones(41) )
291plothist = np.append(plothist,120*np.ones(19) )
292plothist = np.append(plothist,170*np.ones(9)  )
293plothist = np.append(plothist,240*np.ones(2)  )
294plothist = np.append(plothist,340*np.ones(1)  )
295
296#zebins = [50.,71.,100.,141.,200.,282.,400.]
297#plothist = []
298#plothist = np.append(plothist,50. *np.ones(36))
299#plothist = np.append(plothist,71. *np.ones(18) )
300#plothist = np.append(plothist,100. *np.ones(12) )
301#plothist = np.append(plothist,141. *np.ones(6) )
302#plothist = np.append(plothist,200.*np.ones(4) )
303#plothist = np.append(plothist,282.*np.ones(1)  )
304#plothist = np.append(plothist,400.*np.ones(2)  )
305
306exit()
307
308plt.figure(3)
309[alpha, xmin, L] = plfit.plfit(plothist,'xmin',30)#50.)
310print alpha,xmin
311#a = plpva.plpva(plothist,30,'xmin',30)
312h = plplot.plplot(plothist,xmin,alpha)
313plt.loglog(h[0], h[1], 'k--',linewidth=2)
314plt.hist(    plothist,\
315             log=True,\
316             bins=zebins,\
317#             cumulative=-1,\
318             normed=True,\
319        )
320plt.xscale('log')
321plt.xlabel('Pressure (micro Pa)')
322plt.ylabel('Population (normalized)')
323myp.makeplotres("data",res=200,disp=False)
324
325#plt.figure(4)
326#[alpha, xmin, L] = plfit.plfit(plothist,'xmin',50.) #,'xmin',30.)
327#print alpha,xmin
328#h = plplot.plplot(plothist,xmin,alpha)
329#myp.makeplotres("datafit",res=200,disp=False)
Note: See TracBrowser for help on using the repository browser.