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

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

PYTHON.

Added new option: --yint, to switch the way vertical coordinates are treated when --vert z1,z2 is specified.

Without specifying --yint, a simple average of the field is done along the vertical coordinate:
mean(vert[z1,z2])

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