Context Navigation

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

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

GRAPHICS: PYTHON:
updated movie capabilities with HTML support.
corrected problems with cat, winds, ...
added avi and html possibility to "save" keyword.

-S avi is simply equivalent to --rate 8.
-S html creates automatically webpage with nice animations and controls.

2 examples are commented in README.PP :
http://www.lmd.jussieu.fr/~aslmd/EXAMPLES/LMD_GCM_movie_tsurf_UV/anim.html
http://www.lmd.jussieu.fr/~aslmd/EXAMPLES/LMD_MMM_d1_10km_movie_QDUST_-1000m-AMR_lat_-3_Ls134.8/anim.html

File size: 43.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: