Changeset 876
- Timestamp:
- Feb 11, 2013, 2:22:32 PM (12 years ago)
- Location:
- trunk/UTIL/PYTHON
- Files:
-
- 4 edited
-
myplot.py (modified) (42 diffs)
-
myscript.py (modified) (1 diff)
-
planetoplot.py (modified) (48 diffs)
-
pp.py (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/UTIL/PYTHON/myplot.py
r874 r876 1 2 import numpy as np 3 import netCDF4 4 1 5 ## Author: AS 2 6 def errormess(text,printvar=None): … … 8 12 ## Author: AS 9 13 def adjust_length (tab, zelen): 10 from numpy import ones11 14 if tab is None: 12 outtab = ones(zelen) * -99999915 outtab = np.ones(zelen) * -999999 13 16 else: 14 17 if zelen != len(tab): 15 18 print "not enough or too much values... setting same values all variables" 16 outtab = ones(zelen) * tab[0]19 outtab = np.ones(zelen) * tab[0] 17 20 else: 18 21 outtab = tab … … 31 34 ## Author: AS + AC 32 35 def localtime(time,lon,namefile): # lon is the mean longitude of the plot, not of the domain. central lon of domain is taken from cen_lon 33 import numpy as np34 from netCDF4 import Dataset35 36 ## THIS IS FOR MESOSCALE 36 nc = Dataset(namefile)37 nc = netCDF4.Dataset(namefile) 37 38 ## get start date and intervals 38 39 dt_hour=1. ; start=0. … … 93 94 94 95 ## Author: AS 96 def getfieldred (nc,var,indexlon,indexlat,indexvert,indextime): 97 dimension = len(nc.variables[var].dimensions) 98 ## this allows to get much faster and use much less memory esp. with large datasets 99 if dimension == 2: field = nc.variables[var][indextime,indexlon] 100 elif dimension == 3: field = nc.variables[var][indextime,indexlat,indexlon] 101 elif dimension == 4: field = nc.variables[var][indextime,indexvert,indexlat,indexlon] 102 elif dimension == 1: field = nc.variables[var][indextime] 103 return field 104 105 ## Author: AS 95 106 def getfield (nc,var): 107 dimension = len(nc.variables[var].dimensions) 96 108 ## this allows to get much faster (than simply referring to nc.variables[var]) 97 import numpy as np 98 dimension = len(nc.variables[var].dimensions) 99 #print " Opening variable with", dimension, "dimensions ..." 109 print " Opening variable",var," with", dimension, "dimensions ..." 100 110 if dimension == 2: field = nc.variables[var][:,:] 101 111 elif dimension == 3: field = nc.variables[var][:,:,:] … … 142 152 # The corresponding variable to call is UV or uvmet (to use api) 143 153 def windamplitude (nc,mode): 144 import numpy as np145 154 varinfile = nc.variables.keys() 146 155 if "U" in varinfile: zu=getfield(nc,'U') … … 177 186 # this only requires co2col so that you can concat.nc at low cost 178 187 def enrichment_factor(nc,lon,lat,time): 179 import numpy as np180 188 from myplot import reducefield 181 189 varinfile = nc.variables.keys() … … 213 221 # The corresponding variable to call is SLOPEXY 214 222 def slopeamplitude (nc): 215 import numpy as np216 223 varinfile = nc.variables.keys() 217 224 if "slopex" in varinfile: zu=getfield(nc,'slopex') … … 233 240 # The corresponding variable to call is DELTAT 234 241 def deltat0t1 (nc): 235 import numpy as np236 242 varinfile = nc.variables.keys() 237 243 if "tsurf" in varinfile: zu=getfield(nc,'tsurf') … … 251 257 ### it would be actually better to name d4 d3 d2 d1 as t z y x 252 258 ### ... note, anomaly is only computed over d1 and d2 for the moment 253 import numpy as np254 259 from mymath import max,mean,min,sum,getmask 255 260 csmooth = 12 ## a fair amount of grid points (too high results in high computation time) … … 444 449 from mymath import max,mean 445 450 from scipy import integrate 446 import numpy as np447 451 if yint and vert is not None and indice is not None: 448 452 if type(input).__name__=='MaskedArray': … … 483 487 rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. ) 484 488 elif numplot <= 6: 485 subv = 2 486 subh = 3 487 #fig.subplots_adjust(wspace = 0.4, hspace = 0.0) 488 fig.subplots_adjust(wspace = 0.5, hspace = 0.3) 489 subv = 3#2 490 subh = 2#3 491 ##fig.subplots_adjust(wspace = 0.4, hspace = 0.0) 492 #fig.subplots_adjust(wspace = 0.5, hspace = 0.3) 493 fig.subplots_adjust(wspace = 0.3, hspace = 0.5) 489 494 rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) 490 495 elif numplot <= 8: … … 531 536 ## Author: AS 532 537 def getlschar ( namefile, getaxis=False ): 533 from netCDF4 import Dataset534 538 from timestuff import sol2ls 535 from numpy import array536 539 from string import rstrip 537 540 import os as daos … … 542 545 namefile = namefiletest 543 546 #### we assume that wrfout is next to wrfout_z and wrfout_zabg 544 nc = Dataset(namefile)547 nc = netCDF4.Dataset(namefile) 545 548 zetime = None 546 549 days_in_month = [61, 66, 66, 65, 60, 54, 50, 46, 47, 47, 51, 56] … … 548 551 if 'Times' in nc.variables: 549 552 zetime = nc.variables['Times'][0] 550 shape = array(nc.variables['Times']).shape553 shape = np.array(nc.variables['Times']).shape 551 554 if shape[0] < 2: zetime = None 552 555 if zetime is not None \ … … 621 624 ## Author: AS 622 625 def polarinterv (lon2d,lat2d): 623 import numpy as np624 626 wlon = [np.min(lon2d),np.max(lon2d)] 625 627 ind = np.array(lat2d).shape[0] / 2 ## to get a good boundlat and to get the pole … … 629 631 ## Author: AS 630 632 def simplinterv (lon2d,lat2d): 631 import numpy as np632 633 return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]] 633 634 … … 679 680 ## Author: AS + AC 680 681 def getcoorddef ( nc ): 681 import numpy as np682 682 ## getcoord2d for predefined types 683 683 typefile = whatkindfile(nc) … … 712 712 ## Author: AS 713 713 def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False): 714 import numpy as np715 714 if nlon == "nothing" or nlat == "nothing": 716 715 print "NO LAT LON FIELDS. I AM TRYING MY BEST. I ASSUME GLOBAL FIELD." … … 740 739 ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth 741 740 def smooth1d(x,window_len=11,window='hanning'): 742 import numpy743 741 """smooth the data using a window with requested size. 744 742 This method is based on the convolution of a scaled window with the signal. … … 762 760 TODO: the window parameter could be the window itself if an array instead of a string 763 761 """ 764 x = n umpy.array(x)762 x = np.array(x) 765 763 if x.ndim != 1: 766 764 raise ValueError, "smooth only accepts 1 dimension arrays." … … 771 769 if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']: 772 770 raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'" 773 s=n umpy.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]771 s=np.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]] 774 772 #print(len(s)) 775 773 if window == 'flat': #moving average 776 w=n umpy.ones(window_len,'d')774 w=np.ones(window_len,'d') 777 775 else: 778 776 w=eval('numpy.'+window+'(window_len)') 779 y=n umpy.convolve(w/w.sum(),s,mode='valid')777 y=np.convolve(w/w.sum(),s,mode='valid') 780 778 return y 781 779 … … 789 787 ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth 790 788 def gauss_kern(size, sizey=None): 791 import numpy as np792 789 # Returns a normalized 2D gauss kernel array for convolutions 793 790 size = int(size) … … 832 829 ## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs. 833 830 import matplotlib.pyplot as plt 834 import numpy as np835 831 #posx = np.min(x) - np.std(x) / 10. 836 832 #posy = np.min(y) - np.std(y) / 10. … … 874 870 def define_proj (char,wlon,wlat,back=None,blat=None,blon=None): 875 871 from mpl_toolkits.basemap import Basemap 876 import numpy as np877 872 import matplotlib as mpl 878 873 from mymath import max … … 978 973 ## Author: AS 979 974 def calculate_bounds(field,vmin=None,vmax=None): 980 import numpy as np981 975 from mymath import max,min,mean 982 976 ind = np.where(field < 9e+35) … … 1022 1016 ## Author : AC 1023 1017 def hole_bounds(what_I_plot,zevmin,zevmax): 1024 import numpy as np1025 1018 zi=0 1026 1019 for i in what_I_plot: … … 1264 1257 ## Author: TN 1265 1258 def separatenames (name): 1266 from numpy import concatenate1267 1259 # look for comas in the input name to separate different names (files, variables,etc ..) 1268 1260 if name is None: … … 1279 1271 else: 1280 1272 name1 = currentname[0:indexvir] 1281 names = concatenate((names,[name1]))1273 names = np.concatenate((names,[name1])) 1282 1274 currentname = currentname[indexvir+1:len(currentname)] 1283 1275 return names … … 1286 1278 ## Author: TN 1287 1279 def readslices(saxis): 1288 from numpy import empty1289 1280 if saxis == None: 1290 1281 zesaxis = None 1291 1282 else: 1292 zesaxis = empty((len(saxis),2))1283 zesaxis = np.empty((len(saxis),2)) 1293 1284 for i in range(len(saxis)): 1294 1285 a = separatenames(saxis[i]) … … 1304 1295 def readdata(data,datatype,coord1,coord2): 1305 1296 ## Read sparse data 1306 from numpy import empty1307 1297 if datatype == 'txt': 1308 1298 if len(data[coord1].shape) == 1: … … 1320 1310 ## Author: AS 1321 1311 def bidimfind(lon2d,lat2d,vlon,vlat,file=None): 1322 import numpy as np1323 1312 import matplotlib.pyplot as mpl 1324 1313 if vlat is None: array = (lon2d - vlon)**2 … … 1342 1331 # input : all the desired slices and the good index 1343 1332 # output : all indexes to be taken into account for reducing field 1344 import numpy as np1345 1333 if ( np.array(axis).ndim == 2): 1346 1334 axis = axis[:,0] … … 1364 1352 # x axis priority: 1/time 2/lon 3/lat 4/vertical 1365 1353 # To be improved !!!... 1366 from numpy import array,swapaxes1367 1354 x = None 1368 1355 y = None 1369 1356 count = 0 1370 what_I_plot = array(what_I_plot)1357 what_I_plot = np.array(what_I_plot) 1371 1358 shape = what_I_plot.shape 1372 1359 if indextime is None and len(time) > 1: … … 1394 1381 else: y = vert 1395 1382 count = count+1 1396 x = array(x)1397 y = array(y)1383 x = np.array(x) 1384 y = np.array(y) 1398 1385 print "CHECK SHAPE: what_I_plot, x, y", what_I_plot.shape, x.shape, y.shape 1399 1386 if len(shape) == 1: … … 1402 1389 elif len(shape) == 2: 1403 1390 if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]: 1404 print "INFO: swapaxes: ",what_I_plot.shape,shape ; what_I_plot = swapaxes(what_I_plot,0,1)1391 print "INFO: swapaxes: ",what_I_plot.shape,shape ; what_I_plot = np.swapaxes(what_I_plot,0,1) 1405 1392 else: 1406 1393 if shape[0] != len(y): print "WARNING: shape[0] != len(y). Correcting." ; what_I_plot = what_I_plot[0:len(y),:] … … 1445 1432 vecx=None,vecy=None,stride=2 ): 1446 1433 ### an easy way to map a field over lat/lon grid 1447 import numpy as np1448 1434 import matplotlib.pyplot as mpl 1449 1435 from matplotlib.cm import get_cmap … … 1496 1482 specificname_gcm = ['enfact'] 1497 1483 1484 zncvar = znc.variables 1485 1498 1486 ## Check for variable in file: 1499 1487 if mode == 'check': 1500 1488 varname = zvarname 1501 varinfile=znc .variables.keys()1502 logical_novarname = zvarname not in znc .variables1489 varinfile=zncvar.keys() 1490 logical_novarname = zvarname not in zncvar 1503 1491 logical_nospecificname_meso = not ((ztypefile in ['meso']) and (zvarname in specificname_meso)) 1504 1492 logical_nospecificname_gcm = not ((ztypefile in ['gcm']) and (zvarname in specificname_gcm)) … … 1519 1507 all_var=get_tsat(yalt,tinput,zlon=ylon,zlat=ylat,zalt=yalt,ztime=ytime) 1520 1508 ### ----------- 2. wind amplitude 1521 elif ((zvarname in ['UV','uv','uvmet']) and (ztypefile in ['meso']) and (zvarname not in znc .variables)):1509 elif ((zvarname in ['UV','uv','uvmet']) and (ztypefile in ['meso']) and (zvarname not in zncvar)): 1522 1510 all_var=windamplitude(znc,'amplitude') 1523 elif ((zvarname in ['hodograph','hodograph_2']) and (ztypefile in ['meso']) and (zvarname not in znc .variables)):1511 elif ((zvarname in ['hodograph','hodograph_2']) and (ztypefile in ['meso']) and (zvarname not in zncvar)): 1524 1512 plot_x, plot_y = windamplitude(znc,zvarname) 1525 1513 if plot_x is not None: all_var=plot_x # dummy 1526 1514 else: all_var=plot_y ; plot_x = None ; plot_y = None # Hodograph type 2 is not 'xy' mode 1527 elif ((zvarname in ['slopexy','SLOPEXY']) and (ztypefile in ['meso']) and (zvarname not in znc .variables)):1515 elif ((zvarname in ['slopexy','SLOPEXY']) and (ztypefile in ['meso']) and (zvarname not in zncvar)): 1528 1516 all_var=slopeamplitude(znc) 1529 1517 ### ------------ 3. Near surface instability 1530 elif ((zvarname in ['DELTAT','deltat']) and (ztypefile in ['meso']) and (zvarname not in znc .variables)):1518 elif ((zvarname in ['DELTAT','deltat']) and (ztypefile in ['meso']) and (zvarname not in zncvar)): 1531 1519 all_var=deltat0t1(znc) 1532 1520 ### ------------ 4. Enrichment factor … … 1534 1522 all_var=enrichment_factor(znc,ylon,ylat,ytime) 1535 1523 ### ------------ 5. teta -> temp 1536 elif ((ztypefile in ['meso']) and (zvarname in ['tk']) and ('tk' not in znc .variables.keys())):1524 elif ((ztypefile in ['meso']) and (zvarname in ['tk']) and ('tk' not in zncvar.keys())): 1537 1525 all_var=teta_to_tk(znc) 1538 1526 else: … … 1549 1537 # (which is not efficient but it is still ok) and then, make the average (if the user wants to) 1550 1538 def spectrum(var,time,vert,lat,lon): 1551 import numpy as np1552 1539 fft=np.fft.fft(var,axis=1) 1553 1540 N=len(vert) … … 1565 1552 # Computes temperature from potential temperature for mesoscale files, without the need to use API, i.e. using natural vertical grid 1566 1553 def teta_to_tk(nc): 1567 import numpy as np1568 1554 varinfile = nc.variables.keys() 1569 1555 p0=610. … … 1588 1574 # 6/ in this slab, find the point at which the surface pressure is minimum 1589 1575 def find_devil(nc,indextime): 1590 import numpy as np1591 1576 from scipy import ndimage 1592 1577 from mymath import array2image,image2array -
trunk/UTIL/PYTHON/myscript.py
r865 r876 3 3 ### I/O 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 OK: use -f "foo*" DONT FORGET QUOTES "" !!!!') 5 parser.add_option('-L', '--large', action='store_true',dest='monster', default=False, help='speedy version for large files (EXPERIMENTAL)') 5 6 parser.add_option('--seevar', action='store_true',dest='seevar', default=False, help='display the list of variables in the file') 6 7 parser.add_option('-t', '--target', action='store',dest='tgt', type="string", default=None, help='destination folder') -
trunk/UTIL/PYTHON/planetoplot.py
r874 r876 78 78 streamflag=False,\ 79 79 nocolorb=False,\ 80 analysis=None): 80 analysis=None,\ 81 monster=False): 81 82 82 83 #################################################################################################################### … … 85 86 ################################# 86 87 ### Load librairies and functions 87 from netCDF4 import Dataset 88 from myplot import getcoord2d,define_proj,makeplotres,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\ 89 fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\ 90 zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\ 91 getname,localtime,check_localtime,polarinterv,getsindex,define_axis,determineplot,readslices,bidimfind,getlschar,hole_bounds,\ 92 getdimfromvar,select_getfield,find_devil 93 from mymath import deg,max,min,mean,writeascii,fig2data,fig2img 88 import netCDF4 89 94 90 import matplotlib as mpl 95 from matplotlib.pyplot import contour,contourf,hist, text,subplot, figure, rcParams, savefig, colorbar, \ 96 pcolor, show, plot, clabel, title, close, legend, xlabel, axis, ylabel, \ 97 subplots_adjust, axes, clabel 98 from matplotlib.cm import get_cmap 91 import matplotlib.pyplot 92 import matplotlib.cm 99 93 from mpl_toolkits.basemap import cm 94 if mrate is not None: from videosink import VideoSink 95 100 96 import numpy as np 101 97 from numpy.core.defchararray import find 102 from scipy.stats import gaussian_kde,describe 103 from videosink import VideoSink 98 import scipy 99 if analysis in ['laplace']: import scipy.ndimage.laplace as laplace 100 101 import itertools 102 import os 103 104 import myplot 105 from mymath import deg,max,min,mean,writeascii,fig2data,fig2img 104 106 from times import sol2ls 105 import subprocess106 107 #from singlet import singlet 107 from itertools import cycle108 import os109 from scipy import ndimage110 108 111 109 … … 117 115 print "********** WELCOME TO PLANETOPLOT **********" 118 116 print "********************************************" 117 if monster: 118 print "********* SPEED MODE. EXPERIMENTAL. ********" 119 print "********************************************" 119 120 ### we ensure namefiles and var are arrays 120 121 if not isinstance(namefiles, np.ndarray): namefiles = [namefiles] … … 132 133 if trycol: clb = ["Greys","Blues","YlOrRd","jet","spectral","hot","RdBu","RdYlBu","Paired"] ; zetitle = clb ; var = [var[0]]*9 133 134 ### we avoid specific cases not yet implemented 134 if mrate is not None and len(var) > 1: errormess("multivar not allowed in movies. should be fixed soon!")135 if mrate is not None and len(var) > 1: myplot.errormess("multivar not allowed in movies. should be fixed soon!") 135 136 ### we prepare number of plot fields [zelen] and number of plot instances [numplot] according to user choices 136 137 ### --> we support multifile and multivar plots : files + vars separated by commas are on the same figure 137 nlon, nlat, nvert, ntime, mapmode, nslices = determineplot(slon, slat, svert, stime, redope)138 nlon, nlat, nvert, ntime, mapmode, nslices = myplot.determineplot(slon, slat, svert, stime, redope) 138 139 zelen = len(namefiles)*len(var) 140 if (nslices > 1 and monster): myplot.errormess("multislice + monster not supported yet. to be done soon") 139 141 ### we have a special mode obtained by -p noproj in which lat/lon plots are just flat 2D plots 140 142 if proj == "noproj": mapmode = 0 … … 157 159 ### LOOP OVER PLOT FIELDS ### 158 160 ############################# 159 161 160 162 for nnn in range(len(namefiles)): 161 163 for vvv in range(len(var)): … … 163 165 ### we load each NETCDF objects in namefiles 164 166 namefile = namefiles[nnn] 165 nc = Dataset(namefile)167 nc = netCDF4.Dataset(namefile) 166 168 ### we explore the variables in the file 167 169 varinfile = nc.variables.keys() 168 170 if seevar: print varinfile ; exit() 169 171 ### we define the type of file we have (gcm, meso, etc...) 170 typefile = whatkindfile(nc)172 typefile = myplot.whatkindfile(nc) 171 173 if firstfile: typefile0 = typefile 172 elif typefile != typefile0: errormess("Not the same kind of files !", [typefile0, typefile])174 elif typefile != typefile0: myplot.errormess("Not the same kind of files !", [typefile0, typefile]) 173 175 ### we care for input file being 1D simulations 174 176 is1d=999 175 if "longitude" in nc.dimensions and "latitude" in nc.dimensions: is1d = len(nc. variables["longitude"][:])*len(nc.variables["latitude"][:])176 elif "lon" in nc.dimensions and "lat" in nc.dimensions: is1d = len(nc. variables["lon"][:])*len(nc.variables["lat"][:])177 if "longitude" in nc.dimensions and "latitude" in nc.dimensions: is1d = len(nc.dimensions["longitude"])*len(nc.dimensions["latitude"]) 178 elif "lon" in nc.dimensions and "lat" in nc.dimensions: is1d = len(nc.dimensions["lon"])*len(nc.dimensions["lat"]) 177 179 if typefile in ['gcm','earthgcm'] and is1d == 1: mapmode=0 ; winds=False 178 180 ### we create default vert and time prescriptions if not here in case mapping mode is on (lat/lon maps) 179 181 if redope is None and mapmode == 1: 180 if svert is None: svert = readslices(str(level)) ; nvert=1182 if svert is None: svert = myplot.readslices(str(level)) ; nvert=1 181 183 if stime is None and mrate is None: 182 stime = readslices(str(0)) ; ntime=1 ## this is a default choice184 stime = myplot.readslices(str(0)) ; ntime=1 ## this is a default choice 183 185 print "WELL... nothing about time axis. I took default: first time reference stored in file." 184 186 ### we get the names of variables to be read. in case only one available, we choose this one. 185 187 ### (we take out of the test specific names e.g. UV is not in the file but used to ask a wind speed computation) 186 varname = select_getfield(zvarname=var[vvv],znc=nc,ztypefile=typefile,mode='check')188 varname = myplot.select_getfield(zvarname=var[vvv],znc=nc,ztypefile=typefile,mode='check') 187 189 ### we get the names of wind variables to be read (if any) 188 190 if winds: 189 [uchar,vchar,metwind] = getwinddef(nc)191 [uchar,vchar,metwind] = myplot.getwinddef(nc) 190 192 if uchar == 'not found': winds = False 191 193 ### we tell the user that either no var or no wind is not acceptable 192 if not varname and not winds: errormess("please set at least winds or var",printvar=nc.variables)194 if not varname and not winds: myplot.errormess("please set at least winds or var",printvar=nc.variables) 193 195 ### we get the coordinates lat/lon to be used 194 [lon2d,lat2d] = getcoorddef(nc)196 [lon2d,lat2d] = myplot.getcoorddef(nc) 195 197 ### we get an adapted map projection if none is provided by the user 196 if proj == None: proj = getproj(nc)198 if proj == None: proj = myplot.getproj(nc) 197 199 ### we define plot boundaries given projection or user choices 198 200 if firstfile: 199 if proj in ["npstere","spstere"]: [wlon,wlat] = polarinterv(lon2d,lat2d)200 elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d)201 else: [wlon,wlat] = simplinterv(lon2d,lat2d)202 if zoom: [wlon,wlat] = zoomset(wlon,wlat,zoom)203 elif zarea is not None: [wlon,wlat] = latinterv(area=zarea)201 if proj in ["npstere","spstere"]: [wlon,wlat] = myplot.polarinterv(lon2d,lat2d) 202 elif proj in ["lcc","laea"]: [wlon,wlat] = myplot.wrfinterv(lon2d,lat2d) 203 else: [wlon,wlat] = myplot.simplinterv(lon2d,lat2d) 204 if zoom: [wlon,wlat] = myplot.zoomset(wlon,wlat,zoom) 205 elif zarea is not None: [wlon,wlat] = myplot.latinterv(area=zarea) 204 206 205 207 ############################################################# … … 228 230 ### --> add a line here if your reference is not present 229 231 else: 230 dadim = getdimfromvar(nc) ; print "No altitude found. Try to build a simple axis.",dadim232 dadim = myplot.getdimfromvar(nc) ; print "No altitude found. Try to build a simple axis.",dadim 231 233 if len(dadim) == 4: print "-- 4D field. Assume z is dim 2." ; vert = np.arange(dadim[-3]) 232 234 elif len(dadim) == 3: print "-- 3D field. Assume z is dim 1." ; vert = [0.] 233 235 else: vert = [0.] 234 236 ### we define time vector. either it is referenced or it is guessed based on last variable's dimensions. 235 if "Time" in nc.variables: time = nc.variables["Time"][:]236 elif "time_counter" in nc.variables: time = nc.variables["time_counter"][:]/86400. #### convert from s to days237 elif "time" in nc.variables: time = nc.variables["time"][:]237 if "Time" in nc.variables: letime = "Time" 238 elif "time_counter" in nc.variables: letime = "time_counter" 239 elif "time" in nc.variables: letime = "time" 238 240 ### --> add a line here if your reference is not present 241 else: letime = None 242 if letime is not None: 243 if monster: 244 ### very long simulation... we just retrieve 3 values for time 245 timelength = len(nc.dimensions[letime]) 246 dafirst = nc.variables[letime][0] 247 dalast = nc.variables[letime][timelength-1] 248 step = nc.variables[letime][1] - dafirst 249 time = np.arange(dafirst,dalast,step) 250 else: 251 ### if the time record is not too long, what follows is pretty quick 252 time = nc.variables[letime][:] 239 253 else: 240 254 print "No time found. Try to build a simple axis. Assume t is dim 1." 241 dadim = getdimfromvar(nc)255 dadim = myplot.getdimfromvar(nc) 242 256 if len(dadim) == 4: time = np.arange(dadim[-4]) 243 257 elif len(dadim) == 3: time = np.arange(dadim[-3]) 244 else: time = [0.] #errormess("no time axis found.") 258 else: time = [0.] #myplot.errormess("no time axis found.") 259 ### (SPECIFIC?) 260 if "time_counter" in nc.variables: time = time / 86400. #### convert from s to days 245 261 ### (SPECIFIC. convert to Ls for Martian GCM files.) 246 262 if axtime in ["ls"]: … … 256 272 ### (simply ask for subscript) 257 273 if axtime in ["ind"]: 258 dadim = getdimfromvar(nc)274 dadim = myplot.getdimfromvar(nc) 259 275 if len(dadim) == 4: time = np.arange(dadim[-4]) 260 276 elif len(dadim) == 3: time = np.arange(dadim[-3]) … … 280 296 if slon is not None: vlon = slon[0][iii] ### note: slon[:][0] does not work 281 297 if slat is not None: vlat = slat[0][jjj] ### note: slon[:][0] does not work 282 indices[iii,jjj,:] = bidimfind(lon2d,lat2d,vlon,vlat,file=iwantawhereplot)298 indices[iii,jjj,:] = myplot.bidimfind(lon2d,lat2d,vlon,vlat,file=iwantawhereplot) 283 299 lonp,latp = ( lon2d[indices[iii,jjj,0],indices[iii,jjj,1]] , lat2d[indices[iii,jjj,0],indices[iii,jjj,1]] ) 284 300 ### possible bug here if several --lat … … 290 306 ### we get rid of boundary relaxation zone for plots. important to do that now and not before. 291 307 if (typefile in ['meso'] and mapmode == 1): 292 if '9999' not in getattr(nc,'START_DATE'): lon2d = dumpbdy(lon2d,6) ; lat2d =dumpbdy(lat2d,6)308 if '9999' not in getattr(nc,'START_DATE'): lon2d = myplot.dumpbdy(lon2d,6) ; lat2d = myplot.dumpbdy(lat2d,6) 293 309 ### we read the keyword for vertical dimension. we take care for vertical staggering. 294 310 if varname in ['PHTOT','W']: vertdim='BOTTOM-TOP_PATCH_END_STAG' … … 304 320 ### we define the time axis and take care of various specificities (lt, ls, sol) or issues (concatenation) 305 321 if axtime in ["ls","sol"]: 306 lstab, soltab, lttab = getlschar ( namefile, getaxis = True )322 lstab, soltab, lttab = myplot.getlschar ( namefile, getaxis = True ) 307 323 if axtime == "ls": time = lstab 308 324 elif axtime == "sol": time = soltab … … 315 331 if axtime in ["lt"]: 316 332 ftime = np.zeros(len(time)) 317 for i in range(len(time)): ftime[i] = localtime ( time[i], slon , namefile )318 time=ftime ; time =check_localtime(time)333 for i in range(len(time)): ftime[i] = myplot.localtime ( time[i], slon , namefile ) 334 time=ftime ; time = myplot.check_localtime(time) 319 335 print "LOCAL TIMES.... ", time 320 336 ### we define the vertical axis (or lack thereof) and cover possibilities for it to be altitude, pressure, geopotential. quite SPECIFIC. 321 if typefile in ['geo']: vert = [0.] ; stime = readslices(str(0))337 if typefile in ['geo']: vert = [0.] ; stime = myplot.readslices(str(0)) 322 338 else: 323 339 if vertmode is None: vertmode=0 … … 332 348 #################################################################### 333 349 350 334 351 ### we fill the arrays of varname, namefile, time, colorbar at the current step considered (NB: why use both k and nnn ?) 335 352 all_varname[k] = varname 336 353 all_namefile[k] = namefile 337 all_time[k] = time338 all_vert[k] = vert339 all_lat[k] = lat340 all_lon[k] = lon341 354 all_colorb[k] = clb[vvv] 342 if var2: all_var2[k], plot_x[k], plot_y[k] = select_getfield(zvarname=var2,znc=nc,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis) 343 if winds: all_windu[k] = getfield(nc,uchar) ; all_windv[k] = getfield(nc,vchar) 344 ### we fill the arrays of fields to be plotted at the current step considered 345 all_var[k], plot_x[k], plot_y[k] = select_getfield(zvarname=all_varname[k],znc=nc,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis) 355 356 ############################################################################## 357 ##### LARGE FILE. WE'D BETTER NOT FILL THE MEMORY WITH THE STUPID WHOLE ARRAY! 358 if monster: 359 #################################################################### 360 ## get all indexes to be taken into account for this subplot and then reduce field 361 ## We plot 1) all lon slices 2) all lat slices 3) all vert slices 4) all time slices and then go to the next slice 362 if ope is not None: 363 if fileref is not None: index_f = (k//(nlon*nlat*nvert*ntime))%(3*len(namefiles)) ## OK only 1 var, see test in the beginning 364 elif "var" in ope: index_f = (k//(nlon*nlat*nvert*ntime))%(len(var)+1) ## OK only 1 file, see test in the beginning 365 elif "cat" in ope: index_f = 0 366 elif not (k == 0): 367 #if len(namefiles) > 1 and len(var) == 1 and which == "unidim": pass ## TROUVER UN MOYEN POUR unidim 368 if len(namefiles) > 1 and len(var) == 1: pass 369 else: yeah = len(namefiles)*len(var) ; index_f = (k//(nlon*nlat*nvert*ntime))%yeah 370 else: yeah = len(namefiles)*len(var) ; index_f = (k//(nlon*nlat*nvert*ntime))%yeah 371 372 ilon = myplot.getsindex(sslon,k%nlon,lon) 373 ilat = myplot.getsindex(sslat,(k//nlon)%nlat,lat) 374 ivert = myplot.getsindex(svert,(k//(nlon*nlat))%nvert,vert) 375 376 if mrate is not None: itime = None 377 else: itime = myplot.getsindex(stime,(k//(nlon*nlat*nvert))%ntime,time) 378 ltst = None 379 if typefile in ['meso'] and itime is not None and len(itime) < 2: ltst = myplot.localtime ( itime, slon , all_namefile[index_f] ) 380 381 if ilat is not None: print "********** LAT : INDEX",ilat[0], ilat[-1], "VALUE",lat[ilat[0]], lat[ilat[-1]] 382 else: ilat = range(len(lat)) 383 if ilon is not None: print "********** LON : INDEX",ilon[0], ilon[-1], "VALUE",lon[ilon[0]], lon[ilon[-1]] 384 else: ilon = range(len(lon)) 385 if ivert is not None: print "********** VERT: INDEX",ivert[0], ivert[-1], "VALUE",vert[ivert[0]], vert[ivert[-1]] 386 else: ivert = range(len(vert)) 387 if itime is not None: print "********** TIME: INDEX",itime[0], itime[-1], "VALUE",time[itime[0]], time[itime[-1]] 388 else: itime = range(len(time)) 389 390 all_time[k] = time[itime] 391 all_vert[k] = vert[ivert] 392 all_lat[k] = lat[ilat] 393 all_lon[k] = lon[ilon] 394 395 all_var[k] = myplot.getfieldred(nc,all_varname[k],ilon,ilat,ivert,itime) 396 if var2: all_var2[k] = myplot.getfieldred(nc,var2,ilon,ilat,ivert,itime) 397 if winds: 398 all_windu[k] = myplot.getfield(nc,uchar,ilon,ilat,ivert,itime) 399 all_windv[k] = myplot.getfield(nc,vchar,ilon,ilat,ivert,itime) 400 plot_x[k] = None ; plot_y[k] = None 401 402 else: 403 ## regular stuff: not large array. 404 all_time[k] = time 405 all_vert[k] = vert 406 all_lat[k] = lat 407 all_lon[k] = lon 408 if var2: all_var2[k], plot_x[k], plot_y[k] = myplot.select_getfield(zvarname=var2,znc=nc,ztypefile=typefile,\ 409 mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis) 410 if winds: all_windu[k] = myplot.getfield(nc,uchar) ; all_windv[k] = myplot.getfield(nc,vchar) 411 ### we fill the arrays of fields to be plotted at the current step considered 412 all_var[k], plot_x[k], plot_y[k] = myplot.select_getfield(zvarname=all_varname[k],znc=nc,ztypefile=typefile,\ 413 mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis) 414 #################################################################### 346 415 347 416 # last line of for namefile in namefiles … … 354 423 if ope is not None and "var" in ope: 355 424 print "********** OPERATION: ",ope 356 if len(namefiles) > 1: errormess("for this operation... please set only one file !")357 if len(var) > 2: errormess("not sure this works for more than 2 vars... please check.")425 if len(namefiles) > 1: myplot.errormess("for this operation... please set only one file !") 426 if len(var) > 2: myplot.errormess("not sure this works for more than 2 vars... please check.") 358 427 if "div_var" in ope: all_var[k] = all_var[k-2] / all_var[k-1] ; insert = '_div_' 359 428 elif "mul_var" in ope: all_var[k] = all_var[k-2] * all_var[k-1] ; insert = '_mul_' 360 429 elif "add_var" in ope: all_var[k] = all_var[k-2] + all_var[k-1] ; insert = '_add_' 361 430 elif "sub_var" in ope: all_var[k] = all_var[k-2] - all_var[k-1] ; insert = '_sub_' 362 else: errormess(ope+" : non-implemented operation. Check pp.py --help")431 else: myplot.errormess(ope+" : non-implemented operation. Check pp.py --help") 363 432 plot_x[k] = None ; plot_y[k] = None 364 433 all_time[k] = all_time[k-1] ; all_vert[k] = all_vert[k-1] ; all_lat[k] = all_lat[k-1] ; all_lon[k] = all_lon[k-1] ; all_namefile[k] = all_namefile[k-1] … … 383 452 print "********** OPERATION: ",ope 384 453 for k in np.arange(zelen): 385 if len(var) > 1: errormess("for this operation... please set only one var !")454 if len(var) > 1: myplot.errormess("for this operation... please set only one var !") 386 455 if ope in ["-","+","-%","-_only","+_only","-%_only","-_histo"]: 387 if fileref is None: errormess("fileref is missing!")456 if fileref is None: myplot.errormess("fileref is missing!") 388 457 else:ncref = Dataset(fileref) 389 458 … … 399 468 print "SETTING REFERENCE PLOT" 400 469 all_varname[k] = all_varname[k-1] ; all_time[k] = all_time[k-1] ; all_vert[k] = all_vert[k-1] ; all_lat[k] = all_lat[k-1] ; all_lon[k] = all_lon[k-1] ; all_namefile[k] = all_namefile[k-1] ; all_var2[k] = all_var2[k-1] ; all_colorb[k] = all_colorb[k-1] 401 all_var[k], plot_x[k], plot_y[k] = select_getfield(zvarname=all_varname[k-1],znc=ncref,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis)402 if winds: all_windu[k] = getfield(ncref,uchar) ; all_windv[k] =getfield(ncref,vchar)470 all_var[k], plot_x[k], plot_y[k] = myplot.select_getfield(zvarname=all_varname[k-1],znc=ncref,ztypefile=typefile,mode='getvar',ztsat=tsat,ylon=all_lon[k],ylat=all_lat[k],yalt=all_vert[k],ytime=all_time[k],analysis=analysis) 471 if winds: all_windu[k] = myplot.getfield(ncref,uchar) ; all_windv[k] = myplot.getfield(ncref,vchar) 403 472 404 473 if (k+1)%3==0: ## operation plots … … 434 503 all_time[0] = np.array(tabtime) ; all_var[0] = np.array(tab) ; numplot = 1 435 504 if var2: all_var2[0] = np.array(tab2) 436 else: errormess(ope+" : non-implemented operation. Check pp.py --help")505 else: myplot.errormess(ope+" : non-implemented operation. Check pp.py --help") 437 506 if "only" in ope: 438 507 numplot = 1 ; all_var[0] = all_var[k] … … 442 511 ################################## 443 512 ### Open a figure and set subplots 444 fig = figure()445 subv,subh = definesubplot( numplot, fig, ipreferline = (mapmode == 1) )513 fig = mpl.pyplot.figure() 514 subv,subh = myplot.definesubplot( numplot, fig, ipreferline = (mapmode == 1) ) 446 515 if ope in ['-','-%','-_histo'] and len(namefiles) ==1 : subv,subh = 2,2 447 516 elif ope in ['-','-%'] and len(namefiles)>1 : subv, subh = len(namefiles), 3 … … 450 519 ### Time loop for plotting device 451 520 nplot = 1 ; error = False ; firstpass = True 452 if lstyle is not None: linecycler = cycle(lstyle)453 else: linecycler = cycle(["-","--","-.",":"])521 if lstyle is not None: linecycler = itertools.cycle(lstyle) 522 else: linecycler = itertools.cycle(["-","--","-.",":"]) 454 523 print "********************************************" 455 524 while error is False: … … 470 539 else: yeah = len(namefiles)*len(var) ; index_f = ((nplot-1)//(nlon*nlat*nvert*ntime))%yeah 471 540 time = all_time[index_f] ; vert = all_vert[index_f] ; lat = all_lat[index_f] ; lon = all_lon[index_f] 472 indexlon = getsindex(sslon,(nplot-1)%nlon,lon)473 indexlat = getsindex(sslat,((nplot-1)//nlon)%nlat,lat)474 indexvert = getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert)541 indexlon = myplot.getsindex(sslon,(nplot-1)%nlon,lon) 542 indexlat = myplot.getsindex(sslat,((nplot-1)//nlon)%nlat,lat) 543 indexvert = myplot.getsindex(svert,((nplot-1)//(nlon*nlat))%nvert,vert) 475 544 plotx=plot_x[index_f] ; ploty=plot_y[index_f] 476 545 if mrate is not None: indextime = None 477 else: indextime = getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time)546 else: indextime = myplot.getsindex(stime,((nplot-1)//(nlon*nlat*nvert))%ntime,time) 478 547 ltst = None 479 if typefile in ['meso'] and indextime is not None and len(indextime) < 2: ltst = localtime ( indextime, slon , all_namefile[index_f]) 480 481 482 #if mapmode == 0: 483 # if xlab is None: 484 # if indexlon is None: xlab = "Longitude" 485 # elif indexlat is None: xlab = "Latitude" 486 # if ylab is None: 487 # if indexvert is None: ylab = "Altitude" 488 489 490 print "********** INDEX LON:",indexlon," LAT:",indexlat," VERT:",indexvert," TIME:",indextime 548 if typefile in ['meso'] and indextime is not None and len(indextime) < 2: ltst = myplot.localtime ( indextime, slon , all_namefile[index_f] ) 549 550 if not monster: 551 if indexlat is not None: print "********** LAT : INDEX",indexlat[0], indexlat[-1], "VALUE",lat[indexlat[0]], lat[indexlat[-1]] 552 if indexlon is not None: print "********** LON : INDEX",indexlon[0], indexlon[-1], "VALUE",lon[indexlon[0]], lon[indexlon[-1]] 553 if indexvert is not None: print "********** VERT: INDEX",indexvert[0], indexvert[-1], "VALUE",vert[indexvert[0]], vert[indexvert[-1]] 554 if indextime is not None: print "********** TIME: INDEX",indextime[0], indextime[-1], "VALUE",time[indextime[0]], time[indextime[-1]] 555 491 556 ##var = nc.variables["phisinit"][:,:] 492 ##contourf(np.transpose(var),30,cmap = get_cmap(name="Greys_r") ) ; axis('off') ; plot(indexlat,indexlon,'mx',mew=4.0,ms=20.0)493 ## show()557 ##contourf(np.transpose(var),30,cmap = get_cmap(name="Greys_r") ) ; mpl.pyplot.axis('off') ; plot(indexlat,indexlon,'mx',mew=4.0,ms=20.0) 558 ##mpl.pyplot.show() 494 559 ##exit() 495 560 #truc = True … … 503 568 which='' 504 569 if varname: ### what is shaded. 505 what_I_plot, error = reducefield( all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \570 what_I_plot, error = myplot.reducefield( all_var[index_f], d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \ 506 571 yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area, unidim=is1d) 507 572 if add != 0.: what_I_plot = what_I_plot + add … … 510 575 511 576 if var2: ### what is contoured. 512 what_I_plot_contour, error = reducefield( all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert, \577 what_I_plot_contour, error = myplot.reducefield( all_var2[index_f], d4=indextime, d1=indexlon, d2=indexlat , d3=indexvert, \ 513 578 yint=yintegral, alt=vert, redope=redope ) 514 579 if winds: ### what is plot as vectors. 515 vecx, error = reducefield( all_windu[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert)516 vecy, error = reducefield( all_windv[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert)580 vecx, error = myplot.reducefield( all_windu[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert) 581 vecy, error = myplot.reducefield( all_windv[index_f], d4=indextime, d3=indexvert, yint=yintegral, alt=vert) 517 582 if varname in [uchar,vchar]: what_I_plot = np.sqrt( np.square(vecx) + np.square(vecy) ) ; varname = "wind" 518 583 519 584 if plotx is not None: 520 plotx, error = reducefield( plotx, d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \585 plotx, error = myplot.reducefield( plotx, d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \ 521 586 yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area, unidim=is1d) 522 ploty, error = reducefield( ploty, d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \587 ploty, error = myplot.reducefield( ploty, d4=indextime, d1=indexlon, d2=indexlat, d3=indexvert, \ 523 588 yint=yintegral, alt=vert, anomaly=anomaly, redope=redope, mesharea=area, unidim=is1d) 524 589 which='xy' … … 526 591 #if truc: 527 592 # nx = what_I_plot.shape[2] ; ny = what_I_plot.shape[1] ; nz = what_I_plot.shape[0] 528 # for k in range(nz): print k,' over ',nz ; what_I_plot[k,:,:] = what_I_plot[k,:,:] / smooth(what_I_plot[k,:,:],12)593 # for k in range(nz): print k,' over ',nz ; what_I_plot[k,:,:] = what_I_plot[k,:,:] / myplot.smooth(what_I_plot[k,:,:],12) 529 594 # for iii in range(nx): 530 595 # for jjj in range(ny): … … 534 599 # if np.abs(what_I_plot[0,jjj,iii]) > 1.5: 535 600 # print iii,jjj,what_I_plot[0,jjj,iii],int(abs(1.-mx)*100.),int(abs(1.-mn)*100.) 536 # plot(rel)537 # show()601 # mpl.pyplot.plot(rel) 602 # mpl.pyplot.show() 538 603 # anomaly = True ### pour avoir les bons reglages plots 539 604 # what_I_plot = what_I_plot[0,:,:] … … 543 608 ### General plot settings 544 609 changesubplot = (numplot > 1) and (len(what_I_plot.shape) != 1) and (which != "xy") ## default for 1D plots: superimposed. to be reworked for better flexibility. 545 if changesubplot: subplot(subv,subh,nplot) #;subplots_adjust(wspace=0,hspace=0)610 if changesubplot: mpl.pyplot.subplot(subv,subh,nplot) #; mpl.pyplot.subplots_adjust(wspace=0,hspace=0) 546 611 colorb = all_colorb[index_f] 547 612 #################################################################### 548 613 if error: 549 errormess("There is an error in reducing field !")614 myplot.errormess("There is an error in reducing field !") 550 615 else: 551 616 ticks = ndiv + 1 … … 565 630 if slon is not None or proj == "noproj": latyeah = lat2d[:,milieux] 566 631 if slat is not None or proj == "noproj": lonyeah = lon2d[milieuy,:] 567 what_I_plot, x, y = define_axis(lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\632 what_I_plot, x, y = myplot.define_axis(lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\ 568 633 itime,what_I_plot, len(all_var[index_f].shape),vertmode,redope) 569 634 ### 570 if analysis in ['laplace']: what_I_plot = ndimage.laplace(what_I_plot)635 if analysis in ['laplace']: what_I_plot = laplace(what_I_plot) 571 636 ### 572 if (fileref is not None) and ((index_f+1)%3 == 0): zevmin, zevmax = calculate_bounds(what_I_plot,vmin=minop,vmax=maxop)573 else: zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax)637 if (fileref is not None) and ((index_f+1)%3 == 0): zevmin, zevmax = myplot.calculate_bounds(what_I_plot,vmin=minop,vmax=maxop) 638 else: zevmin, zevmax = myplot.calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax) 574 639 #if (fileref is not None) and (index_f == numplot-1): colorb = "RdBu_r" 575 if colorb in ["def","nobar","onebar"]: palette = get_cmap(name=defcolorb(fvar.upper()))640 if colorb in ["def","nobar","onebar"]: palette = mpl.cm.get_cmap(name=myplot.defcolorb(fvar.upper())) 576 641 elif colorb == "relief": palette = cm.GMT_relief 577 642 elif colorb == "haxby": palette = cm.GMT_haxby 578 else: palette = get_cmap(name=colorb)643 else: palette = mpl.cm.get_cmap(name=colorb) 579 644 #palette = cm.GMT_split 580 645 #palette = cm.GMT_globe … … 585 650 elif len(what_I_plot.shape) >= 4: 586 651 print "WARNING!!! ",len(what_I_plot.shape),"-D PLOT NOT SUPPORTED !!! dimensions: ",what_I_plot.shape 587 errormess("Are you sure you did not forget to prescribe a dimension ?")652 myplot.errormess("Are you sure you did not forget to prescribe a dimension ?") 588 653 ##### 2. HANDLE simple 1D/2D field and movies of 1D/2D fields 589 654 else: … … 595 660 if var2 and which == '': which = "contour" ## have to start with contours rather than shading 596 661 elif which == '': which = "regular" 597 if mrate is None: errormess("3D field. Use --rate RATE for movie or specify --time TIME. Exit.")662 if mrate is None: myplot.errormess("3D field. Use --rate RATE for movie or specify --time TIME. Exit.") 598 663 elif len(what_I_plot.shape) == 2: 599 664 if var2 and which == '': which = "contour" ## have to start with contours rather than shading … … 621 686 #if mrate is not None: 622 687 if mapmode == 1: 623 m = define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon) ## this is dirty, defined above but out of imov loop688 m = myplot.define_proj(proj,wlon,wlat,back=back,blat=blat,blon=blon) ## this is dirty, defined above but out of imov loop 624 689 x, y = m(lon2d, lat2d) ## this is dirty, defined above but out of imov loop 625 690 if (typefile in ['meso'] and mapmode == 1): 626 if '9999' not in getattr(nc,'START_DATE'): what_I_plot_frame = dumpbdy(what_I_plot_frame,6,condition=True)627 # if typefile in ['mesoideal']: what_I_plot_frame = dumpbdy(what_I_plot_frame,0,stag='W',condition=dumped_vert_stag)691 if '9999' not in getattr(nc,'START_DATE'): what_I_plot_frame = myplot.dumpbdy(what_I_plot_frame,6,condition=True) 692 # if typefile in ['mesoideal']: what_I_plot_frame = myplot.dumpbdy(what_I_plot_frame,0,stag='W',condition=dumped_vert_stag) 628 693 629 694 if which == "unidim": … … 645 710 else: zemarker = None 646 711 this_is_a_regular_plot = (indexvert is not None) or (indextime is None) or (indexlat is None) or (indexlon is None) 647 if this_is_a_regular_plot: plot(x,what_I_plot_frame,zeline,label=lbl,marker=zemarker) ## vertical profile648 else: plot(what_I_plot_frame,x,zeline,label=lbl,marker=zemarker) ## regular plot712 if this_is_a_regular_plot: mpl.pyplot.plot(x,what_I_plot_frame,zeline,label=lbl,marker=zemarker) ## vertical profile 713 else: mpl.pyplot.plot(what_I_plot_frame,x,zeline,label=lbl,marker=zemarker) ## regular plot 649 714 mpl.pyplot.grid(True) 650 if nplot > 1: legend(loc='best')651 if indextime is None and axtime is not None and xlab is None: xlabel(axtime.upper()) ## define the right label715 if nplot > 1: mpl.pyplot.legend(loc='best') 716 if indextime is None and axtime is not None and xlab is None: mpl.pyplot.xlabel(axtime.upper()) ## define the right label 652 717 if save == 'txt': writeascii(np.transpose([x,np.transpose(what_I_plot)]),'profile'+str(nplot*1000+imov)+'.txt') 653 718 if axtime == "lt" and indextime is None: … … 682 747 else: 683 748 plotx = np.linspace(min(ploty.flatten()),max(ploty.flatten()),1000) 684 density = gaussian_kde(ploty.flatten())749 density = scipy.stats.gaussian_kde(ploty.flatten()) 685 750 # density.covariance_factor = lambda : .25 # adjust the covariance factor to change the bandwidth if needed 686 751 # density._compute_covariance() 687 752 # display the mean and variance of the kde: 688 753 sample = density.resample(size=20000) 689 n, (smin, smax), sm, sv, ss, sk = describe(sample[0])754 n, (smin, smax), sm, sv, ss, sk = scipy.stats.describe(sample[0]) 690 755 mpl.pyplot.plot(plotx,density(plotx), label = lbl+'\nmean: '+str(sm)[0:5]+' std: '+str(np.sqrt(sv))[0:5]+'\nskewness: '+str(ss)[0:5]+' kurtosis: '+str(sk)[0:5]) 691 756 if analysis == 'histodensity': # plot both histo and density (to assess the rightness of the kernel density estimate for exemple) and display the estimated variance … … 693 758 if index_f == zelen-1: mpl.pyplot.legend() ; mpl.pyplot.title(zetitle) 694 759 else: 695 plot(plotx,ploty,label = lbl)760 mpl.pyplot.plot(plotx,ploty,label = lbl) 696 761 if index_f == zelen-1: mpl.pyplot.legend() ; mpl.pyplot.title(zetitle) 697 762 mpl.pyplot.grid(True) 698 763 else: 699 plot(plotx,ploty,label = lbl)764 mpl.pyplot.plot(plotx,ploty,label = lbl) 700 765 if index_f == zelen-1: mpl.pyplot.legend() ; mpl.pyplot.title(zetitle) 701 766 if varname == 'hodograph': … … 715 780 print 'INFO: Using stream file',streamfile, 'for stream lines' 716 781 ncstream = Dataset(streamfile) 717 psi = getfield(ncstream,'psi')782 psi = myplot.getfield(ncstream,'psi') 718 783 psi = psi[0,:,:,0] # time and longitude are dummy dimensions 719 784 if psi.shape[1] != len(x) or psi.shape[0] != len(y): 720 errormess('stream file does not fit! Dimensions: '+str(psi.shape)+' '+str(x.shape)+' '+str(y.shape))785 myplot.errormess('stream file does not fit! Dimensions: '+str(psi.shape)+' '+str(x.shape)+' '+str(y.shape)) 721 786 zelevels = np.arange(-1.e10,1.e10,1.e9) 722 787 zemin = np.min(abs(zelevels)) 723 788 zemax = np.max(abs(zelevels)) 724 789 zewidth = (abs(zelevels)-zemin)*(5.- 0.5)/(zemax - zemin) + 0.5 # linewidth ranges from 5 to 0.5 725 cs = contour( x,y,psi, zelevels, colors='k', linewidths = zewidth)726 clabel(cs, inline=True, fontsize = 4.*rcParams['font.size']/5., fmt="%1.1e")790 cs = mpl.pyplot.contour( x,y,psi, zelevels, colors='k', linewidths = zewidth) 791 mpl.pyplot.clabel(cs, inline=True, fontsize = 4.*mpl.pyplot.rcParams['font.size']/5., fmt="%1.1e") 727 792 else: 728 793 print 'WARNING: STREAM FILE',streamfile, 'DOES NOT EXIST !' 729 794 730 if hole: what_I_plot_frame = hole_bounds(what_I_plot_frame,zevmin,zevmax)731 else: what_I_plot_frame = bounds(what_I_plot_frame,zevmin,zevmax)732 if flagnolow: what_I_plot_frame = nolow(what_I_plot_frame)795 if hole: what_I_plot_frame = myplot.hole_bounds(what_I_plot_frame,zevmin,zevmax) 796 else: what_I_plot_frame = myplot.bounds(what_I_plot_frame,zevmin,zevmax) 797 if flagnolow: what_I_plot_frame = myplot.nolow(what_I_plot_frame) 733 798 if not tile: 734 799 #zelevels = np.linspace(zevmin*(1. + 1.e-7),zevmax*(1. - 1.e-7)) #,num=20) 735 800 zelevels = np.linspace(zevmin,zevmax,num=ticks) 736 #what_I_plot_frame = smooth(what_I_plot_frame,100)801 #what_I_plot_frame = myplot.smooth(what_I_plot_frame,100) 737 802 if mapmode == 1: m.contourf( x, y, what_I_plot_frame, zelevels, cmap = palette, alpha=trans) 738 elif mapmode == 0: contourf( x, y, what_I_plot_frame, zelevels, cmap = palette, alpha=trans)803 elif mapmode == 0: mpl.pyplot.contourf( x, y, what_I_plot_frame, zelevels, cmap = palette, alpha=trans) 739 804 else: 740 805 if mapmode == 1: m.pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans) 741 elif mapmode == 0: pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans)806 elif mapmode == 0: mpl.pyplot.pcolor( x, y, what_I_plot_frame, cmap = palette, vmin=zevmin, vmax=zevmax, alpha=trans) 742 807 743 808 if (cross is not None or markdevil) and mapmode == 1: … … 748 813 mpl.pyplot.plot([idx],[idy],'+k',mew=0.5,ms=18) 749 814 elif markdevil: 750 idx,idy= find_devil(nc,indextime)815 idx,idy=myplot.find_devil(nc,indextime) 751 816 idx,idy=x[idx,idy],y[idx,idy] 752 817 mpl.pyplot.plot([idx],[idy],'+k',mew=0.5,ms=18) … … 756 821 if (fileref is not None) and ((index_f+1)%3 == 0): daformat = "%.3f" 757 822 elif mult != 1: daformat = "%.1f" 758 else: daformat = fmtvar(fvar.upper())823 else: daformat = myplot.fmtvar(fvar.upper()) 759 824 if proj in ['moll']: zeorientation="horizontal" ; zepad = 0.07 ; zefrac = 0.1 #zepad=0.05 760 825 elif proj in ['cyl']: zeorientation="vertical" ; zepad = 0.03 ; zefrac = 0.023 761 826 else: zeorientation="vertical" ; zepad = 0.03 ; zefrac = 0.05 762 827 if mapmode == 0: zefrac = 0.1 763 zecb = colorbar( fraction=zefrac,pad=zepad,format=daformat,orientation=zeorientation,\828 zecb = mpl.pyplot.colorbar( fraction=zefrac,pad=zepad,format=daformat,orientation=zeorientation,\ 764 829 ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,21])),extend='neither',spacing='proportional' ) 765 830 #if zeorientation == "horizontal": … … 769 834 if winds: 770 835 if typefile in ['meso']: 771 if '9999' not in getattr(nc,'START_DATE') : [vecx_frame,vecy_frame] = [ dumpbdy(vecx_frame,6,stag=uchar,condition=True),dumpbdy(vecy_frame,6,stag=vchar,condition=True)]836 if '9999' not in getattr(nc,'START_DATE') : [vecx_frame,vecy_frame] = [myplot.dumpbdy(vecx_frame,6,stag=uchar,condition=True), myplot.dumpbdy(vecy_frame,6,stag=vchar,condition=True)] 772 837 key = True 773 838 if fvar in ['UV','uv','uvmet']: key = False … … 775 840 key = False 776 841 if metwind and mapmode == 1: [vecx_frame,vecy_frame] = m.rotate_vector(vecx_frame, vecy_frame, lon2d, lat2d) 777 if trans != 0.0: colorvec = definecolorvec(colorb)778 else: colorvec = definecolorvec(back)779 vectorfield(vecx_frame, vecy_frame, x, y, stride=stride, csmooth=2,\842 if trans != 0.0: colorvec = myplot.definecolorvec(colorb) 843 else: colorvec = myplot.definecolorvec(back) 844 myplot.vectorfield(vecx_frame, vecy_frame, x, y, stride=stride, csmooth=2,\ 780 845 #scale=15., factor=300., color=colorvec, key=key) 781 846 scale=20., factor=250./facwind, color=colorvec, key=key) … … 783 848 ### THIS IS A QUITE SPECIFIC PIECE (does not work for mesoscale files) 784 849 if ope == '-_histo' and nplot == numplot: # this should work as long as ope is '-' guarantees 3 plots for 4 panels without contour 785 subplot(subv,subh,nplot+1)786 rcParams["legend.fontsize"] = 'xx-large'850 mpl.pyplot.subplot(subv,subh,nplot+1) 851 mpl.pyplot.rcParams["legend.fontsize"] = 'xx-large' 787 852 if indexlat is None: 788 853 latmin = -50.; latmax = 50. # latitude range for histogram of difference … … 794 859 toplot = np.ravel(what_I_plot_frame[np.isnan(what_I_plot_frame)==False]) 795 860 zebins = np.linspace(zevmin,zevmax,num=30) 796 hist(toplot,bins=zebins,histtype='step',linewidth=2,color='k',normed=True)861 mpl.pyplot.hist(toplot,bins=zebins,histtype='step',linewidth=2,color='k',normed=True) 797 862 zestd = np.std(toplot);zemean = np.mean(toplot) 798 863 zebins = np.linspace(zevmin,zevmax,num=300) 799 864 zegauss = (1./(zestd * np.sqrt(2 * np.pi)) * np.exp( - (zebins - zemean)**2 / (2 * zestd**2) ) ) 800 plot(zebins, zegauss, linewidth=1, color='r',label="mean: "+str(zemean)[0:5]+"\nstd: "+str(zestd)[0:5])801 legend(loc=0,frameon=False)802 subplot(subv,subh,nplot) # go back to last plot for title of contour difference865 mpl.pyplot.plot(zebins, zegauss, linewidth=1, color='r',label="mean: "+str(zemean)[0:5]+"\nstd: "+str(zestd)[0:5]) 866 mpl.pyplot.legend(loc=0,frameon=False) 867 mpl.pyplot.subplot(subv,subh,nplot) # go back to last plot for title of contour difference 803 868 if ope is not None and "only" not in ope: title("fig(1) "+ope+" fig(2)") 804 869 elif ope is not None and "only" in ope: title("fig(1) "+ope[0]+" fig(2)") 805 870 806 871 elif which == "contour": 807 rcParams['contour.negative_linestyle'] = 'solid' # no dashed line for negative values808 zevminc, zevmaxc = calculate_bounds(what_I_plot_frame, vmin=min(what_I_plot_frame), vmax=max(what_I_plot_frame))872 mpl.pyplot.rcParams['contour.negative_linestyle'] = 'solid' # no dashed line for negative values 873 zevminc, zevmaxc = myplot.calculate_bounds(what_I_plot_frame, vmin=min(what_I_plot_frame), vmax=max(what_I_plot_frame)) 809 874 zelevels = np.linspace(zevminc,zevmaxc,ticks/2) #20) 810 875 ### another dirty specific stuff in the wall … … 814 879 elif var2 == 'wstar': zelevels = np.arange(0,10,1.0) 815 880 elif var2 == 'zmax_th': zelevels = np.arange(0,10,2.0) ; what_I_plot_frame = what_I_plot_frame / 1000. 881 elif var2 in ['u']: zelevels = np.arange(-400.,400.,10.) 816 882 ### 817 883 if mapmode == 0: 818 what_I_plot_frame, x, y = define_axis( lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\884 what_I_plot_frame, x, y = myplot.define_axis( lonyeah,latyeah,vert,time,indexlon,indexlat,indexvert,\ 819 885 itime,what_I_plot_frame, len(all_var2[index_f].shape),vertmode,redope) 820 886 ## is this needed? only if len(all_var2[index_f].shape) != len(all_var[index_f].shape) 821 cs = contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33)#, alpha=0.5, linestyles='solid')822 ##cs = contour( x,y,what_I_plot_frame, zelevels, colors='w', linewidths = 0.5)#, alpha=0.5, linestyles='solid')823 # clabel(cs, inline=1, fontsize = 4.*rcParams['font.size']/5., fmt=fmtvar(var2.upper()))887 cs = mpl.pyplot.contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33)#, alpha=0.5, linestyles='solid') 888 ##cs = mpl.pyplot.contour( x,y,what_I_plot_frame, zelevels, colors='w', linewidths = 0.5)#, alpha=0.5, linestyles='solid') 889 #mpl.pyplot.clabel(cs, inline=1, fontsize = 4.*rcParams['font.size']/5., fmt=fmtvar(var2.upper())) 824 890 elif mapmode == 1: 825 891 cs = m.contour( x,y,what_I_plot_frame, zelevels, colors='k', linewidths = 0.33)#, alpha=0.5, linestyles='solid') 826 # clabel(cs, inline=0, fontsize = rcParams['font.size'], fmt="%.0f") #fmtvar(var2.upper()))892 #mpl.pyplot.clabel(cs, inline=0, fontsize = mpl.pyplot.rcParams['font.size'], fmt="%.0f") #myplot.fmtvar(var2.upper())) 827 893 if which in ["regular","unidim","xy"]: 828 894 … … 842 908 mpl.pyplot.yscale('log') 843 909 if invert_y: ax = mpl.pyplot.gca() ; ax.set_ylim(ax.get_ylim()[::-1]) 844 if xlab is not None: xlabel(xlab)845 if ylab is not None: ylabel(ylab)910 if xlab is not None: mpl.pyplot.xlabel(xlab) 911 if ylab is not None: mpl.pyplot.ylabel(ylab) 846 912 847 913 if mrate is not None: … … 855 921 moviename='movie' ;W,H = figframe.canvas.get_width_height() 856 922 video = VideoSink((H,W), moviename, rate=mrate, byteorder="rgba") 857 video.run(mframe) ; close()923 video.run(mframe) ; mpl.pyplot.close() 858 924 if imov == iend: video.close() 859 925 ### THIS IS A WEBPAGE MOVIE 860 926 else: 861 927 nameframe = "image"+str(1000+imov) 862 m akeplotres(nameframe,res=100.,disp=False) ;close()928 myplot.makeplotres(nameframe,res=100.,disp=False) ; mpl.pyplot.close() 863 929 if imov == 0: myfile = open("zepics", 'w') 864 930 myfile.write("modImages["+str(imov)+"] = '"+nameframe+"_100.png';"+ '\n') … … 875 941 zevarname = varname 876 942 if redope is not None: zevarname = zevarname + "_" + redope 877 basename = getname(var=zevarname,var2=var2,winds=winds,anomaly=anomaly)943 basename = myplot.getname(var=zevarname,var2=var2,winds=winds,anomaly=anomaly) 878 944 if len(what_I_plot.shape) > 3: 879 basename = basename + getstralt(nc,level)945 basename = basename + myplot.getstralt(nc,level) 880 946 if mrate is not None: basename = "movie_" + basename 881 947 if typefile in ['meso']: … … 884 950 plottitle = basename 885 951 ### dans le nouveau systeme time=ls,sol,lt cette ligne pourrait ne servir a rien (ou deplacer au dessus) 886 if addchar and indextime is not None: [addchar,gogol,gogol2] = getlschar ( all_namefile[index_f] ) ; plottitle = plottitle + addchar952 if addchar and indextime is not None: [addchar,gogol,gogol2] = myplot.getlschar ( all_namefile[index_f] ) ; plottitle = plottitle + addchar 887 953 ### en fait redope is None doit etre remplace par : n'est ni maxt ni mint 888 954 if redope is None and ltst is not None and ( (mapmode == 0) or (proj in ["lcc","laea","merc","nsper"]) ): plottitle = plottitle + "_LT" + str(ltst) … … 905 971 # if indexvert is not None: plottitle = plottitle + " vert: " + str(min(vert[indexvert])) +" "+ str(max(vert[indexvert])) 906 972 # if indextime is not None: plottitle = plottitle + " time: " + str(min(time[indextime])) +" "+ str(max(time[indextime])) 907 if colorb != "onebar": title( plottitle )973 if colorb != "onebar": mpl.pyplot.title( plottitle ) 908 974 if nplot >= numplot: error = True 909 975 nplot += 1 … … 913 979 914 980 if colorb == "onebar": 915 cax = axes([0.1, 0.2, 0.8, 0.03]) # a ameliorer916 zecb = colorbar(cax=cax, orientation="horizontal", format=fmtvar(fvar.upper()),\981 cax = mpl.pyplot.axes([0.1, 0.2, 0.8, 0.03]) # a ameliorer 982 zecb = mpl.pyplot.colorbar(cax=cax, orientation="horizontal", format=myplot.fmtvar(fvar.upper()),\ 917 983 ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,21])),extend='neither',spacing='proportional') 918 984 if zetitle[0] != "fill": zecb.ax.set_xlabel(zetitle[index_f]) ; zetitle[index_f]="" … … 922 988 ### Save the figure in a file in the data folder or an user-defined folder 923 989 if outputname is None: 924 if typefile in ['meso']: prefix = getprefix(nc)990 if typefile in ['meso']: prefix = myplot.getprefix(nc) 925 991 elif typefile in ['gcm']: prefix = 'LMD_GCM_' 926 992 else: prefix = '' … … 942 1008 print "********** SAVE ", save 943 1009 if save == 'png': 944 if display: m akeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True) ## a miniature945 m akeplotres(zeplot,res=resolution,pad_inches_value=pad_inches_value,disp=False)946 elif save in ['eps','svg','pdf']: m akeplotres(zeplot,pad_inches_value=pad_inches_value,disp=False,ext=save)947 elif save == 'gui': show()1010 if display: myplot.makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True) ## a miniature 1011 myplot.makeplotres(zeplot,res=resolution,pad_inches_value=pad_inches_value,disp=False) 1012 elif save in ['eps','svg','pdf']: myplot.makeplotres(zeplot,pad_inches_value=pad_inches_value,disp=False,ext=save) 1013 elif save == 'gui': mpl.pyplot.show() 948 1014 elif save == 'donothing': pass 949 1015 elif save == 'txt': print "Saved results in txt file." 950 1016 else: 951 1017 print "INFO: save mode not supported. using gui instead." 952 show()1018 mpl.pyplot.show() 953 1019 954 1020 ################################### -
trunk/UTIL/PYTHON/pp.py
r865 r876 8 8 if __name__ == "__main__": 9 9 import sys 10 import myplot 11 import os 12 import numpy as np 13 import netCDF4 14 import glob 15 10 16 from optparse import OptionParser ### to be replaced by argparse 11 17 from api_wrapper import api_onelevel 12 18 from gcm_transformations import call_zrecast 13 from netCDF4 import Dataset14 from myplot import getlschar, separatenames, readslices, adjust_length, whatkindfile, errormess15 from os import system16 19 from planetoplot import planetoplot 17 20 from myscript import getparseroptions 18 import glob19 import numpy as np20 21 21 22 … … 23 24 ### Get options and variables 24 25 parser = OptionParser() ; getparseroptions(parser) ; (opt,args) = parser.parse_args() 25 if opt.file is None: errormess("I want to eat one file at least ! Use pp.py -f name_of_my_file. Or type pp.py -h")26 if opt.var is None and opt.anomaly is True: errormess("Cannot ask to compute anomaly if no variable is set")27 if opt.fref is not None and opt.operat is None: errormess("you must specify an operation when using a reference file")28 if opt.operat in ["+","-"] and opt.fref is None: errormess("you must specifiy a reference file when using inter-file operations")26 if opt.file is None: myplot.errormess("I want to eat one file at least ! Use pp.py -f name_of_my_file. Or type pp.py -h") 27 if opt.var is None and opt.anomaly is True: myplot.errormess("Cannot ask to compute anomaly if no variable is set") 28 if opt.fref is not None and opt.operat is None: myplot.errormess("you must specify an operation when using a reference file") 29 if opt.operat in ["+","-"] and opt.fref is None: myplot.errormess("you must specifiy a reference file when using inter-file operations") 29 30 if opt.fref is not None and opt.operat is not None and opt.itp is not None: interpref=True 30 31 else: interpref=False … … 35 36 ############################# 36 37 ### Get infos about slices 37 zeslat = readslices(opt.slat) ; zeslon = readslices(opt.slon) ; zesvert = readslices(opt.svert) ; zestime =readslices(opt.stime)38 zeslat = myplot.readslices(opt.slat) ; zeslon = myplot.readslices(opt.slon) ; zesvert = myplot.readslices(opt.svert) ; zestime = myplot.readslices(opt.stime) 38 39 39 40 reffile = opt.fref … … 51 52 for i in range(len(opt.file)): 52 53 53 zefiles = separatenames(opt.file[i])54 zefiles = myplot.separatenames(opt.file[i]) 54 55 55 typefile = whatkindfile(Dataset(zefiles[0])) ; stralt = None56 typefile = myplot.whatkindfile(netCDF4.Dataset(zefiles[0])) ; stralt = None 56 57 if typefile in ["meso"]: 57 [lschar,zehour,zehourin] = getlschar ( zefiles[0] )58 [lschar,zehour,zehourin] = myplot.getlschar ( zefiles[0] ) 58 59 if opt.var is None: opt.var = ["HGT"] ; opt.nocolorb = True 59 60 elif typefile in ["geo"]: … … 77 78 for j in range(len(opt.var)): 78 79 79 zevars = separatenames(opt.var[j])80 zevars = myplot.separatenames(opt.var[j]) 80 81 81 inputnvert = separatenames(opt.lvl)82 inputnvert = myplot.separatenames(opt.lvl) 82 83 if np.array(inputnvert).size == 1: 83 84 zelevel = float(inputnvert[0]) … … 104 105 if typefile in ["meso"]: 105 106 if zelevel == 0. and opt.itp == 4: zelevel = 0.010 106 if np.array(inputnvert).size == 1: zesvert = readslices([str(zelevel)])107 if np.array(inputnvert).size == 1: zesvert = myplot.readslices([str(zelevel)]) 107 108 ### winds or no winds 108 109 if opt.winds : zefields = 'uvmet' … … 165 166 name = planetoplot (zefiles,level=int(zelevel),vertmode=opt.itp,\ 166 167 proj=opt.proj,back=opt.back,target=opt.tgt,stride=opt.ste,var=zevars,\ 167 clb= separatenames(opt.clb),winds=opt.winds,\168 clb=myplot.separatenames(opt.clb),winds=opt.winds,\ 168 169 addchar=lschar,vmin=zevmin,vmax=zevmax,\ 169 170 tile=opt.tile,zoom=opt.zoom,display=opt.display,\ 170 171 hole=opt.hole,save=opt.save,\ 171 172 anomaly=opt.anomaly,var2=opt.var2,ndiv=opt.ndiv,\ 172 mult=opt.mult,add=opt.add,zetitle= separatenames(opt.zetitle),\173 mult=opt.mult,add=opt.add,zetitle=myplot.separatenames(opt.zetitle),\ 173 174 slon=zeslon,slat=zeslat,svert=zesvert,stime=zestime,\ 174 175 outputname=opt.out,resolution=opt.res,\ … … 177 178 blat=opt.blat,blon=opt.blon,tsat=opt.tsat,flagnolow=opt.nolow,\ 178 179 mrate=opt.rate,mquality=opt.quality,trans=opt.trans,zarea=opt.area,axtime=opt.axtime,\ 179 redope=opt.redope,seevar=opt.seevar,xlab=opt.xlab,ylab=opt.ylab,lbls= separatenames(opt.labels),\180 lstyle= separatenames(opt.linestyle),cross=readslices(opt.mark),markdevil=opt.mdevil,facwind=opt.facwind,\181 trycol=opt.trycol,streamflag=opt.stream,nocolorb=opt.nocolorb,analysis=opt.analysis )180 redope=opt.redope,seevar=opt.seevar,xlab=opt.xlab,ylab=opt.ylab,lbls=myplot.separatenames(opt.labels),\ 181 lstyle=myplot.separatenames(opt.linestyle),cross=myplot.readslices(opt.mark),markdevil=opt.mdevil,facwind=opt.facwind,\ 182 trycol=opt.trycol,streamflag=opt.stream,nocolorb=opt.nocolorb,analysis=opt.analysis,monster=opt.monster) 182 183 print 'DONE: '+name 183 system("rm -f to_be_erased")184 os.system("rm -f to_be_erased") 184 185 185 186 ######################################################### … … 189 190 #if typefile not in ["meso","mesoapi"]: name = 'pycommand' 190 191 if opt.save == "gui": name = 'pycommand' 191 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)192 elif opt.save == "html": os.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) 192 193 f = open(name+'.sh', 'w') 193 194 f.write(command)
Note: See TracChangeset
for help on using the changeset viewer.
