source: trunk/UTIL/PYTHON/myplot.py @ 430

Last change on this file since 430 was 430, checked in by acolaitis, 13 years ago


PYTHON

Added... movies !

Example :

pp.py -f wrfout_d01_9999-01-01_07:11:40 --var upward --lat 50 --time 5 -m 0 -M 10 -c Spectral_r --div 20 --ymax 200

will generate a 2D altitude-latitude plot of upward tracer concentration, whereas :

pp.py -f wrfout_d01_9999-01-01_07:11:40 --var upward --lat 50 -m 0 -M 10 -c Spectral_r --div 20 --ymax 200

will generate a movie of the previous plot, along the time axis (default axis for movies for now).

Requirements :

The way the movie maker is written is very light-weight and fast, as no intermediary pictures are saved. The data is buffered and streamed right into an encoder.
I could not make this work with ffmpeg... So you MUST have mencoder installed. It is not on the farm... It is on ciclad, and I know that Aymeric installed it on his machine (at least before getting Ulrich).

Known problems :

I spent litteraly all day making this work, but the option is still basic for now, i.e. it is not yea possible to use -w with movies, or 1D plot movies. (several other stuff are not possible yet, like title specification, axis limits (ylim,xlim), etc... I know how to make it work, I just need some time.

Default name for the output movie is "test.avi" for now.

Tested with LES data without projections. Not test on meso or gcm in mapmode 1.

Details :

A videosink.py
----------------- New class that contains some routines for video stuff. This could either be transformed in a function in an other .py, or other functions could be added to this class to make it more usefull.

M 427 mymath.py
----------------- New routines to convert figures to RGBA data (4-dimensions arrays of pixel data) and figures to PIL images.

M * 428 planetoplot.py
----------------- Added support for movies of 2D plots without overlines

M * 428 myplot.py
----------------- Minor stuff about LES and dumpbdy for W. This could be made in a much cleaner way... maybe later.
File size: 41.0 KB
Line 
1## Author: AS
2def errormess(text,printvar=None):
3    print text
4    if printvar is not None: print printvar
5    exit()
6    return
7
8## Author: AS
9def adjust_length (tab, zelen):
10    from numpy import ones
11    if tab is None:
12        outtab = ones(zelen) * -999999
13    else:
14        if zelen != len(tab):
15            print "not enough or too much values... setting same values all variables"
16            outtab = ones(zelen) * tab[0]
17        else:
18            outtab = tab
19    return outtab
20
21## Author: AS
22def getname(var=False,winds=False,anomaly=False):
23    if var and winds:     basename = var + '_UV'
24    elif var:             basename = var
25    elif winds:           basename = 'UV'
26    else:                 errormess("please set at least winds or var",printvar=nc.variables)
27    if anomaly:           basename = 'd' + basename
28    return basename
29
30## Author: AS
31def localtime(utc,lon):
32    ltst = utc + lon / 15.
33    ltst = int (ltst * 10) / 10.
34    ltst = ltst % 24
35    return ltst
36
37## Author: AS, AC
38def whatkindfile (nc):
39    if 'controle' in nc.variables:             typefile = 'gcm'
40    elif 'phisinit' in nc.variables:           typefile = 'gcm'
41    elif '9999' in getattr(nc,'START_DATE') :  typefile = 'mesoideal'
42    elif 'vert' in nc.variables:               typefile = 'mesoapi'
43    elif 'U' in nc.variables:                  typefile = 'meso'
44    elif 'HGT_M' in nc.variables:              typefile = 'geo'
45    #else:                            errormess("whatkindfile: typefile not supported.")
46    else:                                      typefile = 'gcm' # for lslin-ed files from gcm
47    return typefile
48
49## Author: AS
50def getfield (nc,var):
51    ## this allows to get much faster (than simply referring to nc.variables[var])
52    import numpy as np
53    dimension = len(nc.variables[var].dimensions)
54    #print "   Opening variable with", dimension, "dimensions ..."
55    if dimension == 2:    field = nc.variables[var][:,:]
56    elif dimension == 3:  field = nc.variables[var][:,:,:]
57    elif dimension == 4:  field = nc.variables[var][:,:,:,:]
58    # if there are NaNs in the ncdf, they should be loaded as a masked array which will be
59    # recasted as a regular array later in reducefield
60    if (np.isnan(np.sum(field)) and (type(field).__name__ not in 'MaskedArray')):
61       print "Warning: netcdf as nan values but is not loaded as a Masked Array."
62       print "recasting array type"
63       out=np.ma.masked_invalid(field)
64       out.set_fill_value([np.NaN])
65    else:
66       out=field
67    return out
68
69## Author: AS + TN + AC
70def reducefield (input,d4=None,d3=None,d2=None,d1=None,yint=False,alt=None,anomaly=False):
71    ### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x"
72    ### it would be actually better to name d4 d3 d2 d1 as t z y x
73    ### ... note, anomaly is only computed over d1 and d2 for the moment
74    import numpy as np
75    from mymath import max,mean
76    csmooth = 12 ## a fair amount of grid points (too high results in high computation time)
77    dimension = np.array(input).ndim
78    shape = np.array(input).shape
79    #print 'd1,d2,d3,d4: ',d1,d2,d3,d4
80    print 'IN  REDUCEFIELD dim,shape: ',dimension,shape
81    if anomaly: print 'ANOMALY ANOMALY'
82    output = input
83    error = False
84    #### this is needed to cope the case where d4,d3,d2,d1 are single integers and not arrays
85    if d4 is not None and not isinstance(d4, np.ndarray): d4=[d4]
86    if d3 is not None and not isinstance(d3, np.ndarray): d3=[d3]
87    if d2 is not None and not isinstance(d2, np.ndarray): d2=[d2]
88    if d1 is not None and not isinstance(d1, np.ndarray): d1=[d1]
89    ### now the main part
90    if dimension == 2:
91        if   d2 >= shape[0]: error = True
92        elif d1 >= shape[1]: error = True
93        elif d1 is not None and d2 is not None:  output = mean(input[d2,:],axis=0); output = mean(output[d1],axis=0)
94        elif d1 is not None:         output = mean(input[:,d1],axis=1)
95        elif d2 is not None:         output = mean(input[d2,:],axis=0)
96    elif dimension == 3:
97        if   max(d4) >= shape[0]: error = True
98        elif max(d2) >= shape[1]: error = True
99        elif max(d1) >= shape[2]: error = True
100        elif d4 is not None and d2 is not None and d1 is not None: 
101            output = mean(input[d4,:,:],axis=0); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0)
102        elif d4 is not None and d2 is not None:    output = mean(input[d4,:,:],axis=0); output=mean(output[d2,:],axis=0)
103        elif d4 is not None and d1 is not None:    output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1)
104        elif d2 is not None and d1 is not None:    output = mean(input[:,d2,:],axis=1); output=mean(output[:,d1],axis=1)
105        elif d1 is not None:                       output = mean(input[:,:,d1],axis=2)
106        elif d2 is not None:                       output = mean(input[:,d2,:],axis=1)
107        elif d4 is not None:                       output = mean(input[d4,:,:],axis=0)
108    elif dimension == 4:
109        if   max(d4) >= shape[0]: error = True
110        elif max(d3) >= shape[1]: error = True
111        elif max(d2) >= shape[2]: error = True
112        elif max(d1) >= shape[3]: error = True
113        elif d4 is not None and d3 is not None and d2 is not None and d1 is not None:
114             output = mean(input[d4,:,:,:],axis=0)
115             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
116             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
117             output = mean(output[d2,:],axis=0)
118             output = mean(output[d1],axis=0)
119        elif d4 is not None and d3 is not None and d2 is not None: 
120             output = mean(input[d4,:,:,:],axis=0)
121             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
122             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) 
123             output = mean(output[d2,:],axis=0)
124        elif d4 is not None and d3 is not None and d1 is not None: 
125             output = mean(input[d4,:,:,:],axis=0)
126             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
127             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) 
128             output = mean(output[:,d1],axis=1)
129        elif d4 is not None and d2 is not None and d1 is not None: 
130             output = mean(input[d4,:,:,:],axis=0)
131             if anomaly:
132                 for k in range(output.shape[0]):  output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
133             output = mean(output[:,d2,:],axis=1)
134             output = mean(output[:,d1],axis=1)
135             #noperturb = smooth1d(output,window_len=7)
136             #lenlen = len(output) ; output = output[1:lenlen-7] ; yeye = noperturb[4:lenlen-4]
137             #plot(output) ; plot(yeye) ; show() ; plot(output-yeye) ; show()
138        elif d3 is not None and d2 is not None and d1 is not None: 
139             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3) 
140             if anomaly:
141                 for k in range(output.shape[0]):  output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
142             output = mean(output[:,d2,:],axis=1)
143             output = mean(output[:,d1],axis=1) 
144        elif d4 is not None and d3 is not None: 
145             output = mean(input[d4,:,:,:],axis=0) 
146             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
147             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) 
148        elif d4 is not None and d2 is not None: 
149             output = mean(input[d4,:,:,:],axis=0) 
150             output = mean(output[:,d2,:],axis=1)
151        elif d4 is not None and d1 is not None: 
152             output = mean(input[d4,:,:,:],axis=0) 
153             output = mean(output[:,:,d1],axis=2)
154        elif d3 is not None and d2 is not None: 
155             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
156             output = mean(output[:,d2,:],axis=1)
157        elif d3 is not None and d1 is not None: 
158             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
159             output = mean(output[:,:,d1],axis=0)
160        elif d2 is not None and d1 is not None: 
161             output = mean(input[:,:,d2,:],axis=2)
162             output = mean(output[:,:,d1],axis=2)
163        elif d1 is not None:        output = mean(input[:,:,:,d1],axis=3)
164        elif d2 is not None:        output = mean(input[:,:,d2,:],axis=2)
165        elif d3 is not None:        output = reduce_zaxis(input[:,d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
166        elif d4 is not None:        output = mean(input[d4,:,:,:],axis=0)
167    dimension = np.array(output).ndim
168    shape = np.array(output).shape
169    print 'OUT REDUCEFIELD dim,shape: ',dimension,shape
170    return output, error
171
172## Author: AC + AS
173def reduce_zaxis (input,ax=None,yint=False,vert=None,indice=None):
174    from mymath import max,mean
175    from scipy import integrate
176    if yint and vert is not None and indice is not None:
177       if type(input).__name__=='MaskedArray':
178          input.set_fill_value([np.NaN])
179          output = integrate.trapz(input.filled(),x=vert[indice],axis=ax)
180       else:
181          output = integrate.trapz(input,x=vert[indice],axis=ax)
182    else:
183       output = mean(input,axis=ax)
184    return output
185
186## Author: AS + TN
187def definesubplot ( numplot, fig ):
188    from matplotlib.pyplot import rcParams
189    rcParams['font.size'] = 12. ## default (important for multiple calls)
190    if numplot <= 0:
191        subv = 99999
192        subh = 99999
193    elif numplot == 1: 
194        subv = 99999
195        subh = 99999
196    elif numplot == 2:
197        subv = 1
198        subh = 2
199        fig.subplots_adjust(wspace = 0.35)
200        rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
201    elif numplot == 3:
202        subv = 2
203        subh = 2
204        fig.subplots_adjust(wspace = 0.5)
205        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
206    elif   numplot == 4:
207        subv = 2
208        subh = 2
209        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
210        rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
211    elif numplot <= 6:
212        subv = 2
213        subh = 3
214        fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
215        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
216    elif numplot <= 8:
217        subv = 2
218        subh = 4
219        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
220        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
221    elif numplot <= 9:
222        subv = 3
223        subh = 3
224        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
225        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
226    elif numplot <= 12:
227        subv = 3
228        subh = 4
229        fig.subplots_adjust(wspace = 0.1, hspace = 0.1)
230        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
231    elif numplot <= 16:
232        subv = 4
233        subh = 4
234        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
235        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
236    else:
237        print "number of plot supported: 1 to 16"
238        exit()
239    return subv,subh
240
241## Author: AS
242def getstralt(nc,nvert):
243    typefile = whatkindfile(nc)
244    if typefile is 'meso':                     
245        stralt = "_lvl" + str(nvert)
246    elif typefile is 'mesoapi':
247        zelevel = int(nc.variables['vert'][nvert])
248        if abs(zelevel) < 10000.:   strheight=str(zelevel)+"m"
249        else:                       strheight=str(int(zelevel/1000.))+"km"
250        if 'altitude'       in nc.dimensions:   stralt = "_"+strheight+"-AMR"
251        elif 'altitude_abg' in nc.dimensions:   stralt = "_"+strheight+"-ALS"
252        elif 'bottom_top'   in nc.dimensions:   stralt = "_"+strheight
253        elif 'pressure'     in nc.dimensions:   stralt = "_"+str(zelevel)+"Pa"
254        else:                                   stralt = ""
255    else:
256        stralt = ""
257    return stralt
258
259## Author: AS
260def getlschar ( namefile ):
261    from netCDF4 import Dataset
262    from timestuff import sol2ls
263    from numpy import array
264    from string import rstrip
265    namefile = rstrip( rstrip( rstrip( namefile, chars="_z"), chars="_zabg"), chars="_p")
266    #### we assume that wrfout is next to wrfout_z and wrfout_zabg
267    nc  = Dataset(namefile)
268    zetime = None
269    if 'Times' in nc.variables:
270        zetime = nc.variables['Times'][0]
271        shape = array(nc.variables['Times']).shape
272        if shape[0] < 2: zetime = None
273    if zetime is not None \
274       and 'vert' not in nc.variables:
275        #### strangely enough this does not work for api or ncrcat results!
276        zetimestart = getattr(nc, 'START_DATE')
277        zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9])
278        if zeday < 0:    lschar=""  ## might have crossed a month... fix soon
279        else:            lschar="_Ls"+str( int( 10. * sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) / 10. )
280        ###
281        zetime2 = nc.variables['Times'][1]
282        one  = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37.
283        next = int(zetime2[11]+zetime2[12]) + int(zetime2[14]+zetime2[15])/37. 
284        zehour    = one
285        zehourin  = abs ( next - one )
286    else:
287        lschar=""
288        zehour = 0
289        zehourin = 1 
290    return lschar, zehour, zehourin
291
292## Author: AS
293def getprefix (nc):
294    prefix = 'LMD_MMM_'
295    prefix = prefix + 'd'+str(getattr(nc,'GRID_ID'))+'_'
296    prefix = prefix + str(int(getattr(nc,'DX')/1000.))+'km_'
297    return prefix
298
299## Author: AS
300def getproj (nc):
301    typefile = whatkindfile(nc)
302    if typefile in ['mesoapi','meso','geo']:
303        ### (il faudrait passer CEN_LON dans la projection ?)
304        map_proj = getattr(nc, 'MAP_PROJ')
305        cen_lat  = getattr(nc, 'CEN_LAT')
306        if map_proj == 2:
307            if cen_lat > 10.:   
308                proj="npstere"
309                #print "NP stereographic polar domain"
310            else:           
311                proj="spstere"
312                #print "SP stereographic polar domain"
313        elif map_proj == 1: 
314            #print "lambert projection domain"
315            proj="lcc"
316        elif map_proj == 3: 
317            #print "mercator projection"
318            proj="merc"
319        else:
320            proj="merc"
321    elif typefile in ['gcm']:  proj="cyl"    ## pb avec les autres (de trace derriere la sphere ?)
322    else:                      proj="ortho"
323    return proj   
324
325## Author: AS
326def ptitle (name):
327    from matplotlib.pyplot import title
328    title(name)
329    print name
330
331## Author: AS
332def polarinterv (lon2d,lat2d):
333    import numpy as np
334    wlon = [np.min(lon2d),np.max(lon2d)]
335    ind = np.array(lat2d).shape[0] / 2  ## to get a good boundlat and to get the pole
336    wlat = [np.min(lat2d[ind,:]),np.max(lat2d[ind,:])]
337    return [wlon,wlat]
338
339## Author: AS
340def simplinterv (lon2d,lat2d):
341    import numpy as np
342    return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]]
343
344## Author: AS
345def wrfinterv (lon2d,lat2d):
346    nx = len(lon2d[0,:])-1
347    ny = len(lon2d[:,0])-1
348    lon1 = lon2d[0,0] 
349    lon2 = lon2d[nx,ny] 
350    lat1 = lat2d[0,0] 
351    lat2 = lat2d[nx,ny] 
352    if abs(0.5*(lat1+lat2)) > 60.:  wider = 0.5 * (abs(lon1)+abs(lon2)) * 0.1
353    else:                           wider = 0.
354    if lon1 < lon2:  wlon = [lon1, lon2 + wider] 
355    else:            wlon = [lon2, lon1 + wider]
356    if lat1 < lat2:  wlat = [lat1, lat2]
357    else:            wlat = [lat2, lat1]
358    return [wlon,wlat]
359
360## Author: AS
361def makeplotres (filename,res=None,pad_inches_value=0.25,folder='',disp=True,ext='png',erase=False):
362    import  matplotlib.pyplot as plt
363    from os import system
364    addstr = ""
365    if res is not None:
366        res = int(res)
367        addstr = "_"+str(res)
368    name = filename+addstr+"."+ext
369    if folder != '':      name = folder+'/'+name
370    plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value)
371    if disp:              display(name)
372    if ext in ['eps','ps','svg']:   system("tar czvf "+name+".tar.gz "+name+" ; rm -f "+name)
373    if erase:   system("mv "+name+" to_be_erased")             
374    return
375
376## Author: AS + AC
377def dumpbdy (field,n,stag=None):
378    nx = len(field[0,:])-1
379    ny = len(field[:,0])-1
380    if stag == 'U': nx = nx-1
381    if stag == 'V': ny = ny-1
382    if stag == 'W': nx = nx+1 #special les case when we dump stag on W
383    return field[n:ny-n,n:nx-n]
384
385## Author: AS
386def getcoorddef ( nc ):   
387    import numpy as np
388    ## getcoord2d for predefined types
389    typefile = whatkindfile(nc)
390    if typefile in ['mesoapi','meso']:
391        [lon2d,lat2d] = getcoord2d(nc)
392        lon2d = dumpbdy(lon2d,6)
393        lat2d = dumpbdy(lat2d,6)
394    elif typefile in ['gcm']: 
395        [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
396    elif typefile in ['geo']:
397        [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M')
398    elif typefile in ['mesoideal']:
399        nx=getattr(nc,'WEST-EAST_GRID_DIMENSION')
400        ny=getattr(nc,'SOUTH-NORTH_GRID_DIMENSION')
401        [lon2d,lat2d] = np.meshgrid(np.arange(nx),np.arange(ny))
402    return lon2d,lat2d   
403
404## Author: AS
405def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False):
406    import numpy as np
407    if is1d:
408        lat = nc.variables[nlat][:]
409        lon = nc.variables[nlon][:]
410        [lon2d,lat2d] = np.meshgrid(lon,lat)
411    else:
412        lat = nc.variables[nlat][0,:,:]
413        lon = nc.variables[nlon][0,:,:]
414        [lon2d,lat2d] = [lon,lat]
415    return lon2d,lat2d
416
417## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
418def smooth1d(x,window_len=11,window='hanning'):
419    import numpy
420    """smooth the data using a window with requested size.
421    This method is based on the convolution of a scaled window with the signal.
422    The signal is prepared by introducing reflected copies of the signal
423    (with the window size) in both ends so that transient parts are minimized
424    in the begining and end part of the output signal.
425    input:
426        x: the input signal
427        window_len: the dimension of the smoothing window; should be an odd integer
428        window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
429            flat window will produce a moving average smoothing.
430    output:
431        the smoothed signal
432    example:
433    t=linspace(-2,2,0.1)
434    x=sin(t)+randn(len(t))*0.1
435    y=smooth(x)
436    see also:
437    numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
438    scipy.signal.lfilter
439    TODO: the window parameter could be the window itself if an array instead of a string   
440    """
441    x = numpy.array(x)
442    if x.ndim != 1:
443        raise ValueError, "smooth only accepts 1 dimension arrays."
444    if x.size < window_len:
445        raise ValueError, "Input vector needs to be bigger than window size."
446    if window_len<3:
447        return x
448    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
449        raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
450    s=numpy.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
451    #print(len(s))
452    if window == 'flat': #moving average
453        w=numpy.ones(window_len,'d')
454    else:
455        w=eval('numpy.'+window+'(window_len)')
456    y=numpy.convolve(w/w.sum(),s,mode='valid')
457    return y
458
459## Author: AS
460def smooth (field, coeff):
461        ## actually blur_image could work with different coeff on x and y
462        if coeff > 1:   result = blur_image(field,int(coeff))
463        else:           result = field
464        return result
465
466## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
467def gauss_kern(size, sizey=None):
468        import numpy as np
469        # Returns a normalized 2D gauss kernel array for convolutions
470        size = int(size)
471        if not sizey:
472                sizey = size
473        else:
474                sizey = int(sizey)
475        x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
476        g = np.exp(-(x**2/float(size)+y**2/float(sizey)))
477        return g / g.sum()
478
479## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
480def blur_image(im, n, ny=None) :
481        from scipy.signal import convolve
482        # blurs the image by convolving with a gaussian kernel of typical size n.
483        # The optional keyword argument ny allows for a different size in the y direction.
484        g = gauss_kern(n, sizey=ny)
485        improc = convolve(im, g, mode='same')
486        return improc
487
488## Author: AS
489def getwinddef (nc):   
490    ## getwinds for predefined types
491    typefile = whatkindfile(nc)
492    ###
493    if typefile is 'mesoapi':                  [uchar,vchar] = ['Um','Vm']
494    elif typefile is 'gcm':                    [uchar,vchar] = ['u','v']
495    elif typefile in ['meso','mesoideal']:     [uchar,vchar] = ['U','V']
496    else:                                      [uchar,vchar] = ['not found','not found']
497    ###
498    if typefile in ['meso']:     metwind = False ## geometrical (wrt grid)
499    else:                        metwind = True  ## meteorological (zon/mer)
500    if metwind is False:         print "Not using meteorological winds. You trust numerical grid as being (x,y)"
501    ###
502    return uchar,vchar,metwind
503
504## Author: AS
505def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True):
506    ## scale regle la reference du vecteur
507    ## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs.
508    import  matplotlib.pyplot               as plt
509    import  numpy                           as np
510    posx = np.min(x) - np.std(x) / 10.
511    posy = np.min(y) - np.std(y) / 10.
512    u = smooth(u,csmooth)
513    v = smooth(v,csmooth)
514    widthvec = 0.003 #0.005 #0.003
515    q = plt.quiver( x[::stride,::stride],\
516                    y[::stride,::stride],\
517                    u[::stride,::stride],\
518                    v[::stride,::stride],\
519                    angles='xy',color=color,pivot='middle',\
520                    scale=factor,width=widthvec )
521    if color in ['white','yellow']:     kcolor='black'
522    else:                               kcolor=color
523    if key: p = plt.quiverkey(q,posx,posy,scale,\
524                   str(int(scale)),coordinates='data',color=kcolor,labelpos='S',labelsep = 0.03)
525    return 
526
527## Author: AS
528def display (name):
529    from os import system
530    system("display "+name+" > /dev/null 2> /dev/null &")
531    return name
532
533## Author: AS
534def findstep (wlon):
535    steplon = int((wlon[1]-wlon[0])/4.)  #3
536    step = 120.
537    while step > steplon and step > 15. :       step = step / 2.
538    if step <= 15.:
539        while step > steplon and step > 5.  :   step = step - 5.
540    if step <= 5.:
541        while step > steplon and step > 1.  :   step = step - 1.
542    if step <= 1.:
543        step = 1. 
544    return step
545
546## Author: AS
547def define_proj (char,wlon,wlat,back=None,blat=None):
548    from    mpl_toolkits.basemap            import Basemap
549    import  numpy                           as np
550    import  matplotlib                      as mpl
551    from mymath import max
552    meanlon = 0.5*(wlon[0]+wlon[1])
553    meanlat = 0.5*(wlat[0]+wlat[1])
554    if blat is None:
555        ortholat=meanlat
556        if   wlat[0] >= 80.:   blat =  40. 
557        elif wlat[1] <= -80.:  blat = -40.
558        elif wlat[1] >= 0.:    blat = wlat[0] 
559        elif wlat[0] <= 0.:    blat = wlat[1]
560    else:  ortholat=blat
561    #print "blat ", blat
562    h = 50.  ## en km
563    radius = 3397200.
564    if   char == "cyl":     m = Basemap(rsphere=radius,projection='cyl',\
565                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
566    elif char == "moll":    m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon)
567    elif char == "ortho":   m = Basemap(rsphere=radius,projection='ortho',lon_0=meanlon,lat_0=ortholat)
568    elif char == "lcc":     m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\
569                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
570    elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.)
571    elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=180.)
572    elif char == "nplaea":  m = Basemap(rsphere=radius,projection='nplaea', boundinglat=wlat[0], lon_0=meanlon)
573    elif char == "laea":    m = Basemap(rsphere=radius,projection='laea',lon_0=meanlon,lat_0=meanlat,lat_ts=meanlat,\
574                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
575    elif char == "nsper":   m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.)
576    elif char == "merc":    m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\
577                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
578    fontsizemer = int(mpl.rcParams['font.size']*3./4.)
579    if char in ["cyl","lcc","merc","nsper","laea"]:   step = findstep(wlon)
580    else:                                             step = 10.
581    steplon = step*2.
582    #if back in ["geolocal"]:                         
583    #    step = np.min([5.,step])
584    #    steplon = step
585    m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color='grey', fontsize=fontsizemer)
586    m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer)
587    if back: m.warpimage(marsmap(back),scale=0.75)
588            #if not back:
589            #    if not var:                                        back = "mola"    ## if no var:         draw mola
590            #    elif typefile in ['mesoapi','meso','geo'] \
591            #       and proj not in ['merc','lcc','nsper','laea']:  back = "molabw"  ## if var but meso:   draw molabw
592            #    else:                                              pass             ## else:              draw None
593    return m
594
595## Author: AS
596#### test temporaire
597def putpoints (map,plot):
598    #### from http://www.scipy.org/Cookbook/Matplotlib/Maps
599    # lat/lon coordinates of five cities.
600    lats = [18.4]
601    lons = [-134.0]
602    points=['Olympus Mons']
603    # compute the native map projection coordinates for cities.
604    x,y = map(lons,lats)
605    # plot filled circles at the locations of the cities.
606    map.plot(x,y,'bo')
607    # plot the names of those five cities.
608    wherept = 0 #1000 #50000
609    for name,xpt,ypt in zip(points,x,y):
610       plot.text(xpt+wherept,ypt+wherept,name)
611    ## le nom ne s'affiche pas...
612    return
613
614## Author: AS
615def calculate_bounds(field,vmin=None,vmax=None):
616    import numpy as np
617    from mymath import max,min,mean
618    ind = np.where(field < 9e+35)
619    fieldcalc = field[ ind ] # la syntaxe compacte ne marche si field est un tuple
620    ###
621    dev = np.std(fieldcalc)*3.0
622    ###
623    if vmin is None:
624        zevmin = mean(fieldcalc) - dev
625    else:             zevmin = vmin
626    ###
627    if vmax is None:  zevmax = mean(fieldcalc) + dev
628    else:             zevmax = vmax
629    if vmin == vmax:
630                      zevmin = mean(fieldcalc) - dev  ### for continuity
631                      zevmax = mean(fieldcalc) + dev  ### for continuity           
632    ###
633    if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0.
634    print "BOUNDS field ", min(fieldcalc), max(fieldcalc)
635    print "BOUNDS adopted ", zevmin, zevmax
636    return zevmin, zevmax
637
638## Author: AS
639def bounds(what_I_plot,zevmin,zevmax):
640    from mymath import max,min,mean
641    ### might be convenient to add the missing value in arguments
642    #what_I_plot[ what_I_plot < zevmin ] = zevmin#*(1. + 1.e-7)
643    if zevmin < 0: what_I_plot[ what_I_plot < zevmin*(1. - 1.e-7) ] = zevmin*(1. - 1.e-7)
644    else:          what_I_plot[ what_I_plot < zevmin*(1. + 1.e-7) ] = zevmin*(1. + 1.e-7)
645    print "NEW MIN ", min(what_I_plot)
646    what_I_plot[ what_I_plot > 9e+35  ] = -9e+35
647    what_I_plot[ what_I_plot > zevmax ] = zevmax
648    print "NEW MAX ", max(what_I_plot)
649    return what_I_plot
650
651## Author: AS
652def nolow(what_I_plot):
653    from mymath import max,min
654    lim = 0.15*0.5*(abs(max(what_I_plot))+abs(min(what_I_plot)))
655    print "NO PLOT BELOW VALUE ", lim
656    what_I_plot [ abs(what_I_plot) < lim ] = 1.e40 
657    return what_I_plot
658
659
660## Author : AC
661def hole_bounds(what_I_plot,zevmin,zevmax):
662    import numpy as np
663    zi=0
664    for i in what_I_plot:
665        zj=0
666        for j in i:
667            if ((j < zevmin) or (j > zevmax)):what_I_plot[zi,zj]=np.NaN
668            zj=zj+1
669        zi=zi+1
670
671    return what_I_plot
672
673## Author: AS
674def zoomset (wlon,wlat,zoom):
675    dlon = abs(wlon[1]-wlon[0])/2.
676    dlat = abs(wlat[1]-wlat[0])/2.
677    [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\
678                    [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ]
679    print "ZOOM %",zoom,wlon,wlat
680    return wlon,wlat
681
682## Author: AS
683def fmtvar (whichvar="def"):
684    fmtvar    =     { \
685             "TK":           "%.0f",\
686             # Variables from TES ncdf format
687             "T_NADIR_DAY":  "%.0f",\
688             "T_NADIR_NIT":  "%.0f",\
689             # Variables from tes.py ncdf format
690             "TEMP_DAY":     "%.0f",\
691             "TEMP_NIGHT":   "%.0f",\
692             # Variables from MCS and mcs.py ncdf format
693             "DTEMP":        "%.0f",\
694             "NTEMP":        "%.0f",\
695             "DNUMBINTEMP":  "%.0f",\
696             "NNUMBINTEMP":  "%.0f",\
697             # other stuff
698             "TPOT":         "%.0f",\
699             "TSURF":        "%.0f",\
700             "def":          "%.1e",\
701             "PTOT":         "%.0f",\
702             "HGT":          "%.1e",\
703             "USTM":         "%.2f",\
704             "HFX":          "%.0f",\
705             "ICETOT":       "%.1e",\
706             "TAU_ICE":      "%.2f",\
707             "VMR_ICE":      "%.1e",\
708             "MTOT":         "%.1f",\
709             "ANOMALY":      "%.1f",\
710             "W":            "%.1f",\
711             "WMAX_TH":      "%.1f",\
712             "QSURFICE":     "%.0f",\
713             "UM":           "%.0f",\
714             "ALBBARE":      "%.2f",\
715             "TAU":          "%.1f",\
716             "CO2":          "%.2f",\
717             #### T.N.
718             "TEMP":         "%.0f",\
719             "VMR_H2OICE":   "%.0f",\
720             "VMR_H2OVAP":   "%.0f",\
721             "TAUTES":       "%.2f",\
722             "TAUTESAP":     "%.2f",\
723
724                    }
725    if whichvar not in fmtvar:
726        whichvar = "def"
727    return fmtvar[whichvar]
728
729## Author: AS
730####################################################################################################################
731### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
732def defcolorb (whichone="def"):
733    whichcolorb =    { \
734             "def":          "spectral",\
735             "HGT":          "spectral",\
736             "HGT_M":        "spectral",\
737             "TK":           "gist_heat",\
738             "TPOT":         "Paired",\
739             "TSURF":        "RdBu_r",\
740             "QH2O":         "PuBu",\
741             "USTM":         "YlOrRd",\
742             #"T_nadir_nit":  "RdBu_r",\
743             #"T_nadir_day":  "RdBu_r",\
744             "HFX":          "RdYlBu",\
745             "ICETOT":       "YlGnBu_r",\
746             #"MTOT":         "PuBu",\
747             "CCNQ":         "YlOrBr",\
748             "CCNN":         "YlOrBr",\
749             "TEMP":         "Jet",\
750             "TAU_ICE":      "Blues",\
751             "VMR_ICE":      "Blues",\
752             "W":            "jet",\
753             "WMAX_TH":      "spectral",\
754             "ANOMALY":      "RdBu_r",\
755             "QSURFICE":     "hot_r",\
756             "ALBBARE":      "spectral",\
757             "TAU":          "YlOrBr_r",\
758             "CO2":          "YlOrBr_r",\
759             #### T.N.
760             "MTOT":         "Jet",\
761             "H2O_ICE_S":    "RdBu",\
762             "VMR_H2OICE":   "PuBu",\
763             "VMR_H2OVAP":   "PuBu",\
764                     }
765#W --> spectral ou jet
766#spectral BrBG RdBu_r
767    #print "predefined colorbars"
768    if whichone not in whichcolorb:
769        whichone = "def"
770    return whichcolorb[whichone]
771
772## Author: AS
773def definecolorvec (whichone="def"):
774        whichcolor =    { \
775                "def":          "black",\
776                "vis":          "yellow",\
777                "vishires":     "yellow",\
778                "molabw":       "yellow",\
779                "mola":         "black",\
780                "gist_heat":    "white",\
781                "hot":          "tk",\
782                "gist_rainbow": "black",\
783                "spectral":     "black",\
784                "gray":         "red",\
785                "PuBu":         "black",\
786                        }
787        if whichone not in whichcolor:
788                whichone = "def"
789        return whichcolor[whichone]
790
791## Author: AS
792def marsmap (whichone="vishires"):
793        from os import uname
794        mymachine = uname()[1]
795        ### not sure about speed-up with this method... looks the same
796        if "lmd.jussieu.fr" in mymachine: domain = "/u/aslmd/WWW/maps/"
797        else:                             domain = "http://www.lmd.jussieu.fr/~aslmd/maps/"
798        whichlink =     { \
799                #"vis":         "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\
800                #"vishires":    "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\
801                #"geolocal":    "http://dl.dropbox.com/u/11078310/geolocal.jpg",\
802                #"mola":        "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\
803                #"molabw":      "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\
804                "vis":         domain+"mar0kuu2.jpg",\
805                "vishires":    domain+"MarsMap_2500x1250.jpg",\
806                "geolocal":    domain+"geolocal.jpg",\
807                "mola":        domain+"mars-mola-2k.jpg",\
808                "molabw":      domain+"MarsElevation_2500x1250.jpg",\
809                "clouds":      "http://www.johnstonsarchive.net/spaceart/marswcloudmap.jpg",\
810                "jupiter":     "http://www.mmedia.is/~bjj/data/jupiter_css/jupiter_css.jpg",\
811                "jupiter_voy": "http://www.mmedia.is/~bjj/data/jupiter/jupiter_vgr2.jpg",\
812                "bw":          "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg",\
813                "contrast":    "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg",\
814                "nice":        "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg",\
815                "blue":        "http://eoimages.gsfc.nasa.gov/ve/2430/land_ocean_ice_2048.jpg",\
816                "blueclouds":  "http://eoimages.gsfc.nasa.gov/ve/2431/land_ocean_ice_cloud_2048.jpg",\
817                "justclouds":  "http://eoimages.gsfc.nasa.gov/ve/2432/cloud_combined_2048.jpg",\
818                        }
819        ### see http://www.mmedia.is/~bjj/planetary_maps.html
820        if whichone not in whichlink: 
821                print "marsmap: choice not defined... you'll get the default one... "
822                whichone = "vishires" 
823        return whichlink[whichone]
824
825#def earthmap (whichone):
826#       if   whichone == "contrast":    whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg"
827#       elif whichone == "bw":          whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg"
828#       elif whichone == "nice":        whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg"
829#       return whichlink
830
831## Author: AS
832def latinterv (area="Whole"):
833    list =    { \
834        "Europe":                [[ 20., 80.],[- 50.,  50.]],\
835        "Central_America":       [[-10., 40.],[ 230., 300.]],\
836        "Africa":                [[-20., 50.],[- 50.,  50.]],\
837        "Whole":                 [[-90., 90.],[-180., 180.]],\
838        "Southern_Hemisphere":   [[-90., 60.],[-180., 180.]],\
839        "Northern_Hemisphere":   [[-60., 90.],[-180., 180.]],\
840        "Tharsis":               [[-30., 60.],[-170.,- 10.]],\
841        "Whole_No_High":         [[-60., 60.],[-180., 180.]],\
842        "Chryse":                [[-60., 60.],[- 60.,  60.]],\
843        "North_Pole":            [[ 50., 90.],[-180., 180.]],\
844        "Close_North_Pole":      [[ 75., 90.],[-180., 180.]],\
845        "Far_South_Pole":        [[-90.,-40.],[-180., 180.]],\
846        "South_Pole":            [[-90.,-50.],[-180., 180.]],\
847        "Close_South_Pole":      [[-90.,-75.],[-180., 180.]],\
848              }
849    if area not in list:   area = "Whole"
850    [olat,olon] = list[area]
851    return olon,olat
852
853## Author: TN
854def separatenames (name):
855  from numpy import concatenate
856  # look for comas in the input name to separate different names (files, variables,etc ..)
857  if name is None:
858     names = None
859  else:
860    names = []
861    stop = 0
862    currentname = name
863    while stop == 0:
864      indexvir = currentname.find(',')
865      if indexvir == -1:
866        stop = 1
867        name1 = currentname
868      else:
869        name1 = currentname[0:indexvir]
870      names = concatenate((names,[name1]))
871      currentname = currentname[indexvir+1:len(currentname)]
872  return names
873
874## Author: TN [old]
875def getopmatrix (kind,n):
876  import numpy as np
877  # compute matrix of operations between files
878  # the matrix is 'number of files'-square
879  # 1: difference (row minus column), 2: add
880  # not 0 in diag : just plot
881  if n == 1:
882    opm = np.eye(1)
883  elif kind == 'basic':
884    opm = np.eye(n)
885  elif kind == 'difference1': # show differences with 1st file
886    opm = np.zeros((n,n))
887    opm[0,:] = 1
888    opm[0,0] = 0
889  elif kind == 'difference2': # show differences with 1st file AND show 1st file
890    opm = np.zeros((n,n))
891    opm[0,:] = 1
892  else:
893    opm = np.eye(n)
894  return opm
895 
896## Author: TN [old]
897def checkcoherence (nfiles,nlat,nlon,ntime):
898  if (nfiles > 1):
899     if (nlat > 1):
900        errormess("what you asked is not possible !")
901  return 1
902
903## Author: TN
904def readslices(saxis):
905  from numpy import empty
906  if saxis == None:
907     zesaxis = None
908  else:
909     zesaxis = empty((len(saxis),2))
910     for i in range(len(saxis)):
911        a = separatenames(saxis[i])
912        if len(a) == 1:
913           zesaxis[i,:] = float(a[0])
914        else:
915           zesaxis[i,0] = float(a[0])
916           zesaxis[i,1] = float(a[1])
917           
918  return zesaxis
919
920## Author: AS
921def bidimfind(lon2d,lat2d,vlon,vlat):
922   import numpy as np
923   if vlat is None:    array = (lon2d - vlon)**2
924   elif vlon is None:  array = (lat2d - vlat)**2
925   else:               array = (lon2d - vlon)**2 + (lat2d - vlat)**2
926   idy,idx = np.unravel_index( np.argmin(array), lon2d.shape )
927   if vlon is not None:
928       #print lon2d[idy,idx],vlon
929       if (np.abs(lon2d[idy,idx]-vlon)) > 5: errormess("longitude not found ",printvar=lon2d)
930   if vlat is not None:
931       #print lat2d[idy,idx],vlat
932       if (np.abs(lat2d[idy,idx]-vlat)) > 5: errormess("latitude not found ",printvar=lat2d)
933   return idx,idy
934
935## Author: TN
936def getsindex(saxis,index,axis):
937# input  : all the desired slices and the good index
938# output : all indexes to be taken into account for reducing field
939  import numpy as np
940  ### added by AS to treat the case of stime = - LT
941  if saxis is not None:
942      if saxis[0][0] < 0: saxis = - saxis
943  ###
944  if ( np.array(axis).ndim == 2):
945      axis = axis[:,0]
946  if saxis is None:
947      zeindex = None
948  else:
949      aaa = int(np.argmin(abs(saxis[index,0] - axis)))
950      bbb = int(np.argmin(abs(saxis[index,1] - axis)))
951      [imin,imax] = np.sort(np.array([aaa,bbb]))
952      zeindex = np.array(range(imax-imin+1))+imin
953      # because -180 and 180 are the same point in longitude,
954      # we get rid of one for averaging purposes.
955      if axis[imin] == -180 and axis[imax] == 180:
956         zeindex = zeindex[0:len(zeindex)-1]
957         print "INFO: whole longitude averaging asked, so last point is not taken into account."
958  return zeindex
959     
960## Author: TN
961def define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,indextime,what_I_plot,dim0,vertmode):
962# Purpose of define_axis is to find x and y axis scales in a smart way
963# x axis priority: 1/time 2/lon 3/lat 4/vertical
964# To be improved !!!...
965   from numpy import array,swapaxes
966   x = None
967   y = None
968   count = 0
969   what_I_plot = array(what_I_plot)
970   shape = what_I_plot.shape
971   if indextime is None:
972      print "AXIS is time"
973      x = time
974      count = count+1
975   if indexlon is None:
976      print "AXIS is lon"
977      if count == 0: x = lon
978      else: y = lon
979      count = count+1
980   if indexlat is None:
981      print "AXIS is lat"
982      if count == 0: x = lat
983      else: y = lat
984      count = count+1
985   if indexvert is None and dim0 is 4:
986      print "AXIS is vert"
987      if vertmode == 0: # vertical axis is as is (GCM grid)
988         if count == 0: x=range(len(vert))
989         else: y=range(len(vert))
990         count = count+1
991      else: # vertical axis is in kms
992         if count == 0: x = vert
993         else: y = vert
994         count = count+1
995   x = array(x)
996   y = array(y)
997   print "CHECK: what_I_plot.shape", what_I_plot.shape
998   print "CHECK: x.shape, y.shape", x.shape, y.shape
999   if len(shape) == 1:
1000      if shape[0] != len(x):
1001         print "WARNING HERE !!!"
1002         x = y
1003   elif len(shape) == 2:
1004      if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]:
1005         what_I_plot = swapaxes(what_I_plot,0,1)
1006         print "INFO: swapaxes", what_I_plot.shape, shape
1007   return what_I_plot,x,y
1008
1009# Author: TN + AS
1010def determineplot(slon, slat, svert, stime):
1011    nlon = 1 # number of longitudinal slices -- 1 is None
1012    nlat = 1
1013    nvert = 1
1014    ntime = 1
1015    nslices = 1
1016    if slon is not None:
1017        nslices = nslices*len(slon)
1018        nlon = len(slon)
1019    if slat is not None:
1020        nslices = nslices*len(slat)
1021        nlat = len(slat)
1022    if svert is not None:
1023        nslices = nslices*len(svert)
1024        nvert = len(svert)
1025    if stime is not None:
1026        nslices = nslices*len(stime)
1027        ntime = len(stime)
1028    #else:
1029    #    nslices = 2 
1030    mapmode = 0
1031    if slon is None and slat is None:
1032       mapmode = 1 # in this case we plot a map, with the given projection
1033
1034    return nlon, nlat, nvert, ntime, mapmode, nslices
Note: See TracBrowser for help on using the repository browser.