Context Navigation

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

Last change on this file since 1017 was 984, checked in by lfita, 9 years ago
Fixing explanation of `freqv' in 'draw_lines'
File size: 159.8 KB

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

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

Download in other formats: