source: trunk/UTIL/PYTHON/mcd/myplot.py @ 937

Last change on this file since 937 was 937, checked in by aslmd, 12 years ago

UTIL PYTHON MCDonline. small bug fixes and mcd5

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