Context Navigation

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

Last change on this file since 338 was 328, checked in by lfita, 10 years ago
Adding kg kg-1 s-1 to 'units_lunits'
File size: 176.5 KB

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

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

Download in other formats: