Context Navigation

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

Last change on this file since 938 was 936, checked in by aslmd, 12 years ago
UTIL PYTHON planetoplot_v2. fixed bug with minimum close to 0. number of ticks in 1D plots can be set with self.div. better labelling (actual values). better handling of time in case this is Ls.
File size: 28.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"
	27	mpl.rcParams['axes.color_cycle'] = "r,g,b,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
	183	# a function to define subplot
	184	# ... user can change settings in set_multiplot.txt read above
	185	# -------------------------------
	186	def definesubplot(numplot, fig):
	187	try:
	188	mpl.rcParams['font.size'] = font_t[numplot]
	189	except:
	190	mpl.rcParams['font.size'] = 18
	191	try:
	192	fig.subplots_adjust(wspace = wspace_t[numplot], hspace = hspace_t[numplot])
	193	subv, subh = subv_t[numplot],subh_t[numplot]
	194	except:
	195	print "!! WARNING !! no preset found from set_multiplot.txt, or this setting file was not found."
	196	subv = 1 ; subh = numplot
	197	return subv,subh
	198
	199	# a function to calculate automatically bounds (or simply prescribe those)
	200	# -------------------------------
	201	def calculate_bounds(field,vmin=None,vmax=None):
	202	# prescribed cases first
	203	zevmin = vmin
	204	zevmax = vmax
	205	# computed cases
	206	if zevmin is None or zevmax is None:
	207	# select values
	208	ind = np.where(field < 9e+35)
	209	fieldcalc = field[ ind ] # field must be a numpy array
	210	# calculate stdev and mean
	211	dev = np.std(fieldcalc)*how_many_sigma
	212	damean = ppcompute.mean(fieldcalc)
	213	# fill min/max if needed
	214	if vmin is None: zevmin = damean - dev
	215	if vmax is None: zevmax = damean + dev
	216	# special case: negative values with stddev while field is positive
[933]	217	if zevmin < 0. and ppcompute.min(fieldcalc) >= 0.: zevmin = 0.
[920]	218	# check that bounds are not too tight given the field
	219	amin = ppcompute.min(field)
	220	amax = ppcompute.max(field)
[936]	221	if np.abs(amin) < 1.e-15:
	222	zevmin = 0.
	223	cmin = 0.
	224	else:
	225	cmin = 100.*np.abs((amin - zevmin)/amin)
[920]	226	cmax = 100.*np.abs((amax - zevmax)/amax)
	227	if cmin > 150. or cmax > 150.:
	228	print "!! WARNING !! Bounds are a bit too tight. Might need to reconsider those."
	229	print "!! WARNING !! --> actual",amin,amax,"adopted",zevmin,zevmax
[910]	230	return zevmin, zevmax
	231	#### treat vmin = vmax for continuity
	232	#if vmin == vmax: zevmin = damean - dev ; zevmax = damean + dev
	233
	234	# a function to solve the problem with blank bounds !
	235	# -------------------------------
	236	def bounds(what_I_plot,zevmin,zevmax,miss=9e+35):
	237	small_enough = 1.e-7
	238	if zevmin < 0: what_I_plot[ what_I_plot < zevmin(1.-small_enough) ] = zevmin(1.-small_enough)
	239	else: what_I_plot[ what_I_plot < zevmin(1.+small_enough) ] = zevmin(1.+small_enough)
	240	what_I_plot[ what_I_plot > miss ] = -miss
	241	what_I_plot[ what_I_plot > zevmax ] = zevmax*(1.-small_enough)
	242	return what_I_plot
	243
	244	# a generic function to show (GUI) or save a plot (PNG,EPS,PDF,...)
	245	# -------------------------------
[917]	246	def save(mode="gui",filename="plot",folder="./",includedate=True,res=150,custom=False):
[933]	247	if mode != "nothing":
	248	# a few settings
	249	possiblesave = ['eps','ps','svg','png','jpg','pdf']
	250	# now the main instructions
	251	if mode == "gui":
	252	mpl.show()
	253	elif mode in possiblesave:
	254	## name of plot
	255	name = folder+'/'+filename
	256	if includedate:
	257	for ttt in timelib.gmtime():
	258	name = name + "_" + str(ttt)
	259	name = name +"."+mode
	260	## save file
	261	print "**** Saving file in "+mode+" format... Please wait."
	262	if not custom:
	263	# ... regular plots
	264	mpl.savefig(name,dpi=res,pad_inches=pad_inches_value,bbox_inches='tight')
	265	else:
	266	# ... mapping mode, adapted space for labels etc...
	267	mpl.savefig(name,dpi=res)
	268	else:
	269	print "!! ERROR !! File format not supported. Supported: ",possiblesave
[910]	270
	271	##################################
	272	# a generic class to make a plot #
	273	##################################
	274	class plot():
	275
	276	# print out a help string when help is invoked on the object
	277	# -------------------------------
	278	def __repr__(self):
	279	whatprint = 'plot object. \"help(plot)\" for more information\n'
	280	return whatprint
	281
	282	# default settings
	283	# -- user can define settings by two methods.
	284	# -- 1. yeah = plot2d(title="foo")
	285	# -- 2. yeah = pp() ; yeah.title = "foo"
	286	# -------------------------------
	287	def __init__(self,\
	288	var=None,\
	289	field=None,\
	290	absc=None,\
	291	xlabel="",\
	292	ylabel="",\
	293	div=20,\
	294	logx=False,\
	295	logy=False,\
	296	swap=False,\
	297	swaplab=True,\
	298	invert=False,\
	299	xcoeff=1.,\
	300	ycoeff=1.,\
[933]	301	fmt=None,\
	302	units="",\
[910]	303	title=""):
	304	## what could be defined by the user
	305	self.var = var
	306	self.field = field
	307	self.absc = absc
	308	self.xlabel = xlabel
	309	self.ylabel = ylabel
	310	self.title = title
	311	self.div = div
	312	self.logx = logx
	313	self.logy = logy
	314	self.swap = swap
	315	self.swaplab = swaplab # NB: swaplab only used if swap=True
	316	self.invert = invert
	317	self.xcoeff = xcoeff
	318	self.ycoeff = ycoeff
[933]	319	self.fmt = fmt
	320	self.units = units
[910]	321	## other useful arguments
	322	## ... not used here in ppplot but need to be attached to plot object
	323	self.axisbg = "white"
	324	self.superpose = False
	325
	326	# check
	327	# -------------------------------
	328	def check(self):
	329	if self.field is None: print "!! ERROR !! Please define a field to be plotted" ; exit()
	330
	331	# define_from_var
	332	# ... this uses settings in set_var.txt
	333	# -------------------------------
	334	def define_from_var(self):
	335	if self.var is not None:
	336	if self.var.upper() in vl.keys():
	337	self.title = vl[self.var.upper()]
[933]	338	self.units = vu[self.var.upper()]
[910]	339
	340	# make
	341	# this is generic to all plots
	342	# -------------------------------
	343	def make(self):
	344	self.check()
[928]	345	# labels, title, etc...
[910]	346	mpl.xlabel(self.xlabel)
	347	mpl.ylabel(self.ylabel)
	348	if self.swap:
	349	if self.swaplab:
	350	mpl.xlabel(self.ylabel)
	351	mpl.ylabel(self.xlabel)
	352	mpl.title(self.title)
	353	# if masked array, set masked values to filled values (e.g. np.nan) for plotting purposes
	354	if type(self.field).__name__ in 'MaskedArray':
	355	self.field[self.field.mask] = self.field.fill_value
	356
	357	################################
	358	# a subclass to make a 1D plot #
	359	################################
	360	class plot1d(plot):
	361
	362	# print out a help string when help is invoked on the object
	363	# -------------------------------
	364	def __repr__(self):
	365	whatprint = 'plot1d object. \"help(plot1d)\" for more information\n'
	366	return whatprint
	367
	368	# default settings
	369	# -- user can define settings by two methods.
	370	# -- 1. yeah = plot1d(title="foo")
	371	# -- 2. yeah = pp() ; yeah.title = "foo"
	372	# -------------------------------
	373	def __init__(self,\
	374	lstyle='-',\
	375	color='b',\
[923]	376	marker='x',\
	377	label=None):
[910]	378	## get initialization from parent class
	379	plot.__init__(self)
	380	## what could be defined by the user
	381	self.lstyle = lstyle
	382	self.color = color
	383	self.marker = marker
[923]	384	self.label = label
[910]	385
	386	# define_from_var
	387	# ... this uses settings in set_var.txt
	388	# -------------------------------
	389	def define_from_var(self):
	390	# get what is done in the parent class
	391	plot.define_from_var(self)
	392	# add specific stuff
	393	if self.var is not None:
	394	if self.var.upper() in vl.keys():
	395	self.ylabel = vl[self.var.upper()]
	396	self.title = ""
[933]	397	self.fmt = vf[self.var.upper()]
[910]	398
	399	# make
	400	# -------------------------------
	401	def make(self):
	402	# get what is done in the parent class
	403	plot.make(self)
[933]	404	if self.fmt is None: self.fmt = '%.0f'
[910]	405	# add specific stuff
	406	mpl.grid(color='grey')
	407	if self.lstyle == "": self.lstyle = " " # to allow for no line at all with ""
	408	# swapping if requested
	409	if self.swap: x = self.field ; y = self.absc
	410	else: x = self.absc ; y = self.field
	411	# coefficients on axis
	412	x=xself.xcoeff ; y=yself.ycoeff
	413	# check axis
	414	if x.size != y.size:
	415	print "!! ERROR !! x and y sizes don't match on 1D plot.", x.size, y.size
	416	exit()
	417	# make the 1D plot
	418	# either request linestyle or let matplotlib decide
	419	if self.lstyle is not None and self.color is not None:
[923]	420	mpl.plot(x,y,self.color+self.lstyle,marker=self.marker,label=self.label)
[910]	421	else:
[923]	422	mpl.plot(x,y,marker=self.marker,label=self.label)
[910]	423	# make log axes and/or invert ordinate
	424	# ... this must be after plot so that axis bounds are well-defined
	425	# ... also inverting must be after making the thing logarithmic
	426	if self.logx: mpl.semilogx()
	427	if self.logy: mpl.semilogy()
	428	if self.invert: ax = mpl.gca() ; ax.set_ylim(ax.get_ylim()[::-1])
[923]	429	# add a label for line(s)
	430	if self.label is not None:
	431	if self.label != "":
	432	mpl.legend(loc="best",fancybox=True)
[933]	433	# AXES
	434	ax = mpl.gca()
	435	# format labels
[936]	436	if self.swap: ax.xaxis.set_major_formatter(FormatStrFormatter(self.fmt))
	437	else: ax.yaxis.set_major_formatter(FormatStrFormatter(self.fmt))
[933]	438	# plot limits: ensure that no curve would be outside the window
	439	# x-axis
	440	x1, x2 = ax.get_xbound()
	441	xmin = ppcompute.min(x)
	442	xmax = ppcompute.max(x)
	443	if xmin < x1: x1 = xmin
	444	if xmax > x2: x2 = xmax
	445	ax.set_xbound(lower=x1,upper=x2)
	446	# y-axis
	447	y1, y2 = ax.get_ybound()
	448	ymin = ppcompute.min(y)
	449	ymax = ppcompute.max(y)
	450	if ymin < y1: y1 = ymin
	451	if ymax > y2: y2 = ymax
	452	ax.set_ybound(lower=y1,upper=y2)
[936]	453	## set with .div the number of ticks. (is it better than automatic?)
	454	ax.xaxis.set_major_locator(MaxNLocator(self.div/2))
[910]	455
	456	################################
	457	# a subclass to make a 2D plot #
	458	################################
	459	class plot2d(plot):
	460
	461	# print out a help string when help is invoked on the object
	462	# -------------------------------
	463	def __repr__(self):
	464	whatprint = 'plot2d object. \"help(plot2d)\" for more information\n'
	465	return whatprint
	466
	467	# default settings
	468	# -- user can define settings by two methods.
	469	# -- 1. yeah = plot2d(title="foo")
	470	# -- 2. yeah = pp() ; yeah.title = "foo"
	471	# -------------------------------
	472	def __init__(self,\
	473	ordi=None,\
	474	mapmode=False,\
	475	proj="cyl",\
	476	back=None,\
	477	colorb="jet",\
	478	trans=1.0,\
	479	addvecx=None,\
	480	addvecy=None,\
	481	addcontour=None,\
	482	blon=None,\
	483	blat=None,\
	484	area=None,\
	485	vmin=None,\
	486	vmax=None,\
	487	colorvec="black"):
	488	## get initialization from parent class
	489	plot.__init__(self)
	490	## what could be defined by the user
	491	self.ordi = ordi
	492	self.mapmode = mapmode
	493	self.proj = proj
	494	self.back = back
	495	self.colorb = colorb
	496	self.trans = trans
	497	self.addvecx = addvecx
	498	self.addvecy = addvecy
	499	self.colorvec = colorvec
	500	self.addcontour = addcontour
	501	self.blon = blon ; self.blat = blat
	502	self.area = area
	503	self.vmin = vmin ; self.vmax = vmax
	504
	505	# define_from_var
	506	# ... this uses settings in set_var.txt
	507	# -------------------------------
	508	def define_from_var(self):
	509	# get what is done in the parent class
	510	plot.define_from_var(self)
	511	# add specific stuff
	512	if self.var is not None:
	513	if self.var.upper() in vl.keys():
	514	self.colorb = vc[self.var.upper()]
	515	self.fmt = vf[self.var.upper()]
	516
	517	# make
	518	# -------------------------------
	519	def make(self):
	520	# get what is done in the parent class...
	521	plot.make(self)
[933]	522	if self.fmt is None: self.fmt = "%.2e"
[910]	523	# ... then add specific stuff
	524	############################################################################################
	525	### PRE-SETTINGS
	526	############################################################################################
	527	# transposing if necessary
	528	shape = self.field.shape
	529	if shape[0] != shape[1]:
	530	if len(self.absc) == shape[0] and len(self.ordi) == shape[1]:
	531	print "!! WARNING !! Transposing axes"
	532	self.field = np.transpose(self.field)
	533	# bound field
	534	zevmin, zevmax = calculate_bounds(self.field,vmin=self.vmin,vmax=self.vmax)
	535	what_I_plot = bounds(self.field,zevmin,zevmax)
	536	# define contour field levels. define color palette
	537	ticks = self.div + 1
	538	zelevels = np.linspace(zevmin,zevmax,ticks)
	539	palette = get_cmap(name=self.colorb)
	540	# do the same thing for possible contourline entries
	541	if self.addcontour is not None:
	542	# if masked array, set masked values to filled values (e.g. np.nan) for plotting purposes
	543	if type(self.addcontour).__name__ in 'MaskedArray':
	544	self.addcontour[self.addcontour.mask] = self.addcontour.fill_value
	545	zevminc, zevmaxc = calculate_bounds(self.addcontour)
	546	what_I_contour = bounds(self.addcontour,zevminc,zevmaxc)
	547	ticks = self.div + 1
	548	zelevelsc = np.linspace(zevminc,zevmaxc,ticks)
	549	############################################################################################
	550	### MAIN PLOT
	551	### NB: contour lines are done before contour shades otherwise colorar error
	552	############################################################################################
	553	if not self.mapmode:
	554	## A SIMPLE 2D PLOT
	555	###################
	556	# swapping if requested
	557	if self.swap: x = self.ordi ; y = self.absc
	558	else: x = self.absc ; y = self.ordi
	559	# coefficients on axis
	560	x=xself.xcoeff ; y=yself.ycoeff
	561	# make shaded and line contours
	562	if self.addcontour is not None:
	563	mpl.contour(x, y, what_I_contour, \
	564	zelevelsc, colors = ccol, linewidths = cline)
	565	mpl.contourf(x, y, \
	566	self.field, \
	567	zelevels, cmap=palette)
	568	# make log axes and/or invert ordinate
	569	if self.logx: mpl.semilogx()
	570	if self.logy: mpl.semilogy()
	571	if self.invert: ax = mpl.gca() ; ax.set_ylim(ax.get_ylim()[::-1])
	572	else:
	573	## A 2D MAP USING PROJECTIONS (basemap)
	574	#######################################
	575	mpl.xlabel("") ; mpl.ylabel("")
	576	# additional security in case self.proj is None here
	577	# ... we set cylindrical projection (the simplest one)
	578	if self.proj is None: self.proj = "cyl"
	579	# get lon and lat in 2D version.
	580	# (but first ensure we do have 2D coordinates)
	581	if self.absc.ndim == 1: [self.absc,self.ordi] = np.meshgrid(self.absc,self.ordi)
	582	elif self.absc.ndim > 2: print "!! ERROR !! lon and lat arrays must be 1D or 2D"
	583	# get lat lon intervals and associated settings
	584	wlon = [np.min(self.absc),np.max(self.absc)]
	585	wlat = [np.min(self.ordi),np.max(self.ordi)]
	586	# -- area presets are in set_area.txt
	587	if self.area is not None:
	588	if self.area in area.keys():
	589	wlon, wlat = area[self.area]
	590	# -- settings for meridians and parallels
[936]	591	steplon = int(abs(wlon[1]-wlon[0])/6.)
	592	steplat = int(abs(wlat[1]-wlat[0])/3.)
	593	#mertab = np.r_[wlon[0]:wlon[1]:steplon] ; merlab = [0,0,0,1]
	594	#partab = np.r_[wlat[0]:wlat[1]:steplat] ; parlab = [1,0,0,0]
	595	mertab = np.r_[-360.:360.:steplon] ; merlab = [0,0,0,1]
	596	partab = np.r_[-90.:90.:steplat] ; parlab = [1,0,0,0]
[910]	597	format = '%.1f'
	598	# -- center of domain and bounding lats
	599	lon_0 = 0.5*(wlon[0]+wlon[1])
	600	lat_0 = 0.5*(wlat[0]+wlat[1])
	601	# some tests, bug fixes, and good-looking settings
	602	# ... cyl is good for global and regional
	603	if self.proj == "cyl":
[917]	604	format = '%.0f'
[910]	605	# ... global projections
	606	elif self.proj in ["ortho","moll","robin"]:
	607	wlat[0] = None ; wlat[1] = None ; wlon[0] = None ; wlon[1] = None
	608	steplon = 30. ; steplat = 30.
[923]	609	if self.proj in ["robin","moll"]: steplon = 60.
[910]	610	mertab = np.r_[-360.:360.:steplon]
	611	partab = np.r_[-90.:90.:steplat]
	612	if self.proj == "ortho":
	613	merlab = [0,0,0,0] ; parlab = [0,0,0,0]
	614	# in ortho projection, blon and blat can be used to set map center
	615	if self.blon is not None: lon_0 = self.blon
	616	if self.blat is not None: lat_0 = self.blat
	617	elif self.proj == "moll":
	618	merlab = [0,0,0,0]
	619	format = '%.0f'
	620	# ... regional projections
	621	elif self.proj in ["lcc","laea","merc"]:
	622	if self.proj in ["lcc","laea"] and wlat[0] == -wlat[1]:
	623	print "!! ERROR !! with Lambert lat1 must be different than lat2" ; exit()
	624	if wlat[0] < -80. and wlat[1] > 80.:
	625	print "!! ERROR !! set an area (not global)" ; exit()
	626	format = '%.0f'
	627	elif self.proj in ["npstere","spstere"]:
	628	# in polar projections, blat gives the bounding lat
	629	# if not set, set something reasonable
	630	if self.blat is None: self.blat = 60.
	631	# help the user who forgets self.blat would better be negative in spstere
	632	# (this actually serves for the default setting just above)
	633	if self.proj == "spstere" and self.blat > 0: self.blat = -self.blat
[932]	634	# labels
	635	mertab = np.r_[-360.:360.:15.]
	636	partab = np.r_[-90.:90.:5.]
[910]	637	# ... unsupported projections
	638	else:
	639	print "!! ERROR !! unsupported projection. supported: "+\
	640	"cyl, npstere, spstere, ortho, moll, robin, lcc, laea, merc"
	641	# finally define projection
	642	m = Basemap(projection=self.proj,\
	643	lat_0=lat_0,lon_0=lon_0,\
	644	boundinglat=self.blat,\
	645	llcrnrlat=wlat[0],urcrnrlat=wlat[1],\
	646	llcrnrlon=wlon[0],urcrnrlon=wlon[1])
	647	# draw meridians and parallels
	648	ft = int(mpl.rcParams['font.size']*3./4.)
[932]	649	zelatmax = 85.
	650	m.drawmeridians(mertab,labels=merlab,color='grey',linewidth=0.75,fontsize=ft,fmt=format,latmax=zelatmax)
	651	m.drawparallels(partab,labels=parlab,color='grey',linewidth=0.75,fontsize=ft,fmt=format,latmax=zelatmax)
[910]	652	# define background (see set_back.txt)
	653	if self.back is not None:
	654	if self.back in back.keys():
	655	print "**** info: loading a background, please wait.",self.back
[933]	656	if self.back not in ["coast","sea"]: m.warpimage(back[self.back],scale=0.75)
	657	elif self.back == "coast": m.drawcoastlines()
	658	elif self.back == "sea": m.drawlsmask(land_color='white',ocean_color='aqua')
[910]	659	else:
	660	print "!! ERROR !! requested background not defined. change name or fill in set_back.txt" ; exit()
	661	# define x and y given the projection
	662	x, y = m(self.absc, self.ordi)
	663	# contour field. first line contour then shaded contour.
	664	if self.addcontour is not None:
	665	m.contour(x, y, what_I_contour, \
	666	zelevelsc, colors = ccol, linewidths = cline)
	667	m.contourf(x, y, what_I_plot, zelevels, cmap = palette, alpha = self.trans)
	668	############################################################################################
	669	### COLORBAR
	670	############################################################################################
	671	if self.trans > 0.:
	672	## draw colorbar. settings are different with projections. or if not mapmode.
[933]	673	if not self.mapmode: orientation=zeorientation ; frac = 0.075 ; pad = 0.03 ; lu = 0.5
	674	elif self.proj in ['moll']: orientation="horizontal" ; frac = 0.075 ; pad = 0.03 ; lu = 1.0
	675	elif self.proj in ['cyl']: orientation="vertical" ; frac = 0.023 ; pad = 0.03 ; lu = 0.5
	676	else: orientation = zeorientation ; frac = zefrac ; pad = 0.03 ; lu = 0.5
[910]	677	zelevpal = np.linspace(zevmin,zevmax,num=min([ticks/2+1,21]))
[917]	678	zecb = mpl.colorbar(fraction=frac,pad=pad,\
[910]	679	format=self.fmt,orientation=orientation,\
	680	ticks=zelevpal,\
[933]	681	extend='neither',spacing='proportional')
[910]	682	if zeorientation == "horizontal": zecb.ax.set_xlabel(self.title) ; self.title = ""
[933]	683	# colorbar title --> units
	684	if self.units not in ["dimless",""]:
	685	zecb.ax.set_title("["+self.units+"]",fontsize=3.*mpl.rcParams['font.size']/4.,x=lu,y=1.025)
	686
[910]	687	############################################################################################
	688	### VECTORS. must be after the colorbar. we could also leave possibility for streamlines.
	689	############################################################################################
	690	### not expecting NaN in self.addvecx and self.addvecy. masked arrays is just enough.
	691	stride = 3
	692	if self.addvecx is not None and self.addvecy is not None and self.mapmode:
[932]	693	### for metwinds only ???
[910]	694	[vecx,vecy] = m.rotate_vector(self.addvecx,self.addvecy,self.absc,self.ordi)
[932]	695	#vecx,vecy = self.addvecx,self.addvecy
[910]	696	# reference vector is scaled
	697	zescale = ppcompute.mean(np.sqrt(self.addvecxself.addvecx+self.addvecyself.addvecy))
	698	# make vector field
	699	q = m.quiver( x[::stride,::stride],y[::stride,::stride],\
	700	vecx[::stride,::stride],vecy[::stride,::stride],\
	701	angles='xy',color=self.colorvec,pivot='middle',\
	702	scale=zescale*reducevec,width=widthvec )
[930]	703	# make vector key.
[933]	704	#keyh = 1.025 ; keyv = 1.05 # upper right corner over colorbar
	705	keyh = 0.98 ; keyv = 1.06
[930]	706	#keyh = -0.03 ; keyv = 1.08 # upper left corner
[923]	707	p = mpl.quiverkey(q,keyh,keyv,\
[910]	708	zescale,str(int(zescale)),\
	709	color='black',labelpos='S',labelsep = 0.07)

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

Download in other formats: