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

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

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