source: trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/find_devils.py @ 296

#! /usr/bin/env python

def detsize( xx, res=1, thres=3, loga=False ):
    import numpy as np
    import math
    size = []
    sizecalc = 1
    diff = np.asarray( np.roll(xx,-1) - xx )
    for i in diff:
        if abs(i) > 1:
            if sizecalc >= thres: 
                if loga: addthis = math.log(sizecalc*res)
                else:    addthis = sizecalc*res
                size.append(addthis)
            sizecalc = 1
        else:
            sizecalc += 1
    return size

def getsize(filename):

    import numpy as np
    from scipy.ndimage.measurements import minimum_position
    from netCDF4 import Dataset
    
    ### LOAD NETCDF DATA
    filename = "psfc.nc"
    nc = Dataset(filename) 
    psfc = nc.variables["PSFC"]
    
    ### LOOP on TIME
    ### NB: a same event could be counted several times...
    shape = np.array(psfc).shape
    allsizesx = []
    allsizesy = []
    for i in range(0,shape[0]):
    #for i in range(0,2):
    
        print i, ' on ', shape[0]
        psfc2d = np.array ( psfc [ i, : , : ] )
    
        ############### CRITERION
        ave = np.mean(psfc2d,dtype=np.float64) ## dtype otherwise inaccuracy
        where = np.where(psfc2d - ave < -0.3)  ## comme le papier Phoenix
    
        std = np.std(psfc2d,dtype=np.float64) ## dtype otherwise inaccuracy
        fac = 2.5  ## not too low, otherwise vortices are not caught. 4 good choice.
        lim = ave - fac*std
        print lim, ave, std
        where = np.where(psfc2d < lim)
        ############### END CRITERION
    
        lab = np.zeros(np.array(psfc2d).shape) ## points to be treated by the minimum_position routine
        lab[where] = 1.  ## do not treat points close to 'mean' (background) pressure
    
        xx = []
        yy = []
        while 1 in lab:
            p = minimum_position(psfc2d,labels=lab)
            lab[p] = 0 ## once a minimum has been found in a grid point, do not search here again.
            if p[0] not in xx: xx.append(p[0]) ## if x coordinate not yet in the list add it
            if p[1] not in yy: yy.append(p[1]) ## if y coordinate not yet in the list add it
        xx.sort()
        yy.sort()
        ### now xx and yy are sorted arrays containing grid points with pressure minimum
       
        ######## DETERMINE SIZE OF STRUCTURES
        ######## this is rather brute-force...
        sizex = detsize( xx, res = 10, loga=False, thres=2 )
        sizey = detsize( yy, res = 10, loga=False, thres=2 )
        ###
        allsizesx = np.append(allsizesx,sizex)
        allsizesy = np.append(allsizesy,sizey)
        ########
    
    allsizesx.sort()
    allsizesy.sort()
    
    return allsizesx, allsizesy

def writeascii ( tab, filename ):
    mydata = tab
    myfile = open(filename, 'w')
    for line in mydata:
        myfile.write(str(line) + '\n')
    myfile.close()
    return

import matplotlib.pyplot as plt
import pickle
import numpy as np
import matplotlib.mlab as mlab

save = True
#save = False
if save:
    allsizesx, allsizesy = getsize("psfc.nc")
    writeascii(allsizesx,'allsizex.txt')
    writeascii(allsizesy,'allsizey.txt')
    
    myfile = open('allsizex.bin', 'wb')
    pickle.dump(allsizesx, myfile)
    myfile.close()
    
    myfile = open('allsizey.bin', 'wb')
    pickle.dump(allsizesy, myfile)
    myfile.close()


myfile = open('allsizex.bin', 'r')
allsizesx = pickle.load(myfile)
myfile = open('allsizey.bin', 'r')
allsizesy = pickle.load(myfile)

plothist = np.append(allsizesx,allsizesy)
plothist.sort()

nbins=100
zebins = [2.0]
for i in range(0,nbins):  zebins.append(zebins[i]*np.sqrt(2))
zebins = np.array(zebins)
print zebins

zebins = zebins [ zebins > 10. ] 
zebins = zebins [ zebins < 1000. ]
print zebins

#### HISTOGRAM
plt.hist(    plothist,\
             log=True,\
             bins=zebins,\
             #cumulative=-1,\
        )
plt.xscale('log')

plt.show()