Context Navigation

myplot.py @ 395

Last change on this file since 395 was 395, checked in by acolaitis, 13 years ago

PYTHON

M 393 mymath.py
^{----------------- Cosmetic change}

M 393 make_netcdf.py
^{----------------- Cosmetic change}

M 393 planetoplot.py
^{----------------- corrected bug with varname and --tsat}

M 393 myplot.py
^{----------------- added possibility for the script to read ncdf files with NaN value and treat them correctly without messing with the mathematical operations.}

File size: 36.2 KB

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

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

Context Navigation

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

Download in other formats: