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

Last change on this file since 610 was 610, checked in by aslmd, 13 years ago

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