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

Last change on this file since 448 was 448, checked in by aslmd, 13 years ago

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