Context Navigation

source: lmdz_wrf/trunk/tools/drawing_tools.py @ 255

Last change on this file since 255 was 254, checked in by lfita, 10 years ago
Adding 'K s-1' in 'units_lunits'
File size: 167.2 KB

Rev	Line
[192]	1	# -- coding: iso-8859-15 --
	2	#import pylab as plt
	3	# From http://stackoverflow.com/questions/13336823/matplotlib-python-error
	4	import numpy as np
	5	import matplotlib as mpl
	6	mpl.use('Agg')
	7	import matplotlib.pyplot as plt
	8	from mpl_toolkits.basemap import Basemap
	9	import os
	10	from netCDF4 import Dataset as NetCDFFile
	11
	12	errormsg = 'ERROR -- error -- ERROR -- error'
	13	warnmsg = 'WARNING -- waring -- WARNING -- warning'
	14
	15	fillValue = 1.e20
	16
	17	####### Funtions
	18	# searchInlist:
	19	# datetimeStr_datetime:
	20	# dateStr_date:
	21	# numVector_String:
	22	# timeref_datetime:
	23	# slice_variable:
	24	# interpolate_locs:
	25	# datetimeStr_conversion:
	26	# percendone:
	27	# netCDFdatetime_realdatetime:
	28	# file_nlines:
	29	# variables_values:
	30	# check_colorBar:
	31	# units_lunits:
	32	# ASCII_LaTeX:
	33	# pretty_int:
	34	# DegGradSec_deg:
	35	# intT2dt:
	36	# lonlat_values:
	37	# date_CFtime:
	38	# pot_values:
	39	# CFtimes_plot:
	40	# color_lines:
	41	# output_kind:
	42	# check_arguments:
	43	# Str_Bool:
	44	# plot_points:
	45	# plot_2Dfield:
	46	# plot_2Dfield_easy:
	47	# plot_topo_geogrid:
	48	# plot_topo_geogrid_boxes:
	49	# plot_2D_shadow:
	50	# plot_2D_shadow_time: Plotting a 2D field with one of the axes being time
	51	# plot_Neighbourghood_evol:Plotting neighbourghood evolution# plot_Trajectories
	52	# plot_2D_shadow_contour:
	53	# plot_2D_shadow_contour_time:
	54	# dxdy_lonlat: Function to provide lon/lat 2D lilke-matrices from any sort of dx,dy values
	55	# plot_2D_shadow_line:
	56	# plot_lines: Function to plot a collection of lines
	57
	58	# From nc_var_tools.py
	59	def searchInlist(listname, nameFind):
	60	""" Function to search a value within a list
	61	listname = list
	62	nameFind = value to find
	63	>>> searInlist(['1', '2', '3', '5'], '5')
	64	True
	65	"""
	66	for x in listname:
	67	if x == nameFind:
	68	return True
	69	return False
	70
	71	def datetimeStr_datetime(StringDT):
	72	""" Function to transform a string date ([YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format) to a date object
	73	>>> datetimeStr_datetime('1976-02-17_00:00:00')
	74	1976-02-17 00:00:00
	75	"""
	76	import datetime as dt
	77
	78	fname = 'datetimeStr_datetime'
	79
	80	dateD = np.zeros((3), dtype=int)
	81	timeT = np.zeros((3), dtype=int)
	82
	83	dateD[0] = int(StringDT[0:4])
	84	dateD[1] = int(StringDT[5:7])
	85	dateD[2] = int(StringDT[8:10])
	86
	87	trefT = StringDT.find(':')
	88	if not trefT == -1:
	89	# print ' ' + fname + ': refdate with time!'
	90	timeT[0] = int(StringDT[11:13])
	91	timeT[1] = int(StringDT[14:16])
	92	timeT[2] = int(StringDT[17:19])
	93
	94	if int(dateD[0]) == 0:
	95	print warnmsg
	96	print ' ' + fname + ': 0 reference year!! changing to 1'
	97	dateD[0] = 1
	98
	99	newdatetime = dt.datetime(dateD[0], dateD[1], dateD[2], timeT[0], timeT[1], timeT[2])
	100
	101	return newdatetime
	102
	103	def dateStr_date(StringDate):
	104	""" Function to transform a string date ([YYYY]-[MM]-[DD] format) to a date object
	105	>>> dateStr_date('1976-02-17')
	106	1976-02-17
	107	"""
	108	import datetime as dt
	109
	110	dateD = StringDate.split('-')
	111	if int(dateD[0]) == 0:
	112	print warnmsg
	113	print ' dateStr_date: 0 reference year!! changing to 1'
	114	dateD[0] = 1
	115	newdate = dt.date(int(dateD[0]), int(dateD[1]), int(dateD[2]))
	116	return newdate
	117
	118	def numVector_String(vec,char):
	119	""" Function to transform a vector of numbers to a single string [char] separated
	120	numVector_String(vec,char)
	121	vec= vector with the numerical values
	122	char= single character to split the values
	123	>>> print numVector_String(np.arange(10),' ')
	124	0 1 2 3 4 5 6 7 8 9
	125	"""
	126	fname = 'numVector_String'
	127
	128	if vec == 'h':
	129	print fname + '_____________________________________________________________'
	130	print numVector_String.__doc__
	131	quit()
	132
	133	Nvals = len(vec)
	134
	135	string=''
	136	for i in range(Nvals):
	137	if i == 0:
	138	string = str(vec[i])
	139	else:
	140	string = string + char + str(vec[i])
	141
	142	return string
	143
	144	def timeref_datetime(refd, timeval, tu):
	145	""" Function to transform from a [timeval] in [tu] units from the time referece [tref] to datetime object
	146	refd: time of reference (as datetime object)
	147	timeval: time value (as [tu] from [tref])
	148	tu: time units
	149	>>> timeref = date(1949,12,1,0,0,0)
	150	>>> timeref_datetime(timeref, 229784.36, hours)
	151	1976-02-17 08:21:36
	152	"""
	153	import datetime as dt
	154	import numpy as np
	155
	156	## Not in timedelta
	157	# if tu == 'years':
	158	# realdate = refdate + dt.timedelta(years=float(timeval))
	159	# elif tu == 'months':
	160	# realdate = refdate + dt.timedelta(months=float(timeval))
	161	if tu == 'weeks':
	162	realdate = refd + dt.timedelta(weeks=float(timeval))
	163	elif tu == 'days':
	164	realdate = refd + dt.timedelta(days=float(timeval))
	165	elif tu == 'hours':
	166	realdate = refd + dt.timedelta(hours=float(timeval))
	167	elif tu == 'minutes':
	168	realdate = refd + dt.timedelta(minutes=float(timeval))
	169	elif tu == 'seconds':
	170	realdate = refd + dt.timedelta(seconds=float(timeval))
	171	elif tu == 'milliseconds':
	172	realdate = refd + dt.timedelta(milliseconds=float(timeval))
	173	else:
	174	print errormsg
	175	print ' timeref_datetime: time units "' + tu + '" not ready!!!!'
	176	quit(-1)
	177
	178	return realdate
	179
	180	def slice_variable(varobj, dimslice):
	181	""" Function to return a slice of a given variable according to values to its
	182	dimensions
	183	slice_variable(varobj, dims)
	184	varobj= object wit the variable
	185	dimslice= [[dimname1]:[value1]\|[[dimname2]:[value2], ...] pairs of dimension
	186	[value]:
	187	* [integer]: which value of the dimension
	188	* -1: all along the dimension
	189	* [beg]:[end] slice from [beg] to [end]
	190	"""
	191	fname = 'slice_variable'
	192
	193	if varobj == 'h':
	194	print fname + '_____________________________________________________________'
	195	print slice_variable.__doc__
	196	quit()
	197
	198	vardims = varobj.dimensions
	199	Ndimvar = len(vardims)
	200
	201	Ndimcut = len(dimslice.split('\|'))
	202	dimsl = dimslice.split('\|')
	203
	204	varvalsdim = []
	205	dimnslice = []
	206
	207	for idd in range(Ndimvar):
	208	for idc in range(Ndimcut):
	209	dimcutn = dimsl[idc].split(':')[0]
	210	dimcutv = dimsl[idc].split(':')[1]
	211	if vardims[idd] == dimcutn:
	212	posfrac = dimcutv.find(':')
	213	if posfrac != -1:
	214	inifrac = int(dimcutv.split(':')[0])
	215	endfrac = int(dimcutv.split(':')[1])
	216	varvalsdim.append(slice(inifrac,endfrac))
	217	dimnslice.append(vardims[idd])
	218	else:
	219	if int(dimcutv) == -1:
	220	varvalsdim.append(slice(0,varobj.shape[idd]))
	221	dimnslice.append(vardims[idd])
	222	else:
	223	varvalsdim.append(int(dimcutv))
	224	break
	225
	226	varvalues = varobj[tuple(varvalsdim)]
	227
	228	return varvalues, dimnslice
	229
	230	def interpolate_locs(locs,coords,kinterp):
	231	""" Function to provide interpolate locations on a given axis
	232	interpolate_locs(locs,axis,kinterp)
	233	locs= locations to interpolate
	234	coords= axis values with the reference of coordinates
	235	kinterp: kind of interpolation
	236	'lin': linear
	237	>>> coordinates = np.arange((10), dtype=np.float)
	238	>>> values = np.array([-1.2, 2.4, 5.6, 7.8, 12.0])
	239	>>> interpolate_locs(values,coordinates,'lin')
	240	[ -1.2 2.4 5.6 7.8 13. ]
	241	>>> coordinates[0] = 0.5
	242	>>> coordinates[2] = 2.5
	243	>>> interpolate_locs(values,coordinates,'lin')
	244	[ -3.4 1.93333333 5.6 7.8 13. ]
	245	"""
	246
	247	fname = 'interpolate_locs'
	248
	249	if locs == 'h':
	250	print fname + '_____________________________________________________________'
	251	print interpolate_locs.__doc__
	252	quit()
	253
	254	Nlocs = locs.shape[0]
	255	Ncoords = coords.shape[0]
	256
	257	dcoords = coords[Ncoords-1] - coords[0]
	258
	259	intlocs = np.zeros((Nlocs), dtype=np.float)
	260	minc = np.min(coords)
	261	maxc = np.max(coords)
	262
	263	for iloc in range(Nlocs):
	264	for icor in range(Ncoords-1):
	265	if locs[iloc] < minc and dcoords > 0.:
	266	a = 0.
	267	b = 1. / (coords[1] - coords[0])
	268	c = coords[0]
	269	elif locs[iloc] > maxc and dcoords > 0.:
	270	a = (Ncoords-1)*1.
	271	b = 1. / (coords[Ncoords-1] - coords[Ncoords-2])
	272	c = coords[Ncoords-2]
	273	elif locs[iloc] < minc and dcoords < 0.:
	274	a = (Ncoords-1)*1.
	275	b = 1. / (coords[Ncoords-1] - coords[Ncoords-2])
	276	c = coords[Ncoords-2]
	277	elif locs[iloc] > maxc and dcoords < 0.:
	278	a = 0.
	279	b = 1. / (coords[1] - coords[0])
	280	c = coords[0]
	281	elif locs[iloc] >= coords[icor] and locs[iloc] < coords[icor+1] and dcoords > 0.:
	282	a = icor*1.
	283	b = 1. / (coords[icor+1] - coords[icor])
	284	c = coords[icor]
	285	print coords[icor], locs[iloc], coords[icor+1], ':', icor, '->', a, b
	286	elif locs[iloc] <= coords[icor] and locs[iloc] > coords[icor+1] and dcoords < 0.:
	287	a = icor*1.
	288	b = 1. / (coords[icor+1] - coords[icor])
	289	c = coords[icor]
	290
	291	if kinterp == 'lin':
	292	intlocs[iloc] = a + (locs[iloc] - c)*b
	293	else:
	294	print errormsg
	295	print ' ' + fname + ": interpolation kind '" + kinterp + "' not ready !!!!!"
	296	quit(-1)
	297
	298	return intlocs
	299
	300	def datetimeStr_conversion(StringDT,typeSi,typeSo):
	301	""" Function to transform a string date to an another date object
	302	StringDT= string with the date and time
	303	typeSi= type of datetime string input
	304	typeSo= type of datetime string output
	305	[typeSi/o]
	306	'cfTime': [time],[units]; ]time in CF-convention format [units] = [tunits] since [refdate]
	307	'matYmdHMS': numerical vector with [[YYYY], [MM], [DD], [HH], [MI], [SS]]
	308	'YmdHMS': [YYYY][MM][DD][HH][MI][SS] format
	309	'Y-m-d_H:M:S': [YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format
	310	'Y-m-d H:M:S': [YYYY]-[MM]-[DD] [HH]:[MI]:[SS] format
	311	'Y/m/d H-M-S': [YYYY]/[MM]/[DD] [HH]-[MI]-[SS] format
	312	'WRFdatetime': [Y], [Y], [Y], [Y], '-', [M], [M], '-', [D], [D], '_', [H],
	313	[H], ':', [M], [M], ':', [S], [S]
	314	>>> datetimeStr_conversion('1976-02-17_08:32:05','Y-m-d_H:M:S','matYmdHMS')
	315	[1976 2 17 8 32 5]
	316	>>> datetimeStr_conversion(str(137880)+',minutes since 1979-12-01_00:00:00','cfTime','Y/m/d H-M-S')
	317	1980/03/05 18-00-00
	318	"""
	319	import datetime as dt
	320
	321	fname = 'datetimeStr_conversion'
	322
	323	if StringDT[0:1] == 'h':
	324	print fname + '_____________________________________________________________'
	325	print datetimeStr_conversion.__doc__
	326	quit()
	327
	328	if typeSi == 'cfTime':
	329	timeval = np.float(StringDT.split(',')[0])
	330	tunits = StringDT.split(',')[1].split(' ')[0]
	331	Srefdate = StringDT.split(',')[1].split(' ')[2]
	332
	333	# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
	334	##
	335	yrref=Srefdate[0:4]
	336	monref=Srefdate[5:7]
	337	dayref=Srefdate[8:10]
	338
	339	trefT = Srefdate.find(':')
	340	if not trefT == -1:
	341	# print ' ' + fname + ': refdate with time!'
	342	horref=Srefdate[11:13]
	343	minref=Srefdate[14:16]
	344	secref=Srefdate[17:19]
	345	refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref + \
	346	'_' + horref + ':' + minref + ':' + secref)
	347	else:
	348	refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref + \
	349	+ '_00:00:00')
	350
	351	if tunits == 'weeks':
	352	newdate = refdate + dt.timedelta(weeks=float(timeval))
	353	elif tunits == 'days':
	354	newdate = refdate + dt.timedelta(days=float(timeval))
	355	elif tunits == 'hours':
	356	newdate = refdate + dt.timedelta(hours=float(timeval))
	357	elif tunits == 'minutes':
	358	newdate = refdate + dt.timedelta(minutes=float(timeval))
	359	elif tunits == 'seconds':
	360	newdate = refdate + dt.timedelta(seconds=float(timeval))
	361	elif tunits == 'milliseconds':
	362	newdate = refdate + dt.timedelta(milliseconds=float(timeval))
	363	else:
	364	print errormsg
	365	print ' timeref_datetime: time units "' + tunits + '" not ready!!!!'
	366	quit(-1)
	367
	368	yr = newdate.year
	369	mo = newdate.month
	370	da = newdate.day
	371	ho = newdate.hour
	372	mi = newdate.minute
	373	se = newdate.second
	374	elif typeSi == 'matYmdHMS':
	375	yr = StringDT[0]
	376	mo = StringDT[1]
	377	da = StringDT[2]
	378	ho = StringDT[3]
	379	mi = StringDT[4]
	380	se = StringDT[5]
	381	elif typeSi == 'YmdHMS':
	382	yr = int(StringDT[0:4])
	383	mo = int(StringDT[4:6])
	384	da = int(StringDT[6:8])
	385	ho = int(StringDT[8:10])
	386	mi = int(StringDT[10:12])
	387	se = int(StringDT[12:14])
	388	elif typeSi == 'Y-m-d_H:M:S':
	389	dateDT = StringDT.split('_')
	390	dateD = dateDT[0].split('-')
	391	timeT = dateDT[1].split(':')
	392	yr = int(dateD[0])
	393	mo = int(dateD[1])
	394	da = int(dateD[2])
	395	ho = int(timeT[0])
	396	mi = int(timeT[1])
	397	se = int(timeT[2])
	398	elif typeSi == 'Y-m-d H:M:S':
	399	dateDT = StringDT.split(' ')
	400	dateD = dateDT[0].split('-')
	401	timeT = dateDT[1].split(':')
	402	yr = int(dateD[0])
	403	mo = int(dateD[1])
	404	da = int(dateD[2])
	405	ho = int(timeT[0])
	406	mi = int(timeT[1])
	407	se = int(timeT[2])
	408	elif typeSi == 'Y/m/d H-M-S':
	409	dateDT = StringDT.split(' ')
	410	dateD = dateDT[0].split('/')
	411	timeT = dateDT[1].split('-')
	412	yr = int(dateD[0])
	413	mo = int(dateD[1])
	414	da = int(dateD[2])
	415	ho = int(timeT[0])
	416	mi = int(timeT[1])
	417	se = int(timeT[2])
	418	elif typeSi == 'WRFdatetime':
	419	yr = int(StringDT[0])1000 + int(StringDT[1])100 + int(StringDT[2])*10 + \
	420	int(StringDT[3])
	421	mo = int(StringDT[5])*10 + int(StringDT[6])
	422	da = int(StringDT[8])*10 + int(StringDT[9])
	423	ho = int(StringDT[11])*10 + int(StringDT[12])
	424	mi = int(StringDT[14])*10 + int(StringDT[15])
	425	se = int(StringDT[17])*10 + int(StringDT[18])
	426	else:
	427	print errormsg
	428	print ' ' + fname + ': type of String input date "' + typeSi + \
	429	'" not ready !!!!'
	430	quit(-1)
	431
	432	if typeSo == 'matYmdHMS':
	433	dateYmdHMS = np.zeros((6), dtype=int)
	434	dateYmdHMS[0] = yr
	435	dateYmdHMS[1] = mo
	436	dateYmdHMS[2] = da
	437	dateYmdHMS[3] = ho
	438	dateYmdHMS[4] = mi
	439	dateYmdHMS[5] = se
	440	elif typeSo == 'YmdHMS':
	441	dateYmdHMS = str(yr).zfill(4) + str(mo).zfill(2) + str(da).zfill(2) + \
	442	str(ho).zfill(2) + str(mi).zfill(2) + str(se).zfill(2)
	443	elif typeSo == 'Y-m-d_H:M:S':
	444	dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' + \
	445	str(da).zfill(2) + '_' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \
	446	str(se).zfill(2)
	447	elif typeSo == 'Y-m-d H:M:S':
	448	dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' + \
	449	str(da).zfill(2) + ' ' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \
	450	str(se).zfill(2)
	451	elif typeSo == 'Y/m/d H-M-S':
	452	dateYmdHMS = str(yr).zfill(4) + '/' + str(mo).zfill(2) + '/' + \
	453	str(da).zfill(2) + ' ' + str(ho).zfill(2) + '-' + str(mi).zfill(2) + '-' + \
	454	str(se).zfill(2)
	455	elif typeSo == 'WRFdatetime':
	456	dateYmdHMS = []
	457	yM = yr/1000
	458	yC = (yr-yM*1000)/100
	459	yD = (yr-yM1000-yC100)/10
	460	yU = yr-yM1000-yC100-yD*10
	461
	462	mD = mo/10
	463	mU = mo-mD*10
	464
	465	dD = da/10
	466	dU = da-dD*10
	467
	468	hD = ho/10
	469	hU = ho-hD*10
	470
	471	miD = mi/10
	472	miU = mi-miD*10
	473
	474	sD = se/10
	475	sU = se-sD*10
	476
	477	dateYmdHMS.append(str(yM))
	478	dateYmdHMS.append(str(yC))
	479	dateYmdHMS.append(str(yD))
	480	dateYmdHMS.append(str(yU))
	481	dateYmdHMS.append('-')
	482	dateYmdHMS.append(str(mD))
	483	dateYmdHMS.append(str(mU))
	484	dateYmdHMS.append('-')
	485	dateYmdHMS.append(str(dD))
	486	dateYmdHMS.append(str(dU))
	487	dateYmdHMS.append('_')
	488	dateYmdHMS.append(str(hD))
	489	dateYmdHMS.append(str(hU))
	490	dateYmdHMS.append(':')
	491	dateYmdHMS.append(str(miD))
	492	dateYmdHMS.append(str(miU))
	493	dateYmdHMS.append(':')
	494	dateYmdHMS.append(str(sD))
	495	dateYmdHMS.append(str(sU))
	496	else:
	497	print errormsg
	498	print ' ' + fname + ': type of output date "' + typeSo + '" not ready !!!!'
	499	quit(-1)
	500
	501	return dateYmdHMS
	502
	503	def percendone(nvals,tot,percen,msg):
	504	""" Function to provide the percentage of an action across the matrix
	505	nvals=number of values
	506	tot=total number of values
	507	percen=percentage frequency for which the message is wanted
	508	msg= message
	509	"""
	510	from sys import stdout
	511
	512	num = int(tot * percen/100)
	513	if (nvals%num == 0):
	514	print '\r ' + msg + '{0:8.3g}'.format(nvals*100./tot) + ' %',
	515	stdout.flush()
	516
	517	return ''
	518
	519	def netCDFdatetime_realdatetime(units, tcalendar, times):
	520	""" Function to transfrom from netCDF CF-compilant times to real time
	521	"""
	522	import datetime as dt
	523
	524	txtunits = units.split(' ')
	525	tunits = txtunits[0]
	526	Srefdate = txtunits[len(txtunits) - 1]
	527
	528	# Calendar type
	529	##
	530	is360 = False
	531	if tcalendar is not None:
	532	print ' netCDFdatetime_realdatetime: There is a calendar attribute'
	533	if tcalendar == '365_day' or tcalendar == 'noleap':
	534	print ' netCDFdatetime_realdatetime: No leap years!'
	535	isleapcal = False
	536	elif tcalendar == 'proleptic_gregorian' or tcalendar == 'standard' or tcalendar == 'gregorian':
	537	isleapcal = True
	538	elif tcalendar == '360_day':
	539	is360 = True
	540	isleapcal = False
	541	else:
	542	print errormsg
	543	print ' netCDFdatetime_realdatetime: Calendar "' + tcalendar + '" not prepared!'
	544	quit(-1)
	545
	546	# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
	547	##
	548	timeval = Srefdate.find(':')
	549
	550	if not timeval == -1:
	551	print ' netCDFdatetime_realdatetime: refdate with time!'
	552	refdate = datetimeStr_datetime(Srefdate)
	553	else:
	554	refdate = dateStr_date(Srefdate + '_00:00:00')
	555
	556	dimt = len(times)
	557	# datetype = type(dt.datetime(1972,02,01))
	558	# realdates = np.array(dimt, datetype)
	559	# print realdates
	560
	561	## Not in timedelta
	562	# if tunits == 'years':
	563	# for it in range(dimt):
	564	# realdate = refdate + dt.timedelta(years=float(times[it]))
	565	# realdates[it] = int(realdate.year)
	566	# elif tunits == 'months':
	567	# for it in range(dimt):
	568	# realdate = refdate + dt.timedelta(months=float(times[it]))
	569	# realdates[it] = int(realdate.year)
	570	# realdates = []
	571	realdates = np.zeros((dimt, 6), dtype=int)
	572	if tunits == 'weeks':
	573	for it in range(dimt):
	574	realdate = refdate + dt.timedelta(weeks=float(times[it]))
	575	realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
	576	elif tunits == 'days':
	577	for it in range(dimt):
	578	realdate = refdate + dt.timedelta(days=float(times[it]))
	579	realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
	580	elif tunits == 'hours':
	581	for it in range(dimt):
	582	realdate = refdate + dt.timedelta(hours=float(times[it]))
	583	# if not isleapcal:
	584	# Nleapdays = cal.leapdays(int(refdate.year), int(realdate.year))
	585	# realdate = realdate - dt.timedelta(days=Nleapdays)
	586	# if is360:
	587	# Nyears360 = int(realdate.year) - int(refdate.year) + 1
	588	# realdate = realdate -dt.timedelta(days=Nyears360*5)
	589	# realdates[it] = realdate
	590	# realdates = refdate + dt.timedelta(hours=float(times))
	591	realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
	592	elif tunits == 'minutes':
	593	for it in range(dimt):
	594	realdate = refdate + dt.timedelta(minutes=float(times[it]))
	595	realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
	596	elif tunits == 'seconds':
	597	for it in range(dimt):
	598	realdate = refdate + dt.timedelta(seconds=float(times[it]))
	599	realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
	600	elif tunits == 'milliseconds':
	601	for it in range(dimt):
	602	realdate = refdate + dt.timedelta(milliseconds=float(times[it]))
	603	realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
	604	elif tunits == 'microseconds':
	605	for it in range(dimt):
	606	realdate = refdate + dt.timedelta(microseconds=float(times[it]))
	607	realdates[it,:]=[realdate.year, realdate.month, realdate.day, realdate.hour, realdate.minute, realdate.second]
	608	else:
	609	print errormsg
	610	print ' netCDFdatetime_realdatetime: time units "' + tunits + '" is not ready!!!'
	611	quit(-1)
	612
	613	return realdates
	614
	615	def file_nlines(filen):
	616	""" Function to provide the number of lines of a file
	617	filen= name of the file
	618	>>> file_nlines('trajectory.dat')
	619	49
	620	"""
	621	fname = 'file_nlines'
	622
	623	if not os.path.isfile(filen):
	624	print errormsg
	625	print ' ' + fname + ' file: "' + filen + '" does not exist !!'
	626	quit(-1)
	627
	628	fo = open(filen,'r')
	629
	630	nlines=0
	631	for line in fo: nlines = nlines + 1
	632
	633	fo.close()
	634
	635	return nlines
	636
	637	def realdatetime1_CFcompilant(time, Srefdate, tunits):
	638	""" Function to transform a matrix with a real time value ([year, month, day,
	639	hour, minute, second]) to a netCDF one
	640	time= matrix with time
	641	Srefdate= reference date ([YYYY][MM][DD][HH][MI][SS] format)
	642	tunits= units of time respect to Srefdate
	643	>>> realdatetime1_CFcompilant([1976, 2, 17, 8, 20, 0], '19491201000000', 'hours')
	644	229784.33333333
	645	"""
	646
	647	import datetime as dt
	648	yrref=int(Srefdate[0:4])
	649	monref=int(Srefdate[4:6])
	650	dayref=int(Srefdate[6:8])
	651	horref=int(Srefdate[8:10])
	652	minref=int(Srefdate[10:12])
	653	secref=int(Srefdate[12:14])
	654
	655	refdate=dt.datetime(yrref, monref, dayref, horref, minref, secref)
	656
	657	if tunits == 'weeks':
	658	cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5])-refdate
	659	cfdates = (cfdate.days + cfdate.seconds/(3600.*24.))/7.
	660	elif tunits == 'days':
	661	cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
	662	cfdates = cfdate.days + cfdate.seconds/(3600.*24.)
	663	elif tunits == 'hours':
	664	cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
	665	cfdates = cfdate.days*24. + cfdate.seconds/3600.
	666	elif tunits == 'minutes':
	667	cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
	668	cfdates = cfdate.days24.60. + cfdate.seconds/60.
	669	elif tunits == 'seconds':
	670	cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
	671	cfdates = cfdate.days24.3600. + cfdate.seconds
	672	elif tunits == 'milliseconds':
	673	cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],time[5]) - refdate
	674	cfdates = cfdate.days1000.24.3600. + cfdate.seconds1000.
	675	elif tunits == 'microseconds':
	676	cfdate = dt.datetime(time[0],time[1],time[2],time[3],time[4],times[5]) - refdate
	677	cfdates = cfdate.days1000000.24.3600. + cfdate.seconds1000000.
	678	else:
	679	print errormsg
	680	print ' ' + fname + ': time units "' + tunits + '" is not ready!!!'
	681	quit(-1)
	682
	683	return cfdates
	684
	685	def basicvardef(varobj, vstname, vlname, vunits):
	686	""" Function to give the basic attributes to a variable
	687	varobj= netCDF variable object
	688	vstname= standard name of the variable
	689	vlname= long name of the variable
	690	vunits= units of the variable
	691	"""
	692	attr = varobj.setncattr('standard_name', vstname)
	693	attr = varobj.setncattr('long_name', vlname)
	694	attr = varobj.setncattr('units', vunits)
	695
	696	return
	697
	698	def variables_values(varName):
	699	""" Function to provide values to plot the different variables
	700	variables_values(varName)
	701	[varName]= name of the variable
	702	return: [var name], [std name], [minimum], [maximum],
	703	[long name]('\|' for spaces), [units], [color palette] (following:
	704	http://matplotlib.org/1.3.1/examples/color/colormaps_reference.html)
	705	[varn]: original name of the variable
	706	NOTE: It might be better doing it with an external ASII file. But then we
	707	got an extra dependency...
	708	>>> variables_values('WRFght')
	709	['z', 'geopotential_height', 0.0, 80000.0, 'geopotential\|height', 'm2s-2', 'rainbow']
	710	"""
	711	fname='variables_values'
	712
	713	if varName == 'h':
	714	print fname + '_____________________________________________________________'
	715	print variables_values.__doc__
	716	quit()
	717
	718	# Variable name might come with a statistical surname...
	719	stats=['min','max','mean','stdv']
	720
	721	ifst = False
	722	for st in stats:
	723	if varName.find(st) > -1:
	724	print ' '+ fname + ": varibale '" + varName + "' with a statistical "+\
	725	" surname: '",st,"' !!"
	726	Lst = len(st)
	727	LvarName = len(varName)
	728	varn = varName[0:LvarName - Lst]
	729	ifst = True
	730	break
	731	if not ifst:
	732	varn = varName
	733
	734	if varn[0:6] == 'varDIM':
	735	# Variable from a dimension (all with 'varDIM' prefix)
	736	Lvarn = len(varn)
	737	varvals = [varn[6:Lvarn+1], varn[6:Lvarn+1], 0., 1., \
	738	"variable\|from\|size\|of\|dimension\|'" + varn[6:Lvarn+1] + "'", '1', 'rainbox']
	739	elif varn == 'a_tht' or varn == 'LA_THT':
[238]	740	varvals = ['ath', 'total_thermal_plume_cover', 0., 1., \
[192]	741	'total\|column\|thermal\|plume\|cover', '1', 'YlGnBu']
	742	elif varn == 'bils' or varn == 'LBILS':
	743	varvals = ['bils', 'surface_total_heat_flux', -100., 100., \
	744	'surface\|total\|heat\|flux', 'Wm-2', 'seismic']
	745	elif varn == 'landcat' or varn == 'category':
	746	varvals = ['landcat', 'land_categories', 0., 22., 'land\|categories', '1', \
	747	'rainbow']
	748	elif varn == 'c' or varn == 'QCLOUD' or varn == 'oliq' or varn == 'OLIQ':
	749	varvals = ['c', 'condensed_water_mixing_ratio', 0., 3.e-4, \
	750	'condensed\|water\|mixing\|ratio', 'kgkg-1', 'BuPu']
[238]	751	elif varn == 'ci' or varn == 'iwcon' or varn == 'LIWCON':
	752	varvals = ['ci', 'cloud_iced_water_mixing_ratio', 0., 0.0003, \
	753	'cloud\|iced\|water\|mixing\|ratio', 'kgkg-1', 'Purples']
	754	elif varn == 'cl' or varn == 'lwcon' or varn == 'LLWCON':
	755	varvals = ['cl', 'cloud_liquidwater_mixing_ratio', 0., 0.0003, \
	756	'cloud\|liquid\|water\|mixing\|ratio', 'kgkg-1', 'Blues']
[233]	757	elif varn == 'cld' or varn == 'CLDFRA' or varn == 'rneb' or varn == 'lrneb' or \
	758	varn == 'LRNEB':
	759	varvals = ['cld', 'cloud_area_fraction', 0., 1., 'cloud\|fraction', '1', \
	760	'gist_gray']
[234]	761	elif varn == 'cldc' or varn == 'rnebcon' or varn == 'lrnebcon' or \
	762	varn == 'LRNEBCON':
[233]	763	varvals = ['cldc', 'convective_cloud_area_fraction', 0., 1., \
	764	'convective\|cloud\|fraction', '1', 'gist_gray']
	765	elif varn == 'cldl' or varn == 'rnebls' or varn == 'lrnebls' or varn == 'LRNEBLS':
	766	varvals = ['cldl', 'large_scale_cloud_area_fraction', 0., 1., \
	767	'large\|scale\|cloud\|fraction', '1', 'gist_gray']
[192]	768	elif varn == 'clt' or varn == 'CLT' or varn == 'cldt' or \
	769	varn == 'Total cloudiness':
	770	varvals = ['clt', 'cloud_area_fraction', 0., 1., 'total\|cloud\|cover', '1', \
	771	'gist_gray']
	772	elif varn == 'cll' or varn == 'cldl' or varn == 'LCLDL' or \
	773	varn == 'Low-level cloudiness':
	774	varvals = ['cll', 'low_level_cloud_area_fraction', 0., 1., \
	775	'low\|level\|(p\|>\|680\|hPa)\|cloud\|fraction', '1', 'gist_gray']
	776	elif varn == 'clm' or varn == 'cldm' or varn == 'LCLDM' or \
	777	varn == 'Mid-level cloudiness':
	778	varvals = ['clm', 'mid_level_cloud_area_fraction', 0., 1., \
	779	'medium\|level\|(440\|<\|p\|<\|680\|hPa)\|cloud\|fraction', '1', 'gist_gray']
	780	elif varn == 'clh' or varn == 'cldh' or varn == 'LCLDH' or \
	781	varn == 'High-level cloudiness':
	782	varvals = ['clh', 'high_level_cloud_area_fraction', 0., 1., \
	783	'high\|level\|(p\|<\|440\|hPa)\|cloud\|fraction', '1', 'gist_gray']
[249]	784	elif varn == 'clmf' or varn == 'fbase' or varn == 'LFBASE':
	785	varvals = ['clmf', 'cloud_base_max_flux', -0.3, 0.3, 'cloud\|base\|max\|flux', \
	786	'kgm-2s-1', 'seismic']
	787	elif varn == 'clp' or varn == 'pbase' or varn == 'LPBASE':
	788	varvals = ['clp', 'cloud_base_pressure', -0.3, 0.3, 'cloud\|base\|pressure', \
	789	'Pa', 'Reds']
	790	elif varn == 'cpt' or varn == 'ptconv' or varn == 'LPTCONV':
	791	varvals = ['cpt', 'convective_point', 0., 1., 'convective\|point', '1', \
	792	'seismic']
[192]	793	elif varn == 'dqajs' or varn == 'LDQAJS':
	794	varvals = ['dqajs', 'dry_adjustment_water_vapor_tendency', -0.0003, 0.0003, \
	795	'dry\|adjustment\|water\|vapor\|tendency', 'kg/kg/s', 'seismic']
	796	elif varn == 'dqcon' or varn == 'LDQCON':
	797	varvals = ['dqcon', 'convective_water_vapor_tendency', -3e-8, 3.e-8, \
	798	'convective\|water\|vapor\|tendency', 'kg/kg/s', 'seismic']
	799	elif varn == 'dqdyn' or varn == 'LDQDYN':
	800	varvals = ['dqdyn', 'dynamics_water_vapor_tendency', -3.e-7, 3.e-7, \
	801	'dynamics\|water\|vapor\|tendency', 'kg/kg/s', 'seismic']
	802	elif varn == 'dqeva' or varn == 'LDQEVA':
	803	varvals = ['dqeva', 'evaporation_water_vapor_tendency', -3.e-6, 3.e-6, \
	804	'evaporation\|water\|vapor\|tendency', 'kg/kg/s', 'seismic']
	805	elif varn == 'dqlscst' or varn == 'LDQLSCST':
	806	varvals = ['dqlscst', 'stratocumulus_water_vapor_tendency', -3.e-7, 3.e-7, \
	807	'stratocumulus\|water\|vapor\|tendency', 'kg/kg/s', 'seismic']
	808	elif varn == 'dqlscth' or varn == 'LDQLSCTH':
	809	varvals = ['dqlscth', 'thermals_water_vapor_tendency', -3.e-7, 3.e-7, \
	810	'thermal\|plumes\|water\|vapor\|tendency', 'kg/kg/s', 'seismic']
	811	elif varn == 'dqlsc' or varn == 'LDQLSC':
	812	varvals = ['dqlsc', 'condensation_water_vapor_tendency', -3.e-6, 3.e-6, \
	813	'condensation\|water\|vapor\|tendency', 'kg/kg/s', 'seismic']
	814	elif varn == 'dqphy' or varn == 'LDQPHY':
	815	varvals = ['dqphy', 'physics_water_vapor_tendency', -3.e-7, 3.e-7, \
	816	'physics\|water\|vapor\|tendency', 'kg/kg/s', 'seismic']
	817	elif varn == 'dqthe' or varn == 'LDQTHE':
	818	varvals = ['dqthe', 'thermals_water_vapor_tendency', -3.e-7, 3.e-7, \
	819	'thermal\|plumes\|water\|vapor\|tendency', 'kg/kg/s', 'seismic']
	820	elif varn == 'dqvdf' or varn == 'LDQVDF':
	821	varvals = ['dqvdf', 'vertical_difussion_water_vapor_tendency', -3.e-8, 3.e-8,\
	822	'vertical\|difussion\|water\|vapor\|tendency', 'kg/kg/s', 'seismic']
	823	elif varn == 'dqwak' or varn == 'LDQWAK':
	824	varvals = ['dqwak', 'wake_water_vapor_tendency', -3.e-7, 3.e-7, \
	825	'wake\|water\|vapor\|tendency', 'kg/kg/s', 'seismic']
[249]	826	elif varn == 'dta' or varn == 'tnt' or varn == 'LTNT':
	827	varvals = ['dta', 'tendency_air_temperature', -3.e-3, 3.e-3, \
	828	'tendency\|of\|air\|temperature', 'K/s', 'seismic']
	829	elif varn == 'dtac' or varn == 'tntc' or varn == 'LTNTC':
[250]	830	varvals = ['dtac', 'moist_convection_tendency_air_temperature', -3.e-3, \
	831	3.e-3, 'moist\|convection\|tendency\|of\|air\|temperature', 'K/s', 'seismic']
[249]	832	elif varn == 'dtar' or varn == 'tntr' or varn == 'LTNTR':
	833	varvals = ['dtar', 'radiative_heating_tendency_air_temperature', -3.e-3, \
	834	3.e-3, 'radiative\|heating\|tendency\|of\|air\|temperature', 'K/s', 'seismic']
	835	elif varn == 'dtascpbl' or varn == 'tntscpbl' or varn == 'LTNTSCPBL':
	836	varvals = ['dtascpbl', \
	837	'stratiform_cloud_precipitation_BL_mixing_tendency_air_temperature', \
	838	-3.e-6, 3.e-6, \
	839	'stratiform\|cloud\|precipitation\|Boundary\|Layer\|mixing\|tendency\|air\|' +
	840	'temperature', 'K/s', 'seismic']
[192]	841	elif varn == 'dtajs' or varn == 'LDTAJS':
	842	varvals = ['dtajs', 'dry_adjustment_thermal_tendency', -3.e-5, 3.e-5, \
	843	'dry\|adjustment\|thermal\|tendency', 'K/s', 'seismic']
	844	elif varn == 'dtcon' or varn == 'LDTCON':
	845	varvals = ['dtcon', 'convective_thermal_tendency', -3.e-5, 3.e-5, \
	846	'convective\|thermal\|tendency', 'K/s', 'seismic']
	847	elif varn == 'dtdyn' or varn == 'LDTDYN':
	848	varvals = ['dtdyn', 'dynamics_thermal_tendency', -3.e-4, 3.e-4, \
	849	'dynamics\|thermal\|tendency', 'K/s', 'seismic']
	850	elif varn == 'dteva' or varn == 'LDTEVA':
	851	varvals = ['dteva', 'evaporation_thermal_tendency', -3.e-3, 3.e-3, \
	852	'evaporation\|thermal\|tendency', 'K/s', 'seismic']
	853	elif varn == 'dtlscst' or varn == 'LDTLSCST':
	854	varvals = ['dtlscst', 'stratocumulus_thermal_tendency', -3.e-4, 3.e-4, \
	855	'stratocumulus\|thermal\|tendency', 'K/s', 'seismic']
	856	elif varn == 'dtlscth' or varn == 'LDTLSCTH':
	857	varvals = ['dtlscth', 'thermals_thermal_tendency', -3.e-4, 3.e-4, \
	858	'thermal\|plumes\|thermal\|tendency', 'K/s', 'seismic']
	859	elif varn == 'dtlsc' or varn == 'LDTLSC':
	860	varvals = ['dtlsc', 'condensation_thermal_tendency', -3.e-3, 3.e-3, \
	861	'condensation\|thermal\|tendency', 'K/s', 'seismic']
[249]	862	elif varn == 'dtlwr' or varn == 'LDTLWR':
	863	varvals = ['dtlwr', 'long_wave_thermal_tendency', -3.e-3, 3.e-3, \
	864	'long\|wave\|radiation\|thermal\|tendency', 'K/s', 'seismic']
[192]	865	elif varn == 'dtphy' or varn == 'LDTPHY':
	866	varvals = ['dtphy', 'physics_thermal_tendency', -3.e-4, 3.e-4, \
	867	'physics\|thermal\|tendency', 'K/s', 'seismic']
[249]	868	elif varn == 'dtsw0' or varn == 'LDTSW0':
	869	varvals = ['dtsw0', 'cloudy_sky_short_wave_thermal_tendency', -3.e-4, 3.e-4, \
	870	'cloudy\|sky\|short\|wave\|radiation\|thermal\|tendency', 'K/s', 'seismic']
[192]	871	elif varn == 'dtthe' or varn == 'LDTTHE':
	872	varvals = ['dtthe', 'thermals_thermal_tendency', -3.e-4, 3.e-4, \
	873	'thermal\|plumes\|thermal\|tendency', 'K/s', 'seismic']
	874	elif varn == 'dtvdf' or varn == 'LDTVDF':
	875	varvals = ['dtvdf', 'vertical_difussion_thermal_tendency', -3.e-5, 3.e-5, \
	876	'vertical\|difussion\|thermal\|tendency', 'K/s', 'seismic']
	877	elif varn == 'dtwak' or varn == 'LDTWAK':
	878	varvals = ['dtwak', 'wake_thermal_tendency', -3.e-4, 3.e-4, \
	879	'wake\|thermal\|tendency', 'K/s', 'seismic']
[206]	880	elif varn == 'ducon' or varn == 'LDUCON':
	881	varvals = ['ducon', 'convective_eastward_wind_tendency', -3.e-3, 3.e-3, \
	882	'convective\|eastward\|wind\|tendency', 'ms-2', 'seismic']
	883	elif varn == 'dudyn' or varn == 'LDUDYN':
	884	varvals = ['dudyn', 'dynamics_eastward_wind_tendency', -3.e-3, 3.e-3, \
	885	'dynamics\|eastward\|wind\|tendency', 'ms-2', 'seismic']
	886	elif varn == 'duvdf' or varn == 'LDUVDF':
	887	varvals = ['duvdf', 'vertical_difussion_eastward_wind_tendency', -3.e-3, \
	888	3.e-3, 'vertical\|difussion\|eastward\|wind\|tendency', 'ms-2', 'seismic']
	889	elif varn == 'dvcon' or varn == 'LDVCON':
	890	varvals = ['dvcon', 'convective_difussion_northward_wind_tendency', -3.e-3, \
	891	3.e-3, 'convective\|northward\|wind\|tendency', 'ms-2', 'seismic']
	892	elif varn == 'dvdyn' or varn == 'LDVDYN':
	893	varvals = ['dvdyn', 'dynamics_northward_wind_tendency', -3.e-3, \
	894	3.e-3, 'dynamics\|difussion\|northward\|wind\|tendency', 'ms-2', 'seismic']
	895	elif varn == 'dvvdf' or varn == 'LDVVDF':
	896	varvals = ['dvvdf', 'vertical_difussion_northward_wind_tendency', -3.e-3, \
	897	3.e-3, 'vertical\|difussion\|northward\|wind\|tendency', 'ms-2', 'seismic']
[192]	898	elif varn == 'evspsbl' or varn == 'LEVAP' or varn == 'evap':
	899	varvals = ['evspsbl', 'water_evaporation_flux', 0., 1.5e-4, \
	900	'water\|evaporation\|flux', 'kgm-2s-1', 'Blues']
[228]	901	elif varn == 'g' or varn == 'QGRAUPEL':
	902	varvals = ['g', 'grauepl_mixing_ratio', 0., 0.0003, 'graupel\|mixing\|ratio', \
	903	'kgkg-1', 'Purples']
[192]	904	elif varn == 'h2o' or varn == 'LH2O':
	905	varvals = ['h2o', 'water_mass_fraction', 0., 3.e-2, \
	906	'mass\|fraction\|of\|water', '1', 'Blues']
[228]	907	elif varn == 'h' or varn == 'QHAIL':
	908	varvals = ['h', 'hail_mixing_ratio', 0., 0.0003, 'hail\|mixing\|ratio', \
	909	'kgkg-1', 'Purples']
[192]	910	elif varn == 'hfls' or varn == 'LH' or varn == 'LFLAT' or varn == 'flat':
	911	varvals = ['hfls', 'surface_upward_latent_heat_flux', -400., 400., \
	912	'upward\|latnt\|heat\|flux\|at\|the\|surface', 'Wm-2', 'seismic']
	913	elif varn == 'hfss' or varn == 'LSENS' or varn == 'sens':
	914	varvals = ['hfss', 'surface_upward_sensible_heat_flux', -150., 150., \
	915	'upward\|sensible\|heat\|flux\|at\|the\|surface', 'Wm-2', 'seismic']
	916	elif varn == 'hus' or varn == 'WRFrh' or varn == 'LMDZrh' or varn == 'rhum' or \
	917	varn == 'LRHUM':
	918	varvals = ['hus', 'specific_humidity', 0., 1., 'specific\|humidty', '1', \
	919	'BuPu']
	920	elif varn == 'huss' or varn == 'WRFrhs' or varn == 'LMDZrhs' or varn == 'rh2m' or\
	921	varn == 'LRH2M':
	922	varvals = ['huss', 'specific_humidity', 0., 1., 'specific\|humidty\|at\|2m', \
	923	'1', 'BuPu']
[238]	924	elif varn == 'i' or varn == 'QICE':
[228]	925	varvals = ['i', 'iced_water_mixing_ratio', 0., 0.0003, \
[206]	926	'iced\|water\|mixing\|ratio', 'kgkg-1', 'Purples']
[192]	927	elif varn == 'lat' or varn == 'XLAT' or varn == 'XLAT_M' or varn == 'latitude':
	928	varvals = ['lat', 'latitude', -90., 90., 'latitude', 'degrees North', \
	929	'seismic']
[233]	930	elif varn == 'lcl' or varn == 's_lcl' or varn == 'ls_lcl' or varn == 'LS_LCL':
	931	varvals = ['lcl', 'condensation_level', 0., 2500., 'level\|of\|condensation', \
	932	'm', 'Greens']
[249]	933	elif varn == 'lambdath' or varn == 'lambda_th' or varn == 'LLAMBDA_TH':
	934	varvals = ['lambdath', 'thermal_plume_vertical_velocity', -30., 30., \
[252]	935	'thermal\|plume\|vertical\|velocity', 'm/s', 'seismic']
[249]	936	elif varn == 'lmaxth' or varn == 'LLMAXTH':
	937	varvals = ['lmaxth', 'upper_level_thermals', 0., 100., 'upper\|level\|thermals'\
	938	, '1', 'Greens']
[192]	939	elif varn == 'lon' or varn == 'XLONG' or varn == 'XLONG_M':
	940	varvals = ['lon', 'longitude', -180., 180., 'longitude', 'degrees East', \
	941	'seismic']
	942	elif varn == 'longitude':
	943	varvals = ['lon', 'longitude', 0., 360., 'longitude', 'degrees East', \
	944	'seismic']
	945	elif varn == 'orog' or varn == 'HGT' or varn == 'HGT_M':
	946	varvals = ['orog', 'orography', 0., 3000., 'surface\|altitude', 'm','terrain']
[249]	947	elif varn == 'pfc' or varn == 'plfc' or varn == 'LPLFC':
	948	varvals = ['pfc', 'pressure_free_convection', 100., 1100., \
	949	'pressure\|free\|convection', 'hPa', 'BuPu']
	950	elif varn == 'plcl' or varn == 'LPLCL':
	951	varvals = ['plcl', 'pressure_lifting_condensation_level', 700., 1100., \
	952	'pressure\|lifting\|condensation\|level', 'hPa', 'BuPu']
[192]	953	elif varn == 'pr' or varn == 'RAINTOT' or varn == 'precip' or \
	954	varn == 'LPRECIP' or varn == 'Precip Totale liq+sol':
	955	varvals = ['pr', 'precipitation_flux', 0., 1.e-4, 'precipitation\|flux', \
	956	'kgm-2s-1', 'BuPu']
[220]	957	elif varn == 'prprof' or varn == 'vprecip' or varn == 'LVPRECIP':
[227]	958	varvals = ['prprof', 'precipitation_profile', 0., 1.e-3, \
	959	'precipitation\|profile', 'kg/m2/s', 'BuPu']
	960	elif varn == 'prprofci' or varn == 'pr_con_i' or varn == 'LPR_CON_I':
	961	varvals = ['prprofci', 'precipitation_profile_convective_i', 0., 1.e-3, \
	962	'precipitation\|profile\|convective\|i', 'kg/m2/s', 'BuPu']
	963	elif varn == 'prprofcl' or varn == 'pr_con_l' or varn == 'LPR_CON_L':
	964	varvals = ['prprofcl', 'precipitation_profile_convective_l', 0., 1.e-3, \
	965	'precipitation\|profile\|convective\|l', 'kg/m2/s', 'BuPu']
	966	elif varn == 'prprofli' or varn == 'pr_lsc_i' or varn == 'LPR_LSC_I':
	967	varvals = ['prprofli', 'precipitation_profile_large_scale_i', 0., 1.e-3, \
	968	'precipitation\|profile\|large\|scale\|i', 'kg/m2/s', 'BuPu']
	969	elif varn == 'prprofll' or varn == 'pr_lsc_l' or varn == 'LPR_LSC_L':
	970	varvals = ['prprofll', 'precipitation_profile_large_scale_l', 0., 1.e-3, \
	971	'precipitation\|profile\|large\|scale\|l', 'kg/m2/s', 'BuPu']
[192]	972	elif varn == 'pracc' or varn == 'ACRAINTOT':
	973	varvals = ['pracc', 'precipitation_amount', 0., 100., \
	974	'accumulated\|precipitation', 'kgm-2', 'BuPu']
	975	elif varn == 'prc' or varn == 'LPLUC' or varn == 'pluc' or varn == 'WRFprc':
	976	varvals = ['prc', 'convective_precipitation_flux', 0., 2.e-4, \
	977	'convective\|precipitation\|flux', 'kgm-2s-1', 'Blues']
[206]	978	elif varn == 'prci' or varn == 'pr_con_i' or varn == 'LPR_CON_I':
	979	varvals = ['prci', 'convective_ice_precipitation_flux', 0., 0.003, \
	980	'convective\|ice\|precipitation\|flux', 'kgm-2s-1', 'Purples']
	981	elif varn == 'prcl' or varn == 'pr_con_l' or varn == 'LPR_CON_L':
	982	varvals = ['prcl', 'convective_liquid_precipitation_flux', 0., 0.003, \
	983	'convective\|liquid\|precipitation\|flux', 'kgm-2s-1', 'Blues']
[220]	984	elif varn == 'pres' or varn == 'presnivs' or varn == 'pressure' or \
	985	varn == 'lpres' or varn == 'LPRES':
[192]	986	varvals = ['pres', 'air_pressure', 0., 103000., 'air\|pressure', 'Pa', \
	987	'Blues']
	988	elif varn == 'prls' or varn == 'WRFprls' or varn == 'LPLUL' or varn == 'plul':
	989	varvals = ['prls', 'large_scale_precipitation_flux', 0., 2.e-4, \
	990	'large\|scale\|precipitation\|flux', 'kgm-2s-1', 'Blues']
	991	elif varn == 'prsn' or varn == 'SNOW' or varn == 'snow' or varn == 'LSNOW':
	992	varvals = ['prsn', 'snowfall', 0., 1.e-4, 'snowfall\|flux', 'kgm-2s-1', 'BuPu']
	993	elif varn == 'prw' or varn == 'WRFprh':
	994	varvals = ['prw', 'atmosphere_water_vapor_content', 0., 10., \
	995	'water\|vapor"path', 'kgm-2', 'Blues']
	996	elif varn == 'ps' or varn == 'psfc' or varn =='PSFC' or varn == 'psol' or \
	997	varn == 'Surface Pressure':
	998	varvals=['ps', 'surface_air_pressure', 85000., 105400., 'surface\|pressure', \
	999	'hPa', 'cool']
	1000	elif varn == 'psl' or varn == 'mslp' or varn =='WRFmslp':
	1001	varvals=['psl', 'air_pressure_at_sea_level', 85000., 104000., \
	1002	'mean\|sea\|level\|pressure', 'Pa', 'Greens']
[249]	1003	elif varn == 'qth' or varn == 'q_th' or varn == 'LQ_TH':
[238]	1004	varvals = ['qth', 'thermal_plume_total_water_content', 0., 25., \
[192]	1005	'total\|water\|cotent\|in\|thermal\|plume', 'mm', 'YlOrRd']
	1006	elif varn == 'r' or varn == 'QVAPOR' or varn == 'ovap' or varn == 'LOVAP':
	1007	varvals = ['r', 'water_mixing_ratio', 0., 0.03, 'water\|mixing\|ratio', \
	1008	'kgkg-1', 'BuPu']
	1009	elif varn == 'rsds' or varn == 'SWdnSFC' or varn == 'SWdn at surface' or \
	1010	varn == 'SWDOWN':
	1011	varvals=['rsds', 'surface_downwelling_shortwave_flux_in_air', 0., 1200., \
	1012	'downward\|SW\|surface\|radiation', 'Wm-2' ,'Reds']
	1013	elif varn == 'rsdsacc':
	1014	varvals=['rsdsacc', 'accumulated_surface_downwelling_shortwave_flux_in_air', \
	1015	0., 1200., 'accumulated\|downward\|SW\|surface\|radiation', 'Wm-2' ,'Reds']
	1016	elif varn == 'rvor' or varn == 'WRFrvor':
	1017	varvals = ['rvor', 'air_relative_vorticity', -2.5E-3, 2.5E-3, \
	1018	'air\|relative\|vorticity', 's-1', 'seismic']
	1019	elif varn == 'rvors' or varn == 'WRFrvors':
	1020	varvals = ['rvors', 'surface_air_relative_vorticity', -2.5E-3, 2.5E-3, \
	1021	'surface\|air\|relative\|vorticity', 's-1', 'seismic']
[228]	1022	elif varn == 's' or varn == 'QSNOW':
	1023	varvals = ['s', 'snow_mixing_ratio', 0., 0.0003, 'snow\|mixing\|ratio', \
	1024	'kgkg-1', 'Purples']
[238]	1025	elif varn == 'stherm' or varn == 'LS_THERM':
[249]	1026	varvals = ['stherm', 'thermals_excess', 0., 0.8, 'thermals\|excess', 'K', \
[192]	1027	'Reds']
	1028	elif varn == 'ta' or varn == 'WRFt' or varn == 'temp' or varn == 'LTEMP' or \
	1029	varn == 'Air temperature':
	1030	varvals = ['ta', 'air_temperature', 195., 320., 'air\|temperature', 'K', \
	1031	'YlOrRd']
[249]	1032	elif varn == 'tah' or varn == 'theta' or varn == 'LTHETA':
	1033	varvals = ['tah', 'potential_air_temperature', 195., 320., \
	1034	'potential\|air\|temperature', 'K', 'YlOrRd']
[192]	1035	elif varn == 'tas' or varn == 'T2' or varn == 't2m' or varn == 'T2M' or \
	1036	varn == 'Temperature 2m':
	1037	varvals = ['tas', 'air_temperature', 240., 310., 'air\|temperature\|at\|2m', ' \
	1038	K', 'YlOrRd']
	1039	elif varn == 'tds' or varn == 'TH2':
	1040	varvals = ['tds', 'air_dew_point_temperature', 240., 310., \
	1041	'air\|dew\|point\|temperature\|at\|2m', 'K', 'YlGnBu']
[233]	1042	elif varn == 'tke' or varn == 'TKE' or varn == 'tke' or varn == 'LTKE':
	1043	varvals = ['tke', 'turbulent_kinetic_energy', 0., 0.003, \
	1044	'turbulent\|kinetic\|energy', 'm2/s2', 'Reds']
[192]	1045	elif varn == 'time'or varn == 'time_counter':
	1046	varvals = ['time', 'time', 0., 1000., 'time', \
	1047	'hours\|since\|1949/12/01\|00:00:00', 'Reds']
[233]	1048	elif varn == 'tmla' or varn == 's_pblt' or varn == 'LS_PBLT':
	1049	varvals = ['tmla', 'atmosphere_top_boundary_layer_temperature', 250., 330., \
	1050	'atmosphere\|top\|boundary\|layer\|temperature', 'K', 'Reds']
[192]	1051	elif varn == 'ua' or varn == 'vitu' or varn == 'U' or varn == 'Zonal wind' or \
	1052	varn == 'LVITU':
	1053	varvals = ['ua', 'eastward_wind', -30., 30., 'eastward\|wind', 'ms-1', \
	1054	'seismic']
	1055	elif varn == 'uas' or varn == 'u10m' or varn == 'U10' or varn =='Vent zonal 10m':
	1056	varvals = ['uas', 'eastward_wind', -30., 30., 'eastward\|2m\|wind', \
	1057	'ms-1', 'seismic']
	1058	elif varn == 'va' or varn == 'vitv' or varn == 'V' or varn == 'Meridional wind' \
	1059	or varn == 'LVITV':
	1060	varvals = ['va', 'northward_wind', -30., 30., 'northward\|wind', 'ms-1', \
	1061	'seismic']
	1062	elif varn == 'vas' or varn == 'v10m' or varn == 'V10' or \
	1063	varn =='Vent meridien 10m':
	1064	varvals = ['vas', 'northward_wind', -30., 30., 'northward\|2m\|wind', 'ms-1', \
	1065	'seismic']
[238]	1066	elif varn == 'wakedeltaq' or varn == 'wake_deltaq' or varn == 'lwake_deltaq' or \
	1067	varn == 'LWAKE_DELTAQ':
	1068	varvals = ['wakedeltaq', 'wake_delta_vapor', -0.003, 0.003, \
[233]	1069	'wake\|delta\|mixing\|ratio', '-', 'seismic']
[238]	1070	elif varn == 'wakedeltat' or varn == 'wake_deltat' or varn == 'lwake_deltat' or \
	1071	varn == 'LWAKE_DELTAT':
	1072	varvals = ['wakedeltat', 'wake_delta_temp', -0.003, 0.003, \
[233]	1073	'wake\|delta\|temperature', '-', 'seismic']
[238]	1074	elif varn == 'wakeh' or varn == 'wake_h' or varn == 'LWAKE_H':
	1075	varvals = ['wakeh', 'wake_height', 0., 1000., 'height\|of\|the\|wakes', 'm', \
[192]	1076	'YlOrRd']
[238]	1077	elif varn == 'wakeomg' or varn == 'wake_omg' or varn == 'lwake_omg' or \
	1078	varn == 'LWAKE_OMG':
	1079	varvals = ['wakeomg', 'wake_omega', 0., 3., 'wake\|omega', \
[233]	1080	'-', 'BuGn']
[238]	1081	elif varn == 'wakes' or varn == 'wake_s' or varn == 'LWAKE_S':
	1082	varvals = ['wakes', 'wake_area_fraction', 0., 0.5, 'wake\|spatial\|fraction', \
[192]	1083	'1', 'BuGn']
	1084	elif varn == 'wa' or varn == 'W' or varn == 'Vertical wind':
	1085	varvals = ['wa', 'upward_wind', -10., 10., 'upward\|wind', 'ms-1', \
	1086	'seismic']
	1087	elif varn == 'wap' or varn == 'vitw' or varn == 'LVITW':
	1088	varvals = ['wap', 'upward_wind', -3.e-10, 3.e-10, 'upward\|wind', 'mPa-1', \
	1089	'seismic']
	1090	elif varn == 'wss' or varn == 'SPDUV':
	1091	varvals = ['wss', 'air_velocity', 0., 30., 'surface\|horizontal\|wind\|speed', \
	1092	'ms-1', 'Reds']
	1093	elif varn == 'xtime' or varn == 'XTIME':
	1094	varvals = ['xtime', 'time', 0., 1.e5, 'time', \
	1095	'minutes\|since\|simulation\|start', 'Reds']
	1096	elif varn == 'x' or varn == 'X':
	1097	varvals = ['x', 'x', 0., 100., 'x', '-', 'Reds']
	1098	elif varn == 'y' or varn == 'Y':
	1099	varvals = ['y', 'y', 0., 100., 'y', '-', 'Blues']
	1100	elif varn == 'z' or varn == 'Z':
	1101	varvals = ['z', 'z', 0., 100., 'z', '-', 'Greens']
	1102	elif varn == 'zg' or varn == 'WRFght' or varn == 'Geopotential height' or \
	1103	varn == 'geop' or varn == 'LGEOP':
	1104	varvals = ['zg', 'geopotential_height', 0., 80000., 'geopotential\|height', \
	1105	'm2s-2', 'rainbow']
[238]	1106	elif varn == 'zmaxth' or varn == 'zmax_th' or varn == 'LZMAX_TH':
	1107	varvals = ['zmaxth', 'thermal_plume_height', 0., 4000., \
[192]	1108	'maximum\|thermals\|plume\|height', 'm', 'YlOrRd']
	1109	elif varn == 'zmla' or varn == 's_pblh' or varn == 'LS_PBLH':
	1110	varvals = ['zmla', 'atmosphere_boundary_layer_thickness', 0., 2500., \
	1111	'atmosphere\|boundary\|layer\|thickness', 'm', 'Blues']
	1112	else:
	1113	print errormsg
	1114	print ' ' + fname + ': variable ' + varn + ' not defined !!!'
	1115	quit(-1)
	1116
	1117	return varvals
	1118
	1119
	1120	####### ####### ####### ####### ####### ####### ####### ####### ####### #######
	1121
	1122	def check_colorBar(cbarn):
	1123	""" Check if the given colorbar exists in matplotlib
	1124	"""
	1125	fname = 'check_colorBar'
	1126
	1127	# Possible color bars
	1128	colorbars = ['binary', 'Blues', 'BuGn', 'BuPu', 'gist_yarg', 'GnBu', 'Greens', \
	1129	'Greys', 'Oranges', 'OrRd', 'PuBu', 'PuBuGn', 'PuRd', 'Purples', 'RdPu', \
	1130	'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd', 'afmhot', 'autumn', 'bone', \
	1131	'cool', 'copper', 'gist_gray', 'gist_heat', 'gray', 'hot', 'pink', 'spring', \
	1132	'summer', 'winter', 'BrBG', 'bwr', 'coolwarm', 'PiYG', 'PRGn', 'PuOr', 'RdBu', \
	1133	'RdGy', 'RdYlBu', 'RdYlGn', 'seismic', 'Accent', 'Dark2', 'hsv', 'Paired', \
	1134	'Pastel1', 'Pastel2', 'Set1', 'Set2', 'Set3', 'spectral', 'gist_earth', \
	1135	'gist_ncar', 'gist_rainbow', 'gist_stern', 'jet', 'brg', 'CMRmap', 'cubehelix',\
	1136	'gnuplot', 'gnuplot2', 'ocean', 'rainbow', 'terrain', 'flag', 'prism']
	1137
	1138	if not searchInlist(colorbars,cbarn):
	1139	print warnmsg
	1140	print ' ' + fname + ' color bar: "' + cbarn + '" does not exist !!'
	1141	print ' a standard one will be use instead !!'
	1142
	1143	return
	1144
	1145	def units_lunits(u):
	1146	""" Fucntion to provide LaTeX equivalences from a given units
	1147	u= units to transform
	1148	>>> units_lunits('kgkg-1')
	1149	'$kgkg^{-1}$'
	1150	"""
	1151	fname = 'units_lunits'
	1152
	1153	if u == 'h':
	1154	print fname + '_____________________________________________________________'
	1155	print units_lunits.__doc__
	1156	quit()
	1157
	1158	# Units which does not change
	1159	same = ['1', 'category', 'day', 'degrees East', 'degrees Nord', 'degrees North', \
	1160	'g', 'hour', 'hPa', 'K', 'Km', 'kg', 'km', 'm', 'minute', 'mm', 'month', 'Pa', \
[208]	1161	's', 'second', 'um', 'year', '-']
[192]	1162
	1163	if searchInlist(same,u):
	1164	lu = '$' + u + '$'
	1165	elif len(u.split(' ')) > 1 and u.split(' ')[1] == 'since':
	1166	uparts = u.split(' ')
	1167	ip=0
	1168	for up in uparts:
	1169	if ip == 0:
	1170	lu = '$' + up
	1171	else:
	1172	lu = lu + '\ ' + up
	1173	ip=ip+1
	1174	lu = lu + '$'
	1175	else:
[208]	1176	if u == '': lu='-'
	1177	elif u == 'C': lu='$^{\circ}C$'
[192]	1178	elif u == 'days': lu='$day$'
	1179	elif u == 'degrees_east': lu='$degrees\ East$'
	1180	elif u == 'degree_east': lu='$degrees\ East$'
	1181	elif u == 'degrees longitude': lu='$degrees\ East$'
	1182	elif u == 'degrees latitude': lu='$degrees\ North$'
	1183	elif u == 'degrees_north': lu='$degrees\ North$'
	1184	elif u == 'degree_north': lu='$degrees\ North$'
	1185	elif u == 'deg C': lu='$^{\circ}C$'
	1186	elif u == 'degC': lu='$^{\circ}C$'
	1187	elif u == 'deg K': lu='$K$'
	1188	elif u == 'degK': lu='$K$'
	1189	elif u == 'hours': lu='$hour$'
	1190	elif u == 'J/kg': lu='$Jkg^{-1}$'
	1191	elif u == 'Jkg-1': lu='$Jkg^{-1}$'
	1192	elif u == 'K/m': lu='$Km^{-1}$'
	1193	elif u == 'Km-1': lu='$Km^{-1}$'
	1194	elif u == 'K/s': lu='$Ks^{-1}$'
	1195	elif u == 'Ks-1': lu='$Ks^{-1}$'
[254]	1196	elif u == 'K s-1': lu='$Ks^{-1}$'
[192]	1197	elif u == 'kg/kg': lu='$kgkg^{-1}$'
	1198	elif u == 'kgkg-1': lu='$kgkg^{-1}$'
	1199	elif u == 'kg kg-1': lu='$kgkg^{-1}$'
	1200	elif u == '(kg/kg)/s': lu='$kgkg^{-1}s^{-1}$'
	1201	elif u == 'kgkg-1s-1': lu='$kgkg^{-1}s^{-1}$'
	1202	elif u == 'kg/m2': lu='$kgm^{-2}$'
	1203	elif u == 'kgm-2': lu='$kgm^{-2}$'
	1204	elif u == 'kg m-2': lu='$kgm^{-2}$'
	1205	elif u == 'Kg m-2': lu='$kgm^{-2}$'
	1206	elif u == 'kg/m2/s': lu='$kgm^{-2}s^{-1}$'
	1207	elif u == 'kg/(m2*s)': lu='$kgm^{-2}s^{-1}$'
	1208	elif u == 'kg/(s*m2)': lu='$kgm^{-2}s^{-1}$'
	1209	elif u == 'kgm-2s-1': lu='$kgm^{-2}s^{-1}$'
	1210	elif u == '1/m': lu='$m^{-1}$'
	1211	elif u == 'm-1': lu='$m^{-1}$'
	1212	elif u == 'm2/s': lu='$m2s^{-1}$'
	1213	elif u == 'm2s-1': lu='$m2s^{-1}$'
	1214	elif u == 'm2/s2': lu='$m2s^{-2}$'
	1215	elif u == 'm/s': lu='$ms^{-1}$'
	1216	elif u == 'mmh-3': lu='$mmh^{-3}$'
	1217	elif u == 'ms-1': lu='$ms^{-1}$'
	1218	elif u == 'm s-1': lu='$ms^{-1}$'
	1219	elif u == 'm/s2': lu='$ms^{-2}$'
	1220	elif u == 'ms-2': lu='$ms^{-2}$'
	1221	elif u == 'minutes': lu='$minute$'
	1222	elif u == 'Pa/s': lu='$Pas^{-1}$'
	1223	elif u == 'Pas-1': lu='$Pas^{-1}$'
	1224	elif u == 'W m-2': lu='$Wm^{-2}$'
	1225	elif u == 'Wm-2': lu='$Wm^{-2}$'
	1226	elif u == 'W/m2': lu='$Wm^{-2}$'
	1227	elif u == '1/s': lu='$s^{-1}$'
	1228	elif u == 's-1': lu='$s^{-1}$'
	1229	elif u == 'seconds': lu='$second$'
	1230	elif u == '%': lu='\%'
	1231	else:
	1232	print errormsg
	1233	print ' ' + fname + ': units "' + u + '" not ready!!!!'
	1234	quit(-1)
	1235
	1236	return lu
	1237
	1238	def ASCII_LaTeX(ln):
	1239	""" Function to transform from an ASCII line to LaTeX codification
	1240	>>> ASCII_LaTeX('Laboratoire de MÃ©tÃ©orologie Dynamique perÃ² HovmÃ¶ller')
	1241	Laboratoire de M\'et\'eorologie Dynamique per\`o Hovm\"oller
	1242	"""
	1243	fname='ASCII_LaTeX'
	1244
	1245	if ln == 'h':
	1246	print fname + '_____________________________________________________________'
	1247	print ASCII_LaTeX.__doc__
	1248	quit()
	1249
	1250	newln = ln.replace('\\', '\\textbackslash')
	1251
	1252	newln = newln.replace('Ã¡', "\\'a")
	1253	newln = newln.replace('Ã©', "\\'e")
	1254	newln = newln.replace('Ã', "\\'i")
	1255	newln = newln.replace('Ã³', "\\'o")
	1256	newln = newln.replace('Ãº', "\\'u")
	1257
	1258	newln = newln.replace('Ã ', "\\`a")
	1259	newln = newln.replace('Ãš', "\\`e")
	1260	newln = newln.replace('Ã¬', "\\`i")
	1261	newln = newln.replace('Ã²', "\\`o")
	1262	newln = newln.replace('Ã¹', "\\`u")
	1263
	1264	newln = newln.replace('Ã¢', "\\^a")
	1265	newln = newln.replace('Ãª', "\\^e")
	1266	newln = newln.replace('Ã®', "\\^i")
	1267	newln = newln.replace('ÃŽ', "\\^o")
	1268	newln = newln.replace('Ã»', "\\^u")
	1269
	1270	newln = newln.replace('Ã€', '\\"a')
	1271	newln = newln.replace('Ã«', '\\"e')
	1272	newln = newln.replace('Ã¯', '\\"i')
	1273	newln = newln.replace('Ã¶', '\\"o')
	1274	newln = newln.replace('ÃŒ', '\\"u')
	1275
	1276	newln = newln.replace('Ã§', '\c{c}')
	1277	newln = newln.replace('Ã±', '\~{n}')
	1278
	1279	newln = newln.replace('Ã', "\\'A")
	1280	newln = newln.replace('Ã', "\\'E")
	1281	newln = newln.replace('Ã', "\\'I")
	1282	newln = newln.replace('Ã', "\\'O")
	1283	newln = newln.replace('Ã', "\\'U")
	1284
	1285	newln = newln.replace('Ã', "\\`A")
	1286	newln = newln.replace('Ã', "\\`E")
	1287	newln = newln.replace('Ã', "\\`I")
	1288	newln = newln.replace('Ã', "\\`O")
	1289	newln = newln.replace('Ã', "\\`U")
	1290
	1291	newln = newln.replace('Ã', "\\^A")
	1292	newln = newln.replace('Ã', "\\^E")
	1293	newln = newln.replace('Ã', "\\^I")
	1294	newln = newln.replace('Ã', "\\^O")
	1295	newln = newln.replace('Ã', "\\^U")
	1296
	1297	newln = newln.replace('Ã', '\\"A')
	1298	newln = newln.replace('Ã', '\\"E')
	1299	newln = newln.replace('Ã', '\\"I')
	1300	newln = newln.replace('Ã', '\\"O')
	1301	newln = newln.replace('Ã', '\\"U')
	1302
	1303	newln = newln.replace('Ã', '\\c{C}')
	1304	newln = newln.replace('Ã', '\\~{N}')
	1305
	1306	newln = newln.replace('Â¡', '!`')
	1307	newln = newln.replace('Â¿', 'Â¿`')
	1308	newln = newln.replace('%', '\\%')
	1309	newln = newln.replace('#', '\\#')
	1310	newln = newln.replace('&', '\\&')
	1311	newln = newln.replace('$', '\\$')
	1312	newln = newln.replace('_', '\\_')
	1313	newln = newln.replace('Â·', '\\textperiodcentered')
	1314	newln = newln.replace('<', '$<$')
	1315	newln = newln.replace('>', '$>$')
	1316	newln = newln.replace('ï', '*')
	1317	# newln = newln.replace('Âº', '$^{\\circ}$')
	1318	newln = newln.replace('Âª', '$^{a}$')
	1319	newln = newln.replace('Âº', '$^{o}$')
	1320	newln = newln.replace('Â°', '$^{\\circ}$')
	1321	newln = newln.replace('\n', '\\\\\n')
	1322	newln = newln.replace('\t', '\\medskip')
	1323
	1324	return newln
	1325
	1326	def pretty_int(minv,maxv,Nint):
	1327	""" Function to plot nice intervals
	1328	minv= minimum value
	1329	maxv= maximum value
	1330	Nint= number of intervals
	1331	>>> pretty_int(23.50,67.21,5)
	1332	[ 25. 30. 35. 40. 45. 50. 55. 60. 65.]
	1333	>>> pretty_int(-23.50,67.21,15)
	1334	[ 0. 20. 40. 60.]
	1335	pretty_int(14.75,25.25,5)
	1336	[ 16. 18. 20. 22. 24.]
	1337	"""
	1338	fname = 'pretty_int'
	1339	nice_int = [1,2,5]
	1340
	1341	# print 'minv: ',minv,'maxv:',maxv,'Nint:',Nint
	1342
	1343	interval = np.abs(maxv - minv)
	1344
	1345	potinterval = np.log10(interval)
	1346
	1347	Ipotint = int(potinterval)
	1348	intvalue = np.float(interval / np.float(Nint))
	1349
	1350	# new
	1351	potinterval = np.log10(intvalue)
	1352	Ipotint = int(potinterval)
	1353
	1354	# print 'interval:', interval, 'intavlue:', intvalue, 'potinterval:', potinterval, \
	1355	# 'Ipotint:', Ipotint, 'intvalue:', intvalue
	1356
	1357	mindist = 10.e15
	1358	for inice in nice_int:
	1359	# print inice,':',inice10.Ipotint,np.abs(inice10.**Ipotint - intvalue),mindist
	1360	if np.abs(inice10.*Ipotint - intvalue) < mindist:
	1361	mindist = np.abs(inice10.*Ipotint - intvalue)
	1362	closestint = inice
	1363
	1364	Ibeg = int(minv / (closestint10.*Ipotint))
	1365
	1366	values = []
	1367	val = closestint(Ibeg)10.**(Ipotint)
	1368
	1369	# print 'closestint:',closestint,'Ibeg:',Ibeg,'val:',val
	1370
	1371	while val < maxv:
	1372	values.append(val)
	1373	val = val + closestint10.*Ipotint
	1374
	1375	return np.array(values, dtype=np.float)
	1376
	1377	def DegGradSec_deg(grad,deg,sec):
	1378	""" Function to transform from a coordinate in grad deg sec to degrees (decimal)
	1379	>>> DegGradSec_deg(39.,49.,26.)
	1380	39.8238888889
	1381	"""
	1382	fname = 'DegGradSec_deg'
	1383
	1384	if grad == 'h':
	1385	print fname + '_____________________________________________________________'
	1386	print DegGradSec_deg.__doc__
	1387	quit()
	1388
	1389	deg = grad + deg/60. + sec/3600.
	1390
	1391	return deg
	1392
	1393	def intT2dt(intT,tu):
	1394	""" Function to provide an 'timedelta' object from a given interval value
	1395	intT= interval value
	1396	tu= interval units, [tu]= 'd': day, 'w': week, 'h': hour, 'i': minute, 's': second,
	1397	'l': milisecond
	1398
	1399	>>> intT2dt(3.5,'s')
	1400	0:00:03.500000
	1401
	1402	>>> intT2dt(3.5,'w')
	1403	24 days, 12:00:00
	1404	"""
	1405	import datetime as dt
	1406
	1407	fname = 'intT2dt'
	1408
	1409	if tu == 'w':
	1410	dtv = dt.timedelta(weeks=np.float(intT))
	1411	elif tu == 'd':
	1412	dtv = dt.timedelta(days=np.float(intT))
	1413	elif tu == 'h':
	1414	dtv = dt.timedelta(hours=np.float(intT))
	1415	elif tu == 'i':
	1416	dtv = dt.timedelta(minutes=np.float(intT))
	1417	elif tu == 's':
	1418	dtv = dt.timedelta(seconds=np.float(intT))
	1419	elif tu == 'l':
	1420	dtv = dt.timedelta(milliseconds=np.float(intT))
	1421	else:
	1422	print errormsg
	1423	print ' ' + fname + ': time units "' + tu + '" not ready!!!!'
	1424	quit(-1)
	1425
	1426	return dtv
	1427
	1428	def lonlat_values(mapfile,lonvn,latvn):
	1429	""" Function to obtain the lon/lat matrices from a given netCDF file
	1430	lonlat_values(mapfile,lonvn,latvn)
	1431	[mapfile]= netCDF file name
	1432	[lonvn]= variable name with the longitudes
	1433	[latvn]= variable name with the latitudes
	1434	"""
	1435
	1436	fname = 'lonlat_values'
	1437
	1438	if mapfile == 'h':
	1439	print fname + '_____________________________________________________________'
	1440	print lonlat_values.__doc__
	1441	quit()
	1442
	1443	if not os.path.isfile(mapfile):
	1444	print errormsg
	1445	print ' ' + fname + ": map file '" + mapfile + "' does not exist !!"
	1446	quit(-1)
	1447
	1448	ncobj = NetCDFFile(mapfile, 'r')
	1449	lonobj = ncobj.variables[lonvn]
	1450	latobj = ncobj.variables[latvn]
	1451
	1452	if len(lonobj.shape) == 3:
	1453	lonv = lonobj[0,:,:]
	1454	latv = latobj[0,:,:]
	1455	elif len(lonobj.shape) == 2:
	1456	lonv = lonobj[:,:]
	1457	latv = latobj[:,:]
	1458	elif len(lonobj.shape) == 1:
	1459	lon0 = lonobj[:]
	1460	lat0 = latobj[:]
	1461	lonv = np.zeros( (len(lat0),len(lon0)), dtype=np.float )
	1462	latv = np.zeros( (len(lat0),len(lon0)), dtype=np.float )
	1463	for iy in range(len(lat0)):
	1464	lonv[iy,:] = lon0
	1465	for ix in range(len(lon0)):
	1466	latv[:,ix] = lat0
	1467	else:
	1468	print errormsg
	1469	print ' ' + fname + ': lon/lat variables shape:',lonobj.shape,'not ready!!'
	1470	quit(-1)
	1471
	1472	return lonv, latv
	1473
	1474	def date_CFtime(ind,refd,tunits):
	1475	""" Function to transform from a given date object a CF-convention time
	1476	ind= date object to transform
	1477	refd= reference date
	1478	tunits= units for time
	1479	>>> date_CFtime(dt.datetime(1976,02,17,08,30,00), dt.datetime(1949,12,01,00,00,00), 'seconds')
	1480	827224200.0
	1481	"""
	1482	import datetime as dt
	1483
	1484	fname = 'date_CFtime'
	1485
	1486	dt = ind - refd
	1487
	1488	if tunits == 'weeks':
	1489	value = dt.days/7. + dt.seconds/(3600.24.7.)
	1490	elif tunits == 'days':
	1491	value = dt.days + dt.seconds/(3600.*24.)
	1492	elif tunits == 'hours':
	1493	value = dt.days*24. + dt.seconds/(3600.)
	1494	elif tunits == 'minutes':
	1495	value = dt.days24.60. + dt.seconds/(60.)
	1496	elif tunits == 'seconds':
	1497	value = dt.days24.3600. + dt.seconds
	1498	elif tunits == 'milliseconds':
	1499	value = dt.days24.3600.1000. + dt.seconds1000.
	1500	else:
	1501	print errormsg
	1502	print ' ' + fname + ': reference time units "' + trefu + '" not ready!!!!'
	1503	quit(-1)
	1504
	1505	return value
	1506
	1507	def pot_values(values, uvals):
	1508	""" Function to modify a seies of values by their potency of 10
	1509	pot_values(values, uvals)
	1510	values= values to modify
	1511	uvals= units of the values
	1512	>>> vals = np.sin(np.arange(20)np.pi/5.+0.01)10.e-5
	1513	>>> pot_values(vals,'ms-1')
	1514	(array([ 0.00000000e+00, 5.87785252e-01, 9.51056516e-01,
	1515	9.51056516e-01, 5.87785252e-01, 1.22464680e-16,
	1516	-5.87785252e-01, -9.51056516e-01, -9.51056516e-01,
	1517	-5.87785252e-01, -2.44929360e-16, 5.87785252e-01,
	1518	9.51056516e-01, 9.51056516e-01, 5.87785252e-01,
	1519	3.67394040e-16, -5.87785252e-01, -9.51056516e-01,
	1520	-9.51056516e-01, -5.87785252e-01]), -4, 'x10e-4 ms-1', 'x10e-4')
	1521	"""
	1522
	1523	fname = 'pot_values'
	1524
	1525	if np.min(values) != 0.:
	1526	potmin = int( np.log10( np.abs(np.min(values)) ) )
	1527	else:
	1528	potmin = 0
	1529
	1530	if np.max(values) != 0.:
	1531	potmax = int( np.log10( np.abs(np.max(values)) ) )
	1532	else:
	1533	potmax = 0
	1534
	1535	if potmin * potmax > 9:
	1536	potval = -np.min([np.abs(potmin), np.abs(potmax)]) * np.abs(potmin) / potmin
	1537
	1538	newvalues = values10.*potval
	1539	potvalue = - potval
	1540	potS = 'x10e' + str(potvalue)
	1541	newunits = potS + ' ' + uvals
	1542	else:
	1543	newvalues = values
	1544	potvalue = None
	1545	potS = ''
	1546	newunits = uvals
	1547
	1548	return newvalues, potvalue, newunits, potS
	1549
	1550	def CFtimes_plot(timev,units,kind,tfmt):
	1551	""" Function to provide a list of string values from a CF time values in order
	1552	to use them in a plot, according to the series of characteristics.
	1553	String outputs will be suited to the 'human-like' output
	1554	timev= time values (real values)
	1555	units= units string according to CF conventions ([tunits] since
	1556	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]])
	1557	kind= kind of output
	1558	'Nval': according to a given number of values as 'Nval',[Nval]
	1559	'exct': according to an exact time unit as 'exct',[tunit];
	1560	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
	1561	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
	1562	'l': milisecond
	1563	tfmt= desired format
	1564	>>> CFtimes_plot(np.arange(100)*1.,'hours since 1979-12-01 00:00:00', 'Nval,5',"%Y/%m/%d %H:%M:%S")
	1565	0.0 1979/12/01 00:00:00
	1566	24.75 1979/12/02 00:45:00
	1567	49.5 1979/12/03 01:30:00
	1568	74.25 1979/12/04 02:15:00
	1569	99.0 1979/12/05 03:00:00
	1570	>>> CFtimes_plot(np.arange(100)*1.,'hours since 1979-12-01 00:00:00', 'exct,2,d',"%Y/%m/%d %H:%M:%S")
	1571	0.0 1979/12/01 00:00:00
	1572	48.0 1979/12/03 00:00:00
	1573	96.0 1979/12/05 00:00:00
	1574	144.0 1979/12/07 00:00:00
	1575	"""
	1576	import datetime as dt
	1577
	1578	# Seconds between 0001 and 1901 Jan - 01
	1579	secs0001_1901=59958144000.
	1580
	1581	fname = 'CFtimes_plot'
	1582
	1583	if timev == 'h':
	1584	print fname + '_____________________________________________________________'
	1585	print CFtimes_plot.__doc__
	1586	quit()
	1587
	1588	# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
	1589	##
	1590	trefT = units.find(':')
	1591	txtunits = units.split(' ')
	1592	Ntxtunits = len(txtunits)
	1593
	1594	if Ntxtunits == 3:
	1595	Srefdate = txtunits[Ntxtunits - 1]
	1596	else:
	1597	Srefdate = txtunits[Ntxtunits - 2]
	1598
	1599	if not trefT == -1:
	1600	# print ' ' + fname + ': refdate with time!'
	1601	if Ntxtunits == 3:
	1602	refdate = datetimeStr_datetime(Srefdate)
	1603	else:
	1604	refdate = datetimeStr_datetime(Srefdate + '_' + txtunits[Ntxtunits - 1])
	1605	else:
	1606	refdate = datetimeStr_datetime(Srefdate + '_00:00:00')
	1607
	1608	trefunits=units.split(' ')[0]
	1609	if trefunits == 'weeks':
	1610	trefu = 'w'
	1611	elif trefunits == 'days':
	1612	trefu = 'd'
	1613	elif trefunits == 'hours':
	1614	trefu = 'h'
	1615	elif trefunits == 'minutes':
	1616	trefu = 'm'
	1617	elif trefunits == 'seconds':
	1618	trefu = 's'
	1619	elif trefunits == 'milliseconds':
	1620	trefu = 'l'
	1621	else:
	1622	print errormsg
	1623	print ' ' + fname + ': reference time units "' + trefu + '" not ready!!!!'
	1624	quit(-1)
	1625
	1626	okind=kind.split(',')[0]
	1627	dtv = len(timev)
	1628
	1629	if refdate.year == 1:
	1630	print warnmsg
	1631	print ' ' + fname + ': changing reference date: ',refdate, \
	1632	'to 1901-01-01_00:00:00 !!!'
	1633	refdate = datetimeStr_datetime('1901-01-01_00:00:00')
	1634	if trefu == 'w': timev = timev - secs0001_1901/(7.24.3600.)
	1635	if trefu == 'd': timev = timev - secs0001_1901/(24.*3600.)
	1636	if trefu == 'h': timev = timev - secs0001_1901/(3600.)
	1637	if trefu == 'm': timev = timev - secs0001_1901/(60.)
	1638	if trefu == 's': timev = timev - secs0001_1901
	1639	if trefu == 'l': timev = timev - secs0001_1901*1000.
	1640
	1641	firstT = timev[0]
	1642	lastT = timev[dtv-1]
	1643
	1644	# First and last times as datetime objects
	1645	firstTdt = timeref_datetime(refdate, firstT, trefunits)
	1646	lastTdt = timeref_datetime(refdate, lastT, trefunits)
	1647
	1648	# First and last times as [year, mon, day, hour, minut, second] vectors
	1649	firstTvec = np.zeros((6), dtype= np.float)
	1650	lastTvec = np.zeros((6), dtype= np.float)
	1651	chTvec = np.zeros((6), dtype= bool)
	1652
	1653	firstTvec = np.array([firstTdt.year, firstTdt.month, firstTdt.day, firstTdt.hour,\
	1654	firstTdt.minute, firstTdt.second])
	1655	lastTvec = np.array([lastTdt.year, lastTdt.month, lastTdt.day, lastTdt.hour, \
	1656	lastTdt.minute, lastTdt.second])
	1657
	1658	chdate= lastTvec - firstTvec
	1659	chTvec = np.where (chdate != 0., True, False)
	1660
	1661	timeout = []
	1662	if okind == 'Nval':
	1663	nvalues = int(kind.split(',')[1])
	1664	intervT = (lastT - firstT)/(nvalues-1)
	1665	dtintervT = intT2dt(intervT, trefu)
	1666
	1667	for it in range(nvalues):
	1668	timeout.append(firstTdt + dtintervT*it)
	1669	elif okind == 'exct':
	1670	Nunits = int(kind.split(',')[1])
	1671	tu = kind.split(',')[2]
	1672
	1673	# Generic incremental dt [seconds] according to all the possibilities ['c', 'y', 'm',
	1674	# 'w', 'd', 'h', 'i', 's', 'l'], some of them approximated (because they are not
	1675	# already necessary!)
	1676	basedt = np.zeros((9), dtype=np.float)
	1677	basedt[0] = (365.100. + 25.)24.*3600.
	1678	basedt[1] = 365.24.3600.
	1679	basedt[2] = 31.24.3600.
	1680	basedt[3] = 7.24.3600.
	1681	basedt[4] = 24.*3600.
	1682	basedt[5] = 3600.
	1683	basedt[6] = 60.
	1684	basedt[7] = 1.
	1685	basedt[8] = 1000.
	1686
	1687	# Increment according to the units of the CF dates
	1688	if trefunits == 'weeks':
	1689	basedt = basedt/(7.24.3600.)
	1690	elif trefunits == 'days':
	1691	basedt = basedt/(24.*3600.)
	1692	elif trefunits == 'hours':
	1693	basedt = basedt/(3600.)
	1694	elif trefunits == 'minutes':
	1695	basedt = basedt/(60.)
	1696	elif trefunits == 'seconds':
	1697	basedt = basedt
	1698	elif trefunits == 'milliseconds':
	1699	basedt = basedt*1000.
	1700
	1701	if tu == 'c':
	1702	ti = firstTvec[0]
	1703	tf = lastTvec[0]
	1704	centi = firstTvec[0] / 100
	1705
	1706	for it in range((tf - ti)/(Nunits*100) + 1):
	1707	timeout.append(dt.datetime(centi+itNunits100, 1, 1, 0, 0, 0))
	1708	elif tu == 'y':
	1709	ti = firstTvec[0]
	1710	tf = lastTvec[0]
	1711	yeari = firstTvec[0]
	1712
	1713	for it in range((tf - ti)/(Nunits) + 1):
	1714	timeout.append(dt.datetime(yeari+it*Nunits, 1, 1, 0, 0, 0))
	1715	elif tu == 'm':
	1716	ti = firstTvec[1]
	1717	tf = lastTvec[1]
	1718	yr = firstTvec[0]
	1719	mon = firstTvec[1]
	1720
	1721	for it in range((tf - ti)/(Nunits) + 1):
	1722	mon = mon+it*Nunits
	1723	if mon > 12:
	1724	yr = yr + 1
	1725	mon = 1
	1726
	1727	timeout.append(dt.datetime(yr, mon, 1, 0, 0, 0))
	1728	elif tu == 'w':
	1729	datev=firstTdt
	1730	it=0
	1731	while datev <= lastTdt:
	1732	datev = firstTdt + dt.timedelta(days=7Nunitsit)
	1733	timeout.append(datev)
	1734	it = it + 1
	1735	elif tu == 'd':
	1736	# datev=firstTdt
	1737	yr = firstTvec[0]
	1738	mon = firstTvec[1]
	1739	day = firstTvec[2]
	1740
	1741	if np.sum(firstTvec[2:5]) > 0:
	1742	firstTdt = dt.datetime(yr, mon, day+1, 0, 0, 0)
	1743	datev = dt.datetime(yr, mon, day+1, 0, 0, 0)
	1744	else:
	1745	firstTdt = dt.datetime(yr, mon, day, 0, 0, 0)
	1746	datev = dt.datetime(yr, mon, day, 0, 0, 0)
	1747
	1748	it=0
	1749	while datev <= lastTdt:
	1750	datev = firstTdt + dt.timedelta(days=Nunits*it)
	1751	timeout.append(datev)
	1752	it = it + 1
	1753	elif tu == 'h':
	1754	datev=firstTdt
	1755	yr = firstTvec[0]
	1756	mon = firstTvec[1]
	1757	day = firstTvec[2]
	1758	hour = firstTvec[3]
	1759
	1760	if np.sum(firstTvec[4:5]) > 0 or np.mod(hour,Nunits) != 0:
	1761	tadvance = 2*Nunits
	1762	if tadvance >= 24:
	1763	firstTdt = dt.datetime(yr, mon, day+1, 0, 0, 0)
	1764	datev = dt.datetime(yr, mon, day+1, 0, 0, 0)
	1765	else:
	1766	firstTdt = dt.datetime(yr, mon, day, Nunits, 0, 0)
	1767	datev = dt.datetime(yr, mon, day, Nunits, 0, 0)
	1768	else:
	1769	firstTdt = dt.datetime(yr, mon, day, hour, 0, 0)
	1770	datev = dt.datetime(yr, mon, day, hour, 0, 0)
	1771
	1772	it=0
	1773	while datev <= lastTdt:
	1774	datev = firstTdt + dt.timedelta(seconds=Nunits3600it)
	1775	timeout.append(datev)
	1776	it = it + 1
	1777	elif tu == 'i':
	1778	datev=firstTdt
	1779	yr = firstTvec[0]
	1780	mon = firstTvec[1]
	1781	day = firstTvec[2]
	1782	hour = firstTvec[3]
	1783	minu = firstTvec[4]
	1784
	1785	if firstTvec[5] > 0 or np.mod(minu,Nunits) != 0:
	1786	tadvance = 2*Nunits
	1787	if tadvance >= 60:
	1788	firstTdt = dt.datetime(yr, mon, day, hour, 0, 0)
	1789	datev = dt.datetime(yr, mon, day, hour, 0, 0)
	1790	else:
	1791	firstTdt = dt.datetime(yr, mon, day, hour, Nunits, 0)
	1792	datev = dt.datetime(yr, mon, day, hour, Nunits, 0)
	1793	else:
	1794	firstTdt = dt.datetime(yr, mon, day, hour, minu, 0)
	1795	datev = dt.datetime(yr, mon, day, hour, minu, 0)
	1796	it=0
	1797	while datev <= lastTdt:
	1798	datev = firstTdt + dt.timedelta(seconds=Nunits60it)
	1799	timeout.append(datev)
	1800	it = it + 1
	1801	elif tu == 's':
	1802	datev=firstTdt
	1803	it=0
	1804	while datev <= lastTdt:
	1805	datev = firstTdt + dt.timedelta(seconds=Nunits)
	1806	timeout.append(datev)
	1807	it = it + 1
	1808	elif tu == 'l':
	1809	datev=firstTdt
	1810	it=0
	1811	while datev <= lastTdt:
	1812	datev = firstTdt + dt.timedelta(seconds=Nunits*it/1000.)
	1813	timeout.append(datev)
	1814	it = it + 1
	1815	else:
	1816	print errormsg
	1817	print ' ' + fname + ': exact units "' + tu + '" not ready!!!!!'
	1818	quit(-1)
	1819
	1820	else:
	1821	print errormsg
	1822	print ' ' + fname + ': output kind "' + okind + '" not ready!!!!'
	1823	quit(-1)
	1824
	1825	dtout = len(timeout)
	1826
	1827	timeoutS = []
	1828	timeoutv = np.zeros((dtout), dtype=np.float)
	1829
	1830	for it in range(dtout):
	1831	timeoutS.append(timeout[it].strftime(tfmt))
	1832	timeoutv[it] = date_CFtime(timeout[it], refdate, trefunits)
	1833
	1834	# print it,':',timeoutv[it], timeoutS[it]
	1835
	1836	return timeoutv, timeoutS
	1837
	1838	def color_lines(Nlines):
	1839	""" Function to provide a color list to plot lines
	1840	color_lines(Nlines)
	1841	Nlines= number of lines
	1842	"""
	1843
	1844	fname = 'color_lines'
	1845
	1846	colors = ['r', 'b', 'g', 'p', 'g']
	1847
	1848	colorv = []
	1849
	1850	colorv.append('k')
	1851	for icol in range(Nlines):
	1852	colorv.append(colors[icol])
	1853
	1854
	1855	return colorv
	1856
	1857	def output_kind(kindf, namef, close):
	1858	""" Function to generate the output of the figure
	1859	kindf= kind of the output
	1860	null: show in screen
	1861	[jpg/pdf/png/ps]: standard output types
	1862	namef= name of the figure (without extension)
	1863	close= if the graph has to be close or not [True/False]
	1864	"""
	1865	fname = 'output_kind'
	1866
	1867	if kindf == 'h':
	1868	print fname + '_____________________________________________________________'
	1869	print output_kind.__doc__
	1870	quit()
	1871
	1872	if kindf == 'null':
	1873	print 'showing figure...'
	1874	plt.show()
	1875	elif kindf == 'gif':
	1876	plt.savefig(namef + ".gif")
	1877	if close: print "Successfully generation of figure '" + namef + ".jpg' !!!"
	1878	elif kindf == 'jpg':
	1879	plt.savefig(namef + ".jpg")
	1880	if close: print "Successfully generation of figure '" + namef + ".jpg' !!!"
	1881	elif kindf == 'pdf':
	1882	plt.savefig(namef + ".pdf")
	1883	if close: print "Successfully generation of figure '" + namef + ".pdf' !!!"
	1884	elif kindf == 'png':
	1885	plt.savefig(namef + ".png")
	1886	if close: print "Successfully generation of figure '" + namef + ".png' !!!"
	1887	elif kindf == 'ps':
	1888	plt.savefig(namef + ".ps")
	1889	if close: print "Successfully generation of figure '" + namef + ".ps' !!!"
	1890	else:
	1891	print errormsg
	1892	print ' ' + fname + ' output format: "' + kindf + '" not ready !!'
	1893	print errormsg
	1894	quit(-1)
	1895
	1896	if close:
	1897	plt.close()
	1898
	1899	return
	1900
	1901	def check_arguments(funcname,Nargs,args,char,expectargs):
	1902	""" Function to check the number of arguments if they are coincident
	1903	check_arguments(funcname,Nargs,args,char)
	1904	funcname= name of the function/program to check
	1905	Nargs= theoretical number of arguments
	1906	args= passed arguments
	1907	char= character used to split the arguments
	1908	"""
	1909
	1910	fname = 'check_arguments'
	1911
	1912	Nvals = len(args.split(char))
	1913	if Nvals != Nargs:
	1914	print errormsg
	1915	print ' ' + fname + ': wrong number of arguments:',Nvals," passed to '", \
	1916	funcname, "' which requires:",Nargs,'!!'
	1917	print ' given arguments:',args.split(char)
	1918	print ' expected arguments:',expectargs
	1919	quit(-1)
	1920
	1921	return
	1922
	1923	def Str_Bool(val):
	1924	""" Function to transform from a String value to a boolean one
	1925	>>> Str_Bool('True')
	1926	True
	1927	>>> Str_Bool('0')
	1928	False
	1929	>>> Str_Bool('no')
	1930	False
	1931	"""
	1932
	1933	fname = 'Str_Bool'
	1934
	1935	if val == 'True' or val == '1' or val == 'yes':
	1936	boolv = True
	1937	elif val == 'False' or val == '0' or val== 'no':
	1938	boolv = False
	1939	else:
	1940	print errormsg
	1941	print ' ' + fname + ": value '" + val + "' not ready!!"
	1942	quit(-1)
	1943
	1944	return boolv
	1945
	1946	def plot_TimeSeries(valtimes, vunits, tunits, hfileout, vtit, ttit, tkind, tformat, \
	1947	tit, linesn, lloc, kfig):
	1948	""" Function to draw time-series
	1949	valtimes= list of arrays to plot [vals1[1values, 1times], [...,valsM[Mvals,Mtimes]])
	1950	vunits= units of the values
	1951	tunits= units of the times
	1952	hfileout= header of the output figure. Final name: [hfileout]_[vtit].[kfig]
	1953	vtit= variable title to be used in the graph
	1954	ttit= time title to be used in the graph
	1955	tkind= kind of time values to appear in the x-axis
	1956	'Nval': according to a given number of values as 'Nval',[Nval]
	1957	'exct': according to an exact time unit as 'exct',[tunit];
	1958	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
	1959	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
	1960	'l': milisecond
	1961	tformat= desired format of times
	1962	tit= title of the graph
	1963	linesn= list of values fot the legend
	1964	lloc= location of the legend (-1, autmoatic)
	1965	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
	1966	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
	1967	9: 'upper center', 10: 'center'
	1968	kfig= type of figure: jpg, png, pds, ps
	1969	"""
	1970	fname = 'plot_TimeSeries'
	1971
	1972	if valtimes == 'h':
	1973	print fname + '_____________________________________________________________'
	1974	print plot_TimeSeries.__doc__
	1975	quit()
	1976
	1977
	1978	# Canging line kinds every 7 lines (end of standard colors)
	1979	linekinds=['.-','x-','o-']
	1980
	1981	Nlines = len(valtimes)
	1982
	1983	Nvalues = []
	1984	Ntimes = []
	1985
	1986	for il in range(Nlines):
	1987	array = valtimes[il]
	1988
	1989	if Nlines == 1:
	1990	print warnmsg
	1991	print ' ' + fname + ': drawing only one line!'
	1992
	1993	Nvalues.append(array.shape[1])
	1994	Ntimes.append(array.shape[0])
	1995	tmin = np.min(array[1])
	1996	tmax = np.max(array[1])
	1997	vmin = np.min(array[0])
	1998	vmax = np.max(array[0])
	1999	else:
	2000	Nvalues.append(array.shape[1])
	2001	Ntimes.append(array.shape[0])
	2002	tmin = np.min(array[1,:])
	2003	tmax = np.max(array[1,:])
	2004	vmin = np.min(array[0,:])
	2005	vmax = np.max(array[0,:])
	2006
	2007	if il == 0:
	2008	xmin = tmin
	2009	xmax = tmax
	2010	ymin = vmin
	2011	ymax = vmax
	2012	else:
	2013	if tmin < xmin: xmin = tmin
	2014	if tmax > xmax: xmax = tmax
	2015	if vmin < ymin: ymin = vmin
	2016	if vmax > ymax: ymax = vmax
	2017
	2018	dx = np.max(Ntimes)
	2019	dy = np.min(Nvalues)
	2020
	2021	plt.rc('text', usetex=True)
	2022
	2023	print vtit
	2024	if vtit == 'ps':
	2025	plt.ylim(98000.,ymax)
	2026	else:
	2027	plt.ylim(ymin,ymax)
	2028
	2029	plt.xlim(xmin,xmax)
	2030	# print 'x lim:',xmin,xmax
	2031	# print 'y lim:',ymin,ymax
	2032
	2033	N7lines=0
	2034	for il in range(Nlines):
	2035	array = valtimes[il]
	2036	if vtit == 'ps':
	2037	array[0,:] = np.where(array[0,:] < 98000., None, array[0,:])
	2038	plt.plot(array[1,:],array[0,:], linekinds[N7lines], label= linesn[il])
	2039	if il == 6: N7lines = N7lines + 1
	2040
	2041	timevals = np.arange(xmin,xmax)*1.
	2042
	2043	tpos, tlabels = CFtimes_plot(timevals, tunits, tkind, tformat)
	2044
	2045	if len(tpos) > 10:
	2046	print warnmsg
	2047	print ' ' + fname + ': with "' + tkind + '" there are', len(tpos), 'xticks !'
	2048
	2049	plt.xticks(tpos, tlabels)
	2050	# plt.Axes.set_xticklabels(tlabels)
	2051
	2052	plt.legend(loc=lloc)
	2053	plt.xlabel(ttit)
	2054	plt.ylabel(vtit + " (" + vunits + ")")
	2055	plt.title(tit)
	2056
	2057	figname = hfileout + '_' + vtit
	2058
	2059	output_kind(kfig, figname, True)
	2060
	2061	return
	2062
	2063	#Nt = 50
	2064	#Nlines = 3
	2065
	2066	#vtvalsv = []
	2067
	2068	#valsv = np.zeros((2,Nt), dtype=np.float)
	2069	## First
	2070	#valsv[0,:] = np.arange(Nt)
	2071	#valsv[1,:] = np.arange(Nt)*180.
	2072	#vtvalsv.append(valsv)
	2073	#del(valsv)
	2074
	2075	#valsv = np.zeros((2,Nt/2), dtype=np.float)
	2076	## Second
	2077	#valsv[0,:] = np.arange(Nt/2)
	2078	#valsv[1,:] = np.arange(Nt/2)180.2.
	2079	#vtvalsv.append(valsv)
	2080	#del(valsv)
	2081
	2082	#valsv = np.zeros((2,Nt/4), dtype=np.float)
	2083	## Third
	2084	#valsv[0,:] = np.arange(Nt/4)
	2085	#valsv[1,:] = np.arange(Nt/4)180.4.
	2086	#vtvalsv.append(valsv)
	2087	#del(valsv)
	2088
	2089	#varu='mm'
	2090	#timeu='seconds'
	2091
	2092	#title='test'
	2093	#linesname = ['line 1', 'line 2', 'line 3']
	2094
	2095	#plot_TimeSeries(vtvalsv, units_lunits(varu), timeu, 'test', 'vartest', 'time', title, linesname, 'png')
	2096	#quit()
	2097
	2098	def plot_points(xval, yval, ifile, vtit, kfig, mapv):
	2099	""" plotting points
	2100	[x/yval]: x,y values to plot
	2101	vn,vm= minmum and maximum values to plot
	2102	unit= units of the variable
	2103	ifile= name of the input file
	2104	vtit= title of the variable
	2105	kfig= kind of figure (jpg, pdf, png)
	2106	mapv= map characteristics: [proj],[res]
	2107	see full documentation: http://matplotlib.org/basemap/
	2108	[proj]: projection
	2109	* 'cyl', cilindric
	2110	[res]: resolution:
	2111	* 'c', crude
	2112	* 'l', low
	2113	* 'i', intermediate
	2114	* 'h', high
	2115	* 'f', full
	2116	"""
	2117	fname = 'plot_points'
	2118
	2119	dx=xval.shape[0]
	2120	dy=yval.shape[0]
	2121
	2122	plt.rc('text', usetex=True)
	2123
	2124	if not mapv is None:
	2125	lon00 = np.where(xval[:] < 0., 360. + olon[:], olon[:])
	2126	lat00 = yval[:]
	2127	lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
	2128	lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
	2129
	2130	for iy in range(len(lat00)):
	2131	lon0[iy,:] = lon00
	2132	for ix in range(len(lon00)):
	2133	lat0[:,ix] = lat00
	2134
	2135	map_proj=mapv.split(',')[0]
	2136	map_res=mapv.split(',')[1]
	2137
	2138	nlon = lon0[0,0]
	2139	xlon = lon0[dy-1,dx-1]
	2140	nlat = lat0[0,0]
	2141	xlat = lat0[dy-1,dx-1]
	2142
	2143	lon2 = lon0[dy/2,dx/2]
	2144	lat2 = lat0[dy/2,dx/2]
	2145
	2146	print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \
	2147	xlon, ',', xlat
	2148
	2149	if map_proj == 'cyl':
	2150	m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \
	2151	urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	2152	elif map_proj == 'lcc':
	2153	m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
	2154	llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	2155
	2156	lons, lats = np.meshgrid(olon[:], olat[:])
	2157	lons = np.where(lons < 0., lons + 360., lons)
	2158
	2159	x,y = m(lons,lats)
	2160	plt.plot(x,y)
	2161
	2162	m.drawcoastlines()
	2163
	2164	meridians = pretty_int(nlon,xlon,5)
	2165	m.drawmeridians(meridians,labels=[True,False,False,True])
	2166
	2167	parallels = pretty_int(nlat,xlat,5)
	2168	m.drawparallels(parallels,labels=[False,True,True,False])
	2169	else:
	2170	# plt.xlim(0,dx-1)
	2171	# plt.ylim(0,dy-1)
	2172
	2173	plt.plot(xval, yval, '.')
	2174
	2175	figname = ifile.replace('.','_') + '_' + vtit
	2176	graphtit = vtit.replace('_','\_')
	2177
	2178	plt.title(graphtit)
	2179
	2180	output_kind(kfig, figname, True)
	2181
	2182	return
	2183
	2184	def plot_2Dfield(varv,dimn,colorbar,vn,vx,unit,olon,olat,ifile,vtit,zvalue,time,tk, \
	2185	tkt,tobj,tvals,tind,kfig,mapv,reva):
	2186	""" Adding labels and other staff to the graph
	2187	varv= 2D values to plot
	2188	dimn= dimension names to plot
	2189	colorbar= name of the color bar to use
	2190	vn,vm= minmum and maximum values to plot
	2191	unit= units of the variable
	2192	olon,olat= longitude, latitude objects
	2193	ifile= name of the input file
	2194	vtit= title of the variable
	2195	zvalue= value on the z axis
	2196	time= value on the time axis
	2197	tk= kind of time (WRF)
	2198	tkt= kind of time taken
	2199	tobj= tim object
	2200	tvals= values of the time variable
	2201	tind= time index
	2202	kfig= kind of figure (jpg, pdf, png)
	2203	mapv= map characteristics: [proj],[res]
	2204	see full documentation: http://matplotlib.org/basemap/
	2205	[proj]: projection
	2206	* 'cyl', cilindric
	2207	[res]: resolution:
	2208	* 'c', crude
	2209	* 'l', low
	2210	* 'i', intermediate
	2211	* 'h', high
	2212	* 'f', full
	2213	reva= reverse the axes (x-->y, y-->x)
	2214	"""
	2215	## import matplotlib as mpl
	2216	## mpl.use('Agg')
	2217	## import matplotlib.pyplot as plt
	2218
	2219	if reva:
	2220	print ' reversing the axes of the figure (x-->y, y-->x)!!'
	2221	varv = np.transpose(varv)
	2222	dimn0 = []
	2223	dimn0.append(dimn[1] + '')
	2224	dimn0.append(dimn[0] + '')
	2225	dimn = dimn0
	2226
	2227	fname = 'plot_2Dfield'
	2228	dx=varv.shape[1]
	2229	dy=varv.shape[0]
	2230
	2231	plt.rc('text', usetex=True)
	2232	# plt.rc('font', family='serif')
	2233
	2234	if not mapv is None:
	2235	if len(olon[:].shape) == 3:
	2236	lon0 = np.where(olon[0,] < 0., 360. + olon[0,], olon[0,])
	2237	lat0 = olat[0,]
	2238	elif len(olon[:].shape) == 2:
	2239	lon0 = np.where(olon[:] < 0., 360. + olon[:], olon[:])
	2240	lat0 = olat[:]
	2241	elif len(olon[:].shape) == 1:
	2242	lon00 = np.where(olon[:] < 0., 360. + olon[:], olon[:])
	2243	lat00 = olat[:]
	2244	lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
	2245	lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
	2246
	2247	for iy in range(len(lat00)):
	2248	lon0[iy,:] = lon00
	2249	for ix in range(len(lon00)):
	2250	lat0[:,ix] = lat00
	2251
	2252	map_proj=mapv.split(',')[0]
	2253	map_res=mapv.split(',')[1]
	2254
	2255	nlon = lon0[0,0]
	2256	xlon = lon0[dy-1,dx-1]
	2257	nlat = lat0[0,0]
	2258	xlat = lat0[dy-1,dx-1]
	2259
	2260	lon2 = lon0[dy/2,dx/2]
	2261	lat2 = lat0[dy/2,dx/2]
	2262
	2263	print ' lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \
	2264	xlon, ',', xlat
	2265
	2266	if map_proj == 'cyl':
	2267	m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \
	2268	urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	2269	elif map_proj == 'lcc':
	2270	m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
	2271	llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	2272
	2273	if len(olon[:].shape) == 1:
	2274	lons, lats = np.meshgrid(olon[:], olat[:])
	2275	else:
	2276	lons = olon[0,:]
	2277	lats = olat[:,0]
	2278
	2279	lons = np.where(lons < 0., lons + 360., lons)
	2280
	2281	x,y = m(lons,lats)
	2282	plt.pcolormesh(x,y,varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx)
	2283	cbar = plt.colorbar()
	2284
	2285	m.drawcoastlines()
	2286	# if (nlon > 180. or xlon > 180.):
	2287	# nlon0 = nlon
	2288	# xlon0 = xlon
	2289	# if (nlon > 180.): nlon0 = nlon - 360.
	2290	# if (xlon > 180.): xlon0 = xlon - 360.
	2291	# meridians = pretty_int(nlon0,xlon0,5)
	2292	# meridians = np.where(meridians < 0., meridians + 360., meridians)
	2293	# else:
	2294	# meridians = pretty_int(nlon,xlon,5)
	2295
	2296	meridians = pretty_int(nlon,xlon,5)
	2297
	2298	m.drawmeridians(meridians,labels=[True,False,False,True])
	2299	parallels = pretty_int(nlat,xlat,5)
	2300	m.drawparallels(parallels,labels=[False,True,True,False])
	2301
	2302	else:
	2303	plt.xlim(0,dx-1)
	2304	plt.ylim(0,dy-1)
	2305
	2306	plt.pcolormesh(varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx)
	2307	cbar = plt.colorbar()
	2308
	2309	plt.xlabel(dimn[1].replace('_','\_'))
	2310	plt.ylabel(dimn[0].replace('_','\_'))
	2311
	2312	# set the limits of the plot to the limits of the data
	2313	# plt.axis([x.min(), x.max(), y.min(), y.max()])
	2314
	2315	# plt.plot(varv)
	2316	cbar.set_label(unit)
	2317
	2318	figname = ifile.replace('.','_') + '_' + vtit
	2319	graphtit = vtit.replace('_','\_')
	2320
	2321	if zvalue != 'null':
	2322	graphtit = graphtit + ' at z= ' + zvalue
	2323	figname = figname + '_z' + zvalue
	2324	if tkt == 'tstep':
	2325	graphtit = graphtit + ' at time-step= ' + time.split(',')[1]
	2326	figname = figname + '_t' + time.split(',')[1].zfill(4)
	2327	elif tkt == 'CFdate':
	2328	graphtit = graphtit + ' at ' + tobj.strfmt("%Y%m%d%H%M%S")
	2329	figname = figname + '_t' + tobj.strfmt("%Y%m%d%H%M%S")
	2330
	2331	if tk == 'WRF':
	2332	# datev = str(timevals[timeind][0:9])
	2333	datev = tvals[tind][0] + tvals[tind][1] + tvals[tind][2] + \
	2334	timevals[timeind][3] + timevals[timeind][4] + timevals[timeind][5] + \
	2335	timevals[timeind][6] + timevals[timeind][7] + timevals[timeind][8] + \
	2336	timevals[timeind][9]
	2337	# timev = str(timevals[timeind][11:18])
	2338	timev = timevals[timeind][11] + timevals[timeind][12] + \
	2339	timevals[timeind][13] + timevals[timeind][14] + timevals[timeind][15] + \
	2340	timevals[timeind][16] + timevals[timeind][17] + timevals[timeind][18]
	2341	graphtit = vtit.replace('_','\_') + ' (' + datev + ' ' + timev + ')'
	2342
	2343	plt.title(graphtit)
	2344
	2345	output_kind(kfig, figname, True)
	2346
	2347	return
	2348
	2349	def plot_2Dfield_easy(varv,dimxv,dimyv,dimn,colorbar,vn,vx,unit,ifile,vtit,kfig,reva):
	2350	""" Adding labels and other staff to the graph
	2351	varv= 2D values to plot
	2352	dim[x/y]v = values at the axes of x and y
	2353	dimn= dimension names to plot
	2354	colorbar= name of the color bar to use
	2355	vn,vm= minmum and maximum values to plot
	2356	unit= units of the variable
	2357	ifile= name of the input file
	2358	vtit= title of the variable
	2359	kfig= kind of figure (jpg, pdf, png)
	2360	reva= reverse the axes (x-->y, y-->x)
	2361	"""
	2362	## import matplotlib as mpl
	2363	## mpl.use('Agg')
	2364	## import matplotlib.pyplot as plt
	2365	fname = 'plot_2Dfield'
	2366
	2367	if reva:
	2368	print ' reversing the axes of the figure (x-->y, y-->x)!!'
	2369	varv = np.transpose(varv)
	2370	dimn0 = []
	2371	dimn0.append(dimn[1] + '')
	2372	dimn0.append(dimn[0] + '')
	2373	dimn = dimn0
	2374	if len(dimyv.shape) == 2:
	2375	x = np.transpose(dimyv)
	2376	else:
	2377	if len(dimxv.shape) == 2:
	2378	ddx = len(dimyv)
	2379	ddy = dimxv.shape[1]
	2380	else:
	2381	ddx = len(dimyv)
	2382	ddy = len(dimxv)
	2383
	2384	x = np.zeros((ddy,ddx), dtype=np.float)
	2385	for j in range(ddy):
	2386	x[j,:] = dimyv
	2387
	2388	if len(dimxv.shape) == 2:
	2389	y = np.transpose(dimxv)
	2390	else:
	2391	if len(dimyv.shape) == 2:
	2392	ddx = dimyv.shape[0]
	2393	ddy = len(dimxv)
	2394	else:
	2395	ddx = len(dimyv)
	2396	ddy = len(dimxv)
	2397
	2398	y = np.zeros((ddy,ddx), dtype=np.float)
	2399	for i in range(ddx):
	2400	y[:,i] = dimxv
	2401	else:
	2402	if len(dimxv.shape) == 2:
	2403	x = dimxv
	2404	else:
	2405	x = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
	2406	for j in range(len(dimyv)):
	2407	x[j,:] = dimxv
	2408
	2409	if len(dimyv.shape) == 2:
	2410	y = dimyv
	2411	else:
	2412	y = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
	2413	for i in range(len(dimxv)):
	2414	x[:,i] = dimyv
	2415
	2416	dx=varv.shape[1]
	2417	dy=varv.shape[0]
	2418
	2419	plt.rc('text', usetex=True)
	2420	plt.xlim(0,dx-1)
	2421	plt.ylim(0,dy-1)
	2422
	2423	plt.pcolormesh(x, y, varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx)
	2424	# plt.pcolormesh(varv, cmap=plt.get_cmap(colorbar), vmin=vn, vmax=vx)
	2425	cbar = plt.colorbar()
	2426
	2427	plt.xlabel(dimn[1].replace('_','\_'))
	2428	plt.ylabel(dimn[0].replace('_','\_'))
	2429
	2430	# set the limits of the plot to the limits of the data
	2431	plt.axis([x.min(), x.max(), y.min(), y.max()])
	2432	# if varv.shape[1] / varv.shape[0] > 10:
	2433	# plt.axes().set_aspect(0.001)
	2434	# else:
	2435	# plt.axes().set_aspect(np.float(varv.shape[0])/np.float(varv.shape[1]))
	2436
	2437	cbar.set_label(unit)
	2438
	2439	figname = ifile.replace('.','_') + '_' + vtit
	2440	graphtit = vtit.replace('_','\_')
	2441
	2442	plt.title(graphtit)
	2443
	2444	output_kind(kfig, figname, True)
	2445
	2446	return
	2447
	2448	def plot_Trajectories(lonval, latval, linesn, olon, olat, lonlatLims, gtit, kfig, \
	2449	mapv, obsname):
	2450	""" plotting points
	2451	[lon/latval]= lon,lat values to plot (as list of vectors)
	2452	linesn: name of the lines
	2453	o[lon/lat]= object with the longitudes and the latitudes of the map to plot
	2454	lonlatLims: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax]
	2455	gtit= title of the graph
	2456	kfig= kind of figure (jpg, pdf, png)
	2457	mapv= map characteristics: [proj],[res]
	2458	see full documentation: http://matplotlib.org/basemap/
	2459	[proj]: projection
	2460	* 'cyl', cilindric
	2461	* 'lcc', lambert conformal
	2462	[res]: resolution:
	2463	* 'c', crude
	2464	* 'l', low
	2465	* 'i', intermediate
	2466	* 'h', high
	2467	* 'f', full
	2468	obsname= name of the observations in graph (can be None for without).
	2469	Observational trajectory would be the last one
	2470	"""
	2471	fname = 'plot_Trajectories'
	2472
	2473	if lonval == 'h':
	2474	print fname + '_____________________________________________________________'
	2475	print plot_Trajectories.__doc__
	2476	quit()
	2477
	2478	# Canging line kinds every 7 lines (end of standard colors)
	2479	linekinds=['.-','x-','o-']
	2480
	2481	Ntraj = len(lonval)
	2482
	2483	if obsname is not None:
	2484	Ntraj = Ntraj - 1
	2485
	2486	N7lines = 0
	2487
	2488	plt.rc('text', usetex=True)
	2489
	2490	if not mapv is None:
	2491	if len(olon[:].shape) == 3:
	2492	# lon0 = np.where(olon[0,] < 0., 360. + olon[0,], olon[0,])
	2493	lon0 = olon[0,]
	2494	lat0 = olat[0,]
	2495	elif len(olon[:].shape) == 2:
	2496	# lon0 = np.where(olon[:] < 0., 360. + olon[:], olon[:])
	2497	lon0 = olon[:]
	2498	lat0 = olat[:]
	2499	elif len(olon[:].shape) == 1:
	2500	# lon00 = np.where(olon[:] < 0., 360. + olon[:], olon[:])
	2501	lon00 = olon[:]
	2502	lat00 = olat[:]
	2503	lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
	2504	lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
	2505
	2506	for iy in range(len(lat00)):
	2507	lon0[iy,:] = lon00
	2508	for ix in range(len(lon00)):
	2509	lat0[:,ix] = lat00
	2510
	2511	map_proj=mapv.split(',')[0]
	2512	map_res=mapv.split(',')[1]
	2513
	2514	dx = lon0.shape[1]
	2515	dy = lon0.shape[0]
	2516
	2517	nlon = lon0[0,0]
	2518	xlon = lon0[dy-1,dx-1]
	2519	nlat = lat0[0,0]
	2520	xlat = lat0[dy-1,dx-1]
	2521
	2522	lon2 = lon0[dy/2,dx/2]
	2523	lat2 = lat0[dy/2,dx/2]
	2524
	2525	if lonlatLims is not None:
	2526	plt.xlim(lonlatLims[0], lonlatLims[2])
	2527	plt.ylim(lonlatLims[1], lonlatLims[3])
	2528	if map_proj == 'cyl':
	2529	nlon = lonlatLims[0]
	2530	nlat = lonlatLims[1]
	2531	xlon = lonlatLims[2]
	2532	xlat = lonlatLims[3]
	2533
	2534	print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \
	2535	xlon, ',', xlat
	2536
	2537	if map_proj == 'cyl':
	2538	m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \
	2539	urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	2540	elif map_proj == 'lcc':
	2541	m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
	2542	llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	2543
	2544	if len(olon.shape) == 3:
	2545	# lons, lats = np.meshgrid(olon[0,:,:], olat[0,:,:])
	2546	lons = olon[0,:,:]
	2547	lats = olat[0,:,:]
	2548
	2549	elif len(olon.shape) == 2:
	2550	# lons, lats = np.meshgrid(olon[:,:], olat[:,:])
	2551	lons = olon[:,:]
	2552	lats = olat[:,:]
	2553	else:
	2554	print errormsg
	2555	print ' ' + fname + ': shapes of lon/lat objects', olon.shape, \
	2556	'not ready!!!'
	2557
	2558	for il in range(Ntraj):
	2559	plt.plot(lonval[il], latval[il], linekinds[N7lines], label= linesn[il])
	2560	if il == 6: N7lines = N7lines + 1
	2561
	2562	m.drawcoastlines()
	2563
	2564	meridians = pretty_int(nlon,xlon,5)
	2565	m.drawmeridians(meridians,labels=[True,False,False,True])
	2566
	2567	parallels = pretty_int(nlat,xlat,5)
	2568	m.drawparallels(parallels,labels=[False,True,True,False])
	2569
	2570	plt.xlabel('W-E')
	2571	plt.ylabel('S-N')
	2572
	2573	else:
	2574	if len(olon.shape) == 3:
	2575	dx = olon.shape[2]
	2576	dy = olon.shape[1]
	2577	elif len(olon.shape) == 2:
	2578	dx = olon.shape[1]
	2579	dy = olon.shape[0]
	2580	else:
	2581	print errormsg
	2582	print ' ' + fname + ': shapes of lon/lat objects', olon.shape, \
	2583	'not ready!!!'
	2584
	2585	if lonlatLims is not None:
	2586	plt.xlim(lonlatLims[0], lonlatLims[2])
	2587	plt.ylim(lonlatLims[1], lonlatLims[3])
	2588	else:
	2589	plt.xlim(np.min(olon[:]),np.max(olon[:]))
	2590	plt.ylim(np.min(olat[:]),np.max(olat[:]))
	2591
	2592	for il in range(Ntraj):
	2593	plt.plot(lonval[il], latval[il], linekinds[N7lines], label= linesn[il])
	2594	if il == 6: N7lines = N7lines + 1
	2595
	2596	plt.xlabel('x-axis')
	2597	plt.ylabel('y-axis')
	2598
	2599	figname = 'trajectories'
	2600	graphtit = gtit
	2601
	2602	if obsname is not None:
	2603	plt.plot(lonval[Ntraj], latval[Ntraj], linestyle='-', color='k', \
	2604	linewidth=3, label= obsname)
	2605
	2606	plt.title(graphtit)
	2607	plt.legend()
	2608
	2609	output_kind(kfig, figname, True)
	2610
	2611	return
	2612
	2613	def plot_topo_geogrid(varv, olon, olat, mint, maxt, lonlatLims, gtit, kfig, mapv, \
	2614	closeif):
	2615	""" plotting geo_em.d[nn].nc topography from WPS files
	2616	plot_topo_geogrid(domf, mint, maxt, gtit, kfig, mapv)
	2617	varv= topography values
	2618	o[lon/lat]= longitude and latitude objects
	2619	[min/max]t: minimum and maximum values of topography to draw
	2620	lonlatLims: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax]
	2621	gtit= title of the graph
	2622	kfig= kind of figure (jpg, pdf, png)
	2623	mapv= map characteristics: [proj],[res]
	2624	see full documentation: http://matplotlib.org/basemap/
	2625	[proj]: projection
	2626	* 'cyl', cilindric
	2627	* 'lcc', lamvbert conformal
	2628	[res]: resolution:
	2629	* 'c', crude
	2630	* 'l', low
	2631	* 'i', intermediate
	2632	* 'h', high
	2633	* 'f', full
	2634	closeif= Boolean value if the figure has to be closed
	2635	"""
	2636	fname = 'plot_topo_geogrid'
	2637
	2638	if varv == 'h':
	2639	print fname + '_____________________________________________________________'
	2640	print plot_topo_geogrid.__doc__
	2641	quit()
	2642
	2643	dx=varv.shape[1]
	2644	dy=varv.shape[0]
	2645
	2646	plt.rc('text', usetex=True)
	2647	# plt.rc('font', family='serif')
	2648
	2649	if not mapv is None:
	2650	if len(olon[:].shape) == 3:
	2651	lon0 = olon[0,]
	2652	lat0 = olat[0,]
	2653	elif len(olon[:].shape) == 2:
	2654	lon0 = olon[:]
	2655	lat0 = olat[:]
	2656	elif len(olon[:].shape) == 1:
	2657	lon00 = olon[:]
	2658	lat00 = olat[:]
	2659	lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
	2660	lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
	2661
	2662	for iy in range(len(lat00)):
	2663	lon0[iy,:] = lon00
	2664	for ix in range(len(lon00)):
	2665	lat0[:,ix] = lat00
	2666
	2667	map_proj=mapv.split(',')[0]
	2668	map_res=mapv.split(',')[1]
	2669	dx = lon0.shape[1]
	2670	dy = lon0.shape[0]
	2671
	2672	if lonlatLims is not None:
	2673	print ' ' + fname + ': cutting the domain to plot !!!!'
	2674	print ' limits: W-E', lonlatLims[0], lonlatLims[2]
	2675	print ' limits: N-S', lonlatLims[1], lonlatLims[3]
	2676	nlon = lonlatLims[0]
	2677	xlon = lonlatLims[2]
	2678	nlat = lonlatLims[1]
	2679	xlat = lonlatLims[3]
	2680
	2681	if map_proj == 'lcc':
	2682	lon2 = (lonlatLims[0] + lonlatLims[2])/2.
	2683	lat2 = (lonlatLims[1] + lonlatLims[3])/2.
	2684	else:
	2685	nlon = lon0[0,0]
	2686	xlon = lon0[dy-1,dx-1]
	2687	nlat = lat0[0,0]
	2688	xlat = lat0[dy-1,dx-1]
	2689	lon2 = lon0[dy/2,dx/2]
	2690	lat2 = lat0[dy/2,dx/2]
	2691
	2692	plt.xlim(nlon, xlon)
	2693	plt.ylim(nlat, xlat)
	2694	print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \
	2695	xlon, ',', xlat
	2696
	2697	if map_proj == 'cyl':
	2698	m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \
	2699	urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	2700	elif map_proj == 'lcc':
	2701	m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
	2702	llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	2703	else:
	2704	print errormsg
	2705	print ' ' + fname + ": map projection '" + map_proj + "' not ready !!"
	2706	quit(-1)
	2707
	2708	if len(olon[:].shape) == 1:
	2709	lons, lats = np.meshgrid(olon[:], olat[:])
	2710	else:
	2711	if len(olon[:].shape) == 3:
	2712	lons = olon[0,:,:]
	2713	lats = olat[0,:,:]
	2714	else:
	2715	lons = olon[:]
	2716	lats = olat[:]
	2717
	2718	x,y = m(lons,lats)
	2719
	2720	plt.pcolormesh(x,y,varv, cmap=plt.get_cmap('terrain'), vmin=mint, vmax=maxt)
	2721	cbar = plt.colorbar()
	2722
	2723	m.drawcoastlines()
	2724
	2725	meridians = pretty_int(nlon,xlon,5)
	2726	m.drawmeridians(meridians,labels=[True,False,False,True])
	2727
	2728	parallels = pretty_int(nlat,xlat,5)
	2729	m.drawparallels(parallels,labels=[False,True,True,False])
	2730
	2731	plt.xlabel('W-E')
	2732	plt.ylabel('S-N')
	2733	else:
	2734	print emsg
	2735	print ' ' + fname + ': A projection parameter is needed None given !!'
	2736	quit(-1)
	2737
	2738	figname = 'domain'
	2739	graphtit = gtit.replace('_','\_')
	2740	cbar.set_label('height ($m$)')
	2741
	2742	plt.title(graphtit)
	2743
	2744	output_kind(kfig, figname, closeif)
	2745
	2746	return
	2747
	2748	def plot_topo_geogrid_boxes(varv, boxesX, boxesY, boxlabels, olon, olat, mint, maxt, \
	2749	lonlatLims, gtit, kfig, mapv, closeif):
	2750	""" plotting geo_em.d[nn].nc topography from WPS files
	2751	plot_topo_geogrid(domf, mint, maxt, gtit, kfig, mapv)
	2752	varv= topography values
	2753	boxesX/Y= 4-line sets to draw the boxes
	2754	boxlabels= labels for the legend of the boxes
	2755	o[lon/lat]= longitude and latitude objects
	2756	[min/max]t: minimum and maximum values of topography to draw
	2757	lonlatLims: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax]
	2758	gtit= title of the graph
	2759	kfig= kind of figure (jpg, pdf, png)
	2760	mapv= map characteristics: [proj],[res]
	2761	see full documentation: http://matplotlib.org/basemap/
	2762	[proj]: projection
	2763	* 'cyl', cilindric
	2764	* 'lcc', lamvbert conformal
	2765	[res]: resolution:
	2766	* 'c', crude
	2767	* 'l', low
	2768	* 'i', intermediate
	2769	* 'h', high
	2770	* 'f', full
	2771	closeif= Boolean value if the figure has to be closed
	2772	"""
	2773	fname = 'plot_topo_geogrid'
	2774
	2775	if varv == 'h':
	2776	print fname + '_____________________________________________________________'
	2777	print plot_topo_geogrid.__doc__
	2778	quit()
	2779
	2780	cols = color_lines(len(boxlabels))
	2781
	2782	dx=varv.shape[1]
	2783	dy=varv.shape[0]
	2784
	2785	plt.rc('text', usetex=True)
	2786	# plt.rc('font', family='serif')
	2787
	2788	if not mapv is None:
	2789	if len(olon[:].shape) == 3:
	2790	lon0 = olon[0,]
	2791	lat0 = olat[0,]
	2792	elif len(olon[:].shape) == 2:
	2793	lon0 = olon[:]
	2794	lat0 = olat[:]
	2795	elif len(olon[:].shape) == 1:
	2796	lon00 = olon[:]
	2797	lat00 = olat[:]
	2798	lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
	2799	lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
	2800
	2801	for iy in range(len(lat00)):
	2802	lon0[iy,:] = lon00
	2803	for ix in range(len(lon00)):
	2804	lat0[:,ix] = lat00
	2805
	2806	map_proj=mapv.split(',')[0]
	2807	map_res=mapv.split(',')[1]
	2808	dx = lon0.shape[1]
	2809	dy = lon0.shape[0]
	2810
	2811	if lonlatLims is not None:
	2812	print ' ' + fname + ': cutting the domain to plot !!!!'
	2813	print ' limits: W-E', lonlatLims[0], lonlatLims[2]
	2814	print ' limits: N-S', lonlatLims[1], lonlatLims[3]
	2815	nlon = lonlatLims[0]
	2816	xlon = lonlatLims[2]
	2817	nlat = lonlatLims[1]
	2818	xlat = lonlatLims[3]
	2819
	2820	if map_proj == 'lcc':
	2821	lon2 = (lonlatLims[0] + lonlatLims[2])/2.
	2822	lat2 = (lonlatLims[1] + lonlatLims[3])/2.
	2823	else:
	2824	nlon = np.min(lon0)
	2825	xlon = np.max(lon0)
	2826	nlat = np.min(lat0)
	2827	xlat = np.max(lat0)
	2828	lon2 = lon0[dy/2,dx/2]
	2829	lat2 = lat0[dy/2,dx/2]
	2830
	2831	plt.xlim(nlon, xlon)
	2832	plt.ylim(nlat, xlat)
	2833	print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \
	2834	xlon, ',', xlat
	2835
	2836	if map_proj == 'cyl':
	2837	m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \
	2838	urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	2839	elif map_proj == 'lcc':
	2840	m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
	2841	llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	2842
	2843	if len(olon[:].shape) == 1:
	2844	lons, lats = np.meshgrid(olon[:], olat[:])
	2845	else:
	2846	if len(olon[:].shape) == 3:
	2847	lons = olon[0,:,:]
	2848	lats = olat[0,:,:]
	2849	else:
	2850	lons = olon[:]
	2851	lats = olat[:]
	2852
	2853	x,y = m(lons,lats)
	2854
	2855	plt.pcolormesh(x,y,varv, cmap=plt.get_cmap('terrain'), vmin=mint, vmax=maxt)
	2856	cbar = plt.colorbar()
	2857
	2858	Nboxes = len(boxesX)/4
	2859	for ibox in range(Nboxes):
	2860	plt.plot(boxesX[ibox4], boxesY[ibox4], linestyle='-', linewidth=3, \
	2861	label=boxlabels[ibox], color=cols[ibox])
	2862	plt.plot(boxesX[ibox4+1], boxesY[ibox4+1], linestyle='-', linewidth=3, \
	2863	color=cols[ibox])
	2864	plt.plot(boxesX[ibox4+2], boxesY[ibox4+2], linestyle='-', linewidth=3, \
	2865	color=cols[ibox])
	2866	plt.plot(boxesX[ibox4+3], boxesY[ibox4+3], linestyle='-', linewidth=3, \
	2867	color=cols[ibox])
	2868
	2869	m.drawcoastlines()
	2870
	2871	meridians = pretty_int(nlon,xlon,5)
	2872	m.drawmeridians(meridians,labels=[True,False,False,True])
	2873
	2874	parallels = pretty_int(nlat,xlat,5)
	2875	m.drawparallels(parallels,labels=[False,True,True,False])
	2876
	2877	plt.xlabel('W-E')
	2878	plt.ylabel('S-N')
	2879	else:
	2880	print emsg
	2881	print ' ' + fname + ': A projection parameter is needed None given !!'
	2882	quit(-1)
	2883
	2884	figname = 'domain_boxes'
	2885	graphtit = gtit.replace('_','\_')
	2886	cbar.set_label('height ($m$)')
	2887
	2888	plt.title(graphtit)
	2889	plt.legend()
	2890
	2891	output_kind(kfig, figname, closeif)
	2892
	2893	return
	2894
	2895	def plot_2D_shadow(varsv,vnames,dimxv,dimyv,dimxu,dimyu,dimn, \
	2896	colorbar,vs,uts,vtit,kfig,reva,mapv,ifclose):
	2897	""" Adding labels and other staff to the graph
	2898	varsv= 2D values to plot with shading
	2899	vnames= variable names for the figure
	2900	dim[x/y]v = values at the axes of x and y
	2901	dim[x/y]u = units at the axes of x and y
	2902	dimn= dimension names to plot
	2903	colorbar= name of the color bar to use
	2904	vs= minmum and maximum values to plot in shadow or:
	2905	'Srange': for full range
	2906	'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
	2907	'Saroundminmax@val': for minval,maxval
	2908	'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
	2909	percentile_(100-val)-median)
	2910	'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min@val,max@val-mean)
	2911	'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
	2912	'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
	2913	percentile_(100-val)-median)
	2914	uts= units of the variable to shadow
	2915	vtit= title of the variable
	2916	kfig= kind of figure (jpg, pdf, png)
	2917	reva= ('\|' for combination)
	2918	* 'transpose': reverse the axes (x-->y, y-->x)
	2919	* 'flip'@[x/y]: flip the axis x or y
	2920	mapv= map characteristics: [proj],[res]
	2921	see full documentation: http://matplotlib.org/basemap/
	2922	[proj]: projection
	2923	* 'cyl', cilindric
	2924	* 'lcc', lambert conformal
	2925	[res]: resolution:
	2926	* 'c', crude
	2927	* 'l', low
	2928	* 'i', intermediate
	2929	* 'h', high
	2930	* 'f', full
	2931	ifclose= boolean value whether figure should be close (finish) or not
	2932	"""
	2933	## import matplotlib as mpl
	2934	## mpl.use('Agg')
	2935	## import matplotlib.pyplot as plt
	2936	fname = 'plot_2D_shadow'
	2937
	2938	# print dimyv[73,21]
	2939	# dimyv[73,21] = -dimyv[73,21]
	2940	# print 'Lluis dimsv: ',np.min(dimxv), np.max(dimxv), ':', np.min(dimyv), np.max(dimyv)
	2941
	2942	if varsv == 'h':
	2943	print fname + '_____________________________________________________________'
	2944	print plot_2D_shadow.__doc__
	2945	quit()
	2946
	2947	if len(varsv.shape) != 2:
	2948	print errormsg
	2949	print ' ' + fname + ': wrong variable shape:',varv.shape,'is has to be 2D!!'
	2950	quit(-1)
	2951
	2952	reva0 = ''
	2953	if reva.find('\|') != 0:
	2954	revas = reva.split('\|')
	2955	else:
	2956	revas = [reva]
	2957	reva0 = reva
	2958
	2959	for rev in revas:
	2960	if reva[0:4] == 'flip':
	2961	reva0 = 'flip'
	2962	if len(reva.split('@')) != 2:
	2963	print errormsg
	2964	print ' ' + fname + ': flip is given', reva, 'but not axis!'
	2965	quit(-1)
	2966
	2967	if rev == 'transpose':
	2968	print ' reversing the axes of the figure (x-->y, y-->x)!!'
	2969	varsv = np.transpose(varsv)
	2970	dxv = dimyv
	2971	dyv = dimxv
	2972	dimxv = dxv
	2973	dimyv = dyv
	2974
	2975	if len(dimxv[:].shape) == 3:
	2976	xdims = '1,2'
	2977	elif len(dimxv[:].shape) == 2:
	2978	xdims = '0,1'
	2979	elif len(dimxv[:].shape) == 1:
	2980	xdims = '0'
	2981
	2982	if len(dimyv[:].shape) == 3:
	2983	ydims = '1,2'
	2984	elif len(dimyv[:].shape) == 2:
	2985	ydims = '0,1'
	2986	elif len(dimyv[:].shape) == 1:
	2987	ydims = '0'
	2988
	2989	lon0, lat0 = dxdy_lonlat(dimxv,dimyv, xdims, ydims)
	2990
	2991	if not mapv is None:
	2992	map_proj=mapv.split(',')[0]
	2993	map_res=mapv.split(',')[1]
	2994
	2995	dx = lon0.shape[1]
	2996	dy = lat0.shape[0]
	2997
	2998	nlon = lon0[0,0]
	2999	xlon = lon0[dy-1,dx-1]
	3000	nlat = lat0[0,0]
	3001	xlat = lat0[dy-1,dx-1]
	3002
	3003	# Thats too much! :)
	3004	# if lonlatLims is not None:
	3005	# print ' ' + fname + ': cutting the domain to plot !!!!'
	3006	# plt.xlim(lonlatLims[0], lonlatLims[2])
	3007	# plt.ylim(lonlatLims[1], lonlatLims[3])
	3008	# print ' limits: W-E', lonlatLims[0], lonlatLims[2]
	3009	# print ' limits: N-S', lonlatLims[1], lonlatLims[3]
	3010
	3011	# if map_proj == 'cyl':
	3012	# nlon = lonlatLims[0]
	3013	# nlat = lonlatLims[1]
	3014	# xlon = lonlatLims[2]
	3015	# xlat = lonlatLims[3]
	3016	# elif map_proj == 'lcc':
	3017	# lon2 = (lonlatLims[0] + lonlatLims[2])/2.
	3018	# lat2 = (lonlatLims[1] + lonlatLims[3])/2.
	3019	# nlon = lonlatLims[0]
	3020	# xlon = lonlatLims[2]
	3021	# nlat = lonlatLims[1]
	3022	# xlat = lonlatLims[3]
	3023
	3024	lon2 = lon0[dy/2,dx/2]
	3025	lat2 = lat0[dy/2,dx/2]
	3026
	3027	print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \
	3028	xlon, ',', xlat
	3029
	3030	if map_proj == 'cyl':
	3031	m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \
	3032	urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	3033	elif map_proj == 'lcc':
	3034	m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
	3035	llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	3036	else:
	3037	print errormsg
	3038	print ' ' + fname + ": map projection '" + map_proj + "' not defined!!!"
	3039	print ' available: cyl, lcc'
	3040	quit(-1)
	3041
	3042	x,y = m(lon0,lat0)
	3043
	3044	else:
	3045	x = lon0
	3046	y = lat0
	3047
	3048	vsend = np.zeros((2), dtype=np.float)
	3049	# Changing limits of the colors
	3050	if type(vs[0]) != type(np.float(1.)):
	3051	if vs[0] == 'Srange':
	3052	vsend[0] = np.min(varsv)
	3053	elif vs[0][0:11] == 'Saroundmean':
	3054	meanv = np.mean(varsv)
	3055	permean = np.float(vs[0].split('@')[1])
	3056	minv = np.min(varsv)*permean
	3057	maxv = np.max(varsv)*permean
	3058	minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
	3059	vsend[0] = meanv-minextrm
	3060	vsend[1] = meanv+minextrm
	3061	elif vs[0][0:13] == 'Saroundminmax':
	3062	permean = np.float(vs[0].split('@')[1])
	3063	minv = np.min(varsv)*permean
	3064	maxv = np.max(varsv)*permean
	3065	vsend[0] = minv
	3066	vsend[1] = maxv
	3067	elif vs[0][0:17] == 'Saroundpercentile':
	3068	medianv = np.median(varsv)
	3069	valper = np.float(vs[0].split('@')[1])
	3070	minv = np.percentile(varsv, valper)
	3071	maxv = np.percentile(varsv, 100.-valper)
	3072	minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
	3073	vsend[0] = medianv-minextrm
	3074	vsend[1] = medianv+minextrm
	3075	elif vs[0][0:5] == 'Smean':
	3076	meanv = np.mean(varsv)
	3077	permean = np.float(vs[0].split('@')[1])
	3078	minv = np.min(varsv)*permean
	3079	maxv = np.max(varsv)*permean
	3080	minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
	3081	vsend[0] = -minextrm
	3082	vsend[1] = minextrm
	3083	elif vs[0][0:7] == 'Smedian':
	3084	medianv = np.median(varsv)
	3085	permedian = np.float(vs[0].split('@')[1])
	3086	minv = np.min(varsv)*permedian
	3087	maxv = np.max(varsv)*permedian
	3088	minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
	3089	vsend[0] = -minextrm
	3090	vsend[1] = minextrm
	3091	elif vs[0][0:11] == 'Spercentile':
	3092	medianv = np.median(varsv)
	3093	valper = np.float(vs[0].split('@')[1])
	3094	minv = np.percentile(varsv, valper)
	3095	maxv = np.percentile(varsv, 100.-valper)
	3096	minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
	3097	vsend[0] = -minextrm
	3098	vsend[1] = minextrm
	3099	else:
	3100	print errormsg
	3101	print ' ' + fname + ": range '" + vs[0] + "' not ready!!!"
	3102	quit(-1)
	3103	print ' ' + fname + ': modified shadow min,max:',vsend
	3104	else:
	3105	vsend[0] = vs[0]
	3106
	3107	if type(vs[0]) != type(np.float(1.)):
	3108	if vs[1] == 'range':
	3109	vsend[1] = np.max(varsv)
	3110	else:
	3111	vsend[1] = vs[1]
	3112
	3113	plt.rc('text', usetex=True)
	3114
	3115	plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(colorbar), vmin=vsend[0], vmax=vsend[1])
	3116	cbar = plt.colorbar()
	3117
	3118	if not mapv is None:
	3119	m.drawcoastlines()
	3120
	3121	meridians = pretty_int(nlon,xlon,5)
	3122	m.drawmeridians(meridians,labels=[True,False,False,True])
	3123	parallels = pretty_int(nlat,xlat,5)
	3124	m.drawparallels(parallels,labels=[False,True,True,False])
	3125
	3126	plt.xlabel('W-E')
	3127	plt.ylabel('S-N')
	3128	else:
	3129	plt.xlabel(variables_values(dimn[1])[0].replace('_','\_') + ' (' + \
	3130	units_lunits(dimxu) + ')')
	3131	plt.ylabel(variables_values(dimn[0])[0].replace('_','\_') + ' (' + \
	3132	units_lunits(dimyu) + ')')
	3133
	3134	txpos = pretty_int(x.min(),x.max(),5)
	3135	typos = pretty_int(y.min(),y.max(),5)
	3136	txlabels = list(txpos)
	3137	for i in range(len(txlabels)): txlabels[i] = str(txlabels[i])
	3138	tylabels = list(typos)
	3139	for i in range(len(tylabels)): tylabels[i] = str(tylabels[i])
	3140
	3141	# set the limits of the plot to the limits of the data
	3142
	3143	if searchInlist(revas,'transpose'):
	3144	x0 = y
	3145	y0 = x
	3146	x = x0
	3147	y = y0
	3148	# print 'Lluis reva0:',reva0,'x min,max:',x.min(),x.max(),' y min,max:',y.min(),y.max()
	3149
	3150	if reva0 == 'flip':
	3151	if reva.split('@')[1] == 'x':
	3152	plt.axis([x.max(), x.min(), y.min(), y.max()])
	3153	elif reva.split('@')[1] == 'y':
	3154	plt.axis([x.min(), x.max(), y.max(), y.min()])
	3155	else:
	3156	plt.axis([x.max(), x.min(), y.max(), y.min()])
	3157	else:
	3158	plt.axis([x.min(), x.max(), y.min(), y.max()])
	3159
	3160	if mapv is None:
	3161	plt.xticks(txpos, txlabels)
	3162	plt.yticks(typos, tylabels)
	3163
	3164	# units labels
	3165	cbar.set_label(vnames.replace('_','\_') + ' (' + units_lunits(uts) + ')')
	3166
	3167	figname = '2Dfields_shadow'
	3168	graphtit = vtit.replace('_','\_').replace('&','\&')
	3169
	3170	plt.title(graphtit)
	3171
	3172	output_kind(kfig, figname, ifclose)
	3173
	3174	return
	3175
	3176	#Nvals=50
	3177	#vals1 = np.zeros((Nvals,Nvals), dtype= np.float)
	3178	#for j in range(Nvals):
	3179	# for i in range(Nvals):
	3180	# vals1[j,i]=np.sqrt((j-Nvals/2)2. + (i-Nvals/2)2.)
	3181
	3182	#plot_2D_shadow(vals1, 'var1', np.arange(50)1., np.arange(50)1., 'ms-1', \
	3183	# 'm', ['lat','lon'], 'rainbow', [0, Nvals], 'ms-1', 'test var1', 'pdf', 'None', \
	3184	# None, True)
	3185	#quit()
	3186
	3187	def plot_2D_shadow_time(varsv,vnames,dimxv,dimyv,dimxu,dimyu,dimn,colorbar,vs,uts, \
	3188	vtit,kfig,reva,taxis,tpos,tlabs,ifclose):
	3189	""" Plotting a 2D field with one of the axes being time
	3190	varsv= 2D values to plot with shading
	3191	vnames= shading variable name for the figure
	3192	dim[x/y]v= values at the axes of x and y
	3193	dim[x/y]u= units at the axes of x and y
	3194	dimn= dimension names to plot
	3195	colorbar= name of the color bar to use
	3196	vs= minmum and maximum values to plot in shadow or:
	3197	'Srange': for full range
	3198	'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
	3199	'Saroundminmax@val': for minval,maxval
	3200	'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
	3201	percentile_(100-val)-median)
	3202	'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min@val,max@val-mean)
	3203	'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
	3204	'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
	3205	percentile_(100-val)-median)
	3206	uts= units of the variable to shadow
	3207	vtit= title of the variable
	3208	kfig= kind of figure (jpg, pdf, png)
	3209	reva=
	3210	* 'transpose': reverse the axes (x-->y, y-->x)
	3211	* 'flip'@[x/y]: flip the axis x or y
	3212	taxis= Which is the time-axis
	3213	tpos= positions of the time ticks
	3214	tlabs= labels of the time ticks
	3215	ifclose= boolean value whether figure should be close (finish) or not
	3216	"""
	3217	fname = 'plot_2D_shadow_time'
	3218
	3219	if varsv == 'h':
	3220	print fname + '_____________________________________________________________'
	3221	print plot_2D_shadow_time.__doc__
	3222	quit()
	3223
	3224	# Definning ticks labels
	3225	if taxis == 'x':
	3226	txpos = tpos
	3227	txlabels = tlabs
	3228	plxlabel = dimxu
	3229	typos = pretty_int(np.min(dimyv),np.max(dimyv),10)
	3230	tylabels = list(typos)
	3231	for i in range(len(tylabels)): tylabels[i] = str(tylabels[i])
	3232	plylabel = variables_values(dimn[0])[0].replace('_','\_') + ' (' + \
	3233	units_lunits(dimyu) + ')'
	3234	else:
	3235	txpos = pretty_int(np.min(dimxv),np.max(dimxv),10)
	3236	txlabels = list(txpos)
	3237	plxlabel = variables_values(dimn[1])[0].replace('_','\_') + ' (' + \
	3238	units_lunits(dimxu) + ')'
	3239	typos = tpos
	3240	tylabels = tlabs
	3241	plylabel = dimyu
	3242
	3243	# Transposing/flipping axis
	3244	if reva.find('\|') != 0:
	3245	revas = reva.split('\|')
	3246	reva0 = ''
	3247	else:
	3248	revas = [reva]
	3249	reva0 = reva
	3250
	3251	for rev in revas:
	3252	if rev[0:4] == 'flip':
	3253	reva0 = 'flip'
	3254	if len(reva.split('@')) != 2:
	3255	print errormsg
	3256	print ' ' + fname + ': flip is given', reva, 'but not axis!'
	3257	quit(-1)
	3258	else:
	3259	print " flipping '" + rev.split('@')[1] + "' axis !"
	3260
	3261	if rev == 'transpose':
	3262	print ' reversing the axes of the figure (x-->y, y-->x)!!'
	3263	# Flipping values of variable
	3264	varsv = np.transpose(varsv)
	3265	dxv = dimyv
	3266	dyv = dimxv
	3267	dimxv = dxv
	3268	dimyv = dyv
	3269
	3270	if len(dimxv.shape) == 3:
	3271	dxget='1,2'
	3272	elif len(dimxv.shape) == 2:
	3273	dxget='0,1'
	3274	elif len(dimxv.shape) == 1:
	3275	dxget='0'
	3276	else:
	3277	print errormsg
	3278	print ' ' + fname + ': shape of x-values:',dimxv.shape,'not ready!!'
	3279	quit(-1)
	3280
	3281	if len(dimyv.shape) == 3:
	3282	dyget='1,2'
	3283	elif len(dimyv.shape) == 2:
	3284	dyget='0,1'
	3285	elif len(dimyv.shape) == 1:
	3286	dyget='0'
	3287	else:
	3288	print errormsg
	3289	print ' ' + fname + ': shape of y-values:',dimyv.shape,'not ready!!'
	3290	quit(-1)
	3291
	3292	x,y = dxdy_lonlat(dimxv,dimyv,dxget,dyget)
	3293
	3294	plt.rc('text', usetex=True)
	3295
	3296	vsend = np.zeros((2), dtype=np.float)
	3297	# Changing limits of the colors
	3298	if type(vs[0]) != type(np.float(1.)):
	3299	if vs[0] == 'Srange':
	3300	vsend[0] = np.min(varsv)
	3301	elif vs[0][0:11] == 'Saroundmean':
	3302	meanv = np.mean(varsv)
	3303	permean = np.float(vs[0].split('@')[1])
	3304	minv = np.min(varsv)*permean
	3305	maxv = np.max(varsv)*permean
	3306	minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
	3307	vsend[0] = meanv-minextrm
	3308	vsend[1] = meanv+minextrm
	3309	elif vs[0][0:13] == 'Saroundminmax':
	3310	permean = np.float(vs[0].split('@')[1])
	3311	minv = np.min(varsv)*permean
	3312	maxv = np.max(varsv)*permean
	3313	vsend[0] = minv
	3314	vsend[1] = maxv
	3315	elif vs[0][0:17] == 'Saroundpercentile':
	3316	medianv = np.median(varsv)
	3317	valper = np.float(vs[0].split('@')[1])
	3318	minv = np.percentile(varsv, valper)
	3319	maxv = np.percentile(varsv, 100.-valper)
	3320	minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
	3321	vsend[0] = medianv-minextrm
	3322	vsend[1] = medianv+minextrm
	3323	elif vs[0][0:5] == 'Smean':
	3324	meanv = np.mean(varsv)
	3325	permean = np.float(vs[0].split('@')[1])
	3326	minv = np.min(varsv)*permean
	3327	maxv = np.max(varsv)*permean
	3328	minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
	3329	vsend[0] = -minextrm
	3330	vsend[1] = minextrm
	3331	elif vs[0][0:7] == 'Smedian':
	3332	medianv = np.median(varsv)
	3333	permedian = np.float(vs[0].split('@')[1])
	3334	minv = np.min(varsv)*permedian
	3335	maxv = np.max(varsv)*permedian
	3336	minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
	3337	vsend[0] = -minextrm
	3338	vsend[1] = minextrm
	3339	elif vs[0][0:11] == 'Spercentile':
	3340	medianv = np.median(varsv)
	3341	valper = np.float(vs[0].split('@')[1])
	3342	minv = np.percentile(varsv, valper)
	3343	maxv = np.percentile(varsv, 100.-valper)
	3344	minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
	3345	vsend[0] = -minextrm
	3346	vsend[1] = minextrm
	3347	else:
	3348	print errormsg
	3349	print ' ' + fname + ": range '" + vs[0] + "' not ready!!!"
	3350	quit(-1)
	3351	print ' ' + fname + ': modified shadow min,max:',vsend
	3352	else:
	3353	vsend[0] = vs[0]
	3354
	3355	if type(vs[0]) != type(np.float(1.)):
	3356	if vs[1] == 'range':
	3357	vsend[1] = np.max(varsv)
	3358	else:
	3359	vsend[1] = vs[1]
	3360
	3361	plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(colorbar), vmin=vsend[0], vmax=vsend[1])
	3362	cbar = plt.colorbar()
	3363
	3364	# print 'Lluis reva0:',reva0,'x min,max:',x.min(),x.max(),' y min,max:',y.min(),y.max()
	3365
	3366	# set the limits of the plot to the limits of the data
	3367	if reva0 == 'flip':
	3368	if reva.split('@')[1] == 'x':
	3369	plt.axis([x.max(), x.min(), y.min(), y.max()])
	3370	elif reva.split('@')[1] == 'y':
	3371	plt.axis([x.min(), x.max(), y.max(), y.min()])
	3372	else:
	3373	plt.axis([x.max(), x.min(), y.max(), y.min()])
	3374	else:
	3375	plt.axis([x.min(), x.max(), y.min(), y.max()])
	3376
	3377	if searchInlist(revas, 'transpose'):
	3378	plt.xticks(typos, tylabels)
	3379	plt.yticks(txpos, txlabels)
	3380	plt.xlabel(plylabel)
	3381	plt.ylabel(plxlabel)
	3382	else:
	3383	plt.xticks(txpos, txlabels)
	3384	plt.yticks(typos, tylabels)
	3385	plt.xlabel(plxlabel)
	3386	plt.ylabel(plylabel)
	3387
	3388	# units labels
	3389	cbar.set_label(vnames.replace('_','\_') + ' (' + units_lunits(uts) + ')')
	3390
	3391	figname = '2Dfields_shadow_time'
	3392	graphtit = vtit.replace('_','\_').replace('&','\&')
	3393
	3394	plt.title(graphtit)
	3395
	3396	output_kind(kfig, figname, ifclose)
	3397
	3398	return
	3399
	3400	def plot_2D_shadow_contour(varsv,varcv,vnames,dimxv,dimyv,dimxu,dimyu,dimn, \
	3401	colorbar,ckind,clabfmt,vs,vc,uts,vtit,kfig,reva,mapv):
	3402	""" Adding labels and other staff to the graph
	3403	varsv= 2D values to plot with shading
	3404	varcv= 2D values to plot with contours
	3405	vnames= variable names for the figure
	3406	dim[x/y]v = values at the axes of x and y
	3407	dim[x/y]u = units at the axes of x and y
	3408	dimn= dimension names to plot
	3409	colorbar= name of the color bar to use
	3410	ckind= contour kind
	3411	'cmap': as it gets from colorbar
	3412	'fixc,[colname]': fixed color [colname], all stright lines
	3413	'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed line
	3414	clabfmt= format of the labels in the contour plot (None, no labels)
	3415	vs= minmum and maximum values to plot in shadow
	3416	vc= vector with the levels for the contour
	3417	uts= units of the variable [u-shadow, u-contour]
	3418	vtit= title of the variable
	3419	kfig= kind of figure (jpg, pdf, png)
	3420	reva=
	3421	* 'transpose': reverse the axes (x-->y, y-->x)
	3422	* 'flip'@[x/y]: flip the axis x or y
	3423	mapv= map characteristics: [proj],[res]
	3424	see full documentation: http://matplotlib.org/basemap/
	3425	[proj]: projection
	3426	* 'cyl', cilindric
	3427	* 'lcc', lamvbert conformal
	3428	[res]: resolution:
	3429	* 'c', crude
	3430	* 'l', low
	3431	* 'i', intermediate
	3432	* 'h', high
	3433	* 'f', full
	3434	"""
	3435	## import matplotlib as mpl
	3436	## mpl.use('Agg')
	3437	## import matplotlib.pyplot as plt
	3438	fname = 'plot_2D_shadow_contour'
	3439
	3440	if varsv == 'h':
	3441	print fname + '_____________________________________________________________'
	3442	print plot_2D_shadow_contour.__doc__
	3443	quit()
	3444
	3445	if reva[0:4] == 'flip':
	3446	reva0 = 'flip'
	3447	if len(reva.split('@')) != 2:
	3448	print errormsg
	3449	print ' ' + fname + ': flip is given', reva, 'but not axis!'
	3450	quit(-1)
	3451	else:
	3452	reva0 = reva
	3453
	3454	if reva0 == 'transpose':
	3455	print ' reversing the axes of the figure (x-->y, y-->x)!!'
	3456	varsv = np.transpose(varsv)
	3457	varcv = np.transpose(varcv)
	3458	dxv = dimyv
	3459	dyv = dimxv
	3460	dimxv = dxv
	3461	dimyv = dyv
	3462
	3463	if not mapv is None:
	3464	if len(dimxv[:].shape) == 3:
	3465	lon0 = dimxv[0,]
	3466	lat0 = dimyv[0,]
	3467	elif len(dimxv[:].shape) == 2:
	3468	lon0 = dimxv[:]
	3469	lat0 = dimyv[:]
	3470	elif len(dimxv[:].shape) == 1:
	3471	lon00 = dimxv[:]
	3472	lat00 = dimyv[:]
	3473	lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
	3474	lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
	3475
	3476	for iy in range(len(lat00)):
	3477	lon0[iy,:] = lon00
	3478	for ix in range(len(lon00)):
	3479	lat0[:,ix] = lat00
	3480
	3481	map_proj=mapv.split(',')[0]
	3482	map_res=mapv.split(',')[1]
	3483
	3484	dx = lon0.shape[1]
	3485	dy = lon0.shape[0]
	3486
	3487	nlon = lon0[0,0]
	3488	xlon = lon0[dy-1,dx-1]
	3489	nlat = lat0[0,0]
	3490	xlat = lat0[dy-1,dx-1]
	3491
	3492	# Thats too much! :)
	3493	# if lonlatLims is not None:
	3494	# print ' ' + fname + ': cutting the domain to plot !!!!'
	3495	# plt.xlim(lonlatLims[0], lonlatLims[2])
	3496	# plt.ylim(lonlatLims[1], lonlatLims[3])
	3497	# print ' limits: W-E', lonlatLims[0], lonlatLims[2]
	3498	# print ' limits: N-S', lonlatLims[1], lonlatLims[3]
	3499
	3500	# if map_proj == 'cyl':
	3501	# nlon = lonlatLims[0]
	3502	# nlat = lonlatLims[1]
	3503	# xlon = lonlatLims[2]
	3504	# xlat = lonlatLims[3]
	3505	# elif map_proj == 'lcc':
	3506	# lon2 = (lonlatLims[0] + lonlatLims[2])/2.
	3507	# lat2 = (lonlatLims[1] + lonlatLims[3])/2.
	3508	# nlon = lonlatLims[0]
	3509	# xlon = lonlatLims[2]
	3510	# nlat = lonlatLims[1]
	3511	# xlat = lonlatLims[3]
	3512
	3513	lon2 = lon0[dy/2,dx/2]
	3514	lat2 = lat0[dy/2,dx/2]
	3515
	3516	print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \
	3517	xlon, ',', xlat
	3518
	3519	if map_proj == 'cyl':
	3520	m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \
	3521	urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	3522	elif map_proj == 'lcc':
	3523	m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
	3524	llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	3525
	3526	if len(dimxv.shape) == 1:
	3527	lons, lats = np.meshgrid(dimxv, dimyv)
	3528	else:
	3529	if len(dimxv.shape) == 3:
	3530	lons = dimxv[0,:,:]
	3531	lats = dimyv[0,:,:]
	3532	else:
	3533	lons = dimxv[:]
	3534	lats = dimyv[:]
	3535
	3536	x,y = m(lons,lats)
	3537
	3538	else:
	3539	if len(dimxv.shape) == 2:
	3540	x = dimxv
	3541	else:
	3542	if len(dimyv.shape) == 1:
	3543	x = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
	3544	for j in range(len(dimyv)):
	3545	x[j,:] = dimxv
	3546	else:
	3547	x = np.zeros((dimyv.shape), dtype=np.float)
	3548	if x.shape[0] == dimxv.shape[0]:
	3549	for j in range(x.shape[1]):
	3550	x[:,j] = dimxv
	3551	else:
	3552	for j in range(x.shape[0]):
	3553	x[j,:] = dimxv
	3554
	3555	if len(dimyv.shape) == 2:
	3556	y = dimyv
	3557	else:
	3558	if len(dimxv.shape) == 1:
	3559	y = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
	3560	for i in range(len(dimxv)):
	3561	y[:,i] = dimyv
	3562	else:
	3563	y = np.zeros((dimxv.shape), dtype=np.float)
	3564	if y.shape[0] == dimyv.shape[0]:
	3565	for i in range(y.shape[1]):
	3566	y[:,i] = dimyv
	3567	else:
	3568	for i in range(y.shape[0]):
	3569	y[i,:] = dimyv
	3570
[218]	3571	# plt.rc('text', usetex=True)
[192]	3572
	3573	plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(colorbar), vmin=vs[0], vmax=vs[1])
	3574	cbar = plt.colorbar()
	3575
	3576	# contour
	3577	##
	3578	contkind = ckind.split(',')[0]
	3579	if contkind == 'cmap':
	3580	cplot = plt.contour(x, y, varcv, levels=vc)
	3581	elif contkind == 'fixc':
	3582	plt.rcParams['contour.negative_linestyle'] = 'solid'
	3583	coln = ckind.split(',')[1]
	3584	cplot = plt.contour(x, y, varcv, levels=vc, colors=coln)
	3585	elif contkind == 'fixsigc':
	3586	coln = ckind.split(',')[1]
	3587	cplot = plt.contour(x, y, varcv, levels=vc, colors=coln)
	3588	else:
	3589	print errormsg
	3590	print ' ' + fname + ': contour kind "' + contkind + '" not defined !!!!!'
	3591	quit(-1)
	3592
	3593	if clabfmt is not None:
	3594	plt.clabel(cplot, fmt=clabfmt)
	3595	mincntS = format(vc[0], clabfmt[1:len(clabfmt)])
	3596	maxcntS = format(vc[len(vc)-1], clabfmt[1:len(clabfmt)])
	3597	else:
	3598	mincntS = '{:g}'.format(vc[0])
	3599	maxcntS = '{:g}'.format(vc[len(vc)-1])
	3600
	3601	if not mapv is None:
	3602	m.drawcoastlines()
	3603
	3604	meridians = pretty_int(nlon,xlon,5)
	3605	m.drawmeridians(meridians,labels=[True,False,False,True])
	3606	parallels = pretty_int(nlat,xlat,5)
	3607	m.drawparallels(parallels,labels=[False,True,True,False])
	3608
	3609	plt.xlabel('W-E')
	3610	plt.ylabel('S-N')
	3611	else:
	3612	plt.xlabel(variables_values(dimn[1])[0] + ' (' + units_lunits(dimxu) + ')')
	3613	plt.ylabel(variables_values(dimn[0])[0] + ' (' + units_lunits(dimyu) + ')')
	3614
	3615	txpos = pretty_int(x.min(),x.max(),10)
	3616	typos = pretty_int(y.min(),y.max(),10)
	3617	txlabels = list(txpos)
	3618	for i in range(len(txlabels)): txlabels[i] = str(txlabels[i])
	3619	tylabels = list(typos)
	3620	for i in range(len(tylabels)): tylabels[i] = str(tylabels[i])
	3621
	3622	# set the limits of the plot to the limits of the data
	3623	if reva0 == 'flip':
	3624	if reva.split('@')[1] == 'x':
	3625	plt.axis([x.max(), x.min(), y.min(), y.max()])
	3626	else:
	3627	plt.axis([x.min(), x.max(), y.max(), y.min()])
	3628	else:
	3629	plt.axis([x.min(), x.max(), y.min(), y.max()])
	3630
	3631	plt.xticks(txpos, txlabels)
	3632	plt.yticks(typos, tylabels)
	3633
	3634	# units labels
	3635	cbar.set_label(vnames[0].replace('_','\_') + ' (' + units_lunits(uts[0]) + ')')
	3636	plt.annotate(vnames[1].replace('_','\_') +' (' + units_lunits(uts[1]) + ') [' + \
	3637	mincntS + ', ' + maxcntS + ']', xy=(0.55,0.04), xycoords='figure fraction', \
	3638	color=coln)
	3639
	3640	figname = '2Dfields_shadow-contour'
	3641	graphtit = vtit.replace('_','\_').replace('&','\&')
	3642
	3643	plt.title(graphtit)
	3644
	3645	output_kind(kfig, figname, True)
	3646
	3647	return
	3648
	3649	#Nvals=50
	3650	#vals1 = np.zeros((Nvals,Nvals), dtype= np.float)
	3651	#vals2 = np.zeros((Nvals,Nvals), dtype= np.float)
	3652	#for j in range(Nvals):
	3653	# for i in range(Nvals):
	3654	# vals1[j,i]=np.sqrt((j-Nvals/2)2. + (i-Nvals/2)2.)
	3655	# vals2[j,i]=np.sqrt((j-Nvals/2)2. + (i-Nvals/2)2.) - Nvals/2
	3656
	3657	#prettylev=pretty_int(-Nvals/2,Nvals/2,10)
	3658
	3659	#plot_2D_shadow_contour(vals1, vals2, ['var1', 'var2'], np.arange(50)*1., \
	3660	# np.arange(50)*1., ['x-axis','y-axis'], 'rainbow', 'fixc,b', "%.2f", [0, Nvals], \
	3661	# prettylev, ['$ms^{-1}$','$kJm^{-1}s^{-1}$'], 'test var1 & var2', 'pdf', False)
	3662
	3663	def plot_2D_shadow_contour_time(varsv,varcv,vnames,valv,timv,timpos,timlab,valu, \
	3664	timeu,axist,dimn,colorbar,ckind,clabfmt,vs,vc,uts,vtit,kfig,reva,mapv):
	3665	""" Adding labels and other staff to the graph
	3666	varsv= 2D values to plot with shading
	3667	varcv= 2D values to plot with contours
	3668	vnames= variable names for the figure
	3669	valv = values at the axes which is not time
	3670	timv = values for the axis time
	3671	timpos = positions at the axis time
	3672	timlab = labes at the axis time
	3673	valu = units at the axes which is not time
	3674	timeu = units at the axes which is not time
	3675	axist = which is the axis time
	3676	dimn= dimension names to plot
	3677	colorbar= name of the color bar to use
	3678	ckind= contour kind
	3679	'cmap': as it gets from colorbar
	3680	'fixc,[colname]': fixed color [colname], all stright lines
	3681	'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed line
	3682	clabfmt= format of the labels in the contour plot (None, no labels)
	3683	vs= minmum and maximum values to plot in shadow
	3684	vc= vector with the levels for the contour
	3685	uts= units of the variable [u-shadow, u-contour]
	3686	vtit= title of the variable
	3687	kfig= kind of figure (jpg, pdf, png)
	3688	reva=
	3689	* 'transpose': reverse the axes (x-->y, y-->x)
	3690	* 'flip'@[x/y]: flip the axis x or y
	3691	mapv= map characteristics: [proj],[res]
	3692	see full documentation: http://matplotlib.org/basemap/
	3693	[proj]: projection
	3694	* 'cyl', cilindric
	3695	* 'lcc', lamvbert conformal
	3696	[res]: resolution:
	3697	* 'c', crude
	3698	* 'l', low
	3699	* 'i', intermediate
	3700	* 'h', high
	3701	* 'f', full
	3702	"""
	3703	## import matplotlib as mpl
	3704	## mpl.use('Agg')
	3705	## import matplotlib.pyplot as plt
	3706	fname = 'plot_2D_shadow_contour'
	3707
	3708	if varsv == 'h':
	3709	print fname + '_____________________________________________________________'
	3710	print plot_2D_shadow_contour.__doc__
	3711	quit()
	3712
	3713	if axist == 'x':
	3714	dimxv = timv.copy()
	3715	dimyv = valv.copy()
	3716	else:
	3717	dimxv = valv.copy()
	3718	dimyv = timv.copy()
	3719
	3720	if reva[0:4] == 'flip':
	3721	reva0 = 'flip'
	3722	if len(reva.split('@')) != 2:
	3723	print errormsg
	3724	print ' ' + fname + ': flip is given', reva, 'but not axis!'
	3725	quit(-1)
	3726	else:
	3727	reva0 = reva
	3728
	3729	if reva0 == 'transpose':
	3730	if axist == 'x':
	3731	axist = 'y'
	3732	else:
	3733	axist = 'x'
	3734
	3735	if not mapv is None:
	3736	if len(dimxv[:].shape) == 3:
	3737	lon0 = dimxv[0,]
	3738	lat0 = dimyv[0,]
	3739	elif len(dimxv[:].shape) == 2:
	3740	lon0 = dimxv[:]
	3741	lat0 = dimyv[:]
	3742	elif len(dimxv[:].shape) == 1:
	3743	lon00 = dimxv[:]
	3744	lat00 = dimyv[:]
	3745	lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
	3746	lat0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
	3747
	3748	for iy in range(len(lat00)):
	3749	lon0[iy,:] = lon00
	3750	for ix in range(len(lon00)):
	3751	lat0[:,ix] = lat00
	3752	if reva0 == 'transpose':
	3753	print ' reversing the axes of the figure (x-->y, y-->x)!!'
	3754	varsv = np.transpose(varsv)
	3755	varcv = np.transpose(varcv)
	3756	lon0 = np.transpose(lon0)
	3757	lat0 = np.transpose(lat0)
	3758
	3759	map_proj=mapv.split(',')[0]
	3760	map_res=mapv.split(',')[1]
	3761
	3762	dx = lon0.shape[1]
	3763	dy = lon0.shape[0]
	3764
	3765	nlon = lon0[0,0]
	3766	xlon = lon0[dy-1,dx-1]
	3767	nlat = lat0[0,0]
	3768	xlat = lat0[dy-1,dx-1]
	3769
	3770	# Thats too much! :)
	3771	# if lonlatLims is not None:
	3772	# print ' ' + fname + ': cutting the domain to plot !!!!'
	3773	# plt.xlim(lonlatLims[0], lonlatLims[2])
	3774	# plt.ylim(lonlatLims[1], lonlatLims[3])
	3775	# print ' limits: W-E', lonlatLims[0], lonlatLims[2]
	3776	# print ' limits: N-S', lonlatLims[1], lonlatLims[3]
	3777
	3778	# if map_proj == 'cyl':
	3779	# nlon = lonlatLims[0]
	3780	# nlat = lonlatLims[1]
	3781	# xlon = lonlatLims[2]
	3782	# xlat = lonlatLims[3]
	3783	# elif map_proj == 'lcc':
	3784	# lon2 = (lonlatLims[0] + lonlatLims[2])/2.
	3785	# lat2 = (lonlatLims[1] + lonlatLims[3])/2.
	3786	# nlon = lonlatLims[0]
	3787	# xlon = lonlatLims[2]
	3788	# nlat = lonlatLims[1]
	3789	# xlat = lonlatLims[3]
	3790
	3791	lon2 = lon0[dy/2,dx/2]
	3792	lat2 = lat0[dy/2,dx/2]
	3793
	3794	print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \
	3795	xlon, ',', xlat
	3796
	3797	if map_proj == 'cyl':
	3798	m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \
	3799	urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	3800	elif map_proj == 'lcc':
	3801	m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
	3802	llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	3803
	3804	if len(dimxv.shape) == 1:
	3805	lons, lats = np.meshgrid(dimxv, dimyv)
	3806	else:
	3807	if len(dimxv.shape) == 3:
	3808	lons = dimxv[0,:,:]
	3809	lats = dimyv[0,:,:]
	3810	else:
	3811	lons = dimxv[:]
	3812	lats = dimyv[:]
	3813
	3814	x,y = m(lons,lats)
	3815
	3816	else:
	3817	if reva0 == 'transpose':
	3818	print ' reversing the axes of the figure (x-->y, y-->x)!!'
	3819	varsv = np.transpose(varsv)
	3820	varcv = np.transpose(varcv)
	3821	dimn0 = []
	3822	dimn0.append(dimn[1] + '')
	3823	dimn0.append(dimn[0] + '')
	3824	dimn = dimn0
	3825	if len(dimyv.shape) == 2:
	3826	x = np.transpose(dimyv)
	3827	else:
	3828	if len(dimxv.shape) == 2:
	3829	ddx = len(dimyv)
	3830	ddy = dimxv.shape[1]
	3831	else:
	3832	ddx = len(dimyv)
	3833	ddy = len(dimxv)
	3834
	3835	x = np.zeros((ddy,ddx), dtype=np.float)
	3836	for j in range(ddy):
	3837	x[j,:] = dimyv
	3838
	3839	if len(dimxv.shape) == 2:
	3840	y = np.transpose(dimxv)
	3841	else:
	3842	if len(dimyv.shape) == 2:
	3843	ddx = dimyv.shape[0]
	3844	ddy = len(dimxv)
	3845	else:
	3846	ddx = len(dimyv)
	3847	ddy = len(dimxv)
	3848
	3849	y = np.zeros((ddy,ddx), dtype=np.float)
	3850	for i in range(ddx):
	3851	y[:,i] = dimxv
	3852	else:
	3853	if len(dimxv.shape) == 2:
	3854	x = dimxv
	3855	else:
	3856	if len(dimyv.shape) == 1:
	3857	x = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
	3858	for j in range(len(dimyv)):
	3859	x[j,:] = dimxv
	3860	else:
	3861	x = np.zeros((dimyv.shape), dtype=np.float)
	3862	if x.shape[0] == dimxv.shape[0]:
	3863	for j in range(x.shape[1]):
	3864	x[:,j] = dimxv
	3865	else:
	3866	for j in range(x.shape[0]):
	3867	x[j,:] = dimxv
	3868
	3869	if len(dimyv.shape) == 2:
	3870	y = dimyv
	3871	else:
	3872	if len(dimxv.shape) == 1:
	3873	y = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
	3874	for i in range(len(dimxv)):
	3875	y[:,i] = dimyv
	3876	else:
	3877	y = np.zeros((dimxv.shape), dtype=np.float)
	3878	if y.shape[0] == dimyv.shape[0]:
	3879	for i in range(y.shape[1]):
	3880	y[:,i] = dimyv
	3881	else:
	3882	for i in range(y.shape[0]):
	3883	y[i,:] = dimyv
	3884
	3885	dx=varsv.shape[1]
	3886	dy=varsv.shape[0]
	3887
	3888	plt.rc('text', usetex=True)
	3889
	3890	if axist == 'x':
	3891	valpos = pretty_int(y.min(),y.max(),10)
	3892	vallabels = list(valpos)
	3893	for i in range(len(vallabels)): vallabels[i] = str(vallabels[i])
	3894	else:
	3895	valpos = pretty_int(x.min(),x.max(),10)
	3896	vallabels = list(valpos)
	3897	for i in range(len(vallabels)): vallabels[i] = str(vallabels[i])
	3898
	3899	if reva0 == 'flip':
	3900	if reva.split('@')[1] == 'x':
	3901	varsv[:,0:dx-1] = varsv[:,dx-1:0:-1]
	3902	varcv[:,0:dx-1] = varcv[:,dx-1:0:-1]
	3903	plt.xticks(valpos, vallabels[::-1])
	3904	else:
	3905	varsv[0:dy-1,:] = varsv[dy-1:0:-1,:]
	3906	varcv[0:dy-1,:] = varcv[dy-1:0:-1,:]
	3907	plt.yticks(valpos, vallabels[::-1])
	3908	else:
	3909	plt.xlim(0,dx-1)
	3910	plt.ylim(0,dy-1)
	3911
	3912	plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(colorbar), vmin=vs[0], vmax=vs[1])
	3913	cbar = plt.colorbar()
	3914
	3915	# contour
	3916	##
	3917	contkind = ckind.split(',')[0]
	3918	if contkind == 'cmap':
	3919	cplot = plt.contour(x, y, varcv, levels=vc)
	3920	elif contkind == 'fixc':
	3921	plt.rcParams['contour.negative_linestyle'] = 'solid'
	3922	coln = ckind.split(',')[1]
	3923	cplot = plt.contour(x, y, varcv, levels=vc, colors=coln)
	3924	elif contkind == 'fixsigc':
	3925	coln = ckind.split(',')[1]
	3926	cplot = plt.contour(x, y, varcv, levels=vc, colors=coln)
	3927	else:
	3928	print errormsg
	3929	print ' ' + fname + ': contour kind "' + contkind + '" not defined !!!!!'
	3930	quit(-1)
	3931
	3932	if clabfmt is not None:
	3933	plt.clabel(cplot, fmt=clabfmt)
	3934	mincntS = format(vc[0], clabfmt[1:len(clabfmt)])
	3935	maxcntS = format(vc[len(vc)-1], clabfmt[1:len(clabfmt)])
	3936	else:
	3937	mincntS = '{:g}'.format(vc[0])
	3938	maxcntS = '{:g}'.format(vc[len(vc)-1])
	3939
	3940	if not mapv is None:
	3941	m.drawcoastlines()
	3942
	3943	meridians = pretty_int(nlon,xlon,5)
	3944	m.drawmeridians(meridians,labels=[True,False,False,True])
	3945	parallels = pretty_int(nlat,xlat,5)
	3946	m.drawparallels(parallels,labels=[False,True,True,False])
	3947
	3948	plt.xlabel('W-E')
	3949	plt.ylabel('S-N')
	3950	else:
	3951	if axist == 'x':
	3952	plt.xlabel(timeu)
	3953	plt.xticks(timpos, timlab)
	3954	plt.ylabel(variables_values(dimn[0])[0] + ' (' + units_lunits(valu) + ')')
	3955	plt.yticks(valpos, vallabels)
	3956	else:
	3957	plt.xlabel(variables_values(dimn[1])[0] + ' (' + units_lunits(valu) + ')')
	3958	plt.xticks(valpos, vallabels)
	3959	plt.ylabel(timeu)
	3960	plt.yticks(timpos, timlab)
	3961
	3962	# set the limits of the plot to the limits of the data
	3963	plt.axis([x.min(), x.max(), y.min(), y.max()])
	3964
	3965	# units labels
	3966	cbar.set_label(vnames[0].replace('_','\_') + ' (' + units_lunits(uts[0]) + ')')
	3967	plt.annotate(vnames[1].replace('_','\_') +' (' + units_lunits(uts[1]) + ') [' + \
	3968	mincntS + ', ' + maxcntS + ']', xy=(0.55,0.04), xycoords='figure fraction', \
	3969	color=coln)
	3970
	3971	figname = '2Dfields_shadow-contour'
	3972	graphtit = vtit.replace('_','\_').replace('&','\&')
	3973
	3974	plt.title(graphtit)
	3975
	3976	output_kind(kfig, figname, True)
	3977
	3978	return
	3979
	3980	def dxdy_lonlat(dxv,dyv,ddx,ddy):
	3981	""" Function to provide lon/lat 2D lilke-matrices from any sort of dx,dy values
	3982	dxdy_lonlat(dxv,dyv,Lv,lv)
	3983	dx: values for the x
	3984	dy: values for the y
	3985	ddx: ',' list of which dimensions to use from values along x
	3986	ddy: ',' list of which dimensions to use from values along y
	3987	"""
	3988
	3989	fname = 'dxdy_lonlat'
	3990
	3991	if ddx.find(',') > -1:
	3992	dxk = 2
	3993	ddxv = ddx.split(',')
	3994	ddxy = int(ddxv[0])
	3995	ddxx = int(ddxv[1])
	3996	else:
	3997	dxk = 1
	3998	ddxy = int(ddx)
	3999	ddxx = int(ddx)
	4000
	4001	if ddy.find(',') > -1:
	4002	dyk = 2
	4003	ddyv = ddy.split(',')
	4004	ddyy = int(ddyv[0])
	4005	ddyx = int(ddyv[1])
	4006	else:
	4007	dyk = 1
	4008	ddyy = int(ddy)
	4009	ddyx = int(ddy)
	4010
	4011	ddxxv = dxv.shape[ddxx]
	4012	ddxyv = dxv.shape[ddxy]
	4013	ddyxv = dyv.shape[ddyx]
	4014	ddyyv = dyv.shape[ddyy]
	4015
	4016	slicex = []
	4017	if len(dxv.shape) > 1:
	4018	for idim in range(len(dxv.shape)):
	4019	if idim == ddxx or idim == ddxy:
	4020	slicex.append(slice(0,dxv.shape[idim]))
	4021	else:
	4022	slicex.append(0)
	4023	else:
	4024	slicex.append(slice(0,len(dxv)))
	4025
	4026	slicey = []
	4027	if len(dyv.shape) > 1:
	4028	for idim in range(len(dyv.shape)):
	4029	if idim == ddyx or idim == ddyy:
	4030	slicey.append(slice(0,dyv.shape[idim]))
	4031	else:
	4032	slicey.append(0)
	4033	else:
	4034	slicey.append(slice(0,len(dyv)))
	4035
	4036	# print ' ' + fname + ' Lluis shapes dxv:',dxv.shape,'dyv:',dyv.shape
	4037	# print ' ' + fname + ' Lluis slicex:',slicex,'slicey:',slicey
	4038
	4039	if dxk == 2 and dyk == 2:
	4040	if ddxxv != ddyxv:
	4041	print errormsg
	4042	print ' ' + fname + ': wrong dx dimensions! ddxx=',ddxxv,'ddyx=',ddyxv
	4043	print ' choose another for x:',dxv.shape,'or y:',dyv.shape
	4044	quit(-1)
	4045	if ddxyv != ddyyv:
	4046	print errormsg
	4047	print ' ' + fname + ': wrong dy dimensions! ddxy=',ddxyv,'ddyy=',ddyv
	4048	print ' choose another for x:',dxv.shape,'or y:',dyv.shape
	4049	quit(-1)
	4050	dx = ddxxv
	4051	dy = ddxyv
	4052
	4053	print ' ' + fname + ': final dimension 2D lon/lat-like matrices:',dy,',',dx
	4054	lonv = np.zeros((dy,dx), dtype=np.float)
	4055	latv = np.zeros((dy,dx), dtype=np.float)
	4056
	4057
	4058	lonv = dxv[tuple(slicex)]
	4059	latv = dyv[tuple(slicey)]
	4060
	4061	elif dxk == 2 and dyk == 1:
	4062	if not ddxxv == ddyxv and not ddxyv == ddyyv:
	4063	print errormsg
	4064	print ' ' + fname + ': wrong dimensions! ddxx=',ddxxv,'ddyx=',ddyxv, \
	4065	'ddyx=',ddyxv,'ddyy=',ddyyv
	4066	print ' choose another for x:',dxv.shape,'or y:',dyv.shape
	4067	quit(-1)
	4068	dx = ddxvv
	4069	dy = ddxyv
	4070
	4071	print ' ' + fname + ': final dimension 2D lon/lat-like matrices:',dy,',',dx
	4072	lonv = np.zeros((dy,dx), dtype=np.float)
	4073	latv = np.zeros((dy,dx), dtype=np.float)
	4074	lonv = dxv[tuple(slicex)]
	4075
	4076	if ddxxv == ddyxv:
	4077	for iy in range(dy):
	4078	latv[iy,:] = dyv[tuple(slicey)]
	4079	else:
	4080	for ix in range(dx):
	4081	latv[:,ix] = dyv[tuple(slicey)]
	4082
	4083	elif dxk == 1 and dyk == 2:
	4084	if not ddxxv == ddyxv and not ddxyv == ddyyv:
	4085	print errormsg
	4086	print ' ' + fname + ': wrong dimensions! ddxx=',ddxxv,'ddyx=',ddyxv, \
	4087	'ddyx=',ddyxv,'ddyy=',ddyyv
	4088	print ' choose another for x:',dxv.shape,'or y:',dyv.shape
	4089	quit(-1)
	4090	dx = ddyxv
	4091	dy = ddyyv
	4092
	4093	print ' ' + fname + ': final dimension 2D lon/lat-like matrices:',dy,',',dx
	4094	lonv = np.zeros((dy,dx), dtype=np.float)
	4095	latv = np.zeros((dy,dx), dtype=np.float)
	4096
	4097	latv = dyv[tuple(slicey)]
	4098
	4099	if ddyxv == ddxxv:
	4100	for iy in range(dy):
	4101	lonv[iy,:] = dxv[tuple(slicex)]
	4102	else:
	4103	for ix in range(dx):
	4104	lonv[:,ix] = dxv[tuple(slicex)]
	4105
	4106
	4107	elif dxk == 1 and dyk == 1:
	4108	dx = ddxxv
	4109	dy = ddyyv
	4110
	4111	# print 'dx:',dx,'dy:',dy
	4112
	4113	lonv = np.zeros((dy,dx), dtype=np.float)
	4114	latv = np.zeros((dy,dx), dtype=np.float)
	4115
	4116	for iy in range(dy):
	4117	lonv[iy,:] = dxv[tuple(slicex)]
	4118	for ix in range(dx):
	4119	latv[:,ix] = dyv[tuple(slicey)]
	4120
	4121	return lonv,latv
	4122
	4123	def plot_2D_shadow_line(varsv,varlv,vnames,vnamel,dimxv,dimyv,dimxu,dimyu,dimn, \
	4124	colorbar,colln,vs,uts,utl,vtit,kfig,reva,mapv,ifclose):
	4125	""" Plotting a 2D field with shadows and another one with a line
	4126	varsv= 2D values to plot with shading
	4127	varlv= 1D values to plot with line
	4128	vnames= variable names for the shadow variable in the figure
	4129	vnamel= variable names for the line varibale in the figure
	4130	dim[x/y]v = values at the axes of x and y
	4131	dim[x/y]u = units at the axes of x and y
	4132	dimn= dimension names to plot
	4133	colorbar= name of the color bar to use
	4134	colln= color for the line
	4135	vs= minmum and maximum values to plot in shadow
	4136	uts= units of the variable to shadow
	4137	utl= units of the variable to line
	4138	vtit= title of the variable
	4139	kfig= kind of figure (jpg, pdf, png)
	4140	reva=
	4141	* 'transpose': reverse the axes (x-->y, y-->x)
	4142	* 'flip'@[x/y]: flip the axis x or y
	4143	mapv= map characteristics: [proj],[res]
	4144	see full documentation: http://matplotlib.org/basemap/
	4145	[proj]: projection
	4146	* 'cyl', cilindric
	4147	* 'lcc', lambert conformal
	4148	[res]: resolution:
	4149	* 'c', crude
	4150	* 'l', low
	4151	* 'i', intermediate
	4152	* 'h', high
	4153	* 'f', full
	4154	ifclose= boolean value whether figure should be close (finish) or not
	4155	"""
	4156	## import matplotlib as mpl
	4157	## mpl.use('Agg')
	4158	## import matplotlib.pyplot as plt
	4159	fname = 'plot_2D_shadow_line'
	4160
	4161	if varsv == 'h':
	4162	print fname + '_____________________________________________________________'
	4163	print plot_2D_shadow_line.__doc__
	4164	quit()
	4165
	4166	if reva[0:4] == 'flip':
	4167	reva0 = 'flip'
	4168	if len(reva.split('@')) != 2:
	4169	print errormsg
	4170	print ' ' + fname + ': flip is given', reva, 'but not axis!'
	4171	quit(-1)
	4172	else:
	4173	reva0 = reva
	4174
	4175	if reva0 == 'transpose':
	4176	print ' reversing the axes of the figure (x-->y, y-->x)!!'
	4177	varsv = np.transpose(varsv)
	4178	dxv = dimyv
	4179	dyv = dimxv
	4180	dimxv = dxv
	4181	dimyv = dyv
	4182
	4183	if len(dimxv[:].shape) == 3:
	4184	lon0 = dimxv[0,]
	4185	elif len(dimxv[:].shape) == 2:
	4186	lon0 = dimxv[:]
	4187
	4188	if len(dimyv[:].shape) == 3:
	4189	lat0 = dimyv[0,]
	4190	elif len(dimyv[:].shape) == 2:
	4191	lat0 = dimyv[:]
	4192
	4193	if len(dimxv[:].shape) == 1 and len(dimyv[:].shape) == 1:
	4194	lon00 = dimxv[:]
	4195	lon0 = np.zeros( (len(lat00),len(lon00)), dtype=np.float )
	4196
	4197	for iy in range(len(lat00)):
	4198	lon0[iy,:] = lon00
	4199	for ix in range(len(lon00)):
	4200	lat0[:,ix] = lat00
	4201
	4202	if not mapv is None:
	4203	map_proj=mapv.split(',')[0]
	4204	map_res=mapv.split(',')[1]
	4205
	4206	dx = lon0.shape[1]
	4207	dy = lat0.shape[0]
	4208
	4209	nlon = lon0[0,0]
	4210	xlon = lon0[dy-1,dx-1]
	4211	nlat = lat0[0,0]
	4212	xlat = lat0[dy-1,dx-1]
	4213
	4214	# Thats too much! :)
	4215	# if lonlatLims is not None:
	4216	# print ' ' + fname + ': cutting the domain to plot !!!!'
	4217	# plt.xlim(lonlatLims[0], lonlatLims[2])
	4218	# plt.ylim(lonlatLims[1], lonlatLims[3])
	4219	# print ' limits: W-E', lonlatLims[0], lonlatLims[2]
	4220	# print ' limits: N-S', lonlatLims[1], lonlatLims[3]
	4221
	4222	# if map_proj == 'cyl':
	4223	# nlon = lonlatLims[0]
	4224	# nlat = lonlatLims[1]
	4225	# xlon = lonlatLims[2]
	4226	# xlat = lonlatLims[3]
	4227	# elif map_proj == 'lcc':
	4228	# lon2 = (lonlatLims[0] + lonlatLims[2])/2.
	4229	# lat2 = (lonlatLims[1] + lonlatLims[3])/2.
	4230	# nlon = lonlatLims[0]
	4231	# xlon = lonlatLims[2]
	4232	# nlat = lonlatLims[1]
	4233	# xlat = lonlatLims[3]
	4234
	4235	lon2 = lon0[dy/2,dx/2]
	4236	lat2 = lat0[dy/2,dx/2]
	4237
	4238	print 'lon2:', lon2, 'lat2:', lat2, 'SW pt:', nlon, ',', nlat, 'NE pt:', \
	4239	xlon, ',', xlat
	4240
	4241	if map_proj == 'cyl':
	4242	m = Basemap(projection=map_proj, llcrnrlon=nlon, llcrnrlat=nlat, \
	4243	urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	4244	elif map_proj == 'lcc':
	4245	m = Basemap(projection=map_proj, lat_0=lat2, lon_0=lon2, llcrnrlon=nlon, \
	4246	llcrnrlat=nlat, urcrnrlon=xlon, urcrnrlat= xlat, resolution=map_res)
	4247	else:
	4248	print errormsg
	4249	print ' ' + fname + ": map projection '" + map_proj + "' not defined!!!"
	4250	print ' available: cyl, lcc'
	4251	quit(-1)
	4252
	4253	if len(dimxv.shape) == 1:
	4254	lons, lats = np.meshgrid(dimxv, dimyv)
	4255	else:
	4256	if len(dimxv.shape) == 3:
	4257	lons = dimxv[0,:,:]
	4258	else:
	4259	lons = dimxv[:]
	4260
	4261	if len(dimyv.shape) == 3:
	4262	lats = dimyv[0,:,:]
	4263	else:
	4264	lats = dimyv[:]
	4265
	4266	x,y = m(lons,lats)
	4267
	4268	else:
	4269	if len(dimxv.shape) == 3:
	4270	x = dimxv[0,:,:]
	4271	elif len(dimxv.shape) == 2:
	4272	x = dimxv
	4273	else:
	4274	# Attempt of simplier way...
	4275	# x = np.zeros((lon0.shape), dtype=np.float)
	4276	# for j in range(lon0.shape[0]):
	4277	# x[j,:] = dimxv
	4278
	4279	## This way is too complicated and maybe not necessary ? (assuming dimxv.shape == dimyv.shape)
	4280	if len(dimyv.shape) == 1:
	4281	x = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
	4282	for j in range(len(dimxv)):
	4283	x[j,:] = dimxv
	4284	else:
	4285	x = np.zeros((dimyv.shape), dtype=np.float)
	4286	if x.shape[0] == dimxv.shape[0]:
	4287	for j in range(x.shape[1]):
	4288	x[:,j] = dimxv
	4289	else:
	4290	for j in range(x.shape[0]):
	4291	x[j,:] = dimxv
	4292
	4293	if len(dimyv.shape) == 3:
	4294	y = dimyv[0,:,:]
	4295	elif len(dimyv.shape) == 2:
	4296	y = dimyv
	4297	else:
	4298	# y = np.zeros((lat0.shape), dtype=np.float)
	4299	# for i in range(lat0.shape[1]):
	4300	# x[:,i] = dimyv
	4301
	4302	# Idem
	4303	if len(dimxv.shape) == 1:
	4304	y = np.zeros((len(dimyv),len(dimxv)), dtype=np.float)
	4305	for i in range(len(dimxv)):
	4306	y[:,i] = dimyv
	4307	else:
	4308	y = np.zeros((dimxv.shape), dtype=np.float)
	4309	if y.shape[0] == dimyv.shape[0]:
	4310	for i in range(y.shape[1]):
	4311	y[:,i] = dimyv
	4312	else:
	4313	for j in range(y.shape[0]):
	4314	y[j,:] = dimyv
	4315
	4316	plt.rc('text', usetex=True)
	4317
	4318	plt.pcolormesh(x, y, varsv, cmap=plt.get_cmap(colorbar), vmin=vs[0], vmax=vs[1])
	4319	cbar = plt.colorbar()
	4320
	4321	if not mapv is None:
	4322	m.drawcoastlines()
	4323
	4324	meridians = pretty_int(nlon,xlon,5)
	4325	m.drawmeridians(meridians,labels=[True,False,False,True])
	4326	parallels = pretty_int(nlat,xlat,5)
	4327	m.drawparallels(parallels,labels=[False,True,True,False])
	4328
	4329	plt.xlabel('W-E')
	4330	plt.ylabel('S-N')
	4331	else:
	4332	plt.xlabel(variables_values(dimn[1])[0] + ' (' + units_lunits(dimxu) + ')')
	4333	plt.ylabel(variables_values(dimn[0])[0] + ' (' + units_lunits(dimyu) + ')')
	4334
	4335	# Line
	4336	##
	4337
	4338	if reva0 == 'flip' and reva.split('@')[1] == 'y':
	4339	b=-np.max(y[0,:])/np.max(varlv)
	4340	a=np.max(y[0,:])
	4341	else:
	4342	b=np.max(y[0,:])/np.max(varlv)
	4343	a=0.
	4344
	4345	newlinv = varlv*b+a
	4346	if reva0 == 'transpose':
	4347	plt.plot(newlinv, x[0,:], '-', color=colln, linewidth=2)
	4348	else:
	4349	plt.plot(x[0,:], newlinv, '-', color=colln, linewidth=2)
	4350
	4351	txpos = pretty_int(x.min(),x.max(),10)
	4352	typos = pretty_int(y.min(),y.max(),10)
	4353	txlabels = list(txpos)
	4354	for i in range(len(txlabels)): txlabels[i] = str(txlabels[i])
	4355	tylabels = list(typos)
	4356	for i in range(len(tylabels)): tylabels[i] = str(tylabels[i])
	4357
	4358	tllabels = pretty_int(np.min(varlv),np.max(varlv),len(txlabels))
	4359	for it in range(len(tllabels)):
	4360	yval = (tllabels[it]*b+a)
	4361	plt.plot([x.max()*0.97, x.max()], [yval, yval], '-', color='k')
	4362	plt.annotate(tllabels[it], xy=(1.01,tllabels[it]/np.max(varlv)), \
	4363	xycoords='axes fraction')
	4364
	4365	# set the limits of the plot to the limits of the data
	4366	if reva0 == 'flip':
	4367	if reva.split('@')[1] == 'x':
	4368	plt.axis([x.max(), x.min(), y.min(), y.max()])
	4369	else:
	4370	plt.axis([x.min(), x.max(), y.max(), y.min()])
	4371	else:
	4372	plt.axis([x.min(), x.max(), y.min(), y.max()])
	4373
	4374	plt.tick_params(axis='y',right='off')
	4375	if mapv is None:
	4376	plt.xticks(txpos, txlabels)
	4377	plt.yticks(typos, tylabels)
	4378
	4379	tllabels = pretty_int(np.min(varlv),np.max(varlv),len(txlabels))
	4380	for it in range(len(tllabels)):
	4381	plt.annotate(tllabels[it], xy=(1.01,tllabels[it]/np.max(varlv)), xycoords='axes fraction')
	4382
	4383	# units labels
	4384	cbar.set_label(vnames.replace('_','\_') + ' (' + units_lunits(uts) + ')')
	4385
	4386	plt.annotate(vnamel +' (' + units_lunits(utl) + ')', xy=(0.75,0.04),
	4387	xycoords='figure fraction', color=colln)
	4388	figname = '2Dfields_shadow_line'
	4389	graphtit = vtit.replace('_','\_').replace('&','\&')
	4390
	4391	plt.title(graphtit)
	4392
	4393	output_kind(kfig, figname, ifclose)
	4394
	4395	return
	4396
	4397	def plot_Neighbourghood_evol(varsv, dxv, dyv, vnames, ttits, tpos, tlabels, colorbar, \
	4398	Nng, vs, uts, gtit, kfig, ifclose):
	4399	""" Plotting neighbourghood evolution
	4400	varsv= 2D values to plot with shading
	4401	vnames= shading variable name for the figure
	4402	d[x/y]v= values at the axes of x and y
	4403	ttits= titles of both time axis
	4404	tpos= positions of the time ticks
	4405	tlabels= labels of the time ticks
	4406	colorbar= name of the color bar to use
	4407	Nng= Number of grid points of the full side of the box (odd value)
	4408	vs= minmum and maximum values to plot in shadow or:
	4409	'Srange': for full range
	4410	'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
	4411	'Saroundminmax@val': for minval,maxval
	4412	'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
	4413	percentile_(100-val)-median)
	4414	'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min@val,max@val-mean)
	4415	'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
	4416	'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
	4417	percentile_(100-val)-median)
	4418	uts= units of the variable to shadow
	4419	gtit= title of the graph
	4420	kfig= kind of figure (jpg, pdf, png)
	4421	ifclose= boolean value whether figure should be close (finish) or not
	4422	"""
	4423	import numpy.ma as ma
	4424
	4425	fname = 'plot_Neighbourghood_evol'
	4426
	4427	if varsv == 'h':
	4428	print fname + '_____________________________________________________________'
	4429	print plot_Neighbourghood_evol.__doc__
	4430	quit()
	4431
	4432	if len(varsv.shape) != 2:
	4433	print errormsg
	4434	print ' ' + fname + ': wrong number of dimensions of the values: ', \
	4435	varsv.shape
	4436	quit(-1)
	4437
	4438	varsvmask = ma.masked_equal(varsv,fillValue)
	4439
	4440	vsend = np.zeros((2), dtype=np.float)
	4441	# Changing limits of the colors
	4442	if type(vs[0]) != type(np.float(1.)):
	4443	if vs[0] == 'Srange':
	4444	vsend[0] = np.min(varsvmask)
	4445	elif vs[0][0:11] == 'Saroundmean':
	4446	meanv = np.mean(varsvmask)
	4447	permean = np.float(vs[0].split('@')[1])
	4448	minv = np.min(varsvmask)*permean
	4449	maxv = np.max(varsvmask)*permean
	4450	minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
	4451	vsend[0] = meanv-minextrm
	4452	vsend[1] = meanv+minextrm
	4453	elif vs[0][0:13] == 'Saroundminmax':
	4454	permean = np.float(vs[0].split('@')[1])
	4455	minv = np.min(varsvmask)*permean
	4456	maxv = np.max(varsvmask)*permean
	4457	vsend[0] = minv
	4458	vsend[1] = maxv
	4459	elif vs[0][0:17] == 'Saroundpercentile':
	4460	medianv = np.median(varsvmask)
	4461	valper = np.float(vs[0].split('@')[1])
	4462	minv = np.percentile(varsvmask, valper)
	4463	maxv = np.percentile(varsvmask, 100.-valper)
	4464	minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
	4465	vsend[0] = medianv-minextrm
	4466	vsend[1] = medianv+minextrm
	4467	elif vs[0][0:5] == 'Smean':
	4468	meanv = np.mean(varsvmask)
	4469	permean = np.float(vs[0].split('@')[1])
	4470	minv = np.min(varsvmask)*permean
	4471	maxv = np.max(varsvmask)*permean
	4472	minextrm = np.min([np.abs(meanv-minv), np.abs(maxv-meanv)])
	4473	vsend[0] = -minextrm
	4474	vsend[1] = minextrm
	4475	elif vs[0][0:7] == 'Smedian':
	4476	medianv = np.median(varsvmask)
	4477	permedian = np.float(vs[0].split('@')[1])
	4478	minv = np.min(varsvmask)*permedian
	4479	maxv = np.max(varsvmask)*permedian
	4480	minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
	4481	vsend[0] = -minextrm
	4482	vsend[1] = minextrm
	4483	elif vs[0][0:11] == 'Spercentile':
	4484	medianv = np.median(varsvmask)
	4485	valper = np.float(vs[0].split('@')[1])
	4486	minv = np.percentile(varsvmask, valper)
	4487	maxv = np.percentile(varsvmask, 100.-valper)
	4488	minextrm = np.min([np.abs(medianv-minv), np.abs(maxv-medianv)])
	4489	vsend[0] = -minextrm
	4490	vsend[1] = minextrm
	4491	else:
	4492	print errormsg
	4493	print ' ' + fname + ": range '" + vs[0] + "' not ready!!!"
	4494	quit(-1)
	4495	print ' ' + fname + ': modified shadow min,max:',vsend
	4496	else:
	4497	vsend[0] = vs[0]
	4498
	4499	if type(vs[0]) != type(np.float(1.)):
	4500	if vs[1] == 'range':
	4501	vsend[1] = np.max(varsv)
	4502	else:
	4503	vsend[1] = vs[1]
	4504
	4505	plt.rc('text', usetex=True)
	4506
	4507	# plt.pcolormesh(dxv, dyv, varsv, cmap=plt.get_cmap(colorbar), vmin=vsend[0], vmax=vsend[1])
	4508	plt.pcolormesh(varsvmask, cmap=plt.get_cmap(colorbar), vmin=vsend[0], vmax=vsend[1])
	4509	cbar = plt.colorbar()
	4510
	4511	newtposx = (tpos[0][:] - np.min(dxv)) * len(dxv) * Nng / (np.max(dxv) - np.min(dxv))
	4512	newtposy = (tpos[1][:] - np.min(dyv)) * len(dyv) * Nng / (np.max(dyv) - np.min(dyv))
	4513
	4514	plt.xticks(newtposx, tlabels[0])
	4515	plt.yticks(newtposy, tlabels[1])
	4516	plt.xlabel(ttits[0])
	4517	plt.ylabel(ttits[1])
	4518
	4519	plt.axes().set_aspect('equal')
	4520	# From: http://stackoverflow.com/questions/14406214/moving-x-axis-to-the-top-of-a-plot-in-matplotlib
	4521	plt.axes().xaxis.tick_top
	4522	plt.axes().xaxis.set_ticks_position('top')
	4523
	4524	# units labels
	4525	cbar.set_label(vnames.replace('_','\_') + ' (' + units_lunits(uts) + ')')
	4526
	4527	figname = 'Neighbourghood_evol'
	4528	graphtit = gtit.replace('_','\_').replace('&','\&')
	4529
	4530	plt.title(graphtit, position=(0.5,1.05))
	4531
	4532	output_kind(kfig, figname, ifclose)
	4533
	4534	return
	4535
[218]	4536	def plot_lines(vardv, varvv, vaxis, dtit, linesn, vtit, vunit, gtit, gloc, kfig):
[192]	4537	""" Function to plot a collection of lines
	4538	vardv= list of set of dimension values
	4539	varvv= list of set of values
	4540	vaxis= which axis will be used for the values ('x', or 'y')
	4541	dtit= title for the common dimension
	4542	linesn= names for the legend
	4543	vtit= title for the vaxis
	4544	vunit= units of the vaxis
	4545	gtit= main title
[218]	4546	gloc= location of the legend (-1, autmoatic)
	4547	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
	4548	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
	4549	9: 'upper center', 10: 'center'
[192]	4550	kfig= kind of figure
	4551	plot_lines([np.arange(10)], [np.sin(np.arange(10)*np.pi/2.5)], 'y', 'time (s)', \
	4552	['2.5'], 'sin', '-', 'sinus frequency dependency', 'pdf')
	4553	"""
	4554	fname = 'plot_lines'
	4555
	4556	if vardv == 'h':
	4557	print fname + '_____________________________________________________________'
	4558	print plot_lines.__doc__
	4559	quit()
	4560
	4561	# Canging line kinds every 7 lines (end of standard colors)
	4562	linekinds=['.-','x-','o-']
	4563
	4564	Ntraj = len(vardv)
	4565
	4566	N7lines = 0
	4567
	4568	plt.rc('text', usetex=True)
	4569
	4570	if vaxis == 'x':
	4571	for il in range(Ntraj):
	4572	plt.plot(varvv[il], vardv[il], linekinds[N7lines], label= linesn[il])
	4573	if il == 6: N7lines = N7lines + 1
	4574
	4575	plt.xlabel(vtit + ' (' + vunit + ')')
	4576	plt.ylabel(dtit)
	4577	plt.xlim(np.min(varvv[:]),np.max(varvv[:]))
	4578	plt.ylim(np.min(vardv[:]),np.max(vardv[:]))
	4579
	4580	else:
	4581	for il in range(Ntraj):
	4582	plt.plot(vardv[il], varvv[il], linekinds[N7lines], label= linesn[il])
	4583	if il == 6: N7lines = N7lines + 1
	4584
	4585	plt.xlabel(dtit)
	4586	plt.ylabel(vtit + ' (' + vunit + ')')
	4587	plt.xlim(np.min(vardv[:]),np.max(vardv[:]))
	4588	plt.ylim(np.min(varvv[:]),np.max(varvv[:]))
	4589
	4590	figname = 'lines'
	4591	graphtit = gtit
	4592
	4593	plt.title(graphtit)
[218]	4594	plt.legend(loc=gloc)
[192]	4595
	4596	output_kind(kfig, figname, True)
	4597
	4598	return
[194]	4599
	4600	def plot_lines_time(vardv, varvv, vaxis, dtit, linesn, vtit, vunit, tpos, tlabs, \
[204]	4601	gtit, gloc, kfig):
[194]	4602	""" Function to plot a collection of lines with a time axis
	4603	vardv= list of set of dimension values
	4604	varvv= list of set of values
	4605	vaxis= which axis will be used for the time values ('x', or 'y')
	4606	dtit= title for the common dimension
	4607	linesn= names for the legend
	4608	vtit= title for the vaxis
	4609	vunit= units of the vaxis
	4610	tpos= positions of the time ticks
	4611	tlabs= labels of the time ticks
	4612	gtit= main title
[204]	4613	gloc= location of the legend (-1, autmoatic)
	4614	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
	4615	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
	4616	9: 'upper center', 10: 'center'
[194]	4617	kfig= kind of figure
	4618	plot_lines([np.arange(10)], [np.sin(np.arange(10)*np.pi/2.5)], 'y', 'time (s)', \
	4619	['2.5'], 'sin', '-', 'sinus frequency dependency', 'pdf')
	4620	"""
	4621	fname = 'plot_lines'
	4622
	4623	if vardv == 'h':
	4624	print fname + '_____________________________________________________________'
	4625	print plot_lines.__doc__
	4626	quit()
	4627
	4628	# Canging line kinds every 7 lines (end of standard colors)
	4629	linekinds=['.-','x-','o-']
	4630
	4631	Ntraj = len(vardv)
	4632
	4633	N7lines = 0
	4634
	4635	plt.rc('text', usetex=True)
[196]	4636	varTvv = []
	4637	varTdv = []
[194]	4638
	4639	if vaxis == 'x':
	4640	for il in range(Ntraj):
	4641	plt.plot(varvv[il], vardv[il], linekinds[N7lines], label= linesn[il])
[197]	4642	varTvv = varTvv + list(varvv[il])
	4643	varTdv = varTdv + list(vardv[il])
[194]	4644	if il == 6: N7lines = N7lines + 1
	4645
	4646	plt.xlabel(vtit + ' (' + vunit + ')')
	4647	plt.ylabel(dtit)
[197]	4648	plt.xlim(np.min(varTvv),np.max(varTvv))
	4649	plt.ylim(np.min(varTdv),np.max(varTdv))
[194]	4650	plt.yticks(tpos, tlabs)
	4651	else:
	4652	for il in range(Ntraj):
	4653	plt.plot(vardv[il], varvv[il], linekinds[N7lines], label= linesn[il])
[197]	4654	varTvv = varTvv + list(varvv[il])
	4655	varTdv = varTdv + list(vardv[il])
[194]	4656	if il == 6: N7lines = N7lines + 1
	4657
	4658	plt.xlabel(dtit)
	4659	plt.ylabel(vtit + ' (' + vunit + ')')
[197]	4660
	4661	plt.xlim(np.min(varTdv),np.max(varTdv))
	4662	plt.ylim(np.min(varTvv),np.max(varTvv))
[194]	4663	plt.xticks(tpos, tlabs)
	4664
	4665	figname = 'lines_time'
	4666	graphtit = gtit
	4667
	4668	plt.title(graphtit)
[204]	4669	plt.legend(loc=gloc)
[194]	4670
	4671	output_kind(kfig, figname, True)
	4672
	4673	return

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

Download in other formats: