- Timestamp:
- Mar 1, 2012, 12:50:08 AM (13 years ago)
- Location:
- trunk/UTIL/PYTHON
- Files:
-
- 2 edited
-
myplot.py (modified) (8 diffs)
-
planetoplot.py (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/UTIL/PYTHON/myplot.py
r525 r548 41 41 elif 'phisinit' in nc.variables: typefile = 'gcm' 42 42 elif 'time_counter' in nc.variables: typefile = 'earthgcm' 43 elif hasattr(nc,'START_DATE'): 44 if '9999' in getattr(nc,'START_DATE') : typefile = 'mesoideal' 45 elif 'vert' in nc.variables: typefile = 'mesoapi' 46 elif 'U' in nc.variables: typefile = 'meso' 47 #elif 'HGT_M' in nc.variables: typefile = 'geo' 48 elif 'HGT' in nc.variables: typefile = 'meso' 43 elif hasattr(nc,'START_DATE'): typefile = 'meso' 49 44 elif 'HGT_M' in nc.variables: typefile = 'geo' 50 45 else: typefile = 'gcm' # for lslin-ed files from gcm … … 275 270 ## Author: AS 276 271 def getstralt(nc,nvert): 277 typefile = whatkindfile(nc)278 if typefile is 'meso':272 varinfile = nc.variables.keys() 273 if 'vert' not in varinfile: 279 274 stralt = "_lvl" + str(nvert) 280 el if typefile is 'mesoapi':275 else: 281 276 zelevel = int(nc.variables['vert'][nvert]) 282 277 if abs(zelevel) < 10000.: strheight=str(zelevel)+"m" … … 287 282 elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa" 288 283 else: stralt = "" 289 else:290 stralt = ""291 284 return stralt 292 285 … … 347 340 def getproj (nc): 348 341 typefile = whatkindfile(nc) 349 if typefile in ['meso api','meso','geo']:342 if typefile in ['meso','geo']: 350 343 ### (il faudrait passer CEN_LON dans la projection ?) 351 344 map_proj = getattr(nc, 'MAP_PROJ') … … 366 359 else: 367 360 proj="merc" 368 elif typefile in ['gcm']: proj="cyl" ## pb avec les autres (de trace derriere la sphere ?)369 else: proj="ortho"361 elif typefile in ['gcm']: proj="cyl" ## pb avec les autres (de trace derriere la sphere ?) 362 else: proj="ortho" 370 363 return proj 371 364 … … 439 432 ## getcoord2d for predefined types 440 433 typefile = whatkindfile(nc) 441 if typefile in ['mesoapi','meso']: 442 [lon2d,lat2d] = getcoord2d(nc) 434 if typefile in ['meso']: 435 if '9999' not in getattr(nc,'START_DATE') : 436 ## regular mesoscale 437 [lon2d,lat2d] = getcoord2d(nc) 438 else: 439 ## idealized mesoscale 440 nx=getattr(nc,'WEST-EAST_GRID_DIMENSION') 441 ny=getattr(nc,'SOUTH-NORTH_GRID_DIMENSION') 442 dlat=getattr(nc,'DX') 443 ## this is dirty because Martian-specific 444 # ... but this just intended to get "lat-lon" like info 445 falselon = np.arange(-dlat*(nx-1)/2.,dlat*(nx-1)/2.,dlat)/60000. 446 falselat = np.arange(-dlat*(ny-1)/2.,dlat*(ny-1)/2.,dlat)/60000. 447 [lon2d,lat2d] = np.meshgrid(falselon,falselat) ## dummy coordinates 448 print "WARNING: domain plot artificially centered on lat,lon 0,0" 443 449 elif typefile in ['gcm']: 444 450 [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True) … … 447 453 elif typefile in ['geo']: 448 454 [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M') 449 elif typefile in ['mesoideal']:450 nx=getattr(nc,'WEST-EAST_GRID_DIMENSION')451 ny=getattr(nc,'SOUTH-NORTH_GRID_DIMENSION')452 [lon2d,lat2d] = np.meshgrid(np.arange(nx),np.arange(ny))453 455 return lon2d,lat2d 454 456 … … 539 541 ## Author: AS 540 542 def getwinddef (nc): 541 ## getwinds for predefined types542 typefile = whatkindfile(nc)543 543 ### 544 if typefile is 'mesoapi': [uchar,vchar] = ['Um','Vm'] 545 elif typefile is 'gcm': [uchar,vchar] = ['u','v'] 546 elif typefile in ['meso','mesoideal']: [uchar,vchar] = ['U','V'] 547 else: [uchar,vchar] = ['not found','not found'] 544 varinfile = nc.variables.keys() 545 if 'Um' in varinfile: [uchar,vchar] = ['Um','Vm'] #; print "this is API meso file" 546 elif 'U' in varinfile: [uchar,vchar] = ['U','V'] #; print "this is RAW meso file" 547 elif 'u' in varinfile: [uchar,vchar] = ['u','v'] #; print "this is GCM file" 548 ### you can add choices here ! 549 else: [uchar,vchar] = ['not found','not found'] 548 550 ### 549 if typefile in ['meso']:metwind = False ## geometrical (wrt grid)550 else: 551 if metwind is False: 551 if uchar in ['U']: metwind = False ## geometrical (wrt grid) 552 else: metwind = True ## meteorological (zon/mer) 553 if metwind is False: print "Not using meteorological winds. You trust numerical grid as being (x,y)" 552 554 ### 553 555 return uchar,vchar,metwind -
trunk/UTIL/PYTHON/planetoplot.py
r525 r548 126 126 ### ... TYPEFILE 127 127 typefile = whatkindfile(nc) 128 if typefile in ['mesoideal']: mapmode=0 129 elif typefile in ['gcm'] and len(nc.variables["longitude"][:]) is 1 and len(nc.variables["latitude"][:]) is 1: mapmode=0 130 elif typefile in ['earthgcm'] and len(nc.variables["lon"][:]) is 1 and len(nc.variables["lat"][:]) is 1: mapmode=0 131 if mapmode == 0: winds=False 132 elif mapmode == 1 and redope is None: 128 if typefile in ['gcm'] and len(nc.variables["longitude"][:]) is 1 and len(nc.variables["latitude"][:]) is 1: mapmode=0 ; winds=False 129 elif typefile in ['earthgcm'] and len(nc.variables["lon"][:]) is 1 and len(nc.variables["lat"][:]) is 1: mapmode=0 ; winds=False 130 if redope is None and mapmode == 1: 133 131 if svert is None: svert = readslices(str(level)) ; nvert=1 134 132 if stime is None and mrate is None: … … 152 150 [lon2d,lat2d] = getcoorddef(nc) 153 151 ### ... PROJECTION 154 if ((proj == None) and (typefile not in ['mesoideal'])): proj = getproj(nc)152 if proj == None: proj = getproj(nc) 155 153 156 154 if firstfile: … … 184 182 time = (initime+time*24)%24 185 183 print "LOCAL TIMES.... ", time 186 elif typefile in ['meso',' mesoapi','geo','mesoideal']:184 elif typefile in ['meso','geo']: 187 185 area = None ## not active for the moment 188 186 ###### STUFF TO GET THE CORRECT LAT/LON FROM MESOSCALE FILES WITH 2D LAT/LON ARRAYS … … 208 206 lonp,latp = ( lon2d[indices[0,0,0],indices[0,0,1]] , lat2d[indices[0,0,0],indices[0,0,1]] ) 209 207 ###### 210 if typefile in ['meso api','meso'] and mapmode == 1: lon2d = dumpbdy(lon2d,6) ; lat2d = dumpbdy(lat2d,6) ### important to do that now and not before208 if typefile in ['meso'] and mapmode == 1: lon2d = dumpbdy(lon2d,6) ; lat2d = dumpbdy(lat2d,6) ### important to do that now and not before 211 209 ###### 212 210 if varname in ['PHTOT','W']: vertdim='BOTTOM-TOP_PATCH_END_STAG' … … 333 331 else: indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time) 334 332 ltst = None 335 if typefile in ['meso api','meso'] and indextime is not None: ltst = localtime ( interv[0]+indextime*interv[1], 0.5*(wlon[0]+wlon[1]) )333 if typefile in ['meso'] and indextime is not None: ltst = localtime ( interv[0]+indextime*interv[1], 0.5*(wlon[0]+wlon[1]) ) 336 334 print "********** INDEX LON:",indexlon," LAT:",indexlat," VERT:",indexvert," TIME:",indextime 337 335 ##var = nc.variables["phisinit"][:,:] … … 443 441 m = define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon) ## this is dirty, defined above but out of imov loop 444 442 x, y = m(lon2d, lat2d) ## this is dirty, defined above but out of imov loop 445 if typefile in ['meso api','meso'] and mapmode == 1: what_I_plot_frame = dumpbdy(what_I_plot_frame,6,condition=True)443 if typefile in ['meso'] and mapmode == 1: what_I_plot_frame = dumpbdy(what_I_plot_frame,6,condition=True) 446 444 # if typefile in ['mesoideal']: what_I_plot_frame = dumpbdy(what_I_plot_frame,0,stag='W',condition=dumped_vert_stag) 447 445 … … 471 469 #zelevels = np.linspace(zevmin*(1. + 1.e-7),zevmax*(1. - 1.e-7)) #,num=20) 472 470 zelevels = np.linspace(zevmin,zevmax,num=ticks) 471 #what_I_plot_frame = smooth(what_I_plot_frame,100) 473 472 if mapmode == 1: m.contourf( x, y, what_I_plot_frame, zelevels, cmap = palette, alpha=trans) 474 473 elif mapmode == 0: contourf( x, y, what_I_plot_frame, zelevels, cmap = palette, alpha=trans) … … 487 486 if zeorientation == "horizontal" and zetitle != "fill": zecb.ax.set_xlabel(zetitle) ; zetitle="" 488 487 if winds: 489 if typefile in ['meso api','meso']:488 if typefile in ['meso']: 490 489 [vecx_frame,vecy_frame] = [dumpbdy(vecx_frame,6,stag=uchar,condition=True), dumpbdy(vecy_frame,6,stag=vchar,condition=True)] 491 490 key = True 492 491 elif typefile in ['gcm']: 493 492 key = False 494 if metwind :[vecx_frame,vecy_frame] = m.rotate_vector(vecx_frame, vecy_frame, lon2d, lat2d)493 if metwind and mapmode == 1: [vecx_frame,vecy_frame] = m.rotate_vector(vecx_frame, vecy_frame, lon2d, lat2d) 495 494 if var: colorvec = definecolorvec(back) 496 495 else: colorvec = definecolorvec(colorb) … … 559 558 basename = basename + getstralt(nc,level) 560 559 if mrate is not None: basename = "movie_" + basename 561 if typefile in ['meso api','meso']:560 if typefile in ['meso']: 562 561 if slon is not None: basename = basename + "_lon_" + str(int(round(lonp))) 563 562 if slat is not None: basename = basename + "_lat_" + str(int(round(latp))) … … 598 597 ### Save the figure in a file in the data folder or an user-defined folder 599 598 if outputname is None: 600 if typefile in ['meso' ,'mesoapi']: prefix = getprefix(nc)599 if typefile in ['meso']: prefix = getprefix(nc) 601 600 elif typefile in ['gcm']: prefix = 'LMD_GCM_' 602 601 else: prefix = ''
Note: See TracChangeset
for help on using the changeset viewer.
