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

Last change on this file since 782 was 782, checked in by acolaitis, 12 years ago

UTIL PYTHON fixed mesoscale local time bugs introduced at revision 763. now working fine for LMD and MRAMS-generated wrfout files. plus polar domain now OK

File size: 67.9 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 + AC
32def localtime(time,lon,namefile): # lon is the mean longitude of the plot, not of the domain. central lon of domain is taken from cen_lon
33    import numpy as np
34    from netCDF4 import Dataset
35    ## THIS IS FOR MESOSCALE
36    nc = Dataset(namefile)
37    ## get start date and intervals
38    dt_hour=1. ; start=0.
39    if hasattr(nc,'TITLE'):
40        title=getattr(nc, 'TITLE')
41        if hasattr(nc,'DT') and hasattr(nc,'START_DATE') and 'MRAMS' in title: 
42            ## we must adapt what is done in getlschar to MRAMS (outputs from ic.py)
43            dt_hour=getattr(nc, 'DT')/60. 
44            start_date=getattr(nc, 'START_DATE')
45            start_hour=np.float(start_date[11:13])
46            start_minute=np.float(start_date[14:16])/60.
47            start=start_hour+start_minute # start is the local time of simu at longitude 0
48            #LMD MMM is 1 output/hour (and not 1 output/timestep)
49            #MRAMS is 1 output/timestep, unless stride is added in ic.py
50        elif 'WRF' in title: 
51            [dummy,start,dt_hour] = getlschar ( namefile ) # get start hour and interval hour
52    ## get longitude
53    if lon is not None:
54       if lon[0,1]!=lon[0,0]: mean_lon_plot = 0.5*(lon[0,1]-lon[0,0])
55       else: mean_lon_plot=lon[0,0]
56    elif hasattr(nc, 'CEN_LON'): mean_lon_plot=getattr(nc, 'CEN_LON')
57    else: mean_lon_plot=0.
58    ## calculate local time
59    ltst = start + time*dt_hour + mean_lon_plot / 15.
60    ltst = int (ltst * 10) / 10.
61    ltst = ltst % 24
62    return ltst
63
64## Author: AC
65def check_localtime(time):
66    a=-1
67    print time
68    for i in range(len(time)-1):
69       if (time[i] > time[i+1]): a=i
70    if a >= 0 and a < (len(time)-1)/2.:
71       print "Sorry, time axis is not regular."
72       print "Contourf needs regular axis... recasting"
73       for i in range(a+1):
74          time[i]=time[i]-24.
75    if a >= 0 and a >= (len(time)-1)/2.:
76       print "Sorry, time axis is not regular."
77       print "Contourf needs regular axis... recasting"
78       for i in range((len(time)-1) - a):
79          time[a+1+i]=time[a+1+i]+24.
80    return time
81
82## Author: AS, AC, JL
83def whatkindfile (nc):
84    typefile = 'gcm' # default
85    if 'controle' in nc.variables:             typefile = 'gcm'
86    elif 'phisinit' in nc.variables:           typefile = 'gcm'
87    elif 'phis' in nc.variables:               typefile = 'gcm'
88    elif 'time_counter' in nc.variables:       typefile = 'earthgcm'
89    elif hasattr(nc,'START_DATE'):             typefile = 'meso' 
90    elif 'HGT_M' in nc.variables:              typefile = 'geo'
91    elif hasattr(nc,'institution'):
92      if "European Centre" in getattr(nc,'institution'):  typefile = 'ecmwf'
93    return typefile
94
95## Author: AS
96def getfield (nc,var):
97    ## this allows to get much faster (than simply referring to nc.variables[var])
98    import numpy as np
99    dimension = len(nc.variables[var].dimensions)
100    #print "   Opening variable with", dimension, "dimensions ..."
101    if dimension == 2:    field = nc.variables[var][:,:]
102    elif dimension == 3:  field = nc.variables[var][:,:,:]
103    elif dimension == 4:  field = nc.variables[var][:,:,:,:]
104    elif dimension == 1:  field = nc.variables[var][:]
105    # if there are NaNs in the ncdf, they should be loaded as a masked array which will be
106    # recasted as a regular array later in reducefield
107    if (np.isnan(np.sum(field)) and (type(field).__name__ not in 'MaskedArray')):
108       print "Warning: netcdf as nan values but is not loaded as a Masked Array."
109       print "recasting array type"
110       out=np.ma.masked_invalid(field)
111       out.set_fill_value([np.NaN])
112    else:
113    # missing values from zrecast or hrecast are -1e-33
114       masked=np.ma.masked_where(field < -1e30,field)
115       masked2=np.ma.masked_where(field > 1e35,field)
116       masked.set_fill_value([np.NaN]) ; masked2.set_fill_value([np.NaN])
117       mask = np.ma.getmask(masked) ; mask2 = np.ma.getmask(masked2)
118       if (True in np.array(mask)):
119          out=masked
120          print "Masked array... Missing value is NaN"
121       elif (True in np.array(mask2)):
122          out=masked2
123          print "Masked array... Missing value is NaN"
124#       else:
125#       # missing values from api are 1e36
126#          masked=np.ma.masked_where(field > 1e35,field)
127#          masked.set_fill_value([np.NaN])
128#          mask = np.ma.getmask(masked)
129#          if (True in np.array(mask)):out=masked
130#          else:out=field
131       else:
132#       # missing values from MRAMS files are 0.100E+32
133          masked=np.ma.masked_where(field > 1e30,field)
134          masked.set_fill_value([np.NaN])
135          mask = np.ma.getmask(masked)
136          if (True in np.array(mask)):out=masked
137          else:out=field
138#       else:out=field
139    return out
140
141## Author: AC
142# Compute the norm of the winds or return an hodograph
143# The corresponding variable to call is UV or uvmet (to use api)
144def windamplitude (nc,mode):
145    import numpy as np
146    varinfile = nc.variables.keys()
147    if "U" in varinfile: zu=getfield(nc,'U')
148    elif "Um" in varinfile: zu=getfield(nc,'Um')
149    else: errormess("you need slopex or U or Um in your file.")
150    if "V" in varinfile: zv=getfield(nc,'V')
151    elif "Vm" in varinfile: zv=getfield(nc,'Vm')
152    else: errormess("you need V or Vm in your file.")
153    znt,znz,zny,znx = np.array(zu).shape
154    if hasattr(nc,'WEST-EAST_PATCH_END_UNSTAG'):znx=getattr(nc, 'WEST-EAST_PATCH_END_UNSTAG')
155    zuint = np.zeros([znt,znz,zny,znx])
156    zvint = np.zeros([znt,znz,zny,znx])
157    if "U" in varinfile:
158       if hasattr(nc,'SOUTH-NORTH_PATCH_END_STAG'): zny_stag=getattr(nc, 'SOUTH-NORTH_PATCH_END_STAG')
159       if hasattr(nc,'WEST-EAST_PATCH_END_STAG'): znx_stag=getattr(nc, 'WEST-EAST_PATCH_END_STAG')
160       if zny_stag == zny: zvint=zv
161       else:
162          for yy in np.arange(zny):      zvint[:,:,yy,:] = (zv[:,:,yy,:] + zv[:,:,yy+1,:])/2.
163       if znx_stag == znx: zuint=zu
164       else:
165          for xx in np.arange(znx):      zuint[:,:,:,xx] = (zu[:,:,:,xx] + zu[:,:,:,xx+1])/2.
166    else:
167       zuint=zu
168       zvint=zv
169    if mode=='amplitude': return np.sqrt(zuint**2 + zvint**2)
170    if mode=='hodograph': return zuint,zvint
171    if mode=='hodograph_2': return None, 360.*np.arctan(zvint/zuint)/(2.*np.pi)
172
173## Author: AC
174# Compute the enrichment factor of non condensible gases
175# The corresponding variable to call is enfact
176# enrichment factor is computed as in Yuan Lian et al. 2012
177# i.e. you need to have VL2 site at LS 135 in your data
178# this only requires co2col so that you can concat.nc at low cost
179def enrichment_factor(nc,lon,lat,time):
180    import numpy as np
181    from myplot import reducefield
182    varinfile = nc.variables.keys()
183    if "co2col" in varinfile: co2col=getfield(nc,'co2col')
184    else: print "error, you need co2col var in your file"
185    if "ps" in varinfile: ps=getfield(nc,'ps')
186    else: print "error, you need ps var in your file"
187    dimension = len(nc.variables['co2col'].dimensions)
188    if dimension == 2: 
189      zny,znx = np.array(co2col).shape
190      znt=1
191    elif dimension == 3: znt,zny,znx = np.array(co2col).shape
192    mmrarcol = np.zeros([znt,zny,znx])
193    enfact = np.zeros([znt,zny,znx])
194    grav=3.72
195    mmrarcol[:,:,:] = 1. - grav*co2col[:,:,:]/ps[:,:,:]
196# Computation with reference argon mmr at VL2 Ls 135 (as in Yuan Lian et al 2012)
197    lonvl2=np.zeros([1,2])
198    latvl2=np.zeros([1,2])
199    timevl2=np.zeros([1,2])
200    lonvl2[0,0]=-180
201    lonvl2[0,1]=180
202    latvl2[:,:]=48.16
203    timevl2[:,:]=135.
204    indexlon  = getsindex(lonvl2,0,lon)
205    indexlat  = getsindex(latvl2,0,lat)
206    indextime = getsindex(timevl2,0,time)
207    mmrvl2, error = reducefield( mmrarcol, d4=indextime, d1=indexlon, d2=indexlat)
208    print "VL2 Ls 135 mmr arcol:", mmrvl2
209    enfact[:,:,:] = mmrarcol[:,:,:]/mmrvl2
210    return enfact
211
212## Author: AC
213# Compute the norm of the slope angles
214# The corresponding variable to call is SLOPEXY
215def slopeamplitude (nc):
216    import numpy as np
217    varinfile = nc.variables.keys()
218    if "slopex" in varinfile: zu=getfield(nc,'slopex')
219    elif "SLOPEX" in varinfile: zu=getfield(nc,'SLOPEX')
220    else: errormess("you need slopex or SLOPEX in your file.") 
221    if "slopey" in varinfile: zv=getfield(nc,'slopey')
222    elif "SLOPEY" in varinfile: zv=getfield(nc,'SLOPEY')
223    else: errormess("you need slopey or SLOPEY in your file.")
224    znt,zny,znx = np.array(zu).shape
225    zuint = np.zeros([znt,zny,znx])
226    zvint = np.zeros([znt,zny,znx])
227    zuint=zu
228    zvint=zv
229    return np.sqrt(zuint**2 + zvint**2)
230
231## Author: AC
232# Compute the temperature difference between surface and first level.
233# API is automatically called to get TSURF and TK.
234# The corresponding variable to call is DELTAT
235def deltat0t1 (nc):
236    import numpy as np
237    varinfile = nc.variables.keys()
238    if "tsurf" in varinfile: zu=getfield(nc,'tsurf')
239    elif "TSURF" in varinfile: zu=getfield(nc,'TSURF')
240    else: errormess("You need tsurf or TSURF in your file")
241    if "tk" in varinfile: zv=getfield(nc,'tk')
242    elif "TK" in varinfile: zv=getfield(nc,'TK')
243    else: errormess("You need tk or TK in your file. (might need to use API. try to add -i 4 -l XXX)")
244    znt,zny,znx = np.array(zu).shape
245    zuint = np.zeros([znt,zny,znx])
246    zuint=zu - zv[:,0,:,:]
247    return zuint
248
249## Author: AS + TN + AC
250def reducefield (input,d4=None,d3=None,d2=None,d1=None,yint=False,alt=None,anomaly=False,redope=None,mesharea=None,unidim=999):
251    ### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x"
252    ### it would be actually better to name d4 d3 d2 d1 as t z y x
253    ### ... note, anomaly is only computed over d1 and d2 for the moment
254    import numpy as np
255    from mymath import max,mean,min,sum,getmask
256    csmooth = 12 ## a fair amount of grid points (too high results in high computation time)
257    if redope is not None:
258       if   redope == "mint":     input = min(input,axis=0) ; d1 = None
259       elif redope == "maxt":     input = max(input,axis=0) ; d1 = None
260       elif redope == "edge_y1":  input = input[:,:,0,:]    ; d2 = None
261       elif redope == "edge_y2":  input = input[:,:,-1,:]   ; d2 = None
262       elif redope == "edge_x1":  input = input[:,:,:,0]    ; d1 = None
263       elif redope == "edge_x2":  input = input[:,:,:,-1]   ; d1 = None
264       else:                      errormess("not supported. but try lines in reducefield beforehand.")
265       #elif redope == "minz":     input = min(input,axis=1) ; d2 = None
266       #elif redope == "maxz":     input = max(input,axis=1) ; d2 = None
267       #elif redope == "miny":     input = min(input,axis=2) ; d3 = None
268       #elif redope == "maxy":     input = max(input,axis=2) ; d3 = None
269       #elif redope == "minx":     input = min(input,axis=3) ; d4 = None
270       #elif redope == "maxx":     input = max(input,axis=3) ; d4 = None
271    dimension = np.array(input).ndim
272    shape = np.array(np.array(input).shape)
273    #print 'd1,d2,d3,d4: ',d1,d2,d3,d4
274    if anomaly: print 'ANOMALY ANOMALY'
275    output = input
276    error = False
277    #### this is needed to cope the case where d4,d3,d2,d1 are single integers and not arrays
278    if d4 is not None and not isinstance(d4, np.ndarray): d4=[d4]
279    if d3 is not None and not isinstance(d3, np.ndarray): d3=[d3]
280    if d2 is not None and not isinstance(d2, np.ndarray): d2=[d2]
281    if d1 is not None and not isinstance(d1, np.ndarray): d1=[d1]
282    ### now the main part
283    if dimension == 2:
284        #### this is needed for 1d-type files (where dim=2 but axes are time-vert and not lat-lon)
285        if unidim==1: d2=d4 ; d1=d3 ; d4 = None ; d3 = None
286        if mesharea is None: mesharea=np.ones(shape)
287        if   max(d2) >= shape[0]: error = True
288        elif max(d1) >= shape[1]: error = True
289        elif d1 is not None and d2 is not None:
290          try:
291            totalarea = np.ma.masked_where(getmask(output),mesharea)
292            totalarea = mean(totalarea[d2,:],axis=0);totalarea = mean(totalarea[d1])
293          except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
294          output = output*mesharea; output = mean(output[d2,:],axis=0); output = mean(output[d1])/totalarea
295        elif d1 is not None:         output = mean(input[:,d1],axis=1)
296        elif d2 is not None:
297          try:
298            totalarea = np.ma.masked_where(getmask(output),mesharea)
299            totalarea = mean(totalarea[d2,:],axis=0)
300          except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
301          output = output*mesharea; output = mean(output[d2,:],axis=0)/totalarea
302    elif dimension == 3:
303        if mesharea is None: mesharea=np.ones(shape[[1,2]])
304        if   max(d4) >= shape[0]: error = True
305        elif max(d2) >= shape[1]: error = True
306        elif max(d1) >= shape[2]: error = True
307        elif d4 is not None and d2 is not None and d1 is not None:
308          output = mean(input[d4,:,:],axis=0)
309          try:
310            totalarea = np.ma.masked_where(getmask(output),mesharea)
311            totalarea = mean(totalarea[d2,:],axis=0);totalarea = mean(totalarea[d1])
312          except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
313          output = output*mesharea; output = mean(output[d2,:],axis=0); output = mean(output[d1])/totalarea
314        elif d4 is not None and d2 is not None:
315          output = mean(input[d4,:,:],axis=0)
316          try:
317            totalarea = np.ma.masked_where(getmask(output),mesharea)
318            totalarea = mean(totalarea[d2,:],axis=0)
319          except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
320          output = output*mesharea; output = mean(output[d2,:],axis=0)/totalarea
321        elif d4 is not None and d1 is not None:    output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1)
322        elif d2 is not None and d1 is not None:
323          try:
324            totalarea = np.tile(mesharea,(output.shape[0],1,1))
325            totalarea = np.ma.masked_where(getmask(output),totalarea)
326            totalarea = mean(totalarea[:,d2,:],axis=1);totalarea = mean(totalarea[:,d1],axis=1)
327          except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
328          output = output*mesharea; output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)/totalarea
329        elif d1 is not None:   output = mean(input[:,:,d1],axis=2)
330        elif d2 is not None:   
331          try:
332            totalarea = np.tile(mesharea,(output.shape[0],1,1))
333            totalarea = np.ma.masked_where(getmask(output),totalarea)
334            totalarea = mean(totalarea[:,d2,:],axis=1)
335          except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
336          output = output*mesharea; output = mean(output[:,d2,:],axis=1)/totalarea
337        elif d4 is not None:   output = mean(input[d4,:,:],axis=0)
338    elif dimension == 4:
339        if mesharea is None: mesharea=np.ones(shape[[2,3]]) # mesharea=np.random.random_sample(shape[[2,3]])*5. + 2. # pour tester
340        if   max(d4) >= shape[0]: error = True
341        elif max(d3) >= shape[1]: error = True
342        elif max(d2) >= shape[2]: error = True
343        elif max(d1) >= shape[3]: error = True
344        elif d4 is not None and d3 is not None and d2 is not None and d1 is not None:
345             output = mean(input[d4,:,:,:],axis=0)
346             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
347             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
348             try:
349               totalarea = np.ma.masked_where(np.isnan(output),mesharea)
350               totalarea = mean(totalarea[d2,:],axis=0); totalarea = mean(totalarea[d1])
351             except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
352             output = output*mesharea; output = mean(output[d2,:],axis=0); output = mean(output[d1])/totalarea
353        elif d4 is not None and d3 is not None and d2 is not None: 
354             output = mean(input[d4,:,:,:],axis=0)
355             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
356             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
357             try:
358               totalarea = np.ma.masked_where(np.isnan(output),mesharea)
359               totalarea = mean(totalarea[d2,:],axis=0)
360             except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
361             output = output*mesharea; output = mean(output[d2,:],axis=0)/totalarea
362        elif d4 is not None and d3 is not None and d1 is not None: 
363             output = mean(input[d4,:,:,:],axis=0)
364             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
365             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) 
366             output = mean(output[:,d1],axis=1)
367        elif d4 is not None and d2 is not None and d1 is not None: 
368             output = mean(input[d4,:,:,:],axis=0)
369             if anomaly:
370                 for k in range(output.shape[0]):  output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
371             try:
372               totalarea = np.tile(mesharea,(output.shape[0],1,1))
373               totalarea = np.ma.masked_where(getmask(output),mesharea)
374               totalarea = mean(totalarea[:,d2,:],axis=1); totalarea = mean(totalarea[:,d1],axis=1)
375             except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
376             output = output*mesharea; output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)/totalarea
377             #noperturb = smooth1d(output,window_len=7)
378             #lenlen = len(output) ; output = output[1:lenlen-7] ; yeye = noperturb[4:lenlen-4]
379             #plot(output) ; plot(yeye) ; show() ; plot(output-yeye) ; show()
380        elif d3 is not None and d2 is not None and d1 is not None:
381             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
382             if anomaly:
383                 for k in range(output.shape[0]):  output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
384             try:
385               totalarea = np.tile(mesharea,(output.shape[0],1,1))
386               totalarea = np.ma.masked_where(getmask(output),totalarea)
387               totalarea = mean(totalarea[:,d2,:],axis=1); totalarea = mean(totalarea[:,d1],axis=1)
388             except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
389             output = output*mesharea; output = mean(output[:,d2,:],axis=1); output = mean(output[:,d1],axis=1)/totalarea
390        elif d4 is not None and d3 is not None: 
391             output = mean(input[d4,:,:,:],axis=0) 
392             output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
393             if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.) 
394        elif d4 is not None and d2 is not None: 
395             output = mean(input[d4,:,:,:],axis=0)
396             try:
397               totalarea = np.tile(mesharea,(output.shape[0],1,1))
398               totalarea = np.ma.masked_where(getmask(output),mesharea)
399               totalarea = mean(totalarea[:,d2,:],axis=1)
400             except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
401             output = output*mesharea; output = mean(output[:,d2,:],axis=1)/totalarea
402        elif d4 is not None and d1 is not None: 
403             output = mean(input[d4,:,:,:],axis=0) 
404             output = mean(output[:,:,d1],axis=2)
405        elif d3 is not None and d2 is not None:
406             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
407             try:
408               totalarea = np.tile(mesharea,(output.shape[0],1,1))
409               totalarea = np.ma.masked_where(getmask(output),mesharea)
410               totalarea = mean(totalarea[:,d2,:],axis=1)
411             except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
412             output = output*mesharea; output = mean(output[:,d2,:],axis=1)/totalarea
413        elif d3 is not None and d1 is not None: 
414             output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
415             output = mean(output[:,:,d1],axis=2)
416        elif d2 is not None and d1 is not None:
417             try:
418               totalarea = np.tile(mesharea,(output.shape[0],output.shape[1],1,1))
419               totalarea = np.ma.masked_where(getmask(output),totalarea)
420               totalarea = mean(totalarea[:,:,d2,:],axis=2); totalarea = mean(totalarea[:,:,d1],axis=1)
421             except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
422             output = output*mesharea; output = mean(output[:,:,d2,:],axis=2); output = mean(output[:,:,d1],axis=2)/totalarea
423        elif d1 is not None:        output = mean(input[:,:,:,d1],axis=3)
424        elif d2 is not None:
425             try:
426               totalarea = np.tile(mesharea,(output.shape[0],output.shape[1],1,output.shape[3]))
427               totalarea = np.ma.masked_where(getmask(output),totalarea)
428               totalarea = mean(totalarea[:,:,d2,:],axis=2)
429             except: print "(problem with areas. I skip this)" ; mesharea = 1. ; totalarea = 1.
430             output = output*mesharea; output = mean(output[:,:,d2,:],axis=2)/totalarea
431        elif d3 is not None:        output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
432        elif d4 is not None:        output = mean(input[d4,:,:,:],axis=0)
433    dimension2 = np.array(output).ndim
434    shape2 = np.array(output).shape
435    print 'REDUCEFIELD dim,shape: ',dimension,shape,' >>> ',dimension2,shape2
436    return output, error
437
438## Author: AC + AS
439def reduce_zaxis (input,ax=None,yint=False,vert=None,indice=None):
440    from mymath import max,mean
441    from scipy import integrate
442    import numpy as np
443    if yint and vert is not None and indice is not None:
444       if type(input).__name__=='MaskedArray':
445          input.set_fill_value([np.NaN])
446          output = integrate.trapz(input.filled(),x=vert[indice],axis=ax)
447       else:
448          output = integrate.trapz(input,x=vert[indice],axis=ax)
449    else:
450       output = mean(input,axis=ax)
451    return output
452
453## Author: AS + TN
454def definesubplot ( numplot, fig ):
455    from matplotlib.pyplot import rcParams
456    rcParams['font.size'] = 12. ## default (important for multiple calls)
457    if numplot <= 0:
458        subv = 99999
459        subh = 99999
460    elif numplot == 1: 
461        subv = 1
462        subh = 1
463    elif numplot == 2:
464        subv = 1 #2
465        subh = 2 #1
466        fig.subplots_adjust(wspace = 0.35)
467        rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
468    elif numplot == 3:
469        subv = 3
470        subh = 1
471        fig.subplots_adjust(hspace = 0.5)
472        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
473    elif numplot == 4:
474        subv = 2
475        subh = 2
476        #fig.subplots_adjust(wspace = 0.4, hspace = 0.6)
477        fig.subplots_adjust(wspace = 0.4, hspace = 0.3)
478        rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
479    elif numplot <= 6:
480        subv = 2
481        subh = 3
482        #fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
483        fig.subplots_adjust(wspace = 0.5, hspace = 0.3)
484        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
485    elif numplot <= 8:
486        subv = 2
487        subh = 4
488        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
489        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
490    elif numplot <= 9:
491        subv = 3
492        subh = 3
493        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
494        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
495    elif numplot <= 12:
496        subv = 3
497        subh = 4
498        fig.subplots_adjust(wspace = 0, hspace = 0.1)
499        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
500    elif numplot <= 16:
501        subv = 4
502        subh = 4
503        fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
504        rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
505    else:
506        print "number of plot supported: 1 to 16"
507        exit()
508    return subv,subh
509
510## Author: AS
511def getstralt(nc,nvert):
512    varinfile = nc.variables.keys()
513    if 'vert' not in varinfile:
514        stralt = "_lvl" + str(nvert)
515    else:
516        zelevel = int(nc.variables['vert'][nvert])
517        if abs(zelevel) < 10000.:   strheight=str(zelevel)+"m"
518        else:                       strheight=str(int(zelevel/1000.))+"km"
519        if 'altitude'       in nc.dimensions:   stralt = "_"+strheight+"-AMR"
520        elif 'altitude_abg' in nc.dimensions:   stralt = "_"+strheight+"-ALS"
521        elif 'bottom_top'   in nc.dimensions:   stralt = "_"+strheight
522        elif 'pressure'     in nc.dimensions:   stralt = "_"+str(zelevel)+"Pa"
523        else:                                   stralt = ""
524    return stralt
525
526## Author: AS
527def getlschar ( namefile, getaxis=False ):
528    from netCDF4 import Dataset
529    from timestuff import sol2ls
530    from numpy import array
531    from string import rstrip
532    import os as daos
533    namefiletest = rstrip( rstrip( rstrip( namefile, chars="_z"), chars="_zabg"), chars="_p")
534    testexist = daos.path.isfile(namefiletest)
535    zetime = None
536    if testexist: 
537      namefile = namefiletest
538      #### we assume that wrfout is next to wrfout_z and wrfout_zabg
539      nc  = Dataset(namefile)
540      zetime = None
541      days_in_month = [61, 66, 66, 65, 60, 54, 50, 46, 47, 47, 51, 56]
542      plus_in_month = [ 0, 61,127,193,258,318,372,422,468,515,562,613]
543      if 'Times' in nc.variables:
544        zetime = nc.variables['Times'][0]
545        shape = array(nc.variables['Times']).shape
546        if shape[0] < 2: zetime = None
547    if zetime is not None \
548       and 'vert' not in nc.variables:
549        ##### strangely enough this does not work for api or ncrcat results!
550        zesol = plus_in_month[ int(zetime[5]+zetime[6])-1 ] + int(zetime[8]+zetime[9]) - 1 ##les sols GCM commencent a 0
551        dals = int( 10. * sol2ls ( zesol ) ) / 10.
552        ###
553        zetime2 = nc.variables['Times'][1]
554        one  = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37.
555        next = int(zetime2[11]+zetime2[12]) + int(zetime2[14]+zetime2[15])/37. 
556        zehour    = one
557        zehourin  = abs ( next - one )
558        if not getaxis:
559            lschar = "_Ls"+str(dals)
560        else:
561            zelen = len(nc.variables['Times'][:])
562            yeye = range(zelen) ; lsaxis = range(zelen) ; solaxis = range(zelen) ; ltaxis = range(zelen)
563            for iii in yeye:
564               zetime = nc.variables['Times'][iii] 
565               ltaxis[iii] = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37.
566               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
567               lsaxis[iii] = sol2ls ( solaxis[iii] )
568               if ltaxis[iii] < ltaxis[iii-1]: ltaxis[iii] = ltaxis[iii] + 24.
569               #print ltaxis[iii], solaxis[iii], lsaxis[iii], getattr( nc, 'JULDAY' )
570            lschar = lsaxis ; zehour = solaxis ; zehourin = ltaxis
571    else:
572        lschar=""
573        zehour = 0
574        zehourin = 1 
575    return lschar, zehour, zehourin
576
577## Author: AS
578def getprefix (nc):
579    prefix = 'LMD_MMM_'
580    prefix = prefix + 'd'+str(getattr(nc,'GRID_ID'))+'_'
581    prefix = prefix + str(int(getattr(nc,'DX')/1000.))+'km_'
582    return prefix
583
584## Author: AS
585def getproj (nc):
586    typefile = whatkindfile(nc)
587    if typefile in ['meso','geo']:
588        ### (il faudrait passer CEN_LON dans la projection ?)
589        map_proj = getattr(nc, 'MAP_PROJ')
590        cen_lat  = getattr(nc, 'CEN_LAT')
591        if map_proj == 2:
592            if cen_lat > 10.:   
593                proj="npstere"
594                #print "NP stereographic polar domain"
595            else:           
596                proj="spstere"
597                #print "SP stereographic polar domain"
598        elif map_proj == 1: 
599            #print "lambert projection domain"
600            proj="lcc"
601        elif map_proj == 3: 
602            #print "mercator projection"
603            proj="merc"
604        else:
605            proj="merc"
606    elif typefile in ['gcm']:        proj="cyl"    ## pb avec les autres (de trace derriere la sphere ?)
607    else:                            proj="ortho"
608    return proj   
609
610## Author: AS
611def ptitle (name):
612    from matplotlib.pyplot import title
613    title(name)
614    print name
615
616## Author: AS
617def polarinterv (lon2d,lat2d):
618    import numpy as np
619    wlon = [np.min(lon2d),np.max(lon2d)]
620    ind = np.array(lat2d).shape[0] / 2  ## to get a good boundlat and to get the pole
621    wlat = [np.min(lat2d[ind,:]),np.max(lat2d[ind,:])]
622    return [wlon,wlat]
623
624## Author: AS
625def simplinterv (lon2d,lat2d):
626    import numpy as np
627    return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]]
628
629## Author: AS
630def wrfinterv (lon2d,lat2d):
631    nx = len(lon2d[0,:])-1
632    ny = len(lon2d[:,0])-1
633    lon1 = lon2d[0,0] 
634    lon2 = lon2d[nx,ny] 
635    lat1 = lat2d[0,0] 
636    lat2 = lat2d[nx,ny] 
637    if abs(0.5*(lat1+lat2)) > 60.:  wider = 0.5 * (abs(lon1)+abs(lon2)) * 0.1
638    else:                           wider = 0.
639    if lon1 < lon2:  wlon = [lon1, lon2 + wider] 
640    else:            wlon = [lon2, lon1 + wider]
641    if lat1 < lat2:  wlat = [lat1, lat2]
642    else:            wlat = [lat2, lat1]
643    return [wlon,wlat]
644
645## Author: AS
646def makeplotres (filename,res=None,pad_inches_value=0.25,folder='',disp=True,ext='png',erase=False):
647    import  matplotlib.pyplot as plt
648    from os import system
649    addstr = ""
650    if res is not None:
651        res = int(res)
652        addstr = "_"+str(res)
653    name = filename+addstr+"."+ext
654    if folder != '':      name = folder+'/'+name
655    plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value)
656    if disp:              display(name)
657    if ext in ['eps','ps','svg']:   system("tar czvf "+name+".tar.gz "+name+" ; rm -f "+name)
658    if erase:   system("mv "+name+" to_be_erased")             
659    return
660
661## Author: AS + AC
662def dumpbdy (field,n,stag=None,condition=False,onlyx=False,onlyy=False):
663    nx = len(field[0,:])-1
664    ny = len(field[:,0])-1
665    if condition:
666      if stag == 'U': nx = nx-1
667      if stag == 'V': ny = ny-1
668      if stag == 'W': nx = nx+1 #special les case when we dump stag on W
669    if onlyx:    result = field[:,n:nx-n]
670    elif onlyy:  result = field[n:ny-n,:]
671    else:        result = field[n:ny-n,n:nx-n]
672    return result
673
674## Author: AS + AC
675def getcoorddef ( nc ):   
676    import numpy as np
677    ## getcoord2d for predefined types
678    typefile = whatkindfile(nc)
679    if typefile in ['meso']:
680        if '9999' not in getattr(nc,'START_DATE') :   
681            ## regular mesoscale
682            [lon2d,lat2d] = getcoord2d(nc) 
683        else:                     
684            ## idealized mesoscale                     
685            nx=getattr(nc,'WEST-EAST_GRID_DIMENSION')
686            ny=getattr(nc,'SOUTH-NORTH_GRID_DIMENSION')
687            dlat=getattr(nc,'DX')
688            ## this is dirty because Martian-specific
689            # ... but this just intended to get "lat-lon" like info
690            falselon = np.arange(-dlat*(nx-1)/2.,dlat*(nx-1)/2.,dlat)/60000.
691            falselat = np.arange(-dlat*(ny-1)/2.,dlat*(ny-1)/2.,dlat)/60000.
692            [lon2d,lat2d] = np.meshgrid(falselon,falselat) ## dummy coordinates
693            print "WARNING: domain plot artificially centered on lat,lon 0,0"
694    elif typefile in ['gcm','earthgcm','ecmwf']: 
695        #### n est ce pas nc.variables ?
696        if "longitude" in nc.dimensions:  dalon = "longitude"
697        elif "lon" in nc.dimensions:      dalon = "lon"
698        else:                             dalon = "nothing"
699        if "latitude" in nc.dimensions:   dalat = "latitude"
700        elif "lat" in nc.dimensions:      dalat = "lat"
701        else:                             dalat = "nothing"
702        [lon2d,lat2d] = getcoord2d(nc,nlat=dalat,nlon=dalon,is1d=True)
703    elif typefile in ['geo']:
704        [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M')
705    return lon2d,lat2d   
706
707## Author: AS
708def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False):
709    import numpy as np
710    if nlon == "nothing" or nlat == "nothing":
711        print "NO LAT LON FIELDS. I AM TRYING MY BEST. I ASSUME GLOBAL FIELD."
712        lon = np.linspace(-180.,180.,getdimfromvar(nc)[-1])
713        lat = np.linspace(-90.,90.,getdimfromvar(nc)[-2])
714        [lon2d,lat2d] = np.meshgrid(lon,lat)
715    else:
716        if is1d:
717            lat = nc.variables[nlat][:]
718            lon = nc.variables[nlon][:]
719            [lon2d,lat2d] = np.meshgrid(lon,lat)
720        else:
721            lat = nc.variables[nlat][0,:,:]
722            lon = nc.variables[nlon][0,:,:]
723            [lon2d,lat2d] = [lon,lat]
724    return lon2d,lat2d
725
726## Author: AS
727def getdimfromvar (nc):
728    varinfile = nc.variables.keys()
729    dim = nc.variables[varinfile[-1]].shape ## usually the last variable is 4D or 3D
730    return dim
731
732## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
733def smooth1d(x,window_len=11,window='hanning'):
734    import numpy
735    """smooth the data using a window with requested size.
736    This method is based on the convolution of a scaled window with the signal.
737    The signal is prepared by introducing reflected copies of the signal
738    (with the window size) in both ends so that transient parts are minimized
739    in the begining and end part of the output signal.
740    input:
741        x: the input signal
742        window_len: the dimension of the smoothing window; should be an odd integer
743        window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
744            flat window will produce a moving average smoothing.
745    output:
746        the smoothed signal
747    example:
748    t=linspace(-2,2,0.1)
749    x=sin(t)+randn(len(t))*0.1
750    y=smooth(x)
751    see also:
752    numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
753    scipy.signal.lfilter
754    TODO: the window parameter could be the window itself if an array instead of a string   
755    """
756    x = numpy.array(x)
757    if x.ndim != 1:
758        raise ValueError, "smooth only accepts 1 dimension arrays."
759    if x.size < window_len:
760        raise ValueError, "Input vector needs to be bigger than window size."
761    if window_len<3:
762        return x
763    if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
764        raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
765    s=numpy.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
766    #print(len(s))
767    if window == 'flat': #moving average
768        w=numpy.ones(window_len,'d')
769    else:
770        w=eval('numpy.'+window+'(window_len)')
771    y=numpy.convolve(w/w.sum(),s,mode='valid')
772    return y
773
774## Author: AS
775def smooth (field, coeff):
776        ## actually blur_image could work with different coeff on x and y
777        if coeff > 1:   result = blur_image(field,int(coeff))
778        else:           result = field
779        return result
780
781## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
782def gauss_kern(size, sizey=None):
783        import numpy as np
784        # Returns a normalized 2D gauss kernel array for convolutions
785        size = int(size)
786        if not sizey:
787                sizey = size
788        else:
789                sizey = int(sizey)
790        x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
791        g = np.exp(-(x**2/float(size)+y**2/float(sizey)))
792        return g / g.sum()
793
794## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
795def blur_image(im, n, ny=None) :
796        from scipy.signal import convolve
797        # blurs the image by convolving with a gaussian kernel of typical size n.
798        # The optional keyword argument ny allows for a different size in the y direction.
799        g = gauss_kern(n, sizey=ny)
800        improc = convolve(im, g, mode='same')
801        return improc
802
803## Author: AS
804def getwinddef (nc):   
805    ###
806    varinfile = nc.variables.keys()
807    if 'Um' in varinfile:   [uchar,vchar] = ['Um','Vm'] #; print "this is API meso file"
808    elif 'U' in varinfile:  [uchar,vchar] = ['U','V']   #; print "this is RAW meso file"
809    elif 'u' in varinfile:  [uchar,vchar] = ['u','v']   #; print "this is GCM file"
810    elif 'vitu' in varinfile:  [uchar,vchar] = ['vitu','vitv']   #; print "this is GCM v5 file"
811     ### you can add choices here !
812    else:                   [uchar,vchar] = ['not found','not found']
813    ###
814    if uchar in ['U']:         metwind = False ## geometrical (wrt grid)
815    else:                      metwind = True  ## meteorological (zon/mer)
816    if metwind is False:       print "Not using meteorological winds. You trust numerical grid as being (x,y)"
817    ###
818    return uchar,vchar,metwind
819
820## Author: AS
821def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True):
822    ## scale regle la reference du vecteur
823    ## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs.
824    import  matplotlib.pyplot               as plt
825    import  numpy                           as np
826    #posx = np.min(x) - np.std(x) / 10.
827    #posy = np.min(y) - np.std(y) / 10.
828    posx = np.min(x) 
829    posy = np.min(y) - 4.*np.std(y) / 10.
830    u = smooth(u,csmooth)
831    v = smooth(v,csmooth)
832    widthvec = 0.003 #0.005 #0.003
833    q = plt.quiver( x[::stride,::stride],\
834                    y[::stride,::stride],\
835                    u[::stride,::stride],\
836                    v[::stride,::stride],\
837                    angles='xy',color=color,pivot='middle',\
838                    scale=factor,width=widthvec )
839    if color in ['white','yellow']:     kcolor='black'
840    else:                               kcolor=color
841    if key: p = plt.quiverkey(q,posx,posy,scale,\
842                   str(int(scale)),coordinates='data',color=kcolor,labelpos='S',labelsep = 0.03)
843    return 
844
845## Author: AS
846def display (name):
847    from os import system
848    system("display "+name+" > /dev/null 2> /dev/null &")
849    return name
850
851## Author: AS
852def findstep (wlon):
853    steplon = int((wlon[1]-wlon[0])/4.)  #3
854    step = 120.
855    while step > steplon and step > 15. :       step = step / 2.
856    if step <= 15.:
857        while step > steplon and step > 5.  :   step = step - 5.
858    if step <= 5.:
859        while step > steplon and step > 1.  :   step = step - 1.
860    if step <= 1.:
861        step = 1. 
862    return step
863
864## Author: AS
865def define_proj (char,wlon,wlat,back=None,blat=None,blon=None):
866    from    mpl_toolkits.basemap            import Basemap
867    import  numpy                           as np
868    import  matplotlib                      as mpl
869    from mymath import max
870    meanlon = 0.5*(wlon[0]+wlon[1])
871    meanlat = 0.5*(wlat[0]+wlat[1])
872    zewidth = np.abs(wlon[0]-wlon[1])*60000.*np.cos(3.14*meanlat/180.)
873    zeheight = np.abs(wlat[0]-wlat[1])*60000.
874    if blat is None:
875        ortholat=meanlat
876        if   wlat[0] >= 80.:   blat =  -40. 
877        elif wlat[1] <= -80.:  blat = -40.
878        elif wlat[1] >= 0.:    blat = wlat[0] 
879        elif wlat[0] <= 0.:    blat = wlat[1]
880    else:  ortholat=blat
881    if blon is None:  ortholon=meanlon
882    else:             ortholon=blon
883    h = 50.  ## en km
884    radius = 3397200.
885    #print meanlat, meanlon
886    if   char == "cyl":     m = Basemap(rsphere=radius,projection='cyl',\
887                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
888    elif char == "moll":    m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon)
889    elif char == "ortho":   m = Basemap(rsphere=radius,projection='ortho',lon_0=ortholon,lat_0=ortholat)
890    elif char == "lcc":     m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\
891                              width=zewidth,height=zeheight)
892                              #llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
893    elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.)
894    elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=180.)
895    elif char == "nplaea":  m = Basemap(rsphere=radius,projection='nplaea', boundinglat=wlat[0], lon_0=meanlon)
896    elif char == "laea":    m = Basemap(rsphere=radius,projection='laea',lon_0=meanlon,lat_0=meanlat,lat_ts=meanlat,\
897                              width=zewidth,height=zeheight)
898                              #llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
899    elif char == "nsper":   m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.)
900    elif char == "merc":    m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\
901                              llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
902    fontsizemer = int(mpl.rcParams['font.size']*3./4.)
903    if char in ["cyl","lcc","merc","nsper","laea"]:   step = findstep(wlon)
904    else:                                             step = 10.
905    steplon = step*2.
906    zecolor ='grey'
907    zelinewidth = 1
908    zelatmax = 80
909    if meanlat > 75.: zelatmax = 90. ; step = step/2.
910    # to show gcm grid:
911    #zecolor = 'r'
912    #zelinewidth = 1
913    #step = 180./48.
914    #steplon = 360./64.
915    #zelatmax = 90. - step/3
916    if char not in ["moll"]:
917        if wlon[1]-wlon[0] < 2.:  ## LOCAL MODE
918            m.drawmeridians(np.r_[-1.:1.:0.05], labels=[0,0,0,1], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, fmt='%5.2f')
919            m.drawparallels(np.r_[-1.:1.:0.05], labels=[1,0,0,0], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, fmt='%5.2f')
920        else:  ## GLOBAL OR REGIONAL MODE
921            m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, latmax=zelatmax)
922            m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, latmax=zelatmax)
923    if back: m.warpimage(marsmap(back),scale=0.75)
924            #if not back:
925            #    if not var:                                        back = "mola"    ## if no var:         draw mola
926            #    elif typefile in ['mesoapi','meso','geo'] \
927            #       and proj not in ['merc','lcc','nsper','laea']:  back = "molabw"  ## if var but meso:   draw molabw
928            #    else:                                              pass             ## else:              draw None
929    return m
930
931## Author: AS
932#### test temporaire
933def putpoints (map,plot):
934    #### from http://www.scipy.org/Cookbook/Matplotlib/Maps
935    # lat/lon coordinates of five cities.
936    lats = [18.4]
937    lons = [-134.0]
938    points=['Olympus Mons']
939    # compute the native map projection coordinates for cities.
940    x,y = map(lons,lats)
941    # plot filled circles at the locations of the cities.
942    map.plot(x,y,'bo')
943    # plot the names of those five cities.
944    wherept = 0 #1000 #50000
945    for name,xpt,ypt in zip(points,x,y):
946       plot.text(xpt+wherept,ypt+wherept,name)
947    ## le nom ne s'affiche pas...
948    return
949
950## Author: AS
951def calculate_bounds(field,vmin=None,vmax=None):
952    import numpy as np
953    from mymath import max,min,mean
954    ind = np.where(field < 9e+35)
955    fieldcalc = field[ ind ] # la syntaxe compacte ne marche si field est un tuple
956    ###
957    dev = np.std(fieldcalc)*3.0
958    ###
959    if vmin is None:  zevmin = mean(fieldcalc) - dev
960    else:             zevmin = vmin
961    ###
962    if vmax is None:  zevmax = mean(fieldcalc) + dev
963    else:             zevmax = vmax
964    if vmin == vmax:
965                      zevmin = mean(fieldcalc) - dev  ### for continuity
966                      zevmax = mean(fieldcalc) + dev  ### for continuity           
967    ###
968    if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0.
969    print "BOUNDS field ", min(fieldcalc), max(fieldcalc), " //// adopted", zevmin, zevmax
970    return zevmin, zevmax
971
972## Author: AS
973def bounds(what_I_plot,zevmin,zevmax):
974    from mymath import max,min,mean
975    ### might be convenient to add the missing value in arguments
976    #what_I_plot[ what_I_plot < zevmin ] = zevmin#*(1. + 1.e-7)
977    if zevmin < 0: what_I_plot[ what_I_plot < zevmin*(1. - 1.e-7) ] = zevmin*(1. - 1.e-7)
978    else:          what_I_plot[ what_I_plot < zevmin*(1. + 1.e-7) ] = zevmin*(1. + 1.e-7)
979    #print "NEW MIN ", min(what_I_plot)
980    what_I_plot[ what_I_plot > 9e+35  ] = -9e+35
981    what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7)
982    #print "NEW MAX ", max(what_I_plot)
983    return what_I_plot
984
985## Author: AS
986def nolow(what_I_plot):
987    from mymath import max,min
988    lim = 0.15*0.5*(abs(max(what_I_plot))+abs(min(what_I_plot)))
989    print "NO PLOT BELOW VALUE ", lim
990    what_I_plot [ abs(what_I_plot) < lim ] = 1.e40 
991    return what_I_plot
992
993
994## Author : AC
995def hole_bounds(what_I_plot,zevmin,zevmax):
996    import numpy as np
997    zi=0
998    for i in what_I_plot:
999        zj=0
1000        for j in i:
1001            if ((j < zevmin) or (j > zevmax)):what_I_plot[zi,zj]=np.NaN
1002            zj=zj+1
1003        zi=zi+1
1004
1005    return what_I_plot
1006
1007## Author: AS
1008def zoomset (wlon,wlat,zoom):
1009    dlon = abs(wlon[1]-wlon[0])/2.
1010    dlat = abs(wlat[1]-wlat[0])/2.
1011    [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\
1012                    [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ]
1013    print "ZOOM %",zoom,wlon,wlat
1014    return wlon,wlat
1015
1016## Author: AS
1017def fmtvar (whichvar="def"):
1018    fmtvar    =     { \
1019             "MIXED":        "%.0f",\
1020             "UPDRAFT":      "%.0f",\
1021             "DOWNDRAFT":    "%.0f",\
1022             "TK":           "%.0f",\
1023             "T":            "%.0f",\
1024             #"ZMAX_TH":      "%.0f",\
1025             #"WSTAR":        "%.0f",\
1026             # Variables from TES ncdf format
1027             "T_NADIR_DAY":  "%.0f",\
1028             "T_NADIR_NIT":  "%.0f",\
1029             # Variables from tes.py ncdf format
1030             "TEMP_DAY":     "%.0f",\
1031             "TEMP_NIGHT":   "%.0f",\
1032             # Variables from MCS and mcs.py ncdf format
1033             "DTEMP":        "%.0f",\
1034             "NTEMP":        "%.0f",\
1035             "DNUMBINTEMP":  "%.0f",\
1036             "NNUMBINTEMP":  "%.0f",\
1037             # other stuff
1038             "TPOT":         "%.0f",\
1039             "TSURF":        "%.0f",\
1040             "U_OUT1":       "%.0f",\
1041             "T_OUT1":       "%.0f",\
1042             "def":          "%.1e",\
1043             "PTOT":         "%.0f",\
1044             "PSFC":         "%.1f",\
1045             "HGT":          "%.1e",\
1046             "USTM":         "%.2f",\
1047             "HFX":          "%.0f",\
1048             "ICETOT":       "%.1e",\
1049             "TAU_ICE":      "%.2f",\
1050             "TAUICE":       "%.2f",\
1051             "VMR_ICE":      "%.1e",\
1052             "MTOT":         "%.1f",\
1053             "ANOMALY":      "%.1f",\
1054             "W":            "%.2f",\
1055             "WMAX_TH":      "%.1f",\
1056             "WSTAR":        "%.1f",\
1057             "QSURFICE":     "%.0f",\
1058             "UM":           "%.0f",\
1059             "WIND":         "%.0f",\
1060             "UVMET":         "%.0f",\
1061             "UV":         "%.0f",\
1062             "ALBBARE":      "%.2f",\
1063             "TAU":          "%.1f",\
1064             "CO2":          "%.2f",\
1065             "ENFACT":       "%.1f",\
1066             "QDUST":        "%.6f",\
1067             #### T.N.
1068             "TEMP":         "%.0f",\
1069             "VMR_H2OICE":   "%.0f",\
1070             "VMR_H2OVAP":   "%.0f",\
1071             "TAUTES":       "%.2f",\
1072             "TAUTESAP":     "%.2f",\
1073
1074                    }
1075    if "TSURF" in whichvar: whichvar = "TSURF"
1076    if whichvar not in fmtvar:
1077        whichvar = "def"
1078    return fmtvar[whichvar]
1079
1080## Author: AS
1081####################################################################################################################
1082### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
1083def defcolorb (whichone="def"):
1084    whichcolorb =    { \
1085             "def":          "spectral",\
1086             "HGT":          "spectral",\
1087             "HGT_M":        "spectral",\
1088             "TK":           "gist_heat",\
1089             "TPOT":         "Paired",\
1090             "TSURF":        "RdBu_r",\
1091             "QH2O":         "PuBu",\
1092             "USTM":         "YlOrRd",\
1093             "WIND":         "YlOrRd",\
1094             #"T_nadir_nit":  "RdBu_r",\
1095             #"T_nadir_day":  "RdBu_r",\
1096             "HFX":          "RdYlBu",\
1097             "ICETOT":       "YlGnBu_r",\
1098             #"MTOT":         "PuBu",\
1099             "CCNQ":         "YlOrBr",\
1100             "CCNN":         "YlOrBr",\
1101             "TEMP":         "Jet",\
1102             "TAU_ICE":      "Blues",\
1103             "TAUICE":       "Blues",\
1104             "VMR_ICE":      "Blues",\
1105             "W":            "jet",\
1106             "WMAX_TH":      "spectral",\
1107             "ANOMALY":      "RdBu_r",\
1108             "QSURFICE":     "hot_r",\
1109             "ALBBARE":      "spectral",\
1110             "TAU":          "YlOrBr_r",\
1111             "CO2":          "YlOrBr_r",\
1112             "MIXED":        "GnBu",\
1113             #### T.N.
1114             "MTOT":         "spectral",\
1115             "H2O_ICE_S":    "RdBu",\
1116             "VMR_H2OICE":   "PuBu",\
1117             "VMR_H2OVAP":   "PuBu",\
1118             "WATERCAPTAG":  "Blues",\
1119                     }
1120#W --> spectral ou jet
1121#spectral BrBG RdBu_r
1122    #print "predefined colorbars"
1123    if "TSURF" in whichone: whichone = "TSURF"
1124    if whichone not in whichcolorb:
1125        whichone = "def"
1126    return whichcolorb[whichone]
1127
1128## Author: AS
1129def definecolorvec (whichone="def"):
1130        whichcolor =    { \
1131                "def":          "black",\
1132                "vis":          "yellow",\
1133                "vishires":     "green",\
1134                "molabw":       "yellow",\
1135                "mola":         "black",\
1136                "gist_heat":    "white",\
1137                "hot":          "tk",\
1138                "gist_rainbow": "black",\
1139                "spectral":     "black",\
1140                "gray":         "red",\
1141                "PuBu":         "black",\
1142                "titan":        "red",\
1143                        }
1144        if whichone not in whichcolor:
1145                whichone = "def"
1146        return whichcolor[whichone]
1147
1148## Author: AS
1149def marsmap (whichone="vishires"):
1150        from os import uname
1151        mymachine = uname()[1]
1152        ### not sure about speed-up with this method... looks the same
1153        if "lmd.jussieu.fr" in mymachine:   domain = "/u/aslmd/WWW/maps/"
1154        elif "aymeric-laptop" in mymachine: domain = "/home/aymeric/Dropbox/Public/"
1155        else:                               domain = "http://www.lmd.jussieu.fr/~aslmd/maps/"
1156        whichlink =     { \
1157                #"vis":         "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\
1158                #"vishires":    "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\
1159                #"geolocal":    "http://dl.dropbox.com/u/11078310/geolocal.jpg",\
1160                #"mola":        "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\
1161                #"molabw":      "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\
1162                "thermalday":   domain+"thermalday.jpg",\
1163                "thermalnight": domain+"thermalnight.jpg",\
1164                "tesalbedo":    domain+"tesalbedo.jpg",\
1165                "vis":         domain+"mar0kuu2.jpg",\
1166                "vishires":    domain+"MarsMap_2500x1250.jpg",\
1167                "geolocal":    domain+"geolocal.jpg",\
1168                "mola":        domain+"mars-mola-2k.jpg",\
1169                "molabw":      domain+"MarsElevation_2500x1250.jpg",\
1170                "clouds":      "http://www.johnstonsarchive.net/spaceart/marswcloudmap.jpg",\
1171                "jupiter":     "http://www.mmedia.is/~bjj/data/jupiter_css/jupiter_css.jpg",\
1172                "jupiter_voy": "http://www.mmedia.is/~bjj/data/jupiter/jupiter_vgr2.jpg",\
1173                #"bw":          domain+"EarthElevation_2500x1250.jpg",\
1174                "bw":          "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg",\
1175                "contrast":    "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg",\
1176                "nice":        "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg",\
1177                "blue":        "http://eoimages.gsfc.nasa.gov/ve/2430/land_ocean_ice_2048.jpg",\
1178                "blueclouds":  "http://eoimages.gsfc.nasa.gov/ve/2431/land_ocean_ice_cloud_2048.jpg",\
1179                "justclouds":  "http://eoimages.gsfc.nasa.gov/ve/2432/cloud_combined_2048.jpg",\
1180                "pluto":       "http://www.boulder.swri.edu/~buie/pluto/pluto_all.png",\
1181                "triton":      "http://laps.noaa.gov/albers/sos/neptune/triton/triton_rgb_cyl_www.jpg",\
1182                "titan":       "http://laps.noaa.gov/albers/sos/saturn/titan/titan_rgb_cyl_www.jpg",\
1183                #"titan":       "http://laps.noaa.gov/albers/sos/celestia/titan_50.jpg",\
1184                "titanuni":    "http://maps.jpl.nasa.gov/pix/sat6fss1.jpg",\
1185                "venus":       "http://laps.noaa.gov/albers/sos/venus/venus4/venus4_rgb_cyl_www.jpg",\
1186                "cosmic":      "http://laps.noaa.gov/albers/sos/universe/wmap/wmap_rgb_cyl_www.jpg",\
1187                        }
1188        ### see http://www.mmedia.is/~bjj/planetary_maps.html
1189        if whichone not in whichlink: 
1190                print "marsmap: choice not defined... you'll get the default one... "
1191                whichone = "vishires" 
1192        return whichlink[whichone]
1193
1194#def earthmap (whichone):
1195#       if   whichone == "contrast":    whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg"
1196#       elif whichone == "bw":          whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg"
1197#       elif whichone == "nice":        whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg"
1198#       return whichlink
1199
1200## Author: AS
1201def latinterv (area="Whole"):
1202    list =    { \
1203        "Europe":                [[ 20., 80.],[- 50.,  50.]],\
1204        "Central_America":       [[-10., 40.],[ 230., 300.]],\
1205        "Africa":                [[-20., 50.],[- 50.,  50.]],\
1206        "Whole":                 [[-90., 90.],[-180., 180.]],\
1207        "Southern_Hemisphere":   [[-90., 60.],[-180., 180.]],\
1208        "Northern_Hemisphere":   [[-60., 90.],[-180., 180.]],\
1209        "Tharsis":               [[-30., 60.],[-170.,- 10.]],\
1210        "Whole_No_High":         [[-60., 60.],[-180., 180.]],\
1211        "Chryse":                [[-60., 60.],[- 60.,  60.]],\
1212        "North_Pole":            [[ 50., 90.],[-180., 180.]],\
1213        "Close_North_Pole":      [[ 75., 90.],[-180., 180.]],\
1214        "Far_South_Pole":        [[-90.,-40.],[-180., 180.]],\
1215        "South_Pole":            [[-90.,-50.],[-180., 180.]],\
1216        "Close_South_Pole":      [[-90.,-75.],[-180., 180.]],\
1217        "Sirenum_Crater_large":  [[-46.,-34.],[-166.,-151.]],\
1218        "Sirenum_Crater_small":  [[-36.,-26.],[-168.,-156.]],\
1219        "Rupes":                 [[ 72., 90.],[-120.,- 20.]],\
1220        "Xanadu":                [[-40., 20.],[  40., 120.]],\
1221        "Hyperboreae":           [[ 80., 87.],[- 70.,- 10.]],\
1222              }
1223    if area not in list:   area = "Whole"
1224    [olat,olon] = list[area]
1225    return olon,olat
1226
1227## Author: TN
1228def separatenames (name):
1229  from numpy import concatenate
1230  # look for comas in the input name to separate different names (files, variables,etc ..)
1231  if name is None:
1232     names = None
1233  else:
1234    names = []
1235    stop = 0
1236    currentname = name
1237    while stop == 0:
1238      indexvir = currentname.find(',')
1239      if indexvir == -1:
1240        stop = 1
1241        name1 = currentname
1242      else:
1243        name1 = currentname[0:indexvir]
1244      names = concatenate((names,[name1]))
1245      currentname = currentname[indexvir+1:len(currentname)]
1246  return names
1247
1248
1249## Author: TN
1250def readslices(saxis):
1251  from numpy import empty
1252  if saxis == None:
1253     zesaxis = None
1254  else:
1255     zesaxis = empty((len(saxis),2))
1256     for i in range(len(saxis)):
1257        a = separatenames(saxis[i])
1258        if len(a) == 1:
1259           zesaxis[i,:] = float(a[0])
1260        else:
1261           zesaxis[i,0] = float(a[0])
1262           zesaxis[i,1] = float(a[1])
1263           
1264  return zesaxis
1265
1266## Author: TN
1267def readdata(data,datatype,coord1,coord2):
1268## Read sparse data
1269  from numpy import empty
1270  if datatype == 'txt':
1271     if len(data[coord1].shape) == 1:
1272         return data[coord1][:]
1273     elif len(data[coord1].shape) == 2:
1274         return data[coord1][:,int(coord2)-1]
1275     else:
1276         errormess('error in readdata')
1277  elif datatype == 'sav':
1278     return data[coord1][coord2]
1279  else:
1280     errormess(datatype+' type is not supported!')
1281
1282
1283## Author: AS
1284def bidimfind(lon2d,lat2d,vlon,vlat,file=None):
1285   import numpy as np
1286   import matplotlib.pyplot as mpl
1287   if vlat is None:    array = (lon2d - vlon)**2
1288   elif vlon is None:  array = (lat2d - vlat)**2
1289   else:               array = (lon2d - vlon)**2 + (lat2d - vlat)**2
1290   idy,idx = np.unravel_index( np.argmin(array), lon2d.shape )
1291   if vlon is not None:
1292      if (np.abs(lon2d[idy,idx]-vlon)) > 5: errormess("longitude not found ",printvar=lon2d)
1293   if vlat is not None:
1294      if (np.abs(lat2d[idy,idx]-vlat)) > 5: errormess("latitude not found ",printvar=lat2d)
1295   if file is not None:
1296      print idx,idy,lon2d[idy,idx],vlon
1297      print idx,idy,lat2d[idy,idx],vlat
1298      var = file.variables["HGT"][:,:,:]
1299      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)
1300      mpl.show()
1301   return idy,idx
1302
1303## Author: TN
1304def getsindex(saxis,index,axis):
1305# input  : all the desired slices and the good index
1306# output : all indexes to be taken into account for reducing field
1307  import numpy as np
1308  if ( np.array(axis).ndim == 2):
1309      axis = axis[:,0]
1310  if saxis is None:
1311      zeindex = None
1312  else:
1313      aaa = int(np.argmin(abs(saxis[index,0] - axis)))
1314      bbb = int(np.argmin(abs(saxis[index,1] - axis)))
1315      [imin,imax] = np.sort(np.array([aaa,bbb]))
1316      zeindex = np.array(range(imax-imin+1))+imin
1317      # because -180 and 180 are the same point in longitude,
1318      # we get rid of one for averaging purposes.
1319      if axis[imin] == -180 and axis[imax] == 180:
1320         zeindex = zeindex[0:len(zeindex)-1]
1321         print "INFO: whole longitude averaging asked, so last point is not taken into account."
1322  return zeindex
1323     
1324## Author: TN
1325def define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,indextime,what_I_plot,dim0,vertmode,redope):
1326# Purpose of define_axis is to find x and y axis scales in a smart way
1327# x axis priority: 1/time 2/lon 3/lat 4/vertical
1328# To be improved !!!...
1329   from numpy import array,swapaxes
1330   x = None
1331   y = None
1332   count = 0
1333   what_I_plot = array(what_I_plot)
1334   shape = what_I_plot.shape
1335   if indextime is None and len(time) > 1:
1336      print "AXIS is time"
1337      x = time
1338      count = count+1
1339   if indexlon is None and len(lon) > 1 and redope not in ['edge_x1','edge_x2']:
1340      print "AXIS is lon"
1341      if count == 0: x = lon
1342      else: y = lon
1343      count = count+1
1344   if indexlat is None and len(lat) > 1 and redope not in ['edge_y1','edge_y2']:
1345      print "AXIS is lat"
1346      if count == 0: x = lat
1347      else: y = lat
1348      count = count+1
1349   if indexvert is None and ((dim0 == 4) or (y is None)):
1350      print "AXIS is vert"
1351      if vertmode == 0: # vertical axis is as is (GCM grid)
1352         if count == 0: x=range(len(vert))
1353         else: y=range(len(vert))
1354         count = count+1
1355      else: # vertical axis is in kms
1356         if count == 0: x = vert
1357         else: y = vert
1358         count = count+1
1359   x = array(x)
1360   y = array(y)
1361   print "CHECK SHAPE: what_I_plot, x, y", what_I_plot.shape, x.shape, y.shape
1362   if len(shape) == 1:
1363       if shape[0] != len(x):           print "WARNING: shape[0] != len(x). Correcting." ; what_I_plot = what_I_plot[0:len(x)]
1364       if len(y.shape) > 0:             y = ()
1365   elif len(shape) == 2:
1366       if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]:
1367           print "INFO: swapaxes: ",what_I_plot.shape,shape ; what_I_plot = swapaxes(what_I_plot,0,1)
1368       else:
1369           if shape[0] != len(y):       print "WARNING: shape[0] != len(y). Correcting." ; what_I_plot = what_I_plot[0:len(y),:]
1370           elif shape[1] != len(x):     print "WARNING: shape[1] != len(x). Correcting." ; what_I_plot = what_I_plot[:,0:len(x)]
1371   elif len(shape) == 3:
1372       if vertmode < 0: print "not supported. must check array dimensions at some point. not difficult to implement though."
1373   return what_I_plot,x,y
1374
1375# Author: TN + AS + AC
1376def determineplot(slon, slat, svert, stime, redope):
1377    nlon = 1 # number of longitudinal slices -- 1 is None
1378    nlat = 1
1379    nvert = 1
1380    ntime = 1
1381    nslices = 1
1382    if slon is not None:
1383        length=len(slon[:,0])
1384        nslices = nslices*length
1385        nlon = len(slon)
1386    if slat is not None:
1387        length=len(slat[:,0])
1388        nslices = nslices*length
1389        nlat = len(slat)
1390    if svert is not None:
1391        length=len(svert[:,0])
1392        nslices = nslices*length
1393        nvert = len(svert)
1394    if stime is not None:
1395        length=len(stime[:,0])
1396        nslices = nslices*length
1397        ntime = len(stime)
1398    #else:
1399    #    nslices = 2 
1400    mapmode = 0
1401    if slon is None and slat is None and redope not in ['edge_x1','edge_x2','edge_y1','edge_y2']:
1402       mapmode = 1 # in this case we plot a map, with the given projection
1403    return nlon, nlat, nvert, ntime, mapmode, nslices
1404
1405## Author : AS
1406def maplatlon( lon,lat,field,\
1407               proj="cyl",colorb="jet",ndiv=10,zeback="molabw",trans=0.6,title="",\
1408               vecx=None,vecy=None,stride=2 ):
1409    ### an easy way to map a field over lat/lon grid
1410    import numpy as np
1411    import matplotlib.pyplot as mpl
1412    from matplotlib.cm import get_cmap
1413    ## get lon and lat in 2D version. get lat/lon intervals
1414    numdim = len(np.array(lon).shape)
1415    if numdim == 2:     [lon2d,lat2d] = [lon,lat]
1416    elif numdim == 1:   [lon2d,lat2d] = np.meshgrid(lon,lat)
1417    else:               errormess("lon and lat arrays must be 1D or 2D")
1418    [wlon,wlat] = latinterv()
1419    ## define projection and background. define x and y given the projection
1420    m = define_proj(proj,wlon,wlat,back=zeback,blat=None,blon=None)
1421    x, y = m(lon2d, lat2d)
1422    ## define field. bound field.
1423    what_I_plot = np.transpose(field)
1424    zevmin, zevmax = calculate_bounds(what_I_plot)  ## vmin=min(what_I_plot_frame), vmax=max(what_I_plot_frame))
1425    what_I_plot = bounds(what_I_plot,zevmin,zevmax)
1426    ## define contour field levels. define color palette
1427    ticks = ndiv + 1
1428    zelevels = np.linspace(zevmin,zevmax,ticks)
1429    palette = get_cmap(name=colorb)
1430    ## contour field
1431    m.contourf( x, y, what_I_plot, zelevels, cmap = palette, alpha = trans )
1432    ## draw colorbar
1433    if proj in ['moll','cyl']:        zeorientation="horizontal" ; zepad = 0.07
1434    else:                             zeorientation="vertical" ; zepad = 0.03
1435    #daformat = fmtvar(fvar.upper())
1436    daformat = "%.0f"
1437    zecb = mpl.colorbar( fraction=0.05,pad=zepad,format=daformat,orientation=zeorientation,\
1438                 ticks=np.linspace(zevmin,zevmax,num=min([ticks/2+1,21])),extend='neither',spacing='proportional' ) 
1439    ## give a title
1440    if zeorientation == "horizontal": zecb.ax.set_xlabel(title)
1441    else:                             ptitle(title)
1442    ## draw vector
1443    if vecx is not None and vecy is not None:
1444       [vecx_frame,vecy_frame] = m.rotate_vector( np.transpose(vecx), np.transpose(vecy), lon2d, lat2d ) ## for metwinds
1445       vectorfield(vecx_frame, vecy_frame, x, y, stride=stride, csmooth=2,\
1446                                             scale=30., factor=500., color=definecolorvec(colorb), key=True)
1447    ## scale regle la reference du vecteur. factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs.
1448    return
1449## Author : AC
1450## Handles calls to specific computations (e.g. wind norm, enrichment factor...)
1451def select_getfield(zvarname=None,znc=None,ztypefile=None,mode=None,ztsat=None,ylon=None,ylat=None,yalt=None,ytime=None,analysis=None):
1452      from mymath import get_tsat
1453 
1454      ## Specific variables are described here:
1455      # for the mesoscale:
1456      specificname_meso = ['UV','uv','uvmet','slopexy','SLOPEXY','deltat','DELTAT','hodograph','tk','hodograph_2']
1457      # for the gcm:
1458      specificname_gcm = ['enfact']
1459
1460      ## Check for variable in file:
1461      if mode == 'check':     
1462          varname = zvarname
1463          varinfile=znc.variables.keys()
1464          logical_novarname = zvarname not in znc.variables
1465          logical_nospecificname_meso = not ((ztypefile in ['meso']) and (zvarname in specificname_meso))
1466          logical_nospecificname_gcm = not ((ztypefile in ['gcm']) and (zvarname in specificname_gcm))
1467          if ( logical_novarname and logical_nospecificname_meso and logical_nospecificname_gcm ):
1468              if len(varinfile) == 1:   varname = varinfile[0]
1469              else:                     varname = False
1470          ## Return the variable name:
1471          return varname
1472
1473      ## Get the corresponding variable:
1474      if mode == 'getvar':
1475          plot_x = None ; plot_y = None ;
1476          ### ----------- 1. saturation temperature
1477          if zvarname in ["temp","t","T_nadir_nit","T_nadir_day","temp_day","temp_night"] and ztsat:
1478              tt=getfield(znc,zvarname) ; print "computing Tsat-T, I ASSUME Z-AXIS IS PRESSURE"
1479              if type(tt).__name__=='MaskedArray':  tt.set_fill_value([np.NaN]) ; tinput=tt.filled()
1480              else:                                 tinput=tt
1481              all_var=get_tsat(yalt,tinput,zlon=ylon,zlat=ylat,zalt=yalt,ztime=ytime)
1482          ### ----------- 2. wind amplitude
1483          elif ((zvarname in ['UV','uv','uvmet']) and (ztypefile in ['meso']) and (zvarname not in znc.variables)):
1484              all_var=windamplitude(znc,'amplitude')
1485          elif ((zvarname in ['hodograph','hodograph_2']) and (ztypefile in ['meso']) and (zvarname not in znc.variables)):
1486              plot_x, plot_y = windamplitude(znc,zvarname)
1487              if plot_x is not None: all_var=plot_x # dummy
1488              else: all_var=plot_y ; plot_x = None ; plot_y = None # Hodograph type 2 is not 'xy' mode
1489          elif ((zvarname in ['slopexy','SLOPEXY']) and (ztypefile in ['meso']) and (zvarname not in znc.variables)):
1490              all_var=slopeamplitude(znc)
1491          ### ------------ 3. Near surface instability
1492          elif ((zvarname in ['DELTAT','deltat']) and (ztypefile in ['meso']) and (zvarname not in znc.variables)):
1493              all_var=deltat0t1(znc)
1494          ### ------------ 4. Enrichment factor
1495          elif ((ztypefile in ['gcm']) and (zvarname in ['enfact'])):
1496              all_var=enrichment_factor(znc,ylon,ylat,ytime)
1497          ### ------------ 5. teta -> temp
1498          elif ((ztypefile in ['meso']) and (zvarname in ['tk']) and ('tk' not in znc.variables.keys())):
1499              all_var=teta_to_tk(znc)
1500          else:
1501          ### -----------  999. Normal case
1502              all_var = getfield(znc,zvarname)
1503          if analysis is not None:
1504             if analysis in ['histo','density','histodensity']: plot_y=all_var ; plot_x = plot_y
1505             elif analysis == 'fft': plot_y, plot_x = spectrum(all_var,ytime,yalt,ylat,ylon) ; all_var = plot_y
1506          return all_var, plot_x, plot_y
1507
1508# Author : A.C
1509# FFT is computed before reducefield voluntarily, because we dont want to compute
1510# ffts on averaged fields (which would kill all waves). Instead, we take the fft everywhere
1511# (which is not efficient but it is still ok) and then, make the average (if the user wants to)
1512def spectrum(var,time,vert,lat,lon):
1513    import numpy as np
1514    fft=np.fft.fft(var,axis=1)
1515    N=len(vert)
1516    step=(vert[1]-vert[0])*1000.
1517    print "step is: ",step
1518    fftfreq=np.fft.fftfreq(N,d=step)
1519    fftfreq=np.fft.fftshift(fftfreq) # spatial FFT => this is the wavenumber
1520    fft=np.fft.fftshift(fft)
1521    fftfreq = 1./fftfreq # => wavelength (div by 0 expected, don't panic)
1522    fft=np.abs(fft) # => amplitude spectrum
1523#    fft=np.abs(fft)**2 # => power spectrum
1524    return fft,fftfreq
1525
1526# Author : A.C.
1527# Computes temperature from potential temperature for mesoscale files, without the need to use API, i.e. using natural vertical grid
1528def teta_to_tk(nc):
1529    import numpy as np
1530    varinfile = nc.variables.keys() 
1531    p0=610.
1532    t0=220.
1533    r_cp=1./3.89419
1534    if "T" in varinfile: zteta=getfield(nc,'T')
1535    else: errormess("you need T in your file.")
1536    if "PTOT" in varinfile: zptot=getfield(nc,'PTOT')
1537    else: errormess("you need PTOT in your file.")
1538    zt=(zteta+220.)*(zptot/p0)**(r_cp)
1539    return zt
Note: See TracBrowser for help on using the repository browser.