Context Navigation

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

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

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

Download in other formats: