Context Navigation

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

Last change on this file since 245 was 238, checked in by lfita, 10 years ago
Removing all '_' in variables' names on `variables_values'
File size: 164.3 KB

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

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

Download in other formats: