Context Navigation

myplot.py @ 376

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

Added new routine to our python subroutines pool: make_netcdf.py

This routine is an easy way to create a netcdf file conform to the gcm .nc format.
Here is an exemple of a call to the routine:

#####################

from make_netcdf import make_gcm_netcdf

make_gcm_netcdf (zfilename="diagfi_TES.nc", \

zdescription="TES day and night temperature fields", \
zlon=lon, \
zlat=lat, \
zalt=alt, \
ztime=time, \
zvariables=[ps_day,ps_night,aps,bps,temp_day,temp_night], \
znames=["ps_day","ps_night","aps","bps","temp_day","temp_night"])

####################

This will create a file diagfi_TES.nc with variables specified in "zvariables", with names specified in "znames", along the dimensions given in "zlon,zlat,zalt,ztime". One can use less dimensions if the variables to be written are not 4D.

Note: by "gcm .nc format" is implied that dimension names are conform to what is expected when using zrecast/hrecast tools, or gcm.py. One could write a similar routine with variables conform to API and meso.py for mesoscale .nc files.

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",\
[376]	561	"T_NADIR_NIT": "%.0f",\
[204]	562	"tpot": "%.0f",\
[295]	563	"TSURF": "%.0f",\
[204]	564	"def": "%.1e",\
	565	"PTOT": "%.0f",\
	566	"HGT": "%.1e",\
	567	"USTM": "%.2f",\
[225]	568	"HFX": "%.0f",\
[232]	569	"ICETOT": "%.1e",\
[237]	570	"TAU_ICE": "%.2f",\
[252]	571	"VMR_ICE": "%.1e",\
[345]	572	"MTOT": "%.1f",\
[240]	573	"anomaly": "%.1f",\
[241]	574	"W": "%.1f",\
[287]	575	"WMAX_TH": "%.1f",\
	576	"QSURFICE": "%.0f",\
[296]	577	"Um": "%.0f",\
[295]	578	"ALBBARE": "%.2f",\
[317]	579	"TAU": "%.1f",\
[345]	580	#### T.N.
	581	"TEMP": "%.0f",\
	582	"VMR_H2OICE": "%.0f",\
	583	"VMR_H2OVAP": "%.0f",\
	584	"TAUTES": "%.2f",\
	585	"TAUTESAP": "%.2f",\
	586
[204]	587	}
	588	if whichvar not in fmtvar:
	589	whichvar = "def"
	590	return fmtvar[whichvar]
[201]	591
[345]	592	## Author: AS
[233]	593	####################################################################################################################
	594	### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png
[202]	595	def defcolorb (whichone="def"):
[204]	596	whichcolorb = { \
	597	"def": "spectral",\
	598	"HGT": "spectral",\
	599	"tk": "gist_heat",\
[295]	600	"TSURF": "RdBu_r",\
[204]	601	"QH2O": "PuBu",\
	602	"USTM": "YlOrRd",\
[363]	603	#"T_nadir_nit": "RdBu_r",\
	604	#"T_nadir_day": "RdBu_r",\
[225]	605	"HFX": "RdYlBu",\
[310]	606	"ICETOT": "YlGnBu_r",\
[345]	607	#"MTOT": "PuBu",\
	608	"CCNQ": "YlOrBr",\
	609	"CCNN": "YlOrBr",\
	610	"TEMP": "Jet",\
[238]	611	"TAU_ICE": "Blues",\
[252]	612	"VMR_ICE": "Blues",\
[241]	613	"W": "jet",\
[287]	614	"WMAX_TH": "spectral",\
[240]	615	"anomaly": "RdBu_r",\
[287]	616	"QSURFICE": "hot_r",\
[295]	617	"ALBBARE": "spectral",\
[317]	618	"TAU": "YlOrBr_r",\
[345]	619	#### T.N.
	620	"MTOT": "Jet",\
	621	"H2O_ICE_S": "RdBu",\
	622	"VMR_H2OICE": "PuBu",\
	623	"VMR_H2OVAP": "PuBu",\
[204]	624	}
[241]	625	#W --> spectral ou jet
[240]	626	#spectral BrBG RdBu_r
[241]	627	print "predefined colorbars"
[204]	628	if whichone not in whichcolorb:
	629	whichone = "def"
	630	return whichcolorb[whichone]
[202]	631
[345]	632	## Author: AS
[202]	633	def definecolorvec (whichone="def"):
	634	whichcolor = { \
	635	"def": "black",\
	636	"vis": "yellow",\
	637	"vishires": "yellow",\
	638	"molabw": "yellow",\
	639	"mola": "black",\
	640	"gist_heat": "white",\
	641	"hot": "tk",\
	642	"gist_rainbow": "black",\
	643	"spectral": "black",\
	644	"gray": "red",\
	645	"PuBu": "black",\
	646	}
	647	if whichone not in whichcolor:
	648	whichone = "def"
	649	return whichcolor[whichone]
	650
[345]	651	## Author: AS
[180]	652	def marsmap (whichone="vishires"):
[233]	653	from os import uname
	654	mymachine = uname()[1]
	655	### not sure about speed-up with this method... looks the same
	656	if "lmd.jussieu.fr" in mymachine: domain = "/u/aslmd/WWW/maps/"
	657	else: domain = "http://www.lmd.jussieu.fr/~aslmd/maps/"
[180]	658	whichlink = { \
[233]	659	#"vis": "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\
	660	#"vishires": "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\
	661	#"geolocal": "http://dl.dropbox.com/u/11078310/geolocal.jpg",\
	662	#"mola": "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\
	663	#"molabw": "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\
	664	"vis": domain+"mar0kuu2.jpg",\
	665	"vishires": domain+"MarsMap_2500x1250.jpg",\
	666	"geolocal": domain+"geolocal.jpg",\
	667	"mola": domain+"mars-mola-2k.jpg",\
	668	"molabw": domain+"MarsElevation_2500x1250.jpg",\
[238]	669	"clouds": "http://www.johnstonsarchive.net/spaceart/marswcloudmap.jpg",\
	670	"jupiter": "http://www.mmedia.is/~bjj/data/jupiter_css/jupiter_css.jpg",\
	671	"jupiter_voy": "http://www.mmedia.is/~bjj/data/jupiter/jupiter_vgr2.jpg",\
[273]	672	"bw": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg",\
	673	"contrast": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg",\
	674	"nice": "http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg",\
	675	"blue": "http://eoimages.gsfc.nasa.gov/ve/2430/land_ocean_ice_2048.jpg",\
[296]	676	"blueclouds": "http://eoimages.gsfc.nasa.gov/ve/2431/land_ocean_ice_cloud_2048.jpg",\
	677	"justclouds": "http://eoimages.gsfc.nasa.gov/ve/2432/cloud_combined_2048.jpg",\
[180]	678	}
[238]	679	### see http://www.mmedia.is/~bjj/planetary_maps.html
[180]	680	if whichone not in whichlink:
	681	print "marsmap: choice not defined... you'll get the default one... "
	682	whichone = "vishires"
	683	return whichlink[whichone]
	684
[273]	685	#def earthmap (whichone):
	686	# if whichone == "contrast": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg"
	687	# elif whichone == "bw": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg"
	688	# elif whichone == "nice": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg"
	689	# return whichlink
[180]	690
[345]	691	## Author: AS
[241]	692	def latinterv (area="Whole"):
	693	list = { \
	694	"Europe": [[ 20., 80.],[- 50., 50.]],\
	695	"Central_America": [[-10., 40.],[ 230., 300.]],\
	696	"Africa": [[-20., 50.],[- 50., 50.]],\
[273]	697	"Whole": [[-90., 90.],[-180., 180.]],\
	698	"Southern_Hemisphere": [[-90., 60.],[-180., 180.]],\
	699	"Northern_Hemisphere": [[-60., 90.],[-180., 180.]],\
[241]	700	"Tharsis": [[-30., 60.],[-170.,- 10.]],\
	701	"Whole_No_High": [[-60., 60.],[-180., 180.]],\
	702	"Chryse": [[-60., 60.],[- 60., 60.]],\
	703	"North_Pole": [[ 50., 90.],[-180., 180.]],\
	704	"Close_North_Pole": [[ 75., 90.],[-180., 180.]],\
	705	"Far_South_Pole": [[-90.,-40.],[-180., 180.]],\
	706	"South_Pole": [[-90.,-50.],[-180., 180.]],\
	707	"Close_South_Pole": [[-90.,-75.],[-180., 180.]],\
	708	}
	709	if area not in list: area = "Whole"
	710	[olat,olon] = list[area]
	711	return olon,olat
	712
[345]	713	## Author: TN
	714	def separatenames (name):
	715	from numpy import concatenate
	716	# look for comas in the input name to separate different names (files, variables,etc ..)
	717	if name is None:
	718	names = None
	719	else:
	720	names = []
	721	stop = 0
	722	currentname = name
	723	while stop == 0:
	724	indexvir = currentname.find(',')
	725	if indexvir == -1:
	726	stop = 1
	727	name1 = currentname
	728	else:
	729	name1 = currentname[0:indexvir]
	730	names = concatenate((names,[name1]))
	731	currentname = currentname[indexvir+1:len(currentname)]
	732	return names
	733
	734	## Author: TN [old]
	735	def getopmatrix (kind,n):
	736	import numpy as np
	737	# compute matrix of operations between files
	738	# the matrix is 'number of files'-square
	739	# 1: difference (row minus column), 2: add
	740	# not 0 in diag : just plot
	741	if n == 1:
	742	opm = np.eye(1)
	743	elif kind == 'basic':
	744	opm = np.eye(n)
	745	elif kind == 'difference1': # show differences with 1st file
	746	opm = np.zeros((n,n))
	747	opm[0,:] = 1
	748	opm[0,0] = 0
	749	elif kind == 'difference2': # show differences with 1st file AND show 1st file
	750	opm = np.zeros((n,n))
	751	opm[0,:] = 1
	752	else:
	753	opm = np.eye(n)
	754	return opm
	755
	756	## Author: TN [old]
	757	def checkcoherence (nfiles,nlat,nlon,ntime):
	758	if (nfiles > 1):
	759	if (nlat > 1):
	760	errormess("what you asked is not possible !")
	761	return 1
	762
	763	## Author: TN
	764	def readslices(saxis):
	765	from numpy import empty
	766	if saxis == None:
	767	zesaxis = None
	768	else:
	769	zesaxis = empty((len(saxis),2))
	770	for i in range(len(saxis)):
	771	a = separatenames(saxis[i])
	772	if len(a) == 1:
	773	zesaxis[i,:] = float(a[0])
	774	else:
	775	zesaxis[i,0] = float(a[0])
	776	zesaxis[i,1] = float(a[1])
	777
	778	return zesaxis
	779
	780	## Author: TN
	781	def getsindex(saxis,index,axis):
	782	# input : all the desired slices and the good index
	783	# output : all indexes to be taken into account for reducing field
	784	import numpy as np
[349]	785	if ( np.array(axis).ndim == 2):
	786	axis = axis[:,0]
[345]	787	if saxis is None:
	788	zeindex = None
	789	else:
	790	aaa = int(np.argmin(abs(saxis[index,0] - axis)))
	791	bbb = int(np.argmin(abs(saxis[index,1] - axis)))
	792	[imin,imax] = np.sort(np.array([aaa,bbb]))
	793	zeindex = np.array(range(imax-imin+1))+imin
	794	# because -180 and 180 are the same point in longitude,
	795	# we get rid of one for averaging purposes.
	796	if axis[imin] == -180 and axis[imax] == 180:
	797	zeindex = zeindex[0:len(zeindex)-1]
	798	print "whole longitude averaging asked, so last point is not taken into account."
	799	return zeindex
	800
	801	## Author: TN
	802	def define_axis(lon,lat,vert,time,indexlon,indexlat,indexvert,indextime,what_I_plot,dim0,vertmode):
	803	# Purpose of define_axis is to find x and y axis scales in a smart way
	804	# x axis priority: 1/time 2/lon 3/lat 4/vertical
	805	# To be improved !!!...
	806	from numpy import array,swapaxes
	807	x = None
	808	y = None
	809	count = 0
	810	what_I_plot = array(what_I_plot)
	811	shape = what_I_plot.shape
	812	if indextime is None:
[350]	813	print "axis is time"
[345]	814	x = time
	815	count = count+1
	816	if indexlon is None:
[350]	817	print "axis is lon"
[345]	818	if count == 0: x = lon
	819	else: y = lon
	820	count = count+1
	821	if indexlat is None:
[350]	822	print "axis is lat"
[345]	823	if count == 0: x = lat
	824	else: y = lat
	825	count = count+1
	826	if indexvert is None and dim0 is 4:
[350]	827	print "axis is vert"
[345]	828	if vertmode == 0: # vertical axis is as is (GCM grid)
	829	if count == 0: x=range(len(vert))
	830	else: y=range(len(vert))
	831	count = count+1
	832	else: # vertical axis is in kms
	833	if count == 0: x = vert
	834	else: y = vert
	835	count = count+1
	836	x = array(x)
	837	y = array(y)
[350]	838	print "what_I_plot.shape", what_I_plot.shape
	839	print "x.shape, y.shape", x.shape, y.shape
[345]	840	if len(shape) == 1:
[350]	841	print shape[0]
	842	if shape[0] != len(x):
[345]	843	print "WARNING HERE !!!"
	844	x = y
	845	elif len(shape) == 2:
	846	print shape[1], len(y), shape[0], len(x)
	847	if shape[1] == len(y) and shape[0] == len(x) and shape[0] != shape[1]:
	848	what_I_plot = swapaxes(what_I_plot,0,1)
	849	print "swapaxes", what_I_plot.shape, shape
	850	#x0 = x
	851	#x = y
	852	#y = x0
	853	#print "define_axis", x, y
	854	return what_I_plot,x,y
[349]	855
	856	# Author: TN + AS
	857	def determineplot(slon, slat, svert, stime):
	858	nlon = 1 # number of longitudinal slices -- 1 is None
	859	nlat = 1
	860	nvert = 1
	861	ntime = 1
	862	nslices = 1
	863	if slon is not None:
	864	nslices = nslices*len(slon)
	865	nlon = len(slon)
	866	if slat is not None:
	867	nslices = nslices*len(slat)
	868	nlat = len(slat)
	869	if svert is not None:
	870	nslices = nslices*len(svert)
	871	nvert = len(svert)
	872	if stime is not None:
	873	nslices = nslices*len(stime)
	874	ntime = len(stime)
	875	#else:
	876	# nslices = 2
	877
	878	mapmode = 0
	879	if slon is None and slat is None:
	880	mapmode = 1 # in this case we plot a map, with the given projection
	881	#elif proj is not None:
	882	# print "WARNING: you specified a", proj,\
	883	# "projection but asked for slices", slon,"in longitude and",slat,"in latitude"
	884	print "mapmode: ", mapmode
	885
	886	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 @ 376

Download in other formats: