Context Navigation

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

Last change on this file since 394 was 383, checked in by lfita, 10 years ago
Fixing issues on taking values for the dimensions
File size: 93.8 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 tat 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	print ' ' + fname + ': Lluis shapes valshad:',valshad.shape,'valcont:',valcont.shape
584
585	drw.plot_2D_shadow_contour(valshad, valcont, vnamesfig, odimxv, odimyv, odimxu, \
586	odimyu, dimnamesv, colbarn, countkind, countlfmt, shading_nx, levels_cont, \
587	varunits, figtitle, figkind, revals, mapvalue)
588
589	return
590
591	def draw_2D_shad_cont_time(ncfile, values, varn):
592	""" plotting two fields, one with shading and the other with contour lines being
593	one of the dimensions of time characteristics
594	draw_2D_shad_cont(ncfile, values, varn)
595	ncfile= [ncfilevars],[ncfilevarc] files to use (one value, same file)
596	values=[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:[figt]:[kindfig]:[reverse]:[timevals]:[mapv]
597	[vnamefs],[vnamefc]: Name in the figure of the shaded and the contour variables
598	[dimvals/c]: list of [dimname]\|[value] telling at which dimension of the
599	variable a given value is required (no dimension name, all the length)
600	[dimx/yvn]: ',' list with the name of the variables with the values of the dimensions
601	[colorbar]: name of the color bar
602	[ckind]: kind of contours
603	'cmap': as it gets from colorbar
604	'fixc,[colname]': fixed color [colname], all stright lines
605	'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed line
606	[clabfmt]: format of the labels in the contour (None, also possible)
607	[smin/axv]: minimum and maximum value for the shading
608	[sminc]:[smaxv]:[Nlev]: minimum, maximum and number of values for the contour
609	[figt]: title of the figure ('\|' for spaces)
610	[kindfig]: kind of figure
611	[reverse]: modification to the dimensions:
612	'transposed': transpose matrices
613	'flip',[x/y]: flip only the dimension [x] or [y]
614	[timevals]: [timen]\|[units]\|[kind]\|[tfmt]\|[label] time labels characteristics
615	[timen]; name of the time variable
616	[units]; units string according to CF conventions ([tunits] since
617	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
618	[kind]; kind of output
619	'Nval': according to a given number of values as 'Nval',[Nval]
620	'exct': according to an exact time unit as 'exct',[tunit];
621	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
622	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
623	'l': milisecond
624	[tfmt]; desired format
625	[label]; label at the graph ('!' for spaces)
626	[mapv]: map characteristics: [proj],[res]
627	see full documentation: http://matplotlib.org/basemap/
628	[proj]: projection
629	* 'cyl', cilindric
630	* 'lcc', lamvbert conformal
631	[res]: resolution:
632	* 'c', crude
633	* 'l', low
634	* 'i', intermediate
635	* 'h', high
636	* 'f', full
637	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'
638	varn= [varsn],[varcn] name of the variable to plot with shading variable with contour
639	"""
640
641	fname = 'draw_2D_shad_cont_time'
642	if values == 'h':
643	print fname + '_____________________________________________________________'
644	print draw_2D_shad_cont_time.__doc__
645	quit()
646
647	expectargs = ['[vnamefs]', '[dimvals]', '[dimvalc]', '[dimxvn]', '[dimyvn]', \
648	'[colorbar]', '[ckind]', '[clabfmt]', '[sminv],[smaxv]', \
649	'[sminc],[smaxv],[Nlev]', '[figt]', '[kindfig]', '[reverse]', '[timevals]', \
650	'[mapv]']
651
652	drw.check_arguments(fname,len(expectargs),values,':',expectargs)
653
654	vnamesfig = values.split(':')[0].split(',')
655	dimvals= values.split(':')[1].replace('\|',':')
656	dimvalc= values.split(':')[2].replace('\|',':')
657	vdimxn = values.split(':')[3]
658	vdimyn = values.split(':')[4]
659	colbarn = values.split(':')[5]
660	countkind = values.split(':')[6]
661	countlabelfmt = values.split(':')[7]
662	shadminmax = values.split(':')[8]
663	contlevels = values.split(':')[9]
664	figtitle = values.split(':')[10].replace('\|',' ')
665	figkind = values.split(':')[11]
666	revals = values.split(':')[12]
667	timevals = values.split(':')[13]
668	mapvalue = values.split(':')[14]
669
670	if2filenames = ncfile.find(',')
671
672	if if2filenames != -1:
673	ncfiles = ncfile.split(',')[0]
674	ncfilec = ncfile.split(',')[1]
675	else:
676	ncfiles = ncfile
677	ncfilec = ncfile
678
679	if not os.path.isfile(ncfiles):
680	print errormsg
681	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
682	quit(-1)
683
684	if not os.path.isfile(ncfilec):
685	print errormsg
686	print ' ' + fname + ': contour file "' + ncfilec + '" does not exist !!'
687	quit(-1)
688
689	objsf = NetCDFFile(ncfiles, 'r')
690	objcf = NetCDFFile(ncfilec, 'r')
691
692	varns = varn.split(',')[0]
693	varnc = varn.split(',')[1]
694
695	if not objsf.variables.has_key(varns):
696	print errormsg
697	print ' ' + fname + ': shading file "' + ncfiles + \
698	'" does not have variable "' + varns + '" !!'
699	quit(-1)
700
701	if not objcf.variables.has_key(varnc):
702	print errormsg
703	print ' ' + fname + ': contour file "' + ncfilec + \
704	'" does not have variable "' + varnc + '" !!'
705	quit(-1)
706
707	# Variables' values
708	objvars = objsf.variables[varns]
709	objvarc = objcf.variables[varnc]
710
711	if len(objvars.shape) != len(objvarc.shape):
712	print errormsg
713	print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \
714	objvars.shape, 'different than contour variable "' + varnc + '": ', \
715	objvarc.shape,' !!!'
716	quit(-1)
717
718	for idim in range(len(objvars.shape)):
719	if objvars.shape[idim] != objvarc.shape[idim]:
720	print errormsg
721	print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \
722	objvars.shape, 'different than contour variable "' + varnc + '": ', \
723	objvarc.shape,' !!!'
724	quit(-1)
725
726	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
727	valcont, dimscont = drw.slice_variable(objvarc, dimvalc.replace(',','\|'))
728
729	# Dimensions names
730	## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
731	## dimnamesv = []
732	## for idd in range(len(objvars.dimensions)):
733	## cutdim = False
734	## for idc in range(len(dimvals.split(','))):
735	## dimcutn = dimvals.split(',')[idc].split(':')[0]
736	## print objvars.dimensions[idd], dimcutn
737	## if objvars.dimensions[idd] == dimcutn:
738	## cutdim = True
739	## break
740	## if not cutdim: dimnamesv.append(objvars.dimensions[idd])
741	dimnamesv = [vdimyn, vdimxn]
742
743	varunits = []
744	varunits.append(objvars.getncattr('units'))
745	varunits.append(objvarc.getncattr('units'))
746
747	if not objsf.variables.has_key(vdimxn):
748	print errormsg
749	print ' ' + fname + ': shading file "' + ncfiles + \
750	'" does not have dimension variable "' + vdimxn + '" !!'
751	quit(-1)
752	if not objsf.variables.has_key(vdimyn):
753	print errormsg
754	print ' ' + fname + ': shading file "' + ncfiles + \
755	'" does not have dimension variable "' + vdimyn + '" !!'
756	quit(-1)
757
758	timename = timevals.split('\|')[0]
759	timeunit = timevals.split('\|')[1].replace('!',' ')
760	timekind = timevals.split('\|')[2]
761	timefmt = timevals.split('\|')[3]
762	timelabel = timevals.split('\|')[4].replace('!',' ')
763
764	if vdimxn == timename:
765	timevals = objsf.variables[vdimxn][:]
766	dimt = 'x'
767	ovalaxis = objsf.variables[vdimyn]
768	ovalu = ovalaxis.getncattr('units')
769	elif vdimyn == timename:
770	timevals = objsf.variables[vdimyn][:]
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	if len(ovalaxis.shape) <= 2:
782	ovalaxisv = ovalaxis[:]
783
784	elif len(ovalaxis.shape) == 3:
785	ovalaxisv = ovalaxis[0,:]
786	else:
787	print errormsg
788	print ' ' + fname + ': shape of dimension variable:', ovalaxis.shape, \
789	' not ready!!'
790	quit(-1)
791
792	if countlabelfmt == 'None':
793	countlfmt = None
794	else:
795	countlfmt = countlabelfmt
796
797	shading_nx = np.zeros((2), dtype=np.float)
798	shading_nx[0] = np.float(shadminmax.split(',')[0])
799	shading_nx[1] = np.float(shadminmax.split(',')[1])
800
801	clevmin = np.float(contlevels.split(',')[0])
802	clevmax = np.float(contlevels.split(',')[1])
803	Nclevels = int(contlevels.split(',')[2])
804
805	levels_cont = drw.pretty_int(clevmin, clevmax, Nclevels)
806
807	if len(levels_cont) <= 1:
808	print warnmsg
809	print ' ' + fname + ': wrong contour levels:', levels_cont,' !!'
810	del(levels_cont)
811	levels_cont = np.zeros((Nclevels), dtype=np.float)
812	levels_cont = np.arange(7)*(clevmax - clevmin)/(Nclevels-1)
813	print ' generating default ones: ',levels_cont
814
815	if mapvalue == 'None': mapvalue = None
816
817	drw.plot_2D_shadow_contour_time(valshad, valcont, vnamesfig, ovalaxisv, \
818	timevals, timepos, timelabels, ovalu, timelabel, dimt, dimnamesv, colbarn, \
819	countkind, countlfmt, shading_nx, levels_cont, varunits, figtitle, figkind, \
820	revals, mapvalue)
821
822	return
823
824	def draw_2D_shad_line(ncfile, values, varn):
825	""" plotting a fields with shading and another with line
826	draw_2D_shad_line(ncfile, values, varn)
827	ncfile= [ncfiles],[ncfilel] file to use for the shading and for the line
828	values=[vnamefs],[vnamefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar],[colline]:[sminv],[smaxv]:[figt]:
829	[kindfig]:[reverse]:[mapv]:[close]
830	[vnamefs]: Name in the figure of the variable to be shaded
831	[vnamefl]: Name in the figure of the variable to be lined
832	[dimvals]: ',' list of [dimname]\|[value] telling at which dimension of the
833	variable a given value is required (-1, all the length)
834	[dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
835	[colorbar]: name of the color bar
836	[colline]: name of the color for the line
837	[smin/axv]: minimum and maximum value for the shading
838	[figt]: title of the figure ('\|' for spaces)
839	[kindfig]: kind of figure
840	[reverse]: Transformation of the values
841	* 'transpose': reverse the axes (x-->y, y-->x)
842	* 'flip'@[x/y]: flip the axis x or y
843	[mapv]: map characteristics: [proj],[res]
844	see full documentation: http://matplotlib.org/basemap/
845	[proj]: projection
846	* 'cyl', cilindric
847	* 'lcc', lamvbert conformal
848	[res]: resolution:
849	* 'c', crude
850	* 'l', low
851	* 'i', intermediate
852	* 'h', high
853	* 'f', full
854	valules= 'rh:z\|-1,x\|-1:z\|-1,x\|-1:lat:pressure:BuPu:0.,100.:rh:pdf:flip@y:None'
855	varn= [varsn],[varnl] name of the variable to plot with shading and with line
856	"""
857
858	fname = 'draw_2D_shad_line'
859	if values == 'h':
860	print fname + '_____________________________________________________________'
861	print draw_2D_shad_line.__doc__
862	quit()
863
864	farguments = ['[vnamefs],[vnamefl]', '[dimvals]', '[dimxvn]', '[dimyvn]', \
865	'[colorbar],[colline]', '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', \
866	'[mapv]', '[close]']
867	drw.check_arguments(fname,11,values,':',farguments)
868
869	vnamesfig = values.split(':')[0].split(',')[0]
870	dimvals= values.split(':')[1].replace('\|',':')
871	vdimxn = values.split(':')[2]
872	vdimyn = values.split(':')[3]
873	colbarn = values.split(':')[4].split(',')[0]
874	shadminmax = values.split(':')[5]
875	figtitle = values.split(':')[6].replace('\|',' ')
876	figkind = values.split(':')[7]
877	revals = values.split(':')[8]
878	mapvalue = values.split(':')[9]
879	# varn = values.split(':')[10]
880
881	ncfiles = ncfile.split(',')[0]
882
883	if not os.path.isfile(ncfiles):
884	print errormsg
885	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
886	quit(-1)
887
888	objsf = NetCDFFile(ncfiles, 'r')
889
890	varns = varn.split(',')[0]
891
892	if not objsf.variables.has_key(varns):
893	print errormsg
894	print ' ' + fname + ': shading file "' + ncfiles + \
895	'" does not have variable "' + varns + '" !!'
896	quit(-1)
897
898	# Variables' values
899	objvars = objsf.variables[varns]
900
901	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
902
903	# Dimensions names
904	## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
905	## dimnamesv = []
906	## for idd in range(len(objvars.dimensions)):
907	## cutdim = False
908	## for idc in range(len(dimvals.split(','))):
909	## dimcutn = dimvals.split(',')[idc].split(':')[0]
910	## print objvars.dimensions[idd], dimcutn
911	## if objvars.dimensions[idd] == dimcutn:
912	## cutdim = True
913	## break
914	## if not cutdim: dimnamesv.append(objvars.dimensions[idd])
915	dimnamesv = [vdimyn, vdimxn]
916
917	varunits = objvars.getncattr('units')
918
919	if not objsf.variables.has_key(vdimxn):
920	print errormsg
921	print ' ' + fname + ': shading file "' + ncfiles + \
922	'" does not have dimension variable "' + vdimxn + '" !!'
923	quit(-1)
924	if not objsf.variables.has_key(vdimyn):
925	print errormsg
926	print ' ' + fname + ': shading file "' + ncfiles + \
927	'" does not have dimensino variable "' + vdimyn + '" !!'
928	quit(-1)
929
930	objdimx = objsf.variables[vdimxn]
931	objdimy = objsf.variables[vdimyn]
932	odimxu = objdimx.getncattr('units')
933	odimyu = objdimy.getncattr('units')
934
935	if len(objdimx.shape) <= 2:
936	# odimxv = objdimx[valshad.shape]
937	# odimyv = objdimy[valshad.shape]
938	odimxv = objdimx[:]
939	odimyv = objdimy[:]
940
941	elif len(objdimx.shape) == 3:
942	# dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])]
943	# odimxv = objdimx[tuple(dimcut)]
944	# odimyv = objdimy[tuple(dimcut)]
945	odimxv = objdimx[0,:]
946	odimyv = objdimy[0,:]
947	else:
948	print errormsg
949	print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
950	' not ready!!'
951	quit(-1)
952
953	shading_nx = np.zeros((2), dtype=np.float)
954	shading_nx[0] = np.float(shadminmax.split(',')[0])
955	shading_nx[1] = np.float(shadminmax.split(',')[1])
956
957	if mapvalue == 'None': mapvalue = None
958
959	# line plot
960	##
961	ncfilel = ncfile.split(',')[1]
962	vnamelfig = values.split(':')[0].split(',')[1]
963	varnl = varn.split(',')[1]
964	colline = values.split(':')[4].split(',')[1]
965
966	objlf = NetCDFFile(ncfilel,'r')
967	objlvar = objlf.variables[varnl]
968
969	linevals = objlvar[:]
970	varlunits = objlvar.units
971
972	drw.plot_2D_shadow_line(valshad, linevals, vnamesfig, vnamelfig, odimxv, odimyv, \
973	odimxu, odimyu, dimnamesv, colbarn, colline, shading_nx, varunits, varlunits, \
974	figtitle, figkind, revals, mapvalue, True)
975
976	objsf.close()
977	objlf.close()
978
979	return
980
981	def draw_2D_shad_line_time(ncfile, values, varn):
982	""" plotting a fields with shading and a line with time values
983	draw_2D_shad_line(ncfile, values, varn)
984	ncfile= [ncfiles],[ncfilel] files to use to draw with shading and the line
985	values= [vnamefs],[vanemefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
986	[kindfig]:[reverse]:[timevals]:[close]
987	[vnamefs]: Name in the figure of the variable to be shaded
988	[vnamefl]: Name in the figure of the variable to be lined
989	[dimvals]: ',' list of [dimname]\|[value] telling at which dimension of the
990	variable a given value is required (-1, all the length)
991	[dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
992	[colorbar]: name of the color bar
993	[smin/axv]: minimum and maximum value for the shading
994	[figt]: title of the figure ('\|' for spaces)
995	[kindfig]: kind of figure
996	[reverse]: Transformation of the values
997	* 'transpose': reverse the axes (x-->y, y-->x)
998	* 'flip'@[x/y]: flip the axis x or y
999	[timevals]: [timen]\|[units]\|[kind]\|[tfmt]\|[label]\|[timeaxis] time labels characteristics
1000	[timen]; name of the time variable
1001	[units]; units string according to CF conventions ([tunits] since
1002	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
1003	[kind]; kind of output
1004	'Nval': according to a given number of values as 'Nval',[Nval]
1005	'exct': according to an exact time unit as 'exct',[tunit];
1006	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1007	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1008	'l': milisecond
1009	[tfmt]; desired format
1010	[label]; label at the graph ('!' for spaces)
1011	[close]: should figure be closed (finished)
1012	values='dtcon,prc:Time\|-1,bottom_top\|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly\|'
1013	'dtcon:pdf:transpose:time\|hours!since!1949-12-01\|exct,2,d\|%d\|date!([DD])\|x:True
1014	varn= [varsn].[varln] name of the variable to plot with shading and to plot with line
1015	"""
1016	fname = 'draw_2D_shad_line_time'
1017	if values == 'h':
1018	print fname + '_____________________________________________________________'
1019	print draw_2D_shad__line_time.__doc__
1020	quit()
1021
1022	farguments = ['[vnamefs],[vanemefl]', '[dimvals]', '[dimxvn]', '[dimyvn]', \
1023	'[colorbar]', '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', \
1024	'[timevals]', '[close]']
1025	drw.check_arguments(fname,11,values,':',farguments)
1026
1027	vnamesfig = values.split(':')[0].split(',')[0]
1028	dimvals= values.split(':')[1].replace('\|',':')
1029	vdimxn = values.split(':')[2]
1030	vdimyn = values.split(':')[3]
1031	colbarn = values.split(':')[4]
1032	shadminmax = values.split(':')[5]
1033	figtitle = values.split(':')[6].replace('\|',' ')
1034	figkind = values.split(':')[7]
1035	revals = values.split(':')[8]
1036	timevals = values.split(':')[9]
1037	close = values.split(':')[10]
1038
1039	ncfiles = ncfile.split(',')[0]
1040
1041	if not os.path.isfile(ncfiles):
1042	print errormsg
1043	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
1044	quit(-1)
1045
1046	objsf = NetCDFFile(ncfiles, 'r')
1047
1048	varns = varn.split(',')[0]
1049
1050	if not objsf.variables.has_key(varns):
1051	print errormsg
1052	print ' ' + fname + ': shading file "' + ncfiles + \
1053	'" does not have variable "' + varns + '" !!'
1054	quit(-1)
1055
1056	# Variables' values
1057	objvars = objsf.variables[varns]
1058
1059	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
1060
1061	dimnamesv = [vdimyn, vdimxn]
1062
1063	varunits = objvars.getncattr('units')
1064
1065	if not objsf.variables.has_key(vdimxn):
1066	print errormsg
1067	print ' ' + fname + ': shading file "' + ncfiles + \
1068	'" does not have dimension variable "' + vdimxn + '" !!'
1069	quit(-1)
1070	if not objsf.variables.has_key(vdimyn):
1071	print errormsg
1072	print ' ' + fname + ': shading file "' + ncfiles + \
1073	'" does not have dimensino variable "' + vdimyn + '" !!'
1074	quit(-1)
1075
1076	objdimx = objsf.variables[vdimxn]
1077	objdimy = objsf.variables[vdimyn]
1078	odimxu = objdimx.getncattr('units')
1079	odimyu = objdimy.getncattr('units')
1080
1081	if len(objdimx.shape) <= 2:
1082	odimxv = objdimx[:]
1083	odimyv = objdimy[:]
1084
1085	elif len(objdimx.shape) == 3:
1086	odimxv = objdimx[0,:]
1087	odimyv = objdimy[0,:]
1088	else:
1089	print errormsg
1090	print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
1091	' not ready!!'
1092	quit(-1)
1093
1094	timename = timevals.split('\|')[0]
1095	timeunit = timevals.split('\|')[1].replace('!',' ')
1096	timekind = timevals.split('\|')[2]
1097	timefmt = timevals.split('\|')[3]
1098	timelabel = timevals.split('\|')[4].replace('!',' ')
1099
1100	if vdimxn == timename:
1101	odimxv = objsf.variables[vdimxn][:]
1102	odimxu = timelabel
1103	timeaxis = 'x'
1104	odimyv = objsf.variables[vdimyn]
1105	odimyu = odimyv.getncattr('units')
1106	timepos, timelabels = drw.CFtimes_plot(odimxv, timeunit, timekind, timefmt)
1107	elif vdimyn == timename:
1108	odimyv = objsf.variables[vdimyn][:]
1109	odimyu = timelabel
1110	timeaxis = 'y'
1111	odimxv = objsf.variables[vdimxn]
1112	odimxu = odimxv.getncattr('units')
1113	timepos, timelabels = drw.CFtimes_plot(odimyv, timeunit, timekind, timefmt)
1114	else:
1115	print errormsg
1116	print ' ' + fname + ": time variable '" + timename + "' not found!!"
1117	quit(-1)
1118
1119	shading_nx = np.zeros((2), dtype=np.float)
1120	shading_nx[0] = np.float(shadminmax.split(',')[0])
1121	shading_nx[1] = np.float(shadminmax.split(',')[1])
1122
1123	closeval = drw.Str_Bool(close)
1124
1125	drw.plot_2D_shadow_time(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu, \
1126	dimnamesv, colbarn, shading_nx, varunits, figtitle, figkind, revals, timeaxis, \
1127	timepos, timelabels, False)
1128
1129	# Line values
1130	##
1131	ncfilel = ncfile.split(',')[1]
1132
1133	vnamelfig = values.split(':')[0].split(',')[1]
1134	varnl = varn.split(',')[1]
1135
1136	objlf = NetCDFFile(ncfilel,'r')
1137	objlvar = objlf.variables[varnl]
1138
1139	linevals = objlvar[:]
1140	if reva0 == 'tranpose':
1141	plt.plot (linevals, odimxv, '-', color='k')
1142	else:
1143	plt.plot (odimxv, linevals, '-', color='k')
1144
1145	objsf.close()
1146	objsl.close()
1147
1148	return
1149
1150	def draw_topo_geogrid(ncfile, values):
1151	""" plotting geo_em.d[nn].nc topography from WPS files
1152	draw_topo_geogrid(ncfile, values)
1153	ncfile= geo_em.d[nn].nc file to use
1154	values= [minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]
1155	[min/max]Topo: minimum and maximum values of topography to draw
1156	lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] or None
1157	title: title of the graph
1158	graphic_kind: kind of figure (jpg, pdf, png)
1159	mapvalues: map characteristics [proj],[res]
1160	see full documentation: http://matplotlib.org/basemap/
1161	[proj]: projection
1162	* 'cyl', cilindric
1163	* 'lcc', lambert conformal
1164	[res]: resolution:
1165	* 'c', crude
1166	* 'l', low
1167	* 'i', intermediate
1168	* 'h', high
1169	* 'f', full
1170	"""
1171	fname = 'draw_topo_geogrid'
1172
1173	if values == 'h':
1174	print fname + '_____________________________________________________________'
1175	print draw_topo_geogrid.__doc__
1176	quit()
1177
1178	expectargs = ['[minTopo]','[maxTopo]', '[lonlatL]', '[title]', '[graphic_kind]', \
1179	'[mapvalues]']
1180
1181	drw.check_arguments(fname,5,values,':',expectargs)
1182
1183	mintopo = values.split(':')[0].split(',')[0]
1184	maxtopo = values.split(':')[0].split(',')[1]
1185
1186	lonlatLS = values.split(':')[1]
1187	lonlatLv = lonlatLS.split(',')[0]
1188
1189	if lonlatLv == 'None':
1190	lonlatL = None
1191	else:
1192	lonlatL = np.zeros((4), dtype=np.float)
1193	lonlatL[0] = np.float(lonlatLS.split(',')[0])
1194	lonlatL[1] = np.float(lonlatLS.split(',')[1])
1195	lonlatL[2] = np.float(lonlatLS.split(',')[2])
1196	lonlatL[3] = np.float(lonlatLS.split(',')[3])
1197
1198	grtit = values.split(':')[2]
1199	kindfig = values.split(':')[3]
1200	mapvalues = values.split(':')[4]
1201
1202	if not os.path.isfile(ncfile):
1203	print errormsg
1204	print ' ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1205	quit(-1)
1206
1207	objdomf = NetCDFFile(ncfile, 'r')
1208
1209	objhgt = objdomf.variables['HGT_M']
1210	objlon = objdomf.variables['XLONG_M']
1211	objlat = objdomf.variables['XLAT_M']
1212
1213	topography = objhgt[0,:,:]
1214
1215	drw.plot_topo_geogrid(topography, objlon, objlat, mintopo, maxtopo, lonlatL, \
1216	grtit, kindfig, mapvalues, True)
1217
1218	objdomf.close()
1219
1220	return
1221
1222	def draw_topo_geogrid_boxes(ncfiles, values):
1223	""" plotting different geo_em.d[nn].nc topography from WPS files
1224	draw_topo_geogrid_boxes(ncfile, values)
1225	ncfiles= ',' list of geo_em.d[nn].nc files to use (fisrt as topographyc reference)
1226	values= [minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]:[labels]
1227	[min/max]Topo: minimum and maximum values of topography to draw
1228	lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] or None
1229	title: title of the graph
1230	graphic_kind: kind of figure (jpg, pdf, png)
1231	mapvalues: map characteristics [proj],[res]
1232	see full documentation: http://matplotlib.org/basemap/
1233	[proj]: projection
1234	* 'cyl', cilindric
1235	* 'lcc', lambert conformal
1236	[res]: resolution:
1237	* 'c', crude
1238	* 'l', low
1239	* 'i', intermediate
1240	* 'h', high
1241	* 'f', full
1242	labels= labels to write in the graph
1243	"""
1244	# import matplotlib as mpl
1245	# mpl.use('Agg')
1246	import matplotlib.pyplot as plt
1247
1248	fname = 'draw_topo_geogrid_boxes'
1249
1250	if values == 'h':
1251	print fname + '_____________________________________________________________'
1252	print draw_topo_geogrid_boxes.__doc__
1253	quit()
1254
1255	mintopo = values.split(':')[0].split(',')[0]
1256	maxtopo = values.split(':')[0].split(',')[1]
1257
1258	lonlatLS = values.split(':')[1]
1259	lonlatLv = lonlatLS.split(',')[0]
1260
1261	if lonlatLv == 'None':
1262	lonlatL = None
1263	else:
1264	lonlatL = np.zeros((4), dtype=np.float)
1265	lonlatL[0] = np.float(lonlatLS.split(',')[0])
1266	lonlatL[1] = np.float(lonlatLS.split(',')[1])
1267	lonlatL[2] = np.float(lonlatLS.split(',')[2])
1268	lonlatL[3] = np.float(lonlatLS.split(',')[3])
1269
1270	grtit = values.split(':')[2]
1271	kindfig = values.split(':')[3]
1272	mapvalues = values.split(':')[4]
1273	labels = values.split(':')[5]
1274
1275	ncfile = ncfiles.split(',')[0]
1276	if not os.path.isfile(ncfile):
1277	print errormsg
1278	print ' ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1279	quit(-1)
1280
1281	objdomf = NetCDFFile(ncfile, 'r')
1282
1283	objhgt = objdomf.variables['HGT_M']
1284	objlon0 = objdomf.variables['XLONG_M']
1285	objlat0 = objdomf.variables['XLAT_M']
1286
1287	topography = objhgt[0,:,:]
1288
1289	Nfiles = len(ncfiles.split(','))
1290	boxlabels = labels.split(',')
1291
1292	Xboxlines = []
1293	Yboxlines = []
1294
1295	for ifile in range(Nfiles):
1296	ncfile = ncfiles.split(',')[ifile]
1297	# print ifile, ncfile
1298	if not os.path.isfile(ncfile):
1299	print errormsg
1300	print ' ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1301	quit(-1)
1302
1303	objdomfi = NetCDFFile(ncfile, 'r')
1304
1305	objlon = objdomfi.variables['XLONG_M']
1306	objlat = objdomfi.variables['XLAT_M']
1307
1308	dx = objlon.shape[2]
1309	dy = objlon.shape[1]
1310
1311	Xboxlines.append(objlon[0,0,:])
1312	Yboxlines.append(objlat[0,0,:])
1313	Xboxlines.append(objlon[0,dy-1,:])
1314	Yboxlines.append(objlat[0,dy-1,:])
1315	Xboxlines.append(objlon[0,:,0])
1316	Yboxlines.append(objlat[0,:,0])
1317	Xboxlines.append(objlon[0,:,dx-1])
1318	Yboxlines.append(objlat[0,:,dx-1])
1319
1320	objdomfi.close()
1321
1322	drw.plot_topo_geogrid_boxes(topography, Xboxlines, Yboxlines, boxlabels, \
1323	objlon0, objlat0, mintopo, maxtopo, lonlatL, grtit, kindfig, mapvalues, True)
1324
1325	objdomf.close()
1326
1327	return
1328
1329	def movievalslice(origslice, dimmovien, framenum):
1330	""" Function to provide variable slice according to a geneation of a movie
1331	movievals(origslice, dimmovien, framenum)
1332	[origslice]= slice original as [dimname1]\|[val1],[...,[dimnameN]\|[valN]]
1333	([val] = -1, full length)
1334	[dimmovien]= name of the dimension to produce the movie
1335	[framenum]= value of the frame to substitue in [origslice] as
1336	[dimmovien]\|[framenum]
1337	>>> movievalslice('East_West\|-1,North_South\|-1,Time\|2','Time',0)
1338	East_West\|-1,North_South\|-1,Time\|0
1339	"""
1340
1341	fname = 'movievalslice'
1342
1343	if origslice == 'h':
1344	print fname + '_____________________________________________________________'
1345	print movievalslice.__doc__
1346	quit()
1347
1348	dims = origslice.split(',')
1349
1350	movieslice = ''
1351	idim = 0
1352
1353	for dimn in dims:
1354	dn = dimn.split('\|')[0]
1355	if dn == dimmovien:
1356	movieslice = movieslice + dn + '\|' + str(framenum)
1357	else:
1358	movieslice = movieslice + dimn
1359	if idim < len(dims)-1: movieslice = movieslice + ','
1360
1361	idim = idim + 1
1362
1363	return movieslice
1364
1365	class Capturing(list):
1366	""" Class to capture function output as a list
1367	from: http://stackoverflow.com/questions/16571150/how-to-capture-stdout-output-from-a-python-function-call
1368	"""
1369	# from cStringIO import StringIO
1370
1371	def __enter__(self):
1372	self._stdout = sys.stdout
1373	sys.stdout = self._stringio = StringIO()
1374	return self
1375	def __exit__(self, *args):
1376	self.extend(self._stringio.getvalue().splitlines())
1377	sys.stdout = self._stdout
1378
1379	def create_movie(netcdfile, values, variable):
1380	""" Function to create a movie assuming ImageMagick installed!
1381	values= [graph]#[movie_dimension]#[graph_values]
1382	[graph]: which graphic
1383	[movie_dimension]: [dimnmovie]@[dimvmovie]@[moviedelay]@[interval]
1384	[dimnmovie]; name of the dimension from which make the movie
1385	[dimvmovie]; name of the variable with the values of the dimension
1386	[moviedelay]; delay between frames
1387	[interval]; [beg]@[end]@[freq] or -1 (all)
1388	[graph_values]: values to generate the graphic
1389	netcdfile= netCDF file
1390	variable= variable to use (when applicable)
1391	"""
1392	fname = 'create_movie'
1393
1394	if values == 'h':
1395	print fname + '_____________________________________________________________'
1396	print create_movie.__doc__
1397	quit()
1398
1399	graph = values.split('#')[0]
1400	movie_dim = values.split('#')[1]
1401	graph_vals = values.split('#')[2]
1402
1403	ncobj = NetCDFFile(netcdfile, 'r')
1404
1405	# Movie dimension
1406	##
1407	dimnmovie = movie_dim.split('@')[0]
1408	dimvmovie = movie_dim.split('@')[1]
1409	moviedelay = movie_dim.split('@')[2]
1410	moviebeg = int(movie_dim.split('@')[3])
1411
1412	if not drw.searchInlist(ncobj.dimensions.keys(),dimnmovie):
1413	print errormsg
1414	print ' ' + fname + ": file '" + netcdfile + "' has not dimension '" + \
1415	dimnmovie + "' !!!"
1416	quit(-1)
1417
1418	objdmovie = ncobj.dimensions[dimnmovie]
1419	dmovie = len(objdmovie)
1420	if moviebeg != -1:
1421	moviend = int(movie_dim.split('@')[4])
1422	moviefreq = int(movie_dim.split('@')[5])
1423	else:
1424	moviebeg = 0
1425	moviend = dmovie
1426	moviefreq = 1
1427
1428	if dimvmovie == 'WRFTimes':
1429	objvdmovie = ncobj.variables['Times']
1430	vdmovieunits = ''
1431	valsdmovie = []
1432	for it in range(objvdmovie.shape[0]):
1433	valsdmovie.append(drw.datetimeStr_conversion(objvdmovie[it,:], \
1434	'WRFdatetime', 'Y/m/d H-M-S'))
1435	elif dimvmovie == 'CFtime':
1436	objvdmovie = ncobj.variables['time']
1437	vdmovieunits = ''
1438	print objvdmovie.units
1439	valsdmovie0 = drw.netCDFdatetime_realdatetime(objvdmovie.units, 'standard', \
1440	objvdmovie[:])
1441	valsdmovie = []
1442	for it in range(objvdmovie.shape[0]):
1443	valsdmovie.append(drw.datetimeStr_conversion(valsdmovie0[it,:], \
1444	'matYmdHMS', 'Y/m/d H-M-S'))
1445	else:
1446	if not drw.searchInlist(ncobj.variables.keys(),dimvmovie):
1447	print errormsg
1448	print ' ' + fname + ": file '" + netcdfile + "' has not variable '" + \
1449	dimvmovie + "' !!!"
1450	quit(-1)
1451	vdmovieunits = objvdmovie.getncattr('units')
1452	objvdmovie = ncobj.variables[dimvmovie]
1453	if len(objvdmovie.shape) == 1:
1454	vasldmovie = objvdmovie[:]
1455	else:
1456	print errormsg
1457	print ' ' + fname + ': shape', objvdmovie.shape, 'of variable with ' + \
1458	'dimension movie values not ready!!!'
1459	quit(-1)
1460
1461	ncobj.close()
1462	os.system('rm frame_*.png > /dev/null')
1463
1464	# graphic
1465	##
1466	if graph == 'draw_2D_shad':
1467	graphvals = graph_vals.split(':')
1468
1469	for iframe in range(moviebeg,moviend,moviefreq):
1470	iframeS = str(iframe).zfill(4)
1471
1472	drw.percendone((iframe-moviebeg)/moviefreq,(moviend-moviebeg)/moviefreq, \
1473	5, 'frames')
1474	titgraph = dimnmovie + '\|=\|' + str(valsdmovie[iframe]) + '\|' + \
1475	vdmovieunits
1476
1477	graphvals[1] = movievalslice(graphvals[1],dimnmovie,iframe)
1478	graphvals[6] = titgraph
1479	graphvals[7] = 'png'
1480
1481	graphv = drw.numVector_String(graphvals, ":")
1482
1483	with Capturing() as output:
1484	draw_2D_shad(netcdfile, graphv, variable)
1485
1486	os.system('mv 2Dfields_shadow.png frame_' + iframeS + '.png')
1487	else:
1488	print errormsg
1489	print ' ' + fname + ": graphic '" + graph + "' not defined !!!"
1490	quit(-1)
1491
1492	os.system('convert -delay ' + moviedelay + ' -loop 0 frame_*.png create_movie.gif')
1493
1494	print "Succesfuly creation of movie file 'create_movie.gif' !!!"
1495
1496	return
1497
1498	def draw_lines(ncfilens, values, varname):
1499	""" Function to draw different lines at the same time from different files
1500	draw_lines(ncfilens, values, varname):
1501	ncfilens= [filen] ',' separated list of netCDF files
1502	values= [dimvname]:[valuesaxis]:[dimtit]:[leglabels]:[vtit]:[title]:[locleg]:[graphk]
1503	[dimvname]: name of the variable with he values of the common dimension
1504	[valuesaxis]: which axis will be used for the values ('x', or 'y')
1505	[dimtit]: title for the common dimension
1506	[leglabels]: ',' separated list of names for the legend
1507	[vartit]: name of the variable in the graph
1508	[title]: title of the plot ('\|' for spaces)
1509	[locleg]: location of the legend (-1, autmoatic)
1510	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
1511	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
1512	9: 'upper center', 10: 'center'
1513	[graphk]: kind of the graphic
1514	varname= variable to plot
1515	values= 'XLAT:x:latitude:32x32:$wss^{*}$:wss Taylor's turbulence term:pdf'
1516	"""
1517
1518	fname = 'draw_lines'
1519
1520	if values == 'h':
1521	print fname + '_____________________________________________________________'
1522	print draw_lines.__doc__
1523	quit()
1524
1525	expectargs = '[dimvname]:[valuesaxis]:[dimtit]:[leglabels]:[vtit]:[title]:[locleg]:[graphk]'
1526	drw.check_arguments(fname,len(expectargs.split(':')),values,':',expectargs)
1527
1528	ncfiles = ncfilens.split(',')
1529	dimvname = values.split(':')[0]
1530	valuesaxis = values.split(':')[1]
1531	dimtit = values.split(':')[2]
1532	leglabels = values.split(':')[3].replace('_','\_')
1533	vartit = values.split(':')[4]
1534	title = values.split(':')[5].replace('\|',' ')
1535	graphk = values.split(':')[6]
1536
1537	Nfiles = len(ncfiles)
1538
1539	# Getting trajectotries
1540	##
1541
1542	varvalues = []
1543	dimvalues = []
1544
1545	print ' ' + fname
1546	ifn = 0
1547	for ifile in ncfiles:
1548	filen = ifile.split('@')[0]
1549
1550	print ' filen:',filen
1551
1552	if not os.path.isfile(filen):
1553	print errormsg
1554	print ' ' + fname + ": netCDF file '" + filen + "' does not exist !!"
1555	quit(-1)
1556
1557	objfile = NetCDFFile(filen, 'r')
1558
1559	if not objfile.variables.has_key(dimvname):
1560	print errormsg
1561	print ' ' + fname + ": netCDF file '" + filen + \
1562	"' does not have variable '" + dimvname + "' !!"
1563	quit(-1)
1564
1565	if not objfile.variables.has_key(varname):
1566	print errormsg
1567	print ' ' + fname + ": netCDF file '" + filen + \
1568	"' does not have variable '" + varname + "' !!"
1569	quit(-1)
1570
1571	vvobj = objfile.variables[varname]
1572	if len(vvobj.shape) != 1:
1573	print errormsg
1574	print ' ' + fname + ': wrong shape:',vvobj.shape," of variable '" + \
1575	varname + "' !!"
1576	quit(-1)
1577
1578	vdobj = objfile.variables[dimvname]
1579	if len(vdobj.shape) != 1:
1580	print errormsg
1581	print ' ' + fname + ': wrong shape:',vdobj.shape," of variable '" + \
1582	dimvname + "' !!"
1583	quit(-1)
1584
1585	varvalues.append(vvobj[:])
1586	dimvalues.append(vdobj[:])
1587
1588	if ifn == 0:
1589	varunits = vvobj.units
1590
1591	objfile.close()
1592
1593	ifn = ifn + 1
1594
1595
1596	drw.plot_lines(dimvalues, varvalues, valuesaxis, dimtit, leglabels.split(','), \
1597	vartit, varunits, title, graphk)
1598
1599	return
1600
1601	def draw_lines_time(ncfilens, values, varname):
1602	""" Function to draw different lines at the same time from different files with times
1603	draw_lines_time(ncfilens, values, varname):
1604	ncfilens= [filen] ',' separated list of netCDF files
1605	values= [dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];[timevals];[locleg];[graphk]
1606	[dimvname]: name of the variable with he values of the common dimension
1607	[valuesaxis]: which axis will be used for the values ('x', or 'y')
1608	[dimtit]: title for the common dimension
1609	[leglabels]: ',' separated list of names for the legend
1610	[vartit]: name of the variable in the graph
1611	[title]: title of the plot ('\|' for spaces)
1612	[timevals]: [timen]\|[units]\|[kind]\|[tfmt] time labels characteristics
1613	[timen]; name of the time variable
1614	[units]; units string according to CF conventions ([tunits] since
1615	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
1616	[kind]; kind of output
1617	'Nval': according to a given number of values as 'Nval',[Nval]
1618	'exct': according to an exact time unit as 'exct',[tunit];
1619	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1620	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1621	'l': milisecond
1622	[tfmt]; desired format
1623	[locleg]: location of the legend (-1, autmoatic)
1624	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
1625	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
1626	9: 'upper center', 10: 'center'
1627	[graphk]: kind of the graphic
1628	varname= variable to plot
1629	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'
1630	"""
1631
1632	fname = 'draw_lines_time'
1633
1634	if values == 'h':
1635	print fname + '_____________________________________________________________'
1636	print draw_lines_time.__doc__
1637	quit()
1638
1639	expectargs = '[dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];'
1640	expectargs = expectargs + '[timevals];[locleg];[graphk]'
1641	drw.check_arguments(fname,len(expectargs.split(';')),values,';',expectargs)
1642
1643	ncfiles = ncfilens.split(',')
1644	dimvname = values.split(';')[0]
1645	valuesaxis = values.split(';')[1]
1646	dimtit = values.split(';')[2]
1647	leglabels = values.split(';')[3].replace('_','\_')
1648	vartit = values.split(';')[4]
1649	title = values.split(';')[5].replace('\|',' ')
1650	timevals = values.split(';')[6]
1651	locleg = int(values.split(';')[7])
1652	graphk = values.split(';')[8]
1653
1654	Nfiles = len(ncfiles)
1655
1656	# Getting values
1657	##
1658	varvalues = []
1659	dimvalues = []
1660	timvalues = []
1661	timvals0 = timvalues
1662
1663	print ' ' + fname
1664	ifn = 0
1665	mintval = 1.e20
1666	maxtval = -1.e20
1667
1668	for ifile in ncfiles:
1669	filen = ifile.split('@')[0]
1670
1671	print ' filen:',filen
1672
1673	if not os.path.isfile(filen):
1674	print errormsg
1675	print ' ' + fname + ": netCDF file '" + filen + "' does not exist !!"
1676	quit(-1)
1677
1678	objfile = NetCDFFile(filen, 'r')
1679
1680	if not objfile.variables.has_key(dimvname):
1681	print errormsg
1682	print ' ' + fname + ": netCDF file '" + filen + \
1683	"' does not have variable '" + dimvname + "' !!"
1684	quit(-1)
1685
1686	if not objfile.variables.has_key(varname):
1687	print errormsg
1688	print ' ' + fname + ": netCDF file '" + filen + \
1689	"' does not have variable '" + varname + "' !!"
1690	quit(-1)
1691
1692	vvobj = objfile.variables[varname]
1693	if len(vvobj.shape) != 1:
1694	print errormsg
1695	print ' ' + fname + ': wrong shape:',vvobj.shape," of variable '" + \
1696	varname + "' !!"
1697	quit(-1)
1698
1699	vdobj = objfile.variables[dimvname]
1700	if len(vdobj.shape) != 1:
1701	print errormsg
1702	print ' ' + fname + ': wrong shape:',vdobj.shape," of variable '" + \
1703	dimvname + "' !!"
1704	quit(-1)
1705
1706	varvalues.append(vvobj[:])
1707	dimvalues.append(vdobj[:])
1708
1709	mindvals = np.min(vdobj[:])
1710	maxdvals = np.max(vdobj[:])
1711
1712	if mindvals < mintval: mintval = mindvals
1713	if maxdvals > maxtval: maxtval = maxdvals
1714
1715	if ifn == 0:
1716	varunits = drw.units_lunits(vvobj.units)
1717
1718	objfile.close()
1719
1720	ifn = ifn + 1
1721
1722	# Times
1723	timename = timevals.split('\|')[0]
1724	timeunit = timevals.split('\|')[1].replace('!',' ')
1725	timekind = timevals.split('\|')[2]
1726	timefmt = timevals.split('\|')[3]
1727
1728	dtvals = (maxtval - mintval)/5.
1729	tvals = np.arange(mintval, maxtval, dtvals/2.)
1730
1731	timepos, timelabels = drw.CFtimes_plot(tvals, timeunit, timekind, timefmt)
1732
1733	drw.plot_lines_time(dimvalues, varvalues, valuesaxis, dimtit, leglabels.split(','), \
1734	vartit, varunits, timepos, timelabels, title, locleg, graphk)
1735
1736	return
1737
1738	def draw_Neighbourghood_evol(filen, values, variable):
1739	""" Function to draw the temporal evolution of a neighbourghood around a point
1740	draw_Neighbourghood_evol(filen, values, variable)
1741	filen= netCDF file name
1742	values= [gvarname]:[dimsval]:[neigdims]:[Nneig]:[Ncol]:[timetits]:[tkinds]:
1743	[timefmts]:[gtitle]:[shadxtrms]:[cbar]:[gkind]:[ofile]
1744	[dimsval]: [dimn1]\|[val1]\|[dimv1],...,[dimnN]\|[valN]\|[dimvN] dimension names, values to get
1745	(-1, for all; no name/value pair given full length) and variable with values of the dimension
1746	NOTE: when dimsval[X,Y] == neigdims[X,Y], valX,valY --> valX,valY-Nneig/2, valX,valY+Nneig/2
1747	[neigdims]: [dimnX],[dimnY] dimensions mnames along which the neigbourghood should be defined
1748	[Nneig]: Number of grid points of the full side of the box (odd value)
1749	[Ncol]: Number of columns ('auto': square final plot)
1750	[gvarname]: name of the variable to appear in the graph
1751	[timetits]: [titX],[titY] titles of the axes ('\|' for spaces)
1752	[tkinds]: [tkindX]\|[tkindY] kinds of time to appear in the graph
1753	'Nval': according to a given number of values as 'Nval',[Nval]
1754	'exct': according to an exact time unit as 'exct',[tunit];
1755	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1756	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1757	'l': milisecond
1758	[timefmts]: [tfmtX],[tfmtY] format of the time labels
1759	[gtitle]: title of the graphic ('\|' for spaces)
1760	[shadxtrms]: Extremes for the shading
1761	[cbar]: colorbar to use
1762	[gkind]: kind of graphical output
1763	[ofile]: True/False whether the netcdf with data should be created or not
1764	variable= name of the variable
1765	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'
1766	"""
1767
1768	fname = 'draw_Neighbourghood_evol'
1769
1770	if values == 'h':
1771	print fname + '_____________________________________________________________'
1772	print draw_Neighbourghood_evol.__doc__
1773	quit()
1774
1775	expectargs = ['[gvarname]', '[dimsval]', '[neigdims]', '[Nneig]', '[Ncol]', \
1776	'[timetits]', '[tkinds]', '[timefmts]', '[gtitle]', '[shadxtrms]', '[cbar]', \
1777	'[gkind]', '[ofile]']
1778
1779	drw.check_arguments(fname,len(expectargs),values,':',expectargs)
1780
1781	gvarname = values.split(':')[0]
1782	dimsval = values.split(':')[1].split(',')
1783	neigdims = values.split(':')[2].split(',')
1784	Nneig = int(values.split(':')[3])
1785	Ncol0 = values.split(':')[4]
1786	timetits = values.split(':')[5].split(',')
1787	timekinds = values.split(':')[6].split('\|')
1788	timefmts = values.split(':')[7].split(',')
1789	gtitle = values.split(':')[8].replace('\|',' ')
1790	shadxtrms = values.split(':')[9].split(',')
1791	cbar = values.split(':')[10]
1792	gkind = values.split(':')[11]
1793	ofile = values.split(':')[12]
1794
1795	if Ncol0 != 'auto':
1796	Ncol = int(Ncol0)
1797	else:
1798	Ncol = Ncol0
1799
1800	timetits[0] = timetits[0].replace('\|',' ')
1801	timetits[1] = timetits[1].replace('\|',' ')
1802
1803	if np.mod(Nneig,2) == 0:
1804	print errormsg
1805	print ' ' + fname + ": an odd value for 'Nneig':",Nneig,'is required !!!'
1806	quit(-1)
1807
1808	Nneig2 = int(Nneig/2)
1809
1810	# Values to slice the variable
1811	dimvslice = {}
1812	dimvvalues = {}
1813	for dimvs in dimsval:
1814	dimn = dimvs.split('\|')[0]
1815	dimv = int(dimvs.split('\|')[1])
1816	dimnv = dimvs.split('\|')[2]
1817
1818	dimvvalues[dimn] = dimnv
1819	dimvslice[dimn] = dimv
1820
1821	ncobj = NetCDFFile(filen, 'r')
1822
1823	varobj = ncobj.variables[variable]
1824
1825	slicevar = []
1826	newdimn = []
1827	newdimsvar = {}
1828
1829	for dimn in varobj.dimensions:
1830	if not drw.searchInlist(dimvslice.keys(), dimn):
1831	dimsize = len(ncobj.dimensions[dimn])
1832	slicevar.append(slice(0, dimsize+1))
1833	newdimn.append(dimn)
1834	newdimsvar[dimn] = dimseize
1835
1836	for dimslicen in dimvslice.keys():
1837	if dimn == dimslicen:
1838	if dimvslice[dimn] != -1:
1839	if drw.searchInlist(neigdims, dimn):
1840	slicevar.append(slice(dimvslice[dimn]-Nneig2, \
1841	dimvslice[dimn]+Nneig2+1))
1842	newdimn.append(dimn)
1843	newdimsvar[dimn] = Nneig
1844	break
1845	else:
1846	slicevar.append(slice(dimvslice[dimn], dimvslice[dimn]+1))
1847	break
1848	else:
1849	dimsize = len(ncobj.dimensions[dimn])
1850	slicevar.append(slice(0, dimsize+1))
1851	newdimn.append(dimn)
1852	newdimsvar[dimn] = dimsize
1853	break
1854
1855	varv = varobj[tuple(slicevar)]
1856
1857	if len(newdimn) != 3:
1858	print errormsg
1859	print ' ' + fname + ': sliced variable with shape=', varv.shape, \
1860	' must have three dimensions',len(varv.shape),'given !!'
1861	quit(-1)
1862
1863	newdims = []
1864	for nwdims in newdimn:
1865	newdims.append(newdimsvar[nwdims])
1866
1867	# The dimension which is not in the neighbourhood dimensions must be time!
1868	for dim1 in newdimn:
1869	if not drw.searchInlist(neigdims, dim1):
1870	dimt = newdimsvar[dim1]
1871	dimtime = dim1
1872
1873	if Ncol == 'auto':
1874	dimtsqx = int(np.sqrt(dimt)) + 1
1875	dimtsqy = int(np.sqrt(dimt)) + 1
1876	else:
1877	dimtsqx = int(Ncol)
1878	dimtsqy = dimt/dimtsqx + 1
1879
1880	neighbourghood = np.ones((dimtsqyNneig,dimtsqxNneig), dtype=np.float)*fillValue
1881
1882	for it in range(dimt):
1883	ity = int(it/dimtsqx)
1884	itx = it-ity*dimtsqx
1885
1886	itty = (dimtsqy - ity - 1)*Nneig + Nneig2
1887	ittx = itx*Nneig + Nneig2
1888
1889	neighbourghood[itty-Nneig2:itty+Nneig2+1,ittx-Nneig2:ittx+Nneig2+1]= \
1890	varv[it,::-1,:]
1891
1892	variablevals = drw.variables_values(variable)
1893	if drw.searchInlist(varobj.ncattrs(), 'units'):
1894	vunits = varobj.units
1895	else:
1896	vunits = variablevals[5]
1897
1898	# Time values at the X/Y axes
1899	if ncobj.variables[dimvvalues[dimtime]].dtype == '\|S1':
1900	print ' ' + fname + ': WRF time variable!'
1901	refdate = '19491201000000'
1902	tunitsval = 'hours'
1903	dimtvalues = np.zeros((dimt), dtype=np.float)
1904	tvals = ncobj.variables[dimvvalues[dimtime]]
1905	yrref=refdate[0:4]
1906	monref=refdate[4:6]
1907	dayref=refdate[6:8]
1908	horref=refdate[8:10]
1909	minref=refdate[10:12]
1910	secref=refdate[12:14]
1911
1912	refdateS = yrref + '/' + monref + '/' + dayref + '_' + horref + ':' + \
1913	minref + ':' + secref
1914	tunits = tunitsval + ' since ' + refdateS
1915	for it in range(dimt):
1916	wrfdates = drw.datetimeStr_conversion(tvals[it,:],'WRFdatetime', 'matYmdHMS')
1917	dimtvalues[it] = drw.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
1918	else:
1919	dimtvalues = ncobj.variables[dimvvalues[dimtime]][:]
1920	tunits = ncobj.variables[newdimsvar[dimtime]].units
1921
1922	dimxv = dimtvalues[0:dimtsqx]
1923	dimyv = dimtvalues[0:dimt:dimtsqx]
1924
1925	print 'Lluis lens: dimxv, dimyv:',len(dimxv), len(dimyv)
1926
1927	dimn = ['time','time']
1928
1929	if ofile == 'True':
1930	ofilen = 'Neighbourghood_evol.nc'
1931	newnc = NetCDFFile(ofilen, 'w')
1932	# Dimensions
1933	newdim = newnc.createDimension('time',None)
1934	newdim = newnc.createDimension('y',dimtsqy*Nneig)
1935	newdim = newnc.createDimension('x',dimtsqx*Nneig)
1936	# Dimension values
1937	newvar = newnc.createVariable('time','f8',('time'))
1938	newvar[:] = np.arange(dimt)
1939	newattr = drw.basicvardef(newvar, 'time','time',tunits)
1940	# Neighbourhghood variable
1941	newvar = newnc.createVariable(variable + 'neigevol', 'f4', ('y','x'), \
1942	fill_value=fillValue)
1943	newvar[:] = neighbourghood
1944
1945	newnc.sync()
1946	newnc.close()
1947	print fname + ": Successfull generation of file '" + ofilen + "' !!"
1948
1949	# Time ticks
1950	timeposX, timelabelsX = drw.CFtimes_plot(dimxv, tunits, timekinds[0], timefmts[0])
1951	timeposY, timelabelsY = drw.CFtimes_plot(dimyv, tunits, timekinds[1], timefmts[1])
1952
1953	timepos = [timeposX[0:len(timeposX)], timeposY[len(timeposY):0:-1]]
1954	timelabels = [timelabelsX[0:len(timeposX)], timelabelsY[0:len(timeposY)]]
1955
1956	for i in range(2):
1957	if shadxtrms[i][0:1] != 'S':
1958	shadxtrms[i] = np.float(shadxtrms[i])
1959
1960	drw.plot_Neighbourghood_evol(neighbourghood, dimxv, dimyv, gvarname, timetits, \
1961	timepos, timelabels, cbar, Nneig, shadxtrms, vunits, gtitle, gkind, True)
1962
1963	def draw_timeSeries(filen, values, variables):
1964	""" Function to draw a time-series
1965	draw_timeSeries(filen, values, variable):
1966	filen= name of the file
1967	values= [gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[locleg]:[gkind]
1968	[gvarname]: name of the variable to appear in the graph
1969	[timetit]: title of the time axis (assumed x-axis, '\|' for spaces)
1970	[tkind]: kind of time to appear in the graph (assumed x-axis)
1971	'Nval': according to a given number of values as 'Nval',[Nval]
1972	'exct': according to an exact time unit as 'exct',[tunit];
1973	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1974	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1975	'l': milisecond
1976	[timefmt]: format of the time labels
1977	[title]: title of the graphic ('\|' for spaces)
1978	[locleg]: location of the legend (-1, autmoatic)
1979	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
1980	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
1981	9: 'upper center', 10: 'center'
1982	[gkind]: kind of graphical output
1983	variables= [varname],[timename] names of variable and variable with times
1984	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')
1985	"""
1986
1987	fname = 'draw_timeSeries'
1988
1989	if values == 'h':
1990	print fname + '_____________________________________________________________'
1991	print draw_timeSeries.__doc__
1992	quit()
1993
1994	expectargs = ['[gvarname]', '[timetit]', '[tkind]', '[timefmt]', '[title]', \
1995	'[locleg]', '[gkind]']
1996
1997	drw.check_arguments(fname,len(expectargs),values,':',expectargs)
1998
1999	gvarname = values.split(':')[0]
2000	timetit = values.split(':')[1].replace('\|',' ')
2001	tkind = values.split(':')[2]
2002	timefmt = values.split(':')[3]
2003	title = values.split(':')[4].replace('\|',' ')
2004	locleg = int(values.split(':')[5])
2005	gkind = values.split(':')[6]
2006
2007	ncobj = NetCDFFile(filen, 'r')
2008
2009	variable = variables.split(',')[0]
2010	timevar = variables.split(',')[1]
2011
2012	if not ncobj.variables.has_key(variable):
2013	print errormsg
2014	print ' ' + fname + ": file '" + filen + "' does not have variable '" + \
2015	variable + "' !!"
2016	quit(-1)
2017
2018	if not ncobj.variables.has_key(timevar):
2019	print errormsg
2020	print ' ' + fname + ": file '" + filen + "' does not have variable time '" \
2021	+ timevar + "' !!"
2022	quit(-1)
2023
2024	varobj = ncobj.variables[variable]
2025	timeobj = ncobj.variables[timevar]
2026
2027	dimt = len(timeobj[:])
2028	varvals = np.zeros((2,dimt), dtype=np.float)
2029
2030	gunits = varobj.getncattr('units')
2031	tunits = timeobj.getncattr('units')
2032
2033	varvals[0,:], valpot, newgunits, Spot = drw.pot_values(varobj[:].flatten(), gunits)
2034	varvals[1,:] = timeobj[:]
2035
2036	tseriesvals = []
2037	tseriesvals.append(varvals)
2038
2039	drw.plot_TimeSeries(tseriesvals, Spot + drw.units_lunits(gunits), tunits, \
2040	'TimeSeries', gvarname, timetit, tkind, timefmt, title, \
2041	gvarname.replace('_','\_'), locleg, gkind)
2042
2043	return
2044
2045	#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')
2046
2047	def draw_trajectories(trjfilens, values, observations):
2048	""" Function to draw different trajectories at the same time
2049	draw_trajectories(trjfilens, values, observations):
2050	trjfilens= [filen]@[Tint]@[map] ',' separated list of files with trajectories,
2051	time intervals and reference maps (first one will be used to plot)
2052	[filen]: name of the file to use (lines with '#', not readed) as:
2053	[t-step] [x] [y]
2054	[Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2055	[map]: [file]#[lonname]#[latname]
2056	[file]; with the [lon],[lat] matrices
2057	[lonname],[latname]; names of the longitudes and latitudes variables
2058	values=[leglabels]\|[lonlatlims]\|[title]\|[graphk]\|[mapkind]
2059	[leglabels]: ',' separated list of names for the legend
2060	[lonlatlims]: limits of the map [lonmin, latmin, lonmax, latmax] or None
2061	[title]: title of the plot
2062	[graphk]: kind of the graphic
2063	[mapkind]: drawing coastaline ([proj],[res]) or None
2064	[proj]: projection
2065	* 'cyl', cilindric
2066	* 'lcc', lambert conformal
2067	[res]: resolution:
2068	* 'c', crude
2069	* 'l', low
2070	* 'i', intermediate
2071	* 'h', high
2072	* 'f', full
2073	obsevations= [obsfile],[obsname],[Tint],[null]
2074	[obsfile]: name fo the File with the observations as [t-step] [lat] [lon]
2075	[obsname]: name of the observations in the graph
2076	[Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2077	[null]: null value for the observed trajectory
2078	"""
2079
2080	fname = 'draw_trajectories'
2081
2082	if values == 'h':
2083	print fname + '_____________________________________________________________'
2084	print draw_trajectories.__doc__
2085	quit()
2086
2087	trjfiles = trjfilens.split(',')
2088	leglabels = values.split('\|')[0]
2089	lonlatlims = values.split('\|')[1]
2090	title = values.split('\|')[2]
2091	graphk = values.split('\|')[3]
2092	mapkind = values.split('\|')[4]
2093
2094	Nfiles = len(trjfiles)
2095
2096	# Getting trajectotries
2097	##
2098
2099	lontrjvalues = []
2100	lattrjvalues = []
2101
2102	print ' ' + fname
2103	ifn = 0
2104	for ifile in trjfiles:
2105	filen = ifile.split('@')[0]
2106	Tint = ifile.split('@')[1]
2107
2108	print ' trajectory:',filen
2109
2110	if Tint != '-1':
2111	Tbeg = Tint
2112	Tend = ifile.split('@')[2]
2113	mapv = ifile.split('@')[3]
2114	else:
2115	mapv = ifile.split('@')[2]
2116
2117	if not os.path.isfile(filen):
2118	print errormsg
2119	print ' ' + fname + ": trajectory file '" + filen + "' does not exist !!"
2120	quit(-1)
2121
2122	# Charging longitude and latitude values
2123	##
2124	lonvals, latvals = drw.lonlat_values(mapv.split('#')[0], mapv.split('#')[1], \
2125	mapv.split('#')[2])
2126
2127	if ifn == 0: mapref = mapv
2128	ifn = ifn + 1
2129
2130	objfile = open(filen, 'r')
2131	trjtimev = []
2132	trjxv = []
2133	trjyv = []
2134
2135	for line in objfile:
2136	if line[0:1] != '#':
2137	trjtimev.append(int(line.split(' ')[0]))
2138	trjxv.append(int(line.split(' ')[1]))
2139	trjyv.append(int(line.split(' ')[2]))
2140
2141	objfile.close()
2142
2143	if Tint != '-1':
2144	lontrjvalues.append(lonvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
2145	lattrjvalues.append(latvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
2146	else:
2147	lontrjvalues.append(lonvals[trjyv[:], trjxv[:]])
2148	lattrjvalues.append(latvals[trjyv[:], trjxv[:]])
2149
2150	# lonlatlimits
2151	##
2152
2153	if lonlatlims == 'None':
2154	lonlatlimsv = None
2155	else:
2156	lonlatlimsv = np.zeros((4), dtype=np.float)
2157	lonlatlimsv[0] = np.float(lonlatlims.split(',')[0])
2158	lonlatlimsv[1] = np.float(lonlatlims.split(',')[1])
2159	lonlatlimsv[2] = np.float(lonlatlims.split(',')[2])
2160	lonlatlimsv[3] = np.float(lonlatlims.split(',')[3])
2161
2162	# lon/lat objects
2163	##
2164	objnc = NetCDFFile(mapref.split('#')[0])
2165	lonobj = objnc.variables[mapref.split('#')[1]]
2166	latobj = objnc.variables[mapref.split('#')[2]]
2167
2168	# map
2169	##
2170	if mapkind == 'None':
2171	mapkindv = None
2172	else:
2173	mapkindv = mapkind
2174
2175	if observations is None:
2176	obsname = None
2177	else:
2178	obsfile = observations.split(',')[0]
2179	obsname = observations.split(',')[1]
2180	Tint = observations.split(',')[2]
2181	null = np.float(observations.split(',')[3])
2182	print ' observational trajectory:',obsfile
2183
2184	if not os.path.isfile(obsfile):
2185	print errormsg
2186	print ' ' + fname + ": observations trajectory file '" + obsfile + \
2187	"' does not exist !!"
2188	quit(-1)
2189
2190	objfile = open(obsfile, 'r')
2191	obstrjtimev = []
2192	obstrjxv = []
2193	obstrjyv = []
2194
2195	for line in objfile:
2196	if line[0:1] != '#':
2197	lon = np.float(line.split(' ')[2])
2198	lat = np.float(line.split(' ')[1])
2199	if not lon == null and not lat == null:
2200	obstrjtimev.append(int(line.split(' ')[0]))
2201	obstrjxv.append(lon)
2202	obstrjyv.append(lat)
2203	else:
2204	obstrjtimev.append(int(line.split(' ')[0]))
2205	obstrjxv.append(None)
2206	obstrjyv.append(None)
2207
2208	objfile.close()
2209
2210	if Tint != '-1':
2211	Tint = int(observations.split(',')[2].split('@')[0])
2212	Tend = int(observations.split(',')[2].split('@')[1])
2213	lontrjvalues.append(obstrjxv[Tint:Tend+1])
2214	lattrjvalues.append(obstrjyv[Tint:Tend+1])
2215	else:
2216	lontrjvalues.append(obstrjxv[:])
2217	lattrjvalues.append(obstrjyv[:])
2218
2219	drw.plot_Trajectories(lontrjvalues, lattrjvalues, leglabels.split(','), \
2220	lonobj, latobj, lonlatlimsv, title, graphk, mapkindv, obsname)
2221
2222	objnc.close()
2223
2224	return
2225
2226	def draw_vals_trajectories(ncfile, values, variable):
2227	""" Function to draw values from the outputs from 'compute_tevolboxtraj'
2228	draw_vals_trajectories(ncfile, values, variable)
2229	ncfile= [ncfile] ',' list of files to use
2230	values= [statistics]:[Tint]:[labels]@[locleg]:[gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[gkind]
2231	[statistics]: which statistics to use, from: 'center', 'min', 'max', 'mean',
2232	'mean2', 'stdev'
2233	[Tint]: [Tbeg]@[Tend] or None, interval of time to plot or -1 for all the times
2234	[labels]: ',' separated list of labels for the legend
2235	[locleg]: location of the legend (-1, autmoatic)
2236	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2237	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2238	9: 'upper center', 10: 'center'
2239	[gvarname]: name of the variable to appear in the graph
2240	[timetit]: title of the time axis (assumed x-axis, '\|' for spaces)
2241	[tkind]: kind of time to appear in the graph (assumed x-axis)
2242	'Nval': according to a given number of values as 'Nval',[Nval]
2243	'exct': according to an exact time unit as 'exct',[tunit];
2244	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2245	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2246	'l': milisecond
2247	[timefmt]: format of the time labels
2248	[title]: title of the graphic ('\|' for spaces)
2249	[gkind]: kind of graphical output
2250	variable= variable to use
2251	"""
2252
2253	fname = 'draw_vals_trajectories'
2254
2255	if values == 'h':
2256	print fname + '_____________________________________________________________'
2257	print draw_vals_trajectories.__doc__
2258	quit()
2259
2260	sims = ncfile.split(',')
2261
2262	if len(values.split(':')) != 9:
2263	print errormsg
2264	print ' ' + fname + ': wrong number of values!', len(values.split(':')), \
2265	'given 9 needed!!'
2266	print ' ',values.split(':')
2267	quit(-1)
2268
2269	statistics = values.split(':')[0]
2270	Tint = values.split(':')[1]
2271	labels = values.split(':')[2]
2272	gvarname = values.split(':')[3]
2273	timetit = values.split(':')[4].replace('\|',' ')
2274	tkind = values.split(':')[5]
2275	timefmt = values.split(':')[6]
2276	title = values.split(':')[7].replace('\|',' ')
2277	gkind = values.split(':')[8]
2278
2279	leglabels = labels.split('@')[0].split(',')
2280	locleg = int(labels.split('@')[1])
2281
2282	Nsims = len(sims)
2283
2284	if Tint != '-1':
2285	tini = np.float(Tint.split('@')[0])
2286	tend = np.float(Tint.split('@')[1])
2287	else:
2288	tini = -1.
2289	tend = -1.
2290
2291	vartimetrjv = []
2292
2293	print ' ' + fname
2294	for trjfile in sims:
2295	print ' ' + trjfile + ' ...'
2296	if not os.path.isfile(trjfile):
2297	print errormsg
2298	print ' ' + fname + ": trajectory file: '" + trjfile + \
2299	"' does not exist !!"
2300	quit(-1)
2301
2302	trjobj = NetCDFFile(trjfile, 'r')
2303	otim = trjobj.variables['time']
2304	if not trjobj.variables.has_key(statistics + 'box_' + variable):
2305	print errormsg
2306	print ' ' + fname + ": file '" + trjfile + "' does not have variable '"+\
2307	statistics + 'box_' + variable + "' !!"
2308	quit(-1)
2309	ovar = trjobj.variables[statistics + 'box_' + variable]
2310	dimt = otim.shape[0]
2311
2312	if trjfile == sims[0]:
2313	gunits = ovar.getncattr('units')
2314	lname = ovar.getncattr('long_name')
2315	tunits = otim.getncattr('units')
2316
2317	if tini != -1:
2318	tiniid = -1
2319	tendid = -1
2320	for itv in range(dimt):
2321	if otim[itv] <= tini and otim[itv+1] >= tini: tiniid = itv
2322	if otim[itv] <= tend and otim[itv+1] >= tend: tendid = itv
2323
2324	if tiniid == -1 or tendid == -1:
2325	print errormsg
2326	print ' ' + main + ' time interval ', tini,',',tend,' not found: ', \
2327	tendid, ',', tiniid, ' !!'
2328	print ' data interval [',otim[0], otim[dimt-1],']'
2329	quit(-1)
2330	dimt = tendid - tiniid + 1
2331
2332	else:
2333	dimt = otim.shape[0]
2334
2335	valsv = np.zeros((2,dimt), dtype=np.float)
2336
2337	if tini == -1:
2338	valsv[1,:] = otim[:]
2339	valsv[0,:] = ovar[:]
2340	else:
2341	valsv[1,:] = otim[tiniid:tendid+1]
2342	valsv[0,:] = ovar[tiniid:tendid+1]
2343
2344	vartimetrjv.append(valsv)
2345	trjobj.close()
2346
2347	drw.plot_TimeSeries(vartimetrjv, drw.units_lunits(gunits), tunits, \
2348	'val_trajectories_' + statistics, gvarname, timetit, tkind, timefmt, title, \
2349	leglabels, locleg, gkind)
2350
2351	def variable_values(values):
2352	""" Function to give back values for a given variable
2353	values= [varname] name of the variable
2354	"""
2355
2356	fname = 'variable_values'
2357
2358	values = drw.variables_values(values)
2359
2360	print fname,'values:',values
2361	print fname,'variable_name:',values[0]
2362	print fname,'standard_name:',values[1]
2363	print fname,'min,max:',str(values[2]) + ',' + str(values[3])
2364	print fname,'long_name:',values[4]
2365	print fname,'units:',values[5]
2366	print fname,'shad_colors:',values[6]
2367	print fname,'all_values:',drw.numVector_String(values,',')
2368
2369	return
2370
2371	####### ###### ##### #### ### ## #
2372
2373	ngraphics = "'" + drw.numVector_String(namegraphics, "', '") + "'"
2374
2375	### Options
2376	##string_operation="operation to make: " + '\n' + " out, output values -S inidim1,[inidim2,...]:enddim1,[enddim2,...]"
2377	string_operation="""operation to make:
2378	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]
2379	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]
2380	[ckind]:
2381	'cmap': as it gets from colorbar
2382	'fixc,[colname]': fixed color [colname], all stright lines
2383	'fixsignc,[colname]': fixed color [colname], >0 stright, <0 dashed line
2384	"""
2385
2386	#print string_operation
2387
2388	parser = OptionParser()
2389	parser.add_option("-f", "--netCDF_file", dest="ncfile",
2390	help="file to use", metavar="FILE")
2391	parser.add_option("-o", "--operation", type='choice', dest="operation",
2392	choices=namegraphics,
2393	help="operation to make: " + ngraphics, metavar="OPER")
2394	parser.add_option("-S", "--valueS", dest="values",
2395	help="[WHEN APPLICABLE] values to use according to the operation", metavar="VALUES")
2396	parser.add_option("-v", "--variable", dest="varname",
2397	help="[WHEN APPLICABLE] variable to check", metavar="VAR")
2398
2399	(opts, args) = parser.parse_args()
2400
2401	####### #######
2402	## MAIN
2403	#######
2404
2405	# Not checking file operation
2406	Notcheckingfile = ['draw_2D_shad_cont', 'draw_2D_shad_cont_time', \
2407	'draw_2D_shad_line', 'draw_2D_shad_line_time', 'draw_lines', 'draw_lines_time', \
2408	'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories', \
2409	'variable_values']
2410
2411	####### ###### ##### #### ### ## #
2412	errormsg='ERROR -- error -- ERROR -- error'
2413
2414	varn=opts.varname
2415	oper=opts.operation
2416
2417	if opts.ncfile is not None and not os.path.isfile(opts.ncfile) and \
2418	not drw.searchInlist(Notcheckingfile, oper):
2419	print errormsg
2420	print ' ' + main + ': File ' + opts.ncfile + ' does not exist !!'
2421	quit(-1)
2422
2423	if oper == 'create_movie':
2424	create_movie(opts.ncfile, opts.values, opts.varname)
2425	elif oper == 'draw_2D_shad':
2426	draw_2D_shad(opts.ncfile, opts.values, opts.varname)
2427	elif oper == 'draw_2D_shad_time':
2428	draw_2D_shad_time(opts.ncfile, opts.values, opts.varname)
2429	elif oper == 'draw_2D_shad_cont':
2430	draw_2D_shad_cont(opts.ncfile, opts.values, opts.varname)
2431	elif oper == 'draw_2D_shad_cont_time':
2432	draw_2D_shad_cont_time(opts.ncfile, opts.values, opts.varname)
2433	elif oper == 'draw_2D_shad_line':
2434	draw_2D_shad_line(opts.ncfile, opts.values, opts.varname)
2435	elif oper == 'draw_2D_shad_line_time':
2436	draw_2D_shad_line_time(opts.ncfile, opts.values, opts.varname)
2437	elif oper == 'draw_Neighbourghood_evol':
2438	draw_Neighbourghood_evol(opts.ncfile, opts.values, opts.varname)
2439	elif oper == 'draw_lines':
2440	draw_lines(opts.ncfile, opts.values, opts.varname)
2441	elif oper == 'draw_lines_time':
2442	draw_lines_time(opts.ncfile, opts.values, opts.varname)
2443	elif oper == 'draw_timeSeries':
2444	draw_timeSeries(opts.ncfile, opts.values, opts.varname)
2445	elif oper == 'draw_topo_geogrid':
2446	draw_topo_geogrid(opts.ncfile, opts.values)
2447	elif oper == 'draw_topo_geogrid_boxes':
2448	draw_topo_geogrid_boxes(opts.ncfile, opts.values)
2449	elif oper == 'draw_trajectories':
2450	draw_trajectories(opts.ncfile, opts.values, opts.varname)
2451	elif oper == 'draw_vals_trajectories':
2452	draw_vals_trajectories(opts.ncfile, opts.values, opts.varname)
2453	elif oper == 'list_graphics':
2454	# From: http://www.diveintopython.net/power_of_introspection/all_together.html
2455	import drawing as myself
2456	object = myself
2457	for opern in namegraphics:
2458	if opern != 'list_graphics':
2459	print opern + '_______ ______ _____ ____ ___ __ _'
2460	print getattr(object, opern).__doc__
2461	elif oper == 'variable_values':
2462	variable_values(opts.values)
2463	else:
2464	print errormsg
2465	print ' ' + main + ": the graphic '" + oper + "' is not ready !!"
2466	print errormsg
2467	quit()

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

Download in other formats: