Context Navigation

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

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

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

Download in other formats: