Context Navigation

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

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

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

Download in other formats: