source: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py @ 346

Last change on this file since 346 was 345, checked in by aslmd, 14 years ago

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