Context Navigation

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

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

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

Download in other formats: