Context Navigation

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

Last change on this file since 392 was 392, checked in by aslmd, 13 years ago
PYTHON: quite a lot of modifs and tests. extended multivar subplot capabilities + cosmetic changes + general cleaning. corrected some stuff that were not working with mesoscale. and there is now a common script for meso and gcm svn status -uq! it is named pp.py [stands for: planetoplot]. put some examples in README file. testing is probably a bit necessary for most complex options.
File size: 35.7 KB

Line
1	## Author: AS
2	def errormess(text,printvar=None):
3	print text
4	if printvar: 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,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	return basename
29
30	## Author: AS
31	def localtime(utc,lon):
32	ltst = utc + lon / 15.
33	ltst = int (ltst * 10) / 10.
34	ltst = ltst % 24
35	return ltst
36
37	## Author: AS
38	def whatkindfile (nc):
39	if 'controle' in nc.variables: typefile = 'gcm'
40	elif 'phisinit' in nc.variables: typefile = 'gcm'
41	elif 'vert' in nc.variables: typefile = 'mesoapi'
42	elif 'U' in nc.variables: typefile = 'meso'
43	elif 'HGT_M' in nc.variables: typefile = 'geo'
44	#else: errormess("whatkindfile: typefile not supported.")
45	else: typefile = 'gcm' # for lslin-ed files from gcm
46	return typefile
47
48	## Author: AS
49	def getfield (nc,var):
50	## this allows to get much faster (than simply referring to nc.variables[var])
51	dimension = len(nc.variables[var].dimensions)
52	#print " Opening variable with", dimension, "dimensions ..."
53	if dimension == 2: field = nc.variables[var][:,:]
54	elif dimension == 3: field = nc.variables[var][:,:,:]
55	elif dimension == 4: field = nc.variables[var][:,:,:,:]
56	return field
57
58	## Author: AS + TN + AC
59	def reducefield (input,d4=None,d3=None,d2=None,d1=None,yint=False,alt=None):
60	### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x"
61	### it would be actually better to name d4 d3 d2 d1 as t z y x
62	import numpy as np
63	from mymath import max,mean
64	dimension = np.array(input).ndim
65	shape = np.array(input).shape
66	#print 'd1,d2,d3,d4: ',d1,d2,d3,d4
67	print 'IN REDUCEFIELD dim,shape: ',dimension,shape
68	output = input
69	error = False
70	#### this is needed to cope the case where d4,d3,d2,d1 are single integers and not arrays
71	if d4 is not None and not isinstance(d4, np.ndarray): d4=[d4]
72	if d3 is not None and not isinstance(d3, np.ndarray): d3=[d3]
73	if d2 is not None and not isinstance(d2, np.ndarray): d2=[d2]
74	if d1 is not None and not isinstance(d1, np.ndarray): d1=[d1]
75	### now the main part
76	if dimension == 2:
77	if d2 >= shape[0]: error = True
78	elif d1 >= shape[1]: error = True
79	elif d1 is not None and d2 is not None: output = mean(input[d2,:],axis=0); output = mean(output[d1],axis=0)
80	elif d1 is not None: output = mean(input[:,d1],axis=1)
81	elif d2 is not None: output = mean(input[d2,:],axis=0)
82	elif dimension == 3:
83	if max(d4) >= shape[0]: error = True
84	elif max(d2) >= shape[1]: error = True
85	elif max(d1) >= shape[2]: error = True
86	elif d4 is not None and d2 is not None and d1 is not None:
87	output = mean(input[d4,:,:],axis=0); output = mean(output[d2,:],axis=0); output = mean(output[d1],axis=0)
88	elif d4 is not None and d2 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[d2,:],axis=0)
89	elif d4 is not None and d1 is not None: output = mean(input[d4,:,:],axis=0); output=mean(output[:,d1],axis=1)
90	elif d2 is not None and d1 is not None: output = mean(input[:,d2,:],axis=1); output=mean(output[:,d1],axis=1)
91	elif d1 is not None: output = mean(input[:,:,d1],axis=2)
92	elif d2 is not None: output = mean(input[:,d2,:],axis=1)
93	elif d4 is not None: output = mean(input[d4,:,:],axis=0)
94	elif dimension == 4:
95	if max(d4) >= shape[0]: error = True
96	elif max(d3) >= shape[1]: error = True
97	elif max(d2) >= shape[2]: error = True
98	elif max(d1) >= shape[3]: error = True
99	elif d4 is not None and d3 is not None and d2 is not None and d1 is not None:
100	output = mean(input[d4,:,:,:],axis=0)
101	output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
102	output = mean(output[d2,:],axis=0)
103	output = mean(output[d1],axis=0)
104	elif d4 is not None and d3 is not None and d2 is not None:
105	output = mean(input[d4,:,:,:],axis=0)
106	output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
107	output = mean(output[d2,:],axis=0)
108	elif d4 is not None and d3 is not None and d1 is not None:
109	output = mean(input[d4,:,:,:],axis=0)
110	output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
111	output = mean(output[:,d1],axis=1)
112	elif d4 is not None and d2 is not None and d1 is not None:
113	output = mean(input[d4,:,:,:],axis=0)
114	output = mean(output[:,d2,:],axis=1)
115	output = mean(output[:,d1],axis=1)
116	elif d3 is not None and d2 is not None and d1 is not None:
117	output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
118	output = mean(output[:,d2,:],axis=1)
119	output = mean(output[:,d1],axis=1)
120	elif d4 is not None and d3 is not None:
121	output = mean(input[d4,:,:,:],axis=0)
122	output = reduce_zaxis(output[d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
123	elif d4 is not None and d2 is not None:
124	output = mean(input[d4,:,:,:],axis=0)
125	output = mean(output[:,d2,:],axis=1)
126	elif d4 is not None and d1 is not None:
127	output = mean(input[d4,:,:,:],axis=0)
128	output = mean(output[:,:,d1],axis=2)
129	elif d3 is not None and d2 is not None:
130	output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
131	output = mean(output[:,d2,:],axis=1)
132	elif d3 is not None and d1 is not None:
133	output = reduce_zaxis(input[:,d3,:,:],ax=1,yint=yint,vert=alt,indice=d3)
134	output = mean(output[:,:,d1],axis=0)
135	elif d2 is not None and d1 is not None:
136	output = mean(input[:,:,d2,:],axis=2)
137	output = mean(output[:,:,d1],axis=2)
138	elif d1 is not None: output = mean(input[:,:,:,d1],axis=3)
139	elif d2 is not None: output = mean(input[:,:,d2,:],axis=2)
140	elif d3 is not None: output = reduce_zaxis(input[:,d3,:,:],ax=0,yint=yint,vert=alt,indice=d3)
141	elif d4 is not None: output = mean(input[d4,:,:,:],axis=0)
142	dimension = np.array(output).ndim
143	shape = np.array(output).shape
144	print 'OUT REDUCEFIELD dim,shape: ',dimension,shape
145	return output, error
146
147	## Author: AC + AS
148	def reduce_zaxis (input,ax=None,yint=False,vert=None,indice=None):
149	from mymath import max,mean
150	from scipy import integrate
151	if yint and vert is not None and indice is not None:
152	if type(input).__name__=='MaskedArray':
153	input.set_fill_value([np.NaN])
154	output = integrate.trapz(input.filled(),x=vert[indice],axis=ax)
155	else:
156	output = integrate.trapz(input.filled(),x=vert[indice],axis=ax)
157	else:
158	output = mean(input,axis=ax)
159	return output
160
161	## Author: AS + TN
162	def definesubplot ( numplot, fig ):
163	from matplotlib.pyplot import rcParams
164	rcParams['font.size'] = 12. ## default (important for multiple calls)
165	if numplot <= 0:
166	subv = 99999
167	subh = 99999
168	elif numplot == 1:
169	subv = 99999
170	subh = 99999
171	elif numplot == 2:
172	subv = 1
173	subh = 2
174	fig.subplots_adjust(wspace = 0.35)
175	rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. )
176	elif numplot == 3:
177	subv = 2
178	subh = 2
179	fig.subplots_adjust(wspace = 0.5)
180	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
181	elif numplot == 4:
182	subv = 2
183	subh = 2
184	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
185	rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. )
186	elif numplot <= 6:
187	subv = 2
188	subh = 3
189	fig.subplots_adjust(wspace = 0.4, hspace = 0.0)
190	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
191	elif numplot <= 8:
192	subv = 2
193	subh = 4
194	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
195	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
196	elif numplot <= 9:
197	subv = 3
198	subh = 3
199	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
200	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
201	elif numplot <= 12:
202	subv = 3
203	subh = 4
204	fig.subplots_adjust(wspace = 0.1, hspace = 0.1)
205	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
206	elif numplot <= 16:
207	subv = 4
208	subh = 4
209	fig.subplots_adjust(wspace = 0.3, hspace = 0.3)
210	rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. )
211	else:
212	print "number of plot supported: 1 to 16"
213	exit()
214	return subv,subh
215
216	## Author: AS
217	def getstralt(nc,nvert):
218	typefile = whatkindfile(nc)
219	if typefile is 'meso':
220	stralt = "_lvl" + str(nvert)
221	elif typefile is 'mesoapi':
222	zelevel = int(nc.variables['vert'][nvert])
223	if abs(zelevel) < 10000.: strheight=str(zelevel)+"m"
224	else: strheight=str(int(zelevel/1000.))+"km"
225	if 'altitude' in nc.dimensions: stralt = "_"+strheight+"-AMR"
226	elif 'altitude_abg' in nc.dimensions: stralt = "_"+strheight+"-ALS"
227	elif 'bottom_top' in nc.dimensions: stralt = "_"+strheight
228	elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa"
229	else: stralt = ""
230	else:
231	stralt = ""
232	return stralt
233
234	## Author: AS
235	def getlschar ( namefile ):
236	from netCDF4 import Dataset
237	from timestuff import sol2ls
238	from numpy import array
239	nc = Dataset(namefile)
240	zetime = None
241	if 'Times' in nc.variables:
242	zetime = nc.variables['Times'][0]
243	shape = array(nc.variables['Times']).shape
244	if shape[0] < 2: zetime = None
245	if zetime is not None \
246	and 'vert' not in nc.variables:
247	#### strangely enough this does not work for api or ncrcat results!
248	zetimestart = getattr(nc, 'START_DATE')
249	zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9])
250	if zeday < 0: lschar="" ## might have crossed a month... fix soon
251	else: lschar="_Ls"+str( int( 10. * sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) / 10. )
252	###
253	zetime2 = nc.variables['Times'][1]
254	one = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37.
255	next = int(zetime2[11]+zetime2[12]) + int(zetime2[14]+zetime2[15])/37.
256	zehour = one
257	zehourin = abs ( next - one )
258	else:
259	lschar=""
260	zehour = 0
261	zehourin = 1
262	return lschar, zehour, zehourin
263
264	## Author: AS
265	def getprefix (nc):
266	prefix = 'LMD_MMM_'
267	prefix = prefix + 'd'+str(getattr(nc,'GRID_ID'))+'_'
268	prefix = prefix + str(int(getattr(nc,'DX')/1000.))+'km_'
269	return prefix
270
271	## Author: AS
272	def getproj (nc):
273	typefile = whatkindfile(nc)
274	if typefile in ['mesoapi','meso','geo']:
275	### (il faudrait passer CEN_LON dans la projection ?)
276	map_proj = getattr(nc, 'MAP_PROJ')
277	cen_lat = getattr(nc, 'CEN_LAT')
278	if map_proj == 2:
279	if cen_lat > 10.:
280	proj="npstere"
281	#print "NP stereographic polar domain"
282	else:
283	proj="spstere"
284	#print "SP stereographic polar domain"
285	elif map_proj == 1:
286	#print "lambert projection domain"
287	proj="lcc"
288	elif map_proj == 3:
289	#print "mercator projection"
290	proj="merc"
291	else:
292	proj="merc"
293	elif typefile in ['gcm']: proj="cyl" ## pb avec les autres (de trace derriere la sphere ?)
294	else: proj="ortho"
295	return proj
296
297	## Author: AS
298	def ptitle (name):
299	from matplotlib.pyplot import title
300	title(name)
301	print name
302
303	## Author: AS
304	def polarinterv (lon2d,lat2d):
305	import numpy as np
306	wlon = [np.min(lon2d),np.max(lon2d)]
307	ind = np.array(lat2d).shape[0] / 2 ## to get a good boundlat and to get the pole
308	wlat = [np.min(lat2d[ind,:]),np.max(lat2d[ind,:])]
309	return [wlon,wlat]
310
311	## Author: AS
312	def simplinterv (lon2d,lat2d):
313	import numpy as np
314	return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]]
315
316	## Author: AS
317	def wrfinterv (lon2d,lat2d):
318	nx = len(lon2d[0,:])-1
319	ny = len(lon2d[:,0])-1
320	lon1 = lon2d[0,0]
321	lon2 = lon2d[nx,ny]
322	lat1 = lat2d[0,0]
323	lat2 = lat2d[nx,ny]
324	if abs(0.5(lat1+lat2)) > 60.: wider = 0.5 (abs(lon1)+abs(lon2)) * 0.1
325	else: wider = 0.
326	if lon1 < lon2: wlon = [lon1, lon2 + wider]
327	else: wlon = [lon2, lon1 + wider]
328	if lat1 < lat2: wlat = [lat1, lat2]
329	else: wlat = [lat2, lat1]
330	return [wlon,wlat]
331
332	## Author: AS
333	def makeplotres (filename,res=None,pad_inches_value=0.25,folder='',disp=True,ext='png',erase=False):
334	import matplotlib.pyplot as plt
335	from os import system
336	addstr = ""
337	if res is not None:
338	res = int(res)
339	addstr = "_"+str(res)
340	name = filename+addstr+"."+ext
341	if folder != '': name = folder+'/'+name
342	plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value)
343	if disp: display(name)
344	if ext in ['eps','ps','svg']: system("tar czvf "+name+".tar.gz "+name+" ; rm -f "+name)
345	if erase: system("mv "+name+" to_be_erased")
346	return
347
348	## Author: AS
349	def dumpbdy (field,n,stag=None):
350	nx = len(field[0,:])-1
351	ny = len(field[:,0])-1
352	if stag == 'U': nx = nx-1
353	if stag == 'V': ny = ny-1
354	return field[n:ny-n,n:nx-n]
355
356	## Author: AS
357	def getcoorddef ( nc ):
358	import numpy as np
359	## getcoord2d for predefined types
360	typefile = whatkindfile(nc)
361	if typefile in ['mesoapi','meso']:
362	[lon2d,lat2d] = getcoord2d(nc)
363	lon2d = dumpbdy(lon2d,6)
364	lat2d = dumpbdy(lat2d,6)
365	elif typefile in ['gcm']:
366	[lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True)
367	elif typefile in ['geo']:
368	[lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M')
369	return lon2d,lat2d
370
371	## Author: AS
372	def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False):
373	import numpy as np
374	if is1d:
375	lat = nc.variables[nlat][:]
376	lon = nc.variables[nlon][:]
377	[lon2d,lat2d] = np.meshgrid(lon,lat)
378	else:
379	lat = nc.variables[nlat][0,:,:]
380	lon = nc.variables[nlon][0,:,:]
381	[lon2d,lat2d] = [lon,lat]
382	return lon2d,lat2d
383
384	## Author: AS
385	def smooth (field, coeff):
386	## actually blur_image could work with different coeff on x and y
387	if coeff > 1: result = blur_image(field,int(coeff))
388	else: result = field
389	return result
390
391	## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
392	def gauss_kern(size, sizey=None):
393	import numpy as np
394	# Returns a normalized 2D gauss kernel array for convolutions
395	size = int(size)
396	if not sizey:
397	sizey = size
398	else:
399	sizey = int(sizey)
400	x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
401	g = np.exp(-(x2/float(size)+y2/float(sizey)))
402	return g / g.sum()
403
404	## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth
405	def blur_image(im, n, ny=None) :
406	from scipy.signal import convolve
407	# blurs the image by convolving with a gaussian kernel of typical size n.
408	# The optional keyword argument ny allows for a different size in the y direction.
409	g = gauss_kern(n, sizey=ny)
410	improc = convolve(im, g, mode='same')
411	return improc
412
413	## Author: AS
414	def getwinddef (nc):
415	## getwinds for predefined types
416	typefile = whatkindfile(nc)
417	###
418	if typefile is 'mesoapi': [uchar,vchar] = ['Um','Vm']
419	elif typefile is 'gcm': [uchar,vchar] = ['u','v']
420	elif typefile is 'meso': [uchar,vchar] = ['U','V']
421	else: [uchar,vchar] = ['not found','not found']
422	###
423	if typefile in ['meso']: metwind = False ## geometrical (wrt grid)
424	else: metwind = True ## meteorological (zon/mer)
425	if metwind is False: print "Not using meteorological winds. You trust numerical grid as being (x,y)"
426	###
427	return uchar,vchar,metwind
428
429	## Author: AS
430	def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True):
431	## scale regle la reference du vecteur
432	## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs.
433	import matplotlib.pyplot as plt
434	import numpy as np
435	posx = np.min(x) - np.std(x) / 10.
436	posy = np.min(y) - np.std(y) / 10.
437	u = smooth(u,csmooth)
438	v = smooth(v,csmooth)
439	widthvec = 0.003 #0.005 #0.003
440	q = plt.quiver( x[::stride,::stride],\
441	y[::stride,::stride],\
442	u[::stride,::stride],\
443	v[::stride,::stride],\
444	angles='xy',color=color,pivot='middle',\
445	scale=factor,width=widthvec )
446	if color in ['white','yellow']: kcolor='black'
447	else: kcolor=color
448	if key: p = plt.quiverkey(q,posx,posy,scale,\
449	str(int(scale)),coordinates='data',color=kcolor,labelpos='S',labelsep = 0.03)
450	return
451
452	## Author: AS
453	def display (name):
454	from os import system
455	system("display "+name+" > /dev/null 2> /dev/null &")
456	return name
457
458	## Author: AS
459	def findstep (wlon):
460	steplon = int((wlon[1]-wlon[0])/4.) #3
461	step = 120.
462	while step > steplon and step > 15. : step = step / 2.
463	if step <= 15.:
464	while step > steplon and step > 5. : step = step - 5.
465	if step <= 5.:
466	while step > steplon and step > 1. : step = step - 1.
467	if step <= 1.:
468	step = 1.
469	return step
470
471	## Author: AS
472	def define_proj (char,wlon,wlat,back=None,blat=None):
473	from mpl_toolkits.basemap import Basemap
474	import numpy as np
475	import matplotlib as mpl
476	from mymath import max
477	meanlon = 0.5*(wlon[0]+wlon[1])
478	meanlat = 0.5*(wlat[0]+wlat[1])
479	if blat is None:
480	if wlat[0] >= 80.: blat = 40.
481	elif wlat[1] <= -80.: blat = -40.
482	elif wlat[1] >= 0.: blat = wlat[0]
483	elif wlat[0] <= 0.: blat = wlat[1]
484	#print "blat ", blat
485	h = 50. ## en km
486	radius = 3397200.
487	if char == "cyl": m = Basemap(rsphere=radius,projection='cyl',\
488	llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
489	elif char == "moll": m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon)
490	elif char == "ortho": m = Basemap(rsphere=radius,projection='ortho',lon_0=meanlon,lat_0=meanlat)
491	elif char == "lcc": m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\
492	llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
493	elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.)
494	elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=0.)
495	elif char == "nplaea": m = Basemap(rsphere=radius,projection='nplaea', boundinglat=wlat[0], lon_0=meanlon)
496	elif char == "laea": m = Basemap(rsphere=radius,projection='laea',lon_0=meanlon,lat_0=meanlat,lat_ts=meanlat,\
497	llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
498	elif char == "nsper": m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.)
499	elif char == "merc": m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\
500	llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1])
501	fontsizemer = int(mpl.rcParams['font.size']*3./4.)
502	if char in ["cyl","lcc","merc","nsper","laea"]: step = findstep(wlon)
503	else: step = 10.
504	steplon = step*2.
505	#if back in ["geolocal"]:
506	# step = np.min([5.,step])
507	# steplon = step
508	m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color='grey', fontsize=fontsizemer)
509	m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer)
510	if back: m.warpimage(marsmap(back),scale=0.75)
511	#if not back:
512	# if not var: back = "mola" ## if no var: draw mola
513	# elif typefile in ['mesoapi','meso','geo'] \
514	# and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw
515	# else: pass ## else: draw None
516	return m
517
518	## Author: AS
519	#### test temporaire
520	def putpoints (map,plot):
521	#### from http://www.scipy.org/Cookbook/Matplotlib/Maps
522	# lat/lon coordinates of five cities.
523	lats = [18.4]
524	lons = [-134.0]
525	points=['Olympus Mons']
526	# compute the native map projection coordinates for cities.
527	x,y = map(lons,lats)
528	# plot filled circles at the locations of the cities.
529	map.plot(x,y,'bo')
530	# plot the names of those five cities.
531	wherept = 0 #1000 #50000
532	for name,xpt,ypt in zip(points,x,y):
533	plot.text(xpt+wherept,ypt+wherept,name)
534	## le nom ne s'affiche pas...
535	return
536
537	## Author: AS
538	def calculate_bounds(field,vmin=None,vmax=None):
539	import numpy as np
540	from mymath import max,min,mean
541	ind = np.where(field < 9e+35)
542	fieldcalc = field[ ind ] # la syntaxe compacte ne marche si field est un tuple
543	###
544	dev = np.std(fieldcalc)*3.0
545	###
546	if vmin is None:
547	zevmin = mean(fieldcalc) - dev
548	else: zevmin = vmin
549	###
550	if vmax is None: zevmax = mean(fieldcalc) + dev
551	else: zevmax = vmax
552	if vmin == vmax:
553	zevmin = mean(fieldcalc) - dev ### for continuity
554	zevmax = mean(fieldcalc) + dev ### for continuity
555	###
556	if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0.
557	print "BOUNDS field ", min(fieldcalc), max(fieldcalc)
558	print "BOUNDS adopted ", zevmin, zevmax
559	return zevmin, zevmax
560
561	## Author: AS
562	def bounds(what_I_plot,zevmin,zevmax):
563	from mymath import max,min,mean
564	### might be convenient to add the missing value in arguments
565	#what_I_plot[ what_I_plot < zevmin ] = zevmin#*(1. + 1.e-7)
566	if zevmin < 0: what_I_plot[ what_I_plot < zevmin(1. - 1.e-7) ] = zevmin(1. - 1.e-7)
567	else: what_I_plot[ what_I_plot < zevmin(1. + 1.e-7) ] = zevmin(1. + 1.e-7)
568	print "NEW MIN ", min(what_I_plot)
569	what_I_plot[ what_I_plot > 9e+35 ] = -9e+35
570	what_I_plot[ what_I_plot > zevmax ] = zevmax
571	print "NEW MAX ", max(what_I_plot)
572	return what_I_plot
573
574	## Author: AS
575	def nolow(what_I_plot):
576	from mymath import max,min
577	lim = 0.150.5(abs(max(what_I_plot))+abs(min(what_I_plot)))
578	print "NO PLOT BELOW VALUE ", lim
579	what_I_plot [ abs(what_I_plot) < lim ] = 1.e40
580	return what_I_plot
581
582	## Author: AS
583	def zoomset (wlon,wlat,zoom):
584	dlon = abs(wlon[1]-wlon[0])/2.
585	dlat = abs(wlat[1]-wlat[0])/2.
586	[wlon,wlat] = [ [wlon[0]+zoomdlon/100.,wlon[1]-zoomdlon/100.],\
587	[wlat[0]+zoomdlat/100.,wlat[1]-zoomdlat/100.] ]
588	print "ZOOM %",zoom,wlon,wlat
589	return wlon,wlat
590
591	## Author: AS
592	def fmtvar (whichvar="def"):
593	fmtvar = { \
594	"tk": "%.0f",\
595	"T_NADIR_DAY": "%.0f",\
596	"T_NADIR_NIT": "%.0f",\
597	"tpot": "%.0f",\
598	"TSURF": "%.0f",\
599	"def": "%.1e",\
600	"PTOT": "%.0f",\
601	"HGT": "%.1e",\
602	"USTM": "%.2f",\
603	"HFX": "%.0f",\
604	"ICETOT": "%.1e",\
605	"TAU_ICE": "%.2f",\
606	"VMR_ICE": "%.1e",\
607	"MTOT": "%.1f",\
608	"anomaly": "%.1f",\
609	"W": "%.1f",\
610	"WMAX_TH": "%.1f",\
611	"QSURFICE": "%.0f",\
612	"Um": "%.0f",\
613	"ALBBARE": "%.2f",\
614	"TAU": "%.1f",\
615	"CO2": "%.2f",\
616	#### T.N.
617	"TEMP": "%.0f",\
618	"VMR_H2OICE": "%.0f",\
619	"VMR_H2OVAP": "%.0f",\
620	"TAUTES": "%.2f",\
621	"TAUTESAP": "%.2f",\
622
623	}
624	if whichvar not in fmtvar:
625	whichvar = "def"
626	return fmtvar[whichvar]
627
628	## Author: AS
629	####################################################################################################################
630	### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
631	def defcolorb (whichone="def"):
632	whichcolorb = { \
633	"def": "spectral",\
634	"HGT": "spectral",\
635	"tk": "gist_heat",\
636	"TSURF": "RdBu_r",\
637	"QH2O": "PuBu",\
638	"USTM": "YlOrRd",\
639	#"T_nadir_nit": "RdBu_r",\
640	#"T_nadir_day": "RdBu_r",\
641	"HFX": "RdYlBu",\
642	"ICETOT": "YlGnBu_r",\
643	#"MTOT": "PuBu",\
644	"CCNQ": "YlOrBr",\
645	"CCNN": "YlOrBr",\
646	"TEMP": "Jet",\
647	"TAU_ICE": "Blues",\
648	"VMR_ICE": "Blues",\
649	"W": "jet",\
650	"WMAX_TH": "spectral",\
651	"anomaly": "RdBu_r",\
652	"QSURFICE": "hot_r",\
653	"ALBBARE": "spectral",\
654	"TAU": "YlOrBr_r",\
655	"CO2": "YlOrBr_r",\
656	#### T.N.
657	"MTOT": "Jet",\
658	"H2O_ICE_S": "RdBu",\
659	"VMR_H2OICE": "PuBu",\
660	"VMR_H2OVAP": "PuBu",\
661	}
662	#W --> spectral ou jet
663	#spectral BrBG RdBu_r
664	#print "predefined colorbars"
665	if whichone not in whichcolorb:
666	whichone = "def"
667	return whichcolorb[whichone]
668
669	## Author: AS
670	def definecolorvec (whichone="def"):
671	whichcolor = { \
672	"def": "black",\
673	"vis": "yellow",\
674	"vishires": "yellow",\
675	"molabw": "yellow",\
676	"mola": "black",\
677	"gist_heat": "white",\
678	"hot": "tk",\
679	"gist_rainbow": "black",\
680	"spectral": "black",\
681	"gray": "red",\
682	"PuBu": "black",\
683	}
684	if whichone not in whichcolor:
685	whichone = "def"
686	return whichcolor[whichone]
687
688	## Author: AS
689	def marsmap (whichone="vishires"):
690	from os import uname
691	mymachine = uname()[1]
692	### not sure about speed-up with this method... looks the same
693	if "lmd.jussieu.fr" in mymachine: domain = "/u/aslmd/WWW/maps/"
694	else: domain = "http://www.lmd.jussieu.fr/~aslmd/maps/"
695	whichlink = { \
696	#"vis": "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\
697	#"vishires": "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\
698	#"geolocal": "http://dl.dropbox.com/u/11078310/geolocal.jpg",\
699	#"mola": "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\
700	#"molabw": "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\
701	"vis": domain+"mar0kuu2.jpg",\
702	"vishires": domain+"MarsMap_2500x1250.jpg",\
703	"geolocal": domain+"geolocal.jpg",\
704	"mola": domain+"mars-mola-2k.jpg",\
705	"molabw": domain+"MarsElevation_2500x1250.jpg",\
706	"clouds": "http://www.johnstonsarchive.net/spaceart/marswcloudmap.jpg",\
707	"jupiter": "http://www.mmedia.is/~bjj/data/jupiter_css/jupiter_css.jpg",\
708	"jupiter_voy": "http://www.mmedia.is/~bjj/data/jupiter/jupiter_vgr2.jpg",\
709	"bw": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg",\
710	"contrast": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg",\
711	"nice": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg",\
712	"blue": "http://eoimages.gsfc.nasa.gov/ve/2430/land_ocean_ice_2048.jpg",\
713	"blueclouds": "http://eoimages.gsfc.nasa.gov/ve/2431/land_ocean_ice_cloud_2048.jpg",\
714	"justclouds": "http://eoimages.gsfc.nasa.gov/ve/2432/cloud_combined_2048.jpg",\
715	}
716	### see http://www.mmedia.is/~bjj/planetary_maps.html
717	if whichone not in whichlink:
718	print "marsmap: choice not defined... you'll get the default one... "
719	whichone = "vishires"
720	return whichlink[whichone]
721
722	#def earthmap (whichone):
723	# if whichone == "contrast": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg"
724	# elif whichone == "bw": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg"
725	# elif whichone == "nice": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg"
726	# return whichlink
727
728	## Author: AS
729	def latinterv (area="Whole"):
730	list = { \
731	"Europe": [[ 20., 80.],[- 50., 50.]],\
732	"Central_America": [[-10., 40.],[ 230., 300.]],\
733	"Africa": [[-20., 50.],[- 50., 50.]],\
734	"Whole": [[-90., 90.],[-180., 180.]],\
735	"Southern_Hemisphere": [[-90., 60.],[-180., 180.]],\
736	"Northern_Hemisphere": [[-60., 90.],[-180., 180.]],\
737	"Tharsis": [[-30., 60.],[-170.,- 10.]],\
738	"Whole_No_High": [[-60., 60.],[-180., 180.]],\
739	"Chryse": [[-60., 60.],[- 60., 60.]],\
740	"North_Pole": [[ 50., 90.],[-180., 180.]],\
741	"Close_North_Pole": [[ 75., 90.],[-180., 180.]],\
742	"Far_South_Pole": [[-90.,-40.],[-180., 180.]],\
743	"South_Pole": [[-90.,-50.],[-180., 180.]],\
744	"Close_South_Pole": [[-90.,-75.],[-180., 180.]],\
745	}
746	if area not in list: area = "Whole"
747	[olat,olon] = list[area]
748	return olon,olat
749
750	## Author: TN
751	def separatenames (name):
752	from numpy import concatenate
753	# look for comas in the input name to separate different names (files, variables,etc ..)
754	if name is None:
755	names = None
756	else:
757	names = []
758	stop = 0
759	currentname = name
760	while stop == 0:
761	indexvir = currentname.find(',')
762	if indexvir == -1:
763	stop = 1
764	name1 = currentname
765	else:
766	name1 = currentname[0:indexvir]
767	names = concatenate((names,[name1]))
768	currentname = currentname[indexvir+1:len(currentname)]
769	return names
770
771	## Author: TN [old]
772	def getopmatrix (kind,n):
773	import numpy as np
774	# compute matrix of operations between files
775	# the matrix is 'number of files'-square
776	# 1: difference (row minus column), 2: add
777	# not 0 in diag : just plot
778	if n == 1:
779	opm = np.eye(1)
780	elif kind == 'basic':
781	opm = np.eye(n)
782	elif kind == 'difference1': # show differences with 1st file
783	opm = np.zeros((n,n))
784	opm[0,:] = 1
785	opm[0,0] = 0
786	elif kind == 'difference2': # show differences with 1st file AND show 1st file
787	opm = np.zeros((n,n))
788	opm[0,:] = 1
789	else:
790	opm = np.eye(n)
791	return opm
792
793	## Author: TN [old]
794	def checkcoherence (nfiles,nlat,nlon,ntime):
795	if (nfiles > 1):
796	if (nlat > 1):
797	errormess("what you asked is not possible !")
798	return 1
799
800	## Author: TN
801	def readslices(saxis):
802	from numpy import empty
803	if saxis == None:
804	zesaxis = None
805	else:
806	zesaxis = empty((len(saxis),2))
807	for i in range(len(saxis)):
808	a = separatenames(saxis[i])
809	if len(a) == 1:
810	zesaxis[i,:] = float(a[0])
811	else:
812	zesaxis[i,0] = float(a[0])
813	zesaxis[i,1] = float(a[1])
814
815	return zesaxis
816
817	## Author: TN
818	def getsindex(saxis,index,axis):
819	# input : all the desired slices and the good index
820	# output : all indexes to be taken into account for reducing field
821	import numpy as np
822	if ( np.array(axis).ndim == 2):
823	axis = axis[:,0]
824	if saxis is None:
825	zeindex = None
826	else:
827	aaa = int(np.argmin(abs(saxis[index,0] - axis)))
828	bbb = int(np.argmin(abs(saxis[index,1] - axis)))
829	[imin,imax] = np.sort(np.array([aaa,bbb]))
830	zeindex = np.array(range(imax-imin+1))+imin
831	# because -180 and 180 are the same point in longitude,
832	# we get rid of one for averaging purposes.
833	if axis[imin] == -180 and axis[imax] == 180:
834	zeindex = zeindex[0:len(zeindex)-1]
835	print "INFO: whole longitude averaging asked, so last point is not taken into account."
836	return zeindex
837
838	## Author: TN
839	def define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,indextime,what_I_plot,dim0,vertmode):
840	# Purpose of define_axis is to find x and y axis scales in a smart way
841	# x axis priority: 1/time 2/lon 3/lat 4/vertical
842	# To be improved !!!...
843	from numpy import array,swapaxes
844	x = None
845	y = None
846	count = 0
847	what_I_plot = array(what_I_plot)
848	shape = what_I_plot.shape
849	if indextime is None:
850	print "AXIS is time"
851	x = time
852	count = count+1
853	if indexlon is None:
854	print "AXIS is lon"
855	if count == 0: x = lon
856	else: y = lon
857	count = count+1
858	if indexlat is None:
859	print "AXIS is lat"
860	if count == 0: x = lat
861	else: y = lat
862	count = count+1
863	if indexvert is None and dim0 is 4:
864	print "AXIS is vert"
865	if vertmode == 0: # vertical axis is as is (GCM grid)
866	if count == 0: x=range(len(vert))
867	else: y=range(len(vert))
868	count = count+1
869	else: # vertical axis is in kms
870	if count == 0: x = vert
871	else: y = vert
872	count = count+1
873	x = array(x)
874	y = array(y)
875	print "CHECK: what_I_plot.shape", what_I_plot.shape
876	print "CHECK: x.shape, y.shape", x.shape, y.shape
877	if len(shape) == 1:
878	print shape[0]
879	if shape[0] != len(x):
880	print "WARNING HERE !!!"
881	x = y
882	elif len(shape) == 2:
883	print shape[1], len(y), shape[0], len(x)
884	if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]:
885	what_I_plot = swapaxes(what_I_plot,0,1)
886	print "INFO: swapaxes", what_I_plot.shape, shape
887	#x0 = x
888	#x = y
889	#y = x0
890	#print "define_axis", x, y
891	return what_I_plot,x,y
892
893	# Author: TN + AS
894	def determineplot(slon, slat, svert, stime):
895	nlon = 1 # number of longitudinal slices -- 1 is None
896	nlat = 1
897	nvert = 1
898	ntime = 1
899	nslices = 1
900	if slon is not None:
901	nslices = nslices*len(slon)
902	nlon = len(slon)
903	if slat is not None:
904	nslices = nslices*len(slat)
905	nlat = len(slat)
906	if svert is not None:
907	nslices = nslices*len(svert)
908	nvert = len(svert)
909	if stime is not None:
910	nslices = nslices*len(stime)
911	ntime = len(stime)
912	#else:
913	# nslices = 2
914
915	mapmode = 0
916	if slon is None and slat is None:
917	mapmode = 1 # in this case we plot a map, with the given projection
918
919	return nlon, nlat, nvert, ntime, mapmode, nslices

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

Download in other formats: