Context Navigation

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

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

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

Download in other formats: