Context Navigation

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

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

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

Download in other formats: