Context Navigation

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

Last change on this file since 453 was 453, checked in by tnavarro, 13 years ago
GRAPHICS: transparent fields and 3 new background images
File size: 43.5 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 = 2
199	subh = 1
200	fig.subplots_adjust(wspace = 0.35)
201	rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
202	elif numplot == 3:
203	subv = 3
204	subh = 1
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, 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	zecolor ='grey'
589	zelinewidth = 1
590	zelatmax = 80
591	# to show gcm grid:
592	#zecolor = 'r'
593	#zelinewidth = 1
594	#step = 5.625
595	#steplon = 5.625
596	#zelatmax = 89.9
597	m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, latmax=zelatmax)
598	m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, latmax=zelatmax)
599	if back: m.warpimage(marsmap(back),scale=0.75)
600	#if not back:
601	# if not var: back = "mola" ## if no var: draw mola
602	# elif typefile in ['mesoapi','meso','geo'] \
603	# and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw
604	# else: pass ## else: draw None
605	return m
606
607	## Author: AS
608	#### test temporaire
609	def putpoints (map,plot):
610	#### from http://www.scipy.org/Cookbook/Matplotlib/Maps
611	# lat/lon coordinates of five cities.
612	lats = [18.4]
613	lons = [-134.0]
614	points=['Olympus Mons']
615	# compute the native map projection coordinates for cities.
616	x,y = map(lons,lats)
617	# plot filled circles at the locations of the cities.
618	map.plot(x,y,'bo')
619	# plot the names of those five cities.
620	wherept = 0 #1000 #50000
621	for name,xpt,ypt in zip(points,x,y):
622	plot.text(xpt+wherept,ypt+wherept,name)
623	## le nom ne s'affiche pas...
624	return
625
626	## Author: AS
627	def calculate_bounds(field,vmin=None,vmax=None):
628	import numpy as np
629	from mymath import max,min,mean
630	ind = np.where(field < 9e+35)
631	fieldcalc = field[ ind ] # la syntaxe compacte ne marche si field est un tuple
632	###
633	dev = np.std(fieldcalc)*3.0
634	###
635	if vmin is None:
636	zevmin = mean(fieldcalc) - dev
637	else: zevmin = vmin
638	###
639	if vmax is None: zevmax = mean(fieldcalc) + dev
640	else: zevmax = vmax
641	if vmin == vmax:
642	zevmin = mean(fieldcalc) - dev ### for continuity
643	zevmax = mean(fieldcalc) + dev ### for continuity
644	###
645	if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0.
646	print "BOUNDS field ", min(fieldcalc), max(fieldcalc)
647	print "BOUNDS adopted ", zevmin, zevmax
648	return zevmin, zevmax
649
650	## Author: AS
651	def bounds(what_I_plot,zevmin,zevmax):
652	from mymath import max,min,mean
653	### might be convenient to add the missing value in arguments
654	#what_I_plot[ what_I_plot < zevmin ] = zevmin#*(1. + 1.e-7)
655	if zevmin < 0: what_I_plot[ what_I_plot < zevmin(1. - 1.e-7) ] = zevmin(1. - 1.e-7)
656	else: what_I_plot[ what_I_plot < zevmin(1. + 1.e-7) ] = zevmin(1. + 1.e-7)
657	#print "NEW MIN ", min(what_I_plot)
658	what_I_plot[ what_I_plot > 9e+35 ] = -9e+35
659	what_I_plot[ what_I_plot > zevmax ] = zevmax
660	#print "NEW MAX ", max(what_I_plot)
661	return what_I_plot
662
663	## Author: AS
664	def nolow(what_I_plot):
665	from mymath import max,min
666	lim = 0.150.5(abs(max(what_I_plot))+abs(min(what_I_plot)))
667	print "NO PLOT BELOW VALUE ", lim
668	what_I_plot [ abs(what_I_plot) < lim ] = 1.e40
669	return what_I_plot
670
671
672	## Author : AC
673	def hole_bounds(what_I_plot,zevmin,zevmax):
674	import numpy as np
675	zi=0
676	for i in what_I_plot:
677	zj=0
678	for j in i:
679	if ((j < zevmin) or (j > zevmax)):what_I_plot[zi,zj]=np.NaN
680	zj=zj+1
681	zi=zi+1
682
683	return what_I_plot
684
685	## Author: AS
686	def zoomset (wlon,wlat,zoom):
687	dlon = abs(wlon[1]-wlon[0])/2.
688	dlat = abs(wlat[1]-wlat[0])/2.
689	[wlon,wlat] = [ [wlon[0]+zoomdlon/100.,wlon[1]-zoomdlon/100.],\
690	[wlat[0]+zoomdlat/100.,wlat[1]-zoomdlat/100.] ]
691	print "ZOOM %",zoom,wlon,wlat
692	return wlon,wlat
693
694	## Author: AS
695	def fmtvar (whichvar="def"):
696	fmtvar = { \
697	"TK": "%.0f",\
698	# Variables from TES ncdf format
699	"T_NADIR_DAY": "%.0f",\
700	"T_NADIR_NIT": "%.0f",\
701	# Variables from tes.py ncdf format
702	"TEMP_DAY": "%.0f",\
703	"TEMP_NIGHT": "%.0f",\
704	# Variables from MCS and mcs.py ncdf format
705	"DTEMP": "%.0f",\
706	"NTEMP": "%.0f",\
707	"DNUMBINTEMP": "%.0f",\
708	"NNUMBINTEMP": "%.0f",\
709	# other stuff
710	"TPOT": "%.0f",\
711	"TSURF": "%.0f",\
712	"def": "%.1e",\
713	"PTOT": "%.0f",\
714	"HGT": "%.1e",\
715	"USTM": "%.2f",\
716	"HFX": "%.0f",\
717	"ICETOT": "%.1e",\
718	"TAU_ICE": "%.2f",\
719	"TAUICE": "%.2f",\
720	"VMR_ICE": "%.1e",\
721	"MTOT": "%.1f",\
722	"ANOMALY": "%.1f",\
723	"W": "%.1f",\
724	"WMAX_TH": "%.1f",\
725	"QSURFICE": "%.0f",\
726	"UM": "%.0f",\
727	"ALBBARE": "%.2f",\
728	"TAU": "%.1f",\
729	"CO2": "%.2f",\
730	#### T.N.
731	"TEMP": "%.0f",\
732	"VMR_H2OICE": "%.0f",\
733	"VMR_H2OVAP": "%.0f",\
734	"TAUTES": "%.2f",\
735	"TAUTESAP": "%.2f",\
736
737	}
738	if whichvar not in fmtvar:
739	whichvar = "def"
740	return fmtvar[whichvar]
741
742	## Author: AS
743	####################################################################################################################
744	### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
745	def defcolorb (whichone="def"):
746	whichcolorb = { \
747	"def": "spectral",\
748	"HGT": "spectral",\
749	"HGT_M": "spectral",\
750	"TK": "gist_heat",\
751	"TPOT": "Paired",\
752	"TSURF": "RdBu_r",\
753	"QH2O": "PuBu",\
754	"USTM": "YlOrRd",\
755	#"T_nadir_nit": "RdBu_r",\
756	#"T_nadir_day": "RdBu_r",\
757	"HFX": "RdYlBu",\
758	"ICETOT": "YlGnBu_r",\
759	#"MTOT": "PuBu",\
760	"CCNQ": "YlOrBr",\
761	"CCNN": "YlOrBr",\
762	"TEMP": "Jet",\
763	"TAU_ICE": "Blues",\
764	"TAUICE": "Blues",\
765	"VMR_ICE": "Blues",\
766	"W": "jet",\
767	"WMAX_TH": "spectral",\
768	"ANOMALY": "RdBu_r",\
769	"QSURFICE": "hot_r",\
770	"ALBBARE": "spectral",\
771	"TAU": "YlOrBr_r",\
772	"CO2": "YlOrBr_r",\
773	#### T.N.
774	"MTOT": "Jet",\
775	"H2O_ICE_S": "RdBu",\
776	"VMR_H2OICE": "PuBu",\
777	"VMR_H2OVAP": "PuBu",\
778	"WATERCAPTAG": "Blues",\
779	}
780	#W --> spectral ou jet
781	#spectral BrBG RdBu_r
782	#print "predefined colorbars"
783	if whichone not in whichcolorb:
784	whichone = "def"
785	return whichcolorb[whichone]
786
787	## Author: AS
788	def definecolorvec (whichone="def"):
789	whichcolor = { \
790	"def": "black",\
791	"vis": "yellow",\
792	"vishires": "yellow",\
793	"molabw": "yellow",\
794	"mola": "black",\
795	"gist_heat": "white",\
796	"hot": "tk",\
797	"gist_rainbow": "black",\
798	"spectral": "black",\
799	"gray": "red",\
800	"PuBu": "black",\
801	}
802	if whichone not in whichcolor:
803	whichone = "def"
804	return whichcolor[whichone]
805
806	## Author: AS
807	def marsmap (whichone="vishires"):
808	from os import uname
809	mymachine = uname()[1]
810	### not sure about speed-up with this method... looks the same
811	if "lmd.jussieu.fr" in mymachine: domain = "/u/aslmd/WWW/maps/"
812	else: domain = "http://www.lmd.jussieu.fr/~aslmd/maps/"
813	whichlink = { \
814	#"vis": "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\
815	#"vishires": "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\
816	#"geolocal": "http://dl.dropbox.com/u/11078310/geolocal.jpg",\
817	#"mola": "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\
818	#"molabw": "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\
819	"thermalday": domain+"thermalday.jpg",\
820	"thermalnight": domain+"thermalnight.jpg",\
821	"tesalbedo": domain+"tesalbedo.jpg",\
822	"vis": domain+"mar0kuu2.jpg",\
823	"vishires": domain+"MarsMap_2500x1250.jpg",\
824	"geolocal": domain+"geolocal.jpg",\
825	"mola": domain+"mars-mola-2k.jpg",\
826	"molabw": domain+"MarsElevation_2500x1250.jpg",\
827	"clouds": "http://www.johnstonsarchive.net/spaceart/marswcloudmap.jpg",\
828	"jupiter": "http://www.mmedia.is/~bjj/data/jupiter_css/jupiter_css.jpg",\
829	"jupiter_voy": "http://www.mmedia.is/~bjj/data/jupiter/jupiter_vgr2.jpg",\
830	"bw": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg",\
831	"contrast": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg",\
832	"nice": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg",\
833	"blue": "http://eoimages.gsfc.nasa.gov/ve/2430/land_ocean_ice_2048.jpg",\
834	"blueclouds": "http://eoimages.gsfc.nasa.gov/ve/2431/land_ocean_ice_cloud_2048.jpg",\
835	"justclouds": "http://eoimages.gsfc.nasa.gov/ve/2432/cloud_combined_2048.jpg",\
836	}
837	### see http://www.mmedia.is/~bjj/planetary_maps.html
838	if whichone not in whichlink:
839	print "marsmap: choice not defined... you'll get the default one... "
840	whichone = "vishires"
841	return whichlink[whichone]
842
843	#def earthmap (whichone):
844	# if whichone == "contrast": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg"
845	# elif whichone == "bw": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg"
846	# elif whichone == "nice": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg"
847	# return whichlink
848
849	## Author: AS
850	def latinterv (area="Whole"):
851	list = { \
852	"Europe": [[ 20., 80.],[- 50., 50.]],\
853	"Central_America": [[-10., 40.],[ 230., 300.]],\
854	"Africa": [[-20., 50.],[- 50., 50.]],\
855	"Whole": [[-90., 90.],[-180., 180.]],\
856	"Southern_Hemisphere": [[-90., 60.],[-180., 180.]],\
857	"Northern_Hemisphere": [[-60., 90.],[-180., 180.]],\
858	"Tharsis": [[-30., 60.],[-170.,- 10.]],\
859	"Whole_No_High": [[-60., 60.],[-180., 180.]],\
860	"Chryse": [[-60., 60.],[- 60., 60.]],\
861	"North_Pole": [[ 50., 90.],[-180., 180.]],\
862	"Close_North_Pole": [[ 75., 90.],[-180., 180.]],\
863	"Far_South_Pole": [[-90.,-40.],[-180., 180.]],\
864	"South_Pole": [[-90.,-50.],[-180., 180.]],\
865	"Close_South_Pole": [[-90.,-75.],[-180., 180.]],\
866	}
867	if area not in list: area = "Whole"
868	[olat,olon] = list[area]
869	return olon,olat
870
871	## Author: TN
872	def separatenames (name):
873	from numpy import concatenate
874	# look for comas in the input name to separate different names (files, variables,etc ..)
875	if name is None:
876	names = None
877	else:
878	names = []
879	stop = 0
880	currentname = name
881	while stop == 0:
882	indexvir = currentname.find(',')
883	if indexvir == -1:
884	stop = 1
885	name1 = currentname
886	else:
887	name1 = currentname[0:indexvir]
888	names = concatenate((names,[name1]))
889	currentname = currentname[indexvir+1:len(currentname)]
890	return names
891
892	## Author: TN [old]
893	def getopmatrix (kind,n):
894	import numpy as np
895	# compute matrix of operations between files
896	# the matrix is 'number of files'-square
897	# 1: difference (row minus column), 2: add
898	# not 0 in diag : just plot
899	if n == 1:
900	opm = np.eye(1)
901	elif kind == 'basic':
902	opm = np.eye(n)
903	elif kind == 'difference1': # show differences with 1st file
904	opm = np.zeros((n,n))
905	opm[0,:] = 1
906	opm[0,0] = 0
907	elif kind == 'difference2': # show differences with 1st file AND show 1st file
908	opm = np.zeros((n,n))
909	opm[0,:] = 1
910	else:
911	opm = np.eye(n)
912	return opm
913
914	## Author: TN [old]
915	def checkcoherence (nfiles,nlat,nlon,ntime):
916	if (nfiles > 1):
917	if (nlat > 1):
918	errormess("what you asked is not possible !")
919	return 1
920
921	## Author: TN
922	def readslices(saxis):
923	from numpy import empty
924	if saxis == None:
925	zesaxis = None
926	else:
927	zesaxis = empty((len(saxis),2))
928	for i in range(len(saxis)):
929	a = separatenames(saxis[i])
930	if len(a) == 1:
931	zesaxis[i,:] = float(a[0])
932	else:
933	zesaxis[i,0] = float(a[0])
934	zesaxis[i,1] = float(a[1])
935
936	return zesaxis
937
938	## Author: AS
939	def bidimfind(lon2d,lat2d,vlon,vlat):
940	import numpy as np
941	if vlat is None: array = (lon2d - vlon)**2
942	elif vlon is None: array = (lat2d - vlat)**2
943	else: array = (lon2d - vlon)2 + (lat2d - vlat)2
944	idy,idx = np.unravel_index( np.argmin(array), lon2d.shape )
945	if vlon is not None:
946	#print lon2d[idy,idx],vlon
947	if (np.abs(lon2d[idy,idx]-vlon)) > 5: errormess("longitude not found ",printvar=lon2d)
948	if vlat is not None:
949	#print lat2d[idy,idx],vlat
950	if (np.abs(lat2d[idy,idx]-vlat)) > 5: errormess("latitude not found ",printvar=lat2d)
951	return idx,idy
952
953	## Author: TN
954	def getsindex(saxis,index,axis):
955	# input : all the desired slices and the good index
956	# output : all indexes to be taken into account for reducing field
957	import numpy as np
958	### added by AS to treat the case of stime = - LT
959	if saxis is not None:
960	if saxis[0][0] < 0: saxis = - saxis
961	###
962	if ( np.array(axis).ndim == 2):
963	axis = axis[:,0]
964	if saxis is None:
965	zeindex = None
966	else:
967	aaa = int(np.argmin(abs(saxis[index,0] - axis)))
968	bbb = int(np.argmin(abs(saxis[index,1] - axis)))
969	[imin,imax] = np.sort(np.array([aaa,bbb]))
970	zeindex = np.array(range(imax-imin+1))+imin
971	# because -180 and 180 are the same point in longitude,
972	# we get rid of one for averaging purposes.
973	if axis[imin] == -180 and axis[imax] == 180:
974	zeindex = zeindex[0:len(zeindex)-1]
975	print "INFO: whole longitude averaging asked, so last point is not taken into account."
976	return zeindex
977
978	## Author: TN
979	def define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,indextime,what_I_plot,dim0,vertmode):
980	# Purpose of define_axis is to find x and y axis scales in a smart way
981	# x axis priority: 1/time 2/lon 3/lat 4/vertical
982	# To be improved !!!...
983	from numpy import array,swapaxes
984	x = None
985	y = None
986	count = 0
987	what_I_plot = array(what_I_plot)
988	shape = what_I_plot.shape
989	if indextime is None:
990	print "AXIS is time"
991	x = time
992	count = count+1
993	if indexlon is None:
994	print "AXIS is lon"
995	if count == 0: x = lon
996	else: y = lon
997	count = count+1
998	if indexlat is None:
999	print "AXIS is lat"
1000	if count == 0: x = lat
1001	else: y = lat
1002	count = count+1
1003	if indexvert is None and dim0 is 4:
1004	print "AXIS is vert"
1005	if vertmode == 0: # vertical axis is as is (GCM grid)
1006	if count == 0: x=range(len(vert))
1007	else: y=range(len(vert))
1008	count = count+1
1009	else: # vertical axis is in kms
1010	if count == 0: x = vert
1011	else: y = vert
1012	count = count+1
1013	x = array(x)
1014	y = array(y)
1015	print "CHECK: what_I_plot.shape", what_I_plot.shape
1016	print "CHECK: x.shape, y.shape", x.shape, y.shape
1017	if len(shape) == 1:
1018	if shape[0] != len(x):
1019	print "WARNING HERE !!!"
1020	x = y
1021	elif len(shape) == 2:
1022	if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]:
1023	what_I_plot = swapaxes(what_I_plot,0,1)
1024	print "INFO: swapaxes", what_I_plot.shape, shape
1025	return what_I_plot,x,y
1026
1027	# Author: TN + AS
1028	def determineplot(slon, slat, svert, stime):
1029	nlon = 1 # number of longitudinal slices -- 1 is None
1030	nlat = 1
1031	nvert = 1
1032	ntime = 1
1033	nslices = 1
1034	if slon is not None:
1035	nslices = nslices*len(slon)
1036	nlon = len(slon)
1037	if slat is not None:
1038	nslices = nslices*len(slat)
1039	nlat = len(slat)
1040	if svert is not None:
1041	nslices = nslices*len(svert)
1042	nvert = len(svert)
1043	if stime is not None:
1044	nslices = nslices*len(stime)
1045	ntime = len(stime)
1046	#else:
1047	# nslices = 2
1048	mapmode = 0
1049	if slon is None and slat is None:
1050	mapmode = 1 # in this case we plot a map, with the given projection
1051
1052	return nlon, nlat, nvert, ntime, mapmode, nslices
1053
1054	## Author: AC
1055	## Reduce complexity of main script by moving the contour part here. Also allow to call it from elsewhere
1056	## which can be usefull
1057	#
1058	#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):
1059	# from matplotlib.pyplot import contour, plot, clabel
1060	# import numpy as np
1061	# #what_I_plot = what_I_plot*mult
1062	# if not error:
1063	# if mapmode == 1:
1064	# if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_plot,6)
1065	# zevmin, zevmax = calculate_bounds(what_I_plot)
1066	# zelevels = np.linspace(zevmin,zevmax,ticks) #20)
1067	# if var2 == 'HGT': zelevels = np.arange(-10000.,30000.,2000.)
1068	# if mapmode == 0:
1069	# #if typefile in ['mesoideal']: what_I_plot = dumpbdy(what_I_plot,0,stag='W')
1070	# itime=indextime
1071	# if len(what_I_plot.shape) is 3: itime=[0]
1072	# what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
1073	# itime,what_I_plot, len(ze_var2.shape),vertmode)
1074	# ### If we plot a 2-D field
1075	# if len(what_I_plot.shape) is 2:
1076	# #zelevels=[1.]
1077	# if mapmode == 0:cs = contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5)
1078	# elif mapmode == 1:cs = m.contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5)
1079	# #clabel(cs,zelevels,inline=3,fmt='%1.1f',fontsize=7)
1080	# ### If we plot a 1-D field
1081	# elif len(what_I_plot.shape) is 1:
1082	# plot(what_I_plot,x)
1083	# else:
1084	# errormess("There is an error in reducing field !")
1085	# return error
1086

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

Download in other formats: