Changeset 451
- Timestamp:
- Dec 5, 2011, 2:42:33 AM (13 years ago)
- Location:
- trunk/UTIL/PYTHON
- Files:
-
- 2 added
- 5 edited
-
README.PP (modified) (2 diffs)
-
body.html (added)
-
header.html (added)
-
myplot.py (modified) (8 diffs)
-
myscript.py (modified) (2 diffs)
-
planetoplot.py (modified) (17 diffs)
-
pp.py (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/UTIL/PYTHON/README.PP
r448 r451 12 12 ************************************** 13 13 14 15 ************************************* 16 MAPMODE 1 17 MAPPING MODE 18 SIMPLE EXAMPLES on a SAMPLE GCM FILE 19 ************************************* 20 14 ***************************************************************** 15 MAPMODE 1... MAPPING MODE... SIMPLE EXAMPLES on a SAMPLE GCM FILE 16 ***************************************************************** 21 17 Goal: The simplest, most minimal example. Mapping topography. 22 18 pp.py -f diagfired.nc … … 73 69 pp.py -f diagfi.nc -v tsurf --time 4,7 74 70 75 --- mesoscale stuff 71 [only mesoscale for the moment] 76 72 Goal: I want to plot results for two different LOCAL times in the file next to one another 77 73 pp.py -f wrfout**** -v TSURF --time -4 -- time -7 78 74 79 Goal: The classic mountain GW plot 75 *********************************************************************************** 76 EXAMPLE : The classic mountain GW plot 77 *********************************************************************************** 80 78 pp.py -f wrfout_d01_9999-09-09_09:00:00 -v W,tpot --lat 60 --time 15 -i 4 -l 30,130,100 --div 50 79 *********************************************************************************** 80 81 *********************************************************************************** 82 COMMENTED EXAMPLE : The globe with surface temperature and winds 83 *********************************************************************************** 84 pp.py -f diagfired.nc -v tsurf -w phisinit -m 120 -M 320 --div 20 -W -s 1 --vert 0 -p ortho --blat 20 --blon -80 -S html -t $W 85 *********************************************************************************** 86 See results here: http://www.lmd.jussieu.fr/~aslmd/EXAMPLES/LMD_GCM_movie_tsurf_UV/anim.html 87 *********************************************************************************** 88 pp.py -f diagfired.nc 89 OK. You probably get that one. 90 -v tsurf -w phisinit 91 Shade surface temperature. Contour topography. 92 -m 120 -M 320 --div 20 93 Surface temperature is shown with bounds 120K to 320K. Use 20 levels for shading. 94 -W -s 1 95 Include wind vectors. Prescribe a stride of 1: vectors are shown at every grid point. 96 --vert 0 97 Show fields in the first (lowermost) level. 98 -p ortho --blat 20 --blon -80 99 Use orthographic projection ('whole sphere' view). Center view on lon -80E and lat 20N. 100 -S html 101 Make nice webpage with animation and controls. 102 -t /u/aslmd/WWW/EXAMPLES 103 Move resulting plot files to the given folder. 104 105 *********************************************************************************** 106 COMMENTED EXAMPLE : The dust storm section movie 107 *********************************************************************************** 108 pp.py -f wrfout_d01_2024-05-30_12:00:00,wrfout_d01_2024-05-30_18:00:00,wrfout_d01_2024-05-31_00\:00\:00 --operation cat -v QDUST --lat -3. -i 3 -l -1,37,100 --div 30 -c Oranges_r -m 0. -M 5.e-5 -t $W --rate 12 --xmin=5 --xmax=115 109 *********************************************************************************** 110 See results here: http://www.lmd.jussieu.fr/~aslmd/EXAMPLES/LMD_MMM_d1_10km_movie_QDUST_-1000m-AMR_lat_-3_Ls134.8/anim.html 111 *********************************************************************************** 112 pp.py 113 OK. You probably get that one. 114 -f wrfout_d01_2024-05-30_12:00:00,wrfout_d01_2024-05-30_18:00:00,wrfout_d01_2024-05-31_00\:00\:00 --operation cat 115 Mesoscale outputs are splitted in several files. In that case, this was 1 file per 6 simulated hours. 116 The above options allow to concatenate files along time axis for 1D time series or animated movies 117 [if --operation cat is omitted, this is a multiplot call, with one subplot per files in -f]. 118 -v QDUST 119 Choose to plot dust mass mixing ratio. 120 -i 3 -l -1,37,100 121 Set a call to vertical interpolator [compiled with f2py, thereby being embedded as a Python routine] for each of the files in the -f instance. 122 -i sets the kind of interpolation, 3 means Above MOLA Reference Altitude. -l sets the range for altitude levels: from -1 km to 37 km with 100 levels. 123 --div 30 124 The number of contours used for shaded plots. Higher value means smoother appearance. 125 -c Oranges_r 126 Choose a colorbar adapted to display a dust storm. 127 -m 0. -M 5.e-5 128 Choose bounds for the plotted field. This one is adapted to show dust mass mixing ratio. 129 -t $W 130 Put resulting figure or movie in another destination folder. 131 Personally I have an environnement variable W which is somewhere in my system where the file automatically appears on the web, 132 hence is easy to see from a remote place. 133 --rate 12 --lat -3. 134 Define prescribed axis. A section in latitude -3°N. A time animation with 12 frame per seconds. 135 So the displayed field will be an altitude/longitude section. Alternative: "-S avi" instead of "--rate 12" creates a default 8 fps movie. 136 Alternative II: "-S html" instead of --rate 12 creates a nice webpage. 137 --xmin=5 --xmax=115 138 Define limits for the displayed section. Here we just want to get rid of transition rows where atmospheric fields are 139 relaxed towards prescribed GCM fields. 81 140 82 141 *********************************************************************************** -
trunk/UTIL/PYTHON/myplot.py
r448 r451 376 376 377 377 ## Author: AS + AC 378 def dumpbdy (field,n,stag=None,condition=False ):378 def dumpbdy (field,n,stag=None,condition=False,onlyx=False,onlyy=False): 379 379 nx = len(field[0,:])-1 380 380 ny = len(field[:,0])-1 … … 383 383 if stag == 'V': ny = ny-1 384 384 if stag == 'W': nx = nx+1 #special les case when we dump stag on W 385 return field[n:ny-n,n:nx-n] 385 if onlyx: result = field[:,n:nx-n] 386 elif onlyy: result = field[n:ny-n,:] 387 else: result = field[n:ny-n,n:nx-n] 388 return result 386 389 387 390 ## Author: AS + AC … … 547 550 548 551 ## Author: AS 549 def define_proj (char,wlon,wlat,back=None,blat=None ):552 def define_proj (char,wlon,wlat,back=None,blat=None,blon=None): 550 553 from mpl_toolkits.basemap import Basemap 551 554 import numpy as np … … 561 564 elif wlat[0] <= 0.: blat = wlat[1] 562 565 else: ortholat=blat 563 #print "blat ", blat 566 if blon is None: ortholon=meanlon 567 else: ortholon=blon 564 568 h = 50. ## en km 565 569 radius = 3397200. … … 567 571 llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) 568 572 elif char == "moll": m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon) 569 elif char == "ortho": m = Basemap(rsphere=radius,projection='ortho',lon_0= meanlon,lat_0=ortholat)573 elif char == "ortho": m = Basemap(rsphere=radius,projection='ortho',lon_0=ortholon,lat_0=ortholat) 570 574 elif char == "lcc": m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\ 571 575 llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) … … 645 649 if zevmin < 0: what_I_plot[ what_I_plot < zevmin*(1. - 1.e-7) ] = zevmin*(1. - 1.e-7) 646 650 else: what_I_plot[ what_I_plot < zevmin*(1. + 1.e-7) ] = zevmin*(1. + 1.e-7) 647 print "NEW MIN ", min(what_I_plot)651 #print "NEW MIN ", min(what_I_plot) 648 652 what_I_plot[ what_I_plot > 9e+35 ] = -9e+35 649 653 what_I_plot[ what_I_plot > zevmax ] = zevmax 650 print "NEW MAX ", max(what_I_plot)654 #print "NEW MAX ", max(what_I_plot) 651 655 return what_I_plot 652 656 … … 707 711 "ICETOT": "%.1e",\ 708 712 "TAU_ICE": "%.2f",\ 713 "TAUICE": "%.2f",\ 709 714 "VMR_ICE": "%.1e",\ 710 715 "MTOT": "%.1f",\ … … 751 756 "TEMP": "Jet",\ 752 757 "TAU_ICE": "Blues",\ 758 "TAUICE": "Blues",\ 753 759 "VMR_ICE": "Blues",\ 754 760 "W": "jet",\ -
trunk/UTIL/PYTHON/myscript.py
r444 r451 4 4 parser.add_option('-f', '--file', action='append',dest='file', type="string", default=None, help='[NEEDED] filename. Append: different figures. Comma-separated: same figure (+ possible --operation). Regex tip: -f `echo foo* | sed s+" "+","+g`') 5 5 parser.add_option('-t', '--target', action='store',dest='tgt', type="string", default=None, help='destination folder') 6 parser.add_option('-S', '--save', action='store',dest='save', type="string", default="gui", help='save mode ( png,eps,svg,pdf,txt or gui) [gui]')6 parser.add_option('-S', '--save', action='store',dest='save', type="string", default="gui", help='save mode (gui,png,eps,svg,pdf,txt,html,avi) [gui]') 7 7 parser.add_option('-d', '--display',action='store_false',dest='display', default=True, help='do not pop up created images') 8 8 parser.add_option('-O','--output', action='store',dest='out', type="string", default=None, help='output file name') … … 37 37 parser.add_option('-W', '--winds', action='store_true',dest='winds', default=False, help='wind vectors [False]') 38 38 parser.add_option('-z', '--zoom', action='store',dest='zoom', type="float", default=None, help='zoom factor in %') 39 parser.add_option('--blat', action='store',dest='blat', type="int", default=None, help='bounding latitude for stereographic plots [computed]') 39 parser.add_option('--blat', action='store',dest='blat', type="int", default=None, help='reference lat (or bounding lat for stere) [computed]') 40 parser.add_option('--blon', action='store',dest='blon', type="int", default=None, help='reference lon [computed]') 40 41 41 42 ### SPECIFIC FOR SLICING [MAPMODE 0] -
trunk/UTIL/PYTHON/planetoplot.py
r448 r451 8 8 ### A. Spiga -- LMD -- 11~12/2011 -- Extended multivar subplot capabilities + cosmetic changes + general cleaning and tests 9 9 ### A. Colaitis -- LMD -- 12/2011 -- Added movie capability [mencoder must be installed] 10 ### A. Spiga -- LMD -- 12/2011 -- Added HTML animated page capability + general tests of consistency [winds, etc...] 10 11 11 12 def planetoplot (namefiles,\ … … 50 51 yintegral=False,\ 51 52 blat=None,\ 53 blon=None,\ 52 54 tsat=False,\ 53 55 flagnolow=False,\ … … 68 70 from mymath import deg,max,min,mean,get_tsat,writeascii,fig2data,fig2img 69 71 import matplotlib as mpl 70 from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel, title 72 from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show, plot, clabel, title, close 71 73 from matplotlib.cm import get_cmap 72 74 import numpy as np … … 97 99 stime = readslices(str(0)) ; ntime=1 ## this is a default choice 98 100 print "WELL... nothing about time axis. I took default: first time reference stored in file." 101 if mrate is not None and len(var) > 1: errormess("multivar not allowed in movies. should be fixed soon!") 99 102 zelen = len(namefiles)*len(var) 100 103 numplot = zelen*nslices … … 103 106 if fileref is not None: zelen = zelen + 2 104 107 elif "var" in ope: zelen = zelen + 1 105 all_var = [[]]*zelen ; all_var2 = [[]]*zelen ; all_title = [[]]*zelen ; all_varname = [[]]*zelen ; all_namefile = [[]]*zelen ; all_time = [[]]*zelen 108 all_var = [[]]*zelen ; all_var2 = [[]]*zelen ; all_title = [[]]*zelen ; all_varname = [[]]*zelen ; all_namefile = [[]]*zelen ; all_time = [[]]*zelen ; all_windu = [[]]*zelen ; all_windv = [[]]*zelen 106 109 107 110 ################################################################################################# … … 191 194 all_time[k] = time 192 195 if var2: all_var2[k] = getfield(nc,var2) 196 if winds: all_windu[k] = getfield(nc,uchar) ; all_windv[k] = getfield(nc,vchar) 193 197 ##### SPECIFIC 194 198 if varname in ["temp","t","T_nadir_nit","T_nadir_day","temp_day","temp_night"] and tsat: … … 219 223 elif ope in ["cat"]: 220 224 tab = all_var[0];k = 1 221 while k != len(namefiles) -1:225 while k != len(namefiles): 222 226 tab = np.append(tab,all_var[k],axis=0) ; k += 1 223 227 all_time[0] = np.arange(0,len(tab),1) ### AS: time reference is too simplistic, should be better … … 252 256 253 257 ### Map projection 254 if mapmode == 1: m = define_proj(proj,wlon,wlat,back=back,blat=blat ) ; x, y = m(lon2d, lat2d)258 if mapmode == 1: m = define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon) ; x, y = m(lon2d, lat2d) 255 259 elif mapmode ==0: m = None ; x = None ; y = None 256 260 … … 264 268 if fileref is not None: index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%(len(namefiles)+2) ## OK only 1 var, see test in the beginning 265 269 elif "var" in ope: index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%(len(var)+1) ## OK only 1 file, see test in the beginning 270 elif "cat" in ope: index_f = 0 266 271 else: yeah = len(namefiles)*len(var) ; index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%yeah 267 272 time = all_time[index_f] … … 284 289 yint=yintegral, alt=vert, anomaly=anomaly ) 285 290 what_I_plot = what_I_plot*mult 286 if var2: ### what is contoured 291 if var2: ### what is contoured. 287 292 what_I_plot_contour, error = reducefield( all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert, \ 288 293 yint=yintegral, alt=vert ) 294 if winds: ### what is plot as vectors. 295 vecx, error = reducefield( all_windu[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert) 296 vecy, error = reducefield( all_windv[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert) 289 297 #################################################################### 290 298 … … 325 333 if mrate is None: what_I_plot_frame = what_I_plot 326 334 else: what_I_plot_frame = what_I_plot[imov,:,:] 335 if winds: 336 if mrate is None: vecx_frame = vecx ; vecy_frame = vecy 337 else: vecx_frame = vecx[imov,:,:] ; vecy_frame = vecy[imov,:,:] 327 338 elif which == "contour": 328 339 if mrate is None or what_I_plot_contour.ndim < 3: what_I_plot_frame = what_I_plot_contour … … 330 341 if mrate is not None: 331 342 if mapmode == 1: 332 m = define_proj(proj,wlon,wlat,back=back,blat=blat ) ## this is dirty, defined above but out of imov loop333 x, y = m(lon2d, lat2d) ## this is dirty, defined above but out of imov loop343 m = define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon) ## this is dirty, defined above but out of imov loop 344 x, y = m(lon2d, lat2d) ## this is dirty, defined above but out of imov loop 334 345 if typefile in ['mesoapi','meso'] and mapmode == 1: what_I_plot_frame = dumpbdy(what_I_plot_frame,6,condition=True) 335 346 if typefile in ['mesoideal']: what_I_plot_frame = dumpbdy(what_I_plot_frame,0,stag='W',condition=dumped_vert_stag) … … 362 373 colorbar( fraction=0.05,pad=0.03,format=daformat,\ 363 374 ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,20])),extend='neither',spacing='proportional' ) 375 if winds: 376 if typefile in ['mesoapi','meso']: 377 [vecx_frame,vecy_frame] = [dumpbdy(vecx_frame,6,stag=uchar,condition=True), dumpbdy(vecy_frame,6,stag=vchar,condition=True)] 378 key = True 379 elif typefile in ['gcm']: 380 key = False 381 if metwind: [vecx_frame,vecy_frame] = m.rotate_vector(vecx_frame, vecy_frame, lon2d, lat2d) 382 if var: colorvec = definecolorvec(back) 383 else: colorvec = definecolorvec(colorb) 384 vectorfield(vecx_frame, vecy_frame, x, y, stride=stride, csmooth=2,\ 385 #scale=15., factor=300., color=colorvec, key=key) 386 scale=20., factor=250., color=colorvec, key=key) 387 #200. ## or csmooth=stride 364 388 elif which == "contour": 365 389 zevminc, zevmaxc = calculate_bounds(what_I_plot_frame) … … 375 399 if which == "regular": 376 400 if mrate is not None: 377 figframe=mpl.pyplot.gcf() 378 if mquality: figframe.set_dpi(600.) 379 else: figframe.set_dpi(200.) 380 mframe=fig2img(figframe) 381 if ((mrate is not None) and (imov is 0)): 382 moviename='test' ;W,H = figframe.canvas.get_width_height() 383 video = VideoSink((H,W), moviename, rate=mrate, byteorder="rgba") 384 video.run(mframe) ; mpl.pyplot.close() 401 ### THIS IS A MENCODER MOVIE 402 if mrate > 0: 403 figframe=mpl.pyplot.gcf() 404 if mquality: figframe.set_dpi(600.) 405 else: figframe.set_dpi(200.) 406 mframe=fig2img(figframe) 407 if imov == 0: 408 moviename='movie' ;W,H = figframe.canvas.get_width_height() 409 video = VideoSink((H,W), moviename, rate=mrate, byteorder="rgba") 410 video.run(mframe) ; close() 411 if imov == iend: video.close() 412 ### THIS IS A WEBPAGE MOVIE 413 else: 414 nameframe = "image"+str(1000+imov) 415 makeplotres(nameframe,res=100.,disp=False) ; close() 416 if imov == 0: myfile = open("zepics", 'w') 417 myfile.write("modImages["+str(imov)+"] = '"+nameframe+"_100.png';"+ '\n') 418 if imov == iend: 419 myfile.write("first_image = 0;"+ '\n') 420 myfile.write("last_image = "+str(iend)+";"+ '\n') 421 myfile.close() 385 422 if var2: which = "contour" 386 423 imov = imov+1 … … 388 425 which = "regular" 389 426 390 if mrate is not None: video.close()391 392 427 ##### 1D field 393 428 elif len(what_I_plot.shape) is 1: … … 395 430 if save == 'txt': writeascii(np.transpose(what_I_plot),'profile'+str(nplot)+'.txt') 396 431 397 #### Other cases: (maybe plot 3-D field one day or movie??)432 #### Other cases: (maybe plot 3-D field one day ??) 398 433 else: 399 434 print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!! dimensions: ",what_I_plot.shape … … 402 437 errormess("There is an error in reducing field !") 403 438 404 ### Vector plot --- a adapter405 if winds:406 vecx, error = reducefield( getfield(nc,uchar), d4=indextime, d3=indexvert , yint=yintegral , alt=vert)407 vecy, error = reducefield( getfield(nc,vchar), d4=indextime, d3=indexvert , yint=yintegral , alt=vert)408 #what_I_plot, error = reducefield(all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert )409 if not error:410 if typefile in ['mesoapi','meso']:411 [vecx,vecy] = [dumpbdy(vecx,6,stag=uchar,condition=True), dumpbdy(vecy,6,stag=vchar,condition=True)]412 key = True413 elif typefile in ['gcm']:414 key = False415 if metwind: [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d)416 if varname == False: colorvec = definecolorvec(back)417 else: colorvec = definecolorvec(colorb)418 vectorfield(vecx, vecy,\419 x, y, stride=stride, csmooth=2,\420 #scale=15., factor=300., color=colorvec, key=key)421 scale=20., factor=250., color=colorvec, key=key)422 #200. ## or csmooth=stride423 424 439 ### Next subplot 425 440 basename = getname(var=varname,winds=winds,anomaly=anomaly) 426 basename = basename + getstralt(nc,level) 441 basename = basename + getstralt(nc,level) 442 if mrate is not None: basename = "movie_" + basename 427 443 if typefile in ['mesoapi','meso']: 428 444 if slon is not None: basename = basename + "_lon_" + str(int(lon2d[indices[1],indices[0]])) -
trunk/UTIL/PYTHON/pp.py
r444 r451 27 27 if opt.fref is not None and opt.operat is not None and opt.itp is not None: interpref=True 28 28 else: interpref=False 29 if opt.rate is not None: opt.save = "avi" 30 elif opt.save == "avi": opt.rate = 8 ## this is a default value for -S avi 31 if opt.save == "html": opt.rate = -1 ## this is convenient because everything is done in planetoplot with mrate 29 32 30 33 ############################# … … 137 140 ope=opt.operat,fileref=reffile,minop=opt.vminope,maxop=opt.vmaxope,titleref=opt.titref,\ 138 141 invert_y=opt.inverty,xaxis=zexaxis,yaxis=zeyaxis,ylog=opt.logy,yintegral=opt.column,\ 139 blat=opt.blat, tsat=opt.tsat,flagnolow=opt.nolow,mrate=opt.rate,mquality=opt.quality)142 blat=opt.blat,blon=opt.blon,tsat=opt.tsat,flagnolow=opt.nolow,mrate=opt.rate,mquality=opt.quality) 140 143 print 'DONE: '+name 141 144 system("rm -f to_be_erased") … … 146 149 for arg in sys.argv: command = command + arg + ' ' 147 150 #if typefile not in ["meso","mesoapi"]: name = 'pycommand' 148 if opt.save is "gui": name = 'pycommand' 151 if opt.save == "gui": name = 'pycommand' 152 elif opt.save == "avi": system("mv -f movie*.avi "+name+".avi") 153 elif opt.save == "html": system("cat $PYTHONPATH/header.html > anim.html ; cat zepics >> anim.html ; cat $PYTHONPATH/body.html >> anim.html ; rm -rf zepics "+name+" ; mkdir "+name+" ; mv anim.html image*png "+name) 149 154 f = open(name+'.sh', 'w') 150 155 f.write(command)
Note: See TracChangeset
for help on using the changeset viewer.
