Context Navigation

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

Last change on this file since 498 was 489, checked in by lfita, 10 years ago
Adding use of `dxdy_lonlatDIMS' on 'shad_contour' Removing some residual debug printings...
File size: 188.1 KB

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

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

Download in other formats: