Context Navigation

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

Last change on this file since 762 was 705, checked in by lfita, 9 years ago
Adding `draw_river_desc' to plot rivers' description after` river_desc.nc' from ORCDHIEE
File size: 142.1 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	## e.g. # drawing.py -o draw_vectors -f wrfout_d01_2001-11-11_00:00:00 -S 'T\|Time\|Times\|2,Y\|south_north\|XLAT\|-1,X\|west_east\|XLONG\|-1:3@3,wind@rainbow,9:10m wind,ms-1:cyl,l:WRF 10 m winds:pdf:winds' -v U10,V10
29	## e.g. # drawing.py -o draw_basins -f routing.py -S 'Y\|y\|nav_lat\|-1,X\|x\|nav_lon\|-1:1@1,rainbow,9:basins,-:cyl,l:ORCDHIEE river-basins:pdf:basins_named' -v nav_lon,nav_lat,trip,basins
30	## e.g. # drawing.py -o draw_river_desc -f diver_desc.nc -S 'Y\|lat\|lat\|-1,X\|lon\|lon\|-1:red,green:Blues:cyl,l:ORCDHIEE rivers:pdf:0:or_rivers -v Amazon
31
32	main = 'drawing.py'
33
34	errormsg = 'ERROR -- error -- ERROR -- error'
35	warnmsg = 'WARNING -- waring -- WARNING -- warning'
36	fillValue=1.e20
37
38	namegraphics = ['create_movie', 'draw_2D_shad', 'draw_2D_shad_time', \
39	'draw_2D_shad_cont', 'draw_2D_shad_cont_time', 'draw_2D_shad_line', \
40	'draw_2D_shad_line_time', 'draw_barbs', 'draw_basins', \
41	'draw_lines', 'draw_lines_time', 'draw_Neighbourghood_evol', \
42	'draw_points', 'draw_points_lonlat', \
43	'draw_ptZvals', 'draw_river_desc', 'draw_timeSeries', \
44	'draw_topo_geogrid', \
45	'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories', \
46	'draw_vectors', 'list_graphics', 'variable_values']
47
48	def draw_2D_shad(ncfile, values, varn):
49	""" plotting a fields with shading
50	draw_2D_shad(ncfile, values, varn)
51	ncfile= file to use
52	values=[vnamefs]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
53	[kindfig]:[reverse]:[mapv]:[close]
54	[vnamefs]: Name in the figure of the variable to be shaded
55	[dimvals]: ',' list of [dimname]\|[value] telling at which dimension of the
56	variable a given value is required (-1, all the length)
57	[dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
58	[colorbar]: name of the color bar
59	[smin/axv]: minimum and maximum value for the shading or:
60	'Srange': for full range
61	'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
62	'Saroundminmax@val': for minval,maxval
63	'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
64	percentile_(100-val)-median)
65	'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min@val,max@val-mean)
66	'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
67	'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
68	percentile_(100-val)-median)
69	[figt]: title of the figure ('\|' for spaces)
70	[kindfig]: kind of figure
71	[reverse]: Transformation of the values
72	* 'transpose': reverse the axes (x-->y, y-->x)
73	* 'flip'@[x/y]: flip the axis x or y
74	[mapv]: map characteristics: [proj],[res]
75	see full documentation: http://matplotlib.org/basemap/
76	[proj]: projection
77	* 'cyl', cilindric
78	* 'lcc', lamvbert conformal
79	[res]: resolution:
80	* 'c', crude
81	* 'l', low
82	* 'i', intermediate
83	* 'h', high
84	* 'f', full
85	valules= 'rh:z\|-1,x\|-1:z\|-1,x\|-1:lat:pressure:BuPu:0.,100.:rh:pdf:flip@y:None'
86	varn= [varsn] name of the variable to plot with shading
87	"""
88
89	fname = 'draw_2D_shad'
90	if values == 'h':
91	print fname + '_____________________________________________________________'
92	print draw_2D_shad.__doc__
93	quit()
94
95	expectargs = '[vnamefs]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]'+\
96	':[figt]:[kindfig]:[reverse]:[mapv]:[close]'
97
98	drw.check_arguments(fname,values,expectargs,':')
99
100	vnamesfig = values.split(':')[0]
101	dimvals= values.split(':')[1].replace('\|',':')
102	vdimxn = values.split(':')[2]
103	vdimyn = values.split(':')[3]
104	colbarn = values.split(':')[4]
105	shadminmax = values.split(':')[5]
106	figtitle = values.split(':')[6].replace('\|',' ')
107	figkind = values.split(':')[7]
108	revals = values.split(':')[8]
109	mapvalue = values.split(':')[9]
110	# varn = values.split(':')[10]
111
112	ncfiles = ncfile
113
114	if not os.path.isfile(ncfiles):
115	print errormsg
116	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
117	quit(-1)
118
119	objsf = NetCDFFile(ncfiles, 'r')
120
121	varns = varn.split(',')[0]
122
123	if not objsf.variables.has_key(varns):
124	print errormsg
125	print ' ' + fname + ': shading file "' + ncfiles + \
126	'" does not have variable "' + varns + '" !!'
127	quit(-1)
128
129	# Variables' values
130	objvars = objsf.variables[varns]
131
132	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
133
134	# Dimensions names
135	## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
136	## dimnamesv = []
137	## for idd in range(len(objvars.dimensions)):
138	## cutdim = False
139	## for idc in range(len(dimvals.split(','))):
140	## dimcutn = dimvals.split(',')[idc].split(':')[0]
141	## print objvars.dimensions[idd], dimcutn
142	## if objvars.dimensions[idd] == dimcutn:
143	## cutdim = True
144	## break
145	## if not cutdim: dimnamesv.append(objvars.dimensions[idd])
146	dimnamesv = [vdimyn, vdimxn]
147
148	if drw.searchInlist(objvars.ncattrs(),'units'):
149	varunits = objvars.getncattr('units')
150	else:
151	print warnmsg
152	print ' ' + fname + ": variable '" + varn + "' without units!!"
153	varunits = '-'
154
155	if not objsf.variables.has_key(vdimxn):
156	print errormsg
157	print ' ' + fname + ': shading file "' + ncfiles + \
158	'" does not have dimension variable "' + vdimxn + '" !!'
159	quit(-1)
160	if not objsf.variables.has_key(vdimyn):
161	print errormsg
162	print ' ' + fname + ': shading file "' + ncfiles + \
163	'" does not have dimension variable "' + vdimyn + '" !!'
164	quit(-1)
165
166	objdimx = objsf.variables[vdimxn]
167	objdimy = objsf.variables[vdimyn]
168	if drw.searchInlist(objdimx.ncattrs(),'units'):
169	odimxu = objdimx.getncattr('units')
170	else:
171	print warnmsg
172	print ' ' + fname + ": variable dimension '" + vdimxn + "' without units!!"
173	odimxu = '-'
174
175	if drw.searchInlist(objdimy.ncattrs(),'units'):
176	odimyu = objdimy.getncattr('units')
177	else:
178	print warnmsg
179	print ' ' + fname + ": variable dimension '" + vdimyn + "' without units!!"
180	odimyu = '-'
181
182	odimxv, odimyv = drw.dxdy_lonlatDIMS(objdimx[:], objdimy[:], objdimx.dimensions, \
183	objdimy.dimensions, dimvals.replace(':','\|').split(','))
184
185
186	# if len(objdimx.shape) <= 2:
187	## odimxv = objdimx[valshad.shape]
188	## odimyv = objdimy[valshad.shape]
189	# odimxv = objdimx[:]
190	# odimyv = objdimy[:]
191
192	# elif len(objdimx.shape) == 3:
193	## dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])]
194	## odimxv = objdimx[tuple(dimcut)]
195	## odimyv = objdimy[tuple(dimcut)]
196	# odimxv = objdimx[0,:]
197	# odimyv = objdimy[0,:]
198	# else:
199	# print errormsg
200	# print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
201	# ' not ready!!'
202	# quit(-1)
203
204	shading_nx = []
205	if shadminmax.split(',')[0][0:1] != 'S':
206	shading_nx.append(np.float(shadminmax.split(',')[0]))
207	else:
208	shading_nx.append(shadminmax.split(',')[0])
209
210	if shadminmax.split(',')[1][0:1] != 'S':
211	shading_nx.append(np.float(shadminmax.split(',')[1]))
212	else:
213	shading_nx.append(shadminmax.split(',')[1])
214
215	if mapvalue == 'None': mapvalue = None
216
217	drw.plot_2D_shadow(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu, dimnamesv,\
218	colbarn, shading_nx, varunits, figtitle, figkind, revals, mapvalue, True)
219
220	return
221
222	def draw_2D_shad_time(ncfile, values, varn):
223	""" plotting a fields with shading with time values
224	draw_2D_shad(ncfile, values, varn)
225	ncfile= file to use
226	values=[vnamefs]~[dimvals]~[dimxvn]~[dimyvn]~[colorbar]~[sminv],[smaxv]~[figt]~
227	[kindfig]~[reverse]~[timevals]~[close]
228	[vnamefs]: Name in the figure of the variable to be shaded
229	[dimvals]: ',' list of [dimname]\|[value] telling at which dimension of the
230	variable a given value is required (-1, all the length, [beg]@[end] for an interval)
231	[dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
232	[colorbar]: name of the color bar
233	[smin/axv]: minimum and maximum value for the shading or:
234	'Srange': for full range
235	'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
236	'Saroundminmax@val': for minval,maxval
237	'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
238	percentile_(100-val)-median)
239	'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min@val,max@val-mean)
240	'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
241	'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
242	percentile_(100-val)-median)
243	[figt]: title of the figure ('\|' for spaces)
244	[kindfig]: kind of figure
245	[reverse]: Transformation of the values
246	* 'transpose': reverse the axes (x-->y, y-->x)
247	* 'flip'@[x/y]: flip the axis x or y
248	[timevals]: [timen]\|[units]\|[kind]\|[tfmt]\|[label]\|[timeaxis] time labels characteristics
249	[timen]; name of the time variable
250	[units]; units string according to CF conventions ([tunits] since
251	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
252	[kind]; kind of output
253	'Nval': according to a given number of values as 'Nval',[Nval]
254	'exct': according to an exact time unit as 'exct',[tunit];
255	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
256	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
257	'l': milisecond
258	[tfmt]; desired format
259	[label]; label at the graph ('!' for spaces)
260	[close]: should figure be closed (finished)
261	values='dtcon~Time\|-1,bottom_top\|-1~presmean~time~seismic~-3.e-6,3.e-6~monthly\|'
262	'dtcon~pdf~transpose~time\|hours!since!1949-12-01\|exct,2,d\|%d\|date!([DD])~True
263	varn= [varsn] name of the variable to plot with shading
264	"""
265	fname = 'draw_2D_shad_time'
266
267	if values == 'h':
268	print fname + '_____________________________________________________________'
269	print draw_2D_shad_time.__doc__
270	quit()
271
272	farguments = '[vnamefs]~[dimvals]~[dimxvn]~[dimyvn]~[colorbar]~[sminv],[smaxv]~'+\
273	'[figt]~[kindfig]~[reverse]~[timevals]~[close]'
274	drw.check_arguments(fname,values,farguments,'~')
275
276	vnamesfig = values.split('~')[0]
277	dimvals= values.split('~')[1].replace('\|',':')
278	vdimxn = values.split('~')[2]
279	vdimyn = values.split('~')[3]
280	colbarn = values.split('~')[4]
281	shadminmax = values.split('~')[5]
282	figtitle = values.split('~')[6].replace('\|',' ')
283	figkind = values.split('~')[7]
284	revals = values.split('~')[8]
285	timevals = values.split('~')[9]
286	close = values.split('~')[10]
287
288	ncfiles = ncfile
289
290	if not os.path.isfile(ncfiles):
291	print errormsg
292	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
293	quit(-1)
294
295	objsf = NetCDFFile(ncfiles, 'r')
296
297	varns = varn.split(',')[0]
298
299	if not objsf.variables.has_key(varns):
300	print errormsg
301	print ' ' + fname + ': shading file "' + ncfiles + \
302	'" does not have variable "' + varns + '" !!'
303	quit(-1)
304
305	# Variables' values
306	objvars = objsf.variables[varns]
307
308	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
309
310	dimnamesv = [vdimyn, vdimxn]
311
312	varunits = objvars.getncattr('units')
313
314	if not objsf.variables.has_key(vdimxn):
315	print errormsg
316	print ' ' + fname + ': shading file "' + ncfiles + \
317	'" does not have dimension variable "' + vdimxn + '" !!'
318	quit(-1)
319	if not objsf.variables.has_key(vdimyn):
320	print errormsg
321	print ' ' + fname + ': shading file "' + ncfiles + \
322	'" does not have dimensino variable "' + vdimyn + '" !!'
323	quit(-1)
324
325	objdimx = objsf.variables[vdimxn]
326	objdimy = objsf.variables[vdimyn]
327	odimxu = objdimx.getncattr('units')
328	odimyu = objdimy.getncattr('units')
329
330	if len(objdimx.shape) <= 2:
331	odimxv0 = objdimx[:]
332	odimyv0 = objdimy[:]
333
334	elif len(objdimx.shape) == 3:
335	odimxv0 = objdimx[0,:]
336	odimyv0 = objdimy[0,:]
337	else:
338	print errormsg
339	print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
340	' not ready!!'
341	quit(-1)
342
343	timename = timevals.split('\|')[0]
344	timeunit = timevals.split('\|')[1].replace('!',' ')
345	timekind = timevals.split('\|')[2]
346	timefmt = timevals.split('\|')[3]
347	timelabel = timevals.split('\|')[4].replace('!',' ')
348
349	# Dimensional values
350	odxv, dimsdxv = drw.slice_variable(objsf.variables[vdimxn], dimvals.replace(',','\|'))
351	odyv, dimsdyv = drw.slice_variable(objsf.variables[vdimyn], dimvals.replace(',','\|'))
352
353	if vdimxn == timename:
354	odimxv = objsf.variables[vdimxn][:]
355	odimxu = timelabel
356	timeaxis = 'x'
357	odimyv = objsf.variables[vdimyn]
358	odimyu = odimyv.getncattr('units')
359	timepos, timelabels = drw.CFtimes_plot(odxv, timeunit, timekind, timefmt)
360	elif vdimyn == timename:
361	odimyv = objsf.variables[vdimyn]
362	odimyu = timelabel
363	timeaxis = 'y'
364	odimxv = objsf.variables[vdimxn]
365	odimxu = odimxv.getncattr('units')
366	timepos, timelabels = drw.CFtimes_plot(odyv, timeunit, timekind, timefmt)
367	else:
368	print errormsg
369	print ' ' + fname + ": time variable '" + timename + "' not found!!"
370	quit(-1)
371
372	shading_nx = []
373	if shadminmax.split(',')[0][0:1] != 'S':
374	shading_nx.append(np.float(shadminmax.split(',')[0]))
375	else:
376	shading_nx.append(shadminmax.split(',')[0])
377
378	if shadminmax.split(',')[1][0:1] != 'S':
379	shading_nx.append(np.float(shadminmax.split(',')[1]))
380	else:
381	shading_nx.append(shadminmax.split(',')[1])
382
383	closeval = drw.Str_Bool(close)
384
385	drw.plot_2D_shadow_time(valshad, vnamesfig, odxv, odyv, odimxu, odimyu, \
386	dimnamesv, colbarn, shading_nx, varunits, figtitle, figkind, revals, timeaxis, \
387	timepos, timelabels, closeval)
388
389	return
390
391	def draw_2D_shad_cont(ncfile, values, varn):
392	""" plotting two fields, one with shading and the other with contour lines
393	draw_2D_shad_cont(ncfile, values, varn)
394	ncfile= [ncfilevars],[ncfilevarc] files to use (one value, same file)
395	values=[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:[figt]:[kindfig]:[reverse]:[mapv]
396	[vnamefs],[vnamefc]: Name in the figure of the shaded and the contour variables
397	[dimvals/c]: list of [dimname]\|[value] telling at which dimension of the
398	variable a given value is required (no dimension name, all the length)
399	[dimx/yvn]: names of the variables with the values of the dimensions for the plot
400	[colorbar]: name of the color bar
401	[ckind]: kind of contours
402	'cmap': as it gets from colorbar
403	'fixc,[colname]': fixed color [colname], all stright lines
404	'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed line
405	[clabfmt]: format of the labels in the contour (None, also possible)
406	[smin/axv]: minimum and maximum value for the shading
407	[sminc]:[smaxv]:[Nlev]: minimum, maximum and number of values for the contour
408	[figt]: title of the figure ('\|' for spaces)
409	[kindfig]: kind of figure
410	[reverse]: does the values be transposed? 'True/False',
411	[mapv]: map characteristics: [proj],[res]
412	see full documentation: http://matplotlib.org/basemap/
413	[proj]: projection
414	* 'cyl', cilindric
415	* 'lcc', lamvbert conformal
416	[res]: resolution:
417	* 'c', crude
418	* 'l', low
419	* 'i', intermediate
420	* 'h', high
421	* 'f', full
422	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'
423	varn= [varsn],[varcn] name of the variable to plot with shading variable with contour
424	"""
425
426	fname = 'draw_2D_shad_cont'
427	if values == 'h':
428	print fname + '_____________________________________________________________'
429	print draw_2D_shad_cont.__doc__
430	quit()
431
432	expectargs = '[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:'
433	expectargs = expectargs + '[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],'
434	expectargs = expectargs + '[Nlev]:[figt]:[kindfig]:[reverse]:[mapv]'
435
436	drw.check_arguments(fname,values,expectargs,':')
437
438	vnamesfig = values.split(':')[0].split(',')
439	dimvals= values.split(':')[1].replace('\|',':')
440	dimvalc= values.split(':')[2].replace('\|',':')
441	vdimxn = values.split(':')[3]
442	vdimyn = values.split(':')[4]
443	colbarn = values.split(':')[5]
444	countkind = values.split(':')[6]
445	countlabelfmt = values.split(':')[7]
446	shadminmax = values.split(':')[8]
447	contlevels = values.split(':')[9]
448	figtitle = values.split(':')[10].replace('\|',' ')
449	figkind = values.split(':')[11]
450	revals = values.split(':')[12]
451	mapvalue = values.split(':')[13]
452
453	if2filenames = ncfile.find(',')
454
455	if if2filenames != -1:
456	ncfiles = ncfile.split(',')[0]
457	ncfilec = ncfile.split(',')[1]
458	else:
459	ncfiles = ncfile
460	ncfilec = ncfile
461
462	if not os.path.isfile(ncfiles):
463	print errormsg
464	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
465	quit(-1)
466
467	if not os.path.isfile(ncfilec):
468	print errormsg
469	print ' ' + fname + ': contour file "' + ncfilec + '" does not exist !!'
470	quit(-1)
471
472	objsf = NetCDFFile(ncfiles, 'r')
473	objcf = NetCDFFile(ncfilec, 'r')
474
475	varns = varn.split(',')[0]
476	varnc = varn.split(',')[1]
477
478	if not objsf.variables.has_key(varns):
479	print errormsg
480	print ' ' + fname + ': shading file "' + ncfiles + \
481	'" does not have variable "' + varns + '" !!'
482	quit(-1)
483
484	if not objcf.variables.has_key(varnc):
485	print errormsg
486	print ' ' + fname + ': contour file "' + ncfilec + \
487	'" does not have variable "' + varnc + '" !!'
488	quit(-1)
489
490	# Variables' values
491	objvars = objsf.variables[varns]
492	objvarc = objcf.variables[varnc]
493
494	if len(objvars.shape) != len(objvarc.shape):
495	print errormsg
496	print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \
497	objvars.shape, 'different than contour variable "' + varnc + '": ', \
498	objvarc.shape,' !!!'
499	quit(-1)
500
501	for idim in range(len(objvars.shape)):
502	if objvars.shape[idim] != objvarc.shape[idim]:
503	print errormsg
504	print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \
505	objvars.shape, 'different than contour variable "' + varnc + '": ', \
506	objvarc.shape,' !!!'
507	quit(-1)
508
509	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
510	valcont, dimscont = drw.slice_variable(objvarc, dimvalc.replace(',','\|'))
511
512	# Dimensions names
513	## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
514	## dimnamesv = []
515	## for idd in range(len(objvars.dimensions)):
516	## cutdim = False
517	## for idc in range(len(dimvals.split(','))):
518	## dimcutn = dimvals.split(',')[idc].split(':')[0]
519	## print objvars.dimensions[idd], dimcutn
520	## if objvars.dimensions[idd] == dimcutn:
521	## cutdim = True
522	## break
523	## if not cutdim: dimnamesv.append(objvars.dimensions[idd])
524	dimnamesv = [vdimyn, vdimxn]
525
526	varunits = []
527	varunits.append(objvars.getncattr('units'))
528	varunits.append(objvarc.getncattr('units'))
529
530	if not objsf.variables.has_key(vdimxn):
531	print errormsg
532	print ' ' + fname + ': shading file "' + ncfiles + \
533	'" does not have dimension variable "' + vdimxn + '" !!'
534	quit(-1)
535	if not objsf.variables.has_key(vdimyn):
536	print errormsg
537	print ' ' + fname + ': shading file "' + ncfiles + \
538	'" does not have dimensino variable "' + vdimyn + '" !!'
539	quit(-1)
540
541	objdimx = objsf.variables[vdimxn]
542	objdimy = objsf.variables[vdimyn]
543	odimxu = objdimx.getncattr('units')
544	odimyu = objdimy.getncattr('units')
545
546	# Getting only that dimensions with coincident names
547	odimxv, odimyv = drw.dxdy_lonlatDIMS(objdimx[:], objdimy[:], objdimx.dimensions, \
548	objdimy.dimensions, dimvals.replace(':','\|').split(','))
549
550	# dimnvx = objdimx.dimensions
551	# cutslice = []
552	# for idimn in objdimx.dimensions:
553	# found = False
554	# for dimsn in dimsshad:
555	# if idimn == dimsn:
556	# cutslice.append(slice(0,len(objsf.dimensions[idimn])))
557	# found = True
558	# if not found: cutslice.append(0)
559	#
560	# odimxv = objdimx[tuple(cutslice)]
561	#
562	# dimnvy = objdimy.dimensions
563	# cutslice = []
564	# for idimn in objdimy.dimensions:
565	# found = False
566	# for dimsn in dimsshad:
567	# if idimn == dimsn:
568	# cutslice.append(slice(0,len(objsf.dimensions[idimn])))
569	# found = True
570	# if not found: cutslice.append(0)
571	#
572	# odimyv = objdimy[tuple(cutslice)]
573
574	# if len(objdimx.shape) <= 2:
575	# odimxv = objdimx[:]
576	# odimyv = objdimy[:]
577	# elif len(objdimx.shape) == 3:
578	# odimxv = objdimx[0,:]
579	# odimyv = objdimy[0,:]
580	# else:
581	# print errormsg
582	# print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
583	# ' not ready!!'
584	# quit(-1)
585
586	if countlabelfmt == 'None':
587	countlfmt = None
588	else:
589	countlfmt = countlabelfmt
590
591	shading_nx = np.zeros((2), dtype=np.float)
592	shading_nx[0] = np.float(shadminmax.split(',')[0])
593	shading_nx[1] = np.float(shadminmax.split(',')[1])
594
595	clevmin = np.float(contlevels.split(',')[0])
596	clevmax = np.float(contlevels.split(',')[1])
597	Nclevels = int(contlevels.split(',')[2])
598
599	levels_cont = drw.pretty_int(clevmin, clevmax, Nclevels)
600
601	if len(levels_cont) <= 1:
602	print warnmsg
603	print ' ' + fname + ': wrong contour levels:', levels_cont,' !!'
604	del(levels_cont)
605	levels_cont = np.zeros((Nclevels), dtype=np.float)
606	levels_cont = np.arange(7)*(clevmax - clevmin)/(Nclevels-1)
607	print ' generating default ones: ',levels_cont
608
609	if mapvalue == 'None': mapvalue = None
610
611	drw.plot_2D_shadow_contour(valshad, valcont, vnamesfig, odimxv, odimyv, odimxu, \
612	odimyu, dimnamesv, colbarn, countkind, countlfmt, shading_nx, levels_cont, \
613	varunits, figtitle, figkind, revals, mapvalue)
614
615	return
616
617	def draw_2D_shad_cont_time(ncfile, values, varn):
618	""" plotting two fields, one with shading and the other with contour lines being
619	one of the dimensions of time characteristics
620	draw_2D_shad_cont(ncfile, values, varn)
621	ncfile= [ncfilevars],[ncfilevarc] files to use (one value, same file)
622	values=[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:[figt]:[kindfig]:[reverse]:[timevals]:[mapv]
623	[vnamefs],[vnamefc]: Name in the figure of the shaded and the contour variables
624	[dimvals/c]: list of [dimname]\|[value] telling at which dimension of the
625	variable a given value is required (no dimension name, all the length)
626	[dimx/yvn]: ',' list with the name of the variables with the values of the dimensions
627	[colorbar]: name of the color bar
628	[ckind]: kind of contours
629	'cmap': as it gets from colorbar
630	'fixc,[colname]': fixed color [colname], all stright lines
631	'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed line
632	[clabfmt]: format of the labels in the contour (None, also possible)
633	[smin/axv]: minimum and maximum value for the shading
634	[sminc]:[smaxv]:[Nlev]: minimum, maximum and number of values for the contour
635	[figt]: title of the figure ('\|' for spaces)
636	[kindfig]: kind of figure
637	[reverse]: modification to the dimensions:
638	'transposed': transpose matrices
639	'flip',[x/y]: flip only the dimension [x] or [y]
640	[timevals]: [timen]\|[units]\|[kind]\|[tfmt]\|[label] time labels characteristics
641	[timen]; name of the time variable
642	[units]; units string according to CF conventions ([tunits] since
643	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
644	[kind]; kind of output
645	'Nval': according to a given number of values as 'Nval',[Nval]
646	'exct': according to an exact time unit as 'exct',[tunit];
647	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
648	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
649	'l': milisecond
650	[tfmt]; desired format
651	[label]; label at the graph ('!' for spaces)
652	[mapv]: map characteristics: [proj],[res]
653	see full documentation: http://matplotlib.org/basemap/
654	[proj]: projection
655	* 'cyl', cilindric
656	* 'lcc', lamvbert conformal
657	[res]: resolution:
658	* 'c', crude
659	* 'l', low
660	* 'i', intermediate
661	* 'h', high
662	* 'f', full
663	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'
664	varn= [varsn],[varcn] name of the variable to plot with shading variable with contour
665	"""
666
667	fname = 'draw_2D_shad_cont_time'
668	if values == 'h':
669	print fname + '_____________________________________________________________'
670	print draw_2D_shad_cont_time.__doc__
671	quit()
672
673	expectargs = '[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:' + \
674	'[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:[figt]:[kindfig]:' + \
675	'[reverse]:[timevals]:[mapv]'
676
677	drw.check_arguments(fname,values,expectargs,':')
678
679	vnamesfig = values.split(':')[0].split(',')
680	dimvals= values.split(':')[1].replace('\|',':')
681	dimvalc= values.split(':')[2].replace('\|',':')
682	vdimxn = values.split(':')[3]
683	vdimyn = values.split(':')[4]
684	colbarn = values.split(':')[5]
685	countkind = values.split(':')[6]
686	countlabelfmt = values.split(':')[7]
687	shadminmax = values.split(':')[8]
688	contlevels = values.split(':')[9]
689	figtitle = values.split(':')[10].replace('\|',' ')
690	figkind = values.split(':')[11]
691	revals = values.split(':')[12]
692	timevals = values.split(':')[13]
693	mapvalue = values.split(':')[14]
694
695	if2filenames = ncfile.find(',')
696
697	if if2filenames != -1:
698	ncfiles = ncfile.split(',')[0]
699	ncfilec = ncfile.split(',')[1]
700	else:
701	ncfiles = ncfile
702	ncfilec = ncfile
703
704	if not os.path.isfile(ncfiles):
705	print errormsg
706	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
707	quit(-1)
708
709	if not os.path.isfile(ncfilec):
710	print errormsg
711	print ' ' + fname + ': contour file "' + ncfilec + '" does not exist !!'
712	quit(-1)
713
714	objsf = NetCDFFile(ncfiles, 'r')
715	objcf = NetCDFFile(ncfilec, 'r')
716
717	varns = varn.split(',')[0]
718	varnc = varn.split(',')[1]
719
720	if not objsf.variables.has_key(varns):
721	print errormsg
722	print ' ' + fname + ': shading file "' + ncfiles + \
723	'" does not have variable "' + varns + '" !!'
724	quit(-1)
725
726	if not objcf.variables.has_key(varnc):
727	print errormsg
728	print ' ' + fname + ': contour file "' + ncfilec + \
729	'" does not have variable "' + varnc + '" !!'
730	quit(-1)
731
732	# Variables' values
733	objvars = objsf.variables[varns]
734	objvarc = objcf.variables[varnc]
735
736	if len(objvars.shape) != len(objvarc.shape):
737	print errormsg
738	print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \
739	objvars.shape, 'different than contour variable "' + varnc + '": ', \
740	objvarc.shape,' !!!'
741	quit(-1)
742
743	for idim in range(len(objvars.shape)):
744	if objvars.shape[idim] != objvarc.shape[idim]:
745	print errormsg
746	print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \
747	objvars.shape, 'different than contour variable "' + varnc + '": ', \
748	objvarc.shape,' !!!'
749	quit(-1)
750
751	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
752	valcont, dimscont = drw.slice_variable(objvarc, dimvalc.replace(',','\|'))
753
754	# Dimensions names
755	## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
756	## dimnamesv = []
757	## for idd in range(len(objvars.dimensions)):
758	## cutdim = False
759	## for idc in range(len(dimvals.split(','))):
760	## dimcutn = dimvals.split(',')[idc].split(':')[0]
761	## print objvars.dimensions[idd], dimcutn
762	## if objvars.dimensions[idd] == dimcutn:
763	## cutdim = True
764	## break
765	## if not cutdim: dimnamesv.append(objvars.dimensions[idd])
766	dimnamesv = [vdimyn, vdimxn]
767
768	varunits = []
769	varunits.append(objvars.getncattr('units'))
770	varunits.append(objvarc.getncattr('units'))
771
772	if not objsf.variables.has_key(vdimxn):
773	print errormsg
774	print ' ' + fname + ': shading file "' + ncfiles + \
775	'" does not have dimension variable "' + vdimxn + '" !!'
776	quit(-1)
777	if not objsf.variables.has_key(vdimyn):
778	print errormsg
779	print ' ' + fname + ': shading file "' + ncfiles + \
780	'" does not have dimension variable "' + vdimyn + '" !!'
781	quit(-1)
782
783	timename = timevals.split('\|')[0]
784	timeunit = timevals.split('\|')[1].replace('!',' ')
785	timekind = timevals.split('\|')[2]
786	timefmt = timevals.split('\|')[3]
787	timelabel = timevals.split('\|')[4].replace('!',' ')
788
789	if vdimxn == timename:
790	timevals = objsf.variables[vdimxn][:]
791	timedims = objsf.variables[vdimxn].dimensions
792	dimt = 'x'
793	ovalaxis = objsf.variables[vdimyn]
794	ovalu = ovalaxis.getncattr('units')
795	elif vdimyn == timename:
796	timevals = objsf.variables[vdimyn][:]
797	timedims = objsf.variables[vdimyn].dimensions
798	dimt = 'y'
799	ovalaxis = objsf.variables[vdimxn]
800	ovalu = ovalaxis.getncattr('units')
801	else:
802	print errormsg
803	print ' ' + fname + ": time variable '" + timename + "' not found!!"
804	quit(-1)
805
806	timepos, timelabels = drw.CFtimes_plot(timevals, timeunit, timekind, timefmt)
807
808	# Getting only that dimensions with coincident names
809	dimnvx = ovalaxis.dimensions
810
811	cutslice = []
812	for idimn in dimsshad:
813	found = False
814	for dimsn in dimnvx:
815	if idimn == dimsn:
816	cutslice.append(slice(0,len(objsf.dimensions[idimn])))
817	found = True
818	if not found: cutslice.append(0)
819
820	ovalaxisv = ovalaxis[tuple(cutslice)]
821
822	## if len(ovalaxis.shape) <= 2:
823	## ovalaxisv = ovalaxis[:]
824
825	## elif len(ovalaxis.shape) == 3:
826	## ovalaxisv = ovalaxis[0,:]
827	## else:
828	## print errormsg
829	## print ' ' + fname + ': shape of dimension variable:', ovalaxis.shape, \
830	## ' not ready!!'
831	## quit(-1)
832
833	if countlabelfmt == 'None':
834	countlfmt = None
835	else:
836	countlfmt = countlabelfmt
837
838	shading_nx = np.zeros((2), dtype=np.float)
839	shading_nx[0] = np.float(shadminmax.split(',')[0])
840	shading_nx[1] = np.float(shadminmax.split(',')[1])
841
842	clevmin = np.float(contlevels.split(',')[0])
843	clevmax = np.float(contlevels.split(',')[1])
844	Nclevels = int(contlevels.split(',')[2])
845
846	levels_cont = drw.pretty_int(clevmin, clevmax, Nclevels)
847
848	if len(levels_cont) <= 1:
849	print warnmsg
850	print ' ' + fname + ': wrong contour levels:', levels_cont,' !!'
851	del(levels_cont)
852	levels_cont = np.zeros((Nclevels), dtype=np.float)
853	levels_cont = np.arange(7)*(clevmax - clevmin)/(Nclevels-1)
854	print ' generating default ones: ',levels_cont
855
856	if mapvalue == 'None': mapvalue = None
857
858	drw.plot_2D_shadow_contour_time(valshad, valcont, vnamesfig, ovalaxisv, \
859	timevals, timepos, timelabels, ovalu, timelabel, dimt, dimnamesv, colbarn, \
860	countkind, countlfmt, shading_nx, levels_cont, varunits, figtitle, figkind, \
861	revals, mapvalue)
862
863	return
864
865	def draw_2D_shad_line(ncfile, values, varn):
866	""" plotting a fields with shading and another with line
867	draw_2D_shad_line(ncfile, values, varn)
868	ncfile= [ncfiles],[ncfilel] file to use for the shading and for the line
869	values=[vnamefs],[vnamefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar],[colline]:[sminv],[smaxv]:[figt]:
870	[kindfig]:[reverse]:[mapv]:[close]
871	[vnamefs]: Name in the figure of the variable to be shaded
872	[vnamefl]: Name in the figure of the variable to be lined
873	[dimvals]: ',' list of [dimname]\|[value] telling at which dimension of the
874	variable a given value is required (-1, all the length)
875	[dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
876	[colorbar]: name of the color bar
877	[colline]: name of the color for the line
878	[smin/axv]: minimum and maximum value for the shading
879	[figt]: title of the figure ('\|' for spaces)
880	[kindfig]: kind of figure
881	[reverse]: Transformation of the values
882	* 'transpose': reverse the axes (x-->y, y-->x)
883	* 'flip'@[x/y]: flip the axis x or y
884	[mapv]: map characteristics: [proj],[res]
885	see full documentation: http://matplotlib.org/basemap/
886	[proj]: projection
887	* 'cyl', cilindric
888	* 'lcc', lamvbert conformal
889	[res]: resolution:
890	* 'c', crude
891	* 'l', low
892	* 'i', intermediate
893	* 'h', high
894	* 'f', full
895	valules= 'rh:z\|-1,x\|-1:z\|-1,x\|-1:lat:pressure:BuPu:0.,100.:rh:pdf:flip@y:None'
896	varn= [varsn],[varnl] name of the variable to plot with shading and with line
897	"""
898
899	fname = 'draw_2D_shad_line'
900	if values == 'h':
901	print fname + '_____________________________________________________________'
902	print draw_2D_shad_line.__doc__
903	quit()
904
905	farguments = '[vnamefs],[vnamefl]:[dimvals]:[dimxvn]:[dimyvn]:' + \
906	'[colorbar],[colline]:[sminv],[smaxv]:[figt]:[kindfig]:[reverse]:' + \
907	'[mapv]:[close]'
908	drw.check_arguments(fname,values,farguments,':')
909
910	vnamesfig = values.split(':')[0].split(',')[0]
911	dimvals= values.split(':')[1].replace('\|',':')
912	vdimxn = values.split(':')[2]
913	vdimyn = values.split(':')[3]
914	colbarn = values.split(':')[4].split(',')[0]
915	shadminmax = values.split(':')[5]
916	figtitle = values.split(':')[6].replace('\|',' ')
917	figkind = values.split(':')[7]
918	revals = values.split(':')[8]
919	mapvalue = values.split(':')[9]
920	# varn = values.split(':')[10]
921
922	ncfiles = ncfile.split(',')[0]
923
924	if not os.path.isfile(ncfiles):
925	print errormsg
926	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
927	quit(-1)
928
929	objsf = NetCDFFile(ncfiles, 'r')
930
931	varns = varn.split(',')[0]
932
933	if not objsf.variables.has_key(varns):
934	print errormsg
935	print ' ' + fname + ': shading file "' + ncfiles + \
936	'" does not have variable "' + varns + '" !!'
937	quit(-1)
938
939	# Variables' values
940	objvars = objsf.variables[varns]
941
942	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
943
944	# Dimensions names
945	## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
946	## dimnamesv = []
947	## for idd in range(len(objvars.dimensions)):
948	## cutdim = False
949	## for idc in range(len(dimvals.split(','))):
950	## dimcutn = dimvals.split(',')[idc].split(':')[0]
951	## print objvars.dimensions[idd], dimcutn
952	## if objvars.dimensions[idd] == dimcutn:
953	## cutdim = True
954	## break
955	## if not cutdim: dimnamesv.append(objvars.dimensions[idd])
956	dimnamesv = [vdimyn, vdimxn]
957
958	varunits = objvars.getncattr('units')
959
960	if not objsf.variables.has_key(vdimxn):
961	print errormsg
962	print ' ' + fname + ': shading file "' + ncfiles + \
963	'" does not have dimension variable "' + vdimxn + '" !!'
964	quit(-1)
965	if not objsf.variables.has_key(vdimyn):
966	print errormsg
967	print ' ' + fname + ': shading file "' + ncfiles + \
968	'" does not have dimensino variable "' + vdimyn + '" !!'
969	quit(-1)
970
971	objdimx = objsf.variables[vdimxn]
972	objdimy = objsf.variables[vdimyn]
973	odimxu = objdimx.getncattr('units')
974	odimyu = objdimy.getncattr('units')
975
976	if len(objdimx.shape) <= 2:
977	# odimxv = objdimx[valshad.shape]
978	# odimyv = objdimy[valshad.shape]
979	odimxv = objdimx[:]
980	odimyv = objdimy[:]
981
982	elif len(objdimx.shape) == 3:
983	# dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])]
984	# odimxv = objdimx[tuple(dimcut)]
985	# odimyv = objdimy[tuple(dimcut)]
986	odimxv = objdimx[0,:]
987	odimyv = objdimy[0,:]
988	else:
989	print errormsg
990	print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
991	' not ready!!'
992	quit(-1)
993
994	shading_nx = np.zeros((2), dtype=np.float)
995	shading_nx[0] = np.float(shadminmax.split(',')[0])
996	shading_nx[1] = np.float(shadminmax.split(',')[1])
997
998	if mapvalue == 'None': mapvalue = None
999
1000	# line plot
1001	##
1002	ncfilel = ncfile.split(',')[1]
1003	vnamelfig = values.split(':')[0].split(',')[1]
1004	varnl = varn.split(',')[1]
1005	colline = values.split(':')[4].split(',')[1]
1006
1007	objlf = NetCDFFile(ncfilel,'r')
1008	objlvar = objlf.variables[varnl]
1009
1010	linevals = objlvar[:]
1011	varlunits = objlvar.units
1012
1013	drw.plot_2D_shadow_line(valshad, linevals, vnamesfig, vnamelfig, odimxv, odimyv, \
1014	odimxu, odimyu, dimnamesv, colbarn, colline, shading_nx, varunits, varlunits, \
1015	figtitle, figkind, revals, mapvalue, True)
1016
1017	objsf.close()
1018	objlf.close()
1019
1020	return
1021
1022	def draw_2D_shad_line_time(ncfile, values, varn):
1023	""" plotting a fields with shading and a line with time values
1024	draw_2D_shad_line(ncfile, values, varn)
1025	ncfile= [ncfiles],[ncfilel] files to use to draw with shading and the line
1026	values= [vnamefs],[vanemefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
1027	[kindfig]:[reverse]:[timevals]:[close]
1028	[vnamefs]: Name in the figure of the variable to be shaded
1029	[vnamefl]: Name in the figure of the variable to be lined
1030	[dimvals]: ',' list of [dimname]\|[value] telling at which dimension of the
1031	variable a given value is required (-1, all the length)
1032	[dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
1033	[colorbar]: name of the color bar
1034	[smin/axv]: minimum and maximum value for the shading
1035	[figt]: title of the figure ('\|' for spaces)
1036	[kindfig]: kind of figure
1037	[reverse]: Transformation of the values
1038	* 'transpose': reverse the axes (x-->y, y-->x)
1039	* 'flip'@[x/y]: flip the axis x or y
1040	[timevals]: [timen]\|[units]\|[kind]\|[tfmt]\|[label]\|[timeaxis] time labels characteristics
1041	[timen]; name of the time variable
1042	[units]; units string according to CF conventions ([tunits] since
1043	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
1044	[kind]; kind of output
1045	'Nval': according to a given number of values as 'Nval',[Nval]
1046	'exct': according to an exact time unit as 'exct',[tunit];
1047	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1048	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1049	'l': milisecond
1050	[tfmt]; desired format
1051	[label]; label at the graph ('!' for spaces)
1052	[close]: should figure be closed (finished)
1053	values='dtcon,prc:Time\|-1,bottom_top\|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly\|'
1054	'dtcon:pdf:transpose:time\|hours!since!1949-12-01\|exct,2,d\|%d\|date!([DD])\|x:True
1055	varn= [varsn].[varln] name of the variable to plot with shading and to plot with line
1056	"""
1057	fname = 'draw_2D_shad_line_time'
1058	if values == 'h':
1059	print fname + '_____________________________________________________________'
1060	print draw_2D_shad__line_time.__doc__
1061	quit()
1062
1063	farguments = '[vnamefs],[vanemefl]:[dimvals]:[dimxvn]:[dimyvn]:' + \
1064	'[colorbar]:[sminv],[smaxv]:[figt]:[kindfig]:[reverse]:[timevals]:[close]'
1065	drw.check_arguments(fname,values,farguments,':')
1066
1067	vnamesfig = values.split(':')[0].split(',')[0]
1068	dimvals= values.split(':')[1].replace('\|',':')
1069	vdimxn = values.split(':')[2]
1070	vdimyn = values.split(':')[3]
1071	colbarn = values.split(':')[4]
1072	shadminmax = values.split(':')[5]
1073	figtitle = values.split(':')[6].replace('\|',' ')
1074	figkind = values.split(':')[7]
1075	revals = values.split(':')[8]
1076	timevals = values.split(':')[9]
1077	close = values.split(':')[10]
1078
1079	ncfiles = ncfile.split(',')[0]
1080
1081	if not os.path.isfile(ncfiles):
1082	print errormsg
1083	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
1084	quit(-1)
1085
1086	objsf = NetCDFFile(ncfiles, 'r')
1087
1088	varns = varn.split(',')[0]
1089
1090	if not objsf.variables.has_key(varns):
1091	print errormsg
1092	print ' ' + fname + ': shading file "' + ncfiles + \
1093	'" does not have variable "' + varns + '" !!'
1094	quit(-1)
1095
1096	# Variables' values
1097	objvars = objsf.variables[varns]
1098
1099	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
1100
1101	dimnamesv = [vdimyn, vdimxn]
1102
1103	varunits = objvars.getncattr('units')
1104
1105	if not objsf.variables.has_key(vdimxn):
1106	print errormsg
1107	print ' ' + fname + ': shading file "' + ncfiles + \
1108	'" does not have dimension variable "' + vdimxn + '" !!'
1109	quit(-1)
1110	if not objsf.variables.has_key(vdimyn):
1111	print errormsg
1112	print ' ' + fname + ': shading file "' + ncfiles + \
1113	'" does not have dimensino variable "' + vdimyn + '" !!'
1114	quit(-1)
1115
1116	objdimx = objsf.variables[vdimxn]
1117	objdimy = objsf.variables[vdimyn]
1118	odimxu = objdimx.getncattr('units')
1119	odimyu = objdimy.getncattr('units')
1120
1121	if len(objdimx.shape) <= 2:
1122	odimxv = objdimx[:]
1123	odimyv = objdimy[:]
1124
1125	elif len(objdimx.shape) == 3:
1126	odimxv = objdimx[0,:]
1127	odimyv = objdimy[0,:]
1128	else:
1129	print errormsg
1130	print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
1131	' not ready!!'
1132	quit(-1)
1133
1134	timename = timevals.split('\|')[0]
1135	timeunit = timevals.split('\|')[1].replace('!',' ')
1136	timekind = timevals.split('\|')[2]
1137	timefmt = timevals.split('\|')[3]
1138	timelabel = timevals.split('\|')[4].replace('!',' ')
1139
1140	if vdimxn == timename:
1141	odimxv = objsf.variables[vdimxn][:]
1142	odimxu = timelabel
1143	timeaxis = 'x'
1144	odimyv = objsf.variables[vdimyn]
1145	odimyu = odimyv.getncattr('units')
1146	timepos, timelabels = drw.CFtimes_plot(odimxv, timeunit, timekind, timefmt)
1147	elif vdimyn == timename:
1148	odimyv = objsf.variables[vdimyn][:]
1149	odimyu = timelabel
1150	timeaxis = 'y'
1151	odimxv = objsf.variables[vdimxn]
1152	odimxu = odimxv.getncattr('units')
1153	timepos, timelabels = drw.CFtimes_plot(odimyv, timeunit, timekind, timefmt)
1154	else:
1155	print errormsg
1156	print ' ' + fname + ": time variable '" + timename + "' not found!!"
1157	quit(-1)
1158
1159	shading_nx = np.zeros((2), dtype=np.float)
1160	shading_nx[0] = np.float(shadminmax.split(',')[0])
1161	shading_nx[1] = np.float(shadminmax.split(',')[1])
1162
1163	closeval = drw.Str_Bool(close)
1164
1165	drw.plot_2D_shadow_time(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu, \
1166	dimnamesv, colbarn, shading_nx, varunits, figtitle, figkind, revals, timeaxis, \
1167	timepos, timelabels, False)
1168
1169	# Line values
1170	##
1171	ncfilel = ncfile.split(',')[1]
1172
1173	vnamelfig = values.split(':')[0].split(',')[1]
1174	varnl = varn.split(',')[1]
1175
1176	objlf = NetCDFFile(ncfilel,'r')
1177	objlvar = objlf.variables[varnl]
1178
1179	linevals = objlvar[:]
1180	if reva0 == 'tranpose':
1181	plt.plot (linevals, odimxv, '-', color='k')
1182	else:
1183	plt.plot (odimxv, linevals, '-', color='k')
1184
1185	objsf.close()
1186	objsl.close()
1187
1188	return
1189
1190	def draw_barbs(ncfile, values, varns):
1191	""" Function to plot wind barbs
1192	values= [dimname]\|[vardimname]\|[value]:[vecvals]:[windlabs]:[mapvalues]:
1193	[gtit]:[kindfig]:[figuren]
1194	'X/Y/Z/T'\|[dimname]\|[vardimname]\|[value]: ',', list for each basic dimension '\|' separated of:
1195	[dimname]: name of the dimension in the file
1196	[vardimname]: name of the variable with the values for the dimension in the file
1197	[value]: which value of the given dimension (-1, all; [ibeg]@[iend], i-range beginning, end)
1198	No value takes all the range of the dimension
1199	[vecvals]= [frequency],[color],[length]
1200	[frequency]: [xfreq]@[yfreq] frequency of values allong each axis ('None', all grid points;
1201	'auto', computed automatically to have 20 vectors along each axis)
1202	[color]: color of the vectors ('auto', for 'red')
1203	[length]: length of the wind barbs ('auto', for 9)
1204	[windlabs]= [windname],[windunits]
1205	[windname]: name of the wind variable in the graph
1206	[windunits]: units of the wind variable in the graph ('None', for the value in the file)
1207	[mapvalues]= map characteristics: [proj],[res]
1208	see full documentation: http://matplotlib.org/basemap/
1209	[proj]: projection
1210	* 'cyl', cilindric
1211	* 'lcc', lambert conformal
1212	[res]: resolution:
1213	* 'c', crude
1214	* 'l', low
1215	* 'i', intermediate
1216	* 'h', high
1217	* 'f', full
1218	gtit= title of the graph ('\|', for spaces)
1219	kindfig= kind of figure
1220	figuren= name of the figure
1221	ncfile= file to use
1222	varns= [uwind],[ywind] ',' list of the name of the variables with the u-wind,y-wind component
1223	"""
1224	fname = 'draw_barbs'
1225
1226	if values == 'h':
1227	print fname + '_____________________________________________________________'
1228	print draw_barbs.__doc__
1229	quit()
1230
1231	expectargs = '[X/Y/Z/T]\|[dimname]\|[vardimname]\|[value]:[vecvals]:[windlabs]:' + \
1232	'[mapvalues]:[gtit]:[kindfig]:[figuren]'
1233
1234	drw.check_arguments(fname,values,expectargs,':')
1235
1236	dimvals = values.split(':')[0]
1237	vecvals = values.split(':')[1]
1238	windlabels = values.split(':')[2]
1239	mapvalues = values.split(':')[3]
1240	gtit = values.split(':')[4]
1241	kindfig = values.split(':')[5]
1242	figuren = values.split(':')[6]
1243
1244	of = NetCDFFile(ncfile,'r')
1245
1246	dims = {}
1247	for dimv in dimvals.split(','):
1248	dns = dimv.split('\|')
1249	dims[dns[0]] = [dns[1], dns[2], dns[3]]
1250
1251	varNs = []
1252	for dn in dims.keys():
1253	if dn == 'X':
1254	varNs.append(dims[dn][1])
1255	dimx = len(of.dimensions[dims[dn][0]])
1256	elif dn == 'Y':
1257	varNs.append(dims[dn][1])
1258	dimy = len(of.dimensions[dims[dn][0]])
1259
1260	ivar = 0
1261	for wvar in varns.split(','):
1262	if not drw.searchInlist(of.variables.keys(), wvar):
1263	print errormsg
1264	print ' ' + fname + ": file does not have variable '" + wvar + "' !!"
1265	quit(-1)
1266	if ivar == 0:
1267	varNs.append(wvar)
1268	else:
1269	varNs.append(wvar)
1270
1271	ivar = 0
1272	for varN in varNs:
1273	varslice = []
1274
1275	ovarN = of.variables[varN]
1276	vard = ovarN.dimensions
1277	for vdn in vard:
1278	found = False
1279	for dd in dims.keys():
1280	if dims[dd][0] == vdn:
1281	if dims[dd][2].find('@') != -1:
1282	rvals = dims[dd][2].split('@')
1283	varslice.append(slice(int(rvals[0]), int(rvals[1])))
1284	elif dims[dd][2] == '-1':
1285	varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
1286	else:
1287	varslice.append(int(dims[dd][2]))
1288
1289	found = True
1290	break
1291	if not found:
1292	varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
1293
1294	if varN == dims['X'][1]:
1295	lonvals0 = np.squeeze(ovarN[tuple(varslice)])
1296	elif varN == dims['Y'][1]:
1297	latvals0 = np.squeeze(ovarN[tuple(varslice)])
1298	elif ivar == 2:
1299	uwvals = np.squeeze(np.array(ovarN[tuple(varslice)]))
1300	elif ivar == 3:
1301	vwvals = np.squeeze(ovarN[tuple(varslice)])
1302
1303	ivar = ivar + 1
1304
1305	# print 'Final shapes:',lonvals0.shape,':',latvals0.shape,':',uwvals.shape,':',
1306	# vwvals.shape
1307
1308	if len(uwvals.shape) != 2 or len(vwvals.shape) != 2:
1309	print errormsg
1310	print ' ' + fname + ': wrong size of the wind fields! they must be ' + \
1311	'2-dimensional!'
1312	print ' u-winds shape:',uwvals.shape,'dims:',of.variables[varNs[2]]
1313	print ' v-winds shape:',vwvals.shape,'dims:',of.variables[varNs[3]]
1314	print ' provide more values for their dimensions!!'
1315	quit(-1)
1316
1317	if len(lonvals0.shape) == 1:
1318	lonvals, latvals = np.meshgrid(lonvals0, latvals0)
1319	else:
1320	lonvals = lonvals0
1321	latvals = latvals0
1322
1323	# Vecor values
1324	if vecvals.split(',')[0] == 'None':
1325	freqv = None
1326	else:
1327	freqv = vecvals.split(',')[0]
1328	colorv = vecvals.split(',')[1]
1329	lengthv = vecvals.split(',')[2]
1330
1331	# Vector labels
1332	windname = windlabels.split(',')[0]
1333	windunits = windlabels.split(',')[1]
1334
1335	drw.plot_barbs(lonvals, latvals, uwvals, vwvals, freqv, colorv, lengthv,
1336	windname, windunits, mapvalues, gtit, kindfig, figuren)
1337
1338	return
1339
1340	def draw_topo_geogrid(ncfile, values):
1341	""" plotting geo_em.d[nn].nc topography from WPS files
1342	draw_topo_geogrid(ncfile, values)
1343	ncfile= geo_em.d[nn].nc file to use
1344	values= [minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]
1345	[min/max]Topo: minimum and maximum values of topography to draw
1346	lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] or None
1347	title: title of the graph
1348	graphic_kind: kind of figure (jpg, pdf, png)
1349	mapvalues: map characteristics [proj],[res]
1350	see full documentation: http://matplotlib.org/basemap/
1351	[proj]: projection
1352	* 'cyl', cilindric
1353	* 'lcc', lambert conformal
1354	[res]: resolution:
1355	* 'c', crude
1356	* 'l', low
1357	* 'i', intermediate
1358	* 'h', high
1359	* 'f', full
1360	"""
1361	fname = 'draw_topo_geogrid'
1362
1363	if values == 'h':
1364	print fname + '_____________________________________________________________'
1365	print draw_topo_geogrid.__doc__
1366	quit()
1367
1368	expectargs = '[minTopo]:[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]'
1369
1370	drw.check_arguments(fname,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,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];[linewidths];[pointsizes];[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	[lines]: ',' list of type of lines, None for automatic, single value all the same
1839	[collines]: ',' list of colors for the lines, None for automatic, single
1840	value all the same
1841	[points]: ',' list of type of points for the lines, None for automatic, single
1842	value all the same
1843	[linewidths]: ',' list of widths for the lines, None for automatic, single
1844	value all the same
1845	[pointsizes]: ',' list of widths for the lines, None for automatic, single
1846	value all the same
1847	[pointfreq]: frequency of point plotting, 'all' for all time steps
1848	[period]: which period to plot
1849	'-1': all period
1850	[beg],[end]: beginning and end of the period in reference time-units of first file
1851	varname0= ',' list of variable names to plot (assuming only 1 variable per file)
1852	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'
1853	"""
1854
1855	fname = 'draw_lines_time'
1856
1857	if values == 'h':
1858	print fname + '_____________________________________________________________'
1859	print draw_lines_time.__doc__
1860	quit()
1861
1862	expectargs = '[dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];'
1863	expectargs = expectargs + '[timevals];[locleg];[graphk];[lines];[collines];[points];'
1864	expectargs = expectargs + '[linewidths];[pointsizes];[pointfreq];[period]'
1865	drw.check_arguments(fname,values,expectargs,';')
1866
1867	ncfiles = ncfilens.split(',')
1868	dimvname0 = values.split(';')[0]
1869	valuesaxis = values.split(';')[1]
1870	dimtit = values.split(';')[2].replace('\|',' ')
1871	leglabels = values.split(';')[3].replace('_','\_')
1872	vartit = values.split(';')[4]
1873	title = values.split(';')[5].replace('\|',' ')
1874	timevals = values.split(';')[6]
1875	locleg = int(values.split(';')[7])
1876	graphk = values.split(';')[8]
1877	lines0 = values.split(';')[9]
1878	collines0 = values.split(';')[10]
1879	points0 = values.split(';')[11]
1880	linewidths0 = values.split(';')[12]
1881	pointsizes0 = values.split(';')[13]
1882	pointfreq0 = values.split(';')[14]
1883	period = values.split(';')[15]
1884
1885	Nfiles = len(ncfiles)
1886
1887	# Multiple variable-dimension names?
1888	if dimvname0.find(',') != -1:
1889	dimvname = dimvname0.split(',')
1890	else:
1891	dimvname = [dimvname0]
1892
1893	# Multiple variables?
1894	if varname0.find(',') != -1:
1895	varname = varname0.split(',')
1896	else:
1897	varname = [varname0]
1898
1899	# Multiple lines types?
1900	if lines0.find(',') != -1:
1901	lines = lines0.split(',')
1902	elif lines0 == 'None':
1903	lines = None
1904	else:
1905	lines = []
1906	for il in range(Nfiles):
1907	lines.append(lines0)
1908
1909	# Multiple color names?
1910	if collines0.find(',') != -1:
1911	collines = collines0.split(',')
1912	elif collines0 == 'None':
1913	collines = None
1914	else:
1915	collines = []
1916	for ip in range(Nfiles):
1917	collines.append(collines0)
1918
1919	# Multiple point types?
1920	if points0.find(',') != -1:
1921	points = points0.split(',')
1922	elif points0 == 'None':
1923	points = None
1924	else:
1925	points = []
1926	for ip in range(Nfiles):
1927	points.append(points0)
1928
1929	# Multiple line sizes?
1930	if linewidths0.find(',') != -1:
1931	linewidths = []
1932	Nlines = len(linewidths0.split(','))
1933	for il in Nlines:
1934	linewidths.append(np.float(linewidths0.split(',')[il]))
1935	elif linewidths0 == 'None':
1936	linewidths = None
1937	else:
1938	linewidths = [np.float(linewidths0)]
1939
1940	# Multiple point sizes?
1941	if pointsizes0.find(',') != -1:
1942	pointsizes = []
1943	Npts = len(pointsizes0.split(','))
1944	for ip in Npts:
1945	pointsizes.append(np.float(pointsizes0.split(',')[ip]))
1946	elif pointsizes0 == 'None':
1947	pointsizes = None
1948	else:
1949	pointsizes = [np.float(pointsizes0)]
1950
1951	# Getting values
1952	##
1953	varvalues = []
1954	dimvalues = []
1955	timvalues = []
1956	timvals0 = timvalues
1957
1958	print ' ' + fname
1959	ifn = 0
1960	mintval = 1.e20
1961	maxtval = -1.e20
1962
1963	for ifile in ncfiles:
1964	filen = ifile.split('@')[0]
1965
1966	print ' filen:',filen
1967
1968	if not os.path.isfile(filen):
1969	print errormsg
1970	print ' ' + fname + ": netCDF file '" + filen + "' does not exist !!"
1971	quit(-1)
1972
1973	objfile = NetCDFFile(filen, 'r')
1974
1975	founddvar = False
1976	for dvar in dimvname:
1977	if objfile.variables.has_key(dvar):
1978	founddvar = True
1979	vdobj = objfile.variables[dvar]
1980	if len(vdobj.shape) != 1:
1981	print errormsg
1982	print ' ' + fname + ': wrong shape:',vdobj.shape," of " + \
1983	"variable '" + dvar + "' !!"
1984	quit(-1)
1985	break
1986	if not founddvar:
1987	print errormsg
1988	print ' ' + fname + ": netCDF file '" + filen + \
1989	"' has any variable '", dimvname, "' !!"
1990	quit(-1)
1991
1992	foundvar = False
1993	for var in varname:
1994	if objfile.variables.has_key(var):
1995	foundvar = True
1996	vvobj = objfile.variables[var]
1997	if len(vvobj.shape) != 1:
1998	print errormsg
1999	print ' ' + fname + ': wrong shape:',vvobj.shape," of " + \
2000	"variable '" + var + "' !!"
2001	quit(-1)
2002
2003	break
2004	if not foundvar:
2005	print errormsg
2006	print ' ' + fname + ": netCDF file '" + filen + \
2007	"' has any variable '", varname, "' !!"
2008	quit(-1)
2009
2010
2011	# Getting period
2012	dimt = len(vdobj[:])
2013
2014	if period == '-1':
2015	varvalues.append(vvobj[:])
2016	dimvalues.append(vdobj[:])
2017	mindvals = np.min(vdobj[:])
2018	maxdvals = np.max(vdobj[:])
2019	else:
2020	ibeg=-1
2021	iend=-1
2022	tbeg = np.float(period.split(',')[0])
2023	tend = np.float(period.split(',')[1])
2024
2025	for it in range(dimt-1):
2026	if vdobj[it] <= tbeg and vdobj[it+1] > tbeg : ibeg = it
2027	if vdobj[it] <= tend and vdobj[it+1] > tend: iend = it + 1
2028	if ibeg != -1 and iend != -1: break
2029
2030	if ibeg == -1 and iend == -1:
2031	print errormsg
2032	print ' ' + fname + ': Period:',tbeg,',',tend,'not found!!'
2033	print ' ibeg:',ibeg,'iend:',iend
2034	print ' period in file:',np.min(vdobj[:]), np.max(vdobj[:])
2035	quit(-1)
2036	elif iend == -1:
2037	iend = dimt
2038	print warnmsg
2039	print ' ' + fname + ': end of Period:',tbeg,',',tend,'not found!!'
2040	print ' getting last available time instead'
2041	print ' ibeg:',ibeg,'iend:',iend
2042	print ' period in file:',np.min(vdobj[:]), np.max(vdobj[:])
2043	elif ibeg == -1:
2044	ibeg = 0
2045	print warnmsg
2046	print ' ' + fname + ': beginning of Period:',tbeg,',',tend, \
2047	'not found!!'
2048	print ' getting first available time instead'
2049	print ' ibeg:',ibeg,'iend:',iend
2050	print ' period in file:',np.min(vdobj[:]), np.max(vdobj[:])
2051
2052	varvalues.append(vvobj[ibeg:iend])
2053	dimvalues.append(vdobj[ibeg:iend])
2054	mindvals = np.min(vdobj[ibeg:iend])
2055	maxdvals = np.max(vdobj[ibeg:iend])
2056
2057	dimt = iend - ibeg
2058
2059	if mindvals < mintval: mintval = mindvals
2060	if maxdvals > maxtval: maxtval = maxdvals
2061
2062	if ifn == 0:
2063	varunits = drw.units_lunits(vvobj.units)
2064
2065	objfile.close()
2066
2067	ifn = ifn + 1
2068
2069	# Times
2070	timename = timevals.split('\|')[0]
2071	timeunit = timevals.split('\|')[1].replace('!',' ')
2072	timekind = timevals.split('\|')[2]
2073	timefmt = timevals.split('\|')[3]
2074
2075	dtvals = (maxtval - mintval)/dimt
2076	tvals = np.arange(mintval, maxtval, dtvals/2.)
2077
2078	timepos, timelabels = drw.CFtimes_plot(tvals, timeunit, timekind, timefmt)
2079
2080	if leglabels != 'None':
2081	legvals = leglabels.split(',')
2082	else:
2083	legvals = None
2084
2085	if pointfreq0 == 'all':
2086	pointfreq = None
2087	else:
2088	pointfreq = int(pointfreq0)
2089
2090	drw.plot_lines_time(dimvalues, varvalues, valuesaxis, dimtit, legvals, vartit, \
2091	varunits, timepos, timelabels, title, locleg, graphk, lines, collines, points, \
2092	linewidths, pointsizes, pointfreq)
2093
2094	return
2095
2096	def draw_Neighbourghood_evol(filen, values, variable):
2097	""" Function to draw the temporal evolution of a neighbourghood around a point
2098	draw_Neighbourghood_evol(filen, values, variable)
2099	filen= netCDF file name
2100	values= [gvarname]:[dimsval]:[neigdims]:[Nneig]:[Ncol]:[timetits]:[tkinds]:
2101	[timefmts]:[gtitle]:[shadxtrms]:[cbar]:[gkind]:[ofile]
2102	[dimsval]: [dimn1]\|[val1]\|[dimv1],...,[dimnN]\|[valN]\|[dimvN] dimension names, values to get
2103	(-1, for all; no name/value pair given full length) and variable with values of the dimension
2104	NOTE: when dimsval[X,Y] == neigdims[X,Y], valX,valY --> valX,valY-Nneig/2, valX,valY+Nneig/2
2105	[neigdims]: [dimnX],[dimnY] dimensions mnames along which the neigbourghood should be defined
2106	[Nneig]: Number of grid points of the full side of the box (odd value)
2107	[Ncol]: Number of columns ('auto': square final plot)
2108	[gvarname]: name of the variable to appear in the graph
2109	[timetits]: [titX],[titY] titles of the axes ('\|' for spaces)
2110	[tkinds]: [tkindX]\|[tkindY] kinds of time to appear in the graph
2111	'Nval': according to a given number of values as 'Nval',[Nval]
2112	'exct': according to an exact time unit as 'exct',[tunit];
2113	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2114	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2115	'l': milisecond
2116	[timefmts]: [tfmtX],[tfmtY] format of the time labels
2117	[gtitle]: title of the graphic ('\|' for spaces)
2118	[shadxtrms]: Extremes for the shading
2119	[cbar]: colorbar to use
2120	[gkind]: kind of graphical output
2121	[ofile]: True/False whether the netcdf with data should be created or not
2122	variable= name of the variable
2123	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'
2124	"""
2125
2126	fname = 'draw_Neighbourghood_evol'
2127
2128	if values == 'h':
2129	print fname + '_____________________________________________________________'
2130	print draw_Neighbourghood_evol.__doc__
2131	quit()
2132
2133	expectargs = '[gvarname]:[dimsval]:[neigdims]:[Nneig]:[Ncol]:' + \
2134	'[timetits]:[tkinds]:[timefmts]:[gtitle]:[shadxtrms]:[cbar]:[gkind]:[ofile]'
2135
2136	drw.check_arguments(fname,values,expectargs,':')
2137
2138	gvarname = values.split(':')[0]
2139	dimsval = values.split(':')[1].split(',')
2140	neigdims = values.split(':')[2].split(',')
2141	Nneig = int(values.split(':')[3])
2142	Ncol0 = values.split(':')[4]
2143	timetits = values.split(':')[5].split(',')
2144	timekinds = values.split(':')[6].split('\|')
2145	timefmts = values.split(':')[7].split(',')
2146	gtitle = values.split(':')[8].replace('\|',' ')
2147	shadxtrms = values.split(':')[9].split(',')
2148	cbar = values.split(':')[10]
2149	gkind = values.split(':')[11]
2150	ofile = values.split(':')[12]
2151
2152	if Ncol0 != 'auto':
2153	Ncol = int(Ncol0)
2154	else:
2155	Ncol = Ncol0
2156
2157	timetits[0] = timetits[0].replace('\|',' ')
2158	timetits[1] = timetits[1].replace('\|',' ')
2159
2160	if np.mod(Nneig,2) == 0:
2161	print errormsg
2162	print ' ' + fname + ": an odd value for 'Nneig':",Nneig,'is required !!!'
2163	quit(-1)
2164
2165	Nneig2 = int(Nneig/2)
2166
2167	# Values to slice the variable
2168	dimvslice = {}
2169	dimvvalues = {}
2170	for dimvs in dimsval:
2171	dimn = dimvs.split('\|')[0]
2172	dimv = int(dimvs.split('\|')[1])
2173	dimnv = dimvs.split('\|')[2]
2174
2175	dimvvalues[dimn] = dimnv
2176	dimvslice[dimn] = dimv
2177
2178	ncobj = NetCDFFile(filen, 'r')
2179
2180	varobj = ncobj.variables[variable]
2181
2182	slicevar = []
2183	newdimn = []
2184	newdimsvar = {}
2185
2186	for dimn in varobj.dimensions:
2187	if not drw.searchInlist(dimvslice.keys(), dimn):
2188	dimsize = len(ncobj.dimensions[dimn])
2189	slicevar.append(slice(0, dimsize+1))
2190	newdimn.append(dimn)
2191	newdimsvar[dimn] = dimseize
2192
2193	for dimslicen in dimvslice.keys():
2194	if dimn == dimslicen:
2195	if dimvslice[dimn] != -1:
2196	if drw.searchInlist(neigdims, dimn):
2197	slicevar.append(slice(dimvslice[dimn]-Nneig2, \
2198	dimvslice[dimn]+Nneig2+1))
2199	newdimn.append(dimn)
2200	newdimsvar[dimn] = Nneig
2201	break
2202	else:
2203	slicevar.append(slice(dimvslice[dimn], dimvslice[dimn]+1))
2204	break
2205	else:
2206	dimsize = len(ncobj.dimensions[dimn])
2207	slicevar.append(slice(0, dimsize+1))
2208	newdimn.append(dimn)
2209	newdimsvar[dimn] = dimsize
2210	break
2211
2212	varv = varobj[tuple(slicevar)]
2213
2214	if len(newdimn) != 3:
2215	print errormsg
2216	print ' ' + fname + ': sliced variable with shape=', varv.shape, \
2217	' must have three dimensions',len(varv.shape),'given !!'
2218	quit(-1)
2219
2220	newdims = []
2221	for nwdims in newdimn:
2222	newdims.append(newdimsvar[nwdims])
2223
2224	# The dimension which is not in the neighbourhood dimensions must be time!
2225	for dim1 in newdimn:
2226	if not drw.searchInlist(neigdims, dim1):
2227	dimt = newdimsvar[dim1]
2228	dimtime = dim1
2229
2230	if Ncol == 'auto':
2231	dimtsqx = int(np.sqrt(dimt)) + 1
2232	dimtsqy = int(np.sqrt(dimt)) + 1
2233	else:
2234	dimtsqx = int(Ncol)
2235	dimtsqy = dimt/dimtsqx + 1
2236
2237	neighbourghood = np.ones((dimtsqyNneig,dimtsqxNneig), dtype=np.float)*fillValue
2238
2239	for it in range(dimt):
2240	ity = int(it/dimtsqx)
2241	itx = it-ity*dimtsqx
2242
2243	itty = (dimtsqy - ity - 1)*Nneig + Nneig2
2244	ittx = itx*Nneig + Nneig2
2245
2246	neighbourghood[itty-Nneig2:itty+Nneig2+1,ittx-Nneig2:ittx+Nneig2+1]= \
2247	varv[it,::-1,:]
2248
2249	variablevals = drw.variables_values(variable)
2250	if drw.searchInlist(varobj.ncattrs(), 'units'):
2251	vunits = varobj.units
2252	else:
2253	vunits = variablevals[5]
2254
2255	# Time values at the X/Y axes
2256	if ncobj.variables[dimvvalues[dimtime]].dtype == '\|S1':
2257	print ' ' + fname + ': WRF time variable!'
2258	refdate = '19491201000000'
2259	tunitsval = 'hours'
2260	dimtvalues = np.zeros((dimt), dtype=np.float)
2261	tvals = ncobj.variables[dimvvalues[dimtime]]
2262	yrref=refdate[0:4]
2263	monref=refdate[4:6]
2264	dayref=refdate[6:8]
2265	horref=refdate[8:10]
2266	minref=refdate[10:12]
2267	secref=refdate[12:14]
2268
2269	refdateS = yrref + '/' + monref + '/' + dayref + '_' + horref + ':' + \
2270	minref + ':' + secref
2271	tunits = tunitsval + ' since ' + refdateS
2272	for it in range(dimt):
2273	wrfdates = drw.datetimeStr_conversion(tvals[it,:],'WRFdatetime', 'matYmdHMS')
2274	dimtvalues[it] = drw.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
2275	else:
2276	dimtvalues = ncobj.variables[dimvvalues[dimtime]][:]
2277	tunits = ncobj.variables[newdimsvar[dimtime]].units
2278
2279	dimxv = dimtvalues[0:dimtsqx]
2280	dimyv = dimtvalues[0:dimt:dimtsqx]
2281
2282	dimn = ['time','time']
2283
2284	if ofile == 'True':
2285	ofilen = 'Neighbourghood_evol.nc'
2286	newnc = NetCDFFile(ofilen, 'w')
2287	# Dimensions
2288	newdim = newnc.createDimension('time',None)
2289	newdim = newnc.createDimension('y',dimtsqy*Nneig)
2290	newdim = newnc.createDimension('x',dimtsqx*Nneig)
2291	# Dimension values
2292	newvar = newnc.createVariable('time','f8',('time'))
2293	newvar[:] = np.arange(dimt)
2294	newattr = drw.basicvardef(newvar, 'time','time',tunits)
2295	# Neighbourhghood variable
2296	newvar = newnc.createVariable(variable + 'neigevol', 'f4', ('y','x'), \
2297	fill_value=fillValue)
2298	newvar[:] = neighbourghood
2299
2300	newnc.sync()
2301	newnc.close()
2302	print fname + ": Successfull generation of file '" + ofilen + "' !!"
2303
2304	# Time ticks
2305	timeposX, timelabelsX = drw.CFtimes_plot(dimxv, tunits, timekinds[0], timefmts[0])
2306	timeposY, timelabelsY = drw.CFtimes_plot(dimyv, tunits, timekinds[1], timefmts[1])
2307
2308	timepos = [timeposX[0:len(timeposX)], timeposY[len(timeposY):0:-1]]
2309	timelabels = [timelabelsX[0:len(timeposX)], timelabelsY[0:len(timeposY)]]
2310
2311	for i in range(2):
2312	if shadxtrms[i][0:1] != 'S':
2313	shadxtrms[i] = np.float(shadxtrms[i])
2314
2315	drw.plot_Neighbourghood_evol(neighbourghood, dimxv, dimyv, gvarname, timetits, \
2316	timepos, timelabels, cbar, Nneig, shadxtrms, vunits, gtitle, gkind, True)
2317
2318	def draw_points(filen, values):
2319	""" Function to plot a series of points
2320	[values]= [ptasciifile]:[gtit]:[mapvalues]:[kindfigure]:[pointcolor]:[pointlabels]:
2321	[locleg]:[figureko]:[figuren]
2322	[ptasciifile]:[file],[comchar],[collon],[collat],[lab]
2323	[file]: column ASCII file with the location of the points
2324	[comchar]: '\|' list of characters for commentaries
2325	[collon]: number of column with the longitude of the points
2326	[collat]: number of column with the latitude of the points
2327	[collab]: number of column with the labels of the points ('None', and will get
2328	the values from the [pointlabels] variable
2329	[gtit]: title of the figure ('\|' for spaces)
2330	[mapvalues]: drawing coastaline ([proj],[res]) or None
2331	[proj]: projection
2332	* 'cyl', cilindric
2333	* 'lcc', lambert conformal
2334	[res]: resolution:
2335	* 'c', crude
2336	* 'l', low
2337	* 'i', intermediate
2338	* 'h', high
2339	* 'f', full
2340	[kindfigure]: kind of figure
2341	'legend': only points in the map with the legend with the names
2342	'labelled',[txtsize],[txtcol]: points with the names and size, color of text
2343	[pointcolor]: color for the points ('auto' for "red")
2344	[pointlabels]: ',' of labels [only used if [collab]='None'] ('None' for no labels)
2345	[locleg]: location of the legend (0, autmoatic)
2346	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2347	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2348	9: 'upper center', 10: 'center'
2349	[figureko]: kind of the output file (pdf, png, ...)
2350	[figuren]: name of the figure
2351	[filen]= [ncfile],[lonvarn],[latvarn][,[varn],[dimvals],[vargn],[min],[max],[cbar],[varu]]
2352	[ncfile]: netCDF to use to geolocalize the points
2353	[lonvarn]: name of the variable with the longitudes
2354	[latvarn]: name of the variable with the latitudes
2355	[varn]: optional variable to add staff into the graph
2356	[dimval]: '@' list of [dimn]\|[dimval] to get the values for [varn]
2357	[dimn]: name of the dimension
2358	[dimval]: value of the dimension (no value all range)
2359	[vargn]: name of the variable in the graph
2360	[min]: minimum value for the extra variable
2361	[max]: maximum value for the extra variable
2362	[cbar]: color bar
2363	[varu]: units of the variable
2364	"""
2365	fname = 'draw_points'
2366
2367	if values == 'h':
2368	print fname + '_____________________________________________________________'
2369	print draw_points.__doc__
2370	quit()
2371
2372	expectargs = '[ptasciifile]:[gtit]:[mapvalues]:[kindfigure]:[pointcolor]:' + \
2373	'[pointlabels]:[locleg]:[figurek]:[figuren]'
2374
2375	drw.check_arguments(fname,values,expectargs,':')
2376
2377	ptasciifile = values.split(':')[0]
2378	gtit = values.split(':')[1]
2379	mapvalues = values.split(':')[2]
2380	kindfigure = values.split(':')[3]
2381	pointcolor = values.split(':')[4]
2382	pointlabels = values.split(':')[5]
2383	locleg = int(values.split(':')[6])
2384	figureko = values.split(':')[7]
2385	figuren = values.split(':')[8]
2386
2387	# Getting station locations
2388	##
2389	filev = ptasciifile.split(',')[0]
2390	comchar = ptasciifile.split(',')[1].split('\|')
2391	collon = int(ptasciifile.split(',')[2])
2392	collat = int(ptasciifile.split(',')[3])
2393	collab = ptasciifile.split(',')[4]
2394
2395	if not os.path.isfile(filev):
2396	print errormsg
2397	print ' ' + fname + ": file '" + filev + "' does not exist!!"
2398	quit(-1)
2399
2400	# Getting points position and labels
2401	oascii = open(filev, 'r')
2402	xptval = []
2403	yptval = []
2404	if collab != 'None':
2405	ptlabels = []
2406	for line in oascii:
2407	if not drw.searchInlist(comchar, line[0:1]):
2408	linevals = drw.reduce_spaces(line)
2409	xptval.append(np.float(linevals[collon].replace('\n','')))
2410	yptval.append(np.float(linevals[collat].replace('\n','')))
2411	ptlabels.append(linevals[int(collab)].replace('\n',''))
2412	else:
2413	ptlabels = None
2414	for line in oascii:
2415	if not drw.searchInlist(comchar, line[0:1]):
2416	linevals = drw.reduce_spaces(line)
2417	xptval.append(np.float(linevals[collon].replace('\n','')))
2418	yptval.append(np.float(linevals[collat].replace('\n','')))
2419
2420	oascii.close()
2421
2422	if pointlabels != 'None' and collab == 'None':
2423	ptlabels = pointlabels.split(',')
2424
2425	# Getting localization of the points
2426	##
2427	filev = filen.split(',')
2428	Nvals = len(filev)
2429
2430	ncfile = filev[0]
2431	lonvarn = filev[1]
2432	latvarn = filev[2]
2433	varn = None
2434	varextrav = None
2435	if Nvals == 10:
2436	varn = filev[3]
2437	dimvals = filev[4]
2438	varextrav = [filev[5], np.float(filev[6]), np.float(filev[7]), filev[8], \
2439	filev[9]]
2440
2441	if not os.path.isfile(ncfile):
2442	print errormsg
2443	print ' ' + fname + ": file '" + ncfile + "' does not exist!!"
2444	quit(-1)
2445
2446	onc = NetCDFFile(ncfile, 'r')
2447
2448	lonv, latv = drw.lonlat2D(onc.variables[lonvarn], onc.variables[latvarn])
2449
2450	if varn is not None:
2451	objextra = onc.variables[varn]
2452	vard = objextra.dimensions
2453	dd = {}
2454	for dn in dimvals.split('@'):
2455	ddn = dn.split('\|')[0]
2456	ddv = dn.split('\|')[1]
2457	dd[ddn] = ddv
2458	slicevar = []
2459	for dv in vard:
2460	found= False
2461	for dn in dd.keys():
2462	if dn == dv:
2463	slicevar.append(int(dd[dn]))
2464	found = True
2465	break
2466	if not found:
2467	slicevar.append(slice(0,len(onc.dimensions[dv])))
2468
2469	varextra = np.squeeze(objextra[tuple(slicevar)])
2470
2471	if mapvalues == 'None':
2472	mapV = None
2473	else:
2474	mapV = mapvalues
2475
2476	drw.plot_points(xptval, yptval, lonv, latv, varextra, varextrav, gtit, mapV, \
2477	kindfigure, pointcolor, ptlabels, locleg, figureko, figuren)
2478
2479	onc.close()
2480
2481	return
2482
2483	def draw_points_lonlat(filen, values):
2484	""" Function to plot a series of lon/lat points
2485	filen= name of the file
2486	values= [lonvarname]:[latvarname]:[gkind]:[gtit]:[ptcolor]:[pttype]:[ptsize]:[labels]:[locleg]:[figureK]
2487	[lonvarname]: name of the variable longitude
2488	[latvarname]: name of the variable latitude
2489	[gkind]: kind of graphical output
2490	[gtit]: graphic title '!' for spaces
2491	[ptcolor]: color of the points ('auto', for "red")
2492	[pttype]: type of point
2493	[ptsize]: size of point
2494	[labels]: ',' list of labels to use
2495	[locleg]: location of the legend (0, automatic)
2496	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2497	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2498	9: 'upper center', 10: 'center'
2499	[figureK]= kind of figure
2500	'legend': only points in the map with the legend with the names
2501	'labelled',[txtsize],[txtcol]: points with the names and size, color of text
2502	"""
2503	fname = 'draw_points_lonlat'
2504
2505	if values == 'h':
2506	print fname + '_____________________________________________________________'
2507	print draw_points_lonlat.__doc__
2508	quit()
2509
2510	expectargs = '[lonvarname]:[latvarname]:[gkind]:[gtit]:[ptcolor]:[pttype]:' + \
2511	'[ptsize]:[labels]:[locleg]:[figureK]'
2512
2513	drw.check_arguments(fname,values,expectargs,':')
2514
2515	lonname = values.split(':')[0]
2516	latname = values.split(':')[1]
2517	kindfigure = values.split(':')[2]
2518	gtit = values.split(':')[3].replace('!',' ')
2519	pointcolor = values.split(':')[4]
2520	pointtype = values.split(':')[5]
2521	pointsize = np.float(values.split(':')[6])
2522	labelsv = values.split(':')[7]
2523	loclegend = values.split(':')[8]
2524	figureK = values.split(':')[9]
2525
2526	fname = 'points_lonlat'
2527
2528	onc = NetCDFFile(filen, 'r')
2529	if not onc.variables.has_key(lonname):
2530	print errormsg
2531	print fname + ": file '" + filen + "' does not have longitudes '" + lonname +\
2532	"' !!"
2533	quit(-1)
2534	if not onc.variables.has_key(lonname):
2535	print errormsg
2536	print fname + ": file '" + filen + "' does not have longitudes '" + lonname +\
2537	"' !!"
2538	quit(-1)
2539
2540	olon = onc.variables[lonname]
2541	olat = onc.variables[latname]
2542
2543	Ndimlon = len(olon.shape)
2544	if Ndimlon == 1:
2545	dx = olon.shape[0]
2546	dy = olat.shape[0]
2547	if dx == dy:
2548	lonvals = olon[:]
2549	latvals = olat[:]
2550	else:
2551	lonvals0 = np.zeros((dy,dx), dtype=np.float)
2552	latvals0 = np.zeros((dy,dx), dtype=np.float)
2553	for iL in range(dy):
2554	lonvals0[iL,:] = olon[:]
2555	for il in range(dx):
2556	latvals0[:,il] = olat[:]
2557	lonvals = lonvals0.flatten()
2558	latvals = latvals0.flatten()
2559
2560	elif Ndimlon == 2:
2561	lonvals = olon[:].flatten()
2562	latvals = olat[:].flatten()
2563	elif Ndimlon == 3:
2564	lonvals = olon[1,:,:].flatten()
2565	latvals = olat[1,:,:].flatten()
2566	# Playing for Anna
2567	# lonvals = olon[:].flatten()
2568	# latvals = olat[:].flatten()
2569	elif Ndimlon == 4:
2570	lonvals = olon[1,0,:,:].flatten()
2571	latvals = olat[1,0,:,:].flatten()
2572	else:
2573	print errormsg
2574	print ' ' + fname + ': longitude size:',len(olon),' not ready!!'
2575	quit(-1)
2576
2577	if labelsv == 'None':
2578	labels = None
2579	else:
2580	labels = labelsv.split(',')
2581
2582	drw.plot_list_points(lonvals, latvals, lonname, latname, gtit, figureK, pointcolor, pointtype, \
2583	pointsize, labels, loclegend, kindfigure, fname)
2584
2585	onc.close()
2586
2587	return
2588
2589	def draw_timeSeries(filen, values, variables):
2590	""" Function to draw a time-series
2591	draw_timeSeries(filen, values, variable):
2592	filen= name of the file
2593	values= [gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[locleg]:[gkind]:[colorlines]:[pointtype]:[pointfreq]
2594	[gvarname]: name of the variable to appear in the graph
2595	[timetit]: title of the time axis (assumed x-axis, '\|' for spaces)
2596	[tkind]: kind of time to appear in the graph (assumed x-axis)
2597	'Nval': according to a given number of values as 'Nval',[Nval]
2598	'exct': according to an exact time unit as 'exct',[tunit];
2599	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2600	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2601	'l': milisecond
2602	[timefmt]: format of the time labels
2603	[title]: title of the graphic ('\|' for spaces)
2604	[locleg]: location of the legend (0, automatic)
2605	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2606	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2607	9: 'upper center', 10: 'center'
2608	[gkind]: kind of graphical output
2609	[colorlines]: ',' list of colors for the lines, None for automatic, single
2610	value all the same
2611	[pointtype]: ',' list of type of points for the lines, None for automatic, single
2612	value all the same
2613	[pointfreq]: frequency of point plotting, 'all' for all time steps
2614	variables= [varname],[timename] names of variable and variable with times
2615	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')
2616	"""
2617
2618	fname = 'draw_timeSeries'
2619
2620	if values == 'h':
2621	print fname + '_____________________________________________________________'
2622	print draw_timeSeries.__doc__
2623	quit()
2624
2625	expectargs = '[gvarname]:[timetit]:[tkind]:[timefmt]:[title]:' + \
2626	'[locleg]:[gkind]:[colorlines]:[pointtype]:[pointfreq]'
2627
2628	drw.check_arguments(fname,values,expectargs,':')
2629
2630	gvarname = values.split(':')[0]
2631	timetit = values.split(':')[1].replace('\|',' ')
2632	tkind = values.split(':')[2]
2633	timefmt = values.split(':')[3]
2634	title = values.split(':')[4].replace('\|',' ')
2635	locleg = int(values.split(':')[5])
2636	gkind = values.split(':')[6]
2637	colorlines = values.split(':')[7]
2638	pointtype = values.split(':')[8]
2639	pointfreq0 = values.split(':')[9]
2640
2641	ncobj = NetCDFFile(filen, 'r')
2642
2643	variable = variables.split(',')[0]
2644	timevar = variables.split(',')[1]
2645
2646	if not ncobj.variables.has_key(variable):
2647	print errormsg
2648	print ' ' + fname + ": file '" + filen + "' does not have variable '" + \
2649	variable + "' !!"
2650	quit(-1)
2651
2652	if not ncobj.variables.has_key(timevar):
2653	print errormsg
2654	print ' ' + fname + ": file '" + filen + "' does not have variable time '" \
2655	+ timevar + "' !!"
2656	quit(-1)
2657
2658	varobj = ncobj.variables[variable]
2659	timeobj = ncobj.variables[timevar]
2660
2661	dimt = len(timeobj[:])
2662	varvals = np.zeros((2,dimt), dtype=np.float)
2663
2664	gunits = varobj.getncattr('units')
2665	tunits = timeobj.getncattr('units')
2666
2667	varvals[0,:], valpot, newgunits, Spot = drw.pot_values(varobj[:].flatten(), gunits)
2668	varvals[1,:] = timeobj[:]
2669
2670	tseriesvals = []
2671	tseriesvals.append(varvals)
2672
2673	if colorlines == 'None':
2674	collines = None
2675	else:
2676	collines = colorlines.split(',')
2677	if pointtype == 'None':
2678	pttype = None
2679	else:
2680	pttype = pointtype.split(',')
2681
2682	if pointfreq0 == 'all':
2683	pointfreq = None
2684	else:
2685	pointfreq = int(pointfreq0)
2686
2687	drw.plot_TimeSeries(tseriesvals, Spot + drw.units_lunits(gunits), tunits, \
2688	'TimeSeries', gvarname, timetit, tkind, timefmt, title, \
2689	gvarname.replace('_','\_'), locleg, gkind, collines, pttype, pointfreq)
2690
2691	return
2692
2693	#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')
2694
2695	def draw_trajectories(trjfilens, values, observations):
2696	""" Function to draw different trajectories at the same time
2697	draw_trajectories(trjfilens, values, observations):
2698	trjfilens= [filen]@[Tint]@[map] ',' separated list of files with trajectories,
2699	time intervals and reference maps (first one will be used to plot)
2700	[filen]: name of the file to use (lines with '#', not readed) as:
2701	[t-step] [x] [y]
2702	[Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2703	[map]: [file]#[lonname]#[latname]
2704	[file]; with the [lon],[lat] matrices
2705	[lonname],[latname]; names of the longitudes and latitudes variables
2706	values=[leglabels]\|[lonlatlims]\|[title]\|[graphk]\|[mapkind]
2707	[leglabels]: ',' separated list of names for the legend
2708	[lonlatlims]: ',' list of limits of the map [lonmin, latmin, lonmax, latmax] or None
2709	[title]: title of the plot ('!' for spaces)
2710	[graphk]: kind of the graphic
2711	[mapkind]: drawing coastaline ([proj],[res]) or None
2712	[proj]: projection
2713	* 'cyl', cilindric
2714	* 'lcc', lambert conformal
2715	[res]: resolution:
2716	* 'c', crude
2717	* 'l', low
2718	* 'i', intermediate
2719	* 'h', high
2720	* 'f', full
2721	obsevations= [obsfile],[obsname],[Tint],[null]
2722	[obsfile]: name fo the File with the observations as [t-step] [lat] [lon]
2723	[obsname]: name of the observations in the graph
2724	[Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2725	[null]: null value for the observed trajectory
2726	"""
2727
2728	fname = 'draw_trajectories'
2729
2730	if values == 'h':
2731	print fname + '_____________________________________________________________'
2732	print draw_trajectories.__doc__
2733	quit()
2734
2735	expectargs = '[leglabels]\|[lonlatlims]\|[title]\|[graphk]\|[mapkind]'
2736
2737	drw.check_arguments(fname,values,expectargs,'\|')
2738
2739	trjfiles = trjfilens.split(',')
2740	leglabels = values.split('\|')[0]
2741	lonlatlims = values.split('\|')[1]
2742	title = values.split('\|')[2].replace('!',' ')
2743	graphk = values.split('\|')[3]
2744	mapkind = values.split('\|')[4]
2745
2746	Nfiles = len(trjfiles)
2747
2748	# Getting trajectotries
2749	##
2750
2751	lontrjvalues = []
2752	lattrjvalues = []
2753
2754	print ' ' + fname
2755	ifn = 0
2756	for ifile in trjfiles:
2757	filen = ifile.split('@')[0]
2758	Tint = ifile.split('@')[1]
2759
2760	print ' trajectory:',filen
2761
2762	if Tint != '-1':
2763	Tbeg = Tint
2764	Tend = ifile.split('@')[2]
2765	mapv = ifile.split('@')[3]
2766	else:
2767	mapv = ifile.split('@')[2]
2768
2769	if not os.path.isfile(filen):
2770	print errormsg
2771	print ' ' + fname + ": trajectory file '" + filen + "' does not exist !!"
2772	quit(-1)
2773
2774	# Charging longitude and latitude values
2775	##
2776	lonvals, latvals = drw.lonlat_values(mapv.split('#')[0], mapv.split('#')[1], \
2777	mapv.split('#')[2])
2778
2779	if ifn == 0: mapref = mapv
2780	ifn = ifn + 1
2781
2782	objfile = open(filen, 'r')
2783	trjtimev = []
2784	trjxv = []
2785	trjyv = []
2786
2787	for line in objfile:
2788	if line[0:1] != '#':
2789	trjtimev.append(int(line.split(' ')[0]))
2790	trjxv.append(int(line.split(' ')[1]))
2791	trjyv.append(int(line.split(' ')[2]))
2792
2793	objfile.close()
2794
2795	if Tint != '-1':
2796	lontrjvalues.append(lonvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
2797	lattrjvalues.append(latvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
2798	else:
2799	lontrjvalues.append(lonvals[trjyv[:], trjxv[:]])
2800	lattrjvalues.append(latvals[trjyv[:], trjxv[:]])
2801
2802	# lonlatlimits
2803	##
2804
2805	if lonlatlims == 'None':
2806	lonlatlimsv = None
2807	else:
2808	lonlatlimsv = np.zeros((4), dtype=np.float)
2809	lonlatlimsv[0] = np.float(lonlatlims.split(',')[0])
2810	lonlatlimsv[1] = np.float(lonlatlims.split(',')[1])
2811	lonlatlimsv[2] = np.float(lonlatlims.split(',')[2])
2812	lonlatlimsv[3] = np.float(lonlatlims.split(',')[3])
2813
2814	# lon/lat objects
2815	##
2816	objnc = NetCDFFile(mapref.split('#')[0])
2817	lonobj = objnc.variables[mapref.split('#')[1]]
2818	latobj = objnc.variables[mapref.split('#')[2]]
2819
2820	# map
2821	##
2822	if mapkind == 'None':
2823	mapkindv = None
2824	else:
2825	mapkindv = mapkind
2826
2827	if observations is None:
2828	obsname = None
2829	else:
2830	obsfile = observations.split(',')[0]
2831	obsname = observations.split(',')[1]
2832	Tint = observations.split(',')[2]
2833	null = np.float(observations.split(',')[3])
2834	print ' observational trajectory:',obsfile
2835
2836	if not os.path.isfile(obsfile):
2837	print errormsg
2838	print ' ' + fname + ": observations trajectory file '" + obsfile + \
2839	"' does not exist !!"
2840	quit(-1)
2841
2842	objfile = open(obsfile, 'r')
2843	obstrjtimev = []
2844	obstrjxv = []
2845	obstrjyv = []
2846
2847	for line in objfile:
2848	if line[0:1] != '#':
2849	lon = np.float(line.split(' ')[2])
2850	lat = np.float(line.split(' ')[1])
2851	if not lon == null and not lat == null:
2852	obstrjtimev.append(int(line.split(' ')[0]))
2853	obstrjxv.append(lon)
2854	obstrjyv.append(lat)
2855	else:
2856	obstrjtimev.append(int(line.split(' ')[0]))
2857	obstrjxv.append(None)
2858	obstrjyv.append(None)
2859
2860	objfile.close()
2861
2862	if Tint != '-1':
2863	Tint = int(observations.split(',')[2].split('@')[0])
2864	Tend = int(observations.split(',')[2].split('@')[1])
2865	lontrjvalues.append(obstrjxv[Tint:Tend+1])
2866	lattrjvalues.append(obstrjyv[Tint:Tend+1])
2867	else:
2868	lontrjvalues.append(obstrjxv[:])
2869	lattrjvalues.append(obstrjyv[:])
2870
2871	drw.plot_Trajectories(lontrjvalues, lattrjvalues, leglabels.split(','), \
2872	lonobj, latobj, lonlatlimsv, title, graphk, mapkindv, obsname)
2873
2874	objnc.close()
2875
2876	return
2877
2878	def draw_vals_trajectories(ncfile, values, variable):
2879	""" Function to draw values from the outputs from 'compute_tevolboxtraj'
2880	draw_vals_trajectories(ncfile, values, variable)
2881	ncfile= [ncfile] ',' list of files to use
2882	values= [statisticskind]:[Tint]:[labels]@[locleg]:[gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[gkind]
2883	[statisticskind]=[statistics][kind]
2884	[statistics]: which statistics to use, from: 'center', 'min', 'max', 'mean',
2885	'mean2', 'stdev'
2886	[kind]: 'box', 'circle' statistics taking the values from a box or a circle
2887	'trj': value following the trajectory
2888	[Tint]: [Tbeg]@[Tend] or None, interval of time to plot or -1 for all the times
2889	[labels]: ',' separated list of labels for the legend
2890	[locleg]: location of the legend (0, automatic)
2891	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2892	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2893	9: 'upper center', 10: 'center'
2894	[gvarname]: name of the variable to appear in the graph
2895	[timetit]: title of the time axis (assumed x-axis, '\|' for spaces)
2896	[tkind]: kind of time to appear in the graph (assumed x-axis)
2897	'Nval': according to a given number of values as 'Nval',[Nval]
2898	'exct': according to an exact time unit as 'exct',[tunit];
2899	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2900	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2901	'l': milisecond
2902	[timefmt]: format of the time labels
2903	[title]: title of the graphic ('\|' for spaces)
2904	[gkind]: kind of graphical output
2905	variable= variable to use
2906	"""
2907
2908	fname = 'draw_vals_trajectories'
2909
2910	if values == 'h':
2911	print fname + '_____________________________________________________________'
2912	print draw_vals_trajectories.__doc__
2913	quit()
2914
2915	sims = ncfile.split(',')
2916
2917	if len(values.split(':')) != 9:
2918	print errormsg
2919	print ' ' + fname + ': wrong number of values!', len(values.split(':')), \
2920	'given 9 needed!!'
2921	print ' ',values.split(':')
2922	quit(-1)
2923
2924	statisticskind = values.split(':')[0]
2925	Tint = values.split(':')[1]
2926	labels = values.split(':')[2]
2927	gvarname = values.split(':')[3]
2928	timetit = values.split(':')[4].replace('\|',' ')
2929	tkind = values.split(':')[5]
2930	timefmt = values.split(':')[6]
2931	title = values.split(':')[7].replace('\|',' ')
2932	gkind = values.split(':')[8]
2933
2934	leglabels = labels.split('@')[0].split(',')
2935	locleg = int(labels.split('@')[1])
2936
2937	Nsims = len(sims)
2938
2939	if Tint != '-1':
2940	tini = np.float(Tint.split('@')[0])
2941	tend = np.float(Tint.split('@')[1])
2942	else:
2943	tini = -1.
2944	tend = -1.
2945
2946	vartimetrjv = []
2947
2948	print ' ' + fname
2949	for trjfile in sims:
2950	print ' ' + trjfile + ' ...'
2951	if not os.path.isfile(trjfile):
2952	print errormsg
2953	print ' ' + fname + ": trajectory file: '" + trjfile + \
2954	"' does not exist !!"
2955	quit(-1)
2956
2957	trjobj = NetCDFFile(trjfile, 'r')
2958	otim = trjobj.variables['time']
2959	if not trjobj.variables.has_key(statisticskind + '_' + variable):
2960	print errormsg
2961	print ' ' + fname + ": file '" + trjfile + "' does not have variable '"+\
2962	statisticskind + '_' + variable + "' !!"
2963	quit(-1)
2964	ovar = trjobj.variables[statisticskind + '_' + variable]
2965	dimt = otim.shape[0]
2966
2967	if trjfile == sims[0]:
2968	gunits = ovar.getncattr('units')
2969	lname = ovar.getncattr('long_name')
2970	tunits = otim.getncattr('units')
2971
2972	if tini != -1:
2973	tiniid = -1
2974	tendid = -1
2975	for itv in range(dimt):
2976	if otim[itv] <= tini and otim[itv+1] >= tini: tiniid = itv
2977	if otim[itv] <= tend and otim[itv+1] >= tend: tendid = itv
2978
2979	if tiniid == -1 or tendid == -1:
2980	print errormsg
2981	print ' ' + main + ' time interval ', tini,',',tend,' not found: ', \
2982	tendid, ',', tiniid, ' !!'
2983	print ' data interval [',otim[0], otim[dimt-1],']'
2984	quit(-1)
2985	dimt = tendid - tiniid + 1
2986
2987	else:
2988	dimt = otim.shape[0]
2989
2990	valsv = np.zeros((2,dimt), dtype=np.float)
2991	# Checking for time consistency
2992	if otim.getncattr('units') != tunits:
2993	print warnmsg
2994	print ' ' + fname + ': different time units in the plot!!'
2995	newtimes = drw.coincident_CFtimes(otim[:], tunits, otim.getncattr('units'))
2996	else:
2997	newtimes = otim[:]
2998
2999	if tini == -1:
3000	valsv[1,:] = newtimes[:]
3001	valsv[0,:] = ovar[:]
3002	else:
3003	valsv[1,:] = newtimes[tiniid:tendid+1]
3004	valsv[0,:] = ovar[tiniid:tendid+1]
3005
3006	vartimetrjv.append(valsv)
3007	trjobj.close()
3008
3009	drw.plot_TimeSeries(vartimetrjv, drw.units_lunits(gunits), tunits, \
3010	'val_trajectories_' + statisticskind, gvarname, timetit, tkind, timefmt, title,\
3011	leglabels, locleg, gkind)
3012
3013	def variable_values(values):
3014	""" Function to give back values for a given variable
3015	values= [varname] name of the variable
3016	"""
3017
3018	fname = 'variable_values'
3019
3020	values = drw.variables_values(values)
3021
3022	print fname,'values:',values
3023	print fname,'variable_name:',values[0]
3024	print fname,'standard_name:',values[1]
3025	print fname,'min,max:',str(values[2]) + ',' + str(values[3])
3026	print fname,'long_name:',values[4]
3027	print fname,'units:',values[5]
3028	print fname,'shad_colors:',values[6]
3029	print fname,'all_values:',drw.numVector_String(values,',')
3030
3031	return
3032
3033	def draw_ptZvals(ncfile, values, variable):
3034	""" Function to plot a given list of points and values
3035	ncfile= netCDF file to use
3036	values= [fvname]:[XYvar]:[pointype]:[pointsize]:[graphlimits]:[nxtype]:
3037	[figuretitle]:[colorbar]:[mapvalue]:[kindfig]
3038	fvname: name of the variable in the graph
3039	XYvar: [lon],[lat] variable names
3040	ptype: type of the point
3041	ptsize: size of the point
3042	graphlimits: minLON,minLAT,maxLON,maxLAT limits of the graph 'None' for the full size
3043	nxtype: minimum and maximum type
3044	'auto': values taken from the extrems of the data
3045	[min],[max]: given minimum and maximum values
3046	figtitle: title of the figure
3047	cbar: color bar
3048	mapv: map characteristics: [proj],[res]
3049	see full documentation: http://matplotlib.org/basemap/
3050	[proj]: projection
3051	* 'cyl', cilindric
3052	* 'lcc', lambert-conformal
3053	[res]: resolution:
3054	* 'c', crude
3055	* 'l', low
3056	* 'i', intermediate
3057	* 'h', high
3058	* 'f', full
3059	kfig: kind of figure
3060	variable= name of the variable to plot
3061	"""
3062	fname = 'draw_ptZvals'
3063	import numpy.ma as ma
3064
3065	if values == 'h':
3066	print fname + '_____________________________________________________________'
3067	print draw_ptZvals.__doc__
3068	quit()
3069
3070	expectargs = '[fvname]:[XYvar]:[pointype]:[pointsize]:[graphlmits]:[nxtype]:' + \
3071	'[figuretit]:[colorbar]:[mapvalue]:[kindfig]'
3072
3073	drw.check_arguments(fname,values,expectargs,':')
3074
3075	fvname = values.split(':')[0]
3076	XYvar = values.split(':')[1]
3077	pointype = values.split(':')[2]
3078	pointsize = int(values.split(':')[3])
3079	graphlimits = values.split(':')[4]
3080	nxtype = values.split(':')[5]
3081	figuretitle = values.split(':')[6].replace('!',' ')
3082	colorbar = values.split(':')[7]
3083	mapvalue = values.split(':')[8]
3084	kindfig = values.split(':')[9]
3085
3086	onc = NetCDFFile(ncfile, 'r')
3087
3088	if not onc.variables.has_key(variable):
3089	print errormsg
3090	print ' ' + fname + ": file '" + ncfile + "' does not have variable '" + \
3091	variable + "' !!"
3092	quit(-1)
3093
3094	# points
3095	lonvarn = XYvar.split(',')[0]
3096	latvarn = XYvar.split(',')[1]
3097
3098	if not onc.variables.has_key(lonvarn):
3099	print errormsg
3100	print ' ' + fname + ": file '" + ncfile + "' does not have longitude " + \
3101	"variable '" + lonvarn + "' !!"
3102	quit(-1)
3103
3104	if not onc.variables.has_key(latvarn):
3105	print errormsg
3106	print ' ' + fname + ": file '" + ncfile + "' does not have latitude " + \
3107	"variable '" + latvarn + "' !!"
3108	quit(-1)
3109
3110	olonvar = onc.variables[lonvarn]
3111	olatvar = onc.variables[latvarn]
3112	ovarvar = onc.variables[variable]
3113
3114	Lpts = len(olonvar[:].flatten())
3115
3116	pointvalues = ma.masked_array(np.zeros((Lpts,3), dtype=np.float))
3117	pointvalues[:,0] = olonvar[:]
3118	pointvalues[:,1] = olatvar[:]
3119	pointvalues[:,2] = ovarvar[:]
3120
3121	varattrs = ovarvar.ncattrs()
3122	if drw.searchInlist(varattrs, 'units'):
3123	fvunits = ovarvar.getncattr('units')
3124	else:
3125	fvunits = drw.variables_values(variable)[5]
3126
3127	# map value
3128	if mapvalue == 'None': mapvalue = None
3129
3130	# Graph limits
3131	if graphlimits != 'None':
3132	graphlts = np.zeros((4), dtype=np.float)
3133	for il in range(4): graphlts[il] = np.float(graphlimits.split(',')[il])
3134	pointvalues[:,0] = ma.masked_outside(pointvalues[:,0], graphlts[0], \
3135	graphlts[2])
3136	pointvalues[:,1] = ma.masked_outside(pointvalues[:,1], graphlts[3], \
3137	graphlts[2])
3138
3139	# for ip in range(Lpts):
3140	# if pointvalues[ip,0] < graphlts[0] or pointvalues[ip,0] > graphlts[2] \
3141	# or pointvalues[ip,1] < graphlts[1] or pointvalues[ip,1] > graphlts[3]:
3142	# print ip,pointvalues[ip,0:2], graphlts
3143	# pointvalues[ip,2] = None
3144	else:
3145	graphlts = None
3146
3147	drw.plot_ptZvals(fvname,fvunits,pointvalues,pointype,pointsize,graphlts, nxtype, \
3148	figuretitle,colorbar,mapvalue,kindfig)
3149
3150	return
3151
3152	#draw_ptZvals('OBSnetcdf.nc', 'pracc:lon,lat:o:80:2,42,7,47,:values!of!values:Blues:cyl,l:pdf', 'pr')
3153
3154	def draw_vectors(ncfile, values, varns):
3155	""" Function to plot wind vectors
3156	values= [dimname]\|[vardimname]\|[value]:[vecvals]:[windlabs]:[mapvalues]:
3157	[gtit]:[kindfig]:[figuren]
3158	'X/Y/Z/T'\|[dimname]\|[vardimname]\|[value]: ',', list for each basic dimension '\|' separated of:
3159	[dimname]: name of the dimension in the file
3160	[vardimname]: name of the variable with the values for the dimension in the file
3161	[value]: which value of the given dimension (-1, all; [ibeg]@[iend], i-range beginning, end)
3162	No value takes all the range of the dimension
3163	[vecvals]= [frequency],[color],[length]
3164	[frequency]: [xfreq]@[yfreq] frequency of values allong each axis ('None', all grid points;
3165	'auto', computed automatically to have 20 vectors along each axis)
3166	[color]: color of the vectors
3167	'singlecol'@[colorn]: all the vectors same color ('auto': for 'red')
3168	'wind'@[colorbar]: color of the vectors according to wind speed sqrt(u^2+v^2) and given [colorbar]
3169	all vectors the same length
3170	'3rdvar'@[colorbar]@[varn]@[units]: color of the vectors according to a 3rd variable (to be added at -v) and given [colorbar]
3171	all vectors the same length
3172	[length]: length of the wind vectors ('auto', for 9)
3173	[windlabs]= [windname],[windunits]
3174	[windname]: name of the wind variable in the graph
3175	[windunits]: units of the wind variable in the graph ('None', for the value in the file)
3176	[mapvalues]= map characteristics: [proj],[res]
3177	see full documentation: http://matplotlib.org/basemap/
3178	[proj]: projection
3179	* 'cyl', cilindric
3180	* 'lcc', lambert conformal
3181	[res]: resolution:
3182	* 'c', crude
3183	* 'l', low
3184	* 'i', intermediate
3185	* 'h', high
3186	* 'f', full
3187	gtit= title of the graph ('\|', for spaces)
3188	kindfig= kind of figure
3189	figuren= name of the figure
3190	ncfile= file to use
3191	varns= [uwind],[ywind] ',' list of the name of the variables with the u-wind,y-wind component
3192	"""
3193	fname = 'draw_vectors'
3194
3195	if values == 'h':
3196	print fname + '_____________________________________________________________'
3197	print draw_vectors.__doc__
3198	quit()
3199
3200	expectargs = '[X/Y/Z/T]\|[dimname]\|[vardimname]\|[value]:[vecvals]:[windlabs]:' + \
3201	'[mapvalues]:[gtit]:[kindfig]:[figuren]'
3202
3203	drw.check_arguments(fname,values,expectargs,':')
3204
3205	dimvals = values.split(':')[0]
3206	vecvals = values.split(':')[1]
3207	windlabels = values.split(':')[2]
3208	mapvalues = values.split(':')[3]
3209	gtit = values.split(':')[4]
3210	kindfig = values.split(':')[5]
3211	figuren = values.split(':')[6]
3212
3213	of = NetCDFFile(ncfile,'r')
3214
3215	dims = {}
3216	for dimv in dimvals.split(','):
3217	dns = dimv.split('\|')
3218	dims[dns[0]] = [dns[1], dns[2], dns[3]]
3219
3220	varNs = []
3221	for dn in dims.keys():
3222	if dn == 'X':
3223	varNs.append(dims[dn][1])
3224	dimx = len(of.dimensions[dims[dn][0]])
3225	elif dn == 'Y':
3226	varNs.append(dims[dn][1])
3227	dimy = len(of.dimensions[dims[dn][0]])
3228
3229	ivar = 0
3230	for wvar in varns.split(','):
3231	if not drw.searchInlist(of.variables.keys(), wvar):
3232	print errormsg
3233	print ' ' + fname + ": file does not have variable '" + wvar + "' !!"
3234	quit(-1)
3235	if ivar == 0:
3236	varNs.append(wvar)
3237	else:
3238	varNs.append(wvar)
3239
3240	ivar = 0
3241	for varN in varNs:
3242	varslice = []
3243
3244	ovarN = of.variables[varN]
3245	vard = ovarN.dimensions
3246	for vdn in vard:
3247	found = False
3248	for dd in dims.keys():
3249	if dims[dd][0] == vdn:
3250	if dims[dd][2].find('@') != -1:
3251	rvals = dims[dd][2].split('@')
3252	varslice.append(slice(int(rvals[0]), int(rvals[1])))
3253	elif dims[dd][2] == '-1':
3254	varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3255	else:
3256	varslice.append(int(dims[dd][2]))
3257
3258	found = True
3259	break
3260	if not found:
3261	varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3262
3263	if varN == dims['X'][1]:
3264	lonvals0 = np.squeeze(ovarN[tuple(varslice)])
3265	elif varN == dims['Y'][1]:
3266	latvals0 = np.squeeze(ovarN[tuple(varslice)])
3267	elif ivar == 2:
3268	uwvals = np.squeeze(np.array(ovarN[tuple(varslice)]))
3269	elif ivar == 3:
3270	vwvals = np.squeeze(ovarN[tuple(varslice)])
3271
3272	ivar = ivar + 1
3273
3274	# print 'Final shapes:',lonvals0.shape,':',latvals0.shape,':',uwvals.shape,':',
3275	# vwvals.shape
3276
3277	if len(uwvals.shape) != 2 or len(vwvals.shape) != 2:
3278	print errormsg
3279	print ' ' + fname + ': wrong size of the wind fields! they must be ' + \
3280	'2-dimensional!'
3281	print ' u-winds shape:',uwvals.shape,'dims:',of.variables[varNs[2]]
3282	print ' v-winds shape:',vwvals.shape,'dims:',of.variables[varNs[3]]
3283	print ' provide more values for their dimensions!!'
3284	quit(-1)
3285
3286	if len(lonvals0.shape) == 1:
3287	lonvals, latvals = np.meshgrid(lonvals0, latvals0)
3288	else:
3289	lonvals = lonvals0
3290	latvals = latvals0
3291
3292	# Vector values
3293	if vecvals.split(',')[0] == 'None':
3294	freqv = None
3295	else:
3296	freqv = vecvals.split(',')[0]
3297	colorvals = vecvals.split(',')[1]
3298	coln = colorvals.split('@')[0]
3299	colv = colorvals.split('@')[1]
3300	if coln == 'singlecol':
3301	colorv = colv
3302	elif coln == 'wind':
3303	colorv = np.sqrt(uwvals2 + vwvals2)
3304	elif coln == '3rdvar':
3305	if len(varn.split(',')) != 3:
3306	print errormsg
3307	print ' ' + fname + ": color of vectors should be according to '" + \
3308	coln + "' but a third varibale is not provided !!"
3309	quit(-1)
3310	ocolvec = of.variables[varNs[4]]
3311	colorv = ocolvec[:]
3312	stdvn, lonvn, unitsvn = drw.var_3desc(ocolvec)
3313	colorvals = colorvals + '@' + stdvn + '@' + unitsvn
3314	else:
3315	print errormsg
3316	print ' ' + fname + ": color type '" + coln + "' not ready !!"
3317	quit(-1)
3318
3319	lengthv = vecvals.split(',')[2]
3320
3321	# Vector labels
3322	windname = windlabels.split(',')[0]
3323	windunits = windlabels.split(',')[1]
3324
3325	drw.plot_vector(lonvals, latvals, uwvals, vwvals, freqv, colorvals, colorv, \
3326	lengthv, windname, windunits, mapvalues, gtit, kindfig, figuren)
3327
3328	of.close()
3329
3330	return
3331
3332	def draw_basins(ncfile, values, varns):
3333	""" Function to plot wind basins
3334	values= [dimname]\|[vardimname]\|[value]:[vecvals]:[windlabs]:[mapvalues]:
3335	[gtit]:[kindfig]:[figuren]
3336	'X/Y/Z/T'\|[dimname]\|[vardimname]\|[value]: ',', list for each basic dimension '\|' separated of:
3337	[dimname]: name of the dimension in the file
3338	[vardimname]: name of the variable with the values for the dimension in the file
3339	[value]: which value of the given dimension (-1, all; [ibeg]@[iend], i-range beginning, end)
3340	No value takes all the range of the dimension
3341	[vecvals]= [frequency],[color],[length]
3342	[frequency]: [xfreq]@[yfreq] frequency of values allong each axis ('None', all grid points;
3343	'auto', computed automatically to have 20 vectors along each axis)
3344	[color]: [colorbar]@[varn]@[units]: color of the vectors according to a 3rd variable (to be added at -v) and given [colorbar]
3345	all vectors the same length
3346	[length]: length of the wind vectors ('auto', for 9)
3347	[windlabs]= [windname],[windunits]
3348	[windname]: name of the wind variable in the graph
3349	[windunits]: units of the wind variable in the graph ('None', for the value in the file)
3350	[mapvalues]= map characteristics: [proj],[res]
3351	see full documentation: http://matplotlib.org/basemap/
3352	[proj]: projection
3353	* 'cyl', cilindric
3354	* 'lcc', lambert conformal
3355	[res]: resolution:
3356	* 'c', crude
3357	* 'l', low
3358	* 'i', intermediate
3359	* 'h', high
3360	* 'f', full
3361	gtit= title of the graph ('\|', for spaces)
3362	kindfig= kind of figure
3363	figuren= name of the figure
3364	ncfile= file to use
3365	varns= [lon],[lat],[outflow],[basinID] ',' list of the name of the variables with the lon,lat, the outflow and the basin ID
3366	"""
3367	fname = 'draw_basins'
3368
3369	if values == 'h':
3370	print fname + '_____________________________________________________________'
3371	print draw_vectors.__doc__
3372	quit()
3373
3374	expectargs = '[X/Y/Z/T]\|[dimname]\|[vardimname]\|[value]:[vecvals]:[windlabs]:' + \
3375	'[mapvalues]:[gtit]:[kindfig]:[figuren]'
3376
3377	drw.check_arguments(fname,values,expectargs,':')
3378
3379	dimvals = values.split(':')[0]
3380	vecvals = values.split(':')[1]
3381	windlabels = values.split(':')[2]
3382	mapvalues = values.split(':')[3]
3383	gtit = values.split(':')[4]
3384	kindfig = values.split(':')[5]
3385	figuren = values.split(':')[6]
3386
3387	of = NetCDFFile(ncfile,'r')
3388
3389	dims = {}
3390	for dimv in dimvals.split(','):
3391	dns = dimv.split('\|')
3392	dims[dns[0]] = [dns[1], dns[2], dns[3]]
3393
3394	varNs = []
3395	for dn in dims.keys():
3396	if dn == 'X':
3397	varNs.append(dims[dn][1])
3398	dimx = len(of.dimensions[dims[dn][0]])
3399	elif dn == 'Y':
3400	varNs.append(dims[dn][1])
3401	dimy = len(of.dimensions[dims[dn][0]])
3402
3403	ivar = 0
3404	for wvar in varns.split(','):
3405	if not drw.searchInlist(of.variables.keys(), wvar):
3406	print errormsg
3407	print ' ' + fname + ": file does not have variable '" + wvar + "' !!"
3408	quit(-1)
3409	if ivar == 0:
3410	varNs.append(wvar)
3411	else:
3412	varNs.append(wvar)
3413
3414	ivar = 0
3415	for varN in varNs:
3416	varslice = []
3417
3418	ovarN = of.variables[varN]
3419	vard = ovarN.dimensions
3420	for vdn in vard:
3421	found = False
3422	for dd in dims.keys():
3423	if dims[dd][0] == vdn:
3424	if dims[dd][2].find('@') != -1:
3425	rvals = dims[dd][2].split('@')
3426	varslice.append(slice(int(rvals[0]), int(rvals[1])))
3427	elif dims[dd][2] == '-1':
3428	varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3429	else:
3430	varslice.append(int(dims[dd][2]))
3431
3432	found = True
3433	break
3434	if not found:
3435	varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3436
3437	if varN == dims['X'][1]:
3438	lonvals0 = np.squeeze(ovarN[tuple(varslice)])
3439	elif varN == dims['Y'][1]:
3440	latvals0 = np.squeeze(ovarN[tuple(varslice)])
3441
3442	ivar = ivar + 1
3443
3444	if len(lonvals0.shape) == 1:
3445	lonvals, latvals = np.meshgrid(lonvals0, latvals0)
3446	else:
3447	lonvals = lonvals0
3448	latvals = latvals0
3449
3450	# Vector values
3451	if vecvals.split(',')[0] == 'None':
3452	freqv = None
3453	else:
3454	freqv = vecvals.split(',')[0]
3455
3456	colorvals = vecvals.split(',')[1]
3457	if len(varn.split(',')) != 3:
3458	print errormsg
3459	print ' ' + fname + ": color of vectors should be according to '" + \
3460	coln + "' but a third varibale is not provided !!"
3461	quit(-1)
3462
3463	ocolvec = of.variables[varNs[3]]
3464	colorv = ocolvec[:]
3465	stdvn, lonvn, unitsvn = drw.var_3desc(ocolvec)
3466	colorvals = colorvals + '@' + stdvn + '@' + unitsvn
3467
3468	lengthv = vecvals.split(',')[2]
3469
3470	# Vector labels
3471	windname = windlabels.split(',')[0]
3472	windunits = windlabels.split(',')[1]
3473
3474	# Vector angles
3475	oflow = ofile.variables[varNs[2]]
3476	angle = (oflow[:] - 1)*np.pi/4
3477	xflow = np.where(oflow[:] < 9, np.float(lengthv)*np.sin(angle), 0.)
3478	yflow = np.where(oflow[:] < 9, np.float(lengthv)*np.cos(angle), 0.)
3479
3480	drw.plot_basins(lonvals, latvals, xflow, yflow, freqv, colorvals, colorv, \
3481	lengthv, windname, windunits, mapvalues, gtit, kindfig, figuren)
3482
3483	of.close()
3484
3485	return
3486
3487	def draw_river_desc(ncfile, values, riverns):
3488	""" Function to plot rivers' description from ORCHIDEE's routing scheme
3489	values= [dimname]\|[vardimname]\|[value]:[basinvals]:[upstreamvals]:[mapvalues]:
3490	[gtit]:[kindfig]:[legloc]:[figuren]
3491	'X/Y'\|[dimname]\|[vardimname]\|[value]: ',', list for each basic dimension '\|' separated of:
3492	[dimname]: name of the dimension in the file
3493	[vardimname]: name of the variable with the values for the dimension in the file
3494	[value]: which value of the given dimension (-1, all; [ibeg]@[iend], i-range beginning, end)
3495	No value takes all the range of the dimension
3496	[basinsvals]= [colorline]
3497	[basincolor]: ',' list of colors of the line to use to mark the basins contours (single value also possible)
3498	[upstreamvals]= [upstreamvarn],[colorbar]
3499	[upstreamcolor]: colorbar to use to plot the basins upstream values
3500	[mapvalues]= map characteristics: [proj],[res]
3501	see full documentation: http://matplotlib.org/basemap/
3502	[proj]: projection
3503	* 'cyl', cilindric
3504	* 'lcc', lambert conformal
3505	[res]: resolution:
3506	* 'c', crude
3507	* 'l', low
3508	* 'i', intermediate
3509	* 'h', high
3510	* 'f', full
3511	gtit= title of the graph ('\|', for spaces)
3512	kindfig= kind of figure
3513	legloc= location of the legend (0, automatic)
3514	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
3515	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
3516	9: 'upper center', 10: 'center' kfig= kind of figure
3517	figuren= name of the figure
3518	ncfile= file to use
3519	riverns= ',' list of the name of the rivers to plot
3520	"""
3521	import numpy.ma as ma
3522	fname = 'draw_river_desc'
3523
3524	if values == 'h':
3525	print fname + '_____________________________________________________________'
3526	print draw_river_desc.__doc__
3527	quit()
3528
3529	expectargs = '[X/Y/Z/T]\|[dimname]\|[vardimname]\|[value]:[basinvals]:' + \
3530	'[upstreamvals]:[mapvalues]:[gtit]:[kindfig]:[legloc]:[figuren]'
3531
3532	drw.check_arguments(fname,values,expectargs,':')
3533
3534	dimvals = values.split(':')[0]
3535	basinvals = values.split(':')[1]
3536	upstreamvals = values.split(':')[2]
3537	mapvals = values.split(':')[3]
3538	gtit = values.split(':')[4]
3539	kindfig = values.split(':')[5]
3540	legloc = int(values.split(':')[6])
3541	figuren = values.split(':')[7]
3542
3543	basincol = basinvals
3544	if basincol.find(',') != 1:
3545	basincolor = basincol.split(',')
3546	else:
3547	basincolor = [basincol]
3548
3549	upstreamcolor = upstreamvals
3550
3551	of = NetCDFFile(ncfile,'r')
3552
3553	dims = {}
3554	for dimv in dimvals.split(','):
3555	dns = dimv.split('\|')
3556	dims[dns[0]] = [dns[1], dns[2], dns[3]]
3557
3558	varNs = []
3559	for dn in dims.keys():
3560	if dn == 'X':
3561	varNs.append(dims[dn][1])
3562	dimx = len(of.dimensions[dims[dn][0]])
3563	elif dn == 'Y':
3564	varNs.append(dims[dn][1])
3565	dimy = len(of.dimensions[dims[dn][0]])
3566
3567	if riverns.find(',') != -1:
3568	riverns = riverns.split(',')
3569	else:
3570	riverns = [riverns]
3571
3572	rivers = []
3573	riversubbasins = {}
3574	riversupstream = {}
3575	riversoutflow = {}
3576	for rivern in riverns:
3577	print rivern
3578
3579	# subBasins
3580	basinvar = rivern + '_coding'
3581	if not drw.searchInlist(of.variables.keys(), basinvar):
3582	print errormsg
3583	print ' ' + fname + ": file does not have variable '" + basinvar + "' !!"
3584	quit(-1)
3585	rivers.append(rivern)
3586	obasin = of.variables[basinvar]
3587	riversubbasins[rivern] = obasin[:]
3588	if rivern == riverns[0]:
3589	finalmask = obasin[:].mask
3590	else:
3591	finalmask = finalmask * obasin[:].mask
3592
3593	# upstream
3594	upstreamvar = rivern + '_upstream'
3595	if not drw.searchInlist(of.variables.keys(), upstreamvar):
3596	print errormsg
3597	print ' ' + fname + ": file does not have variable '" + upstreamvar + "' !!"
3598	quit(-1)
3599	oupstream = of.variables[upstreamvar]
3600	riversupstream[rivern] = oupstream[:]
3601	if rivern == riverns[0]:
3602	uunits = oupstream.getncattr('units')
3603
3604	# River metadata
3605	fracvar = rivern + '_frac'
3606	if not drw.searchInlist(of.variables.keys(), fracvar):
3607	print errormsg
3608	print ' ' + fname + ": file does not have variable '" + fracvar + "' !!"
3609	quit(-1)
3610	ofrac = of.variables[fracvar]
3611	riversoutflow[rivern] = [ofrac.getncattr('Longitude_of_outflow_point'), \
3612	ofrac.getncattr('Latitude_of_outflow_point')]
3613
3614	ivar = 0
3615	for varN in varNs:
3616	varslice = []
3617
3618	ovarN = of.variables[varN]
3619	vard = ovarN.dimensions
3620	for vdn in vard:
3621	found = False
3622	for dd in dims.keys():
3623	if dims[dd][0] == vdn:
3624	if dims[dd][2].find('@') != -1:
3625	rvals = dims[dd][2].split('@')
3626	varslice.append(slice(int(rvals[0]), int(rvals[1])))
3627	elif dims[dd][2] == '-1':
3628	varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3629	else:
3630	varslice.append(int(dims[dd][2]))
3631
3632	found = True
3633	break
3634	if not found:
3635	varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3636
3637	if varN == dims['X'][1]:
3638	lonvals0 = np.squeeze(ovarN[tuple(varslice)])
3639	elif varN == dims['Y'][1]:
3640	latvals0 = np.squeeze(ovarN[tuple(varslice)])
3641
3642	ivar = ivar + 1
3643
3644	if len(lonvals0.shape) == 1:
3645	lonvals, latvals = np.meshgrid(lonvals0, latvals0)
3646	else:
3647	lonvals = lonvals0
3648	latvals = latvals0
3649
3650	# Masking only the lon,lat with rivers
3651	malonvals = ma.masked_array(lonvals, mask=finalmask)
3652	malatvals = ma.masked_array(latvals, mask=finalmask)
3653
3654	if mapvals == 'None':
3655	mapvalues = None
3656	else:
3657	mapvalues = mapvals
3658
3659	drw.plot_river_desc(malonvals, malatvals, rivers, riversubbasins, riversupstream, riversoutflow, \
3660	basincolor, upstreamcolor, uunits, mapvalues, gtit, kindfig, legloc, figuren)
3661
3662	of.close()
3663
3664	#quit()
3665
3666	####### ###### ##### #### ### ## #
3667
3668	ngraphics = "'" + drw.numVector_String(namegraphics, "', '") + "'"
3669
3670	### Options
3671	##string_operation="operation to make: " + '\n' + " out, output values -S inidim1,[inidim2,...]:enddim1,[enddim2,...]"
3672	string_operation="""operation to make:
3673	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]
3674	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]
3675	[ckind]:
3676	'cmap': as it gets from colorbar
3677	'fixc,[colname]': fixed color [colname], all stright lines
3678	'fixsignc,[colname]': fixed color [colname], >0 stright, <0 dashed line
3679	"""
3680
3681	#print string_operation
3682
3683	parser = OptionParser()
3684	parser.add_option("-f", "--netCDF_file", dest="ncfile",
3685	help="file to use", metavar="FILE")
3686	parser.add_option("-o", "--operation", type='choice', dest="operation",
3687	choices=namegraphics,
3688	help="operation to make: " + ngraphics, metavar="OPER")
3689	parser.add_option("-S", "--valueS", dest="values",
3690	help="[WHEN APPLICABLE] values to use according to the operation", metavar="VALUES")
3691	parser.add_option("-v", "--variable", dest="varname",
3692	help="[WHEN APPLICABLE] variable to check", metavar="VAR")
3693
3694	(opts, args) = parser.parse_args()
3695
3696	####### #######
3697	## MAIN
3698	#######
3699
3700	# Not checking file operation
3701	Notcheckingfile = ['draw_2D_shad_cont', 'draw_2D_shad_cont_time', \
3702	'draw_2D_shad_line', 'draw_2D_shad_line_time', 'draw_lines', 'draw_lines_time', \
3703	'draw_points', 'draw_topo_geogrid_boxes', 'draw_trajectories', \
3704	'draw_vals_trajectories', 'variable_values']
3705
3706	####### ###### ##### #### ### ## #
3707	errormsg='ERROR -- error -- ERROR -- error'
3708
3709	varn=opts.varname
3710	oper=opts.operation
3711
3712	if opts.ncfile is not None and not os.path.isfile(opts.ncfile) and \
3713	not drw.searchInlist(Notcheckingfile, oper):
3714	print errormsg
3715	print ' ' + main + ': File ' + opts.ncfile + ' does not exist !!'
3716	quit(-1)
3717
3718	if oper == 'create_movie':
3719	create_movie(opts.ncfile, opts.values, opts.varname)
3720	elif oper == 'draw_2D_shad':
3721	draw_2D_shad(opts.ncfile, opts.values, opts.varname)
3722	elif oper == 'draw_2D_shad_time':
3723	draw_2D_shad_time(opts.ncfile, opts.values, opts.varname)
3724	elif oper == 'draw_2D_shad_cont':
3725	draw_2D_shad_cont(opts.ncfile, opts.values, opts.varname)
3726	elif oper == 'draw_2D_shad_cont_time':
3727	draw_2D_shad_cont_time(opts.ncfile, opts.values, opts.varname)
3728	elif oper == 'draw_2D_shad_line':
3729	draw_2D_shad_line(opts.ncfile, opts.values, opts.varname)
3730	elif oper == 'draw_2D_shad_line_time':
3731	draw_2D_shad_line_time(opts.ncfile, opts.values, opts.varname)
3732	elif oper == 'draw_barbs':
3733	draw_barbs(opts.ncfile, opts.values, opts.varname)
3734	elif oper == 'draw_basins':
3735	draw_basins(opts.ncfile, opts.values, opts.varname)
3736	elif oper == 'draw_Neighbourghood_evol':
3737	draw_Neighbourghood_evol(opts.ncfile, opts.values, opts.varname)
3738	elif oper == 'draw_lines':
3739	draw_lines(opts.ncfile, opts.values, opts.varname)
3740	elif oper == 'draw_lines_time':
3741	draw_lines_time(opts.ncfile, opts.values, opts.varname)
3742	elif oper == 'draw_points':
3743	draw_points(opts.ncfile, opts.values)
3744	elif oper == 'draw_points_lonlat':
3745	draw_points_lonlat(opts.ncfile, opts.values)
3746	elif oper == 'draw_ptZvals':
3747	draw_ptZvals(opts.ncfile, opts.values, opts.varname)
3748	elif oper == 'draw_river_desc':
3749	draw_river_desc(opts.ncfile, opts.values, opts.varname)
3750	elif oper == 'draw_timeSeries':
3751	draw_timeSeries(opts.ncfile, opts.values, opts.varname)
3752	elif oper == 'draw_topo_geogrid':
3753	draw_topo_geogrid(opts.ncfile, opts.values)
3754	elif oper == 'draw_topo_geogrid_boxes':
3755	draw_topo_geogrid_boxes(opts.ncfile, opts.values)
3756	elif oper == 'draw_trajectories':
3757	draw_trajectories(opts.ncfile, opts.values, opts.varname)
3758	elif oper == 'draw_vals_trajectories':
3759	draw_vals_trajectories(opts.ncfile, opts.values, opts.varname)
3760	elif oper == 'draw_vectors':
3761	draw_vectors(opts.ncfile, opts.values, opts.varname)
3762	elif oper == 'list_graphics':
3763	# From: http://www.diveintopython.net/power_of_introspection/all_together.html
3764	import drawing as myself
3765	object = myself
3766	for opern in namegraphics:
3767	if opern != 'list_graphics':
3768	print opern + '_______ ______ _____ ____ ___ __ _'
3769	print getattr(object, opern).__doc__
3770	elif oper == 'variable_values':
3771	variable_values(opts.values)
3772	else:
3773	print errormsg
3774	print ' ' + main + ": the graphic '" + oper + "' is not ready !!"
3775	print errormsg
3776	quit()

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

Download in other formats: