Changeset 647
- Timestamp:
- May 3, 2012, 6:46:55 PM (13 years ago)
- Location:
- trunk/UTIL/PYTHON
- Files:
-
- 5 edited
-
mymath.py (modified) (2 diffs)
-
myplot.py (modified) (11 diffs)
-
myscript.py (modified) (1 diff)
-
planetoplot.py (modified) (8 diffs)
-
pp.py (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/UTIL/PYTHON/mymath.py
r525 r647 58 58 else: 59 59 return np.array(field).sum(axis=axis) 60 61 def getmask (field): 62 import numpy as np 63 if field is None: return None 64 if type(field).__name__=='MaskedArray': 65 return np.ma.getmask(field) 66 else: 67 return np.isnan(field) 68 60 69 61 70 def deg (): … … 66 75 myfile = open(filename, 'w') 67 76 for line in mydata: 68 myfile.write(str(line) + '\n') 77 zeline = str(line) 78 zeline = zeline.replace('[','') 79 zeline = zeline.replace(']','') 80 myfile.write(zeline + '\n') 69 81 myfile.close() 70 82 return -
trunk/UTIL/PYTHON/myplot.py
r638 r647 55 55 ## Author: AS, AC, JL 56 56 def whatkindfile (nc): 57 typefile = 'gcm' # default 57 58 if 'controle' in nc.variables: typefile = 'gcm' 58 59 elif 'phisinit' in nc.variables: typefile = 'gcm' … … 62 63 elif hasattr(nc,'institution'): 63 64 if "European Centre" in getattr(nc,'institution'): typefile = 'ecmwf' 64 else: typefile = 'gcm' # for lslin-ed files from gcm65 65 return typefile 66 66 … … 165 165 ### ... note, anomaly is only computed over d1 and d2 for the moment 166 166 import numpy as np 167 from mymath import max,mean,min,sum 167 from mymath import max,mean,min,sum,getmask 168 168 csmooth = 12 ## a fair amount of grid points (too high results in high computation time) 169 169 if redope is not None: … … 194 194 elif max(d1) >= shape[1]: error = True 195 195 elif d1 is not None and d2 is not None: 196 totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0) 197 output = output*mesharea; output = sum(output[d2,:],axis=0); output = sum(output[d1],axis=0)/totalarea 196 totalarea = np.ma.masked_where(getmask(output),mesharea) 197 totalarea = mean(totalarea[d2,:],axis=0);totalarea = mean(totalarea[d1]) 198 output = output*mesharea; output = mean(output[d2,:],axis=0); output = mean(output[d1])/totalarea 198 199 elif d1 is not None: output = mean(input[:,d1],axis=1) 199 elif d2 is not None: totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[d2,:],axis=0)/totalarea 200 elif d2 is not None: 201 totalarea = np.ma.masked_where(getmask(output),mesharea) 202 totalarea = mean(totalarea[d2,:],axis=0) 203 output = output*mesharea; output = mean(output[d2,:],axis=0)/totalarea 200 204 elif dimension == 3: 201 205 if mesharea is None: mesharea=np.ones(shape[[1,2]]) … … 203 207 elif max(d2) >= shape[1]: error = True 204 208 elif max(d1) >= shape[2]: error = True 205 elif d4 is not None and d2 is not None and d1 is not None: 206 totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0) 209 elif d4 is not None and d2 is not None and d1 is not None: 207 210 output = mean(input[d4,:,:],axis=0) 208 output = output*mesharea; output = sum(output[d2,:],axis=0); output = sum(output[d1],axis=0)/totalarea 211 totalarea = np.ma.masked_where(getmask(output),mesharea) 212 totalarea = mean(totalarea[d2,:],axis=0);totalarea = mean(totalarea[d1]) 213 output = output*mesharea; output = mean(output[d2,:],axis=0); output = mean(output[d1])/totalarea 209 214 elif d4 is not None and d2 is not None: 210 215 output = mean(input[d4,:,:],axis=0) 211 totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[d2,:],axis=0)/totalarea 216 totalarea = np.ma.masked_where(getmask(output),mesharea) 217 totalarea = mean(totalarea[d2,:],axis=0) 218 output = output*mesharea; output = mean(output[d2,:],axis=0)/totalarea 212 219 elif d4 is not None and d1 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1) 213 220 elif d2 is not None and d1 is not None: 214 totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0) 215 output = output*mesharea; output = sum(output[:,d2,:],axis=1); output = sum(output[:,d1],axis=1)/totalarea 221 totalarea = np.tile(mesharea,(output.shape[0],1,1)) 222 totalarea = np.ma.masked_where(getmask(output),totalarea) 223 totalarea = mean(totalarea[:,d2,:],axis=1);totalarea = mean(totalarea[:,d1],axis=1) 224 output = output*mesharea; output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)/totalarea 216 225 elif d1 is not None: output = mean(input[:,:,d1],axis=2) 217 elif d2 is not None: totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,d2,:],axis=1)/totalarea 226 elif d2 is not None: 227 totalarea = np.tile(mesharea,(output.shape[0],1,1)) 228 totalarea = np.ma.masked_where(getmask(output),totalarea) 229 totalarea = mean(totalarea[:,d2,:],axis=1) 230 output = output*mesharea; output = mean(output[:,d2,:],axis=1)/totalarea 218 231 elif d4 is not None: output = mean(input[d4,:,:],axis=0) 219 232 elif dimension == 4: 220 if mesharea is None: mesharea=np.ones(shape[[2,3]]) 233 if mesharea is None: mesharea=np.ones(shape[[2,3]]) # mesharea=np.random.random_sample(shape[[2,3]])*5. + 2. # pour tester 221 234 if max(d4) >= shape[0]: error = True 222 235 elif max(d3) >= shape[1]: error = True … … 227 240 output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3) 228 241 if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) 229 totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0) 230 output = output*mesharea 231 output = sum(output[d2,:],axis=0) 232 output = sum(output[d1],axis=0)/totalarea 242 totalarea = np.ma.masked_where(np.isnan(output),mesharea) 243 totalarea = mean(totalarea[d2,:],axis=0); totalarea = mean(totalarea[d1]) 244 output = output*mesharea; output = mean(output[d2,:],axis=0); output = mean(output[d1])/totalarea 233 245 elif d4 is not None and d3 is not None and d2 is not None: 234 246 output = mean(input[d4,:,:,:],axis=0) 235 247 output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3) 236 248 if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) 237 totalarea=sum(mesharea[d2,:],axis=0) 238 output = output*mesharea; output = sum(output[d2,:],axis=0)/totalarea 249 totalarea = np.ma.masked_where(np.isnan(output),mesharea) 250 totalarea = mean(totalarea[d2,:],axis=0) 251 output = output*mesharea; output = mean(output[d2,:],axis=0)/totalarea 239 252 elif d4 is not None and d3 is not None and d1 is not None: 240 253 output = mean(input[d4,:,:,:],axis=0) … … 246 259 if anomaly: 247 260 for k in range(output.shape[0]): output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.) 248 totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0) 249 output = output*mesharea; output = sum(output[:,d2,:],axis=1); output = sum(output[:,d1],axis=1)/totalarea 261 totalarea = np.tile(mesharea,(output.shape[0],1,1)) 262 totalarea = np.ma.masked_where(getmask(output),mesharea) 263 totalarea = mean(totalarea[:,d2,:],axis=1); totalarea = mean(totalarea[:,d1],axis=1) 264 output = output*mesharea; output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)/totalarea 250 265 #noperturb = smooth1d(output,window_len=7) 251 266 #lenlen = len(output) ; output = output[1:lenlen-7] ; yeye = noperturb[4:lenlen-4] 252 267 #plot(output) ; plot(yeye) ; show() ; plot(output-yeye) ; show() 253 elif d3 is not None and d2 is not None and d1 is not None: 254 output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3) 268 elif d3 is not None and d2 is not None and d1 is not None: 269 output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3) 255 270 if anomaly: 256 271 for k in range(output.shape[0]): output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.) 257 totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0) 258 output = output*mesharea; output = sum(output[:,d2,:],axis=1); output = sum(output[:,d1],axis=1)/totalarea 272 totalarea = np.tile(mesharea,(output.shape[0],1,1)) 273 totalarea = np.ma.masked_where(getmask(output),totalarea) 274 totalarea = mean(totalarea[:,d2,:],axis=1); totalarea = mean(totalarea[:,d1],axis=1) 275 output = output*mesharea; output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)/totalarea 259 276 elif d4 is not None and d3 is not None: 260 277 output = mean(input[d4,:,:,:],axis=0) … … 262 279 if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) 263 280 elif d4 is not None and d2 is not None: 264 output = mean(input[d4,:,:,:],axis=0) 265 totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,d2,:],axis=1)/totalarea 281 output = mean(input[d4,:,:,:],axis=0) 282 totalarea = np.tile(mesharea,(output.shape[0],1,1)) 283 totalarea = np.ma.masked_where(getmask(output),mesharea) 284 totalarea = mean(totalarea[:,d2,:],axis=1) 285 output = output*mesharea; output = mean(output[:,d2,:],axis=1)/totalarea 266 286 elif d4 is not None and d1 is not None: 267 287 output = mean(input[d4,:,:,:],axis=0) 268 288 output = mean(output[:,:,d1],axis=2) 269 elif d3 is not None and d2 is not None: 289 elif d3 is not None and d2 is not None: 270 290 output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3) 271 totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,d2,:],axis=1)/totalarea 291 totalarea = np.tile(mesharea,(output.shape[0],1,1)) 292 totalarea = np.ma.masked_where(getmask(output),mesharea) 293 totalarea = mean(totalarea[:,d2,:],axis=1) 294 output = output*mesharea; output = mean(output[:,d2,:],axis=1)/totalarea 272 295 elif d3 is not None and d1 is not None: 273 296 output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3) 274 297 output = mean(output[:,:,d1],axis=2) 275 elif d2 is not None and d1 is not None: 276 totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0) 277 output = output*mesharea; output = sum(output[:,:,d2,:],axis=2); output = sum(output[:,:,d1],axis=2)/totalarea 298 elif d2 is not None and d1 is not None: 299 totalarea = np.tile(mesharea,(output.shape[0],output.shape[1],1,1)) 300 totalarea = np.ma.masked_where(getmask(output),totalarea) 301 totalarea = mean(totalarea[:,:,d2,:],axis=2); totalarea = mean(totalarea[:,:,d1],axis=1) 302 output = output*mesharea; output = mean(output[:,:,d2,:],axis=2); output = mean(output[:,:,d1],axis=1)/totalarea 278 303 elif d1 is not None: output = mean(input[:,:,:,d1],axis=3) 279 elif d2 is not None: 280 totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,:,d2,:],axis=2)/totalarea 304 elif d2 is not None: 305 totalarea = np.tile(mesharea,(output.shape[0],output.shape[1],1,output.shape[3])) 306 totalarea = np.ma.masked_where(getmask(output),totalarea) 307 totalarea = mean(totalarea[:,:,d2,:],axis=2) 308 output = output*mesharea; output = mean(output[:,:,d2,:],axis=2)/totalarea 281 309 elif d3 is not None: output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3) 282 310 elif d4 is not None: output = mean(input[d4,:,:,:],axis=0) … … 739 767 #zecolor = 'r' 740 768 #zelinewidth = 1 741 #step = 5.625742 #steplon = 5.625743 #zelatmax = 89.9769 #step = 180./48. 770 #steplon = 360./64. 771 #zelatmax = 90. - step/3 744 772 if char not in ["moll"]: 745 773 m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, latmax=zelatmax) … … 932 960 "CO2": "YlOrBr_r",\ 933 961 #### T.N. 934 "MTOT": " Jet",\962 "MTOT": "spectral",\ 935 963 "H2O_ICE_S": "RdBu",\ 936 964 "VMR_H2OICE": "PuBu",\ … … 1056 1084 return names 1057 1085 1058 ## Author: TN [old]1059 def getopmatrix (kind,n):1060 import numpy as np1061 # compute matrix of operations between files1062 # the matrix is 'number of files'-square1063 # 1: difference (row minus column), 2: add1064 # not 0 in diag : just plot1065 if n == 1:1066 opm = np.eye(1)1067 elif kind == 'basic':1068 opm = np.eye(n)1069 elif kind == 'difference1': # show differences with 1st file1070 opm = np.zeros((n,n))1071 opm[0,:] = 11072 opm[0,0] = 01073 elif kind == 'difference2': # show differences with 1st file AND show 1st file1074 opm = np.zeros((n,n))1075 opm[0,:] = 11076 else:1077 opm = np.eye(n)1078 return opm1079 1080 ## Author: TN [old]1081 def checkcoherence (nfiles,nlat,nlon,ntime):1082 if (nfiles > 1):1083 if (nlat > 1):1084 errormess("what you asked is not possible !")1085 return 11086 1086 1087 1087 ## Author: TN -
trunk/UTIL/PYTHON/myscript.py
r637 r647 74 74 ### SPECIAL 75 75 parser.add_option('--tsat', action='store_true',dest='tsat', default=False,help='convert temperature field T in Tsat-T using pressure') 76 parser.add_option('--stream', action='store_true',dest='stream', default=False,help='plot streamlines from streamfunction.e for vert/lat slices.') 76 77 77 78 return parser -
trunk/UTIL/PYTHON/planetoplot.py
r642 r647 66 66 lbls=None,\ 67 67 lstyle=None,\ 68 cross=None): 68 cross=None,\ 69 streamflag=False): 69 70 70 71 #################################################################################################################### … … 91 92 #from singlet import singlet 92 93 from itertools import cycle 94 import os 93 95 94 96 ######################### … … 414 416 if mult != 2718.: what_I_plot = what_I_plot*mult 415 417 else: what_I_plot = np.log10(what_I_plot) ; print "log plot" 418 416 419 if var2: ### what is contoured. 417 420 what_I_plot_contour, error = reducefield( all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert, \ … … 523 526 if indexvert is not None: lbl = lbl + " iz" + str(indexvert[0]) 524 527 if indextime is not None: lbl = lbl + " it" + str(indextime[0]) 525 if lbl == "": lbl = namefiles[index_f]528 if lbl == "": lbl = all_namefile[index_f] 526 529 527 530 if mrate is not None: x = y ## because swapaxes... … … 539 542 540 543 elif which == "regular": 544 545 # plot stream lines if there is a stream file and a vert/lat slice. Might not work with movies ?? 546 if streamflag and sslat is None and svert is None: 547 streamfile = all_namefile[index_f].replace('.nc','_stream.nc') 548 teststream = os.path.exists(streamfile) 549 if teststream: 550 print 'INFO: Using stream file',streamfile, 'for stream lines' 551 ncstream = Dataset(streamfile) 552 psi = getfield(ncstream,'psi') 553 psi = psi[0,:,:,0] # time and longitude are dummy dimensions 554 if psi.shape[1] != len(x) or psi.shape[0] != len(y): 555 errormess('stream file does not fit! Dimensions: '+str(psi.shape)+' '+str(x.shape)+' '+str(y.shape)) 556 zelevels = np.arange(-1.e10,1.e10,1.e9) 557 zemin = np.min(abs(zelevels)) 558 zemax = np.max(abs(zelevels)) 559 zewidth = (abs(zelevels)-zemin)*(5.- 0.5)/(zemax - zemin) + 0.5 # linewidth ranges from 5 to 0.5 560 cs = contour( x,y,psi, zelevels, colors='k', linewidths = zewidth) 561 clabel(cs, inline=True, fontsize = 4.*rcParams['font.size']/5., fmt="%1.1e") 562 else: 563 print 'WARNING: STREAM FILE',streamfile, 'DOES NOT EXIST !' 564 541 565 if hole: what_I_plot_frame = hole_bounds(what_I_plot_frame,zevmin,zevmax) 542 566 else: what_I_plot_frame = bounds(what_I_plot_frame,zevmin,zevmax) … … 554 578 if mapmode == 1: m.pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans) 555 579 elif mapmode == 0: pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans) 580 556 581 557 582 if colorb not in ['nobar','onebar']: … … 591 616 zebins = np.linspace(zevmin,zevmax,num=30) 592 617 hist(toplot,bins=zebins,histtype='step',linewidth=2,color='k',normed=True) 593 zestd = np.std(toplot) 594 zemean = np.mean(toplot) 618 zestd = np.std(toplot);zemean = np.mean(toplot) 595 619 zebins = np.linspace(zevmin,zevmax,num=300) 596 620 zegauss = (1./(zestd * np.sqrt(2 * np.pi)) * np.exp( - (zebins - zemean)**2 / (2 * zestd**2) ) ) … … 675 699 if index_f is numplot-1: plottitle = basename+' '+"fig(1) "+ope+" fig(2)" 676 700 elif index_f is numplot-2: plottitle = basename+' '+fileref 677 else: plottitle = basename+' '+ namefiles[0]#index_f]678 else: plottitle = basename+' '+ namefiles[0]#index_f]701 else: plottitle = basename+' '+all_namefile[index_f] 702 else: plottitle = basename+' '+all_namefile[index_f] 679 703 if mult != 1: plottitle = '{:.0e}'.format(mult) + "*" + plottitle 680 704 if zetitle[0] != "fill": -
trunk/UTIL/PYTHON/pp.py
r612 r647 172 172 mrate=opt.rate,mquality=opt.quality,trans=opt.trans,zarea=opt.area,axtime=opt.axtime,\ 173 173 redope=opt.redope,seevar=opt.seevar,xlab=opt.xlab,ylab=opt.ylab,lbls=separatenames(opt.labels),\ 174 lstyle=separatenames(opt.linestyle),cross=readslices(opt.mark) )174 lstyle=separatenames(opt.linestyle),cross=readslices(opt.mark),streamflag=opt.stream) 175 175 print 'DONE: '+name 176 176 system("rm -f to_be_erased")
Note: See TracChangeset
for help on using the changeset viewer.
