Context Navigation

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

Last change on this file since 628 was 625, checked in by lfita, 10 years ago
Adding radial averaged plots
File size: 198.1 KB

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

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

Download in other formats: