- Timestamp:
- Sep 29, 2011, 3:59:01 PM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/examples/scripts/find_devils.py
r251 r310 23 23 from scipy.ndimage.measurements import minimum_position 24 24 from netCDF4 import Dataset 25 25 import matplotlib.pyplot as plt 26 import myplot as myp 27 26 28 ### LOAD NETCDF DATA 27 29 filename = "psfc.nc" … … 34 36 allsizesx = [] 35 37 allsizesy = [] 36 for i in range(0,shape[0]): 37 #for i in range(0,2): 38 39 print i, ' on ', shape[0] 38 stride = 1 #5 39 for i in range(0,shape[0],stride): 40 40 41 psfc2d = np.array ( psfc [ i, : , : ] ) 41 42 42 43 ############### CRITERION 43 ave = np.mean(psfc2d,dtype=np.float64) ## dtype otherwise inaccuracy 44 where = np.where(psfc2d - ave < -0.3) ## comme le papier Phoenix 44 ave = np.mean(psfc2d,dtype=np.float64) ## dtype otherwise inaccuracy 45 46 #limdp = -0.2 ## on en loupe pas mal 47 #where = np.where(psfc2d - ave < limdp) ## comme le papier Phoenix 45 48 46 std = np.std(psfc2d,dtype=np.float64) ## dtype otherwise inaccuracy 47 fac = 2.5 ## not too low, otherwise vortices are not caught. 4 good choice. 49 std = np.std(psfc2d,dtype=np.float64) ## dtype otherwise inaccuracy 50 fac = 4. ## how many sigmas. not too low, otherwise vortices are not caught. 4 good choice. 51 ## 2.5 clearly too low, 3.5 not too bad, 4 probably good 52 fac = 3.5 53 #fac = 2.5 48 54 lim = ave - fac*std 49 print lim, ave, std50 55 where = np.where(psfc2d < lim) 51 56 ############### END CRITERION 52 57 53 58 lab = np.zeros(np.array(psfc2d).shape) ## points to be treated by the minimum_position routine 54 59 lab[where] = 1. ## do not treat points close to 'mean' (background) pressure 55 60 61 draw = False 62 if draw: 63 ################################################################################## 64 vmin = -0.3 65 vmax = 0.0 66 ndiv = 3 67 palette = plt.get_cmap(name="YlGnBu") 68 what_I_plot = psfc2d-ave 69 zevmin, zevmax = myp.calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax) 70 what_I_plot = myp.bounds(what_I_plot,zevmin,zevmax) 71 zelevels = np.linspace(zevmin,zevmax) 72 fig = plt.figure() 73 sub = myp.definesubplot(2,fig) 74 plt.subplot(sub) 75 plt.contourf(what_I_plot,zelevels,cmap=palette) 76 plt.colorbar(fraction=0.05,pad=0.03,format="%.1f",\ 77 ticks=np.linspace(zevmin,zevmax,ndiv+1),\ 78 extend='both',spacing='proportional') 79 plt.subplot(sub+1) 80 palette = plt.get_cmap(name="binary") 81 plt.contourf(lab,2,cmap=palette) 82 plt.show() 83 ################################################################################## 84 56 85 xx = [] 57 86 yy = [] … … 69 98 sizex = detsize( xx, res = 10, loga=False, thres=2 ) 70 99 sizey = detsize( yy, res = 10, loga=False, thres=2 ) 100 #sizex = detsize( xx, res = 10, loga=False, thres=3 ) 101 #sizey = detsize( yy, res = 10, loga=False, thres=3 ) 71 102 ### 103 #if ( mym.max(sizex) > mym.max(sizey) ): sizey = sizex ### un peu limite dans certains cas 104 if (len(sizex) > len(sizey)) : sizey = sizex ### plus fidele mais petit souci lorsque PBC 105 elif (len(sizex) == len(sizey)) : 106 if ( mym.max(sizex) > mym.max(sizey) ): sizey = sizex 107 else : sizex = sizey 108 else : sizex = sizey 72 109 allsizesx = np.append(allsizesx,sizex) 73 110 allsizesy = np.append(allsizesy,sizey) 111 print i, ' on ', shape[0], ' caught ', len(sizex), ' vortices ', sizex 74 112 ######## 75 113 76 114 allsizesx.sort() 77 115 allsizesy.sort() … … 79 117 return allsizesx, allsizesy 80 118 81 def writeascii ( tab, filename ): 82 mydata = tab 83 myfile = open(filename, 'w') 84 for line in mydata: 85 myfile.write(str(line) + '\n') 86 myfile.close() 87 return 119 ######################################################################### 120 ######################################################################### 88 121 89 122 import matplotlib.pyplot as plt … … 91 124 import numpy as np 92 125 import matplotlib.mlab as mlab 126 import mymath as mym 93 127 94 128 save = True 95 #save = False 129 save = False 130 96 131 if save: 97 132 allsizesx, allsizesy = getsize("psfc.nc") 98 writeascii(allsizesx,'allsizex.txt') 99 writeascii(allsizesy,'allsizey.txt') 100 133 ### sauvegarde texte pour inspection 134 mym.writeascii(allsizesx,'allsizex.txt') 135 mym.writeascii(allsizesy,'allsizey.txt') 136 ### sauvegarde binaire pour utilisation python 101 137 myfile = open('allsizex.bin', 'wb') 102 138 pickle.dump(allsizesx, myfile) 103 139 myfile.close() 104 105 140 myfile = open('allsizey.bin', 'wb') 106 141 pickle.dump(allsizesy, myfile) 107 142 myfile.close() 108 143 109 144 ### load files 110 145 myfile = open('allsizex.bin', 'r') 111 146 allsizesx = pickle.load(myfile) … … 113 148 allsizesy = pickle.load(myfile) 114 149 150 ### append sizes 115 151 plothist = np.append(allsizesx,allsizesy) 152 plothist = allsizesx 153 #plothist = allsizesy 116 154 plothist.sort() 117 155 118 nbins=100 156 ### make bins 157 nbins = 100 119 158 zebins = [2.0] 159 #nbins = 8 160 #zebins = [19.0] 161 #nbins = 15 162 #zebins = [11.] ##12 non mais donne un peu la meme chose 163 #zebins = [20.] ##20 non car trop pres du premier 164 nbins = 100 165 zebins = [2./np.sqrt(2.)] ## ne pas tomber sur une dizaine ronde 166 nbins = 200 167 zebins = [0.2/np.sqrt(2.)] ## ne pas tomber sur une dizaine ronde 168 120 169 for i in range(0,nbins): zebins.append(zebins[i]*np.sqrt(2)) 121 170 zebins = np.array(zebins) 122 print zebins 123 171 ### select reasonable bins for DD 124 172 zebins = zebins [ zebins > 10. ] 125 173 zebins = zebins [ zebins < 1000. ] 126 print zebins 174 #print 'bins ',zebins 175 print 'corrected bins ',zebins 176 print 'max value ',mym.max(plothist),mym.max(allsizesx),mym.max(allsizesy) 177 178 #plothist = [20.,30.,40.,50.,70.,100.,130.,190.,260.,370.,520.] ## pour calibrer 127 179 128 180 #### HISTOGRAM … … 130 182 log=True,\ 131 183 bins=zebins,\ 132 #cumulative=-1,\184 # cumulative=-1,\ 133 185 ) 134 186 plt.xscale('log') 135 136 187 plt.show() 137 188
Note: See TracChangeset
for help on using the changeset viewer.
