Context Navigation

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

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

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

Download in other formats: