Changeset 225 for trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py
- Timestamp:
- Jul 15, 2011, 2:47:21 PM (13 years ago)
- File:
-
- 1 edited
-
trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py
r207 r225 14 14 target=None, 15 15 stride=3,\ 16 numplot= 4,\16 numplot=2,\ 17 17 var=None,\ 18 18 colorb=True,\ … … 52 52 elif 'vert' in nc.variables: typefile = 'mesoapi' 53 53 elif 'U' in nc.variables: typefile = 'meso' 54 elif 'HGT_M' in nc.variables: typefile = 'geo' 54 55 else: 55 56 print "typefile not supported." … … 59 60 ############################################################## 60 61 ### Try to guess the projection from wrfout if not set by user 61 if typefile in ['mesoapi','meso' ]:62 if typefile in ['mesoapi','meso','geo']: 62 63 if proj == None: proj = getproj(nc) 63 64 ### (il faudrait passer CEN_LON dans la projection ?) … … 70 71 if not back: 71 72 if not var: back = "mola" ## if no var: draw mola 72 elif typefile in ['mesoapi','meso' ] \73 and proj not in ['merc','lcc','nsper','laea']:back = "molabw" ## if var but meso: draw molabw73 elif typefile in ['mesoapi','meso','geo'] \ 74 and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw 74 75 else: pass ## else: draw None 75 76 … … 82 83 elif typefile in ['gcm']: 83 84 [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True) 85 elif typefile in ['geo']: 86 [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M') 84 87 if proj == "npstere": [wlon,wlat] = latinterv("North_Pole") 85 88 elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d) 86 89 else: [wlon,wlat] = simplinterv(lon2d,lat2d) 90 print wlon 91 print wlat 87 92 if zoom: 88 93 dlon = abs(wlon[1]-wlon[0])/2. … … 105 110 metwind = False ## geometrical (wrt grid) 106 111 print "Beware ! Not using meteorological winds. You trust numerical grid as being (x,y)." 112 elif typefile is 'geo': 113 winds = None 107 114 108 115 ##################################################### … … 122 129 if vmax is None: zevmax = mean(field) + dev 123 130 else: zevmax = vmax 131 if vmin == vmax: 132 zevmin = mean(field) - dev ### for continuity 133 zevmax = mean(field) + dev ### for continuity 134 print "field ", min(field), max(field) 124 135 print "bounds ", zevmin, zevmax 125 136 ### some already defined colormaps … … 137 148 ######################################### 138 149 ### Name for title and graphics save file 139 if winds: basename = "UV_" 140 else: basename = "" 141 if var: basename = basename + var 142 ### 143 if dimension == 4: 150 if var and winds: basename = var + '_UV' 151 elif var: basename = var 152 elif winds: basename = 'UV' 153 else: exit() 154 ### 155 if dimension == 4 or winds: 144 156 if typefile is 'meso': stralt = "_lvl" + str(nvert) 145 157 elif typefile is 'mesoapi': … … 158 170 ### Open a figure and set subplots 159 171 fig = figure() 172 rcParams['font.size'] = 12. 173 if typefile in ['geo']: numplot = 1 160 174 if numplot > 0: 161 175 if numplot == 4: … … 228 242 if var: 229 243 if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(field[i,:,:]) 244 elif typefile in ['geo']: what_I_plot = field[0,:,:] 230 245 elif typefile in ['gcm']: 231 246 if dimension == 2: what_I_plot = field[:,:] … … 233 248 palette = get_cmap(name=colorb) 234 249 #palette.set_over('b', 1.0) 250 #print np.array(x).shape 251 #print np.array(y).shape 252 #print np.array(what_I_plot).shape 235 253 if not tile: 236 254 zelevels = np.linspace(zevmin,zevmax) 255 what_I_plot[ what_I_plot < zevmin ] = zevmin*(1. + 1.e-7) 256 what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7) 237 257 contourf( x, y, what_I_plot, 10, cmap = palette, levels = zelevels ) 238 258 else: 239 259 pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax ) 240 if var in ['HGT']: 260 if var in ['HGT']: pass 241 261 elif colorb: 242 243 244 245 246 262 ndiv = 10 263 colorbar(fraction=0.05,pad=0.1,format=fmtvar(var),\ 264 ticks=np.linspace(zevmin,zevmax,ndiv+1),\ 265 extend='max',spacing='proportional') 266 # both min max neither 247 267 ### Vector plot 248 268 if winds: … … 257 277 else: colorvec = definecolorvec(colorb) 258 278 vectorfield(vecx, vecy,\ 259 x, y, stride=stride, csmooth= stride,\279 x, y, stride=stride, csmooth=2,\ 260 280 scale=15., factor=300., color=colorvec, key=key) 261 #200. ## or csmooth= 2281 #200. ## or csmooth=stride 262 282 263 283 ### Next subplot 264 284 plottitle = basename 265 if addchar: plottitle = plottitle + addchar + "_LT"+str(ltst) 266 else: plottitle = plottitle + "_LT"+str(ltst) 285 if typefile in ['mesoapi','meso']: 286 if addchar: plottitle = plottitle + addchar + "_LT"+str(ltst) 287 else: plottitle = plottitle + "_LT"+str(ltst) 267 288 ptitle( plottitle ) 268 289 sub += 1 … … 297 318 from netCDF4 import Dataset 298 319 from myplot import getlschar 320 from numpy import ones 299 321 300 322 ############################# 301 323 ### Get options and variables 302 324 parser = OptionParser() 303 parser.add_option('-f', action=' store', dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file')325 parser.add_option('-f', action='append', dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file (append)') 304 326 parser.add_option('-l', action='store', dest='nvert', type="float", default=0, help='vertical level (def=0)(-i 2: p,mbar)(-i 3,4: z,km)') 305 327 parser.add_option('-p', action='store', dest='proj', type="string", default=None, help='projection') … … 307 329 parser.add_option('-t', action='store', dest='target', type="string", default=None, help='destination folder') 308 330 parser.add_option('-s', action='store', dest='stride', type="int", default=3, help='stride vectors (def=3)') 309 parser.add_option('-v', action=' store', dest='var', type="string", default=None, help='variable contoured')310 parser.add_option('-n', action='store', dest='numplot', type="int", default= 4, help='number of plots (def=1)(<0: 1 plot of LT -*numplot*)')331 parser.add_option('-v', action='append', dest='var', type="string", default=None, help='variable contoured (append)') 332 parser.add_option('-n', action='store', dest='numplot', type="int", default=2, help='number of plots (def=1)(<0: 1 plot of LT -*numplot*)') 311 333 parser.add_option('-i', action='store', dest='interp', type="int", default=None, help='interpolation method (2: press, 3: z-amr, 4:z-als)') 312 334 parser.add_option('-c', action='store', dest='colorb', type="string", default=True, help='change colormap') 313 335 parser.add_option('-x', action='store_false', dest='winds', default=True, help='no wind vectors') 314 parser.add_option('-m', action=' store', dest='vmin', type="float", default=None, help='bounding minimum value for color plot')315 parser.add_option('-M', action=' store', dest='vmax', type="float", default=None, help='bounding maximum value for color plot')336 parser.add_option('-m', action='append', dest='vmin', type="float", default=None, help='bounding minimum value (append)') 337 parser.add_option('-M', action='append', dest='vmax', type="float", default=None, help='bounding maximum value (append)') 316 338 parser.add_option('-T', action='store_true', dest='tile', default=False, help='draw a tiled plot (no blank zone)') 317 339 parser.add_option('-z', action='store', dest='zoom', type="float", default=None, help='zoom factor in %') … … 323 345 exit() 324 346 print "Options:", opt 325 zefile = opt.namefile 326 zelevel = opt.nvert 327 stralt = None 328 [lschar,zehour,zehourin] = getlschar ( zefile ) ## getlschar from wrfout (or simply return "" if another file) 329 330 ##################################################### 331 ### Call Fortran routines for vertical interpolations 332 if opt.interp is not None: 333 if opt.nvert is 0 and opt.interp is 4: zelevel = 0.010 334 ### winds or no winds 335 if opt.winds : zefields = 'uvmet' 336 else : zefields = '' 337 ### var or no var 338 if opt.var is None : pass 339 elif zefields == '' : zefields = opt.var 340 else : zefields = zefields + "," + opt.var 341 print zefields 342 zefile = api_onelevel ( path_to_input = '', \ 343 input_name = zefile, \ 344 path_to_output = opt.target, \ 345 fields = zefields, \ 346 interp_method = opt.interp, \ 347 onelevel = zelevel, \ 348 nocall = opt.nocall ) 349 print zefile 350 zelevel = 0 ## so that zelevel could play again the role of nvert 351 352 ############# 353 ### Main call 354 name = winds (zefile,int(zelevel),\ 355 proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\ 356 addchar=lschar,interv=[zehour,zehourin],vmin=opt.vmin,vmax=opt.vmax,tile=opt.tile,zoom=opt.zoom) 357 print 'Done: '+name 358 359 ######################################################### 360 ### Generate a .sh file with the used command saved in it 361 command = "" 362 for arg in sys.argv: command = command + arg + ' ' 363 f = open(name+'.sh', 'w') 364 f.write(command) 347 348 if opt.var is None: 349 pass 350 else: 351 listvar = '' 352 zelen = len(opt.var) 353 if zelen == 1: listvar = opt.var[0] + ',' 354 else : 355 for jjj in range(zelen): listvar += opt.var[jjj] + ',' 356 listvar = listvar[0:len(listvar)-1] 357 if opt.vmin: 358 if zelen != len(opt.vmin): 359 print "not enough or too much vmin values... setting same values all variables" 360 vmintab = ones(zelen) * opt.vmin[0] 361 else: 362 vmintab = opt.vmin 363 else: 364 vmintab = None 365 if opt.vmax: 366 if zelen != len(opt.vmax): 367 print "not enough or too much vmax values... setting same values all variables" 368 vmaxtab = ones(zelen) * opt.vmax[0] 369 else: 370 vmaxtab = opt.vmax 371 else: 372 vmaxtab = None 373 374 for i in range(len(opt.namefile)): 375 376 zefile = opt.namefile[i] 377 print zefile 378 zelevel = opt.nvert 379 stralt = None 380 [lschar,zehour,zehourin] = getlschar ( zefile ) ## getlschar from wrfout (or simply return "" if another file) 381 382 ##################################################### 383 ### Call Fortran routines for vertical interpolations 384 if opt.interp is not None: 385 if opt.nvert is 0 and opt.interp is 4: zelevel = 0.010 386 ### winds or no winds 387 if opt.winds : zefields = 'uvmet' 388 else : zefields = '' 389 ### var or no var 390 if opt.var is None : pass 391 elif zefields == '' : zefields = listvar 392 else : zefields = zefields + "," + listvar 393 print zefields 394 zefile = api_onelevel ( path_to_input = '', \ 395 input_name = zefile, \ 396 #path_to_output = opt.target, \ 397 fields = zefields, \ 398 interp_method = opt.interp, \ 399 onelevel = zelevel, \ 400 nocall = opt.nocall ) 401 print zefile 402 zelevel = 0 ## so that zelevel could play again the role of nvert 403 404 if opt.var is None: 405 ############# 406 ### Main call 407 name = winds (zefile,int(zelevel),\ 408 proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\ 409 addchar=lschar,interv=[zehour,zehourin],vmin=opt.vmin,vmax=opt.vmax,tile=opt.tile,zoom=opt.zoom) 410 print 'Done: '+name 411 else: 412 for jjj in range(len(opt.var)): 413 if vmintab: argvmin = vmintab[jjj] 414 else: argvmin = None 415 if vmaxtab: argvmax = vmaxtab[jjj] 416 else: argvmax = None 417 ############# 418 ### Main call 419 name = winds (zefile,int(zelevel),\ 420 proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var[jjj],numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\ 421 addchar=lschar,interv=[zehour,zehourin],vmin=argvmin,vmax=argvmax,tile=opt.tile,zoom=opt.zoom) 422 print 'Done: '+name 423 424 ######################################################### 425 ### Generate a .sh file with the used command saved in it 426 command = "" 427 for arg in sys.argv: command = command + arg + ' ' 428 f = open(name+'.sh', 'w') 429 f.write(command)
Note: See TracChangeset
for help on using the changeset viewer.
