Context Navigation

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

Last change on this file since 232 was 228, checked in by lfita, 10 years ago
Adding 'QSNOW','QGRAUPEL','QICE','QHAIL' in 'variables_values' Fixing 'QCLOUD' in 'variables_values'
File size: 162.0 KB

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

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

Download in other formats: