Context Navigation

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

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

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

Download in other formats: