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

Last change on this file since 453 was 453, checked in by tnavarro, 13 years ago

GRAPHICS: transparent fields and 3 new background images

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