Context Navigation

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

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

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

Download in other formats: