Context Navigation

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

Last change on this file since 592 was 586, checked in by lfita, 10 years ago
Adding 360 deg transformation for the position of the points on `draw_points'
File size: 197.3 KB

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

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

Download in other formats: