Context Navigation

source: trunk/UTIL/PYTHON/planetoplot_v2/ppplot.py @ 978

Last change on this file since 978 was 977, checked in by aslmd, 12 years ago
UTIL PYTHON planetoplot_v2. Added rainbow spahetti plot based on an idea by Alizee. Simply set a colormap for a 1D plot and you get it. Also added in pp_reload (simple rainbow function in ppplot). Added also xmin and xmax capabilitiy. Fixed various bugs related to dummy axis, logarithmic axes. Removed time_counter from set_ppclass.txt it is useless in most cases.
File size: 31.3 KB

Rev	Line
[910]	1	###############################################
	2	## PLANETOPLOT ##
	3	## --> PPPLOT ##
	4	###############################################
	5	## Author: Aymeric Spiga. 02-03/2013 ##
	6	###############################################
	7	# python built-in librairies
	8	import time as timelib
	9	# added librairies
	10	import numpy as np
	11	import matplotlib.pyplot as mpl
	12	from matplotlib.cm import get_cmap
	13	from mpl_toolkits.basemap import Basemap
[936]	14	from matplotlib.ticker import FormatStrFormatter,MaxNLocator
[910]	15	# personal librairies
	16	import ppcompute
	17	###############################################
	18
	19	#################################
	20	# global variables and settings #
	21	#################################
	22
	23	# matplotlib settings
	24	# http://matplotlib.org/users/customizing.html
	25	# -------------------------------
	26	mpl.rcParams['font.family'] = "serif"
[947]	27	mpl.rcParams['axes.color_cycle'] = "b,r,g,k"
[928]	28	mpl.rcParams['contour.negative_linestyle'] = "dashed" # ou "solid"
[910]	29	mpl.rcParams['verbose.level'] = "silent"
	30	mpl.rcParams['lines.linewidth'] = 1.5
	31	mpl.rcParams['lines.markersize'] = 10
	32	mpl.rcParams['xtick.major.pad'] = 10
	33	mpl.rcParams['ytick.major.pad'] = 10
	34
	35	# global variables
	36	# -------------------------------
	37	# - where settings files are located
	38	# (None means planetoplot_v2 in PYTHONPATH)
	39	whereset = None
	40	# - some good default settings.
	41	# (bounds)
	42	how_many_sigma = 3.0
	43	# (contours)
	44	ccol = 'black'
	45	cline = 0.55
	46	# (vectors)
	47	widthvec = 0.002
	48	reducevec = 30.
	49	# (colorbar)
	50	zeorientation="vertical"
	51	zefrac = 0.05
	52	# (save figures)
	53	pad_inches_value=0.25
[928]	54	# (negative mode)
	55	def_negative = False
[910]	56	###############################################
	57
[928]	58	### settings for 'negative-like' mode
	59	if def_negative:
	60	mpl.rc('figure', facecolor='k', edgecolor='k')
	61	mpl.rcParams['text.color'] = 'w'
	62	mpl.rc('axes',labelcolor='w',facecolor='k',edgecolor='w')
	63	mpl.rcParams['xtick.color'] = 'w'
	64	mpl.rcParams['ytick.color'] = 'w'
	65	mpl.rcParams['grid.color'] = 'w'
	66	mpl.rc('savefig',facecolor='k',edgecolor='k')
	67
[910]	68	##########################
	69	# executed when imported #
	70	##########################
	71	###########################################
	72	# we load user-defined automatic settings #
	73	###########################################
	74	# initialize the warning variable about file not present...
	75	files_not_present = ""
	76	# ... and the whereset variable
	77	whereset = ppcompute.findset(whereset)
	78
	79	# - variable settings
	80	# -------------------------------
	81	zefile = "set_var.txt"
[933]	82	vf = {} ; vc = {} ; vl = {} ; vu = {}
[910]	83	try:
	84	f = open(whereset+zefile, 'r')
	85	for line in f:
	86	if "#" in line: pass
	87	else:
[933]	88	var, format, colorb, label, units = line.strip().split(';')
[910]	89	ind = var.strip()
	90	vf[ind] = format.strip()
	91	vc[ind] = colorb.strip()
	92	vl[ind] = label.strip()
[933]	93	vu[ind] = units.strip()
[910]	94	f.close()
	95	except IOError:
	96	files_not_present = files_not_present + zefile + " "
	97
	98	## - file settings
	99	## -------------------------------
	100	#zefile = "set_file.txt"
	101	#prefix_t = {} ; suffix_t = {} ; name_t = {} ; proj_t = {} ; vecx_t = {} ; vecy_t = {}
	102	#try:
	103	# f = open(whereset+zefile, 'r')
	104	# for line in f:
	105	# if "#" in line: pass
	106	# else:
	107	# prefix, suffix, name, proj, vecx, vecy = line.strip().split(';')
	108	# ind = name.strip()
	109	# prefix_t[ind] = prefix.strip()
	110	# suffix_t[ind] = suffix.strip()
	111	# #name_t[ind] = name.strip()
	112	# proj_t[ind] = proj.strip()
	113	# vecx_t[ind] = vecx.strip()
	114	# vecy_t[ind] = vecy.strip()
	115	# f.close()
	116	#except IOError:
	117	# files_not_present = files_not_present + zefile + " "
	118
	119	# - multiplot settings
	120	# -------------------------------
	121	zefile = "set_multiplot.txt"
	122	subv_t = {} ; subh_t = {} ; wspace_t = {} ; hspace_t = {} ; font_t = {}
	123	try:
	124	f = open(whereset+zefile, 'r')
	125	for line in f:
	126	if "#" in line: pass
	127	else:
	128	num, subv, subh, wspace, hspace, font = line.strip().split(';')
	129	ind = int(num.strip())
	130	subv_t[ind] = int(subv.strip())
	131	subh_t[ind] = int(subh.strip())
	132	wspace_t[ind] = float(wspace.strip())
	133	hspace_t[ind] = float(hspace.strip())
	134	font_t[ind] = float(font.strip())
	135	f.close()
	136	except IOError:
	137	files_not_present = files_not_present + zefile + " "
	138
	139	# - background settings
	140	# -------------------------------
	141	zefile = "set_back.txt"
	142	back = {}
	143	try:
	144	f = open(whereset+zefile, 'r')
	145	for line in f:
	146	if "#" in line: pass
	147	else:
	148	name, link = line.strip().split(';')
	149	ind = name.strip()
	150	back[ind] = link.strip()
	151	f.close()
	152	except IOError:
	153	files_not_present = files_not_present + zefile + " "
	154
	155	# - area settings
	156	# -------------------------------
	157	zefile = "set_area.txt"
	158	area = {}
	159	try:
	160	f = open(whereset+zefile, 'r')
	161	for line in f:
	162	if "#" in line: pass
	163	else:
	164	name, wlat1, wlat2, wlon1, wlon2 = line.strip().split(';')
	165	area[name.strip()] = [[float(wlon1.strip()),float(wlon2.strip())],\
	166	[float(wlat1.strip()),float(wlat2.strip())]]
	167	f.close()
	168	except IOError:
	169	files_not_present = files_not_present + zefile + " "
	170	# A note to the user about missing files
	171	if files_not_present != "":
	172	print "warning: files "+files_not_present+" not in "+whereset+" ; those presets will be missing"
	173
	174	# TBD: should change vector color with colormaps
	175	#"gist_heat": "white",\
	176	#"hot": "white",\
	177	#"gray": "red",\
	178
	179	####################
	180	# useful functions #
	181	####################
	182
[977]	183	# a function to change color lines according to a color map (idea by A. Pottier)
	184	# ... two optional arguments: color maps + a number telling how much intervals are needed
	185	def rainbow(cb="jet",num=8):
	186	ax = mpl.gca()
	187	pal = mpl.cm.get_cmap(name=cb)
	188	ax._get_lines.set_color_cycle([pal(i) for i in np.linspace(0,0.9,num)])
	189
[910]	190	# a function to define subplot
	191	# ... user can change settings in set_multiplot.txt read above
	192	# -------------------------------
	193	def definesubplot(numplot, fig):
	194	try:
	195	mpl.rcParams['font.size'] = font_t[numplot]
	196	except:
	197	mpl.rcParams['font.size'] = 18
	198	try:
	199	fig.subplots_adjust(wspace = wspace_t[numplot], hspace = hspace_t[numplot])
	200	subv, subh = subv_t[numplot],subh_t[numplot]
	201	except:
	202	print "!! WARNING !! no preset found from set_multiplot.txt, or this setting file was not found."
	203	subv = 1 ; subh = numplot
	204	return subv,subh
	205
	206	# a function to calculate automatically bounds (or simply prescribe those)
	207	# -------------------------------
	208	def calculate_bounds(field,vmin=None,vmax=None):
	209	# prescribed cases first
	210	zevmin = vmin
	211	zevmax = vmax
	212	# computed cases
	213	if zevmin is None or zevmax is None:
	214	# select values
	215	ind = np.where(field < 9e+35)
	216	fieldcalc = field[ ind ] # field must be a numpy array
	217	# calculate stdev and mean
	218	dev = np.std(fieldcalc)*how_many_sigma
	219	damean = ppcompute.mean(fieldcalc)
	220	# fill min/max if needed
	221	if vmin is None: zevmin = damean - dev
	222	if vmax is None: zevmax = damean + dev
	223	# special case: negative values with stddev while field is positive
[933]	224	if zevmin < 0. and ppcompute.min(fieldcalc) >= 0.: zevmin = 0.
[920]	225	# check that bounds are not too tight given the field
	226	amin = ppcompute.min(field)
	227	amax = ppcompute.max(field)
[936]	228	if np.abs(amin) < 1.e-15:
	229	cmin = 0.
	230	else:
	231	cmin = 100.*np.abs((amin - zevmin)/amin)
[920]	232	cmax = 100.*np.abs((amax - zevmax)/amax)
	233	if cmin > 150. or cmax > 150.:
	234	print "!! WARNING !! Bounds are a bit too tight. Might need to reconsider those."
	235	print "!! WARNING !! --> actual",amin,amax,"adopted",zevmin,zevmax
[910]	236	return zevmin, zevmax
	237	#### treat vmin = vmax for continuity
	238	#if vmin == vmax: zevmin = damean - dev ; zevmax = damean + dev
	239
	240	# a function to solve the problem with blank bounds !
	241	# -------------------------------
	242	def bounds(what_I_plot,zevmin,zevmax,miss=9e+35):
	243	small_enough = 1.e-7
	244	if zevmin < 0: what_I_plot[ what_I_plot < zevmin(1.-small_enough) ] = zevmin(1.-small_enough)
	245	else: what_I_plot[ what_I_plot < zevmin(1.+small_enough) ] = zevmin(1.+small_enough)
	246	what_I_plot[ what_I_plot > miss ] = -miss
	247	what_I_plot[ what_I_plot > zevmax ] = zevmax*(1.-small_enough)
	248	return what_I_plot
	249
	250	# a generic function to show (GUI) or save a plot (PNG,EPS,PDF,...)
	251	# -------------------------------
[917]	252	def save(mode="gui",filename="plot",folder="./",includedate=True,res=150,custom=False):
[933]	253	if mode != "nothing":
	254	# a few settings
	255	possiblesave = ['eps','ps','svg','png','jpg','pdf']
	256	# now the main instructions
	257	if mode == "gui":
	258	mpl.show()
	259	elif mode in possiblesave:
	260	## name of plot
	261	name = folder+'/'+filename
	262	if includedate:
	263	for ttt in timelib.gmtime():
	264	name = name + "_" + str(ttt)
	265	name = name +"."+mode
	266	## save file
	267	print "**** Saving file in "+mode+" format... Please wait."
	268	if not custom:
	269	# ... regular plots
	270	mpl.savefig(name,dpi=res,pad_inches=pad_inches_value,bbox_inches='tight')
	271	else:
	272	# ... mapping mode, adapted space for labels etc...
	273	mpl.savefig(name,dpi=res)
	274	else:
	275	print "!! ERROR !! File format not supported. Supported: ",possiblesave
[910]	276
	277	##################################
	278	# a generic class to make a plot #
	279	##################################
	280	class plot():
	281
	282	# print out a help string when help is invoked on the object
	283	# -------------------------------
	284	def __repr__(self):
	285	whatprint = 'plot object. \"help(plot)\" for more information\n'
	286	return whatprint
	287
	288	# default settings
	289	# -- user can define settings by two methods.
	290	# -- 1. yeah = plot2d(title="foo")
	291	# -- 2. yeah = pp() ; yeah.title = "foo"
	292	# -------------------------------
	293	def __init__(self,\
	294	var=None,\
	295	field=None,\
	296	absc=None,\
	297	xlabel="",\
	298	ylabel="",\
	299	div=20,\
	300	logx=False,\
	301	logy=False,\
	302	swap=False,\
	303	swaplab=True,\
	304	invert=False,\
	305	xcoeff=1.,\
	306	ycoeff=1.,\
[933]	307	fmt=None,\
[977]	308	colorb="jet",\
[933]	309	units="",\
[910]	310	title=""):
	311	## what could be defined by the user
	312	self.var = var
	313	self.field = field
	314	self.absc = absc
	315	self.xlabel = xlabel
	316	self.ylabel = ylabel
	317	self.title = title
	318	self.div = div
	319	self.logx = logx
	320	self.logy = logy
	321	self.swap = swap
	322	self.swaplab = swaplab # NB: swaplab only used if swap=True
	323	self.invert = invert
	324	self.xcoeff = xcoeff
	325	self.ycoeff = ycoeff
[933]	326	self.fmt = fmt
	327	self.units = units
[977]	328	self.colorb = colorb
[910]	329	## other useful arguments
	330	## ... not used here in ppplot but need to be attached to plot object
	331	self.axisbg = "white"
	332	self.superpose = False
	333
	334	# check
	335	# -------------------------------
	336	def check(self):
	337	if self.field is None: print "!! ERROR !! Please define a field to be plotted" ; exit()
	338
	339	# define_from_var
	340	# ... this uses settings in set_var.txt
	341	# -------------------------------
	342	def define_from_var(self):
	343	if self.var is not None:
	344	if self.var.upper() in vl.keys():
	345	self.title = vl[self.var.upper()]
[933]	346	self.units = vu[self.var.upper()]
[910]	347
	348	# make
	349	# this is generic to all plots
	350	# -------------------------------
	351	def make(self):
	352	self.check()
[928]	353	# labels, title, etc...
[910]	354	mpl.xlabel(self.xlabel)
	355	mpl.ylabel(self.ylabel)
	356	if self.swap:
	357	if self.swaplab:
	358	mpl.xlabel(self.ylabel)
	359	mpl.ylabel(self.xlabel)
	360	mpl.title(self.title)
	361	# if masked array, set masked values to filled values (e.g. np.nan) for plotting purposes
	362	if type(self.field).__name__ in 'MaskedArray':
	363	self.field[self.field.mask] = self.field.fill_value
	364
	365	################################
	366	# a subclass to make a 1D plot #
	367	################################
	368	class plot1d(plot):
	369
	370	# print out a help string when help is invoked on the object
	371	# -------------------------------
	372	def __repr__(self):
	373	whatprint = 'plot1d object. \"help(plot1d)\" for more information\n'
	374	return whatprint
	375
	376	# default settings
	377	# -- user can define settings by two methods.
	378	# -- 1. yeah = plot1d(title="foo")
	379	# -- 2. yeah = pp() ; yeah.title = "foo"
	380	# -------------------------------
	381	def __init__(self,\
[947]	382	lstyle=None,\
	383	color=None,\
[923]	384	marker='x',\
	385	label=None):
[910]	386	## get initialization from parent class
	387	plot.__init__(self)
	388	## what could be defined by the user
	389	self.lstyle = lstyle
	390	self.color = color
	391	self.marker = marker
[923]	392	self.label = label
[977]	393	self.xmin = None
	394	self.ymin = None
	395	self.xmax = None
	396	self.ymax = None
[910]	397
	398	# define_from_var
	399	# ... this uses settings in set_var.txt
	400	# -------------------------------
	401	def define_from_var(self):
	402	# get what is done in the parent class
	403	plot.define_from_var(self)
	404	# add specific stuff
	405	if self.var is not None:
	406	if self.var.upper() in vl.keys():
[943]	407	self.ylabel = vl[self.var.upper()] + " (" + vu[self.var.upper()] + ")"
[910]	408	self.title = ""
[933]	409	self.fmt = vf[self.var.upper()]
[910]	410
	411	# make
	412	# -------------------------------
	413	def make(self):
	414	# get what is done in the parent class
	415	plot.make(self)
[933]	416	if self.fmt is None: self.fmt = '%.0f'
[910]	417	# add specific stuff
	418	mpl.grid(color='grey')
	419	if self.lstyle == "": self.lstyle = " " # to allow for no line at all with ""
[977]	420	# set dummy x axis if not defined
	421	if self.absc is None:
	422	self.absc = np.array(range(self.field.shape[0]))
	423	print "!! WARNING !! dummy coordinates on x axis"
[910]	424	# swapping if requested
	425	if self.swap: x = self.field ; y = self.absc
	426	else: x = self.absc ; y = self.field
	427	# coefficients on axis
	428	x=xself.xcoeff ; y=yself.ycoeff
	429	# check axis
	430	if x.size != y.size:
	431	print "!! ERROR !! x and y sizes don't match on 1D plot.", x.size, y.size
	432	exit()
	433	# make the 1D plot
	434	# either request linestyle or let matplotlib decide
	435	if self.lstyle is not None and self.color is not None:
[923]	436	mpl.plot(x,y,self.color+self.lstyle,marker=self.marker,label=self.label)
[942]	437	elif self.color is not None:
	438	mpl.plot(x,y,color=self.color,marker=self.marker,label=self.label)
	439	elif self.lstyle is not None:
[947]	440	mpl.plot(x,y,linestyle=self.lstyle,marker=self.marker,label=self.label)
[910]	441	else:
[923]	442	mpl.plot(x,y,marker=self.marker,label=self.label)
[910]	443	# make log axes and/or invert ordinate
	444	# ... this must be after plot so that axis bounds are well-defined
	445	# ... also inverting must be after making the thing logarithmic
[977]	446	if self.logx: mpl.xscale("log") # not mpl.semilogx() because excludes log on y
	447	if self.logy: mpl.yscale("log") # not mpl.semilogy() because excludes log on x
[910]	448	if self.invert: ax = mpl.gca() ; ax.set_ylim(ax.get_ylim()[::-1])
[923]	449	# add a label for line(s)
	450	if self.label is not None:
	451	if self.label != "":
	452	mpl.legend(loc="best",fancybox=True)
[933]	453	# AXES
	454	ax = mpl.gca()
	455	# format labels
[936]	456	if self.swap: ax.xaxis.set_major_formatter(FormatStrFormatter(self.fmt))
	457	else: ax.yaxis.set_major_formatter(FormatStrFormatter(self.fmt))
[933]	458	# plot limits: ensure that no curve would be outside the window
	459	# x-axis
	460	x1, x2 = ax.get_xbound()
	461	xmin = ppcompute.min(x)
	462	xmax = ppcompute.max(x)
	463	if xmin < x1: x1 = xmin
	464	if xmax > x2: x2 = xmax
[977]	465	if self.xmin is not None: x1 = self.xmin
	466	if self.xmax is not None: x2 = self.xmax
[933]	467	ax.set_xbound(lower=x1,upper=x2)
	468	# y-axis
	469	y1, y2 = ax.get_ybound()
	470	ymin = ppcompute.min(y)
	471	ymax = ppcompute.max(y)
	472	if ymin < y1: y1 = ymin
	473	if ymax > y2: y2 = ymax
[977]	474	if self.ymin is not None: y1 = self.ymin
	475	if self.ymax is not None: y2 = self.ymax
[933]	476	ax.set_ybound(lower=y1,upper=y2)
[936]	477	## set with .div the number of ticks. (is it better than automatic?)
	478	ax.xaxis.set_major_locator(MaxNLocator(self.div/2))
[910]	479
	480	################################
	481	# a subclass to make a 2D plot #
	482	################################
	483	class plot2d(plot):
	484
	485	# print out a help string when help is invoked on the object
	486	# -------------------------------
	487	def __repr__(self):
	488	whatprint = 'plot2d object. \"help(plot2d)\" for more information\n'
	489	return whatprint
	490
	491	# default settings
	492	# -- user can define settings by two methods.
	493	# -- 1. yeah = plot2d(title="foo")
	494	# -- 2. yeah = pp() ; yeah.title = "foo"
	495	# -------------------------------
	496	def __init__(self,\
	497	ordi=None,\
	498	mapmode=False,\
	499	proj="cyl",\
	500	back=None,\
	501	trans=1.0,\
	502	addvecx=None,\
	503	addvecy=None,\
	504	addcontour=None,\
	505	blon=None,\
	506	blat=None,\
	507	area=None,\
	508	vmin=None,\
	509	vmax=None,\
[964]	510	showcb=True,\
[972]	511	wscale=None,\
[910]	512	colorvec="black"):
	513	## get initialization from parent class
	514	plot.__init__(self)
	515	## what could be defined by the user
	516	self.ordi = ordi
	517	self.mapmode = mapmode
	518	self.proj = proj
	519	self.back = back
	520	self.trans = trans
	521	self.addvecx = addvecx
	522	self.addvecy = addvecy
	523	self.colorvec = colorvec
	524	self.addcontour = addcontour
	525	self.blon = blon ; self.blat = blat
	526	self.area = area
	527	self.vmin = vmin ; self.vmax = vmax
[964]	528	self.showcb = showcb
[972]	529	self.wscale = wscale
[910]	530
	531	# define_from_var
	532	# ... this uses settings in set_var.txt
	533	# -------------------------------
	534	def define_from_var(self):
	535	# get what is done in the parent class
	536	plot.define_from_var(self)
	537	# add specific stuff
	538	if self.var is not None:
	539	if self.var.upper() in vl.keys():
	540	self.colorb = vc[self.var.upper()]
	541	self.fmt = vf[self.var.upper()]
	542
	543	# make
	544	# -------------------------------
	545	def make(self):
	546	# get what is done in the parent class...
	547	plot.make(self)
[977]	548	if self.fmt is None: self.fmt = "%.1e"
[910]	549	# ... then add specific stuff
	550	############################################################################################
	551	### PRE-SETTINGS
	552	############################################################################################
[960]	553	# set dummy xy axis if not defined
	554	if self.absc is None:
	555	self.absc = np.array(range(self.field.shape[0]))
	556	self.mapmode = False
	557	print "!! WARNING !! dummy coordinates on x axis"
	558	if self.ordi is None:
	559	self.ordi = np.array(range(self.field.shape[1]))
	560	self.mapmode = False
	561	print "!! WARNING !! dummy coordinates on y axis"
[977]	562	# transposing if necessary
	563	shape = self.field.shape
	564	if shape[0] != shape[1]:
	565	if len(self.absc) == shape[0] and len(self.ordi) == shape[1]:
	566	print "!! WARNING !! Transposing axes"
	567	self.field = np.transpose(self.field)
[910]	568	# bound field
	569	zevmin, zevmax = calculate_bounds(self.field,vmin=self.vmin,vmax=self.vmax)
	570	what_I_plot = bounds(self.field,zevmin,zevmax)
	571	# define contour field levels. define color palette
	572	ticks = self.div + 1
	573	zelevels = np.linspace(zevmin,zevmax,ticks)
	574	palette = get_cmap(name=self.colorb)
	575	# do the same thing for possible contourline entries
	576	if self.addcontour is not None:
	577	# if masked array, set masked values to filled values (e.g. np.nan) for plotting purposes
	578	if type(self.addcontour).__name__ in 'MaskedArray':
	579	self.addcontour[self.addcontour.mask] = self.addcontour.fill_value
	580	zevminc, zevmaxc = calculate_bounds(self.addcontour)
	581	what_I_contour = bounds(self.addcontour,zevminc,zevmaxc)
	582	ticks = self.div + 1
	583	zelevelsc = np.linspace(zevminc,zevmaxc,ticks)
	584	############################################################################################
	585	### MAIN PLOT
	586	### NB: contour lines are done before contour shades otherwise colorar error
	587	############################################################################################
	588	if not self.mapmode:
	589	## A SIMPLE 2D PLOT
	590	###################
	591	# swapping if requested
	592	if self.swap: x = self.ordi ; y = self.absc
	593	else: x = self.absc ; y = self.ordi
	594	# coefficients on axis
	595	x=xself.xcoeff ; y=yself.ycoeff
	596	# make shaded and line contours
	597	if self.addcontour is not None:
[943]	598	objC = mpl.contour(x, y, what_I_contour, \
[910]	599	zelevelsc, colors = ccol, linewidths = cline)
[943]	600	#mpl.clabel(objC, inline=1, fontsize=10)
[910]	601	mpl.contourf(x, y, \
	602	self.field, \
	603	zelevels, cmap=palette)
	604	# make log axes and/or invert ordinate
	605	if self.logx: mpl.semilogx()
	606	if self.logy: mpl.semilogy()
	607	if self.invert: ax = mpl.gca() ; ax.set_ylim(ax.get_ylim()[::-1])
	608	else:
	609	## A 2D MAP USING PROJECTIONS (basemap)
	610	#######################################
	611	mpl.xlabel("") ; mpl.ylabel("")
	612	# additional security in case self.proj is None here
	613	# ... we set cylindrical projection (the simplest one)
	614	if self.proj is None: self.proj = "cyl"
	615	# get lon and lat in 2D version.
	616	# (but first ensure we do have 2D coordinates)
	617	if self.absc.ndim == 1: [self.absc,self.ordi] = np.meshgrid(self.absc,self.ordi)
	618	elif self.absc.ndim > 2: print "!! ERROR !! lon and lat arrays must be 1D or 2D"
	619	# get lat lon intervals and associated settings
	620	wlon = [np.min(self.absc),np.max(self.absc)]
	621	wlat = [np.min(self.ordi),np.max(self.ordi)]
	622	# -- area presets are in set_area.txt
	623	if self.area is not None:
	624	if self.area in area.keys():
	625	wlon, wlat = area[self.area]
	626	# -- settings for meridians and parallels
[936]	627	steplon = int(abs(wlon[1]-wlon[0])/6.)
	628	steplat = int(abs(wlat[1]-wlat[0])/3.)
	629	#mertab = np.r_[wlon[0]:wlon[1]:steplon] ; merlab = [0,0,0,1]
	630	#partab = np.r_[wlat[0]:wlat[1]:steplat] ; parlab = [1,0,0,0]
[972]	631	mertab = np.r_[-180.:180.:steplon] ; merlab = [0,0,0,1] #-360:360.
[936]	632	partab = np.r_[-90.:90.:steplat] ; parlab = [1,0,0,0]
[910]	633	format = '%.1f'
	634	# -- center of domain and bounding lats
	635	lon_0 = 0.5*(wlon[0]+wlon[1])
	636	lat_0 = 0.5*(wlat[0]+wlat[1])
	637	# some tests, bug fixes, and good-looking settings
	638	# ... cyl is good for global and regional
	639	if self.proj == "cyl":
[917]	640	format = '%.0f'
[910]	641	# ... global projections
	642	elif self.proj in ["ortho","moll","robin"]:
	643	wlat[0] = None ; wlat[1] = None ; wlon[0] = None ; wlon[1] = None
	644	steplon = 30. ; steplat = 30.
[923]	645	if self.proj in ["robin","moll"]: steplon = 60.
[910]	646	mertab = np.r_[-360.:360.:steplon]
	647	partab = np.r_[-90.:90.:steplat]
	648	if self.proj == "ortho":
	649	merlab = [0,0,0,0] ; parlab = [0,0,0,0]
	650	# in ortho projection, blon and blat can be used to set map center
	651	if self.blon is not None: lon_0 = self.blon
	652	if self.blat is not None: lat_0 = self.blat
	653	elif self.proj == "moll":
	654	merlab = [0,0,0,0]
	655	format = '%.0f'
	656	# ... regional projections
	657	elif self.proj in ["lcc","laea","merc"]:
	658	if self.proj in ["lcc","laea"] and wlat[0] == -wlat[1]:
	659	print "!! ERROR !! with Lambert lat1 must be different than lat2" ; exit()
	660	if wlat[0] < -80. and wlat[1] > 80.:
	661	print "!! ERROR !! set an area (not global)" ; exit()
	662	format = '%.0f'
	663	elif self.proj in ["npstere","spstere"]:
	664	# in polar projections, blat gives the bounding lat
	665	# if not set, set something reasonable
	666	if self.blat is None: self.blat = 60.
	667	# help the user who forgets self.blat would better be negative in spstere
	668	# (this actually serves for the default setting just above)
	669	if self.proj == "spstere" and self.blat > 0: self.blat = -self.blat
[932]	670	# labels
	671	mertab = np.r_[-360.:360.:15.]
	672	partab = np.r_[-90.:90.:5.]
[910]	673	# ... unsupported projections
	674	else:
	675	print "!! ERROR !! unsupported projection. supported: "+\
	676	"cyl, npstere, spstere, ortho, moll, robin, lcc, laea, merc"
	677	# finally define projection
	678	m = Basemap(projection=self.proj,\
	679	lat_0=lat_0,lon_0=lon_0,\
	680	boundinglat=self.blat,\
	681	llcrnrlat=wlat[0],urcrnrlat=wlat[1],\
	682	llcrnrlon=wlon[0],urcrnrlon=wlon[1])
	683	# draw meridians and parallels
	684	ft = int(mpl.rcParams['font.size']*3./4.)
[932]	685	zelatmax = 85.
	686	m.drawmeridians(mertab,labels=merlab,color='grey',linewidth=0.75,fontsize=ft,fmt=format,latmax=zelatmax)
	687	m.drawparallels(partab,labels=parlab,color='grey',linewidth=0.75,fontsize=ft,fmt=format,latmax=zelatmax)
[910]	688	# define background (see set_back.txt)
	689	if self.back is not None:
	690	if self.back in back.keys():
	691	print "**** info: loading a background, please wait.",self.back
[933]	692	if self.back not in ["coast","sea"]: m.warpimage(back[self.back],scale=0.75)
	693	elif self.back == "coast": m.drawcoastlines()
	694	elif self.back == "sea": m.drawlsmask(land_color='white',ocean_color='aqua')
[910]	695	else:
	696	print "!! ERROR !! requested background not defined. change name or fill in set_back.txt" ; exit()
	697	# define x and y given the projection
	698	x, y = m(self.absc, self.ordi)
	699	# contour field. first line contour then shaded contour.
	700	if self.addcontour is not None:
[943]	701	objC2 = m.contour(x, y, what_I_contour, \
[910]	702	zelevelsc, colors = ccol, linewidths = cline)
[943]	703	#mpl.clabel(objC2, inline=1, fontsize=10)
[910]	704	m.contourf(x, y, what_I_plot, zelevels, cmap = palette, alpha = self.trans)
	705	############################################################################################
	706	### COLORBAR
	707	############################################################################################
[964]	708	if self.trans > 0. and self.showcb:
[910]	709	## draw colorbar. settings are different with projections. or if not mapmode.
[933]	710	if not self.mapmode: orientation=zeorientation ; frac = 0.075 ; pad = 0.03 ; lu = 0.5
[942]	711	elif self.proj in ['moll']: orientation="horizontal" ; frac = 0.08 ; pad = 0.03 ; lu = 1.0
[933]	712	elif self.proj in ['cyl']: orientation="vertical" ; frac = 0.023 ; pad = 0.03 ; lu = 0.5
	713	else: orientation = zeorientation ; frac = zefrac ; pad = 0.03 ; lu = 0.5
[910]	714	zelevpal = np.linspace(zevmin,zevmax,num=min([ticks/2+1,21]))
[917]	715	zecb = mpl.colorbar(fraction=frac,pad=pad,\
[910]	716	format=self.fmt,orientation=orientation,\
	717	ticks=zelevpal,\
[933]	718	extend='neither',spacing='proportional')
[910]	719	if zeorientation == "horizontal": zecb.ax.set_xlabel(self.title) ; self.title = ""
[933]	720	# colorbar title --> units
	721	if self.units not in ["dimless",""]:
	722	zecb.ax.set_title("["+self.units+"]",fontsize=3.*mpl.rcParams['font.size']/4.,x=lu,y=1.025)
	723
[910]	724	############################################################################################
	725	### VECTORS. must be after the colorbar. we could also leave possibility for streamlines.
	726	############################################################################################
	727	### not expecting NaN in self.addvecx and self.addvecy. masked arrays is just enough.
	728	stride = 3
[942]	729	#stride = 1
[910]	730	if self.addvecx is not None and self.addvecy is not None and self.mapmode:
[932]	731	### for metwinds only ???
[910]	732	[vecx,vecy] = m.rotate_vector(self.addvecx,self.addvecy,self.absc,self.ordi)
[932]	733	#vecx,vecy = self.addvecx,self.addvecy
[910]	734	# reference vector is scaled
[972]	735	if self.wscale is None: self.wscale = ppcompute.mean(np.sqrt(self.addvecxself.addvecx+self.addvecyself.addvecy))
[910]	736	# make vector field
	737	q = m.quiver( x[::stride,::stride],y[::stride,::stride],\
	738	vecx[::stride,::stride],vecy[::stride,::stride],\
	739	angles='xy',color=self.colorvec,pivot='middle',\
[972]	740	scale=self.wscale*reducevec,width=widthvec )
[930]	741	# make vector key.
[933]	742	#keyh = 1.025 ; keyv = 1.05 # upper right corner over colorbar
[972]	743	keyh = 0.97 ; keyv = 1.05
[930]	744	#keyh = -0.03 ; keyv = 1.08 # upper left corner
[923]	745	p = mpl.quiverkey(q,keyh,keyv,\
[972]	746	self.wscale,str(int(self.wscale)),\
	747	fontproperties={'size': 'small'},\
[910]	748	color='black',labelpos='S',labelsep = 0.07)
[942]	749	############################################################################################
	750	### TEXT. ANYWHERE. add_text.txt should be present with lines x ; y ; text ; color
	751	############################################################################################
	752	try:
	753	f = open("add_text.txt", 'r')
	754	for line in f:
	755	if "#" in line: pass
	756	else:
	757	userx, usery, usert, userc = line.strip().split(';')
	758	userc = userc.strip()
	759	usert = usert.strip()
	760	userx = float(userx.strip())
	761	usery = float(usery.strip())
	762	if self.mapmode: userx,usery = m(userx,usery)
	763	mpl.text(userx,usery,usert,\
	764	color = userc,\
	765	horizontalalignment='center',\
	766	verticalalignment='center')
	767	f.close()
	768	except IOError:
	769	pass

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

Download in other formats: