Changeset 233 for trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib
- Timestamp:
- Jul 19, 2011, 2:22:19 AM (14 years ago)
- File:
-
- 1 edited
-
trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/myplot.py (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/myplot.py
r232 r233 1 def errormess(text): 2 print text 3 exit() 4 return 5 6 def whatkindfile (nc): 7 if 'controle' in nc.variables: typefile = 'gcm' 8 elif 'vert' in nc.variables: typefile = 'mesoapi' 9 elif 'U' in nc.variables: typefile = 'meso' 10 elif 'HGT_M' in nc.variables: typefile = 'geo' 11 else: errormess("whatkindfile: typefile not supported.") 12 return typefile 13 14 def getfield (nc,var): 15 ## this allows to get much faster (than simply referring to nc.variables[var]) 16 dimension = len(nc.variables[var].dimensions) 17 if dimension == 2: field = nc.variables[var][:,:] 18 elif dimension == 3: field = nc.variables[var][:,:,:] 19 elif dimension == 4: field = nc.variables[var][:,:,:,:] 20 return field 21 22 def reducefield (input,d4=None,d3=None,d2=None,d1=None): 23 ### do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x" 24 ### it would be actually better to name d4 d3 d2 d1 as t z y x 25 import numpy as np 26 dimension = np.array(input).ndim 27 shape = np.array(input).shape 28 print 'dim,shape: ',dimension,shape 29 output = input 30 error = False 31 if dimension == 2: 32 if d2 >= shape[0]: error = True 33 elif d1 >= shape[1]: error = True 34 elif d1 is not None and d2 is not None: output = input[d2,d1] 35 elif d1 is not None: output = input[:,d1] 36 elif d2 is not None: output = input[d2,:] 37 elif dimension == 3: 38 if d4 >= shape[0]: error = True 39 elif d2 >= shape[1]: error = True 40 elif d1 >= shape[2]: error = True 41 elif d4 is not None and d2 is not None and d1 is not None: output = input[d4,d2,d1] 42 elif d4 is not None and d2 is not None: output = input[d4,d2,:] 43 elif d4 is not None and d1 is not None: output = input[d4,:,d1] 44 elif d2 is not None and d1 is not None: output = input[:,d2,d1] 45 elif d1 is not None: output = input[:,:,d1] 46 elif d2 is not None: output = input[:,d2,:] 47 elif d4 is not None: output = input[d4,:,:] 48 elif dimension == 4: 49 if d4 >= shape[0]: error = True 50 elif d3 >= shape[1]: error = True 51 elif d2 >= shape[2]: error = True 52 elif d1 >= shape[3]: error = True 53 elif d4 is not None and d3 is not None and d2 is not None and d1 is not None: output = input[d4,d3,d2,d1] 54 elif d4 is not None and d3 is not None and d2 is not None: output = input[d4,d3,d2,:] 55 elif d4 is not None and d3 is not None and d1 is not None: output = input[d4,d3,:,d1] 56 elif d4 is not None and d2 is not None and d1 is not None: output = input[d4,:,d2,d1] 57 elif d3 is not None and d2 is not None and d1 is not None: output = input[:,d3,d2,d1] 58 elif d4 is not None and d3 is not None: output = input[d4,d3,:,:] 59 elif d4 is not None and d2 is not None: output = input[d4,:,d2,:] 60 elif d4 is not None and d1 is not None: output = input[d4,:,:,d1] 61 elif d3 is not None and d2 is not None: output = input[:,d3,d2,:] 62 elif d3 is not None and d1 is not None: output = input[:,d3,:,d1] 63 elif d2 is not None and d1 is not None: output = input[:,:,d2,d1] 64 elif d1 is not None: output = input[:,:,:,d1] 65 elif d2 is not None: output = input[:,:,d2,:] 66 elif d3 is not None: output = input[:,d3,:,:] 67 elif d4 is not None: output = input[d4,:,:,:] 68 dimension = np.array(output).ndim 69 shape = np.array(output).shape 70 print 'dim,shape: ',dimension,shape 71 return output, error 72 1 73 def latinterv (area): 2 74 if area == "Europe": … … 35 107 return wlon,wlat 36 108 109 def definesubplot ( numplot, fig ): 110 from matplotlib.pyplot import rcParams 111 rcParams['font.size'] = 12. ## default (important for multiple calls) 112 if numplot == 4: 113 sub = 221 114 fig.subplots_adjust(wspace = 0.3, hspace = 0.3) 115 rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. ) 116 elif numplot == 2: 117 sub = 121 118 fig.subplots_adjust(wspace = 0.35) 119 rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. ) 120 elif numplot == 3: 121 sub = 131 122 fig.subplots_adjust(wspace = 0.5) 123 rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) 124 elif numplot == 6: 125 sub = 231 126 fig.subplots_adjust(wspace = 0.4, hspace = 0.0) 127 rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) 128 elif numplot == 8: 129 sub = 331 #241 130 fig.subplots_adjust(wspace = 0.3, hspace = 0.3) 131 rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) 132 elif numplot == 9: 133 sub = 331 134 fig.subplots_adjust(wspace = 0.3, hspace = 0.3) 135 rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) 136 elif numplot == 1: 137 sub = 99999 138 elif numplot < 0: 139 sub = 99999 140 else: 141 print "supported: 1,2,3,4,6,8,9" 142 exit() 143 return sub 144 145 def getstralt(nc,nvert): 146 typefile = whatkindfile(nc) 147 if typefile is 'meso': 148 stralt = "_lvl" + str(nvert) 149 elif typefile is 'mesoapi': 150 zelevel = int(nc.variables['vert'][nvert]) 151 if abs(zelevel) < 10000.: strheight=str(zelevel)+"m" 152 else: strheight=str(int(zelevel/1000.))+"km" 153 if 'altitude' in nc.dimensions: stralt = "_"+strheight+"-AMR" 154 elif 'altitude_abg' in nc.dimensions: stralt = "_"+strheight+"-ALS" 155 elif 'bottom_top' in nc.dimensions: stralt = "_"+strheight 156 elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa" 157 else: stralt = "" 158 else: 159 stralt = "" 160 return stralt 161 37 162 def getlschar ( namefile ): 38 #### strangely enough this does not work for api or ncrcat results!39 163 from netCDF4 import Dataset 40 164 from timestuff import sol2ls 165 from numpy import array 41 166 nc = Dataset(namefile) 42 if namefile[0]+namefile[1]+namefile[2] != "geo" \ 43 and 'Times' in nc.variables \ 167 if 'Times' in nc.variables: 168 zetime = nc.variables['Times'][0] 169 shape = array(nc.variables['Times']).shape 170 if shape[0] < 2: zetime = None 171 if zetime is not None \ 44 172 and 'vert' not in nc.variables: 45 zetime = nc.variables['Times'][0]173 #### strangely enough this does not work for api or ncrcat results! 46 174 zetimestart = getattr(nc, 'START_DATE') 47 175 zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9]) … … 67 195 68 196 def getproj (nc): 69 map_proj = getattr(nc, 'MAP_PROJ') 70 cen_lat = getattr(nc, 'CEN_LAT') 71 if map_proj == 2: 72 if cen_lat > 10.: 73 proj="npstere" 74 print "NP stereographic polar domain" 75 else: 76 proj="spstere" 77 print "SP stereographic polar domain" 78 elif map_proj == 1: 79 print "lambert projection domain" 80 proj="lcc" 81 elif map_proj == 3: 82 print "mercator projection" 83 proj="merc" 84 else: 85 proj="merc" 197 typefile = whatkindfile(nc) 198 if typefile in ['mesoapi','meso','geo']: 199 ### (il faudrait passer CEN_LON dans la projection ?) 200 map_proj = getattr(nc, 'MAP_PROJ') 201 cen_lat = getattr(nc, 'CEN_LAT') 202 if map_proj == 2: 203 if cen_lat > 10.: 204 proj="npstere" 205 print "NP stereographic polar domain" 206 else: 207 proj="spstere" 208 print "SP stereographic polar domain" 209 elif map_proj == 1: 210 print "lambert projection domain" 211 proj="lcc" 212 elif map_proj == 3: 213 print "mercator projection" 214 proj="merc" 215 else: 216 proj="merc" 217 elif typefile in ['gcm']: proj="cyl" ## pb avec les autres (de trace derriere la sphere ?) 218 else: proj="ortho" 86 219 return proj 87 220 … … 102 235 lat1 = lat2d[0,0] 103 236 lat2 = lat2d[nx,ny] 104 wider = 0.5 * (abs(lon1)+abs(lon2)) * 0.1 105 if lon1 < lon2: wlon = [lon1, lon2 + wider] ## a tester en normal 237 if abs(0.5*(lat1+lat2)) > 60.: wider = 0.5 * (abs(lon1)+abs(lon2)) * 0.1 238 else: wider = 0. 239 if lon1 < lon2: wlon = [lon1, lon2 + wider] 106 240 else: wlon = [lon2, lon1 + wider] 107 241 if lat1 < lat2: wlat = [lat1, lat2] … … 118 252 return 119 253 120 def makeplotpng (filename,pad_inches_value=0.25,minres=100.,folder='' ):121 makeplotpngres(filename,minres, pad_inches_value=pad_inches_value,folder=folder )254 def makeplotpng (filename,pad_inches_value=0.25,minres=100.,folder='',disp=True): 255 makeplotpngres(filename,minres, pad_inches_value=pad_inches_value,folder=folder,disp=disp) 122 256 makeplotpngres(filename,minres+200.,pad_inches_value=pad_inches_value,folder=folder,disp=False) 123 257 return 124 258 125 def dumpbdy (field ):259 def dumpbdy (field,stag=None): 126 260 nx = len(field[0,:])-1 127 261 ny = len(field[:,0])-1 262 if stag == 'U': nx = nx-1 263 if stag == 'V': ny = ny-1 128 264 return field[5:ny-5,5:nx-5] 265 266 def getcoorddef ( nc ): 267 ## getcoord2d for predefined types 268 typefile = whatkindfile(nc) 269 if typefile in ['mesoapi','meso']: 270 [lon2d,lat2d] = getcoord2d(nc) 271 lon2d = dumpbdy(lon2d) 272 lat2d = dumpbdy(lat2d) 273 elif typefile in ['gcm']: 274 [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True) 275 elif typefile in ['geo']: 276 [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M') 277 return lon2d,lat2d 129 278 130 279 def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False): … … 168 317 return improc 169 318 170 def getwinds (nc,charu='Um',charv='Vm'): 171 import numpy as np 172 u = nc.variables[charu] 173 v = nc.variables[charv] 174 if charu == 'U': u = u[:, :, :, 0:len(u[0,0,0,:])-1] 175 if charv == 'V': v = v[:, :, 0:len(v[0,0,:,0])-1, :] 176 ### ou alors prendre les coordonnees speciales 177 return u,v 319 def getwinddef (nc): 320 ## getwinds for predefined types 321 typefile = whatkindfile(nc) 322 ### 323 if typefile is 'mesoapi': [uchar,vchar] = ['Um','Vm'] 324 elif typefile is 'gcm': [uchar,vchar] = ['u','v'] 325 elif typefile is 'meso': [uchar,vchar] = ['U','V'] 326 else: [uchar,vchar] = ['not found','not found'] 327 ### 328 if typefile in ['meso']: metwind = False ## geometrical (wrt grid) 329 else: metwind = True ## meteorological (zon/mer) 330 if metwind is False: print "Not using meteorological winds. You trust numerical grid as being (x,y)" 331 ### 332 return uchar,vchar,metwind 178 333 179 334 def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True): … … 216 371 return step 217 372 218 def define_proj (char,wlon,wlat,back= "."):373 def define_proj (char,wlon,wlat,back=None): 219 374 from mpl_toolkits.basemap import Basemap 220 375 import numpy as np … … 247 402 m.drawmeridians(np.r_[-180.:180.:step*2.], labels=[0,0,0,1], color='grey', fontsize=fontsizemer) 248 403 m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer) 249 if back == ".": m.warpimage(marsmap(),scale=0.75) 250 elif back == None: pass 251 else: m.warpimage(marsmap(back),scale=0.75) 404 if back: m.warpimage(marsmap(back),scale=0.75) 405 #if not back: 406 # if not var: back = "mola" ## if no var: draw mola 407 # elif typefile in ['mesoapi','meso','geo'] \ 408 # and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw 409 # else: pass ## else: draw None 252 410 return m 253 411 … … 270 428 return 271 429 430 def calculate_bounds(field,vmin=None,vmax=None): 431 import numpy as np 432 from mymath import max,min,mean 433 ind = np.where(field < 9e+35) 434 fieldcalc = field[ ind ] # la syntaxe compacte ne marche si field est un tuple 435 ### 436 dev = np.std(fieldcalc)*3.0 437 ### 438 if vmin is None: 439 zevmin = mean(fieldcalc) - dev 440 else: zevmin = vmin 441 ### 442 if vmax is None: zevmax = mean(fieldcalc) + dev 443 else: zevmax = vmax 444 if vmin == vmax: 445 zevmin = mean(fieldcalc) - dev ### for continuity 446 zevmax = mean(fieldcalc) + dev ### for continuity 447 ### 448 if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0. 449 print "field ", min(fieldcalc), max(fieldcalc) 450 print "bounds ", zevmin, zevmax 451 return zevmin, zevmax 452 453 def bounds(what_I_plot,zevmin,zevmax): 454 ### might be convenient to add the missing value in arguments 455 what_I_plot[ what_I_plot < zevmin ] = zevmin*(1. + 1.e-7) 456 what_I_plot[ what_I_plot > 9e+35 ] = -9e+35 457 what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7) 458 return what_I_plot 459 460 def zoomset (wlon,wlat,zoom): 461 dlon = abs(wlon[1]-wlon[0])/2. 462 dlat = abs(wlat[1]-wlat[0])/2. 463 [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\ 464 [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ] 465 print "zoom %",zoom,wlon,wlat 466 return wlon,wlat 272 467 273 468 def fmtvar (whichvar="def"): … … 286 481 return fmtvar[whichvar] 287 482 483 #################################################################################################################### 484 ### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png 288 485 def defcolorb (whichone="def"): 289 486 whichcolorb = { \ … … 319 516 320 517 def marsmap (whichone="vishires"): 518 from os import uname 519 mymachine = uname()[1] 520 ### not sure about speed-up with this method... looks the same 521 if "lmd.jussieu.fr" in mymachine: domain = "/u/aslmd/WWW/maps/" 522 else: domain = "http://www.lmd.jussieu.fr/~aslmd/maps/" 321 523 whichlink = { \ 322 "vis": "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\ 323 "vishires": "http://dl.dropbox.com/u/11078310/MarsMap_2500x1250.jpg",\ 324 "geolocal": "http://dl.dropbox.com/u/11078310/geolocal.jpg",\ 325 "mola": "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\ 326 "molabw": "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\ 524 #"vis": "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\ 525 #"vishires": "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\ 526 #"geolocal": "http://dl.dropbox.com/u/11078310/geolocal.jpg",\ 527 #"mola": "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\ 528 #"molabw": "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\ 529 "vis": domain+"mar0kuu2.jpg",\ 530 "vishires": domain+"MarsMap_2500x1250.jpg",\ 531 "geolocal": domain+"geolocal.jpg",\ 532 "mola": domain+"mars-mola-2k.jpg",\ 533 "molabw": domain+"MarsElevation_2500x1250.jpg",\ 327 534 } 328 535 if whichone not in whichlink:
Note: See TracChangeset
for help on using the changeset viewer.
