Context Navigation

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

Last change on this file since 836 was 833, checked in by lfita, 9 years ago
Adding `basinlinewidth': with of the line to draw the basin
File size: 158.5 KB

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

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

Download in other formats: