Context Navigation

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

Last change on this file since 448 was 448, checked in by aslmd, 13 years ago
GRAPHICS: coded a more rational way for contours and shaded plots. that is, no more need for call_contour and no more (or not too much) duplicates between contour part and shaded part. working well with static or animated plots.
File size: 42.8 KB

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

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

Download in other formats: