Context Navigation

source: lmdz_wrf/trunk/tools/drawing.py @ 470

Last change on this file since 470 was 470, checked in by lfita, 10 years ago
Adding multiple dimvar values on `draw_lines'
File size: 95.3 KB

Line
1	import numpy as np
2	import os
3	from netCDF4 import Dataset as NetCDFFile
4	import drawing_tools as drw
5	from optparse import OptionParser
6	import sys
7	from cStringIO import StringIO
8
9	## e.g. # drawing.py -f /media/data2/etudes/WRF_LMDZ/WL_HyMeX/IIphase/medic950116/wlmdza/wrfout/wrfout_d01_1995-01-13_00:00:00 -o create_movie -S 'draw_2D_shad#Time@WRFTimes@10@95@191@1#tas:East_West\|-1,North_South\|-1,Time\|2:longitude:latitude:Summer:270.,300.:tas\|at\|t=0:pdf:None:None' -v T2
10	## e.g. # drawing.py -f wrfout_d01_1980-03-01_00\:00\:00_Time_B0-E48-I1.nc -o draw_2D_shad -S 'tas:East_West\|-1,North_South\|-1,Time\|2:longitude:latitude:Summer:270.,300.:tas\|at\|t=0:pdf:None:cyl,i' -v T2
11	## e.g. # drawing.py -f ~/etudes/domains/MEDCORDEX/geo_em.d01.nc -o draw_2D_shad_cont -S 'landmask,height:Time\|0:Time\|0:XLONG_M:XLAT_M:rainbow:fixc,k:%.2f:0,1:0.,3000.,10:landmask & height:pdf:False:lcc,i' -v LANDMASK,HGT_M
12	## e.g. # drawing.py -f ~/etudes/domains/MEDCORDEX/geo_em.d01.nc -o draw_2D_shad_cont -S 'height,landmask:Time\|0:Time\|0:XLONG_M:XLAT_M:terrain:fixc,k:None:0.,3000.:0,1,10:MEDCORDEX height & landmask:pdf:False:lcc,i' -v HGT_M,LANDMASK
13	## e.g. # drawing.py -o draw_2D_shad_line -f 'mean_dtcon-pluc-pres_lat.nc,mean_dtcon-pluc-pres_lat.nc' -S 'dtcon,prc:bottom_top\|-1,south_north\|-1:latmean:presmean:seismic,k:-5.,5.:monthly\|dtcon\|&\|prc:pdf:flip@y:None:True' -v 'dtconmean,prcmean'
14	## e.g. # drawing.py -f 'geo_em.d02.nc' -o draw_topo_geogrid -S '0.,3000.:None:FF_3dom d02:png:cyl,i'
15	## e.g. # drawing.py -o draw_topo_geogrid_boxes -S '0.,3000.:None:FF domains:pdf:lcc,i:d01,d02,d03' -f geo_em.d01.nc,geo_em.d02.nc,geo_em.d03.nc
16	## e.g. # drawing.py -o draw_trajectories -f 'WRF/control/trajectory.dat@-1@/home/lluis/etudes/domains/WL_HyMeX_HighRes_C/geo_em.d03.nc#XLONG_M#XLAT_M,WRF_LMDZ/wlmdza/trajectory.dat@-1@/home/lluis/etudes/domains/WL_HyMeX_C/geo_em.d01.nc#XLONG_M#XLAT_M,WRF_LMDZ/wlmdzb/trajectory.dat@-1@/home/lluis/etudes/domains/WL_HyMeX_C/geo_em.d01.nc#XLONG_M#XLAT_M,WRF_LMDZ/wlmdzb_ii/trajectory.dat@-1@/home/lluis/etudes/domains/WL_HyMeX_C/geo_em.d01.nc#XLONG_M#XLAT_M' -S '$WRF_{CRM}$,$LMDZ_{AR4.0}$,$LMDZ_{NPv3.1}$,$LMDZ_{NPv3.1b}$\|None\|medicane trajectories\|pdf\|cyl,i' -v obs/trajectory.dat,satellite,-1
17	## e.g. # drawing.py -o draw_vals_trajectories -f WRF_LMDZ/wlmdza/tevolboxtraj_T2.nc,WRF_LMDZ/wlmdzb/tevolboxtraj_T2.nc,WRF/control/tevolboxtraj_T2.nc -S 'mean:-1:$WRF_{CRM}$,$LMDZ_{AR4.0}$,$LMDZ_{NPv3.1}$@4:tas:time\|($[DD]^[HH]$):exct,6,h:$%d^{%H}$:trajectory\|following\|mean:pdf' -v T2
18	## e.g. # drawing.py -o draw_2D_shad_time -f 'netcdf_concatenated.nc' -S 'dtcon:Time\|-1,bottom_top\|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly\|dtcon:pdf:transpose:time\|hours!since!1949-12-01\|exct,2,d\|%d\|date!([DD])\|x:True' -v 'dtconmean'
19	## e.g. # drawing.py -o variable_values -S PSFC
20	## e.g. # drawing.py -o draw_timeSeries -f wrfout_d01_1979-12-01_00:00:00_bottom_top_B6-E6-I1_south_north_B3-E3-I1_west_east_B26-E26-I1.nc -S 'dt_con:time\|($[DD]^{[HH]}$):exct,12,h:$%d^{%H}$:time\|evolution\|at\|-1\|6\|3\|26:1:pdf' -v 'LDQCON,time'
21	## e.g. # drawing.py -f wrfout_d01_1979-12-01_00:00:00 -o draw_Neighbourghood_evol -S 'q:Time\|-1\|Times,bottom_top\|6\|ZNU,south_north\|3\|XLAT,west_east\|26\|XLONG:south_north,west_east:5:auto:time\|($[DD]^{[HH]}$),time\|($[DD]^{[HH]}$):exct,2,h\|exct,1,d:$%d^{%H}$,$%d^{%H}$:5\|pts\|neighbourghood\|temporal\|evolution:0.0,0.004:BuPu:pdf:True' -v QVAPOR
22	## e.g. # drawing.py -o draw_lines_time -f wrfout_d01_1980-03-01_00:00:00_Time_B0-E48-I1_south_north_B15-E15-I1_west_east_B15-E15-I1.nc -S 'time;y;time ([DD]${[HH]}$);file1;tas;evolution;time\|hours!since!1949-12-01_00:00:00\|exct,12,h\|%d$^{%H}$;pdf' -v T2
23
24	main = 'drawing.py'
25
26	errormsg = 'ERROR -- error -- ERROR -- error'
27	warnmsg = 'WARNING -- waring -- WARNING -- warning'
28	fillValue=1.e20
29
30	namegraphics = ['create_movie', 'draw_2D_shad', 'draw_2D_shad_time', \
31	'draw_2D_shad_cont', 'draw_2D_shad_cont_time', 'draw_2D_shad_line', \
32	'draw_2D_shad_line_time', 'draw_timeSeries', 'draw_topo_geogrid', \
33	'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories', \
34	'draw_lines', 'draw_lines_time', 'draw_Neighbourghood_evol', 'list_graphics', \
35	'variable_values']
36
37	def draw_2D_shad(ncfile, values, varn):
38	""" plotting a fields with shading
39	draw_2D_shad(ncfile, values, varn)
40	ncfile= file to use
41	values=[vnamefs]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
42	[kindfig]:[reverse]:[mapv]:[close]
43	[vnamefs]: Name in the figure of the variable to be shaded
44	[dimvals]: ',' list of [dimname]\|[value] telling at which dimension of the
45	variable a given value is required (-1, all the length)
46	[dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
47	[colorbar]: name of the color bar
48	[smin/axv]: minimum and maximum value for the shading or:
49	'Srange': for full range
50	'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
51	'Saroundminmax@val': for minval,maxval
52	'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
53	percentile_(100-val)-median)
54	'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min@val,max@val-mean)
55	'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
56	'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
57	percentile_(100-val)-median)
58	[figt]: title of the figure ('\|' for spaces)
59	[kindfig]: kind of figure
60	[reverse]: Transformation of the values
61	* 'transpose': reverse the axes (x-->y, y-->x)
62	* 'flip'@[x/y]: flip the axis x or y
63	[mapv]: map characteristics: [proj],[res]
64	see full documentation: http://matplotlib.org/basemap/
65	[proj]: projection
66	* 'cyl', cilindric
67	* 'lcc', lamvbert conformal
68	[res]: resolution:
69	* 'c', crude
70	* 'l', low
71	* 'i', intermediate
72	* 'h', high
73	* 'f', full
74	valules= 'rh:z\|-1,x\|-1:z\|-1,x\|-1:lat:pressure:BuPu:0.,100.:rh:pdf:flip@y:None'
75	varn= [varsn] name of the variable to plot with shading
76	"""
77
78	fname = 'draw_2D_shad'
79	if values == 'h':
80	print fname + '_____________________________________________________________'
81	print draw_2D_shad.__doc__
82	quit()
83
84	expectargs = ['[vnamefs]', '[dimvals]', '[dimxvn]', '[dimyvn]', '[colorbar]', \
85	'[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', '[mapv]', '[close]']
86
87	drw.check_arguments(fname,len(expectargs),values,':',expectargs)
88
89	vnamesfig = values.split(':')[0]
90	dimvals= values.split(':')[1].replace('\|',':')
91	vdimxn = values.split(':')[2]
92	vdimyn = values.split(':')[3]
93	colbarn = values.split(':')[4]
94	shadminmax = values.split(':')[5]
95	figtitle = values.split(':')[6].replace('\|',' ')
96	figkind = values.split(':')[7]
97	revals = values.split(':')[8]
98	mapvalue = values.split(':')[9]
99	# varn = values.split(':')[10]
100
101	ncfiles = ncfile
102
103	if not os.path.isfile(ncfiles):
104	print errormsg
105	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
106	quit(-1)
107
108	objsf = NetCDFFile(ncfiles, 'r')
109
110	varns = varn.split(',')[0]
111
112	if not objsf.variables.has_key(varns):
113	print errormsg
114	print ' ' + fname + ': shading file "' + ncfiles + \
115	'" does not have variable "' + varns + '" !!'
116	quit(-1)
117
118	# Variables' values
119	objvars = objsf.variables[varns]
120
121	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
122
123	# Dimensions names
124	## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
125	## dimnamesv = []
126	## for idd in range(len(objvars.dimensions)):
127	## cutdim = False
128	## for idc in range(len(dimvals.split(','))):
129	## dimcutn = dimvals.split(',')[idc].split(':')[0]
130	## print objvars.dimensions[idd], dimcutn
131	## if objvars.dimensions[idd] == dimcutn:
132	## cutdim = True
133	## break
134	## if not cutdim: dimnamesv.append(objvars.dimensions[idd])
135	dimnamesv = [vdimyn, vdimxn]
136
137	if drw.searchInlist(objvars.ncattrs(),'units'):
138	varunits = objvars.getncattr('units')
139	else:
140	print warnmsg
141	print ' ' + fname + ": variable '" + varn + "' without units!!"
142	varunits = '-'
143
144	if not objsf.variables.has_key(vdimxn):
145	print errormsg
146	print ' ' + fname + ': shading file "' + ncfiles + \
147	'" does not have dimension variable "' + vdimxn + '" !!'
148	quit(-1)
149	if not objsf.variables.has_key(vdimyn):
150	print errormsg
151	print ' ' + fname + ': shading file "' + ncfiles + \
152	'" does not have dimension variable "' + vdimyn + '" !!'
153	quit(-1)
154
155	objdimx = objsf.variables[vdimxn]
156	objdimy = objsf.variables[vdimyn]
157	if drw.searchInlist(objdimx.ncattrs(),'units'):
158	odimxu = objdimx.getncattr('units')
159	else:
160	print warnmsg
161	print ' ' + fname + ": variable dimension '" + vdimxn + "' without units!!"
162	odimxu = '-'
163
164	if drw.searchInlist(objdimy.ncattrs(),'units'):
165	odimyu = objdimy.getncattr('units')
166	else:
167	print warnmsg
168	print ' ' + fname + ": variable dimension '" + vdimyn + "' without units!!"
169	odimyu = '-'
170
171	if len(objdimx.shape) <= 2:
172	# odimxv = objdimx[valshad.shape]
173	# odimyv = objdimy[valshad.shape]
174	odimxv = objdimx[:]
175	odimyv = objdimy[:]
176
177	elif len(objdimx.shape) == 3:
178	# dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])]
179	# odimxv = objdimx[tuple(dimcut)]
180	# odimyv = objdimy[tuple(dimcut)]
181	odimxv = objdimx[0,:]
182	odimyv = objdimy[0,:]
183	else:
184	print errormsg
185	print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
186	' not ready!!'
187	quit(-1)
188
189	shading_nx = []
190	if shadminmax.split(',')[0][0:1] != 'S':
191	shading_nx.append(np.float(shadminmax.split(',')[0]))
192	else:
193	shading_nx.append(shadminmax.split(',')[0])
194
195	if shadminmax.split(',')[1][0:1] != 'S':
196	shading_nx.append(np.float(shadminmax.split(',')[1]))
197	else:
198	shading_nx.append(shadminmax.split(',')[1])
199
200	if mapvalue == 'None': mapvalue = None
201
202	drw.plot_2D_shadow(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu, dimnamesv,\
203	colbarn, shading_nx, varunits, figtitle, figkind, revals, mapvalue, True)
204
205	return
206
207	def draw_2D_shad_time(ncfile, values, varn):
208	""" plotting a fields with shading with time values
209	draw_2D_shad(ncfile, values, varn)
210	ncfile= file to use
211	values=[vnamefs]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
212	[kindfig]:[reverse]:[timevals]:[close]
213	[vnamefs]: Name in the figure of the variable to be shaded
214	[dimvals]: ',' list of [dimname]\|[value] telling at which dimension of the
215	variable a given value is required (-1, all the length)
216	[dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
217	[colorbar]: name of the color bar
218	[smin/axv]: minimum and maximum value for the shading or:
219	'Srange': for full range
220	'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
221	'Saroundminmax@val': for minval,maxval
222	'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
223	percentile_(100-val)-median)
224	'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min@val,max@val-mean)
225	'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
226	'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
227	percentile_(100-val)-median)
228	[figt]: title of the figure ('\|' for spaces)
229	[kindfig]: kind of figure
230	[reverse]: Transformation of the values
231	* 'transpose': reverse the axes (x-->y, y-->x)
232	* 'flip'@[x/y]: flip the axis x or y
233	[timevals]: [timen]\|[units]\|[kind]\|[tfmt]\|[label]\|[timeaxis] time labels characteristics
234	[timen]; name of the time variable
235	[units]; units string according to CF conventions ([tunits] since
236	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
237	[kind]; kind of output
238	'Nval': according to a given number of values as 'Nval',[Nval]
239	'exct': according to an exact time unit as 'exct',[tunit];
240	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
241	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
242	'l': milisecond
243	[tfmt]; desired format
244	[label]; label at the graph ('!' for spaces)
245	[close]: should figure be closed (finished)
246	values='dtcon:Time\|-1,bottom_top\|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly\|'
247	'dtcon:pdf:transpose:time\|hours!since!1949-12-01\|exct,2,d\|%d\|date!([DD])\|x:True
248	varn= [varsn] name of the variable to plot with shading
249	"""
250	fname = 'draw_2D_shad_time'
251
252	if values == 'h':
253	print fname + '_____________________________________________________________'
254	print draw_2D_shad_time.__doc__
255	quit()
256
257	farguments = ['[vnamefs]', '[dimvals]', '[dimxvn]', '[dimyvn]', '[colorbar]', \
258	'[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', '[timevals]', '[close]']
259	drw.check_arguments(fname,11,values,':',farguments)
260
261	vnamesfig = values.split(':')[0]
262	dimvals= values.split(':')[1].replace('\|',':')
263	vdimxn = values.split(':')[2]
264	vdimyn = values.split(':')[3]
265	colbarn = values.split(':')[4]
266	shadminmax = values.split(':')[5]
267	figtitle = values.split(':')[6].replace('\|',' ')
268	figkind = values.split(':')[7]
269	revals = values.split(':')[8]
270	timevals = values.split(':')[9]
271	close = values.split(':')[10]
272
273	ncfiles = ncfile
274
275	if not os.path.isfile(ncfiles):
276	print errormsg
277	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
278	quit(-1)
279
280	objsf = NetCDFFile(ncfiles, 'r')
281
282	varns = varn.split(',')[0]
283
284	if not objsf.variables.has_key(varns):
285	print errormsg
286	print ' ' + fname + ': shading file "' + ncfiles + \
287	'" does not have variable "' + varns + '" !!'
288	quit(-1)
289
290	# Variables' values
291	objvars = objsf.variables[varns]
292
293	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
294
295	dimnamesv = [vdimyn, vdimxn]
296
297	varunits = objvars.getncattr('units')
298
299	if not objsf.variables.has_key(vdimxn):
300	print errormsg
301	print ' ' + fname + ': shading file "' + ncfiles + \
302	'" does not have dimension variable "' + vdimxn + '" !!'
303	quit(-1)
304	if not objsf.variables.has_key(vdimyn):
305	print errormsg
306	print ' ' + fname + ': shading file "' + ncfiles + \
307	'" does not have dimensino variable "' + vdimyn + '" !!'
308	quit(-1)
309
310	objdimx = objsf.variables[vdimxn]
311	objdimy = objsf.variables[vdimyn]
312	odimxu = objdimx.getncattr('units')
313	odimyu = objdimy.getncattr('units')
314
315	if len(objdimx.shape) <= 2:
316	odimxv = objdimx[:]
317	odimyv = objdimy[:]
318
319	elif len(objdimx.shape) == 3:
320	odimxv = objdimx[0,:]
321	odimyv = objdimy[0,:]
322	else:
323	print errormsg
324	print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
325	' not ready!!'
326	quit(-1)
327
328	timename = timevals.split('\|')[0]
329	timeunit = timevals.split('\|')[1].replace('!',' ')
330	timekind = timevals.split('\|')[2]
331	timefmt = timevals.split('\|')[3]
332	timelabel = timevals.split('\|')[4].replace('!',' ')
333
334	if vdimxn == timename:
335	odimxv = objsf.variables[vdimxn][:]
336	odimxu = timelabel
337	timeaxis = 'x'
338	odimyv = objsf.variables[vdimyn]
339	odimyu = odimyv.getncattr('units')
340	timepos, timelabels = drw.CFtimes_plot(odimxv, timeunit, timekind, timefmt)
341	elif vdimyn == timename:
342	odimyv = objsf.variables[vdimyn][:]
343	odimyu = timelabel
344	timeaxis = 'y'
345	odimxv = objsf.variables[vdimxn]
346	odimxu = odimxv.getncattr('units')
347	timepos, timelabels = drw.CFtimes_plot(odimyv, timeunit, timekind, timefmt)
348	else:
349	print errormsg
350	print ' ' + fname + ": time variable '" + timename + "' not found!!"
351	quit(-1)
352
353	shading_nx = []
354	if shadminmax.split(',')[0][0:1] != 'S':
355	shading_nx.append(np.float(shadminmax.split(',')[0]))
356	else:
357	shading_nx.append(shadminmax.split(',')[0])
358
359	if shadminmax.split(',')[1][0:1] != 'S':
360	shading_nx.append(np.float(shadminmax.split(',')[1]))
361	else:
362	shading_nx.append(shadminmax.split(',')[1])
363
364	closeval = drw.Str_Bool(close)
365
366	drw.plot_2D_shadow_time(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu, \
367	dimnamesv, colbarn, shading_nx, varunits, figtitle, figkind, revals, timeaxis, \
368	timepos, timelabels, closeval)
369
370	return
371
372	def draw_2D_shad_cont(ncfile, values, varn):
373	""" plotting two fields, one with shading and the other with contour lines
374	draw_2D_shad_cont(ncfile, values, varn)
375	ncfile= [ncfilevars],[ncfilevarc] files to use (one value, same file)
376	values=[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:[figt]:[kindfig]:[reverse]:[mapv]
377	[vnamefs],[vnamefc]: Name in the figure of the shaded and the contour variables
378	[dimvals/c]: list of [dimname]\|[value] telling at which dimension of the
379	variable a given value is required (no dimension name, all the length)
380	[dimx/yvn]: names of the variables with the values of the dimensions for the plot
381	[colorbar]: name of the color bar
382	[ckind]: kind of contours
383	'cmap': as it gets from colorbar
384	'fixc,[colname]': fixed color [colname], all stright lines
385	'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed line
386	[clabfmt]: format of the labels in the contour (None, also possible)
387	[smin/axv]: minimum and maximum value for the shading
388	[sminc]:[smaxv]:[Nlev]: minimum, maximum and number of values for the contour
389	[figt]: title of the figure ('\|' for spaces)
390	[kindfig]: kind of figure
391	[reverse]: does the values be transposed? 'True/False',
392	[mapv]: map characteristics: [proj],[res]
393	see full documentation: http://matplotlib.org/basemap/
394	[proj]: projection
395	* 'cyl', cilindric
396	* 'lcc', lamvbert conformal
397	[res]: resolution:
398	* 'c', crude
399	* 'l', low
400	* 'i', intermediate
401	* 'h', high
402	* 'f', full
403	valules= 'rh,ta:z\|-1,x\|-1:z\|-1,x\|-1:lat:pressure:BuPu:fixsigc,black:%d:0.,100.:195.,305.,7:Meridonal\|average\|of\|rh\|&\|ta:pdf:flip@y:None'
404	varn= [varsn],[varcn] name of the variable to plot with shading variable with contour
405	"""
406
407	fname = 'draw_2D_shad_cont'
408	if values == 'h':
409	print fname + '_____________________________________________________________'
410	print draw_2D_shad_cont.__doc__
411	quit()
412
413	vnamesfig = values.split(':')[0].split(',')
414	dimvals= values.split(':')[1].replace('\|',':')
415	dimvalc= values.split(':')[2].replace('\|',':')
416	vdimxn = values.split(':')[3]
417	vdimyn = values.split(':')[4]
418	colbarn = values.split(':')[5]
419	countkind = values.split(':')[6]
420	countlabelfmt = values.split(':')[7]
421	shadminmax = values.split(':')[8]
422	contlevels = values.split(':')[9]
423	figtitle = values.split(':')[10].replace('\|',' ')
424	figkind = values.split(':')[11]
425	revals = values.split(':')[12]
426	mapvalue = values.split(':')[13]
427
428	if2filenames = ncfile.find(',')
429
430	if if2filenames != -1:
431	ncfiles = ncfile.split(',')[0]
432	ncfilec = ncfile.split(',')[1]
433	else:
434	ncfiles = ncfile
435	ncfilec = ncfile
436
437	if not os.path.isfile(ncfiles):
438	print errormsg
439	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
440	quit(-1)
441
442	if not os.path.isfile(ncfilec):
443	print errormsg
444	print ' ' + fname + ': contour file "' + ncfilec + '" does not exist !!'
445	quit(-1)
446
447	objsf = NetCDFFile(ncfiles, 'r')
448	objcf = NetCDFFile(ncfilec, 'r')
449
450	varns = varn.split(',')[0]
451	varnc = varn.split(',')[1]
452
453	if not objsf.variables.has_key(varns):
454	print errormsg
455	print ' ' + fname + ': shading file "' + ncfiles + \
456	'" does not have variable "' + varns + '" !!'
457	quit(-1)
458
459	if not objcf.variables.has_key(varnc):
460	print errormsg
461	print ' ' + fname + ': contour file "' + ncfilec + \
462	'" does not have variable "' + varnc + '" !!'
463	quit(-1)
464
465	# Variables' values
466	objvars = objsf.variables[varns]
467	objvarc = objcf.variables[varnc]
468
469	if len(objvars.shape) != len(objvarc.shape):
470	print errormsg
471	print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \
472	objvars.shape, 'different than contour variable "' + varnc + '": ', \
473	objvarc.shape,' !!!'
474	quit(-1)
475
476	for idim in range(len(objvars.shape)):
477	if objvars.shape[idim] != objvarc.shape[idim]:
478	print errormsg
479	print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \
480	objvars.shape, 'different than contour variable "' + varnc + '": ', \
481	objvarc.shape,' !!!'
482	quit(-1)
483
484	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
485	valcont, dimscont = drw.slice_variable(objvarc, dimvalc.replace(',','\|'))
486
487	# Dimensions names
488	## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
489	## dimnamesv = []
490	## for idd in range(len(objvars.dimensions)):
491	## cutdim = False
492	## for idc in range(len(dimvals.split(','))):
493	## dimcutn = dimvals.split(',')[idc].split(':')[0]
494	## print objvars.dimensions[idd], dimcutn
495	## if objvars.dimensions[idd] == dimcutn:
496	## cutdim = True
497	## break
498	## if not cutdim: dimnamesv.append(objvars.dimensions[idd])
499	dimnamesv = [vdimyn, vdimxn]
500
501	varunits = []
502	varunits.append(objvars.getncattr('units'))
503	varunits.append(objvarc.getncattr('units'))
504
505	if not objsf.variables.has_key(vdimxn):
506	print errormsg
507	print ' ' + fname + ': shading file "' + ncfiles + \
508	'" does not have dimension variable "' + vdimxn + '" !!'
509	quit(-1)
510	if not objsf.variables.has_key(vdimyn):
511	print errormsg
512	print ' ' + fname + ': shading file "' + ncfiles + \
513	'" does not have dimensino variable "' + vdimyn + '" !!'
514	quit(-1)
515
516	objdimx = objsf.variables[vdimxn]
517	objdimy = objsf.variables[vdimyn]
518	odimxu = objdimx.getncattr('units')
519	odimyu = objdimy.getncattr('units')
520
521	# Getting only that dimensions with coincident names
522	dimnvx = objdimx.dimensions
523	cutslice = []
524	for idimn in objdimx.dimensions:
525	found = False
526	for dimsn in dimsshad:
527	if idimn == dimsn:
528	cutslice.append(slice(0,len(objsf.dimensions[idimn])))
529	found = True
530	if not found: cutslice.append(0)
531
532	odimxv = objdimx[tuple(cutslice)]
533
534	dimnvy = objdimy.dimensions
535	cutslice = []
536	for idimn in objdimy.dimensions:
537	found = False
538	for dimsn in dimsshad:
539	if idimn == dimsn:
540	cutslice.append(slice(0,len(objsf.dimensions[idimn])))
541	found = True
542	if not found: cutslice.append(0)
543
544	odimyv = objdimy[tuple(cutslice)]
545
546	# if len(objdimx.shape) <= 2:
547	# odimxv = objdimx[:]
548	# odimyv = objdimy[:]
549	# elif len(objdimx.shape) == 3:
550	# odimxv = objdimx[0,:]
551	# odimyv = objdimy[0,:]
552	# else:
553	# print errormsg
554	# print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
555	# ' not ready!!'
556	# quit(-1)
557
558	if countlabelfmt == 'None':
559	countlfmt = None
560	else:
561	countlfmt = countlabelfmt
562
563	shading_nx = np.zeros((2), dtype=np.float)
564	shading_nx[0] = np.float(shadminmax.split(',')[0])
565	shading_nx[1] = np.float(shadminmax.split(',')[1])
566
567	clevmin = np.float(contlevels.split(',')[0])
568	clevmax = np.float(contlevels.split(',')[1])
569	Nclevels = int(contlevels.split(',')[2])
570
571	levels_cont = drw.pretty_int(clevmin, clevmax, Nclevels)
572
573	if len(levels_cont) <= 1:
574	print warnmsg
575	print ' ' + fname + ': wrong contour levels:', levels_cont,' !!'
576	del(levels_cont)
577	levels_cont = np.zeros((Nclevels), dtype=np.float)
578	levels_cont = np.arange(7)*(clevmax - clevmin)/(Nclevels-1)
579	print ' generating default ones: ',levels_cont
580
581	if mapvalue == 'None': mapvalue = None
582
583	drw.plot_2D_shadow_contour(valshad, valcont, vnamesfig, odimxv, odimyv, odimxu, \
584	odimyu, dimnamesv, colbarn, countkind, countlfmt, shading_nx, levels_cont, \
585	varunits, figtitle, figkind, revals, mapvalue)
586
587	return
588
589	def draw_2D_shad_cont_time(ncfile, values, varn):
590	""" plotting two fields, one with shading and the other with contour lines being
591	one of the dimensions of time characteristics
592	draw_2D_shad_cont(ncfile, values, varn)
593	ncfile= [ncfilevars],[ncfilevarc] files to use (one value, same file)
594	values=[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:[figt]:[kindfig]:[reverse]:[timevals]:[mapv]
595	[vnamefs],[vnamefc]: Name in the figure of the shaded and the contour variables
596	[dimvals/c]: list of [dimname]\|[value] telling at which dimension of the
597	variable a given value is required (no dimension name, all the length)
598	[dimx/yvn]: ',' list with the name of the variables with the values of the dimensions
599	[colorbar]: name of the color bar
600	[ckind]: kind of contours
601	'cmap': as it gets from colorbar
602	'fixc,[colname]': fixed color [colname], all stright lines
603	'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed line
604	[clabfmt]: format of the labels in the contour (None, also possible)
605	[smin/axv]: minimum and maximum value for the shading
606	[sminc]:[smaxv]:[Nlev]: minimum, maximum and number of values for the contour
607	[figt]: title of the figure ('\|' for spaces)
608	[kindfig]: kind of figure
609	[reverse]: modification to the dimensions:
610	'transposed': transpose matrices
611	'flip',[x/y]: flip only the dimension [x] or [y]
612	[timevals]: [timen]\|[units]\|[kind]\|[tfmt]\|[label] time labels characteristics
613	[timen]; name of the time variable
614	[units]; units string according to CF conventions ([tunits] since
615	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
616	[kind]; kind of output
617	'Nval': according to a given number of values as 'Nval',[Nval]
618	'exct': according to an exact time unit as 'exct',[tunit];
619	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
620	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
621	'l': milisecond
622	[tfmt]; desired format
623	[label]; label at the graph ('!' for spaces)
624	[mapv]: map characteristics: [proj],[res]
625	see full documentation: http://matplotlib.org/basemap/
626	[proj]: projection
627	* 'cyl', cilindric
628	* 'lcc', lamvbert conformal
629	[res]: resolution:
630	* 'c', crude
631	* 'l', low
632	* 'i', intermediate
633	* 'h', high
634	* 'f', full
635	valules= 'rh,ta:z\|-1,x\|-1:z\|-1,x\|-1:lat:pressure:BuPu:fixsigc,black:%d:0.,100.:195.,305.,7:Meridonal\|average\|of\|rh\|&\|ta:pdf:flip@y:time!hours!since!1949/12/01\|exct,5d\|%d\|date!([DD]):None'
636	varn= [varsn],[varcn] name of the variable to plot with shading variable with contour
637	"""
638
639	fname = 'draw_2D_shad_cont_time'
640	if values == 'h':
641	print fname + '_____________________________________________________________'
642	print draw_2D_shad_cont_time.__doc__
643	quit()
644
645	expectargs = ['[vnamefs]', '[dimvals]', '[dimvalc]', '[dimxvn]', '[dimyvn]', \
646	'[colorbar]', '[ckind]', '[clabfmt]', '[sminv],[smaxv]', \
647	'[sminc],[smaxv],[Nlev]', '[figt]', '[kindfig]', '[reverse]', '[timevals]', \
648	'[mapv]']
649
650	drw.check_arguments(fname,len(expectargs),values,':',expectargs)
651
652	vnamesfig = values.split(':')[0].split(',')
653	dimvals= values.split(':')[1].replace('\|',':')
654	dimvalc= values.split(':')[2].replace('\|',':')
655	vdimxn = values.split(':')[3]
656	vdimyn = values.split(':')[4]
657	colbarn = values.split(':')[5]
658	countkind = values.split(':')[6]
659	countlabelfmt = values.split(':')[7]
660	shadminmax = values.split(':')[8]
661	contlevels = values.split(':')[9]
662	figtitle = values.split(':')[10].replace('\|',' ')
663	figkind = values.split(':')[11]
664	revals = values.split(':')[12]
665	timevals = values.split(':')[13]
666	mapvalue = values.split(':')[14]
667
668	if2filenames = ncfile.find(',')
669
670	if if2filenames != -1:
671	ncfiles = ncfile.split(',')[0]
672	ncfilec = ncfile.split(',')[1]
673	else:
674	ncfiles = ncfile
675	ncfilec = ncfile
676
677	if not os.path.isfile(ncfiles):
678	print errormsg
679	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
680	quit(-1)
681
682	if not os.path.isfile(ncfilec):
683	print errormsg
684	print ' ' + fname + ': contour file "' + ncfilec + '" does not exist !!'
685	quit(-1)
686
687	objsf = NetCDFFile(ncfiles, 'r')
688	objcf = NetCDFFile(ncfilec, 'r')
689
690	varns = varn.split(',')[0]
691	varnc = varn.split(',')[1]
692
693	if not objsf.variables.has_key(varns):
694	print errormsg
695	print ' ' + fname + ': shading file "' + ncfiles + \
696	'" does not have variable "' + varns + '" !!'
697	quit(-1)
698
699	if not objcf.variables.has_key(varnc):
700	print errormsg
701	print ' ' + fname + ': contour file "' + ncfilec + \
702	'" does not have variable "' + varnc + '" !!'
703	quit(-1)
704
705	# Variables' values
706	objvars = objsf.variables[varns]
707	objvarc = objcf.variables[varnc]
708
709	if len(objvars.shape) != len(objvarc.shape):
710	print errormsg
711	print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \
712	objvars.shape, 'different than contour variable "' + varnc + '": ', \
713	objvarc.shape,' !!!'
714	quit(-1)
715
716	for idim in range(len(objvars.shape)):
717	if objvars.shape[idim] != objvarc.shape[idim]:
718	print errormsg
719	print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \
720	objvars.shape, 'different than contour variable "' + varnc + '": ', \
721	objvarc.shape,' !!!'
722	quit(-1)
723
724	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
725	valcont, dimscont = drw.slice_variable(objvarc, dimvalc.replace(',','\|'))
726
727	# Dimensions names
728	## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
729	## dimnamesv = []
730	## for idd in range(len(objvars.dimensions)):
731	## cutdim = False
732	## for idc in range(len(dimvals.split(','))):
733	## dimcutn = dimvals.split(',')[idc].split(':')[0]
734	## print objvars.dimensions[idd], dimcutn
735	## if objvars.dimensions[idd] == dimcutn:
736	## cutdim = True
737	## break
738	## if not cutdim: dimnamesv.append(objvars.dimensions[idd])
739	dimnamesv = [vdimyn, vdimxn]
740
741	varunits = []
742	varunits.append(objvars.getncattr('units'))
743	varunits.append(objvarc.getncattr('units'))
744
745	if not objsf.variables.has_key(vdimxn):
746	print errormsg
747	print ' ' + fname + ': shading file "' + ncfiles + \
748	'" does not have dimension variable "' + vdimxn + '" !!'
749	quit(-1)
750	if not objsf.variables.has_key(vdimyn):
751	print errormsg
752	print ' ' + fname + ': shading file "' + ncfiles + \
753	'" does not have dimension variable "' + vdimyn + '" !!'
754	quit(-1)
755
756	timename = timevals.split('\|')[0]
757	timeunit = timevals.split('\|')[1].replace('!',' ')
758	timekind = timevals.split('\|')[2]
759	timefmt = timevals.split('\|')[3]
760	timelabel = timevals.split('\|')[4].replace('!',' ')
761
762	if vdimxn == timename:
763	timevals = objsf.variables[vdimxn][:]
764	timedims = objsf.variables[vdimxn].dimensions
765	dimt = 'x'
766	ovalaxis = objsf.variables[vdimyn]
767	ovalu = ovalaxis.getncattr('units')
768	elif vdimyn == timename:
769	timevals = objsf.variables[vdimyn][:]
770	timedims = objsf.variables[vdimyn].dimensions
771	dimt = 'y'
772	ovalaxis = objsf.variables[vdimxn]
773	ovalu = ovalaxis.getncattr('units')
774	else:
775	print errormsg
776	print ' ' + fname + ": time variable '" + timename + "' not found!!"
777	quit(-1)
778
779	timepos, timelabels = drw.CFtimes_plot(timevals, timeunit, timekind, timefmt)
780
781	# Getting only that dimensions with coincident names
782	dimnvx = ovalaxis.dimensions
783
784	cutslice = []
785	for idimn in dimsshad:
786	found = False
787	for dimsn in dimnvx:
788	if idimn == dimsn:
789	cutslice.append(slice(0,len(objsf.dimensions[idimn])))
790	found = True
791	if not found: cutslice.append(0)
792
793	ovalaxisv = ovalaxis[tuple(cutslice)]
794
795	## if len(ovalaxis.shape) <= 2:
796	## ovalaxisv = ovalaxis[:]
797
798	## elif len(ovalaxis.shape) == 3:
799	## ovalaxisv = ovalaxis[0,:]
800	## else:
801	## print errormsg
802	## print ' ' + fname + ': shape of dimension variable:', ovalaxis.shape, \
803	## ' not ready!!'
804	## quit(-1)
805
806	if countlabelfmt == 'None':
807	countlfmt = None
808	else:
809	countlfmt = countlabelfmt
810
811	shading_nx = np.zeros((2), dtype=np.float)
812	shading_nx[0] = np.float(shadminmax.split(',')[0])
813	shading_nx[1] = np.float(shadminmax.split(',')[1])
814
815	clevmin = np.float(contlevels.split(',')[0])
816	clevmax = np.float(contlevels.split(',')[1])
817	Nclevels = int(contlevels.split(',')[2])
818
819	levels_cont = drw.pretty_int(clevmin, clevmax, Nclevels)
820
821	if len(levels_cont) <= 1:
822	print warnmsg
823	print ' ' + fname + ': wrong contour levels:', levels_cont,' !!'
824	del(levels_cont)
825	levels_cont = np.zeros((Nclevels), dtype=np.float)
826	levels_cont = np.arange(7)*(clevmax - clevmin)/(Nclevels-1)
827	print ' generating default ones: ',levels_cont
828
829	if mapvalue == 'None': mapvalue = None
830
831	drw.plot_2D_shadow_contour_time(valshad, valcont, vnamesfig, ovalaxisv, \
832	timevals, timepos, timelabels, ovalu, timelabel, dimt, dimnamesv, colbarn, \
833	countkind, countlfmt, shading_nx, levels_cont, varunits, figtitle, figkind, \
834	revals, mapvalue)
835
836	return
837
838	def draw_2D_shad_line(ncfile, values, varn):
839	""" plotting a fields with shading and another with line
840	draw_2D_shad_line(ncfile, values, varn)
841	ncfile= [ncfiles],[ncfilel] file to use for the shading and for the line
842	values=[vnamefs],[vnamefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar],[colline]:[sminv],[smaxv]:[figt]:
843	[kindfig]:[reverse]:[mapv]:[close]
844	[vnamefs]: Name in the figure of the variable to be shaded
845	[vnamefl]: Name in the figure of the variable to be lined
846	[dimvals]: ',' list of [dimname]\|[value] telling at which dimension of the
847	variable a given value is required (-1, all the length)
848	[dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
849	[colorbar]: name of the color bar
850	[colline]: name of the color for the line
851	[smin/axv]: minimum and maximum value for the shading
852	[figt]: title of the figure ('\|' for spaces)
853	[kindfig]: kind of figure
854	[reverse]: Transformation of the values
855	* 'transpose': reverse the axes (x-->y, y-->x)
856	* 'flip'@[x/y]: flip the axis x or y
857	[mapv]: map characteristics: [proj],[res]
858	see full documentation: http://matplotlib.org/basemap/
859	[proj]: projection
860	* 'cyl', cilindric
861	* 'lcc', lamvbert conformal
862	[res]: resolution:
863	* 'c', crude
864	* 'l', low
865	* 'i', intermediate
866	* 'h', high
867	* 'f', full
868	valules= 'rh:z\|-1,x\|-1:z\|-1,x\|-1:lat:pressure:BuPu:0.,100.:rh:pdf:flip@y:None'
869	varn= [varsn],[varnl] name of the variable to plot with shading and with line
870	"""
871
872	fname = 'draw_2D_shad_line'
873	if values == 'h':
874	print fname + '_____________________________________________________________'
875	print draw_2D_shad_line.__doc__
876	quit()
877
878	farguments = ['[vnamefs],[vnamefl]', '[dimvals]', '[dimxvn]', '[dimyvn]', \
879	'[colorbar],[colline]', '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', \
880	'[mapv]', '[close]']
881	drw.check_arguments(fname,11,values,':',farguments)
882
883	vnamesfig = values.split(':')[0].split(',')[0]
884	dimvals= values.split(':')[1].replace('\|',':')
885	vdimxn = values.split(':')[2]
886	vdimyn = values.split(':')[3]
887	colbarn = values.split(':')[4].split(',')[0]
888	shadminmax = values.split(':')[5]
889	figtitle = values.split(':')[6].replace('\|',' ')
890	figkind = values.split(':')[7]
891	revals = values.split(':')[8]
892	mapvalue = values.split(':')[9]
893	# varn = values.split(':')[10]
894
895	ncfiles = ncfile.split(',')[0]
896
897	if not os.path.isfile(ncfiles):
898	print errormsg
899	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
900	quit(-1)
901
902	objsf = NetCDFFile(ncfiles, 'r')
903
904	varns = varn.split(',')[0]
905
906	if not objsf.variables.has_key(varns):
907	print errormsg
908	print ' ' + fname + ': shading file "' + ncfiles + \
909	'" does not have variable "' + varns + '" !!'
910	quit(-1)
911
912	# Variables' values
913	objvars = objsf.variables[varns]
914
915	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
916
917	# Dimensions names
918	## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
919	## dimnamesv = []
920	## for idd in range(len(objvars.dimensions)):
921	## cutdim = False
922	## for idc in range(len(dimvals.split(','))):
923	## dimcutn = dimvals.split(',')[idc].split(':')[0]
924	## print objvars.dimensions[idd], dimcutn
925	## if objvars.dimensions[idd] == dimcutn:
926	## cutdim = True
927	## break
928	## if not cutdim: dimnamesv.append(objvars.dimensions[idd])
929	dimnamesv = [vdimyn, vdimxn]
930
931	varunits = objvars.getncattr('units')
932
933	if not objsf.variables.has_key(vdimxn):
934	print errormsg
935	print ' ' + fname + ': shading file "' + ncfiles + \
936	'" does not have dimension variable "' + vdimxn + '" !!'
937	quit(-1)
938	if not objsf.variables.has_key(vdimyn):
939	print errormsg
940	print ' ' + fname + ': shading file "' + ncfiles + \
941	'" does not have dimensino variable "' + vdimyn + '" !!'
942	quit(-1)
943
944	objdimx = objsf.variables[vdimxn]
945	objdimy = objsf.variables[vdimyn]
946	odimxu = objdimx.getncattr('units')
947	odimyu = objdimy.getncattr('units')
948
949	if len(objdimx.shape) <= 2:
950	# odimxv = objdimx[valshad.shape]
951	# odimyv = objdimy[valshad.shape]
952	odimxv = objdimx[:]
953	odimyv = objdimy[:]
954
955	elif len(objdimx.shape) == 3:
956	# dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])]
957	# odimxv = objdimx[tuple(dimcut)]
958	# odimyv = objdimy[tuple(dimcut)]
959	odimxv = objdimx[0,:]
960	odimyv = objdimy[0,:]
961	else:
962	print errormsg
963	print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
964	' not ready!!'
965	quit(-1)
966
967	shading_nx = np.zeros((2), dtype=np.float)
968	shading_nx[0] = np.float(shadminmax.split(',')[0])
969	shading_nx[1] = np.float(shadminmax.split(',')[1])
970
971	if mapvalue == 'None': mapvalue = None
972
973	# line plot
974	##
975	ncfilel = ncfile.split(',')[1]
976	vnamelfig = values.split(':')[0].split(',')[1]
977	varnl = varn.split(',')[1]
978	colline = values.split(':')[4].split(',')[1]
979
980	objlf = NetCDFFile(ncfilel,'r')
981	objlvar = objlf.variables[varnl]
982
983	linevals = objlvar[:]
984	varlunits = objlvar.units
985
986	drw.plot_2D_shadow_line(valshad, linevals, vnamesfig, vnamelfig, odimxv, odimyv, \
987	odimxu, odimyu, dimnamesv, colbarn, colline, shading_nx, varunits, varlunits, \
988	figtitle, figkind, revals, mapvalue, True)
989
990	objsf.close()
991	objlf.close()
992
993	return
994
995	def draw_2D_shad_line_time(ncfile, values, varn):
996	""" plotting a fields with shading and a line with time values
997	draw_2D_shad_line(ncfile, values, varn)
998	ncfile= [ncfiles],[ncfilel] files to use to draw with shading and the line
999	values= [vnamefs],[vanemefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
1000	[kindfig]:[reverse]:[timevals]:[close]
1001	[vnamefs]: Name in the figure of the variable to be shaded
1002	[vnamefl]: Name in the figure of the variable to be lined
1003	[dimvals]: ',' list of [dimname]\|[value] telling at which dimension of the
1004	variable a given value is required (-1, all the length)
1005	[dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
1006	[colorbar]: name of the color bar
1007	[smin/axv]: minimum and maximum value for the shading
1008	[figt]: title of the figure ('\|' for spaces)
1009	[kindfig]: kind of figure
1010	[reverse]: Transformation of the values
1011	* 'transpose': reverse the axes (x-->y, y-->x)
1012	* 'flip'@[x/y]: flip the axis x or y
1013	[timevals]: [timen]\|[units]\|[kind]\|[tfmt]\|[label]\|[timeaxis] time labels characteristics
1014	[timen]; name of the time variable
1015	[units]; units string according to CF conventions ([tunits] since
1016	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
1017	[kind]; kind of output
1018	'Nval': according to a given number of values as 'Nval',[Nval]
1019	'exct': according to an exact time unit as 'exct',[tunit];
1020	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1021	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1022	'l': milisecond
1023	[tfmt]; desired format
1024	[label]; label at the graph ('!' for spaces)
1025	[close]: should figure be closed (finished)
1026	values='dtcon,prc:Time\|-1,bottom_top\|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly\|'
1027	'dtcon:pdf:transpose:time\|hours!since!1949-12-01\|exct,2,d\|%d\|date!([DD])\|x:True
1028	varn= [varsn].[varln] name of the variable to plot with shading and to plot with line
1029	"""
1030	fname = 'draw_2D_shad_line_time'
1031	if values == 'h':
1032	print fname + '_____________________________________________________________'
1033	print draw_2D_shad__line_time.__doc__
1034	quit()
1035
1036	farguments = ['[vnamefs],[vanemefl]', '[dimvals]', '[dimxvn]', '[dimyvn]', \
1037	'[colorbar]', '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', \
1038	'[timevals]', '[close]']
1039	drw.check_arguments(fname,11,values,':',farguments)
1040
1041	vnamesfig = values.split(':')[0].split(',')[0]
1042	dimvals= values.split(':')[1].replace('\|',':')
1043	vdimxn = values.split(':')[2]
1044	vdimyn = values.split(':')[3]
1045	colbarn = values.split(':')[4]
1046	shadminmax = values.split(':')[5]
1047	figtitle = values.split(':')[6].replace('\|',' ')
1048	figkind = values.split(':')[7]
1049	revals = values.split(':')[8]
1050	timevals = values.split(':')[9]
1051	close = values.split(':')[10]
1052
1053	ncfiles = ncfile.split(',')[0]
1054
1055	if not os.path.isfile(ncfiles):
1056	print errormsg
1057	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
1058	quit(-1)
1059
1060	objsf = NetCDFFile(ncfiles, 'r')
1061
1062	varns = varn.split(',')[0]
1063
1064	if not objsf.variables.has_key(varns):
1065	print errormsg
1066	print ' ' + fname + ': shading file "' + ncfiles + \
1067	'" does not have variable "' + varns + '" !!'
1068	quit(-1)
1069
1070	# Variables' values
1071	objvars = objsf.variables[varns]
1072
1073	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
1074
1075	dimnamesv = [vdimyn, vdimxn]
1076
1077	varunits = objvars.getncattr('units')
1078
1079	if not objsf.variables.has_key(vdimxn):
1080	print errormsg
1081	print ' ' + fname + ': shading file "' + ncfiles + \
1082	'" does not have dimension variable "' + vdimxn + '" !!'
1083	quit(-1)
1084	if not objsf.variables.has_key(vdimyn):
1085	print errormsg
1086	print ' ' + fname + ': shading file "' + ncfiles + \
1087	'" does not have dimensino variable "' + vdimyn + '" !!'
1088	quit(-1)
1089
1090	objdimx = objsf.variables[vdimxn]
1091	objdimy = objsf.variables[vdimyn]
1092	odimxu = objdimx.getncattr('units')
1093	odimyu = objdimy.getncattr('units')
1094
1095	if len(objdimx.shape) <= 2:
1096	odimxv = objdimx[:]
1097	odimyv = objdimy[:]
1098
1099	elif len(objdimx.shape) == 3:
1100	odimxv = objdimx[0,:]
1101	odimyv = objdimy[0,:]
1102	else:
1103	print errormsg
1104	print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
1105	' not ready!!'
1106	quit(-1)
1107
1108	timename = timevals.split('\|')[0]
1109	timeunit = timevals.split('\|')[1].replace('!',' ')
1110	timekind = timevals.split('\|')[2]
1111	timefmt = timevals.split('\|')[3]
1112	timelabel = timevals.split('\|')[4].replace('!',' ')
1113
1114	if vdimxn == timename:
1115	odimxv = objsf.variables[vdimxn][:]
1116	odimxu = timelabel
1117	timeaxis = 'x'
1118	odimyv = objsf.variables[vdimyn]
1119	odimyu = odimyv.getncattr('units')
1120	timepos, timelabels = drw.CFtimes_plot(odimxv, timeunit, timekind, timefmt)
1121	elif vdimyn == timename:
1122	odimyv = objsf.variables[vdimyn][:]
1123	odimyu = timelabel
1124	timeaxis = 'y'
1125	odimxv = objsf.variables[vdimxn]
1126	odimxu = odimxv.getncattr('units')
1127	timepos, timelabels = drw.CFtimes_plot(odimyv, timeunit, timekind, timefmt)
1128	else:
1129	print errormsg
1130	print ' ' + fname + ": time variable '" + timename + "' not found!!"
1131	quit(-1)
1132
1133	shading_nx = np.zeros((2), dtype=np.float)
1134	shading_nx[0] = np.float(shadminmax.split(',')[0])
1135	shading_nx[1] = np.float(shadminmax.split(',')[1])
1136
1137	closeval = drw.Str_Bool(close)
1138
1139	drw.plot_2D_shadow_time(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu, \
1140	dimnamesv, colbarn, shading_nx, varunits, figtitle, figkind, revals, timeaxis, \
1141	timepos, timelabels, False)
1142
1143	# Line values
1144	##
1145	ncfilel = ncfile.split(',')[1]
1146
1147	vnamelfig = values.split(':')[0].split(',')[1]
1148	varnl = varn.split(',')[1]
1149
1150	objlf = NetCDFFile(ncfilel,'r')
1151	objlvar = objlf.variables[varnl]
1152
1153	linevals = objlvar[:]
1154	if reva0 == 'tranpose':
1155	plt.plot (linevals, odimxv, '-', color='k')
1156	else:
1157	plt.plot (odimxv, linevals, '-', color='k')
1158
1159	objsf.close()
1160	objsl.close()
1161
1162	return
1163
1164	def draw_topo_geogrid(ncfile, values):
1165	""" plotting geo_em.d[nn].nc topography from WPS files
1166	draw_topo_geogrid(ncfile, values)
1167	ncfile= geo_em.d[nn].nc file to use
1168	values= [minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]
1169	[min/max]Topo: minimum and maximum values of topography to draw
1170	lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] or None
1171	title: title of the graph
1172	graphic_kind: kind of figure (jpg, pdf, png)
1173	mapvalues: map characteristics [proj],[res]
1174	see full documentation: http://matplotlib.org/basemap/
1175	[proj]: projection
1176	* 'cyl', cilindric
1177	* 'lcc', lambert conformal
1178	[res]: resolution:
1179	* 'c', crude
1180	* 'l', low
1181	* 'i', intermediate
1182	* 'h', high
1183	* 'f', full
1184	"""
1185	fname = 'draw_topo_geogrid'
1186
1187	if values == 'h':
1188	print fname + '_____________________________________________________________'
1189	print draw_topo_geogrid.__doc__
1190	quit()
1191
1192	expectargs = ['[minTopo]','[maxTopo]', '[lonlatL]', '[title]', '[graphic_kind]', \
1193	'[mapvalues]']
1194
1195	drw.check_arguments(fname,5,values,':',expectargs)
1196
1197	mintopo = values.split(':')[0].split(',')[0]
1198	maxtopo = values.split(':')[0].split(',')[1]
1199
1200	lonlatLS = values.split(':')[1]
1201	lonlatLv = lonlatLS.split(',')[0]
1202
1203	if lonlatLv == 'None':
1204	lonlatL = None
1205	else:
1206	lonlatL = np.zeros((4), dtype=np.float)
1207	lonlatL[0] = np.float(lonlatLS.split(',')[0])
1208	lonlatL[1] = np.float(lonlatLS.split(',')[1])
1209	lonlatL[2] = np.float(lonlatLS.split(',')[2])
1210	lonlatL[3] = np.float(lonlatLS.split(',')[3])
1211
1212	grtit = values.split(':')[2]
1213	kindfig = values.split(':')[3]
1214	mapvalues = values.split(':')[4]
1215
1216	if not os.path.isfile(ncfile):
1217	print errormsg
1218	print ' ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1219	quit(-1)
1220
1221	objdomf = NetCDFFile(ncfile, 'r')
1222
1223	objhgt = objdomf.variables['HGT_M']
1224	objlon = objdomf.variables['XLONG_M']
1225	objlat = objdomf.variables['XLAT_M']
1226
1227	topography = objhgt[0,:,:]
1228
1229	drw.plot_topo_geogrid(topography, objlon, objlat, mintopo, maxtopo, lonlatL, \
1230	grtit, kindfig, mapvalues, True)
1231
1232	objdomf.close()
1233
1234	return
1235
1236	def draw_topo_geogrid_boxes(ncfiles, values):
1237	""" plotting different geo_em.d[nn].nc topography from WPS files
1238	draw_topo_geogrid_boxes(ncfile, values)
1239	ncfiles= ',' list of geo_em.d[nn].nc files to use (fisrt as topographyc reference)
1240	values= [minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]:[labels]
1241	[min/max]Topo: minimum and maximum values of topography to draw
1242	lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] or None
1243	title: title of the graph
1244	graphic_kind: kind of figure (jpg, pdf, png)
1245	mapvalues: map characteristics [proj],[res]
1246	see full documentation: http://matplotlib.org/basemap/
1247	[proj]: projection
1248	* 'cyl', cilindric
1249	* 'lcc', lambert conformal
1250	[res]: resolution:
1251	* 'c', crude
1252	* 'l', low
1253	* 'i', intermediate
1254	* 'h', high
1255	* 'f', full
1256	labels= labels to write in the graph
1257	"""
1258	# import matplotlib as mpl
1259	# mpl.use('Agg')
1260	import matplotlib.pyplot as plt
1261
1262	fname = 'draw_topo_geogrid_boxes'
1263
1264	if values == 'h':
1265	print fname + '_____________________________________________________________'
1266	print draw_topo_geogrid_boxes.__doc__
1267	quit()
1268
1269	mintopo = values.split(':')[0].split(',')[0]
1270	maxtopo = values.split(':')[0].split(',')[1]
1271
1272	lonlatLS = values.split(':')[1]
1273	lonlatLv = lonlatLS.split(',')[0]
1274
1275	if lonlatLv == 'None':
1276	lonlatL = None
1277	else:
1278	lonlatL = np.zeros((4), dtype=np.float)
1279	lonlatL[0] = np.float(lonlatLS.split(',')[0])
1280	lonlatL[1] = np.float(lonlatLS.split(',')[1])
1281	lonlatL[2] = np.float(lonlatLS.split(',')[2])
1282	lonlatL[3] = np.float(lonlatLS.split(',')[3])
1283
1284	grtit = values.split(':')[2]
1285	kindfig = values.split(':')[3]
1286	mapvalues = values.split(':')[4]
1287	labels = values.split(':')[5]
1288
1289	ncfile = ncfiles.split(',')[0]
1290	if not os.path.isfile(ncfile):
1291	print errormsg
1292	print ' ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1293	quit(-1)
1294
1295	objdomf = NetCDFFile(ncfile, 'r')
1296
1297	objhgt = objdomf.variables['HGT_M']
1298	objlon0 = objdomf.variables['XLONG_M']
1299	objlat0 = objdomf.variables['XLAT_M']
1300
1301	topography = objhgt[0,:,:]
1302
1303	Nfiles = len(ncfiles.split(','))
1304	boxlabels = labels.split(',')
1305
1306	Xboxlines = []
1307	Yboxlines = []
1308
1309	for ifile in range(Nfiles):
1310	ncfile = ncfiles.split(',')[ifile]
1311	# print ifile, ncfile
1312	if not os.path.isfile(ncfile):
1313	print errormsg
1314	print ' ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1315	quit(-1)
1316
1317	objdomfi = NetCDFFile(ncfile, 'r')
1318
1319	objlon = objdomfi.variables['XLONG_M']
1320	objlat = objdomfi.variables['XLAT_M']
1321
1322	dx = objlon.shape[2]
1323	dy = objlon.shape[1]
1324
1325	Xboxlines.append(objlon[0,0,:])
1326	Yboxlines.append(objlat[0,0,:])
1327	Xboxlines.append(objlon[0,dy-1,:])
1328	Yboxlines.append(objlat[0,dy-1,:])
1329	Xboxlines.append(objlon[0,:,0])
1330	Yboxlines.append(objlat[0,:,0])
1331	Xboxlines.append(objlon[0,:,dx-1])
1332	Yboxlines.append(objlat[0,:,dx-1])
1333
1334	objdomfi.close()
1335
1336	drw.plot_topo_geogrid_boxes(topography, Xboxlines, Yboxlines, boxlabels, \
1337	objlon0, objlat0, mintopo, maxtopo, lonlatL, grtit, kindfig, mapvalues, True)
1338
1339	objdomf.close()
1340
1341	return
1342
1343	def movievalslice(origslice, dimmovien, framenum):
1344	""" Function to provide variable slice according to a geneation of a movie
1345	movievals(origslice, dimmovien, framenum)
1346	[origslice]= slice original as [dimname1]\|[val1],[...,[dimnameN]\|[valN]]
1347	([val] = -1, full length)
1348	[dimmovien]= name of the dimension to produce the movie
1349	[framenum]= value of the frame to substitue in [origslice] as
1350	[dimmovien]\|[framenum]
1351	>>> movievalslice('East_West\|-1,North_South\|-1,Time\|2','Time',0)
1352	East_West\|-1,North_South\|-1,Time\|0
1353	"""
1354
1355	fname = 'movievalslice'
1356
1357	if origslice == 'h':
1358	print fname + '_____________________________________________________________'
1359	print movievalslice.__doc__
1360	quit()
1361
1362	dims = origslice.split(',')
1363
1364	movieslice = ''
1365	idim = 0
1366
1367	for dimn in dims:
1368	dn = dimn.split('\|')[0]
1369	if dn == dimmovien:
1370	movieslice = movieslice + dn + '\|' + str(framenum)
1371	else:
1372	movieslice = movieslice + dimn
1373	if idim < len(dims)-1: movieslice = movieslice + ','
1374
1375	idim = idim + 1
1376
1377	return movieslice
1378
1379	class Capturing(list):
1380	""" Class to capture function output as a list
1381	from: http://stackoverflow.com/questions/16571150/how-to-capture-stdout-output-from-a-python-function-call
1382	"""
1383	# from cStringIO import StringIO
1384
1385	def __enter__(self):
1386	self._stdout = sys.stdout
1387	sys.stdout = self._stringio = StringIO()
1388	return self
1389	def __exit__(self, *args):
1390	self.extend(self._stringio.getvalue().splitlines())
1391	sys.stdout = self._stdout
1392
1393	def create_movie(netcdfile, values, variable):
1394	""" Function to create a movie assuming ImageMagick installed!
1395	values= [graph]#[movie_dimension]#[graph_values]
1396	[graph]: which graphic
1397	[movie_dimension]: [dimnmovie]@[dimvmovie]@[moviedelay]@[interval]
1398	[dimnmovie]; name of the dimension from which make the movie
1399	[dimvmovie]; name of the variable with the values of the dimension
1400	[moviedelay]; delay between frames
1401	[interval]; [beg]@[end]@[freq] or -1 (all)
1402	[graph_values]: values to generate the graphic
1403	netcdfile= netCDF file
1404	variable= variable to use (when applicable)
1405	"""
1406	fname = 'create_movie'
1407
1408	if values == 'h':
1409	print fname + '_____________________________________________________________'
1410	print create_movie.__doc__
1411	quit()
1412
1413	graph = values.split('#')[0]
1414	movie_dim = values.split('#')[1]
1415	graph_vals = values.split('#')[2]
1416
1417	ncobj = NetCDFFile(netcdfile, 'r')
1418
1419	# Movie dimension
1420	##
1421	dimnmovie = movie_dim.split('@')[0]
1422	dimvmovie = movie_dim.split('@')[1]
1423	moviedelay = movie_dim.split('@')[2]
1424	moviebeg = int(movie_dim.split('@')[3])
1425
1426	if not drw.searchInlist(ncobj.dimensions.keys(),dimnmovie):
1427	print errormsg
1428	print ' ' + fname + ": file '" + netcdfile + "' has not dimension '" + \
1429	dimnmovie + "' !!!"
1430	quit(-1)
1431
1432	objdmovie = ncobj.dimensions[dimnmovie]
1433	dmovie = len(objdmovie)
1434	if moviebeg != -1:
1435	moviend = int(movie_dim.split('@')[4])
1436	moviefreq = int(movie_dim.split('@')[5])
1437	else:
1438	moviebeg = 0
1439	moviend = dmovie
1440	moviefreq = 1
1441
1442	if dimvmovie == 'WRFTimes':
1443	objvdmovie = ncobj.variables['Times']
1444	vdmovieunits = ''
1445	valsdmovie = []
1446	for it in range(objvdmovie.shape[0]):
1447	valsdmovie.append(drw.datetimeStr_conversion(objvdmovie[it,:], \
1448	'WRFdatetime', 'Y/m/d H-M-S'))
1449	elif dimvmovie == 'CFtime':
1450	objvdmovie = ncobj.variables['time']
1451	vdmovieunits = ''
1452	print objvdmovie.units
1453	valsdmovie0 = drw.netCDFdatetime_realdatetime(objvdmovie.units, 'standard', \
1454	objvdmovie[:])
1455	valsdmovie = []
1456	for it in range(objvdmovie.shape[0]):
1457	valsdmovie.append(drw.datetimeStr_conversion(valsdmovie0[it,:], \
1458	'matYmdHMS', 'Y/m/d H-M-S'))
1459	else:
1460	if not drw.searchInlist(ncobj.variables.keys(),dimvmovie):
1461	print errormsg
1462	print ' ' + fname + ": file '" + netcdfile + "' has not variable '" + \
1463	dimvmovie + "' !!!"
1464	quit(-1)
1465	vdmovieunits = objvdmovie.getncattr('units')
1466	objvdmovie = ncobj.variables[dimvmovie]
1467	if len(objvdmovie.shape) == 1:
1468	vasldmovie = objvdmovie[:]
1469	else:
1470	print errormsg
1471	print ' ' + fname + ': shape', objvdmovie.shape, 'of variable with ' + \
1472	'dimension movie values not ready!!!'
1473	quit(-1)
1474
1475	ncobj.close()
1476	os.system('rm frame_*.png > /dev/null')
1477
1478	# graphic
1479	##
1480	if graph == 'draw_2D_shad':
1481	graphvals = graph_vals.split(':')
1482
1483	for iframe in range(moviebeg,moviend,moviefreq):
1484	iframeS = str(iframe).zfill(4)
1485
1486	drw.percendone((iframe-moviebeg)/moviefreq,(moviend-moviebeg)/moviefreq, \
1487	5, 'frames')
1488	titgraph = dimnmovie + '\|=\|' + str(valsdmovie[iframe]) + '\|' + \
1489	vdmovieunits
1490
1491	graphvals[1] = movievalslice(graphvals[1],dimnmovie,iframe)
1492	graphvals[6] = titgraph
1493	graphvals[7] = 'png'
1494
1495	graphv = drw.numVector_String(graphvals, ":")
1496
1497	with Capturing() as output:
1498	draw_2D_shad(netcdfile, graphv, variable)
1499
1500	os.system('mv 2Dfields_shadow.png frame_' + iframeS + '.png')
1501	else:
1502	print errormsg
1503	print ' ' + fname + ": graphic '" + graph + "' not defined !!!"
1504	quit(-1)
1505
1506	os.system('convert -delay ' + moviedelay + ' -loop 0 frame_*.png create_movie.gif')
1507
1508	print "Succesfuly creation of movie file 'create_movie.gif' !!!"
1509
1510	return
1511
1512	def draw_lines(ncfilens, values, varname):
1513	""" Function to draw different lines at the same time from different files
1514	draw_lines(ncfilens, values, varname):
1515	ncfilens= [filen] ',' separated list of netCDF files
1516	values= [dimvname]:[valuesaxis]:[dimtit]:[leglabels]:[vtit]:[title]:[locleg]:[graphk]
1517	[dimvname]: ',' list of names of the variable with he values of the common dimension
1518	[valuesaxis]: which axis will be used for the values ('x', or 'y')
1519	[dimtit]: title for the common dimension
1520	[leglabels]: ',' separated list of names for the legend
1521	[vartit]: name of the variable in the graph
1522	[title]: title of the plot ('\|' for spaces)
1523	[locleg]: location of the legend (-1, autmoatic)
1524	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
1525	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
1526	9: 'upper center', 10: 'center'
1527	[graphk]: kind of the graphic
1528	varname= variable to plot
1529	values= 'XLAT:x:latitude:32x32:$wss^{*}$:wss Taylor's turbulence term:pdf'
1530	"""
1531
1532	fname = 'draw_lines'
1533
1534	if values == 'h':
1535	print fname + '_____________________________________________________________'
1536	print draw_lines.__doc__
1537	quit()
1538
1539	expectargs = '[dimvname]:[valuesaxis]:[dimtit]:[leglabels]:[vtit]:[title]:[locleg]:[graphk]'
1540	drw.check_arguments(fname,len(expectargs.split(':')),values,':',expectargs)
1541
1542	ncfiles = ncfilens.split(',')
1543	dimvnames = values.split(':')[0]
1544	valuesaxis = values.split(':')[1]
1545	dimtit = values.split(':')[2]
1546	leglabels = values.split(':')[3].replace('_','\_')
1547	vartit = values.split(':')[4]
1548	title = values.split(':')[5].replace('\|',' ')
1549	graphk = values.split(':')[6]
1550
1551	Nfiles = len(ncfiles)
1552
1553	# Getting trajectotries
1554	##
1555
1556	varvalues = []
1557	dimvalues = []
1558
1559	print ' ' + fname
1560	ifn = 0
1561	for ifile in ncfiles:
1562	filen = ifile.split('@')[0]
1563
1564	print ' filen:',filen
1565
1566	if not os.path.isfile(filen):
1567	print errormsg
1568	print ' ' + fname + ": netCDF file '" + filen + "' does not exist !!"
1569	quit(-1)
1570
1571	objfile = NetCDFFile(filen, 'r')
1572
1573	if dimvnames.find(',') != -1:
1574	dimvname = dimvnames.split(',')
1575	else:
1576	dimvname = [dimvnames]
1577
1578	found = False
1579	for dvn in dimvname:
1580	if objfile.variables.has_key(dvn):
1581	found = True
1582	break
1583
1584	if not found:
1585	print errormsg
1586	print ' ' + fname + ": netCDF file '" + filen + \
1587	"' does not have variables '" + dimvname + "' !!"
1588	quit(-1)
1589
1590	if not objfile.variables.has_key(varname):
1591	print errormsg
1592	print ' ' + fname + ": netCDF file '" + filen + \
1593	"' does not have variable '" + varname + "' !!"
1594	quit(-1)
1595
1596	vvobj = objfile.variables[varname]
1597	if len(vvobj.shape) != 1:
1598	print errormsg
1599	print ' ' + fname + ': wrong shape:',vvobj.shape," of variable '" + \
1600	varname + "' !!"
1601	quit(-1)
1602
1603	for dimvn in dimvname:
1604	if searchInlist(objfile.variables, dimvn):
1605	vdobj = objfile.variables[dimvname]
1606	if len(vdobj.shape) != 1:
1607	print errormsg
1608	print ' ' + fname + ': wrong shape:',vdobj.shape, \
1609	" of variable '" + dimvn + "' !!"
1610	quit(-1)
1611	break
1612
1613	varvalues.append(vvobj[:])
1614	dimvalues.append(vdobj[:])
1615
1616	if ifn == 0:
1617	varunits = vvobj.units
1618
1619	objfile.close()
1620
1621	ifn = ifn + 1
1622
1623
1624	drw.plot_lines(dimvalues, varvalues, valuesaxis, dimtit, leglabels.split(','), \
1625	vartit, varunits, title, graphk)
1626
1627	return
1628
1629	def draw_lines_time(ncfilens, values, varname):
1630	""" Function to draw different lines at the same time from different files with times
1631	draw_lines_time(ncfilens, values, varname):
1632	ncfilens= [filen] ',' separated list of netCDF files
1633	values= [dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];[timevals];[locleg];[graphk]
1634	[dimvname]: name of the variable with he values of the common dimension
1635	[valuesaxis]: which axis will be used for the values ('x', or 'y')
1636	[dimtit]: title for the common dimension
1637	[leglabels]: ',' separated list of names for the legend
1638	[vartit]: name of the variable in the graph
1639	[title]: title of the plot ('\|' for spaces)
1640	[timevals]: [timen]\|[units]\|[kind]\|[tfmt] time labels characteristics
1641	[timen]; name of the time variable
1642	[units]; units string according to CF conventions ([tunits] since
1643	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
1644	[kind]; kind of output
1645	'Nval': according to a given number of values as 'Nval',[Nval]
1646	'exct': according to an exact time unit as 'exct',[tunit];
1647	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1648	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1649	'l': milisecond
1650	[tfmt]; desired format
1651	[locleg]: location of the legend (-1, autmoatic)
1652	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
1653	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
1654	9: 'upper center', 10: 'center'
1655	[graphk]: kind of the graphic
1656	varname= variable to plot
1657	values= 'time;y;time ([DD]${[HH]}$);32x32;$wss^{*}$;wss Taylor's turbulence term;time\|hours!since!1949-12-01_00:00:00;exct,12,h\|%d$^{%H}$;2;pdf'
1658	"""
1659
1660	fname = 'draw_lines_time'
1661
1662	if values == 'h':
1663	print fname + '_____________________________________________________________'
1664	print draw_lines_time.__doc__
1665	quit()
1666
1667	expectargs = '[dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];'
1668	expectargs = expectargs + '[timevals];[locleg];[graphk]'
1669	drw.check_arguments(fname,len(expectargs.split(';')),values,';',expectargs)
1670
1671	ncfiles = ncfilens.split(',')
1672	dimvname = values.split(';')[0]
1673	valuesaxis = values.split(';')[1]
1674	dimtit = values.split(';')[2]
1675	leglabels = values.split(';')[3].replace('_','\_')
1676	vartit = values.split(';')[4]
1677	title = values.split(';')[5].replace('\|',' ')
1678	timevals = values.split(';')[6]
1679	locleg = int(values.split(';')[7])
1680	graphk = values.split(';')[8]
1681
1682	Nfiles = len(ncfiles)
1683
1684	# Getting values
1685	##
1686	varvalues = []
1687	dimvalues = []
1688	timvalues = []
1689	timvals0 = timvalues
1690
1691	print ' ' + fname
1692	ifn = 0
1693	mintval = 1.e20
1694	maxtval = -1.e20
1695
1696	for ifile in ncfiles:
1697	filen = ifile.split('@')[0]
1698
1699	print ' filen:',filen
1700
1701	if not os.path.isfile(filen):
1702	print errormsg
1703	print ' ' + fname + ": netCDF file '" + filen + "' does not exist !!"
1704	quit(-1)
1705
1706	objfile = NetCDFFile(filen, 'r')
1707
1708	if not objfile.variables.has_key(dimvname):
1709	print errormsg
1710	print ' ' + fname + ": netCDF file '" + filen + \
1711	"' does not have variable '" + dimvname + "' !!"
1712	quit(-1)
1713
1714	if not objfile.variables.has_key(varname):
1715	print errormsg
1716	print ' ' + fname + ": netCDF file '" + filen + \
1717	"' does not have variable '" + varname + "' !!"
1718	quit(-1)
1719
1720	vvobj = objfile.variables[varname]
1721	if len(vvobj.shape) != 1:
1722	print errormsg
1723	print ' ' + fname + ': wrong shape:',vvobj.shape," of variable '" + \
1724	varname + "' !!"
1725	quit(-1)
1726
1727	vdobj = objfile.variables[dimvname]
1728	if len(vdobj.shape) != 1:
1729	print errormsg
1730	print ' ' + fname + ': wrong shape:',vdobj.shape," of variable '" + \
1731	dimvname + "' !!"
1732	quit(-1)
1733
1734	varvalues.append(vvobj[:])
1735	dimvalues.append(vdobj[:])
1736
1737	mindvals = np.min(vdobj[:])
1738	maxdvals = np.max(vdobj[:])
1739
1740	if mindvals < mintval: mintval = mindvals
1741	if maxdvals > maxtval: maxtval = maxdvals
1742
1743	if ifn == 0:
1744	varunits = drw.units_lunits(vvobj.units)
1745
1746	objfile.close()
1747
1748	ifn = ifn + 1
1749
1750	# Times
1751	timename = timevals.split('\|')[0]
1752	timeunit = timevals.split('\|')[1].replace('!',' ')
1753	timekind = timevals.split('\|')[2]
1754	timefmt = timevals.split('\|')[3]
1755
1756	dtvals = (maxtval - mintval)/5.
1757	tvals = np.arange(mintval, maxtval, dtvals/2.)
1758
1759	timepos, timelabels = drw.CFtimes_plot(tvals, timeunit, timekind, timefmt)
1760
1761	drw.plot_lines_time(dimvalues, varvalues, valuesaxis, dimtit, leglabels.split(','), \
1762	vartit, varunits, timepos, timelabels, title, locleg, graphk)
1763
1764	return
1765
1766	def draw_Neighbourghood_evol(filen, values, variable):
1767	""" Function to draw the temporal evolution of a neighbourghood around a point
1768	draw_Neighbourghood_evol(filen, values, variable)
1769	filen= netCDF file name
1770	values= [gvarname]:[dimsval]:[neigdims]:[Nneig]:[Ncol]:[timetits]:[tkinds]:
1771	[timefmts]:[gtitle]:[shadxtrms]:[cbar]:[gkind]:[ofile]
1772	[dimsval]: [dimn1]\|[val1]\|[dimv1],...,[dimnN]\|[valN]\|[dimvN] dimension names, values to get
1773	(-1, for all; no name/value pair given full length) and variable with values of the dimension
1774	NOTE: when dimsval[X,Y] == neigdims[X,Y], valX,valY --> valX,valY-Nneig/2, valX,valY+Nneig/2
1775	[neigdims]: [dimnX],[dimnY] dimensions mnames along which the neigbourghood should be defined
1776	[Nneig]: Number of grid points of the full side of the box (odd value)
1777	[Ncol]: Number of columns ('auto': square final plot)
1778	[gvarname]: name of the variable to appear in the graph
1779	[timetits]: [titX],[titY] titles of the axes ('\|' for spaces)
1780	[tkinds]: [tkindX]\|[tkindY] kinds of time to appear in the graph
1781	'Nval': according to a given number of values as 'Nval',[Nval]
1782	'exct': according to an exact time unit as 'exct',[tunit];
1783	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1784	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1785	'l': milisecond
1786	[timefmts]: [tfmtX],[tfmtY] format of the time labels
1787	[gtitle]: title of the graphic ('\|' for spaces)
1788	[shadxtrms]: Extremes for the shading
1789	[cbar]: colorbar to use
1790	[gkind]: kind of graphical output
1791	[ofile]: True/False whether the netcdf with data should be created or not
1792	variable= name of the variable
1793	values = 'q:Time\|-1\|Times,bottom_top\|6\|ZNU,south_north\|3\|XLAT,west_east\|26\|XLONG:south_north,west_east:5:auto:time\|($[DD]^{[HH]}$),time\|($[DD]^{[HH]}$):exct,2,h\|exct,1,d:$%d^{%H}$,$%d^{%H}$:5\|pts\|neighbourghood\|temporal\|evolution:0.0,0.004:BuPu:pdf:True'
1794	"""
1795
1796	fname = 'draw_Neighbourghood_evol'
1797
1798	if values == 'h':
1799	print fname + '_____________________________________________________________'
1800	print draw_Neighbourghood_evol.__doc__
1801	quit()
1802
1803	expectargs = ['[gvarname]', '[dimsval]', '[neigdims]', '[Nneig]', '[Ncol]', \
1804	'[timetits]', '[tkinds]', '[timefmts]', '[gtitle]', '[shadxtrms]', '[cbar]', \
1805	'[gkind]', '[ofile]']
1806
1807	drw.check_arguments(fname,len(expectargs),values,':',expectargs)
1808
1809	gvarname = values.split(':')[0]
1810	dimsval = values.split(':')[1].split(',')
1811	neigdims = values.split(':')[2].split(',')
1812	Nneig = int(values.split(':')[3])
1813	Ncol0 = values.split(':')[4]
1814	timetits = values.split(':')[5].split(',')
1815	timekinds = values.split(':')[6].split('\|')
1816	timefmts = values.split(':')[7].split(',')
1817	gtitle = values.split(':')[8].replace('\|',' ')
1818	shadxtrms = values.split(':')[9].split(',')
1819	cbar = values.split(':')[10]
1820	gkind = values.split(':')[11]
1821	ofile = values.split(':')[12]
1822
1823	if Ncol0 != 'auto':
1824	Ncol = int(Ncol0)
1825	else:
1826	Ncol = Ncol0
1827
1828	timetits[0] = timetits[0].replace('\|',' ')
1829	timetits[1] = timetits[1].replace('\|',' ')
1830
1831	if np.mod(Nneig,2) == 0:
1832	print errormsg
1833	print ' ' + fname + ": an odd value for 'Nneig':",Nneig,'is required !!!'
1834	quit(-1)
1835
1836	Nneig2 = int(Nneig/2)
1837
1838	# Values to slice the variable
1839	dimvslice = {}
1840	dimvvalues = {}
1841	for dimvs in dimsval:
1842	dimn = dimvs.split('\|')[0]
1843	dimv = int(dimvs.split('\|')[1])
1844	dimnv = dimvs.split('\|')[2]
1845
1846	dimvvalues[dimn] = dimnv
1847	dimvslice[dimn] = dimv
1848
1849	ncobj = NetCDFFile(filen, 'r')
1850
1851	varobj = ncobj.variables[variable]
1852
1853	slicevar = []
1854	newdimn = []
1855	newdimsvar = {}
1856
1857	for dimn in varobj.dimensions:
1858	if not drw.searchInlist(dimvslice.keys(), dimn):
1859	dimsize = len(ncobj.dimensions[dimn])
1860	slicevar.append(slice(0, dimsize+1))
1861	newdimn.append(dimn)
1862	newdimsvar[dimn] = dimseize
1863
1864	for dimslicen in dimvslice.keys():
1865	if dimn == dimslicen:
1866	if dimvslice[dimn] != -1:
1867	if drw.searchInlist(neigdims, dimn):
1868	slicevar.append(slice(dimvslice[dimn]-Nneig2, \
1869	dimvslice[dimn]+Nneig2+1))
1870	newdimn.append(dimn)
1871	newdimsvar[dimn] = Nneig
1872	break
1873	else:
1874	slicevar.append(slice(dimvslice[dimn], dimvslice[dimn]+1))
1875	break
1876	else:
1877	dimsize = len(ncobj.dimensions[dimn])
1878	slicevar.append(slice(0, dimsize+1))
1879	newdimn.append(dimn)
1880	newdimsvar[dimn] = dimsize
1881	break
1882
1883	varv = varobj[tuple(slicevar)]
1884
1885	if len(newdimn) != 3:
1886	print errormsg
1887	print ' ' + fname + ': sliced variable with shape=', varv.shape, \
1888	' must have three dimensions',len(varv.shape),'given !!'
1889	quit(-1)
1890
1891	newdims = []
1892	for nwdims in newdimn:
1893	newdims.append(newdimsvar[nwdims])
1894
1895	# The dimension which is not in the neighbourhood dimensions must be time!
1896	for dim1 in newdimn:
1897	if not drw.searchInlist(neigdims, dim1):
1898	dimt = newdimsvar[dim1]
1899	dimtime = dim1
1900
1901	if Ncol == 'auto':
1902	dimtsqx = int(np.sqrt(dimt)) + 1
1903	dimtsqy = int(np.sqrt(dimt)) + 1
1904	else:
1905	dimtsqx = int(Ncol)
1906	dimtsqy = dimt/dimtsqx + 1
1907
1908	neighbourghood = np.ones((dimtsqyNneig,dimtsqxNneig), dtype=np.float)*fillValue
1909
1910	for it in range(dimt):
1911	ity = int(it/dimtsqx)
1912	itx = it-ity*dimtsqx
1913
1914	itty = (dimtsqy - ity - 1)*Nneig + Nneig2
1915	ittx = itx*Nneig + Nneig2
1916
1917	neighbourghood[itty-Nneig2:itty+Nneig2+1,ittx-Nneig2:ittx+Nneig2+1]= \
1918	varv[it,::-1,:]
1919
1920	variablevals = drw.variables_values(variable)
1921	if drw.searchInlist(varobj.ncattrs(), 'units'):
1922	vunits = varobj.units
1923	else:
1924	vunits = variablevals[5]
1925
1926	# Time values at the X/Y axes
1927	if ncobj.variables[dimvvalues[dimtime]].dtype == '\|S1':
1928	print ' ' + fname + ': WRF time variable!'
1929	refdate = '19491201000000'
1930	tunitsval = 'hours'
1931	dimtvalues = np.zeros((dimt), dtype=np.float)
1932	tvals = ncobj.variables[dimvvalues[dimtime]]
1933	yrref=refdate[0:4]
1934	monref=refdate[4:6]
1935	dayref=refdate[6:8]
1936	horref=refdate[8:10]
1937	minref=refdate[10:12]
1938	secref=refdate[12:14]
1939
1940	refdateS = yrref + '/' + monref + '/' + dayref + '_' + horref + ':' + \
1941	minref + ':' + secref
1942	tunits = tunitsval + ' since ' + refdateS
1943	for it in range(dimt):
1944	wrfdates = drw.datetimeStr_conversion(tvals[it,:],'WRFdatetime', 'matYmdHMS')
1945	dimtvalues[it] = drw.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
1946	else:
1947	dimtvalues = ncobj.variables[dimvvalues[dimtime]][:]
1948	tunits = ncobj.variables[newdimsvar[dimtime]].units
1949
1950	dimxv = dimtvalues[0:dimtsqx]
1951	dimyv = dimtvalues[0:dimt:dimtsqx]
1952
1953	dimn = ['time','time']
1954
1955	if ofile == 'True':
1956	ofilen = 'Neighbourghood_evol.nc'
1957	newnc = NetCDFFile(ofilen, 'w')
1958	# Dimensions
1959	newdim = newnc.createDimension('time',None)
1960	newdim = newnc.createDimension('y',dimtsqy*Nneig)
1961	newdim = newnc.createDimension('x',dimtsqx*Nneig)
1962	# Dimension values
1963	newvar = newnc.createVariable('time','f8',('time'))
1964	newvar[:] = np.arange(dimt)
1965	newattr = drw.basicvardef(newvar, 'time','time',tunits)
1966	# Neighbourhghood variable
1967	newvar = newnc.createVariable(variable + 'neigevol', 'f4', ('y','x'), \
1968	fill_value=fillValue)
1969	newvar[:] = neighbourghood
1970
1971	newnc.sync()
1972	newnc.close()
1973	print fname + ": Successfull generation of file '" + ofilen + "' !!"
1974
1975	# Time ticks
1976	timeposX, timelabelsX = drw.CFtimes_plot(dimxv, tunits, timekinds[0], timefmts[0])
1977	timeposY, timelabelsY = drw.CFtimes_plot(dimyv, tunits, timekinds[1], timefmts[1])
1978
1979	timepos = [timeposX[0:len(timeposX)], timeposY[len(timeposY):0:-1]]
1980	timelabels = [timelabelsX[0:len(timeposX)], timelabelsY[0:len(timeposY)]]
1981
1982	for i in range(2):
1983	if shadxtrms[i][0:1] != 'S':
1984	shadxtrms[i] = np.float(shadxtrms[i])
1985
1986	drw.plot_Neighbourghood_evol(neighbourghood, dimxv, dimyv, gvarname, timetits, \
1987	timepos, timelabels, cbar, Nneig, shadxtrms, vunits, gtitle, gkind, True)
1988
1989	def draw_timeSeries(filen, values, variables):
1990	""" Function to draw a time-series
1991	draw_timeSeries(filen, values, variable):
1992	filen= name of the file
1993	values= [gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[locleg]:[gkind]
1994	[gvarname]: name of the variable to appear in the graph
1995	[timetit]: title of the time axis (assumed x-axis, '\|' for spaces)
1996	[tkind]: kind of time to appear in the graph (assumed x-axis)
1997	'Nval': according to a given number of values as 'Nval',[Nval]
1998	'exct': according to an exact time unit as 'exct',[tunit];
1999	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2000	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2001	'l': milisecond
2002	[timefmt]: format of the time labels
2003	[title]: title of the graphic ('\|' for spaces)
2004	[locleg]: location of the legend (-1, autmoatic)
2005	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2006	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2007	9: 'upper center', 10: 'center'
2008	[gkind]: kind of graphical output
2009	variables= [varname],[timename] names of variable and variable with times
2010	draw_timeSeries('wrfout_d01_1979-12-01_00:00:00_bottom_top_B6-E6-I1_south_north_B3-E3-I1_west_east_B26-E26-I1.nc', 'dt_con:time\|($[DD]^{[HH]}$):exct,12,h:$%d^{%H}$:time\|evolution\|at\|-1\|6\|3\|26:1:pdf', 'LDQCON,time')
2011	"""
2012
2013	fname = 'draw_timeSeries'
2014
2015	if values == 'h':
2016	print fname + '_____________________________________________________________'
2017	print draw_timeSeries.__doc__
2018	quit()
2019
2020	expectargs = ['[gvarname]', '[timetit]', '[tkind]', '[timefmt]', '[title]', \
2021	'[locleg]', '[gkind]']
2022
2023	drw.check_arguments(fname,len(expectargs),values,':',expectargs)
2024
2025	gvarname = values.split(':')[0]
2026	timetit = values.split(':')[1].replace('\|',' ')
2027	tkind = values.split(':')[2]
2028	timefmt = values.split(':')[3]
2029	title = values.split(':')[4].replace('\|',' ')
2030	locleg = int(values.split(':')[5])
2031	gkind = values.split(':')[6]
2032
2033	ncobj = NetCDFFile(filen, 'r')
2034
2035	variable = variables.split(',')[0]
2036	timevar = variables.split(',')[1]
2037
2038	if not ncobj.variables.has_key(variable):
2039	print errormsg
2040	print ' ' + fname + ": file '" + filen + "' does not have variable '" + \
2041	variable + "' !!"
2042	quit(-1)
2043
2044	if not ncobj.variables.has_key(timevar):
2045	print errormsg
2046	print ' ' + fname + ": file '" + filen + "' does not have variable time '" \
2047	+ timevar + "' !!"
2048	quit(-1)
2049
2050	varobj = ncobj.variables[variable]
2051	timeobj = ncobj.variables[timevar]
2052
2053	dimt = len(timeobj[:])
2054	varvals = np.zeros((2,dimt), dtype=np.float)
2055
2056	gunits = varobj.getncattr('units')
2057	tunits = timeobj.getncattr('units')
2058
2059	varvals[0,:], valpot, newgunits, Spot = drw.pot_values(varobj[:].flatten(), gunits)
2060	varvals[1,:] = timeobj[:]
2061
2062	tseriesvals = []
2063	tseriesvals.append(varvals)
2064
2065	drw.plot_TimeSeries(tseriesvals, Spot + drw.units_lunits(gunits), tunits, \
2066	'TimeSeries', gvarname, timetit, tkind, timefmt, title, \
2067	gvarname.replace('_','\_'), locleg, gkind)
2068
2069	return
2070
2071	#draw_timeSeries('wrfout_d01_1979-12-01_00:00:00_bottom_top_B6-E6-I1_south_north_B3-E3-I1_west_east_B26-E26-I1.nc', 'dt_con:time\|($[DD]^{[HH]}$):exct,12,h:$%d^{%H}$:time\|evolution\|at\|-1\|6\|3\|26:1:pdf', 'LDQCON,time')
2072
2073	def draw_trajectories(trjfilens, values, observations):
2074	""" Function to draw different trajectories at the same time
2075	draw_trajectories(trjfilens, values, observations):
2076	trjfilens= [filen]@[Tint]@[map] ',' separated list of files with trajectories,
2077	time intervals and reference maps (first one will be used to plot)
2078	[filen]: name of the file to use (lines with '#', not readed) as:
2079	[t-step] [x] [y]
2080	[Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2081	[map]: [file]#[lonname]#[latname]
2082	[file]; with the [lon],[lat] matrices
2083	[lonname],[latname]; names of the longitudes and latitudes variables
2084	values=[leglabels]\|[lonlatlims]\|[title]\|[graphk]\|[mapkind]
2085	[leglabels]: ',' separated list of names for the legend
2086	[lonlatlims]: ',' list of limits of the map [lonmin, latmin, lonmax, latmax] or None
2087	[title]: title of the plot ('!' for spaces)
2088	[graphk]: kind of the graphic
2089	[mapkind]: drawing coastaline ([proj],[res]) or None
2090	[proj]: projection
2091	* 'cyl', cilindric
2092	* 'lcc', lambert conformal
2093	[res]: resolution:
2094	* 'c', crude
2095	* 'l', low
2096	* 'i', intermediate
2097	* 'h', high
2098	* 'f', full
2099	obsevations= [obsfile],[obsname],[Tint],[null]
2100	[obsfile]: name fo the File with the observations as [t-step] [lat] [lon]
2101	[obsname]: name of the observations in the graph
2102	[Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2103	[null]: null value for the observed trajectory
2104	"""
2105
2106	fname = 'draw_trajectories'
2107
2108	if values == 'h':
2109	print fname + '_____________________________________________________________'
2110	print draw_trajectories.__doc__
2111	quit()
2112
2113	expectargs = '[leglabels]\|[lonlatlims]\|[title]\|[graphk]\|[mapkind]'
2114
2115	drw.check_arguments(fname,len(expectargs.split('\|')),values,'\|',expectargs)
2116
2117	trjfiles = trjfilens.split(',')
2118	leglabels = values.split('\|')[0]
2119	lonlatlims = values.split('\|')[1]
2120	title = values.split('\|')[2].replace('!',' ')
2121	graphk = values.split('\|')[3]
2122	mapkind = values.split('\|')[4]
2123
2124	Nfiles = len(trjfiles)
2125
2126	# Getting trajectotries
2127	##
2128
2129	lontrjvalues = []
2130	lattrjvalues = []
2131
2132	print ' ' + fname
2133	ifn = 0
2134	for ifile in trjfiles:
2135	filen = ifile.split('@')[0]
2136	Tint = ifile.split('@')[1]
2137
2138	print ' trajectory:',filen
2139
2140	if Tint != '-1':
2141	Tbeg = Tint
2142	Tend = ifile.split('@')[2]
2143	mapv = ifile.split('@')[3]
2144	else:
2145	mapv = ifile.split('@')[2]
2146
2147	if not os.path.isfile(filen):
2148	print errormsg
2149	print ' ' + fname + ": trajectory file '" + filen + "' does not exist !!"
2150	quit(-1)
2151
2152	# Charging longitude and latitude values
2153	##
2154	lonvals, latvals = drw.lonlat_values(mapv.split('#')[0], mapv.split('#')[1], \
2155	mapv.split('#')[2])
2156
2157	if ifn == 0: mapref = mapv
2158	ifn = ifn + 1
2159
2160	objfile = open(filen, 'r')
2161	trjtimev = []
2162	trjxv = []
2163	trjyv = []
2164
2165	for line in objfile:
2166	if line[0:1] != '#':
2167	trjtimev.append(int(line.split(' ')[0]))
2168	trjxv.append(int(line.split(' ')[1]))
2169	trjyv.append(int(line.split(' ')[2]))
2170
2171	objfile.close()
2172
2173	if Tint != '-1':
2174	lontrjvalues.append(lonvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
2175	lattrjvalues.append(latvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
2176	else:
2177	lontrjvalues.append(lonvals[trjyv[:], trjxv[:]])
2178	lattrjvalues.append(latvals[trjyv[:], trjxv[:]])
2179
2180	# lonlatlimits
2181	##
2182
2183	if lonlatlims == 'None':
2184	lonlatlimsv = None
2185	else:
2186	lonlatlimsv = np.zeros((4), dtype=np.float)
2187	lonlatlimsv[0] = np.float(lonlatlims.split(',')[0])
2188	lonlatlimsv[1] = np.float(lonlatlims.split(',')[1])
2189	lonlatlimsv[2] = np.float(lonlatlims.split(',')[2])
2190	lonlatlimsv[3] = np.float(lonlatlims.split(',')[3])
2191
2192	# lon/lat objects
2193	##
2194	objnc = NetCDFFile(mapref.split('#')[0])
2195	lonobj = objnc.variables[mapref.split('#')[1]]
2196	latobj = objnc.variables[mapref.split('#')[2]]
2197
2198	# map
2199	##
2200	if mapkind == 'None':
2201	mapkindv = None
2202	else:
2203	mapkindv = mapkind
2204
2205	if observations is None:
2206	obsname = None
2207	else:
2208	obsfile = observations.split(',')[0]
2209	obsname = observations.split(',')[1]
2210	Tint = observations.split(',')[2]
2211	null = np.float(observations.split(',')[3])
2212	print ' observational trajectory:',obsfile
2213
2214	if not os.path.isfile(obsfile):
2215	print errormsg
2216	print ' ' + fname + ": observations trajectory file '" + obsfile + \
2217	"' does not exist !!"
2218	quit(-1)
2219
2220	objfile = open(obsfile, 'r')
2221	obstrjtimev = []
2222	obstrjxv = []
2223	obstrjyv = []
2224
2225	for line in objfile:
2226	if line[0:1] != '#':
2227	lon = np.float(line.split(' ')[2])
2228	lat = np.float(line.split(' ')[1])
2229	if not lon == null and not lat == null:
2230	obstrjtimev.append(int(line.split(' ')[0]))
2231	obstrjxv.append(lon)
2232	obstrjyv.append(lat)
2233	else:
2234	obstrjtimev.append(int(line.split(' ')[0]))
2235	obstrjxv.append(None)
2236	obstrjyv.append(None)
2237
2238	objfile.close()
2239
2240	if Tint != '-1':
2241	Tint = int(observations.split(',')[2].split('@')[0])
2242	Tend = int(observations.split(',')[2].split('@')[1])
2243	lontrjvalues.append(obstrjxv[Tint:Tend+1])
2244	lattrjvalues.append(obstrjyv[Tint:Tend+1])
2245	else:
2246	lontrjvalues.append(obstrjxv[:])
2247	lattrjvalues.append(obstrjyv[:])
2248
2249	drw.plot_Trajectories(lontrjvalues, lattrjvalues, leglabels.split(','), \
2250	lonobj, latobj, lonlatlimsv, title, graphk, mapkindv, obsname)
2251
2252	objnc.close()
2253
2254	return
2255
2256	def draw_vals_trajectories(ncfile, values, variable):
2257	""" Function to draw values from the outputs from 'compute_tevolboxtraj'
2258	draw_vals_trajectories(ncfile, values, variable)
2259	ncfile= [ncfile] ',' list of files to use
2260	values= [statisticskind]:[Tint]:[labels]@[locleg]:[gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[gkind]
2261	[statisticskind]=[statistics][kind]
2262	[statistics]: which statistics to use, from: 'center', 'min', 'max', 'mean',
2263	'mean2', 'stdev'
2264	[kind]: 'box', 'circle' statistics taking the values from a box or a circle
2265	'trj': value following the trajectory
2266	[Tint]: [Tbeg]@[Tend] or None, interval of time to plot or -1 for all the times
2267	[labels]: ',' separated list of labels for the legend
2268	[locleg]: location of the legend (-1, autmoatic)
2269	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2270	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2271	9: 'upper center', 10: 'center'
2272	[gvarname]: name of the variable to appear in the graph
2273	[timetit]: title of the time axis (assumed x-axis, '\|' for spaces)
2274	[tkind]: kind of time to appear in the graph (assumed x-axis)
2275	'Nval': according to a given number of values as 'Nval',[Nval]
2276	'exct': according to an exact time unit as 'exct',[tunit];
2277	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2278	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2279	'l': milisecond
2280	[timefmt]: format of the time labels
2281	[title]: title of the graphic ('\|' for spaces)
2282	[gkind]: kind of graphical output
2283	variable= variable to use
2284	"""
2285
2286	fname = 'draw_vals_trajectories'
2287
2288	if values == 'h':
2289	print fname + '_____________________________________________________________'
2290	print draw_vals_trajectories.__doc__
2291	quit()
2292
2293	sims = ncfile.split(',')
2294
2295	if len(values.split(':')) != 9:
2296	print errormsg
2297	print ' ' + fname + ': wrong number of values!', len(values.split(':')), \
2298	'given 9 needed!!'
2299	print ' ',values.split(':')
2300	quit(-1)
2301
2302	statisticskind = values.split(':')[0]
2303	Tint = values.split(':')[1]
2304	labels = values.split(':')[2]
2305	gvarname = values.split(':')[3]
2306	timetit = values.split(':')[4].replace('\|',' ')
2307	tkind = values.split(':')[5]
2308	timefmt = values.split(':')[6]
2309	title = values.split(':')[7].replace('\|',' ')
2310	gkind = values.split(':')[8]
2311
2312	leglabels = labels.split('@')[0].split(',')
2313	locleg = int(labels.split('@')[1])
2314
2315	Nsims = len(sims)
2316
2317	if Tint != '-1':
2318	tini = np.float(Tint.split('@')[0])
2319	tend = np.float(Tint.split('@')[1])
2320	else:
2321	tini = -1.
2322	tend = -1.
2323
2324	vartimetrjv = []
2325
2326	print ' ' + fname
2327	for trjfile in sims:
2328	print ' ' + trjfile + ' ...'
2329	if not os.path.isfile(trjfile):
2330	print errormsg
2331	print ' ' + fname + ": trajectory file: '" + trjfile + \
2332	"' does not exist !!"
2333	quit(-1)
2334
2335	trjobj = NetCDFFile(trjfile, 'r')
2336	otim = trjobj.variables['time']
2337	if not trjobj.variables.has_key(statisticskind + '_' + variable):
2338	print errormsg
2339	print ' ' + fname + ": file '" + trjfile + "' does not have variable '"+\
2340	statisticskind + '_' + variable + "' !!"
2341	quit(-1)
2342	ovar = trjobj.variables[statisticskind + '_' + variable]
2343	dimt = otim.shape[0]
2344
2345	if trjfile == sims[0]:
2346	gunits = ovar.getncattr('units')
2347	lname = ovar.getncattr('long_name')
2348	tunits = otim.getncattr('units')
2349
2350	if tini != -1:
2351	tiniid = -1
2352	tendid = -1
2353	for itv in range(dimt):
2354	if otim[itv] <= tini and otim[itv+1] >= tini: tiniid = itv
2355	if otim[itv] <= tend and otim[itv+1] >= tend: tendid = itv
2356
2357	if tiniid == -1 or tendid == -1:
2358	print errormsg
2359	print ' ' + main + ' time interval ', tini,',',tend,' not found: ', \
2360	tendid, ',', tiniid, ' !!'
2361	print ' data interval [',otim[0], otim[dimt-1],']'
2362	quit(-1)
2363	dimt = tendid - tiniid + 1
2364
2365	else:
2366	dimt = otim.shape[0]
2367
2368	valsv = np.zeros((2,dimt), dtype=np.float)
2369	# Checking for time consistency
2370	if otim.getncattr('units') != tunits:
2371	print warnmsg
2372	print ' ' + fname + ': different time units in the plot!!'
2373	newtimes = drw.coincident_CFtimes(otim[:], tunits, otim.getncattr('units'))
2374	else:
2375	newtimes = otim[:]
2376
2377	if tini == -1:
2378	valsv[1,:] = newtimes[:]
2379	valsv[0,:] = ovar[:]
2380	else:
2381	valsv[1,:] = newtimes[tiniid:tendid+1]
2382	valsv[0,:] = ovar[tiniid:tendid+1]
2383
2384	vartimetrjv.append(valsv)
2385	trjobj.close()
2386
2387	drw.plot_TimeSeries(vartimetrjv, drw.units_lunits(gunits), tunits, \
2388	'val_trajectories_' + statisticskind, gvarname, timetit, tkind, timefmt, title,\
2389	leglabels, locleg, gkind)
2390
2391	def variable_values(values):
2392	""" Function to give back values for a given variable
2393	values= [varname] name of the variable
2394	"""
2395
2396	fname = 'variable_values'
2397
2398	values = drw.variables_values(values)
2399
2400	print fname,'values:',values
2401	print fname,'variable_name:',values[0]
2402	print fname,'standard_name:',values[1]
2403	print fname,'min,max:',str(values[2]) + ',' + str(values[3])
2404	print fname,'long_name:',values[4]
2405	print fname,'units:',values[5]
2406	print fname,'shad_colors:',values[6]
2407	print fname,'all_values:',drw.numVector_String(values,',')
2408
2409	return
2410
2411	####### ###### ##### #### ### ## #
2412
2413	ngraphics = "'" + drw.numVector_String(namegraphics, "', '") + "'"
2414
2415	### Options
2416	##string_operation="operation to make: " + '\n' + " out, output values -S inidim1,[inidim2,...]:enddim1,[enddim2,...]"
2417	string_operation="""operation to make:
2418	draw_topo_geogrid, draws topography from a WPS geo_em.d[nn].nc: -S [minTopo],[maxTopo]:[SW_lon],[SW_lat],[NE_lon],[NE_lat]:[title]:[graphic_kind]:[projection],[res_coastline]
2419	draw_2D_shad_cont, draws two 2D fields, first with shading second with contour lines: -v [varns],[varnc] -S [vnamefs],[vnamefc],[dimxvn],[dimyvn],[colorbar],[ckind],[clabfmt],[sminv]:[smaxv],[sminc]:[smaxv]:[Nlev],[figt],[kindfig],[reverse]
2420	[ckind]:
2421	'cmap': as it gets from colorbar
2422	'fixc,[colname]': fixed color [colname], all stright lines
2423	'fixsignc,[colname]': fixed color [colname], >0 stright, <0 dashed line
2424	"""
2425
2426	#print string_operation
2427
2428	parser = OptionParser()
2429	parser.add_option("-f", "--netCDF_file", dest="ncfile",
2430	help="file to use", metavar="FILE")
2431	parser.add_option("-o", "--operation", type='choice', dest="operation",
2432	choices=namegraphics,
2433	help="operation to make: " + ngraphics, metavar="OPER")
2434	parser.add_option("-S", "--valueS", dest="values",
2435	help="[WHEN APPLICABLE] values to use according to the operation", metavar="VALUES")
2436	parser.add_option("-v", "--variable", dest="varname",
2437	help="[WHEN APPLICABLE] variable to check", metavar="VAR")
2438
2439	(opts, args) = parser.parse_args()
2440
2441	####### #######
2442	## MAIN
2443	#######
2444
2445	# Not checking file operation
2446	Notcheckingfile = ['draw_2D_shad_cont', 'draw_2D_shad_cont_time', \
2447	'draw_2D_shad_line', 'draw_2D_shad_line_time', 'draw_lines', 'draw_lines_time', \
2448	'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories', \
2449	'variable_values']
2450
2451	####### ###### ##### #### ### ## #
2452	errormsg='ERROR -- error -- ERROR -- error'
2453
2454	varn=opts.varname
2455	oper=opts.operation
2456
2457	if opts.ncfile is not None and not os.path.isfile(opts.ncfile) and \
2458	not drw.searchInlist(Notcheckingfile, oper):
2459	print errormsg
2460	print ' ' + main + ': File ' + opts.ncfile + ' does not exist !!'
2461	quit(-1)
2462
2463	if oper == 'create_movie':
2464	create_movie(opts.ncfile, opts.values, opts.varname)
2465	elif oper == 'draw_2D_shad':
2466	draw_2D_shad(opts.ncfile, opts.values, opts.varname)
2467	elif oper == 'draw_2D_shad_time':
2468	draw_2D_shad_time(opts.ncfile, opts.values, opts.varname)
2469	elif oper == 'draw_2D_shad_cont':
2470	draw_2D_shad_cont(opts.ncfile, opts.values, opts.varname)
2471	elif oper == 'draw_2D_shad_cont_time':
2472	draw_2D_shad_cont_time(opts.ncfile, opts.values, opts.varname)
2473	elif oper == 'draw_2D_shad_line':
2474	draw_2D_shad_line(opts.ncfile, opts.values, opts.varname)
2475	elif oper == 'draw_2D_shad_line_time':
2476	draw_2D_shad_line_time(opts.ncfile, opts.values, opts.varname)
2477	elif oper == 'draw_Neighbourghood_evol':
2478	draw_Neighbourghood_evol(opts.ncfile, opts.values, opts.varname)
2479	elif oper == 'draw_lines':
2480	draw_lines(opts.ncfile, opts.values, opts.varname)
2481	elif oper == 'draw_lines_time':
2482	draw_lines_time(opts.ncfile, opts.values, opts.varname)
2483	elif oper == 'draw_timeSeries':
2484	draw_timeSeries(opts.ncfile, opts.values, opts.varname)
2485	elif oper == 'draw_topo_geogrid':
2486	draw_topo_geogrid(opts.ncfile, opts.values)
2487	elif oper == 'draw_topo_geogrid_boxes':
2488	draw_topo_geogrid_boxes(opts.ncfile, opts.values)
2489	elif oper == 'draw_trajectories':
2490	draw_trajectories(opts.ncfile, opts.values, opts.varname)
2491	elif oper == 'draw_vals_trajectories':
2492	draw_vals_trajectories(opts.ncfile, opts.values, opts.varname)
2493	elif oper == 'list_graphics':
2494	# From: http://www.diveintopython.net/power_of_introspection/all_together.html
2495	import drawing as myself
2496	object = myself
2497	for opern in namegraphics:
2498	if opern != 'list_graphics':
2499	print opern + '_______ ______ _____ ____ ___ __ _'
2500	print getattr(object, opern).__doc__
2501	elif oper == 'variable_values':
2502	variable_values(opts.values)
2503	else:
2504	print errormsg
2505	print ' ' + main + ": the graphic '" + oper + "' is not ready !!"
2506	print errormsg
2507	quit()

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

Download in other formats: