Context Navigation

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

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

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

Download in other formats: