Changeset 1197
- Timestamp:
- Mar 4, 2014, 12:00:10 PM (11 years ago)
- Location:
- trunk/UTIL/PYTHON
- Files:
-
- 1 added
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/UTIL/PYTHON/mcd/examples/dustdevil.py
r1007 r1197 9 9 10 10 from ppplot import writeascii 11 12 13 from scipy import stats 11 14 12 15 … … 53 56 ## constant_heights 54 57 #for heights in [20.,100.,500.,1000.,2000.,5000.,7000.,10000.,20000.]: 55 for heights in [10.]:58 #for heights in [10.]: 56 59 57 60 ## multiple_obsheight 58 61 #for heights in [0.1,0.2,0.5,0.75,1.0,1.5,2.0,10.]: 62 for heights in [1.0]: 59 63 60 64 query = mcd() … … 72 76 73 77 ## constant_heights 74 query.xz = heights 75 78 #query.xz = heights 76 79 ## multiple_obsheight 77 #query.xz = ddheight[i-1]*heights80 query.xz = ddheight[i-1]*heights 78 81 79 82 query.update() ; uu = query.zonwind ; vv = query.merwind … … 87 90 if angle < 0.: angle = angle + 360. 88 91 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 93 98 94 99 query.profile() … … 130 135 writeascii(field=mcdwind,absc=ddwind,name=str(heights)+'.txt') 131 136 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) 132 147 133 148 mpl.errorbar(mcdwind, ddwind, yerr=[ddwinderror,ddwinderror], xerr=[mcdwindle,mcdwindhe], fmt='bo') -
trunk/UTIL/PYTHON/powerlaw/ddhist.py
r1070 r1197 4 4 import numpy as np 5 5 import scipy.optimize as sciopt 6 7 ### SAVED8 #namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm2_1.txt" ; drop = False ; typefit=19 #namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm2_1.txt" ; drop = True ; typefit=210 #namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm1_1.txt" ; drop = False ; typefit=111 #namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm1_1.txt" ; drop = True ; typefit=112 #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=138 #drop = True ; typefit=139 #drop = True ; typefit=240 ##########################################################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.544 #nbins = 15 ; www = 2.45 ##nbins = 20 ; www = 1.7546 ##nbins = 30 ; www = 1.547 ##nbins = 50 ; www = 1.248 ##########################################################49 limrest = 4. # restrict limit (multiple of dx)50 #limrest = 2.51 #limrest = 3.52 #limrest = 0.53 ##########################################################54 6 55 7 # -- functions for fitting … … 62 14 def fitfunc2(x,a,b): return a*np.exp(-b*x) 63 15 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 ################################################################################ 25 def histodd(namefile,drop=False,typefit=1,nbins=12,limrest=4,addtitle=""): 68 26 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() 71 139 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 # 75 164 76 # restrictions77 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 bins85 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 histogram92 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 info96 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 superimpose101 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 points108 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.nan112 ind = yeah[0]>0113 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 populations118 total = var.shape[0]119 fit = func*binwidth*total120 pop = yeah[0]*binwidth*total121 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 settings126 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 end133 mpl.show()134 exit()135 136 ################################################################################137 ################################################################################138 ################################################################################139 140 ## UNCOMMENT141 ##import plfit,plpva,plplot142 143 #[alpha, xmin, L] = plfit.plfit(var,'xmin',0.3)144 [alpha, xmin, L] = plfit.plfit(var)145 print alpha,xmin146 147 #a = plpva.plpva(var,0.75,'xmin',0.75)148 #print a149 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 3 3 import numpy as np 4 4 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 6 def statdd(namefile): 14 7 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 6 6 import matplotlib.pyplot as mpl 7 7 import ppcompute 8 9 ## method 3 10 #from skimage import filter,transform,feature 11 #import matplotlib.patches as mpatches 12 ########### 8 13 9 14 ############################################################################### … … 71 76 mean = np.mean(psfc,dtype=np.float64) 72 77 std = np.std(psfc,dtype=np.float64) 78 damax = np.max(psfc) 79 damin = np.min(psfc) 73 80 ## some information about inferred limits 74 81 print "**************************************************************" … … 94 101 95 102 ## MAIN ANALYSIS. LOOP ON FAC. OR METHOD. 96 #for fac in faclist:97 fac = 3.7598 for method in [2,1]:103 for fac in faclist: 104 #fac = 3.75 105 #for method in [2,1]: 99 106 100 107 ## initialize arrays 101 108 tabij = [] ; tabsize = [] ; tabdrop = [] 102 109 tabijcenter = [] ; tabijvortex = [] ; tabijnotconv = [] ; tabdim = [] 103 110 104 111 ################ FIND RELEVANT POINTS TO BE ANALYZED 105 112 ## lab is 1 for points to be treated by minimum_position routine … … 110 117 # method 1: standard deviation 111 118 lab[np.where(psfc2d < mean-fac*std)] = 1 112 el se:119 elif method == 2: 113 120 # method 2: polynomial fit 114 121 # ... tried smooth but too difficult, not accurate and too expensive … … 131 138 limlim = fac*std ## same as method 1 132 139 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 134 196 ## while there are still points to be analyzed... 135 197 while 1 in lab: 136 198 ## ... get the point with the minimum field values 137 if method == 1 :199 if method == 1 or method == 3: 138 200 ij = minimum_position(psfc2d,labels=lab) 139 el se:201 elif method == 2: 140 202 ij = minimum_position(anopsfc2d,labels=lab) 141 203 ## ... store the indexes of the point in tabij … … 153 215 ## --> or even lower as shown by plotting reslab 154 216 reslab = np.zeros(psfc2d.shape) 155 if method == 1 :217 if method == 1 or method == 3: 156 218 reslab[np.where(psfc2d < mean-neighbor_fac*std)] = 1 157 el se:219 elif method == 2: 158 220 reslab[np.where(anopsfc2d < -neighbor_fac*std)] = 1 159 221 … … 227 289 ## OK. this is most likely an actual vortex. we get the drop. 228 290 ## we multiply by mesh area, then square to get approx. size of vortex 229 if method == 1 :291 if method == 1 or method ==3: 230 292 drop = -psfc2d[i,j]+mean 231 293 else: … … 301 363 302 364 #### PLOT PLOT PLOT PLOT 303 if nvortex>0 and plotplot:365 if (nvortex>0 and plotplot) or (nvortex>0 and maxsize > 800.): 304 366 mpl.figure(figsize=(12,8)) 305 367 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. 308 370 myplot.title = str(nvortex)+" vortices found (indicated diameter / pressure drop)" 309 371 myplot.xlabel = "x distance (km)" … … 312 374 #if method == 1: 313 375 #myplot.field = ustm 314 myplot.f ield= psfc2d376 myplot.f = psfc2d 315 377 #myplot.field = polypsfc2d 316 378 myplot.vmin = -2.*std … … 324 386 myplot.fmt = "%.1f" 325 387 myplot.div = 20 326 myplot.colorb = "spectral"388 myplot.colorbar = "spectral" 327 389 myplot.make() 328 390 -
trunk/UTIL/PYTHON/powerlaw/performdd.py
r1070 r1197 2 2 from finddd import finddd 3 3 4 finddd("/planeto/aslmd/LESdata/case_A_50m_145_145_201_12km.nc",dt_out=100.,lt_start=8.)5 4 5 finddd("/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) 6 12 exit() 7 13
Note: See TracChangeset
for help on using the changeset viewer.