Context Navigation

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

Last change on this file since 798 was 792, checked in by lfita, 9 years ago
Adding flattening of the matrices to meet arrays size
File size: 153.4 KB

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

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

Download in other formats: