Context Navigation

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

Last change on this file since 523 was 518, checked in by aslmd, 13 years ago
UTIL/PYTHON: added the possibility to use -c onebar for multiplots as one would easily guess why
File size: 46.1 KB

Line
1	## Author: AS
2	def errormess(text,printvar=None):
3	print text
4	if printvar is not None: print printvar
5	exit()
6	return
7
8	## Author: AS
9	def adjust_length (tab, zelen):
10	from numpy import ones
11	if tab is None:
12	outtab = ones(zelen) * -999999
13	else:
14	if zelen != len(tab):
15	print "not enough or too much values... setting same values all variables"
16	outtab = ones(zelen) * tab[0]
17	else:
18	outtab = tab
19	return outtab
20
21	## Author: AS
22	def getname(var=False,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	elif dimension == 1: field = nc.variables[var][:]
62	# if there are NaNs in the ncdf, they should be loaded as a masked array which will be
63	# recasted as a regular array later in reducefield
64	if (np.isnan(np.sum(field)) and (type(field).__name__ not in 'MaskedArray')):
65	print "Warning: netcdf as nan values but is not loaded as a Masked Array."
66	print "recasting array type"
67	out=np.ma.masked_invalid(field)
68	out.set_fill_value([np.NaN])
69	else:
70	# missing values from zrecast or hrecast are -1e-33
71	masked=np.ma.masked_where(field < -1e30,field)
72	masked.set_fill_value([np.NaN])
73	mask = np.ma.getmask(masked)
74	if (True in np.array(mask)):out=masked
75	else:out=field
76	return out
77
78	## Author: AS + TN + AC
79	def reducefield (input,d4=None,d3=None,d2=None,d1=None,yint=False,alt=None,anomaly=False,redope=None):
80	### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x"
81	### it would be actually better to name d4 d3 d2 d1 as t z y x
82	### ... note, anomaly is only computed over d1 and d2 for the moment
83	import numpy as np
84	from mymath import max,mean,min
85	csmooth = 12 ## a fair amount of grid points (too high results in high computation time)
86	if redope is not None:
87	if redope == "mint": input = min(input,axis=0) ; d1 = None
88	elif redope == "maxt": input = max(input,axis=0) ; d1 = None
89	else: errormess("not supported. but try lines in reducefield beforehand.")
90	#elif redope == "minz": input = min(input,axis=1) ; d2 = None
91	#elif redope == "maxz": input = max(input,axis=1) ; d2 = None
92	#elif redope == "miny": input = min(input,axis=2) ; d3 = None
93	#elif redope == "maxy": input = max(input,axis=2) ; d3 = None
94	#elif redope == "minx": input = min(input,axis=3) ; d4 = None
95	#elif redope == "maxx": input = max(input,axis=3) ; d4 = None
96	dimension = np.array(input).ndim
97	shape = np.array(input).shape
98	#print 'd1,d2,d3,d4: ',d1,d2,d3,d4
99	if anomaly: print 'ANOMALY ANOMALY'
100	output = input
101	error = False
102	#### this is needed to cope the case where d4,d3,d2,d1 are single integers and not arrays
103	if d4 is not None and not isinstance(d4, np.ndarray): d4=[d4]
104	if d3 is not None and not isinstance(d3, np.ndarray): d3=[d3]
105	if d2 is not None and not isinstance(d2, np.ndarray): d2=[d2]
106	if d1 is not None and not isinstance(d1, np.ndarray): d1=[d1]
107	### now the main part
108	if dimension == 2:
109	if d2 >= shape[0]: error = True
110	elif d1 >= shape[1]: error = True
111	elif d1 is not None and d2 is not None: output = mean(input[d2,:],axis=0); output = mean(output[d1],axis=0)
112	elif d1 is not None: output = mean(input[:,d1],axis=1)
113	elif d2 is not None: output = mean(input[d2,:],axis=0)
114	elif dimension == 3:
115	if max(d4) >= shape[0]: error = True
116	elif max(d2) >= shape[1]: error = True
117	elif max(d1) >= shape[2]: error = True
118	elif d4 is not None and d2 is not None and d1 is not None:
119	output = mean(input[d4,:,:],axis=0); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0)
120	elif d4 is not None and d2 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[d2,:],axis=0)
121	elif d4 is not None and d1 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1)
122	elif d2 is not None and d1 is not None: output = mean(input[:,d2,:],axis=1); output=mean(output[:,d1],axis=1)
123	elif d1 is not None: output = mean(input[:,:,d1],axis=2)
124	elif d2 is not None: output = mean(input[:,d2,:],axis=1)
125	elif d4 is not None: output = mean(input[d4,:,:],axis=0)
126	elif dimension == 4:
127	if max(d4) >= shape[0]: error = True
128	elif max(d3) >= shape[1]: error = True
129	elif max(d2) >= shape[2]: error = True
130	elif max(d1) >= shape[3]: error = True
131	elif d4 is not None and d3 is not None and d2 is not None and d1 is not None:
132	output = mean(input[d4,:,:,:],axis=0)
133	output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
134	if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
135	output = mean(output[d2,:],axis=0)
136	output = mean(output[d1],axis=0)
137	elif d4 is not None and d3 is not None and d2 is not None:
138	output = mean(input[d4,:,:,:],axis=0)
139	output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
140	if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
141	output = mean(output[d2,:],axis=0)
142	elif d4 is not None and d3 is not None and d1 is not None:
143	output = mean(input[d4,:,:,:],axis=0)
144	output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
145	if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
146	output = mean(output[:,d1],axis=1)
147	elif d4 is not None and d2 is not None and d1 is not None:
148	output = mean(input[d4,:,:,:],axis=0)
149	if anomaly:
150	for k in range(output.shape[0]): output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
151	output = mean(output[:,d2,:],axis=1)
152	output = mean(output[:,d1],axis=1)
153	#noperturb = smooth1d(output,window_len=7)
154	#lenlen = len(output) ; output = output[1:lenlen-7] ; yeye = noperturb[4:lenlen-4]
155	#plot(output) ; plot(yeye) ; show() ; plot(output-yeye) ; show()
156	elif d3 is not None and d2 is not None and d1 is not None:
157	output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
158	if anomaly:
159	for k in range(output.shape[0]): output[k,:,:] = 100. * ((output[k,:,:] / smooth(output[k,:,:],csmooth)) - 1.)
160	output = mean(output[:,d2,:],axis=1)
161	output = mean(output[:,d1],axis=1)
162	elif d4 is not None and d3 is not None:
163	output = mean(input[d4,:,:,:],axis=0)
164	output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
165	if anomaly: output = 100. * ((output / smooth(output,csmooth)) - 1.)
166	elif d4 is not None and d2 is not None:
167	output = mean(input[d4,:,:,:],axis=0)
168	output = mean(output[:,d2,:],axis=1)
169	elif d4 is not None and d1 is not None:
170	output = mean(input[d4,:,:,:],axis=0)
171	output = mean(output[:,:,d1],axis=2)
172	elif d3 is not None and d2 is not None:
173	output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
174	output = mean(output[:,d2,:],axis=1)
175	elif d3 is not None and d1 is not None:
176	output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
177	output = mean(output[:,:,d1],axis=2)
178	elif d2 is not None and d1 is not None:
179	output = mean(input[:,:,d2,:],axis=2)
180	output = mean(output[:,:,d1],axis=2)
181	elif d1 is not None: output = mean(input[:,:,:,d1],axis=3)
182	elif d2 is not None: output = mean(input[:,:,d2,:],axis=2)
183	elif d3 is not None: output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
184	elif d4 is not None: output = mean(input[d4,:,:,:],axis=0)
185	dimension2 = np.array(output).ndim
186	shape2 = np.array(output).shape
187	print 'REDUCEFIELD dim,shape: ',dimension,shape,' >>> ',dimension2,shape2
188	return output, error
189
190	## Author: AC + AS
191	def reduce_zaxis (input,ax=None,yint=False,vert=None,indice=None):
192	from mymath import max,mean
193	from scipy import integrate
194	if yint and vert is not None and indice is not None:
195	if type(input).__name__=='MaskedArray':
196	input.set_fill_value([np.NaN])
197	output = integrate.trapz(input.filled(),x=vert[indice],axis=ax)
198	else:
199	output = integrate.trapz(input,x=vert[indice],axis=ax)
200	else:
201	output = mean(input,axis=ax)
202	return output
203
204	## Author: AS + TN
205	def definesubplot ( numplot, fig ):
206	from matplotlib.pyplot import rcParams
207	rcParams['font.size'] = 12. ## default (important for multiple calls)
208	if numplot <= 0:
209	subv = 99999
210	subh = 99999
211	elif numplot == 1:
212	subv = 99999
213	subh = 99999
214	elif numplot == 2:
215	subv = 1 #2
216	subh = 2 #1
217	fig.subplots_adjust(wspace = 0.35)
218	rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
219	elif numplot == 3:
220	subv = 3
221	subh = 1
222	fig.subplots_adjust(wspace = 0.5)
223	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
224	elif numplot == 4:
225	subv = 2
226	subh = 2
227	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
228	rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
229	elif numplot <= 6:
230	subv = 2
231	subh = 3
232	fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
233	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
234	elif numplot <= 8:
235	subv = 2
236	subh = 4
237	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
238	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
239	elif numplot <= 9:
240	subv = 3
241	subh = 3
242	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
243	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
244	elif numplot <= 12:
245	subv = 3
246	subh = 4
247	fig.subplots_adjust(wspace = 0, hspace = 0.1)
248	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
249	elif numplot <= 16:
250	subv = 4
251	subh = 4
252	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
253	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
254	else:
255	print "number of plot supported: 1 to 16"
256	exit()
257	return subv,subh
258
259	## Author: AS
260	def getstralt(nc,nvert):
261	typefile = whatkindfile(nc)
262	if typefile is 'meso':
263	stralt = "_lvl" + str(nvert)
264	elif typefile is 'mesoapi':
265	zelevel = int(nc.variables['vert'][nvert])
266	if abs(zelevel) < 10000.: strheight=str(zelevel)+"m"
267	else: strheight=str(int(zelevel/1000.))+"km"
268	if 'altitude' in nc.dimensions: stralt = "_"+strheight+"-AMR"
269	elif 'altitude_abg' in nc.dimensions: stralt = "_"+strheight+"-ALS"
270	elif 'bottom_top' in nc.dimensions: stralt = "_"+strheight
271	elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa"
272	else: stralt = ""
273	else:
274	stralt = ""
275	return stralt
276
277	## Author: AS
278	def getlschar ( namefile, getaxis=False ):
279	from netCDF4 import Dataset
280	from timestuff import sol2ls
281	from numpy import array
282	from string import rstrip
283	namefile = rstrip( rstrip( rstrip( namefile, chars="_z"), chars="_zabg"), chars="_p")
284	#### we assume that wrfout is next to wrfout_z and wrfout_zabg
285	nc = Dataset(namefile)
286	zetime = None
287	days_in_month = [61, 66, 66, 65, 60, 54, 50, 46, 47, 47, 51, 56]
288	plus_in_month = [ 0, 61,127,193,258,318,372,422,468,515,562,613]
289	if 'Times' in nc.variables:
290	zetime = nc.variables['Times'][0]
291	shape = array(nc.variables['Times']).shape
292	if shape[0] < 2: zetime = None
293	if zetime is not None \
294	and 'vert' not in nc.variables:
295	##### strangely enough this does not work for api or ncrcat results!
296	zesol = plus_in_month[ int(zetime[5]+zetime[6])-1 ] + int(zetime[8]+zetime[9]) - 1 ##les sols GCM commencent a 0
297	dals = int( 10. * sol2ls ( zesol ) ) / 10.
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	if not getaxis:
305	lschar = "_Ls"+str(dals)
306	else:
307	zelen = len(nc.variables['Times'][:])
308	yeye = range(zelen) ; lsaxis = range(zelen) ; solaxis = range(zelen) ; ltaxis = range(zelen)
309	for iii in yeye:
310	zetime = nc.variables['Times'][iii]
311	ltaxis[iii] = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37.
312	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
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], getattr( nc, 'JULDAY' )
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	if char not in ["moll"]:
626	m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, latmax=zelatmax)
627	m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color=zecolor, linewidth=zelinewidth, fontsize=fontsizemer, latmax=zelatmax)
628	if back: m.warpimage(marsmap(back),scale=0.75)
629	#if not back:
630	# if not var: back = "mola" ## if no var: draw mola
631	# elif typefile in ['mesoapi','meso','geo'] \
632	# and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw
633	# else: pass ## else: draw None
634	return m
635
636	## Author: AS
637	#### test temporaire
638	def putpoints (map,plot):
639	#### from http://www.scipy.org/Cookbook/Matplotlib/Maps
640	# lat/lon coordinates of five cities.
641	lats = [18.4]
642	lons = [-134.0]
643	points=['Olympus Mons']
644	# compute the native map projection coordinates for cities.
645	x,y = map(lons,lats)
646	# plot filled circles at the locations of the cities.
647	map.plot(x,y,'bo')
648	# plot the names of those five cities.
649	wherept = 0 #1000 #50000
650	for name,xpt,ypt in zip(points,x,y):
651	plot.text(xpt+wherept,ypt+wherept,name)
652	## le nom ne s'affiche pas...
653	return
654
655	## Author: AS
656	def calculate_bounds(field,vmin=None,vmax=None):
657	import numpy as np
658	from mymath import max,min,mean
659	ind = np.where(field < 9e+35)
660	fieldcalc = field[ ind ] # la syntaxe compacte ne marche si field est un tuple
661	###
662	dev = np.std(fieldcalc)*3.0
663	###
664	if vmin is None:
665	zevmin = mean(fieldcalc) - dev
666	else: zevmin = vmin
667	###
668	if vmax is None: zevmax = mean(fieldcalc) + dev
669	else: zevmax = vmax
670	if vmin == vmax:
671	zevmin = mean(fieldcalc) - dev ### for continuity
672	zevmax = mean(fieldcalc) + dev ### for continuity
673	###
674	if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0.
675	print "BOUNDS field ", min(fieldcalc), max(fieldcalc), " //// adopted", zevmin, zevmax
676	return zevmin, zevmax
677
678	## Author: AS
679	def bounds(what_I_plot,zevmin,zevmax):
680	from mymath import max,min,mean
681	### might be convenient to add the missing value in arguments
682	#what_I_plot[ what_I_plot < zevmin ] = zevmin#*(1. + 1.e-7)
683	if zevmin < 0: what_I_plot[ what_I_plot < zevmin(1. - 1.e-7) ] = zevmin(1. - 1.e-7)
684	else: what_I_plot[ what_I_plot < zevmin(1. + 1.e-7) ] = zevmin(1. + 1.e-7)
685	#print "NEW MIN ", min(what_I_plot)
686	what_I_plot[ what_I_plot > 9e+35 ] = -9e+35
687	what_I_plot[ what_I_plot > zevmax ] = zevmax
688	#print "NEW MAX ", max(what_I_plot)
689	return what_I_plot
690
691	## Author: AS
692	def nolow(what_I_plot):
693	from mymath import max,min
694	lim = 0.150.5(abs(max(what_I_plot))+abs(min(what_I_plot)))
695	print "NO PLOT BELOW VALUE ", lim
696	what_I_plot [ abs(what_I_plot) < lim ] = 1.e40
697	return what_I_plot
698
699
700	## Author : AC
701	def hole_bounds(what_I_plot,zevmin,zevmax):
702	import numpy as np
703	zi=0
704	for i in what_I_plot:
705	zj=0
706	for j in i:
707	if ((j < zevmin) or (j > zevmax)):what_I_plot[zi,zj]=np.NaN
708	zj=zj+1
709	zi=zi+1
710
711	return what_I_plot
712
713	## Author: AS
714	def zoomset (wlon,wlat,zoom):
715	dlon = abs(wlon[1]-wlon[0])/2.
716	dlat = abs(wlat[1]-wlat[0])/2.
717	[wlon,wlat] = [ [wlon[0]+zoomdlon/100.,wlon[1]-zoomdlon/100.],\
718	[wlat[0]+zoomdlat/100.,wlat[1]-zoomdlat/100.] ]
719	print "ZOOM %",zoom,wlon,wlat
720	return wlon,wlat
721
722	## Author: AS
723	def fmtvar (whichvar="def"):
724	fmtvar = { \
725	"MIXED": "%.0f",\
726	"UPDRAFT": "%.0f",\
727	"DOWNDRAFT": "%.0f",\
728	"TK": "%.0f",\
729	#"ZMAX_TH": "%.0f",\
730	#"WSTAR": "%.0f",\
731	# Variables from TES ncdf format
732	"T_NADIR_DAY": "%.0f",\
733	"T_NADIR_NIT": "%.0f",\
734	# Variables from tes.py ncdf format
735	"TEMP_DAY": "%.0f",\
736	"TEMP_NIGHT": "%.0f",\
737	# Variables from MCS and mcs.py ncdf format
738	"DTEMP": "%.0f",\
739	"NTEMP": "%.0f",\
740	"DNUMBINTEMP": "%.0f",\
741	"NNUMBINTEMP": "%.0f",\
742	# other stuff
743	"TPOT": "%.0f",\
744	"TSURF": "%.0f",\
745	"def": "%.1e",\
746	"PTOT": "%.0f",\
747	"HGT": "%.1e",\
748	"USTM": "%.2f",\
749	"HFX": "%.0f",\
750	"ICETOT": "%.1e",\
751	"TAU_ICE": "%.2f",\
752	"TAUICE": "%.2f",\
753	"VMR_ICE": "%.1e",\
754	"MTOT": "%.1f",\
755	"ANOMALY": "%.1f",\
756	"W": "%.1f",\
757	"WMAX_TH": "%.1f",\
758	"QSURFICE": "%.0f",\
759	"UM": "%.0f",\
760	"WIND": "%.0f",\
761	"ALBBARE": "%.2f",\
762	"TAU": "%.1f",\
763	"CO2": "%.2f",\
764	#### T.N.
765	"TEMP": "%.0f",\
766	"VMR_H2OICE": "%.0f",\
767	"VMR_H2OVAP": "%.0f",\
768	"TAUTES": "%.2f",\
769	"TAUTESAP": "%.2f",\
770
771	}
772	if "TSURF" in whichvar: whichvar = "TSURF"
773	if whichvar not in fmtvar:
774	whichvar = "def"
775	return fmtvar[whichvar]
776
777	## Author: AS
778	####################################################################################################################
779	### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
780	def defcolorb (whichone="def"):
781	whichcolorb = { \
782	"def": "spectral",\
783	"HGT": "spectral",\
784	"HGT_M": "spectral",\
785	"TK": "gist_heat",\
786	"TPOT": "Paired",\
787	"TSURF": "RdBu_r",\
788	"QH2O": "PuBu",\
789	"USTM": "YlOrRd",\
790	"WIND": "YlOrRd",\
791	#"T_nadir_nit": "RdBu_r",\
792	#"T_nadir_day": "RdBu_r",\
793	"HFX": "RdYlBu",\
794	"ICETOT": "YlGnBu_r",\
795	#"MTOT": "PuBu",\
796	"CCNQ": "YlOrBr",\
797	"CCNN": "YlOrBr",\
798	"TEMP": "Jet",\
799	"TAU_ICE": "Blues",\
800	"TAUICE": "Blues",\
801	"VMR_ICE": "Blues",\
802	"W": "jet",\
803	"WMAX_TH": "spectral",\
804	"ANOMALY": "RdBu_r",\
805	"QSURFICE": "hot_r",\
806	"ALBBARE": "spectral",\
807	"TAU": "YlOrBr_r",\
808	"CO2": "YlOrBr_r",\
809	#### T.N.
810	"MTOT": "Jet",\
811	"H2O_ICE_S": "RdBu",\
812	"VMR_H2OICE": "PuBu",\
813	"VMR_H2OVAP": "PuBu",\
814	"WATERCAPTAG": "Blues",\
815	}
816	#W --> spectral ou jet
817	#spectral BrBG RdBu_r
818	#print "predefined colorbars"
819	if "TSURF" in whichone: whichone = "TSURF"
820	if whichone not in whichcolorb:
821	whichone = "def"
822	return whichcolorb[whichone]
823
824	## Author: AS
825	def definecolorvec (whichone="def"):
826	whichcolor = { \
827	"def": "black",\
828	"vis": "yellow",\
829	"vishires": "yellow",\
830	"molabw": "yellow",\
831	"mola": "black",\
832	"gist_heat": "white",\
833	"hot": "tk",\
834	"gist_rainbow": "black",\
835	"spectral": "black",\
836	"gray": "red",\
837	"PuBu": "black",\
838	}
839	if whichone not in whichcolor:
840	whichone = "def"
841	return whichcolor[whichone]
842
843	## Author: AS
844	def marsmap (whichone="vishires"):
845	from os import uname
846	mymachine = uname()[1]
847	### not sure about speed-up with this method... looks the same
848	if "lmd.jussieu.fr" in mymachine: domain = "/u/aslmd/WWW/maps/"
849	elif "aymeric-laptop" in mymachine: domain = "/home/aymeric/Dropbox/Public/"
850	else: domain = "http://www.lmd.jussieu.fr/~aslmd/maps/"
851	whichlink = { \
852	#"vis": "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\
853	#"vishires": "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\
854	#"geolocal": "http://dl.dropbox.com/u/11078310/geolocal.jpg",\
855	#"mola": "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\
856	#"molabw": "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\
857	"thermalday": domain+"thermalday.jpg",\
858	"thermalnight": domain+"thermalnight.jpg",\
859	"tesalbedo": domain+"tesalbedo.jpg",\
860	"vis": domain+"mar0kuu2.jpg",\
861	"vishires": domain+"MarsMap_2500x1250.jpg",\
862	"geolocal": domain+"geolocal.jpg",\
863	"mola": domain+"mars-mola-2k.jpg",\
864	"molabw": domain+"MarsElevation_2500x1250.jpg",\
865	"clouds": "http://www.johnstonsarchive.net/spaceart/marswcloudmap.jpg",\
866	"jupiter": "http://www.mmedia.is/~bjj/data/jupiter_css/jupiter_css.jpg",\
867	"jupiter_voy": "http://www.mmedia.is/~bjj/data/jupiter/jupiter_vgr2.jpg",\
868	"bw": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg",\
869	"contrast": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg",\
870	"nice": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg",\
871	"blue": "http://eoimages.gsfc.nasa.gov/ve/2430/land_ocean_ice_2048.jpg",\
872	"blueclouds": "http://eoimages.gsfc.nasa.gov/ve/2431/land_ocean_ice_cloud_2048.jpg",\
873	"justclouds": "http://eoimages.gsfc.nasa.gov/ve/2432/cloud_combined_2048.jpg",\
874	}
875	### see http://www.mmedia.is/~bjj/planetary_maps.html
876	if whichone not in whichlink:
877	print "marsmap: choice not defined... you'll get the default one... "
878	whichone = "vishires"
879	return whichlink[whichone]
880
881	#def earthmap (whichone):
882	# if whichone == "contrast": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg"
883	# elif whichone == "bw": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg"
884	# elif whichone == "nice": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg"
885	# return whichlink
886
887	## Author: AS
888	def latinterv (area="Whole"):
889	list = { \
890	"Europe": [[ 20., 80.],[- 50., 50.]],\
891	"Central_America": [[-10., 40.],[ 230., 300.]],\
892	"Africa": [[-20., 50.],[- 50., 50.]],\
893	"Whole": [[-90., 90.],[-180., 180.]],\
894	"Southern_Hemisphere": [[-90., 60.],[-180., 180.]],\
895	"Northern_Hemisphere": [[-60., 90.],[-180., 180.]],\
896	"Tharsis": [[-30., 60.],[-170.,- 10.]],\
897	"Whole_No_High": [[-60., 60.],[-180., 180.]],\
898	"Chryse": [[-60., 60.],[- 60., 60.]],\
899	"North_Pole": [[ 50., 90.],[-180., 180.]],\
900	"Close_North_Pole": [[ 75., 90.],[-180., 180.]],\
901	"Far_South_Pole": [[-90.,-40.],[-180., 180.]],\
902	"South_Pole": [[-90.,-50.],[-180., 180.]],\
903	"Close_South_Pole": [[-90.,-75.],[-180., 180.]],\
904	"Sirenum_Crater_large": [[-46.,-34.],[-166., -151.]],\
905	"Sirenum_Crater_small": [[-36.,-26.],[-168., -156.]],\
906
907	}
908	if area not in list: area = "Whole"
909	[olat,olon] = list[area]
910	return olon,olat
911
912	## Author: TN
913	def separatenames (name):
914	from numpy import concatenate
915	# look for comas in the input name to separate different names (files, variables,etc ..)
916	if name is None:
917	names = None
918	else:
919	names = []
920	stop = 0
921	currentname = name
922	while stop == 0:
923	indexvir = currentname.find(',')
924	if indexvir == -1:
925	stop = 1
926	name1 = currentname
927	else:
928	name1 = currentname[0:indexvir]
929	names = concatenate((names,[name1]))
930	currentname = currentname[indexvir+1:len(currentname)]
931	return names
932
933	## Author: TN [old]
934	def getopmatrix (kind,n):
935	import numpy as np
936	# compute matrix of operations between files
937	# the matrix is 'number of files'-square
938	# 1: difference (row minus column), 2: add
939	# not 0 in diag : just plot
940	if n == 1:
941	opm = np.eye(1)
942	elif kind == 'basic':
943	opm = np.eye(n)
944	elif kind == 'difference1': # show differences with 1st file
945	opm = np.zeros((n,n))
946	opm[0,:] = 1
947	opm[0,0] = 0
948	elif kind == 'difference2': # show differences with 1st file AND show 1st file
949	opm = np.zeros((n,n))
950	opm[0,:] = 1
951	else:
952	opm = np.eye(n)
953	return opm
954
955	## Author: TN [old]
956	def checkcoherence (nfiles,nlat,nlon,ntime):
957	if (nfiles > 1):
958	if (nlat > 1):
959	errormess("what you asked is not possible !")
960	return 1
961
962	## Author: TN
963	def readslices(saxis):
964	from numpy import empty
965	if saxis == None:
966	zesaxis = None
967	else:
968	zesaxis = empty((len(saxis),2))
969	for i in range(len(saxis)):
970	a = separatenames(saxis[i])
971	if len(a) == 1:
972	zesaxis[i,:] = float(a[0])
973	else:
974	zesaxis[i,0] = float(a[0])
975	zesaxis[i,1] = float(a[1])
976
977	return zesaxis
978
979	## Author: AS
980	def bidimfind(lon2d,lat2d,vlon,vlat,file=None):
981	import numpy as np
982	import matplotlib.pyplot as mpl
983	if vlat is None: array = (lon2d - vlon)**2
984	elif vlon is None: array = (lat2d - vlat)**2
985	else: array = (lon2d - vlon)2 + (lat2d - vlat)2
986	idy,idx = np.unravel_index( np.argmin(array), lon2d.shape )
987	if vlon is not None:
988	if (np.abs(lon2d[idy,idx]-vlon)) > 5: errormess("longitude not found ",printvar=lon2d)
989	if vlat is not None:
990	if (np.abs(lat2d[idy,idx]-vlat)) > 5: errormess("latitude not found ",printvar=lat2d)
991	if file is not None:
992	print idx,idy,lon2d[idy,idx],vlon
993	print idx,idy,lat2d[idy,idx],vlat
994	var = file.variables["HGT"][:,:,:]
995	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)
996	mpl.show()
997	return idy,idx
998
999	## Author: TN
1000	def getsindex(saxis,index,axis):
1001	# input : all the desired slices and the good index
1002	# output : all indexes to be taken into account for reducing field
1003	import numpy as np
1004	if ( np.array(axis).ndim == 2):
1005	axis = axis[:,0]
1006	if saxis is None:
1007	zeindex = None
1008	else:
1009	aaa = int(np.argmin(abs(saxis[index,0] - axis)))
1010	bbb = int(np.argmin(abs(saxis[index,1] - axis)))
1011	[imin,imax] = np.sort(np.array([aaa,bbb]))
1012	zeindex = np.array(range(imax-imin+1))+imin
1013	# because -180 and 180 are the same point in longitude,
1014	# we get rid of one for averaging purposes.
1015	if axis[imin] == -180 and axis[imax] == 180:
1016	zeindex = zeindex[0:len(zeindex)-1]
1017	print "INFO: whole longitude averaging asked, so last point is not taken into account."
1018	return zeindex
1019
1020	## Author: TN
1021	def define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,indextime,what_I_plot,dim0,vertmode):
1022	# Purpose of define_axis is to find x and y axis scales in a smart way
1023	# x axis priority: 1/time 2/lon 3/lat 4/vertical
1024	# To be improved !!!...
1025	from numpy import array,swapaxes
1026	x = None
1027	y = None
1028	count = 0
1029	what_I_plot = array(what_I_plot)
1030	shape = what_I_plot.shape
1031	if indextime is None and len(time) > 1:
1032	print "AXIS is time"
1033	x = time
1034	count = count+1
1035	if indexlon is None and len(lon) > 1:
1036	print "AXIS is lon"
1037	if count == 0: x = lon
1038	else: y = lon
1039	count = count+1
1040	if indexlat is None and len(lon) > 1:
1041	print "AXIS is lat"
1042	if count == 0: x = lat
1043	else: y = lat
1044	count = count+1
1045	if indexvert is None and ((dim0 is 4) or (y is None)):
1046	print "AXIS is vert"
1047	if vertmode == 0: # vertical axis is as is (GCM grid)
1048	if count == 0: x=range(len(vert))
1049	else: y=range(len(vert))
1050	count = count+1
1051	else: # vertical axis is in kms
1052	if count == 0: x = vert
1053	else: y = vert
1054	count = count+1
1055	x = array(x)
1056	y = array(y)
1057	print "CHECK SHAPE: what_I_plot, x, y", what_I_plot.shape, x.shape, y.shape
1058	if len(shape) == 1:
1059	if shape[0] != len(x):
1060	print "WARNING HERE !!!"
1061	x = y
1062	elif len(shape) == 2:
1063	if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]:
1064	what_I_plot = swapaxes(what_I_plot,0,1)
1065	print "INFO: swapaxes", what_I_plot.shape, shape
1066	return what_I_plot,x,y
1067
1068	# Author: TN + AS
1069	def determineplot(slon, slat, svert, stime):
1070	nlon = 1 # number of longitudinal slices -- 1 is None
1071	nlat = 1
1072	nvert = 1
1073	ntime = 1
1074	nslices = 1
1075	if slon is not None:
1076	nslices = nslices*len(slon)
1077	nlon = len(slon)
1078	if slat is not None:
1079	nslices = nslices*len(slat)
1080	nlat = len(slat)
1081	if svert is not None:
1082	nslices = nslices*len(svert)
1083	nvert = len(svert)
1084	if stime is not None:
1085	nslices = nslices*len(stime)
1086	ntime = len(stime)
1087	#else:
1088	# nslices = 2
1089	mapmode = 0
1090	if slon is None and slat is None:
1091	mapmode = 1 # in this case we plot a map, with the given projection
1092
1093	return nlon, nlat, nvert, ntime, mapmode, nslices
1094
1095	## Author: AC
1096	## Reduce complexity of main script by moving the contour part here. Also allow to call it from elsewhere
1097	## which can be usefull
1098	#
1099	#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):
1100	# from matplotlib.pyplot import contour, plot, clabel
1101	# import numpy as np
1102	# #what_I_plot = what_I_plot*mult
1103	# if not error:
1104	# if mapmode == 1:
1105	# if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_plot,6)
1106	# zevmin, zevmax = calculate_bounds(what_I_plot)
1107	# zelevels = np.linspace(zevmin,zevmax,ticks) #20)
1108	# if var2 == 'HGT': zelevels = np.arange(-10000.,30000.,2000.)
1109	# if mapmode == 0:
1110	# #if typefile in ['mesoideal']: what_I_plot = dumpbdy(what_I_plot,0,stag='W')
1111	# itime=indextime
1112	# if len(what_I_plot.shape) is 3: itime=[0]
1113	# what_I_plot, x, y = define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,\
1114	# itime,what_I_plot, len(ze_var2.shape),vertmode)
1115	# ### If we plot a 2-D field
1116	# if len(what_I_plot.shape) is 2:
1117	# #zelevels=[1.]
1118	# if mapmode == 0:cs = contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5)
1119	# elif mapmode == 1:cs = m.contour(x,y,what_I_plot, zelevels, colors='k', linewidths = 1 ) #0.33 colors='w' )# , alpha=0.5)
1120	# #clabel(cs,zelevels,inline=3,fmt='%1.1f',fontsize=7)
1121	# ### If we plot a 1-D field
1122	# elif len(what_I_plot.shape) is 1:
1123	# plot(what_I_plot,x)
1124	# else:
1125	# errormess("There is an error in reducing field !")
1126	# return error
1127

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

Download in other formats: