Context Navigation

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

Last change on this file since 680 was 680, checked in by lfita, 9 years ago
Fixing `draw_lines_time', right no-line when including points
File size: 121.8 KB

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

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

Download in other formats: