Context Navigation

myplot.py @ 363

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

M 362 planetoplot.py
M 362 meso.py
M 362 gcm.py
^{----------------- Added --inverty option to force pyplot to inverse the y axis of a contourplot. Usefull when plotting data along pressure vertical axis.}

M 362 zrecast_wrapper.py
^{----------------- Some modif for the pressure interpolation. A following commit will make it cleaner so that interpolation details can be passed to zrecast from gcm.py}

M 362 myplot.py
^{----------------- Minor modifs, added some variables.}

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

Download in other formats: