Context Navigation

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

Last change on this file since 652 was 652, checked in by lfita, 9 years ago
Adding `plot_ptZvals' to plot lon,lat stations at a given time Slicing also dimensions on the '2D_shad' plot Adding monotone removal on `slice_variable'
File size: 205.2 KB

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

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

Download in other formats: