Context Navigation

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

Last change on this file since 372 was 372, checked in by lfita, 10 years ago
Removing of printing of values in 'variables_values'
File size: 180.4 KB

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

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

Download in other formats: