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

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

Added new routine to our python subroutines pool: make_netcdf.py

This routine is an easy way to create a netcdf file conform to the gcm .nc format.
Here is an exemple of a call to the routine:

#####################

from make_netcdf import make_gcm_netcdf

make_gcm_netcdf (zfilename="diagfi_TES.nc", \

zdescription="TES day and night temperature fields", \
zlon=lon, \
zlat=lat, \
zalt=alt, \
ztime=time, \
zvariables=[ps_day,ps_night,aps,bps,temp_day,temp_night], \
znames=["ps_day","ps_night","aps","bps","temp_day","temp_night"])

####################

This will create a file diagfi_TES.nc with variables specified in "zvariables", with names specified in "znames", along the dimensions given in "zlon,zlat,zalt,ztime". One can use less dimensions if the variables to be written are not 4D.

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