Context Navigation

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

Last change on this file since 664 was 664, checked in by lfita, 9 years ago
Changing to float for the min and maximum values for 'topo' plots
File size: 117.3 KB

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

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

Download in other formats: