Context Navigation

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

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

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

Download in other formats: