Context Navigation

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

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

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

Download in other formats: