source: trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/myplot.py @ 251

Last change on this file since 251 was 247, checked in by aslmd, 14 years ago

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