Context Navigation

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

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

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

Download in other formats: