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

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

MESOSCALE:
-- Code in storm scenario corresponds to 'OMEGA' reference case [Julien Faure]
-- Easier settings for dust lifting without the need to recompile [see below]
-- A few modifications to plot and dust-devil detection PYTHON routines

29/09/11 == AS

--> To easily explore sensitivity to lifting thresholds: in dustlift.F, ustar_seuil=sqrt(stress/rho)

and alpha_lift[dust_mass] can be prescribed through an external stress.def parameter file.
--- alpha_lift[dust_number] is computed from alpha_lift[dust_mass] as in initracer.F
--- ustar_seuil is more user-friendly than stress because direct comparison with ustar from model

--> For the moment this is MESOSCALE only, but potentially useful to everyone
File size: 25.8 KB
Line 
1def errormess(text,printvar=None):
2    print text
3    if printvar: print printvar
4    exit()
5    return
6
7def getname(var=False,winds=False,anomaly=False):
8    if var and winds:     basename = var + '_UV'
9    elif var:             basename = var
10    elif winds:           basename = 'UV'
11    else:                 errormess("please set at least winds or var",printvar=nc.variables)
12    if anomaly:           basename = 'd' + basename
13    return basename
14
15def localtime(utc,lon):
16    ltst = utc + lon / 15.
17    ltst = int (ltst * 10) / 10.
18    ltst = ltst % 24
19    return ltst
20
21def whatkindfile (nc):
22    if 'controle' in nc.variables:   typefile = 'gcm'
23    elif 'phisinit' in nc.variables: typefile = 'gcmex'
24    elif 'vert' in nc.variables:     typefile = 'mesoapi'
25    elif 'U' in nc.variables:        typefile = 'meso'
26    elif 'HGT_M' in nc.variables:    typefile = 'geo'
27    else:                            errormess("whatkindfile: typefile not supported.")
28    return typefile
29
30def getfield (nc,var):
31    ## this allows to get much faster (than simply referring to nc.variables[var])
32    dimension = len(nc.variables[var].dimensions)
33    if dimension == 2:    field = nc.variables[var][:,:]
34    elif dimension == 3:  field = nc.variables[var][:,:,:]
35    elif dimension == 4:  field = nc.variables[var][:,:,:,:]
36    return field
37
38def reducefield (input,d4=None,d3=None,d2=None,d1=None):
39    ### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x"
40    ### it would be actually better to name d4 d3 d2 d1 as t z y x
41    import numpy as np
42    dimension = np.array(input).ndim
43    shape = np.array(input).shape
44    print 'dim,shape: ',dimension,shape
45    output = input
46    error = False
47    if dimension == 2:
48        if   d2 >= shape[0]: error = True
49        elif d1 >= shape[1]: error = True
50        elif d1 is not None and d2 is not None:  output = input[d2,d1]
51        elif d1 is not None:         output = input[:,d1]
52        elif d2 is not None:         output = input[d2,:]
53    elif dimension == 3:
54        if   d4 >= shape[0]: error = True
55        elif d2 >= shape[1]: error = True
56        elif d1 >= shape[2]: error = True
57        elif d4 is not None and d2 is not None and d1 is not None:  output = input[d4,d2,d1]
58        elif d4 is not None and d2 is not None:         output = input[d4,d2,:]
59        elif d4 is not None and d1 is not None:         output = input[d4,:,d1]
60        elif d2 is not None and d1 is not None:         output = input[:,d2,d1]
61        elif d1 is not None:                output = input[:,:,d1]
62        elif d2 is not None:                output = input[:,d2,:]
63        elif d4 is not None:                output = input[d4,:,:]
64    elif dimension == 4:
65        if   d4 >= shape[0]: error = True
66        elif d3 >= shape[1]: error = True
67        elif d2 >= shape[2]: error = True
68        elif d1 >= shape[3]: error = True
69        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]
70        elif d4 is not None and d3 is not None and d2 is not None:         output = input[d4,d3,d2,:]
71        elif d4 is not None and d3 is not None and d1 is not None:         output = input[d4,d3,:,d1]
72        elif d4 is not None and d2 is not None and d1 is not None:         output = input[d4,:,d2,d1]
73        elif d3 is not None and d2 is not None and d1 is not None:         output = input[:,d3,d2,d1]
74        elif d4 is not None and d3 is not None:                output = input[d4,d3,:,:]
75        elif d4 is not None and d2 is not None:                output = input[d4,:,d2,:]
76        elif d4 is not None and d1 is not None:                output = input[d4,:,:,d1]
77        elif d3 is not None and d2 is not None:                output = input[:,d3,d2,:]
78        elif d3 is not None and d1 is not None:                output = input[:,d3,:,d1]
79        elif d2 is not None and d1 is not None:                output = input[:,:,d2,d1]
80        elif d1 is not None:                       output = input[:,:,:,d1]
81        elif d2 is not None:                       output = input[:,:,d2,:]
82        elif d3 is not None:                       output = input[:,d3,:,:]
83        elif d4 is not None:                       output = input[d4,:,:,:]
84    dimension = np.array(output).ndim
85    shape = np.array(output).shape
86    print 'dim,shape: ',dimension,shape
87    return output, error
88
89def definesubplot ( numplot, fig ):
90    from matplotlib.pyplot import rcParams
91    rcParams['font.size'] = 12. ## default (important for multiple calls)
92    if   numplot == 4:
93        sub = 221
94        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
95        rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
96    elif numplot == 2:
97        sub = 121
98        fig.subplots_adjust(wspace = 0.35)
99        rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
100    elif numplot == 3:
101        sub = 131
102        fig.subplots_adjust(wspace = 0.5)
103        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
104    elif numplot == 6:
105        sub = 231
106        fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
107        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
108    elif numplot == 8:
109        sub = 331 #241
110        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
111        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
112    elif numplot == 9:
113        sub = 331
114        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
115        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
116    elif numplot == 1:
117        sub = 99999
118    elif numplot <= 0: 
119        sub = 99999
120    else:
121        print "supported: 1,2,3,4,6,8,9"
122        exit()
123    return sub
124
125def getstralt(nc,nvert):
126    typefile = whatkindfile(nc)
127    if typefile is 'meso':                     
128        stralt = "_lvl" + str(nvert)
129    elif typefile is 'mesoapi':
130        zelevel = int(nc.variables['vert'][nvert])
131        if abs(zelevel) < 10000.:   strheight=str(zelevel)+"m"
132        else:                       strheight=str(int(zelevel/1000.))+"km"
133        if 'altitude'       in nc.dimensions:   stralt = "_"+strheight+"-AMR"
134        elif 'altitude_abg' in nc.dimensions:   stralt = "_"+strheight+"-ALS"
135        elif 'bottom_top'   in nc.dimensions:   stralt = "_"+strheight
136        elif 'pressure'     in nc.dimensions:   stralt = "_"+str(zelevel)+"Pa"
137        else:                                   stralt = ""
138    else:
139        stralt = ""
140    return stralt
141
142def getlschar ( namefile ):
143    from netCDF4 import Dataset
144    from timestuff import sol2ls
145    from numpy import array
146    nc  = Dataset(namefile)
147    zetime = None
148    if 'Times' in nc.variables: 
149        zetime = nc.variables['Times'][0]
150        shape = array(nc.variables['Times']).shape
151        if shape[0] < 2: zetime = None
152    if zetime is not None \
153       and 'vert' not in nc.variables:
154       #### strangely enough this does not work for api or ncrcat results!
155        zetimestart = getattr(nc, 'START_DATE')
156        zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9])
157        if zeday < 0:    lschar=""  ## might have crossed a month... fix soon
158        else:            lschar="_Ls"+str( int( 10. * sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) / 10. )
159        ###
160        zetime2 = nc.variables['Times'][1]
161        one  = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37.
162        next = int(zetime2[11]+zetime2[12]) + int(zetime2[14]+zetime2[15])/37. 
163        zehour    = one
164        zehourin  = abs ( next - one )
165    else:
166        lschar=""
167        zehour = 0
168        zehourin = 1 
169    return lschar, zehour, zehourin
170
171def getprefix (nc):
172    prefix = 'LMD_MMM_'
173    prefix = prefix + 'd'+str(getattr(nc,'GRID_ID'))+'_'
174    prefix = prefix + str(int(getattr(nc,'DX')/1000.))+'km_'
175    return prefix
176
177def getproj (nc):
178    typefile = whatkindfile(nc)
179    if typefile in ['mesoapi','meso','geo']:
180        ### (il faudrait passer CEN_LON dans la projection ?)
181        map_proj = getattr(nc, 'MAP_PROJ')
182        cen_lat  = getattr(nc, 'CEN_LAT')
183        if map_proj == 2:
184            if cen_lat > 10.:   
185                proj="npstere"
186                print "NP stereographic polar domain" 
187            else:           
188                proj="spstere"
189                print "SP stereographic polar domain"
190        elif map_proj == 1: 
191            print "lambert projection domain" 
192            proj="lcc"
193        elif map_proj == 3: 
194            print "mercator projection"
195            proj="merc"
196        else:
197            proj="merc"
198    elif typefile in ['gcm']:  proj="cyl"    ## pb avec les autres (de trace derriere la sphere ?)
199    else:                      proj="ortho"
200    return proj   
201
202def ptitle (name):
203    from matplotlib.pyplot import title
204    title(name)
205    print name
206
207def polarinterv (lon2d,lat2d):
208    import numpy as np
209    wlon = [np.min(lon2d),np.max(lon2d)]
210    ind = np.array(lat2d).shape[0] / 2  ## to get a good boundlat and to get the pole
211    wlat = [np.min(lat2d[ind,:]),np.max(lat2d[ind,:])]
212    return [wlon,wlat]
213
214def simplinterv (lon2d,lat2d):
215    import numpy as np
216    return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]]
217
218def wrfinterv (lon2d,lat2d):
219    nx = len(lon2d[0,:])-1
220    ny = len(lon2d[:,0])-1
221    lon1 = lon2d[0,0] 
222    lon2 = lon2d[nx,ny] 
223    lat1 = lat2d[0,0] 
224    lat2 = lat2d[nx,ny] 
225    if abs(0.5*(lat1+lat2)) > 60.:  wider = 0.5 * (abs(lon1)+abs(lon2)) * 0.1
226    else:                           wider = 0.
227    if lon1 < lon2:  wlon = [lon1, lon2 + wider] 
228    else:            wlon = [lon2, lon1 + wider]
229    if lat1 < lat2:  wlat = [lat1, lat2]
230    else:            wlat = [lat2, lat1]
231    return [wlon,wlat]
232
233def makeplotres (filename,res=None,pad_inches_value=0.25,folder='',disp=True,ext='png',erase=False):
234    import  matplotlib.pyplot as plt
235    from os import system
236    addstr = ""
237    if res is not None:
238        res = int(res)
239        addstr = "_"+str(res)
240    name = filename+addstr+"."+ext
241    if folder != '':      name = folder+'/'+name
242    plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value)
243    if disp:              display(name)
244    if ext in ['eps','ps','svg']:   system("tar czvf "+name+".tar.gz "+name+" ; rm -f "+name)
245    if erase:   system("mv "+name+" to_be_erased")             
246    return
247
248def dumpbdy (field,n,stag=None):
249    nx = len(field[0,:])-1
250    ny = len(field[:,0])-1
251    if stag == 'U': nx = nx-1
252    if stag == 'V': ny = ny-1
253    return field[n:ny-n,n:nx-n]
254
255def getcoorddef ( nc ):   
256    ## getcoord2d for predefined types
257    typefile = whatkindfile(nc)
258    if typefile in ['mesoapi','meso']:
259        [lon2d,lat2d] = getcoord2d(nc)
260        lon2d = dumpbdy(lon2d,6)
261        lat2d = dumpbdy(lat2d,6)
262    elif typefile in ['gcm','gcmex']:
263        [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
264    elif typefile in ['geo']:
265        [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M')
266    return lon2d,lat2d   
267
268def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False):
269    import numpy as np
270    if is1d:
271        lat = nc.variables[nlat][:]
272        lon = nc.variables[nlon][:]
273        [lon2d,lat2d] = np.meshgrid(lon,lat)
274    else:
275        lat = nc.variables[nlat][0,:,:]
276        lon = nc.variables[nlon][0,:,:]
277        [lon2d,lat2d] = [lon,lat]
278    return lon2d,lat2d
279
280def smooth (field, coeff):
281        ## actually blur_image could work with different coeff on x and y
282        if coeff > 1:   result = blur_image(field,int(coeff))
283        else:           result = field
284        return result
285
286def gauss_kern(size, sizey=None):
287        import numpy as np
288        ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth     
289        # Returns a normalized 2D gauss kernel array for convolutions
290        size = int(size)
291        if not sizey:
292                sizey = size
293        else:
294                sizey = int(sizey)
295        x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
296        g = np.exp(-(x**2/float(size)+y**2/float(sizey)))
297        return g / g.sum()
298
299def blur_image(im, n, ny=None) :
300        from scipy.signal import convolve
301        ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
302        # blurs the image by convolving with a gaussian kernel of typical size n.
303        # The optional keyword argument ny allows for a different size in the y direction.
304        g = gauss_kern(n, sizey=ny)
305        improc = convolve(im, g, mode='same')
306        return improc
307
308def getwinddef (nc):   
309    ## getwinds for predefined types
310    typefile = whatkindfile(nc)
311    ###
312    if typefile is 'mesoapi':    [uchar,vchar] = ['Um','Vm']
313    elif typefile is 'gcm':      [uchar,vchar] = ['u','v']
314    elif typefile is 'meso':     [uchar,vchar] = ['U','V']
315    else:                        [uchar,vchar] = ['not found','not found']
316    ###
317    if typefile in ['meso']:     metwind = False ## geometrical (wrt grid)
318    else:                        metwind = True  ## meteorological (zon/mer)
319    if metwind is False:         print "Not using meteorological winds. You trust numerical grid as being (x,y)"
320    ###
321    return uchar,vchar,metwind
322
323def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True):
324    ## scale regle la reference du vecteur
325    ## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs.
326    import  matplotlib.pyplot               as plt
327    import  numpy                           as np
328    posx = np.min(x) - np.std(x) / 10.
329    posy = np.min(y) - np.std(y) / 10.
330    u = smooth(u,csmooth)
331    v = smooth(v,csmooth)
332    widthvec = 0.003 #0.005 #0.003
333    q = plt.quiver( x[::stride,::stride],\
334                    y[::stride,::stride],\
335                    u[::stride,::stride],\
336                    v[::stride,::stride],\
337                    angles='xy',color=color,pivot='middle',\
338                    scale=factor,width=widthvec )
339    if color in ['white','yellow']:     kcolor='black'
340    else:                               kcolor=color
341    if key: p = plt.quiverkey(q,posx,posy,scale,\
342                   str(int(scale)),coordinates='data',color=kcolor,labelpos='S',labelsep = 0.03)
343    return 
344
345def display (name):
346    from os import system
347    system("display "+name+" > /dev/null 2> /dev/null &")
348    return name
349
350def findstep (wlon):
351    steplon = int((wlon[1]-wlon[0])/4.)  #3
352    step = 120.
353    while step > steplon and step > 15. :       step = step / 2.
354    if step <= 15.:
355        while step > steplon and step > 5.  :   step = step - 5.
356    if step <= 5.:
357        while step > steplon and step > 1.  :   step = step - 1.
358    if step <= 1.:
359        step = 1. 
360    return step
361
362def define_proj (char,wlon,wlat,back=None):
363    from    mpl_toolkits.basemap            import Basemap
364    import  numpy                           as np
365    import  matplotlib                      as mpl
366    from mymath import max
367    meanlon = 0.5*(wlon[0]+wlon[1])
368    meanlat = 0.5*(wlat[0]+wlat[1])
369    if   wlat[0] >= 80.:   blat =  40. 
370    elif wlat[1] <= -80.:  blat = -40.
371    elif wlat[1] >= 0.:    blat = wlat[0] 
372    elif wlat[0] <= 0.:    blat = wlat[1]
373    print "blat ", blat
374    h = 50.  ## en km
375    radius = 3397200.
376    if   char == "cyl":     m = Basemap(rsphere=radius,projection='cyl',\
377                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
378    elif char == "moll":    m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon)
379    elif char == "ortho":   m = Basemap(rsphere=radius,projection='ortho',lon_0=meanlon,lat_0=meanlat)
380    elif char == "lcc":     m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\
381                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
382    elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.)
383    elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=0.)
384    elif char == "nplaea":  m = Basemap(rsphere=radius,projection='nplaea', boundinglat=wlat[0], lon_0=meanlon)
385    elif char == "laea":    m = Basemap(rsphere=radius,projection='laea',lon_0=meanlon,lat_0=meanlat,lat_ts=meanlat,\
386                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
387    elif char == "nsper":   m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.)
388    elif char == "merc":    m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\
389                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
390    fontsizemer = int(mpl.rcParams['font.size']*3./4.)
391    if char in ["cyl","lcc","merc","nsper","laea"]:   step = findstep(wlon)
392    else:                                             step = 10.
393    steplon = step*2.
394    #if back in ["geolocal"]:                         
395    #    step = np.min([5.,step])
396    #    steplon = step
397    print step
398    m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color='grey', fontsize=fontsizemer)
399    m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer)
400    if back: m.warpimage(marsmap(back),scale=0.75)
401            #if not back:
402            #    if not var:                                        back = "mola"    ## if no var:         draw mola
403            #    elif typefile in ['mesoapi','meso','geo'] \
404            #       and proj not in ['merc','lcc','nsper','laea']:  back = "molabw"  ## if var but meso:   draw molabw
405            #    else:                                              pass             ## else:              draw None
406    return m
407
408#### test temporaire
409def putpoints (map,plot):
410    #### from http://www.scipy.org/Cookbook/Matplotlib/Maps
411    # lat/lon coordinates of five cities.
412    lats = [18.4]
413    lons = [-134.0]
414    points=['Olympus Mons']
415    # compute the native map projection coordinates for cities.
416    x,y = map(lons,lats)
417    # plot filled circles at the locations of the cities.
418    map.plot(x,y,'bo')
419    # plot the names of those five cities.
420    wherept = 0 #1000 #50000
421    for name,xpt,ypt in zip(points,x,y):
422       plot.text(xpt+wherept,ypt+wherept,name)
423    ## le nom ne s'affiche pas...
424    return
425
426def calculate_bounds(field,vmin=None,vmax=None):
427    import numpy as np
428    from mymath import max,min,mean
429    ind = np.where(field < 9e+35)
430    fieldcalc = field[ ind ] # la syntaxe compacte ne marche si field est un tuple
431    ###
432    dev = np.std(fieldcalc)*3.0
433    ###
434    if vmin is None:
435        zevmin = mean(fieldcalc) - dev
436    else:             zevmin = vmin
437    ###
438    if vmax is None:  zevmax = mean(fieldcalc) + dev
439    else:             zevmax = vmax
440    if vmin == vmax:
441                      zevmin = mean(fieldcalc) - dev  ### for continuity
442                      zevmax = mean(fieldcalc) + dev  ### for continuity           
443    ###
444    if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0.
445    print "field ", min(fieldcalc), max(fieldcalc)
446    print "bounds ", zevmin, zevmax
447    return zevmin, zevmax
448
449def bounds(what_I_plot,zevmin,zevmax):
450    from mymath import max,min,mean
451    ### might be convenient to add the missing value in arguments
452    #what_I_plot[ what_I_plot < zevmin ] = zevmin#*(1. + 1.e-7)
453    if zevmin < 0: what_I_plot[ what_I_plot < zevmin*(1. - 1.e-7) ] = zevmin*(1. - 1.e-7)
454    else:          what_I_plot[ what_I_plot < zevmin*(1. + 1.e-7) ] = zevmin*(1. + 1.e-7)
455    print "new min ", min(what_I_plot)
456    what_I_plot[ what_I_plot > 9e+35  ] = -9e+35
457    what_I_plot[ what_I_plot > zevmax ] = zevmax
458    print "new max ", max(what_I_plot)
459   
460    return what_I_plot
461
462def nolow(what_I_plot):
463    from mymath import max,min
464    lim = 0.15*0.5*(abs(max(what_I_plot))+abs(min(what_I_plot)))
465    print "on vire en dessous de ", lim
466    what_I_plot [ abs(what_I_plot) < lim ] = 1.e40 
467    return what_I_plot
468
469def zoomset (wlon,wlat,zoom):
470    dlon = abs(wlon[1]-wlon[0])/2.
471    dlat = abs(wlat[1]-wlat[0])/2.
472    [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\
473                    [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ]
474    print "zoom %",zoom,wlon,wlat
475    return wlon,wlat
476
477def fmtvar (whichvar="def"):
478    fmtvar    =     { \
479             "tk":           "%.0f",\
480             "tpot":         "%.0f",\
481             "TSURF":        "%.0f",\
482             "def":          "%.1e",\
483             "PTOT":         "%.0f",\
484             "HGT":          "%.1e",\
485             "USTM":         "%.2f",\
486             "HFX":          "%.0f",\
487             "ICETOT":       "%.1e",\
488             "TAU_ICE":      "%.2f",\
489             "VMR_ICE":      "%.1e",\
490             "MTOT":         "%.0f",\
491             "anomaly":      "%.1f",\
492             "W":            "%.1f",\
493             "WMAX_TH":      "%.1f",\
494             "QSURFICE":     "%.0f",\
495             "Um":           "%.0f",\
496             "ALBBARE":      "%.2f",\
497                    }
498    if whichvar not in fmtvar:
499        whichvar = "def"
500    return fmtvar[whichvar]
501
502####################################################################################################################
503### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
504def defcolorb (whichone="def"):
505    whichcolorb =    { \
506             "def":          "spectral",\
507             "HGT":          "spectral",\
508             "tk":           "gist_heat",\
509             "TSURF":        "RdBu_r",\
510             "QH2O":         "PuBu",\
511             "USTM":         "YlOrRd",\
512             "HFX":          "RdYlBu",\
513             "ICETOT":       "YlGnBu_r",\
514             "MTOT":         "PuBu",\
515             "TAU_ICE":      "Blues",\
516             "VMR_ICE":      "Blues",\
517             "W":            "jet",\
518             "WMAX_TH":      "spectral",\
519             "anomaly":      "RdBu_r",\
520             "QSURFICE":     "hot_r",\
521             "ALBBARE":      "spectral",\
522                     }
523#W --> spectral ou jet
524#spectral BrBG RdBu_r
525    print "predefined colorbars"
526    if whichone not in whichcolorb:
527        whichone = "def"
528    return whichcolorb[whichone]
529
530def definecolorvec (whichone="def"):
531        whichcolor =    { \
532                "def":          "black",\
533                "vis":          "yellow",\
534                "vishires":     "yellow",\
535                "molabw":       "yellow",\
536                "mola":         "black",\
537                "gist_heat":    "white",\
538                "hot":          "tk",\
539                "gist_rainbow": "black",\
540                "spectral":     "black",\
541                "gray":         "red",\
542                "PuBu":         "black",\
543                        }
544        if whichone not in whichcolor:
545                whichone = "def"
546        return whichcolor[whichone]
547
548def marsmap (whichone="vishires"):
549        from os import uname
550        mymachine = uname()[1]
551        ### not sure about speed-up with this method... looks the same
552        if "lmd.jussieu.fr" in mymachine: domain = "/u/aslmd/WWW/maps/"
553        else:                             domain = "http://www.lmd.jussieu.fr/~aslmd/maps/"
554        whichlink =     { \
555                #"vis":         "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\
556                #"vishires":    "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\
557                #"geolocal":    "http://dl.dropbox.com/u/11078310/geolocal.jpg",\
558                #"mola":        "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\
559                #"molabw":      "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\
560                "vis":         domain+"mar0kuu2.jpg",\
561                "vishires":    domain+"MarsMap_2500x1250.jpg",\
562                "geolocal":    domain+"geolocal.jpg",\
563                "mola":        domain+"mars-mola-2k.jpg",\
564                "molabw":      domain+"MarsElevation_2500x1250.jpg",\
565                "clouds":      "http://www.johnstonsarchive.net/spaceart/marswcloudmap.jpg",\
566                "jupiter":     "http://www.mmedia.is/~bjj/data/jupiter_css/jupiter_css.jpg",\
567                "jupiter_voy": "http://www.mmedia.is/~bjj/data/jupiter/jupiter_vgr2.jpg",\
568                "bw":          "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg",\
569                "contrast":    "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg",\
570                "nice":        "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg",\
571                "blue":        "http://eoimages.gsfc.nasa.gov/ve/2430/land_ocean_ice_2048.jpg",\
572                "blueclouds":  "http://eoimages.gsfc.nasa.gov/ve/2431/land_ocean_ice_cloud_2048.jpg",\
573                "justclouds":  "http://eoimages.gsfc.nasa.gov/ve/2432/cloud_combined_2048.jpg",\
574                        }
575        ### see http://www.mmedia.is/~bjj/planetary_maps.html
576        if whichone not in whichlink: 
577                print "marsmap: choice not defined... you'll get the default one... "
578                whichone = "vishires" 
579        return whichlink[whichone]
580
581#def earthmap (whichone):
582#       if   whichone == "contrast":    whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg"
583#       elif whichone == "bw":          whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg"
584#       elif whichone == "nice":        whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg"
585#       return whichlink
586
587def latinterv (area="Whole"):
588    list =    { \
589        "Europe":                [[ 20., 80.],[- 50.,  50.]],\
590        "Central_America":       [[-10., 40.],[ 230., 300.]],\
591        "Africa":                [[-20., 50.],[- 50.,  50.]],\
592        "Whole":                 [[-90., 90.],[-180., 180.]],\
593        "Southern_Hemisphere":   [[-90., 60.],[-180., 180.]],\
594        "Northern_Hemisphere":   [[-60., 90.],[-180., 180.]],\
595        "Tharsis":               [[-30., 60.],[-170.,- 10.]],\
596        "Whole_No_High":         [[-60., 60.],[-180., 180.]],\
597        "Chryse":                [[-60., 60.],[- 60.,  60.]],\
598        "North_Pole":            [[ 50., 90.],[-180., 180.]],\
599        "Close_North_Pole":      [[ 75., 90.],[-180., 180.]],\
600        "Far_South_Pole":        [[-90.,-40.],[-180., 180.]],\
601        "South_Pole":            [[-90.,-50.],[-180., 180.]],\
602        "Close_South_Pole":      [[-90.,-75.],[-180., 180.]],\
603              }
604    if area not in list:   area = "Whole"
605    [olat,olon] = list[area]
606    return olon,olat
607
Note: See TracBrowser for help on using the repository browser.