Context Navigation

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

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

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

Download in other formats: