Changeset 1197 for trunk/UTIL


Ignore:
Timestamp:
Mar 4, 2014, 12:00:10 PM (11 years ago)
Author:
aslmd
Message:

PYTHON. updated analysis tools for dust devil analysis.

Location:
trunk/UTIL/PYTHON
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/UTIL/PYTHON/mcd/examples/dustdevil.py

    r1007 r1197  
    99
    1010from ppplot import writeascii
     11
     12
     13from scipy import stats
    1114
    1215
    … …  
    5356## constant_heights
    5457#for heights in [20.,100.,500.,1000.,2000.,5000.,7000.,10000.,20000.]:
    55 for heights in [10.]:
     58#for heights in [10.]:
    5659
    5760## multiple_obsheight
    5861#for heights in [0.1,0.2,0.5,0.75,1.0,1.5,2.0,10.]:
     62for heights in [1.0]:
    5963
    6064 query = mcd()
    … …  
    7276
    7377    ## constant_heights
    74     query.xz = heights
    75 
     78    #query.xz = heights
    7679    ## multiple_obsheight
    77     #query.xz = ddheight[i-1]*heights
     80    query.xz = ddheight[i-1]*heights
    7881
    7982    query.update() ; uu = query.zonwind ; vv = query.merwind
    … …  
    8790       if angle < 0.: angle = angle + 360.
    8891
    89        #query.xzs = max(ddheight[i-1] - 2.*ddheighterror[i-1],10.)
    90        #query.xze = ddheight[i-1] + 2.*ddheighterror[i-1]
    91        query.xzs = max(query.xz - 0.1*query.xz,10.)
    92        query.xze = query.xz + 0.1*query.xz
     92       ## multiple_obsheight
     93       query.xzs = max(ddheight[i-1] - 2.*ddheighterror[i-1],10.)
     94       query.xze = ddheight[i-1] + 2.*ddheighterror[i-1]
     95       ## constant_heights
     96       #query.xzs = max(query.xz - 0.1*query.xz,10.)
     97       #query.xze = query.xz + 0.1*query.xz
    9398
    9499       query.profile()
    … …  
    130135 writeascii(field=mcdwind,absc=ddwind,name=str(heights)+'.txt')
    131136
     137 #ind = np.isnan(ddwind)
     138 #y = ddwind
     139 #y[ind] = -9999.
     140 #xi = mcdwind[y > -9999.]
     141 #y = ddwind[y > -9999.]
     142 #slope, intercept, r_value, p_value, std_err = stats.linregress(xi,y)
     143 #linear = slope*xi+intercept
     144 #mpl.plot(xi,linear)
     145 #ecart = y - linear
     146 #print np.std(ecart)
    132147
    133148 mpl.errorbar(mcdwind, ddwind, yerr=[ddwinderror,ddwinderror], xerr=[mcdwindle,mcdwindhe], fmt='bo')

  • trunk/UTIL/PYTHON/powerlaw/ddhist.py

    r1070 r1197  
    44import numpy as np
    55import scipy.optimize as sciopt
    6 
    7 ### SAVED
    8 #namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm2_1.txt" ; drop = False ; typefit=1
    9 #namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm2_1.txt" ; drop = True ; typefit=2
    10 #namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm1_1.txt" ; drop = False ; typefit=1
    11 #namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm1_1.txt" ; drop = True ; typefit=1
    12 #namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc.LMD_LES_MARS.160564.ncm2_1.txt" ; drop = False ; typefit=1 #bof. tous les 10 est mieux.
    13 #
    14 ###########################################################
    15 ##namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm1_1.txt"
    16 ##namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm2_1.txt"
    17 ##namefile = "/home/aymeric/Big_Data/LES_dd/psfc.LMD_LES_MARS.160564.ncm1_1.txt"
    18 ##namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc.LMD_LES_MARS.160564.ncm2_1.txt"
    19 #namefile = "/home/aymeric/Big_Data/LES_dd/press_ustm_exomars.nc1.txt"
    20 #namefile = "/home/aymeric/Big_Data/LES_dd/sav/press_ustm_exomars.ncm2_1.txt"
    21 ##namefile = "/home/aymeric/Big_Data/LES_dd/psfc_oldinsight100m.ncm1_1.txt"
    22 
    23 
    24 
    25 
    26 case = "188324p"
    27 case = "191798"
    28 case = "160564p"
    29 case = "156487"
    30 case = "2007p"
    31 case = "13526p"
    32 case = "172097"
    33 namefile = "/planeto/aslmd/LESdata/"+case+".ncm1_1.txt"
    34 
    35 
    36 
    37 drop = False ; typefit=1
    38 #drop = True ; typefit=1
    39 #drop = True ; typefit=2
    40 ##########################################################
    41 ## define bins (we expect fit do not depend too much on this -- to be checked)
    42 nbins = 7 ; www = 3.
    43 #nbins = 10 ; www = 2.5
    44 #nbins = 15 ; www = 2.
    45 ##nbins = 20 ; www = 1.75
    46 ##nbins = 30 ; www = 1.5
    47 ##nbins = 50 ; www = 1.2
    48 ##########################################################
    49 limrest = 4. # restrict limit (multiple of dx)
    50 #limrest = 2.
    51 #limrest = 3.
    52 #limrest = 0.
    53 ##########################################################
    546
    557# -- functions for fitting
    … …  
    6214def fitfunc2(x,a,b): return a*np.exp(-b*x)
    6315
    64 # load data
    65 data = np.loadtxt(namefile,delimiter=";")
    66 t = data[:,0] ; s = data[:,1] ; d = data[:,2]
    67 i = data[:,3] ; j = data[:,4] # int?
     16################################################################################
     17### HISTODD
     18### namefile: text file to be analyzed
     19### drop: look at size (False) or pressure drop (True)
     20### typefit: use fitfunc above 1, 2, or 3
     21### nbins: bins (we expect fit do not depend too much on this -- to be checked)
     22### limrest: restrict limit (multiple of dx) -- greater equal
     23### addtitle: add a name for the examined case
     24################################################################################
     25def histodd(namefile,drop=False,typefit=1,nbins=12,limrest=4,addtitle=""):
    6826
    69 # a way to guess resolution
    70 dx = np.min(s)/2. ; print "resolution: ", dx
     27    # width adapted to number of bins
     28    widthbin = {7:3.,10:2.5,12:2.2,15:2.0,20:1.75,30:1.5,50:1.2}
     29    www = widthbin[nbins]
     30   
     31    # load data
     32    data = np.loadtxt(namefile,delimiter=";")
     33    t = data[:,0] ; s = data[:,1] ; d = data[:,2]
     34    i = data[:,3] ; j = data[:,4] # int?
     35   
     36    # a way to guess resolution
     37    dx = np.min(s)/2.
     38   
     39    # choose a variable
     40    if drop: var = d
     41    else: var = s
     42   
     43    # restrictions
     44    restrict = (s > 0) # initialization (True everywhere)
     45    restrict = restrict*(s >= limrest*dx) # remove lowest sizes (detection limit)
     46    restrict = restrict*(np.abs(i-j) <= 6.*dx) # condition sur i,j (width,height) pour ~round shape
     47    #restrict = restrict*(d > 0.9) # limit on drop for size
     48                                   # (casse la power law? seulement si keep smaller devils)
     49    #restrict = restrict*(d > 0.5)
     50    #restrict = restrict*(d > 0.3)
     51    poum=False
     52    if poum:
     53      out = var[np.logical_not(restrict)]
     54      outi = i[np.logical_not(restrict)]
     55      outj = j[np.logical_not(restrict)]
     56      print out, outi, outj
     57      print out.size
     58    var = var[restrict]
     59   
     60    ## define bins
     61    zebins = [np.min(var)]
     62    for i in range(0,nbins):  zebins.append(zebins[i]*(www**0.5))
     63    zebins = np.array(zebins)
     64    middle = 0.5*(zebins[1:] + zebins[:-1])
     65    binwidth = zebins[1:] - zebins[:-1]
     66    total = var.shape[0]
     67   
     68    # plot histogram
     69    yeah = mpl.hist(var,log=True,bins=zebins,normed=True,color='white') ; mpl.xscale('log')
     70    #yeah = mpl.hist(var,bins=zebins,normed=True,color='white')
     71   
     72    # add error bars
     73    incertitude_nvortex = 10.
     74    if not drop:
     75      err = incertitude_nvortex/(yeah[0]*binwidth*total+0.0001)
     76      ind = err < 1000. ; err = err[ind] ; y = yeah[0][ind] ; x = middle[ind]
     77      err = y*err
     78      mpl.errorbar(x,y,yerr=[err,err],fmt='k.')
     79   
     80    ## print info
     81    #print "min %5.2e // max %5.2e" % (np.min(var),np.max(var))
     82    #for iii in range(len(zebins)-1):
     83    #    this = zebins[iii] ; next = zebins[iii+1]
     84    #    print "%5.2e in [%5.0f %5.0f] %5.0f" % (yeah[0][iii],this,next,np.abs(next-this))
     85   
     86    # fitting function to superimpose
     87    if typefit == 1: xx = sciopt.curve_fit(fitfunc, middle, yeah[0])
     88    elif typefit == 2: xx = sciopt.curve_fit(fitfunc2, middle, yeah[0])
     89    elif typefit == 3: xx = sciopt.curve_fit(fitfunc3, middle, yeah[0])
     90    print "exponent",xx[0][1],"variance %",100.*xx[1][1][1]/xx[0][1]
     91   
     92    # label
     93    lablab = r"$\alpha=$%4.1f"%(xx[0][1])
     94    titi = addtitle+"LES "+str(int(dx))+"m. "+str(total)+" detected vortices"
     95   
     96    # plot obtained fit along with actual points
     97    if typefit == 1: func = fitfunc(middle,xx[0][0],xx[0][1])
     98    elif typefit == 2: func = fitfunc2(middle,xx[0][0],xx[0][1])
     99    elif typefit == 3: func = fitfunc3(middle,xx[0][0],xx[0][1],xx[0][2])
     100    func[func<0]=np.nan
     101    ind = yeah[0]>0
     102    mpl.plot(middle[ind],yeah[0][ind],'k.')
     103    mpl.plot(middle[ind],func[ind],'r-',label=lablab)
     104    mpl.plot(middle[ind],func[ind],'r.')
     105   
     106    # display bin limits
     107    tata = '[bins:'
     108    yorgl = ind[np.where(ind)].size + 1 # to add the last limit
     109    for el in zebins[0:yorgl]:
     110       tata = tata + "%0.f "%(el)
     111    tata=tata+']'
     112   
     113    # print fit vs. actual populations
     114    fit = func*binwidth*total
     115    pop = yeah[0]*binwidth*total
     116    for iii in range(len(fit)):
     117        if fit[iii] > 0.99 or pop[iii] != 0:
     118            yorgl = 100.*np.abs(fit[iii]-pop[iii])/fit[iii]
     119            yargl = 100.*incertitude_nvortex/fit[iii]
     120            print "fit %4.0f real %4.0f pc %3.0f ipc %3.0f" % (fit[iii],pop[iii],yorgl,yargl)
     121   
     122    # plot settings
     123    ax = mpl.gca() ; ax.set_xbound(lower=np.min(zebins),upper=np.max(zebins))
     124    if drop:
     125        mpl.xlabel("Pressure drop (Pa)")
     126        ax.set_xbound(lower=1.e-1,upper=10.)
     127    else:
     128        mpl.xlabel("Vortex size (m)")
     129        #ax.set_xbound(lower=30.,upper=3000.)
     130        ax.set_xbound(lower=30.,upper=1000.)
     131        ax.set_ybound(lower=5.e-6,upper=5.e-2)
     132    mpl.ylabel('Population density $n / N w_{bin}$')
     133    #mpl.title(titi+'\n '+tata)
     134    mpl.title(titi)
     135    mpl.legend(loc="upper right",fancybox=True)
     136       
     137    # show plot and end
     138    mpl.show()
    71139
    72 # choose a variable
    73 if drop: var = d
    74 else: var = s
     140#################################################################################
     141#################################################################################
     142#################################################################################
     143#
     144### UNCOMMENT
     145###import plfit,plpva,plplot
     146#
     147##[alpha, xmin, L] = plfit.plfit(var,'xmin',0.3)
     148#[alpha, xmin, L] = plfit.plfit(var)
     149#print alpha,xmin
     150#
     151##a = plpva.plpva(var,0.75,'xmin',0.75)
     152##print a
     153#
     154#
     155#h = plplot.plplot(var,xmin,alpha)
     156#
     157#ppplot.save(mode="png")
     158#
     159##mpl.loglog(h[0], h[1], 'k--',linewidth=2)
     160#
     161##mpl.show()
     162#
     163#
    75164
    76 # restrictions
    77 restrict = (s > 0) # initialization (True everywhere)
    78 restrict = restrict*(s >= limrest*dx) # remove lowest sizes (detection limit)
    79 #restrict = restrict*(np.abs(i-j) <= 6.*dx) # condition sur i,j (width,height)
    80 #restrict = restrict*(d > 0.9) # limit on drop for size (casse la power law? seulement si keep smaller devils)
    81 #restrict = restrict*(d > 0.5)
    82 var = var[restrict]
    83 
    84 ## define bins
    85 zebins = [np.min(var)]
    86 for i in range(0,nbins):  zebins.append(zebins[i]*(www**0.5))
    87 zebins = np.array(zebins)
    88 middle = 0.5*(zebins[1:] + zebins[:-1])
    89 binwidth = zebins[1:] - zebins[:-1]
    90 
    91 # plot histogram
    92 yeah = mpl.hist(var,log=True,bins=zebins,normed=True,color='white') ; mpl.xscale('log')
    93 #yeah = mpl.hist(var,bins=zebins)#,normed=True)
    94 
    95 # print info
    96 print "min %5.2e // max %5.2e" % (np.min(var),np.max(var))
    97 for iii in range(len(zebins)-1):
    98     print "%5.2e in [%5.0f %5.0f] %5.0f" % (yeah[0][iii],zebins[iii],zebins[iii+1],np.abs(zebins[iii+1]-zebins[iii]))
    99 
    100 # fitting function to superimpose
    101 if typefit == 1: xx = sciopt.curve_fit(fitfunc, middle, yeah[0])
    102 elif typefit == 2: xx = sciopt.curve_fit(fitfunc2, middle, yeah[0])
    103 elif typefit == 3: xx = sciopt.curve_fit(fitfunc3, middle, yeah[0])
    104 print "exponent",xx[0][1]
    105 print "variance %",100.*xx[1][1][1]/xx[0][1]
    106 
    107 # plot obtained fit along with actual points
    108 if typefit == 1: func = fitfunc(middle,xx[0][0],xx[0][1])
    109 elif typefit == 2: func = fitfunc2(middle,xx[0][0],xx[0][1])
    110 elif typefit == 3: func = fitfunc3(middle,xx[0][0],xx[0][1],xx[0][2])
    111 func[func<0]=np.nan
    112 ind = yeah[0]>0
    113 mpl.plot(middle[ind],yeah[0][ind],'k.')
    114 mpl.plot(middle[ind],func[ind],'r-')
    115 mpl.plot(middle[ind],func[ind],'r.')
    116 
    117 # print fit vs. actual populations
    118 total = var.shape[0]
    119 fit = func*binwidth*total
    120 pop = yeah[0]*binwidth*total
    121 for iii in range(len(fit)):
    122     if fit[iii] > 0.99 or pop[iii] != 0:
    123         print "fit %5.0f actual %5.0f" % (fit[iii],pop[iii])
    124 
    125 # plot settings
    126 ax = mpl.gca() ; ax.set_xbound(lower=np.min(zebins),upper=np.max(zebins))
    127 if drop: mpl.xlabel("Pressure drop (Pa)") ; ax.set_xbound(lower=1.e-1,upper=10.)
    128 else: mpl.xlabel("Vortex size (m)") ; ax.set_xbound(lower=10.,upper=5000.)
    129 mpl.ylabel('Population density $n / N w_{bin}$')
    130 mpl.title("Statistics on "+str(total)+" detected vortices")
    131 
    132 # show plot and end
    133 mpl.show()
    134 exit()
    135 
    136 ################################################################################
    137 ################################################################################
    138 ################################################################################
    139 
    140 ## UNCOMMENT
    141 ##import plfit,plpva,plplot
    142 
    143 #[alpha, xmin, L] = plfit.plfit(var,'xmin',0.3)
    144 [alpha, xmin, L] = plfit.plfit(var)
    145 print alpha,xmin
    146 
    147 #a = plpva.plpva(var,0.75,'xmin',0.75)
    148 #print a
    149 
    150 
    151 h = plplot.plplot(var,xmin,alpha)
    152 
    153 ppplot.save(mode="png")
    154 
    155 #mpl.loglog(h[0], h[1], 'k--',linewidth=2)
    156 
    157 #mpl.show()
    158 
    159 
    160 

  • trunk/UTIL/PYTHON/powerlaw/ddstat.py

    r1070 r1197  
    33import numpy as np
    44
    5 ##########################################################
    6 namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm1_2.txt"
    7 namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm2_2.txt"
    8 ##########################################################
    9 namefile = "/home/aymeric/Big_Data/LES_dd/psfc.LMD_LES_MARS.160564.ncm1_2.txt"
    10 #namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc.LMD_LES_MARS.160564.ncm2_2.txt"
    11 ##########################################################
    12 #namefile = "/home/aymeric/Big_Data/LES_dd/psfc_oldinsight100m.ncm1_2.txt"
    13 ##########################################################
     5### STATDD
     6def statdd(namefile):
    147
    15 
    16  
    17 
    18 #case = "188324p"
    19 #case = "191798"
    20 #case = "160564p"
    21 #case = "156487"
    22 #case = "2007p"
    23 #case = "13526p"
    24 #case = "172097"
    25 
    26 
    27 namefile = "/planeto/aslmd/LESdata/"+case+".ncm1_2.txt"
    28 
    29 
    30 
    31 
    32 
    33 # load data
    34 data = np.loadtxt(namefile,delimiter=";")
    35 t = data[:,0] ; n = data[:,1] ; s = data[:,2] ; d = data[:,3]
    36 
    37 # remove size and drop point when no vortex detected
    38 d[np.where(n==0)] = np.nan
    39 s[np.where(n==0)] = np.nan
    40 
    41 ## PLOTS
    42 
    43 number = ppplot.plot1d()
    44 number.f = n
    45 number.x = t
    46 number.linestyle = ''
    47 number.marker = '.'
    48 number.color = 'b'
    49 number.xlabel = "Local time (hour)"
    50 number.ylabel = "Detected vortices"
    51 number.makeshow()
    52 
    53 drop = ppplot.plot1d()
    54 drop.f = d
    55 drop.x = t
    56 drop.linestyle = ''
    57 drop.marker = '.'
    58 drop.color = 'r'
    59 drop.fmt = "%.1f"
    60 drop.xlabel = "Local time (hour)"
    61 drop.ylabel = "Maximum drop of detected vortices (Pa)"
    62 drop.makeshow()
    63 
    64 size = ppplot.plot1d()
    65 size.f = s
    66 size.x = t
    67 size.linestyle = ''
    68 size.marker = '.'
    69 size.color = 'g'
    70 size.xlabel = "Local time (hour)"
    71 size.ylabel = "Maximum size of detected vortices (m)"
    72 size.makeshow()
     8    # load data
     9    data = np.loadtxt(namefile,delimiter=";")
     10    t = data[:,0] ; n = data[:,1] ; s = data[:,2] ; d = data[:,3]
     11   
     12    # remove size and drop point when no vortex detected
     13    d[np.where(n==0)] = np.nan ; s[np.where(n==0)] = np.nan
     14   
     15    ## PLOTS
     16    number = ppplot.plot1d()
     17    number.f = n
     18    number.x = t
     19    number.linestyle = ''
     20    number.marker = '.'
     21    number.color = 'b'
     22    number.xlabel = "Local time (hour)"
     23    number.ylabel = "Detected vortices"
     24    number.makeshow()
     25   
     26    drop = ppplot.plot1d()
     27    drop.f = d
     28    drop.x = t
     29    drop.linestyle = ''
     30    drop.marker = '.'
     31    drop.color = 'r'
     32    drop.fmt = "%.1f"
     33    drop.xlabel = "Local time (hour)"
     34    drop.ylabel = "Maximum drop of detected vortices (Pa)"
     35    drop.makeshow()
     36   
     37    size = ppplot.plot1d()
     38    size.f = s
     39    size.x = t
     40    size.linestyle = ''
     41    size.marker = '.'
     42    size.color = 'g'
     43    size.xlabel = "Local time (hour)"
     44    size.ylabel = "Maximum size of detected vortices (m)"
     45    size.makeshow()

  • trunk/UTIL/PYTHON/powerlaw/finddd.py

    r1070 r1197  
    66import matplotlib.pyplot as mpl
    77import ppcompute
     8
     9## method 3
     10#from skimage import filter,transform,feature
     11#import matplotlib.patches as mpatches
     12###########
    813
    914###############################################################################
    … …  
    7176    mean = np.mean(psfc,dtype=np.float64)
    7277    std = np.std(psfc,dtype=np.float64)
     78    damax = np.max(psfc)
     79    damin = np.min(psfc)
    7380    ## some information about inferred limits
    7481    print "**************************************************************"
    … …  
    94101
    95102     ## MAIN ANALYSIS. LOOP ON FAC. OR METHOD.
    96      #for fac in faclist:
    97      fac = 3.75
    98      for method in [2,1]:
     103     for fac in faclist:
     104     #fac = 3.75
     105     #for method in [2,1]:
    99106 
    100107      ## initialize arrays
    101108      tabij = [] ; tabsize = [] ; tabdrop = []
    102109      tabijcenter = [] ; tabijvortex = [] ; tabijnotconv = [] ; tabdim = []
    103    
     110
    104111      ################ FIND RELEVANT POINTS TO BE ANALYZED
    105112      ## lab is 1 for points to be treated by minimum_position routine
    … …  
    110117          # method 1: standard deviation
    111118          lab[np.where(psfc2d < mean-fac*std)] = 1
    112       else:
     119      elif method == 2:
    113120          # method 2: polynomial fit
    114121          # ... tried smooth but too difficult, not accurate and too expensive
    … …  
    131138          limlim = fac*std ## same as method 1
    132139          lab[np.where(anopsfc2d < -limlim)] = 1
    133    
     140      elif method == 3:
     141          # method 3 : find centers of circle features using a Hough transform
     142          # initialize the array containing point to be further analyzed
     143          lab = np.zeros(psfc2d.shape)
     144          ## prepare the field to be analyzed by the Hough transform
     145          ## normalize it in an interval [-1,1]
     146          field = psfc2d
     147          # ... test 1. local max / min used for normalization.
     148          mmax = np.max(field) ; mmin = np.min(field) ; dasigma = 2.5
     149          ## ... test 2. global max / min used for normalization.
     150          #mmax = damax ; mmin = damin ; dasigma = 1.0 #1.5 trop restrictif
     151          spec = 2.*((field-mmin)/(mmax-mmin) - 0.5)
     152          # perform an edge detection on the field
     153          # ... returns an array with True on edges and False outside
     154          # http://sciunto.wordpress.com/2013/03/01/detection-de-cercles-par-une-transformation-de-hough-dans-scikit-image/   
     155          edges = filter.canny(filter.sobel(spec),sigma=dasigma)
     156          # initialize plot for checks
     157          if plotplot:
     158            fig, ax = mpl.subplots(ncols=1, nrows=1, figsize=(10,10))
     159            ax.imshow(field, cmap=mpl.cm.gray)
     160          # detect circle with radius 3dx. works well.
     161          # use an Hough circle transform
     162          radii = np.array([3])
     163          hough_res = transform.hough_circle(edges, radii)
     164          # analyze results of the Hough transform
     165          nnn = 0
     166          sigselec = neighbor_fac
     167          #sigselec = 3.
     168          for radius, h in zip(radii, hough_res):
     169            # number of circle features to keep
     170            # ... quite large. but we want to be sure not to miss anything.
     171            nup = 30
     172            maxima = feature.peak_local_max(h, num_peaks=nup)
     173            # loop on detected circle features
     174            for maximum in maxima:
     175              center_x, center_y = maximum - radii.max()
     176              # nup is quite high so there are false positives.
     177              # ... but those are easy to detect
     178              # ... if pressure drop is unclear (or inexistent)
     179              # ... we do not take the point into account for further analysis
     180              # ... NB: for inspection give red vs. green color to displayed circles
     181              diag = field[center_x,center_y] - (mean-sigselec*std)
     182              if diag < 0: 
     183                  col = 'green'
     184                  nnn = nnn + 1
     185                  lab[center_x,center_y] = 1
     186              else:
     187                  col = 'red'
     188              # draw circles
     189              if plotplot:
     190                circ = mpatches.Circle((center_y, center_x), radius,fill=False, edgecolor=col, linewidth=2)
     191                ax.add_patch(circ)
     192          if plotplot:
     193            mpl.title(str(nnn)+" vortices") ; mpl.show()
     194            mpl.close()
     195
    134196      ## while there are still points to be analyzed...
    135197      while 1 in lab:
    136198        ## ... get the point with the minimum field values
    137         if method == 1:
     199        if method == 1 or method == 3:
    138200            ij = minimum_position(psfc2d,labels=lab)
    139         else:
     201        elif method == 2:
    140202            ij = minimum_position(anopsfc2d,labels=lab)
    141203        ## ... store the indexes of the point in tabij
    … …  
    153215      ## --> or even lower as shown by plotting reslab
    154216      reslab = np.zeros(psfc2d.shape)
    155       if method == 1:
     217      if method == 1 or method == 3:
    156218          reslab[np.where(psfc2d < mean-neighbor_fac*std)] = 1
    157       else:
     219      elif method == 2:
    158220          reslab[np.where(anopsfc2d < -neighbor_fac*std)] = 1
    159221     
    … …  
    227289            ## OK. this is most likely an actual vortex. we get the drop.
    228290            ## we multiply by mesh area, then square to get approx. size of vortex
    229             if method == 1:
     291            if method == 1 or method ==3:
    230292                drop = -psfc2d[i,j]+mean
    231293            else:
    … …  
    301363   
    302364      #### PLOT PLOT PLOT PLOT
    303       if nvortex>0 and plotplot:
     365      if (nvortex>0 and plotplot) or (nvortex>0 and maxsize > 800.):
    304366       mpl.figure(figsize=(12,8))
    305367       myplot = plot2d()
    306        myplot.absc = np.array(range(psfc2d.shape[1]))*dx/1000.
    307        myplot.ordi = np.array(range(psfc2d.shape[0]))*dx/1000.
     368       myplot.x = np.array(range(psfc2d.shape[1]))*dx/1000.
     369       myplot.y = np.array(range(psfc2d.shape[0]))*dx/1000.
    308370       myplot.title = str(nvortex)+" vortices found (indicated diameter / pressure drop)"
    309371       myplot.xlabel = "x distance (km)"
    … …  
    312374       #if method == 1:
    313375           #myplot.field = ustm
    314            myplot.field = psfc2d
     376           myplot.f = psfc2d
    315377           #myplot.field = polypsfc2d
    316378           myplot.vmin = -2.*std
    … …  
    324386       myplot.fmt = "%.1f"
    325387       myplot.div = 20
    326        myplot.colorb = "spectral"
     388       myplot.colorbar = "spectral"
    327389       myplot.make()
    328390     

  • trunk/UTIL/PYTHON/powerlaw/performdd.py

    r1070 r1197  
    22from finddd import finddd
    33
    4 finddd("/planeto/aslmd/LESdata/case_A_50m_145_145_201_12km.nc",dt_out=100.,lt_start=8.)
    54
     5finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",method=3)
     6#finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",timelist=[250,385],method=3)
     7
     8#finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",method=2)
     9#finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",method=1)
     10
     11#finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",timelist=range(0,591,20),method=1,save=False,plotplot=True)
    612exit()
    713
Note: See TracChangeset for help on using the changeset viewer.