Context Navigation

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

Last change on this file since 579 was 579, checked in by aslmd, 13 years ago
PYTHON UTIL. improved comments for the first part. corrected a bug: for 3D fields xyt it was necessary to add --vert 0 to have x,y,or t as abscissae, this is no longer the case.
File size: 50.5 KB

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

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

Download in other formats: