Context Navigation

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

Last change on this file since 211 was 208, checked in by lfita, 10 years ago
Chaging from to '-' in 'units_lunits'
File size: 160.1 KB

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

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

Download in other formats: