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

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

PYTHON
======

New options:

i) --mark lon,lat will superimpose a cross at lon,lat (floats) on a longitude latitude plot (mapmode 1)
ii) --lstyle let you customize linestyles for 1D plots. Works in the same way as --labels. Exemple: --lstyle "--r","-r","--b","-b"
iii) One can now ask for the variable -v slopexy or SLOPEXY to plot (or overplot with --var2 slopexy) the amplitude of the slope in mesoscale simulations
iv) One can now ask for the variable -v deltat or DELTAT to plot the temperature difference between surface temperature and an arbitrary model level. This automatically calls API with "tk,tsurf". One should use this in conjunction with -i 4 -l xxx so that the difference is made between 0m (surface) and xxx m.

Minor corrections:

Verification plot showing longitude latitude of a slice in a HGT map now also requires "display" to be true.
Winds can now be overplotted in --operation - maps.
File size: 52.2 KB
Line 
1## Author: AS
2def errormess(text,printvar=None):
3    print text
4    if printvar is not None: 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,var2=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    if var2:              basename = basename + '_' + var2
29    return basename
30
31## Author: AS
32def localtime(utc,lon):
33    ltst = utc + lon / 15.
34    ltst = int (ltst * 10) / 10.
35    ltst = ltst % 24
36    return ltst
37
38## Author: AC
39def check_localtime(time):
40    a=-1
41    for i in range(len(time)-1):
42       if (time[i] > time[i+1]): a=i
43    if a >= 0 and a < (len(time)-1)/2.:
44       print "Sorry, time axis is not regular."
45       print "Contourf needs regular axis... recasting"
46       for i in range(a+1):
47          time[i]=time[i]-24.
48    if a >= 0 and a >= (len(time)-1)/2.:
49       print "Sorry, time axis is not regular."
50       print "Contourf needs regular axis... recasting"
51       for i in range((len(time)-1) - a):
52          time[a+1+i]=time[a+1+i]+24.
53    return time
54
55## Author: AS, AC, JL
56def whatkindfile (nc):
57    if 'controle' in nc.variables:             typefile = 'gcm'
58    elif 'phisinit' in nc.variables:           typefile = 'gcm'
59    elif 'time_counter' in nc.variables:       typefile = 'earthgcm'
60    elif hasattr(nc,'START_DATE'):             typefile = 'meso' 
61    elif 'HGT_M' in nc.variables:              typefile = 'geo'
62    elif hasattr(nc,'institution'):
63      if "European Centre" in getattr(nc,'institution'):  typefile = 'ecmwf'
64    else:                                      typefile = 'gcm' # for lslin-ed files from gcm
65    return typefile
66
67## Author: AS
68def getfield (nc,var):
69    ## this allows to get much faster (than simply referring to nc.variables[var])
70    import numpy as np
71    dimension = len(nc.variables[var].dimensions)
72    #print "   Opening variable with", dimension, "dimensions ..."
73    if dimension == 2:    field = nc.variables[var][:,:]
74    elif dimension == 3:  field = nc.variables[var][:,:,:]
75    elif dimension == 4:  field = nc.variables[var][:,:,:,:]
76    elif dimension == 1:  field = nc.variables[var][:]
77    # if there are NaNs in the ncdf, they should be loaded as a masked array which will be
78    # recasted as a regular array later in reducefield
79    if (np.isnan(np.sum(field)) and (type(field).__name__ not in 'MaskedArray')):
80       print "Warning: netcdf as nan values but is not loaded as a Masked Array."
81       print "recasting array type"
82       out=np.ma.masked_invalid(field)
83       out.set_fill_value([np.NaN])
84    else:
85    # missing values from zrecast or hrecast are -1e-33
86       masked=np.ma.masked_where(field < -1e30,field)
87       masked2=np.ma.masked_where(field > 1e35,field)
88       masked.set_fill_value([np.NaN]) ; masked2.set_fill_value([np.NaN])
89       mask = np.ma.getmask(masked) ; mask2 = np.ma.getmask(masked2)
90       if (True in np.array(mask)):
91          out=masked
92          print "Masked array... Missing value is NaN"
93       elif (True in np.array(mask2)):
94          out=masked2
95          print "Masked array... Missing value is NaN"
96#       else:
97#       # missing values from api are 1e36
98#          masked=np.ma.masked_where(field > 1e35,field)
99#          masked.set_fill_value([np.NaN])
100#          mask = np.ma.getmask(masked)
101#          if (True in np.array(mask)):out=masked
102#          else:out=field
103       else:out=field
104    return out
105
106## Author: AC
107# Compute the norm of the winds
108# The corresponding variable to call is UV or uvmet (to use api)
109def windamplitude (nc):
110    import numpy as np
111    varinfile = nc.variables.keys()
112    if "U" in varinfile: zu=getfield(nc,'U')
113    elif "Um" in varinfile: zu=getfield(nc,'Um')
114    if "V" in varinfile: zv=getfield(nc,'V')
115    elif "Vm" in varinfile: zv=getfield(nc,'Vm')
116    znt,znz,zny,znx = np.array(zu).shape
117    if "U" in varinfile:znx=znx-1
118    zuint = np.zeros([znt,znz,zny,znx])
119    zvint = np.zeros([znt,znz,zny,znx])
120    if "U" in varinfile:
121       for xx in np.arange(znx):      zuint[:,:,:,xx] = (zu[:,:,:,xx] + zu[:,:,:,xx+1])/2.
122       for yy in np.arange(zny):      zvint[:,:,yy,:] = (zv[:,:,yy,:] + zv[:,:,yy+1,:])/2.
123    else:
124       zuint=zu
125       zvint=zv
126    return np.sqrt(zuint**2 + zvint**2)
127
128## Author: AC
129# Compute the norm of the slope angles
130# The corresponding variable to call is SLOPEXY
131def slopeamplitude (nc):
132    import numpy as np
133    varinfile = nc.variables.keys()
134    if "slopex" in varinfile: zu=getfield(nc,'slopex')
135    elif "SLOPEX" in varinfile: zu=getfield(nc,'SLOPEX')
136    if "slopey" in varinfile: zv=getfield(nc,'slopey')
137    elif "SLOPEY" in varinfile: zv=getfield(nc,'SLOPEY')
138    znt,zny,znx = np.array(zu).shape
139    zuint = np.zeros([znt,zny,znx])
140    zvint = np.zeros([znt,zny,znx])
141    zuint=zu
142    zvint=zv
143    return np.sqrt(zuint**2 + zvint**2)
144
145## Author: AC
146# Compute the temperature difference between surface and first level.
147# API is automatically called to get TSURF and TK.
148# The corresponding variable to call is DELTAT
149def deltat0t1 (nc):
150    import numpy as np
151    varinfile = nc.variables.keys()
152    if "tsurf" in varinfile: zu=getfield(nc,'tsurf')
153    elif "TSURF" in varinfile: zu=getfield(nc,'TSURF')
154    if "tk" in varinfile: zv=getfield(nc,'tk')
155    elif "TK" in varinfile: zv=getfield(nc,'TK')
156    znt,zny,znx = np.array(zu).shape
157    zuint = np.zeros([znt,zny,znx])
158    zuint=zu - zv[:,0,:,:]
159    return zuint
160
161## Author: AS + TN + AC
162def reducefield (input,d4=None,d3=None,d2=None,d1=None,yint=False,alt=None,anomaly=False,redope=None,mesharea=None):
163    ### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x"
164    ### it would be actually better to name d4 d3 d2 d1 as t z y x
165    ### ... note, anomaly is only computed over d1 and d2 for the moment
166    import numpy as np
167    from mymath import max,mean,min,sum
168    csmooth = 12 ## a fair amount of grid points (too high results in high computation time)
169    if redope is not None:
170       if   redope == "mint":     input = min(input,axis=0) ; d1 = None
171       elif redope == "maxt":     input = max(input,axis=0) ; d1 = None
172       else:                      errormess("not supported. but try lines in reducefield beforehand.")
173       #elif redope == "minz":     input = min(input,axis=1) ; d2 = None
174       #elif redope == "maxz":     input = max(input,axis=1) ; d2 = None
175       #elif redope == "miny":     input = min(input,axis=2) ; d3 = None
176       #elif redope == "maxy":     input = max(input,axis=2) ; d3 = None
177       #elif redope == "minx":     input = min(input,axis=3) ; d4 = None
178       #elif redope == "maxx":     input = max(input,axis=3) ; d4 = None
179    dimension = np.array(input).ndim
180    shape = np.array(np.array(input).shape)
181    #print 'd1,d2,d3,d4: ',d1,d2,d3,d4
182    if anomaly: print 'ANOMALY ANOMALY'
183    output = input
184    error = False
185    #### this is needed to cope the case where d4,d3,d2,d1 are single integers and not arrays
186    if d4 is not None and not isinstance(d4, np.ndarray): d4=[d4]
187    if d3 is not None and not isinstance(d3, np.ndarray): d3=[d3]
188    if d2 is not None and not isinstance(d2, np.ndarray): d2=[d2]
189    if d1 is not None and not isinstance(d1, np.ndarray): d1=[d1]
190    ### now the main part
191    if dimension == 2:
192        if mesharea is None: mesharea=np.ones(shape)
193        if   max(d2) >= shape[0]: error = True
194        elif max(d1) >= shape[1]: error = True
195        elif d1 is not None and d2 is not None:
196          totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0)
197          output = output*mesharea; output = sum(output[d2,:],axis=0); output = sum(output[d1],axis=0)/totalarea
198        elif d1 is not None:         output = mean(input[:,d1],axis=1)
199        elif d2 is not None:         totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[d2,:],axis=0)/totalarea
200    elif dimension == 3:
201        if mesharea is None: mesharea=np.ones(shape[[1,2]])
202        if   max(d4) >= shape[0]: error = True
203        elif max(d2) >= shape[1]: error = True
204        elif max(d1) >= shape[2]: error = True
205        elif d4 is not None and d2 is not None and d1 is not None: 
206          totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0)
207          output = mean(input[d4,:,:],axis=0)
208          output = output*mesharea; output = sum(output[d2,:],axis=0); output = sum(output[d1],axis=0)/totalarea
209        elif d4 is not None and d2 is not None:
210          output = mean(input[d4,:,:],axis=0)
211          totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[d2,:],axis=0)/totalarea
212        elif d4 is not None and d1 is not None:    output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1)
213        elif d2 is not None and d1 is not None:
214          totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0)
215          output = output*mesharea; output = sum(output[:,d2,:],axis=1); output = sum(output[:,d1],axis=1)/totalarea
216        elif d1 is not None:   output = mean(input[:,:,d1],axis=2)
217        elif d2 is not None:   totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,d2,:],axis=1)/totalarea
218        elif d4 is not None:   output = mean(input[d4,:,:],axis=0)
219    elif dimension == 4:
220        if mesharea is None: mesharea=np.ones(shape[[2,3]])
221        if   max(d4) >= shape[0]: error = True
222        elif max(d3) >= shape[1]: error = True
223        elif max(d2) >= shape[2]: error = True
224        elif max(d1) >= shape[3]: error = True
225        elif d4 is not None and d3 is not None and d2 is not None and d1 is not None:
226             output = mean(input[d4,:,:,:],axis=0)
227             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
228             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
229             totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0)
230             output = output*mesharea
231             output = sum(output[d2,:],axis=0)
232             output = sum(output[d1],axis=0)/totalarea
233        elif d4 is not None and d3 is not None and d2 is not None: 
234             output = mean(input[d4,:,:,:],axis=0)
235             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
236             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
237             totalarea=sum(mesharea[d2,:],axis=0)
238             output = output*mesharea; output = sum(output[d2,:],axis=0)/totalarea
239        elif d4 is not None and d3 is not None and d1 is not None: 
240             output = mean(input[d4,:,:,:],axis=0)
241             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
242             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) 
243             output = mean(output[:,d1],axis=1)
244        elif d4 is not None and d2 is not None and d1 is not None: 
245             output = mean(input[d4,:,:,:],axis=0)
246             if anomaly:
247                 for k in range(output.shape[0]):  output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
248             totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0)
249             output = output*mesharea; output = sum(output[:,d2,:],axis=1); output = sum(output[:,d1],axis=1)/totalarea
250             #noperturb = smooth1d(output,window_len=7)
251             #lenlen = len(output) ; output = output[1:lenlen-7] ; yeye = noperturb[4:lenlen-4]
252             #plot(output) ; plot(yeye) ; show() ; plot(output-yeye) ; show()
253        elif d3 is not None and d2 is not None and d1 is not None: 
254             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3) 
255             if anomaly:
256                 for k in range(output.shape[0]):  output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
257             totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0)
258             output = output*mesharea; output = sum(output[:,d2,:],axis=1); output = sum(output[:,d1],axis=1)/totalarea
259        elif d4 is not None and d3 is not None: 
260             output = mean(input[d4,:,:,:],axis=0) 
261             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
262             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) 
263        elif d4 is not None and d2 is not None: 
264             output = mean(input[d4,:,:,:],axis=0) 
265             totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,d2,:],axis=1)/totalarea
266        elif d4 is not None and d1 is not None: 
267             output = mean(input[d4,:,:,:],axis=0) 
268             output = mean(output[:,:,d1],axis=2)
269        elif d3 is not None and d2 is not None: 
270             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
271             totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,d2,:],axis=1)/totalarea
272        elif d3 is not None and d1 is not None: 
273             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
274             output = mean(output[:,:,d1],axis=2)
275        elif d2 is not None and d1 is not None: 
276             totalarea=sum(mesharea[d2,:],axis=0); totalarea=sum(totalarea[d1],axis=0)
277             output = output*mesharea; output = sum(output[:,:,d2,:],axis=2); output = sum(output[:,:,d1],axis=2)/totalarea
278        elif d1 is not None:        output = mean(input[:,:,:,d1],axis=3)
279        elif d2 is not None: 
280             totalarea=sum(mesharea[d2,:],axis=0); output = output*mesharea; output = sum(output[:,:,d2,:],axis=2)/totalarea
281        elif d3 is not None:        output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
282        elif d4 is not None:        output = mean(input[d4,:,:,:],axis=0)
283    dimension2 = np.array(output).ndim
284    shape2 = np.array(output).shape
285    print 'REDUCEFIELD dim,shape: ',dimension,shape,' >>> ',dimension2,shape2
286    return output, error
287
288## Author: AC + AS
289def reduce_zaxis (input,ax=None,yint=False,vert=None,indice=None):
290    from mymath import max,mean
291    from scipy import integrate
292    if yint and vert is not None and indice is not None:
293       if type(input).__name__=='MaskedArray':
294          input.set_fill_value([np.NaN])
295          output = integrate.trapz(input.filled(),x=vert[indice],axis=ax)
296       else:
297          output = integrate.trapz(input,x=vert[indice],axis=ax)
298    else:
299       output = mean(input,axis=ax)
300    return output
301
302## Author: AS + TN
303def definesubplot ( numplot, fig ):
304    from matplotlib.pyplot import rcParams
305    rcParams['font.size'] = 12. ## default (important for multiple calls)
306    if numplot <= 0:
307        subv = 99999
308        subh = 99999
309    elif numplot == 1: 
310        subv = 1
311        subh = 1
312    elif numplot == 2:
313        subv = 1 #2
314        subh = 2 #1
315        fig.subplots_adjust(wspace = 0.35)
316        rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
317    elif numplot == 3:
318        subv = 3
319        subh = 1
320        fig.subplots_adjust(wspace = 0.5)
321        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
322    elif   numplot == 4:
323        subv = 2
324        subh = 2
325        fig.subplots_adjust(wspace = 0.4, hspace = 0.3)
326        rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
327    elif numplot <= 6:
328        subv = 2
329        subh = 3
330        fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
331        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
332    elif numplot <= 8:
333        subv = 2
334        subh = 4
335        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
336        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
337    elif numplot <= 9:
338        subv = 3
339        subh = 3
340        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
341        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
342    elif numplot <= 12:
343        subv = 3
344        subh = 4
345        fig.subplots_adjust(wspace = 0, hspace = 0.1)
346        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
347    elif numplot <= 16:
348        subv = 4
349        subh = 4
350        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
351        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
352    else:
353        print "number of plot supported: 1 to 16"
354        exit()
355    return subv,subh
356
357## Author: AS
358def getstralt(nc,nvert):
359    varinfile = nc.variables.keys()
360    if 'vert' not in varinfile:
361        stralt = "_lvl" + str(nvert)
362    else:
363        zelevel = int(nc.variables['vert'][nvert])
364        if abs(zelevel) < 10000.:   strheight=str(zelevel)+"m"
365        else:                       strheight=str(int(zelevel/1000.))+"km"
366        if 'altitude'       in nc.dimensions:   stralt = "_"+strheight+"-AMR"
367        elif 'altitude_abg' in nc.dimensions:   stralt = "_"+strheight+"-ALS"
368        elif 'bottom_top'   in nc.dimensions:   stralt = "_"+strheight
369        elif 'pressure'     in nc.dimensions:   stralt = "_"+str(zelevel)+"Pa"
370        else:                                   stralt = ""
371    return stralt
372
373## Author: AS
374def getlschar ( namefile, getaxis=False ):
375    from netCDF4 import Dataset
376    from timestuff import sol2ls
377    from numpy import array
378    from string import rstrip
379    namefile = rstrip( rstrip( rstrip( namefile, chars="_z"), chars="_zabg"), chars="_p")
380    #### we assume that wrfout is next to wrfout_z and wrfout_zabg
381    nc  = Dataset(namefile)
382    zetime = None
383    days_in_month = [61, 66, 66, 65, 60, 54, 50, 46, 47, 47, 51, 56]
384    plus_in_month = [ 0, 61,127,193,258,318,372,422,468,515,562,613]
385    if 'Times' in nc.variables:
386        zetime = nc.variables['Times'][0]
387        shape = array(nc.variables['Times']).shape
388        if shape[0] < 2: zetime = None
389    if zetime is not None \
390       and 'vert' not in nc.variables:
391        ##### strangely enough this does not work for api or ncrcat results!
392        zesol = plus_in_month[ int(zetime[5]+zetime[6])-1 ] + int(zetime[8]+zetime[9]) - 1 ##les sols GCM commencent a 0
393        dals = int( 10. * sol2ls ( zesol ) ) / 10.
394        ###
395        zetime2 = nc.variables['Times'][1]
396        one  = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37.
397        next = int(zetime2[11]+zetime2[12]) + int(zetime2[14]+zetime2[15])/37. 
398        zehour    = one
399        zehourin  = abs ( next - one )
400        if not getaxis:
401            lschar = "_Ls"+str(dals)
402        else:
403            zelen = len(nc.variables['Times'][:])
404            yeye = range(zelen) ; lsaxis = range(zelen) ; solaxis = range(zelen) ; ltaxis = range(zelen)
405            for iii in yeye:
406               zetime = nc.variables['Times'][iii] 
407               ltaxis[iii] = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37.
408               solaxis[iii] = ltaxis[iii] / 24. + plus_in_month[ int(zetime[5]+zetime[6])-1 ] + int(zetime[8]+zetime[9]) - 1 ##les sols GCM commencent a 0
409               lsaxis[iii] = sol2ls ( solaxis[iii] )
410               if ltaxis[iii] < ltaxis[iii-1]: ltaxis[iii] = ltaxis[iii] + 24.
411               #print ltaxis[iii], solaxis[iii], lsaxis[iii], getattr( nc, 'JULDAY' )
412            lschar = lsaxis ; zehour = solaxis ; zehourin = ltaxis
413    else:
414        lschar=""
415        zehour = 0
416        zehourin = 1 
417    return lschar, zehour, zehourin
418
419## Author: AS
420def getprefix (nc):
421    prefix = 'LMD_MMM_'
422    prefix = prefix + 'd'+str(getattr(nc,'GRID_ID'))+'_'
423    prefix = prefix + str(int(getattr(nc,'DX')/1000.))+'km_'
424    return prefix
425
426## Author: AS
427def getproj (nc):
428    typefile = whatkindfile(nc)
429    if typefile in ['meso','geo']:
430        ### (il faudrait passer CEN_LON dans la projection ?)
431        map_proj = getattr(nc, 'MAP_PROJ')
432        cen_lat  = getattr(nc, 'CEN_LAT')
433        if map_proj == 2:
434            if cen_lat > 10.:   
435                proj="npstere"
436                #print "NP stereographic polar domain"
437            else:           
438                proj="spstere"
439                #print "SP stereographic polar domain"
440        elif map_proj == 1: 
441            #print "lambert projection domain"
442            proj="lcc"
443        elif map_proj == 3: 
444            #print "mercator projection"
445            proj="merc"
446        else:
447            proj="merc"
448    elif typefile in ['gcm']:        proj="cyl"    ## pb avec les autres (de trace derriere la sphere ?)
449    else:                            proj="ortho"
450    return proj   
451
452## Author: AS
453def ptitle (name):
454    from matplotlib.pyplot import title
455    title(name)
456    print name
457
458## Author: AS
459def polarinterv (lon2d,lat2d):
460    import numpy as np
461    wlon = [np.min(lon2d),np.max(lon2d)]
462    ind = np.array(lat2d).shape[0] / 2  ## to get a good boundlat and to get the pole
463    wlat = [np.min(lat2d[ind,:]),np.max(lat2d[ind,:])]
464    return [wlon,wlat]
465
466## Author: AS
467def simplinterv (lon2d,lat2d):
468    import numpy as np
469    return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]]
470
471## Author: AS
472def wrfinterv (lon2d,lat2d):
473    nx = len(lon2d[0,:])-1
474    ny = len(lon2d[:,0])-1
475    lon1 = lon2d[0,0] 
476    lon2 = lon2d[nx,ny] 
477    lat1 = lat2d[0,0] 
478    lat2 = lat2d[nx,ny] 
479    if abs(0.5*(lat1+lat2)) > 60.:  wider = 0.5 * (abs(lon1)+abs(lon2)) * 0.1
480    else:                           wider = 0.
481    if lon1 < lon2:  wlon = [lon1, lon2 + wider] 
482    else:            wlon = [lon2, lon1 + wider]
483    if lat1 < lat2:  wlat = [lat1, lat2]
484    else:            wlat = [lat2, lat1]
485    return [wlon,wlat]
486
487## Author: AS
488def makeplotres (filename,res=None,pad_inches_value=0.25,folder='',disp=True,ext='png',erase=False):
489    import  matplotlib.pyplot as plt
490    from os import system
491    addstr = ""
492    if res is not None:
493        res = int(res)
494        addstr = "_"+str(res)
495    name = filename+addstr+"."+ext
496    if folder != '':      name = folder+'/'+name
497    plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value)
498    if disp:              display(name)
499    if ext in ['eps','ps','svg']:   system("tar czvf "+name+".tar.gz "+name+" ; rm -f "+name)
500    if erase:   system("mv "+name+" to_be_erased")             
501    return
502
503## Author: AS + AC
504def dumpbdy (field,n,stag=None,condition=False,onlyx=False,onlyy=False):
505    nx = len(field[0,:])-1
506    ny = len(field[:,0])-1
507    if condition:
508      if stag == 'U': nx = nx-1
509      if stag == 'V': ny = ny-1
510      if stag == 'W': nx = nx+1 #special les case when we dump stag on W
511    if onlyx:    result = field[:,n:nx-n]
512    elif onlyy:  result = field[n:ny-n,:]
513    else:        result = field[n:ny-n,n:nx-n]
514    return result
515
516## Author: AS + AC
517def getcoorddef ( nc ):   
518    import numpy as np
519    ## getcoord2d for predefined types
520    typefile = whatkindfile(nc)
521    if typefile in ['meso']:
522        if '9999' not in getattr(nc,'START_DATE') :   
523            ## regular mesoscale 
524            [lon2d,lat2d] = getcoord2d(nc) 
525        else:                     
526            ## idealized mesoscale                     
527            nx=getattr(nc,'WEST-EAST_GRID_DIMENSION')
528            ny=getattr(nc,'SOUTH-NORTH_GRID_DIMENSION')
529            dlat=getattr(nc,'DX')
530            ## this is dirty because Martian-specific
531            # ... but this just intended to get "lat-lon" like info
532            falselon = np.arange(-dlat*(nx-1)/2.,dlat*(nx-1)/2.,dlat)/60000.
533            falselat = np.arange(-dlat*(ny-1)/2.,dlat*(ny-1)/2.,dlat)/60000.
534            [lon2d,lat2d] = np.meshgrid(falselon,falselat) ## dummy coordinates
535            print "WARNING: domain plot artificially centered on lat,lon 0,0"
536    elif typefile in ['gcm']: 
537        [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
538    elif typefile in ['earthgcm','ecmwf']: 
539        [lon2d,lat2d] = getcoord2d(nc,nlat="lat",nlon="lon",is1d=True)
540    elif typefile in ['geo']:
541        [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M')
542    return lon2d,lat2d   
543
544## Author: AS
545def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False):
546    import numpy as np
547    if is1d:
548        lat = nc.variables[nlat][:]
549        lon = nc.variables[nlon][:]
550        [lon2d,lat2d] = np.meshgrid(lon,lat)
551    else:
552        lat = nc.variables[nlat][0,:,:]
553        lon = nc.variables[nlon][0,:,:]
554        [lon2d,lat2d] = [lon,lat]
555    return lon2d,lat2d
556
557## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
558def smooth1d(x,window_len=11,window='hanning'):
559    import numpy
560    """smooth the data using a window with requested size.
561    This method is based on the convolution of a scaled window with the signal.
562    The signal is prepared by introducing reflected copies of the signal
563    (with the window size) in both ends so that transient parts are minimized
564    in the begining and end part of the output signal.
565    input:
566        x: the input signal
567        window_len: the dimension of the smoothing window; should be an odd integer
568        window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
569            flat window will produce a moving average smoothing.
570    output:
571        the smoothed signal
572    example:
573    t=linspace(-2,2,0.1)
574    x=sin(t)+randn(len(t))*0.1
575    y=smooth(x)
576    see also:
577    numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
578    scipy.signal.lfilter
579    TODO: the window parameter could be the window itself if an array instead of a string   
580    """
581    x = numpy.array(x)
582    if x.ndim != 1:
583        raise ValueError, "smooth only accepts 1 dimension arrays."
584    if x.size < window_len:
585        raise ValueError, "Input vector needs to be bigger than window size."
586    if window_len<3:
587        return x
588    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
589        raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
590    s=numpy.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
591    #print(len(s))
592    if window == 'flat': #moving average
593        w=numpy.ones(window_len,'d')
594    else:
595        w=eval('numpy.'+window+'(window_len)')
596    y=numpy.convolve(w/w.sum(),s,mode='valid')
597    return y
598
599## Author: AS
600def smooth (field, coeff):
601        ## actually blur_image could work with different coeff on x and y
602        if coeff > 1:   result = blur_image(field,int(coeff))
603        else:           result = field
604        return result
605
606## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
607def gauss_kern(size, sizey=None):
608        import numpy as np
609        # Returns a normalized 2D gauss kernel array for convolutions
610        size = int(size)
611        if not sizey:
612                sizey = size
613        else:
614                sizey = int(sizey)
615        x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
616        g = np.exp(-(x**2/float(size)+y**2/float(sizey)))
617        return g / g.sum()
618
619## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
620def blur_image(im, n, ny=None) :
621        from scipy.signal import convolve
622        # blurs the image by convolving with a gaussian kernel of typical size n.
623        # The optional keyword argument ny allows for a different size in the y direction.
624        g = gauss_kern(n, sizey=ny)
625        improc = convolve(im, g, mode='same')
626        return improc
627
628## Author: AS
629def getwinddef (nc):   
630    ###
631    varinfile = nc.variables.keys()
632    if 'Um' in varinfile:   [uchar,vchar] = ['Um','Vm'] #; print "this is API meso file"
633    elif 'U' in varinfile:  [uchar,vchar] = ['U','V']   #; print "this is RAW meso file"
634    elif 'u' in varinfile:  [uchar,vchar] = ['u','v']   #; print "this is GCM file"
635     ### you can add choices here !
636    else:                   [uchar,vchar] = ['not found','not found']
637    ###
638    if uchar in ['U']:         metwind = False ## geometrical (wrt grid)
639    else:                      metwind = True  ## meteorological (zon/mer)
640    if metwind is False:       print "Not using meteorological winds. You trust numerical grid as being (x,y)"
641    ###
642    return uchar,vchar,metwind
643
644## Author: AS
645def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True):
646    ## scale regle la reference du vecteur
647    ## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs.
648    import  matplotlib.pyplot               as plt
649    import  numpy                           as np
650    posx = np.min(x) - np.std(x) / 10.
651    posy = np.min(y) - np.std(y) / 10.
652    u = smooth(u,csmooth)
653    v = smooth(v,csmooth)
654    widthvec = 0.003 #0.005 #0.003
655    q = plt.quiver( x[::stride,::stride],\
656                    y[::stride,::stride],\
657                    u[::stride,::stride],\
658                    v[::stride,::stride],\
659                    angles='xy',color=color,pivot='middle',\
660                    scale=factor,width=widthvec )
661    if color in ['white','yellow']:     kcolor='black'
662    else:                               kcolor=color
663    if key: p = plt.quiverkey(q,posx,posy,scale,\
664                   str(int(scale)),coordinates='data',color=kcolor,labelpos='S',labelsep = 0.03)
665    return 
666
667## Author: AS
668def display (name):
669    from os import system
670    system("display "+name+" > /dev/null 2> /dev/null &")
671    return name
672
673## Author: AS
674def findstep (wlon):
675    steplon = int((wlon[1]-wlon[0])/4.)  #3
676    step = 120.
677    while step > steplon and step > 15. :       step = step / 2.
678    if step <= 15.:
679        while step > steplon and step > 5.  :   step = step - 5.
680    if step <= 5.:
681        while step > steplon and step > 1.  :   step = step - 1.
682    if step <= 1.:
683        step = 1. 
684    return step
685
686## Author: AS
687def define_proj (char,wlon,wlat,back=None,blat=None,blon=None):
688    from    mpl_toolkits.basemap            import Basemap
689    import  numpy                           as np
690    import  matplotlib                      as mpl
691    from mymath import max
692    meanlon = 0.5*(wlon[0]+wlon[1])
693    meanlat = 0.5*(wlat[0]+wlat[1])
694    if blat is None:
695        ortholat=meanlat
696        if   wlat[0] >= 80.:   blat =  -40. 
697        elif wlat[1] <= -80.:  blat = -40.
698        elif wlat[1] >= 0.:    blat = wlat[0] 
699        elif wlat[0] <= 0.:    blat = wlat[1]
700    else:  ortholat=blat
701    if blon is None:  ortholon=meanlon
702    else:             ortholon=blon
703    h = 50.  ## en km
704    radius = 3397200.
705    if   char == "cyl":     m = Basemap(rsphere=radius,projection='cyl',\
706                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
707    elif char == "moll":    m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon)
708    elif char == "ortho":   m = Basemap(rsphere=radius,projection='ortho',lon_0=ortholon,lat_0=ortholat)
709    elif char == "lcc":     m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\
710                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
711    elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.)
712    elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=180.)
713    elif char == "nplaea":  m = Basemap(rsphere=radius,projection='nplaea', boundinglat=wlat[0], lon_0=meanlon)
714    elif char == "laea":    m = Basemap(rsphere=radius,projection='laea',lon_0=meanlon,lat_0=meanlat,lat_ts=meanlat,\
715                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
716    elif char == "nsper":   m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.)
717    elif char == "merc":    m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\
718                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
719    fontsizemer = int(mpl.rcParams['font.size']*3./4.)
720    if char in ["cyl","lcc","merc","nsper","laea"]:   step = findstep(wlon)
721    else:                                             step = 10.
722    steplon = step*2.
723    zecolor ='grey'
724    zelinewidth = 1
725    zelatmax = 80
726    # to show gcm grid:
727    #zecolor = 'r'
728    #zelinewidth = 1
729    #step = 5.625
730    #steplon = 5.625
731    #zelatmax = 89.9
732    if char not in ["moll"]:
733        m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, latmax=zelatmax)
734        m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, latmax=zelatmax)
735    if back: m.warpimage(marsmap(back),scale=0.75)
736            #if not back:
737            #    if not var:                                        back = "mola"    ## if no var:         draw mola
738            #    elif typefile in ['mesoapi','meso','geo'] \
739            #       and proj not in ['merc','lcc','nsper','laea']:  back = "molabw"  ## if var but meso:   draw molabw
740            #    else:                                              pass             ## else:              draw None
741    return m
742
743## Author: AS
744#### test temporaire
745def putpoints (map,plot):
746    #### from http://www.scipy.org/Cookbook/Matplotlib/Maps
747    # lat/lon coordinates of five cities.
748    lats = [18.4]
749    lons = [-134.0]
750    points=['Olympus Mons']
751    # compute the native map projection coordinates for cities.
752    x,y = map(lons,lats)
753    # plot filled circles at the locations of the cities.
754    map.plot(x,y,'bo')
755    # plot the names of those five cities.
756    wherept = 0 #1000 #50000
757    for name,xpt,ypt in zip(points,x,y):
758       plot.text(xpt+wherept,ypt+wherept,name)
759    ## le nom ne s'affiche pas...
760    return
761
762## Author: AS
763def calculate_bounds(field,vmin=None,vmax=None):
764    import numpy as np
765    from mymath import max,min,mean
766    ind = np.where(field < 9e+35)
767    fieldcalc = field[ ind ] # la syntaxe compacte ne marche si field est un tuple
768    ###
769    dev = np.std(fieldcalc)*3.0
770    ###
771    if vmin is None:  zevmin = mean(fieldcalc) - dev
772    else:             zevmin = vmin
773    ###
774    if vmax is None:  zevmax = mean(fieldcalc) + dev
775    else:             zevmax = vmax
776    if vmin == vmax:
777                      zevmin = mean(fieldcalc) - dev  ### for continuity
778                      zevmax = mean(fieldcalc) + dev  ### for continuity           
779    ###
780    if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0.
781    print "BOUNDS field ", min(fieldcalc), max(fieldcalc), " //// adopted", zevmin, zevmax
782    return zevmin, zevmax
783
784## Author: AS
785def bounds(what_I_plot,zevmin,zevmax):
786    from mymath import max,min,mean
787    ### might be convenient to add the missing value in arguments
788    #what_I_plot[ what_I_plot < zevmin ] = zevmin#*(1. + 1.e-7)
789    if zevmin < 0: what_I_plot[ what_I_plot < zevmin*(1. - 1.e-7) ] = zevmin*(1. - 1.e-7)
790    else:          what_I_plot[ what_I_plot < zevmin*(1. + 1.e-7) ] = zevmin*(1. + 1.e-7)
791    #print "NEW MIN ", min(what_I_plot)
792    what_I_plot[ what_I_plot > 9e+35  ] = -9e+35
793    what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7)
794    #print "NEW MAX ", max(what_I_plot)
795    return what_I_plot
796
797## Author: AS
798def nolow(what_I_plot):
799    from mymath import max,min
800    lim = 0.15*0.5*(abs(max(what_I_plot))+abs(min(what_I_plot)))
801    print "NO PLOT BELOW VALUE ", lim
802    what_I_plot [ abs(what_I_plot) < lim ] = 1.e40 
803    return what_I_plot
804
805
806## Author : AC
807def hole_bounds(what_I_plot,zevmin,zevmax):
808    import numpy as np
809    zi=0
810    for i in what_I_plot:
811        zj=0
812        for j in i:
813            if ((j < zevmin) or (j > zevmax)):what_I_plot[zi,zj]=np.NaN
814            zj=zj+1
815        zi=zi+1
816
817    return what_I_plot
818
819## Author: AS
820def zoomset (wlon,wlat,zoom):
821    dlon = abs(wlon[1]-wlon[0])/2.
822    dlat = abs(wlat[1]-wlat[0])/2.
823    [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\
824                    [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ]
825    print "ZOOM %",zoom,wlon,wlat
826    return wlon,wlat
827
828## Author: AS
829def fmtvar (whichvar="def"):
830    fmtvar    =     { \
831             "MIXED":        "%.0f",\
832             "UPDRAFT":      "%.0f",\
833             "DOWNDRAFT":    "%.0f",\
834             "TK":           "%.0f",\
835             #"ZMAX_TH":      "%.0f",\
836             #"WSTAR":        "%.0f",\
837             # Variables from TES ncdf format
838             "T_NADIR_DAY":  "%.0f",\
839             "T_NADIR_NIT":  "%.0f",\
840             # Variables from tes.py ncdf format
841             "TEMP_DAY":     "%.0f",\
842             "TEMP_NIGHT":   "%.0f",\
843             # Variables from MCS and mcs.py ncdf format
844             "DTEMP":        "%.0f",\
845             "NTEMP":        "%.0f",\
846             "DNUMBINTEMP":  "%.0f",\
847             "NNUMBINTEMP":  "%.0f",\
848             # other stuff
849             "TPOT":         "%.0f",\
850             "TSURF":        "%.0f",\
851             "U_OUT1":       "%.0f",\
852             "T_OUT1":       "%.0f",\
853             "def":          "%.1e",\
854             "PTOT":         "%.0f",\
855             "HGT":          "%.1e",\
856             "USTM":         "%.2f",\
857             "HFX":          "%.0f",\
858             "ICETOT":       "%.1e",\
859             "TAU_ICE":      "%.2f",\
860             "TAUICE":       "%.2f",\
861             "VMR_ICE":      "%.1e",\
862             "MTOT":         "%.1f",\
863             "ANOMALY":      "%.1f",\
864             "W":            "%.1f",\
865             "WMAX_TH":      "%.1f",\
866             "WSTAR":        "%.1f",\
867             "QSURFICE":     "%.0f",\
868             "UM":           "%.0f",\
869             "WIND":         "%.0f",\
870             "UVMET":         "%.0f",\
871             "UV":         "%.0f",\
872             "ALBBARE":      "%.2f",\
873             "TAU":          "%.1f",\
874             "CO2":          "%.2f",\
875             #### T.N.
876             "TEMP":         "%.0f",\
877             "VMR_H2OICE":   "%.0f",\
878             "VMR_H2OVAP":   "%.0f",\
879             "TAUTES":       "%.2f",\
880             "TAUTESAP":     "%.2f",\
881
882                    }
883    if "TSURF" in whichvar: whichvar = "TSURF"
884    if whichvar not in fmtvar:
885        whichvar = "def"
886    return fmtvar[whichvar]
887
888## Author: AS
889####################################################################################################################
890### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
891def defcolorb (whichone="def"):
892    whichcolorb =    { \
893             "def":          "spectral",\
894             "HGT":          "spectral",\
895             "HGT_M":        "spectral",\
896             "TK":           "gist_heat",\
897             "TPOT":         "Paired",\
898             "TSURF":        "RdBu_r",\
899             "QH2O":         "PuBu",\
900             "USTM":         "YlOrRd",\
901             "WIND":         "YlOrRd",\
902             #"T_nadir_nit":  "RdBu_r",\
903             #"T_nadir_day":  "RdBu_r",\
904             "HFX":          "RdYlBu",\
905             "ICETOT":       "YlGnBu_r",\
906             #"MTOT":         "PuBu",\
907             "CCNQ":         "YlOrBr",\
908             "CCNN":         "YlOrBr",\
909             "TEMP":         "Jet",\
910             "TAU_ICE":      "Blues",\
911             "TAUICE":       "Blues",\
912             "VMR_ICE":      "Blues",\
913             "W":            "jet",\
914             "WMAX_TH":      "spectral",\
915             "ANOMALY":      "RdBu_r",\
916             "QSURFICE":     "hot_r",\
917             "ALBBARE":      "spectral",\
918             "TAU":          "YlOrBr_r",\
919             "CO2":          "YlOrBr_r",\
920             #### T.N.
921             "MTOT":         "Jet",\
922             "H2O_ICE_S":    "RdBu",\
923             "VMR_H2OICE":   "PuBu",\
924             "VMR_H2OVAP":   "PuBu",\
925             "WATERCAPTAG":  "Blues",\
926                     }
927#W --> spectral ou jet
928#spectral BrBG RdBu_r
929    #print "predefined colorbars"
930    if "TSURF" in whichone: whichone = "TSURF"
931    if whichone not in whichcolorb:
932        whichone = "def"
933    return whichcolorb[whichone]
934
935## Author: AS
936def definecolorvec (whichone="def"):
937        whichcolor =    { \
938                "def":          "black",\
939                "vis":          "yellow",\
940                "vishires":     "yellow",\
941                "molabw":       "yellow",\
942                "mola":         "black",\
943                "gist_heat":    "white",\
944                "hot":          "tk",\
945                "gist_rainbow": "black",\
946                "spectral":     "black",\
947                "gray":         "red",\
948                "PuBu":         "black",\
949                        }
950        if whichone not in whichcolor:
951                whichone = "def"
952        return whichcolor[whichone]
953
954## Author: AS
955def marsmap (whichone="vishires"):
956        from os import uname
957        mymachine = uname()[1]
958        ### not sure about speed-up with this method... looks the same
959        if "lmd.jussieu.fr" in mymachine:   domain = "/u/aslmd/WWW/maps/"
960        elif "aymeric-laptop" in mymachine: domain = "/home/aymeric/Dropbox/Public/"
961        else:                               domain = "http://www.lmd.jussieu.fr/~aslmd/maps/"
962        whichlink =     { \
963                #"vis":         "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\
964                #"vishires":    "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\
965                #"geolocal":    "http://dl.dropbox.com/u/11078310/geolocal.jpg",\
966                #"mola":        "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\
967                #"molabw":      "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\
968                "thermalday":   domain+"thermalday.jpg",\
969                "thermalnight": domain+"thermalnight.jpg",\
970                "tesalbedo":    domain+"tesalbedo.jpg",\
971                "vis":         domain+"mar0kuu2.jpg",\
972                "vishires":    domain+"MarsMap_2500x1250.jpg",\
973                "geolocal":    domain+"geolocal.jpg",\
974                "mola":        domain+"mars-mola-2k.jpg",\
975                "molabw":      domain+"MarsElevation_2500x1250.jpg",\
976                "clouds":      "http://www.johnstonsarchive.net/spaceart/marswcloudmap.jpg",\
977                "jupiter":     "http://www.mmedia.is/~bjj/data/jupiter_css/jupiter_css.jpg",\
978                "jupiter_voy": "http://www.mmedia.is/~bjj/data/jupiter/jupiter_vgr2.jpg",\
979                #"bw":          domain+"EarthElevation_2500x1250.jpg",\
980                "bw":          "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg",\
981                "contrast":    "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg",\
982                "nice":        "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg",\
983                "blue":        "http://eoimages.gsfc.nasa.gov/ve/2430/land_ocean_ice_2048.jpg",\
984                "blueclouds":  "http://eoimages.gsfc.nasa.gov/ve/2431/land_ocean_ice_cloud_2048.jpg",\
985                "justclouds":  "http://eoimages.gsfc.nasa.gov/ve/2432/cloud_combined_2048.jpg",\
986                        }
987        ### see http://www.mmedia.is/~bjj/planetary_maps.html
988        if whichone not in whichlink: 
989                print "marsmap: choice not defined... you'll get the default one... "
990                whichone = "vishires" 
991        return whichlink[whichone]
992
993#def earthmap (whichone):
994#       if   whichone == "contrast":    whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg"
995#       elif whichone == "bw":          whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg"
996#       elif whichone == "nice":        whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg"
997#       return whichlink
998
999## Author: AS
1000def latinterv (area="Whole"):
1001    list =    { \
1002        "Europe":                [[ 20., 80.],[- 50.,  50.]],\
1003        "Central_America":       [[-10., 40.],[ 230., 300.]],\
1004        "Africa":                [[-20., 50.],[- 50.,  50.]],\
1005        "Whole":                 [[-90., 90.],[-180., 180.]],\
1006        "Southern_Hemisphere":   [[-90., 60.],[-180., 180.]],\
1007        "Northern_Hemisphere":   [[-60., 90.],[-180., 180.]],\
1008        "Tharsis":               [[-30., 60.],[-170.,- 10.]],\
1009        "Whole_No_High":         [[-60., 60.],[-180., 180.]],\
1010        "Chryse":                [[-60., 60.],[- 60.,  60.]],\
1011        "North_Pole":            [[ 50., 90.],[-180., 180.]],\
1012        "Close_North_Pole":      [[ 75., 90.],[-180., 180.]],\
1013        "Far_South_Pole":        [[-90.,-40.],[-180., 180.]],\
1014        "South_Pole":            [[-90.,-50.],[-180., 180.]],\
1015        "Close_South_Pole":      [[-90.,-75.],[-180., 180.]],\
1016        "Sirenum_Crater_large":  [[-46.,-34.],[-166., -151.]],\
1017        "Sirenum_Crater_small":  [[-36.,-26.],[-168., -156.]],\
1018
1019              }
1020    if area not in list:   area = "Whole"
1021    [olat,olon] = list[area]
1022    return olon,olat
1023
1024## Author: TN
1025def separatenames (name):
1026  from numpy import concatenate
1027  # look for comas in the input name to separate different names (files, variables,etc ..)
1028  if name is None:
1029     names = None
1030  else:
1031    names = []
1032    stop = 0
1033    currentname = name
1034    while stop == 0:
1035      indexvir = currentname.find(',')
1036      if indexvir == -1:
1037        stop = 1
1038        name1 = currentname
1039      else:
1040        name1 = currentname[0:indexvir]
1041      names = concatenate((names,[name1]))
1042      currentname = currentname[indexvir+1:len(currentname)]
1043  return names
1044
1045## Author: TN [old]
1046def getopmatrix (kind,n):
1047  import numpy as np
1048  # compute matrix of operations between files
1049  # the matrix is 'number of files'-square
1050  # 1: difference (row minus column), 2: add
1051  # not 0 in diag : just plot
1052  if n == 1:
1053    opm = np.eye(1)
1054  elif kind == 'basic':
1055    opm = np.eye(n)
1056  elif kind == 'difference1': # show differences with 1st file
1057    opm = np.zeros((n,n))
1058    opm[0,:] = 1
1059    opm[0,0] = 0
1060  elif kind == 'difference2': # show differences with 1st file AND show 1st file
1061    opm = np.zeros((n,n))
1062    opm[0,:] = 1
1063  else:
1064    opm = np.eye(n)
1065  return opm
1066 
1067## Author: TN [old]
1068def checkcoherence (nfiles,nlat,nlon,ntime):
1069  if (nfiles > 1):
1070     if (nlat > 1):
1071        errormess("what you asked is not possible !")
1072  return 1
1073
1074## Author: TN
1075def readslices(saxis):
1076  from numpy import empty
1077  if saxis == None:
1078     zesaxis = None
1079  else:
1080     zesaxis = empty((len(saxis),2))
1081     for i in range(len(saxis)):
1082        a = separatenames(saxis[i])
1083        if len(a) == 1:
1084           zesaxis[i,:] = float(a[0])
1085        else:
1086           zesaxis[i,0] = float(a[0])
1087           zesaxis[i,1] = float(a[1])
1088           
1089  return zesaxis
1090
1091## Author: TN
1092def readdata(data,datatype,coord1,coord2):
1093## Read sparse data
1094  from numpy import empty
1095  if datatype == 'txt':
1096     if len(data[coord1].shape) == 1:
1097         return data[coord1][:]
1098     elif len(data[coord1].shape) == 2:
1099         return data[coord1][:,int(coord2)-1]
1100     else:
1101         errormess('error in readdata')
1102  elif datatype == 'sav':
1103     return data[coord1][coord2]
1104  else:
1105     errormess(datatype+' type is not supported!')
1106
1107
1108## Author: AS
1109def bidimfind(lon2d,lat2d,vlon,vlat,file=None):
1110   import numpy as np
1111   import matplotlib.pyplot as mpl
1112   if vlat is None:    array = (lon2d - vlon)**2
1113   elif vlon is None:  array = (lat2d - vlat)**2
1114   else:               array = (lon2d - vlon)**2 + (lat2d - vlat)**2
1115   idy,idx = np.unravel_index( np.argmin(array), lon2d.shape )
1116   if vlon is not None:
1117      if (np.abs(lon2d[idy,idx]-vlon)) > 5: errormess("longitude not found ",printvar=lon2d)
1118   if vlat is not None:
1119      if (np.abs(lat2d[idy,idx]-vlat)) > 5: errormess("latitude not found ",printvar=lat2d)
1120   if file is not None:
1121      print idx,idy,lon2d[idy,idx],vlon
1122      print idx,idy,lat2d[idy,idx],vlat
1123      var = file.variables["HGT"][:,:,:]
1124      mpl.contourf(var[0,:,:],30,cmap = mpl.get_cmap(name="Greys_r") ) ; mpl.axis('off') ; mpl.plot(idx,idy,'mx',mew=4.0,ms=20.0)
1125      mpl.show()
1126   return idy,idx
1127
1128## Author: TN
1129def getsindex(saxis,index,axis):
1130# input  : all the desired slices and the good index
1131# output : all indexes to be taken into account for reducing field
1132  import numpy as np
1133  if ( np.array(axis).ndim == 2):
1134      axis = axis[:,0]
1135  if saxis is None:
1136      zeindex = None
1137  else:
1138      aaa = int(np.argmin(abs(saxis[index,0] - axis)))
1139      bbb = int(np.argmin(abs(saxis[index,1] - axis)))
1140      [imin,imax] = np.sort(np.array([aaa,bbb]))
1141      zeindex = np.array(range(imax-imin+1))+imin
1142      # because -180 and 180 are the same point in longitude,
1143      # we get rid of one for averaging purposes.
1144      if axis[imin] == -180 and axis[imax] == 180:
1145         zeindex = zeindex[0:len(zeindex)-1]
1146         print "INFO: whole longitude averaging asked, so last point is not taken into account."
1147  return zeindex
1148     
1149## Author: TN
1150def define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,indextime,what_I_plot,dim0,vertmode):
1151# Purpose of define_axis is to find x and y axis scales in a smart way
1152# x axis priority: 1/time 2/lon 3/lat 4/vertical
1153# To be improved !!!...
1154   from numpy import array,swapaxes
1155   x = None
1156   y = None
1157   count = 0
1158   what_I_plot = array(what_I_plot)
1159   shape = what_I_plot.shape
1160   if indextime is None and len(time) > 1:
1161      print "AXIS is time"
1162      x = time
1163      count = count+1
1164   if indexlon is None and len(lon) > 1:
1165      print "AXIS is lon"
1166      if count == 0: x = lon
1167      else: y = lon
1168      count = count+1
1169   if indexlat is None and len(lat) > 1:
1170      print "AXIS is lat"
1171      if count == 0: x = lat
1172      else: y = lat
1173      count = count+1
1174   if indexvert is None and ((dim0 == 4) or (y is None)):
1175      print "AXIS is vert"
1176      if vertmode == 0: # vertical axis is as is (GCM grid)
1177         if count == 0: x=range(len(vert))
1178         else: y=range(len(vert))
1179         count = count+1
1180      else: # vertical axis is in kms
1181         if count == 0: x = vert
1182         else: y = vert
1183         count = count+1
1184   x = array(x)
1185   y = array(y)
1186   print "CHECK SHAPE: what_I_plot, x, y", what_I_plot.shape, x.shape, y.shape
1187   if len(shape) == 1:
1188       if shape[0] != len(x):           print "WARNING: shape[0] != len(x). Correcting." ; what_I_plot = what_I_plot[0:len(x)]
1189       if len(y.shape) > 0:             y = ()
1190   elif len(shape) == 2:
1191       if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]:
1192           print "INFO: swapaxes: ",what_I_plot.shape,shape ; what_I_plot = swapaxes(what_I_plot,0,1)
1193       else:
1194           if shape[0] != len(y):       print "WARNING: shape[0] != len(y). Correcting." ; what_I_plot = what_I_plot[0:len(y),:]
1195           elif shape[1] != len(x):     print "WARNING: shape[1] != len(x). Correcting." ; what_I_plot = what_I_plot[:,0:len(x)]
1196   elif len(shape) == 3:
1197       if vertmode < 0: print "not supported. must check array dimensions at some point. not difficult to implement though."
1198   return what_I_plot,x,y
1199
1200# Author: TN + AS
1201def determineplot(slon, slat, svert, stime):
1202    nlon = 1 # number of longitudinal slices -- 1 is None
1203    nlat = 1
1204    nvert = 1
1205    ntime = 1
1206    nslices = 1
1207    if slon is not None:
1208        nslices = nslices*len(slon)
1209        nlon = len(slon)
1210    if slat is not None:
1211        nslices = nslices*len(slat)
1212        nlat = len(slat)
1213    if svert is not None:
1214        nslices = nslices*len(svert)
1215        nvert = len(svert)
1216    if stime is not None:
1217        nslices = nslices*len(stime)
1218        ntime = len(stime)
1219    #else:
1220    #    nslices = 2 
1221    mapmode = 0
1222    if slon is None and slat is None:
1223       mapmode = 1 # in this case we plot a map, with the given projection
1224
1225    return nlon, nlat, nvert, ntime, mapmode, nslices
1226
1227## Author: AC
1228## Reduce complexity of main script by moving the contour part here. Also allow to call it from elsewhere
1229## which can be usefull
1230#
1231#def call_contour(what_I_plot,error,x,y,m,lon,lat,vert,time,vertmode,ze_var2,indextime,indexlon,indexlat,indexvert,yintegral,mapmode,typefile,var2,ticks):
1232#      from matplotlib.pyplot import contour, plot, clabel
1233#      import numpy as np
1234#      #what_I_plot = what_I_plot*mult
1235#      if not error:
1236#         if mapmode == 1:
1237#             if typefile in ['mesoapi','meso']:    what_I_plot = dumpbdy(what_I_plot,6)
1238#         zevmin, zevmax = calculate_bounds(what_I_plot)
1239#         zelevels = np.linspace(zevmin,zevmax,ticks) #20)
1240#         if var2 == 'HGT':  zelevels = np.arange(-10000.,30000.,2000.)
1241#         if mapmode == 0:
1242#             #if typefile in ['mesoideal']:    what_I_plot = dumpbdy(what_I_plot,0,stag='W')
1243#             itime=indextime
1244#             if len(what_I_plot.shape) is 3: itime=[0]
1245#             what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
1246#                   itime,what_I_plot, len(ze_var2.shape),vertmode)
1247#         ### If we plot a 2-D field
1248#         if len(what_I_plot.shape) is 2:
1249#             #zelevels=[1.]
1250#             if mapmode == 0:cs = contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5)
1251#             elif mapmode == 1:cs = m.contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5)
1252#             #clabel(cs,zelevels,inline=3,fmt='%1.1f',fontsize=7)
1253#         ### If we plot a 1-D field
1254#         elif len(what_I_plot.shape) is 1:
1255#             plot(what_I_plot,x)
1256#      else:
1257#         errormess("There is an error in reducing field !")
1258#      return error
1259
Note: See TracBrowser for help on using the repository browser.