Context Navigation

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

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

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

Download in other formats: