Context Navigation

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

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

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

Download in other formats: