Context Navigation

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

Last change on this file since 955 was 955, checked in by lfita, 8 years ago
Working version of 'shadow_contour' with 'transform' 'transform' do not flip the locations of the labels of the axis
File size: 158.6 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]:[graphk]
1667	[dimvname]: ',' list of names of the variable with he values of the common dimension
1668	[valuesaxis]: which axis will be used for the values ('x', or 'y')
1669	[dimtit]: title for the common dimension
1670	[leglabels]: ',' separated list of names for the legend
1671	[vartit]: name of the variable in the graph
1672	[title]: title of the plot ('\|' for spaces)
1673	[locleg]: location of the legend (0, autmoatic)
1674	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
1675	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
1676	9: 'upper center', 10: 'center'
1677	[graphk]: kind of the graphic
1678	varname= variable to plot
1679	values= 'XLAT:x:latitude:32x32:$wss^{*}$:wss Taylor's turbulence term:pdf'
1680	"""
1681
1682	fname = 'draw_lines'
1683
1684	if values == 'h':
1685	print fname + '_____________________________________________________________'
1686	print draw_lines.__doc__
1687	quit()
1688
1689	expectargs = '[dimvname]:[valuesaxis]:[dimtit]:[leglabels]:[vtit]:[title]:[locleg]:[graphk]'
1690	drw.check_arguments(fname,values,expectargs,':')
1691
1692	ncfiles = ncfilens.split(',')
1693	dimvnames = values.split(':')[0]
1694	valuesaxis = values.split(':')[1]
1695	dimtit = values.split(':')[2]
1696	leglabels = values.split(':')[3].replace('_','\_')
1697	vartit = values.split(':')[4]
1698	title = values.split(':')[5].replace('\|',' ')
1699	locleg = values.split(':')[6]
1700	graphk = values.split(':')[7]
1701
1702	Nfiles = len(ncfiles)
1703
1704	# Getting trajectotries
1705	##
1706
1707	varvalues = []
1708	dimvalues = []
1709
1710	print ' ' + fname
1711	ifn = 0
1712	for ifile in ncfiles:
1713	filen = ifile.split('@')[0]
1714
1715	print ' filen:',filen
1716
1717	if not os.path.isfile(filen):
1718	print errormsg
1719	print ' ' + fname + ": netCDF file '" + filen + "' does not exist !!"
1720	quit(-1)
1721
1722	objfile = NetCDFFile(filen, 'r')
1723
1724	if dimvnames.find(',') != -1:
1725	dimvname = dimvnames.split(',')
1726	else:
1727	dimvname = [dimvnames]
1728
1729	found = False
1730	for dvn in dimvname:
1731	if objfile.variables.has_key(dvn):
1732	found = True
1733	break
1734
1735	if not found:
1736	print errormsg
1737	print ' ' + fname + ": netCDF file '" + filen + \
1738	"' does not have variables '" + dimvnames + "' !!"
1739	quit(-1)
1740
1741	if not objfile.variables.has_key(varname):
1742	print errormsg
1743	print ' ' + fname + ": netCDF file '" + filen + \
1744	"' does not have variable '" + varname + "' !!"
1745	quit(-1)
1746
1747	vvobj = objfile.variables[varname]
1748	if len(vvobj.shape) != 1:
1749	print errormsg
1750	print ' ' + fname + ': wrong shape:',vvobj.shape," of variable '" + \
1751	varname + "' !!"
1752	quit(-1)
1753
1754	for dimvn in dimvname:
1755	if drw.searchInlist(objfile.variables, dimvn):
1756	vdobj = objfile.variables[dimvn]
1757	if len(vdobj.shape) != 1:
1758	print errormsg
1759	print ' ' + fname + ': wrong shape:',vdobj.shape, \
1760	" of variable '" + dimvn + "' !!"
1761	quit(-1)
1762	break
1763
1764	varvalues.append(vvobj[:])
1765	dimvalues.append(vdobj[:])
1766
1767	if ifn == 0:
1768	varunits = vvobj.units
1769
1770	objfile.close()
1771
1772	ifn = ifn + 1
1773
1774	drw.plot_lines(dimvalues, varvalues, valuesaxis, dimtit, leglabels.split(','), \
1775	vartit, varunits, title, locleg, graphk)
1776
1777	return
1778
1779	def draw_lines_time(ncfilens, values, varname0):
1780	""" Function to draw different lines at the same time from different files with times
1781	draw_lines_time(ncfilens, values, varname):
1782	ncfilens= [filen] ',' separated list of netCDF files
1783	values= [dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];[rangevals];[timevals];
1784	[legvals];[graphk];[collines];[points];[linewidths];[pointsizes];[pointfreq];[period]
1785	[dimvname]: ',' list of names of the variables with he values of the common dimension
1786	[valuesaxis]: which axis will be used for the values ('x', or 'y')
1787	[dimtit]: title for the common dimension ('\|' for spaces)
1788	[leglabels]: ',' separated list of names for the legend ('None', no legend '!' for spaces)
1789	[vartit]: name of the variable in the graph
1790	[title]: title of the plot ('\|' for spaces)
1791	[rangevals]: Range of the axis with the values ('None' for 'auto','auto')
1792	[vmin],[vmax]: minimum and maximum values where [vmNN] can also be:
1793	'auto': the computed minimumm or maximum of the values
1794	[timevals]: [timen]\|[units]\|[kind]\|[tfmt] time labels characteristics
1795	[timen]; name of the time variable
1796	[units]; units string according to CF conventions ([tunits] since
1797	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
1798	[kind]; kind of output
1799	'Nval': according to a given number of values as 'Nval',[Nval]
1800	'exct': according to an exact time unit as 'exct',[tunit];
1801	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1802	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1803	'l': milisecond
1804	[tfmt]; desired format
1805	[legvals]=[locleg]\|[fontsize]:
1806	[locleg]: location of the legend (0, autmoatic)
1807	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
1808	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
1809	9: 'upper center', 10: 'center'
1810	[fontsize]: font size for the legend (auto for 12)
1811	[graphk]: kind of the graphic
1812	[lines]: ',' list of type of lines, None for automatic, single value all the same
1813	[collines]: ',' list of colors for the lines, None for automatic, single
1814	value all the same
1815	[points]: ',' list of type of points for the lines, None for automatic, single
1816	value all the same
1817	[linewidths]: ',' list of widths for the lines, None for automatic, single
1818	value all the same
1819	[pointsizes]: ',' list of widths for the lines, None for automatic, single
1820	value all the same
1821	[pointfreq]: frequency of point plotting, 'all' for all time steps
1822	[period]: which period to plot
1823	'-1': all period
1824	[beg],[end]: beginning and end of the period in reference time-units of first file
1825	varname0= ',' list of variable names to plot (assuming only 1 variable per file)
1826	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'
1827	"""
1828
1829	fname = 'draw_lines_time'
1830
1831	if values == 'h':
1832	print fname + '_____________________________________________________________'
1833	print draw_lines_time.__doc__
1834	quit()
1835
1836	expectargs = '[dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];'
1837	expectargs = expectargs + '[rangevals];[timevals];[legvals];[graphk];[lines];'
1838	expectargs = expectargs + '[collines];[points];[linewidths];[pointsizes];'
1839	expectargs = expectargs + '[pointfreq];[period]'
1840	drw.check_arguments(fname,values,expectargs,';')
1841
1842	ncfiles = ncfilens.split(',')
1843	dimvname0 = values.split(';')[0]
1844	valuesaxis = values.split(';')[1]
1845	dimtit = values.split(';')[2].replace('\|',' ')
1846	leglabels = values.split(';')[3].replace('_','\_').replace('!',' ')
1847	vartit = values.split(';')[4]
1848	title = values.split(';')[5].replace('\|',' ')
1849	rangevals = values.split(';')[6]
1850	timevals = values.split(';')[7]
1851	legvalues = values.split(';')[8]
1852	graphk = values.split(';')[9]
1853	lines0 = values.split(';')[10]
1854	collines0 = values.split(';')[11]
1855	points0 = values.split(';')[12]
1856	linewidths0 = values.split(';')[13]
1857	pointsizes0 = values.split(';')[14]
1858	pointfreq0 = values.split(';')[15]
1859	period = values.split(';')[16]
1860
1861	Nfiles = len(ncfiles)
1862
1863	# Multiple variable-dimension names?
1864	if dimvname0.find(',') != -1:
1865	dimvname = dimvname0.split(',')
1866	else:
1867	dimvname = [dimvname0]
1868
1869	# Multiple variables?
1870	if varname0.find(',') != -1:
1871	varname = varname0.split(',')
1872	else:
1873	varname = [varname0]
1874
1875	# Multiple lines types?
1876	if lines0.find(',') != -1:
1877	lines = lines0.split(',')
1878	elif lines0 == 'None':
1879	lines = None
1880	else:
1881	lines = []
1882	for il in range(Nfiles):
1883	lines.append(lines0)
1884
1885	# Multiple color names?
1886	if collines0.find(',') != -1:
1887	collines = collines0.split(',')
1888	elif collines0 == 'None':
1889	collines = None
1890	else:
1891	collines = []
1892	for ip in range(Nfiles):
1893	collines.append(collines0)
1894
1895	# Multiple point types?
1896	if points0.find(',') != -1:
1897	if len(points0) == 1:
1898	points = []
1899	for ip in range(Nfiles):
1900	points.append(points0)
1901	else:
1902	points = points0.split(',')
1903	elif points0 == 'None':
1904	points = None
1905	else:
1906	points = []
1907	for ip in range(Nfiles):
1908	points.append(points0)
1909
1910	# Multiple line sizes?
1911	if linewidths0.find(',') != -1:
1912	linewidths = []
1913	Nlines = len(linewidths0.split(','))
1914	for il in range(Nlines):
1915	linewidths.append(np.float(linewidths0.split(',')[il]))
1916	elif linewidths0 == 'None':
1917	linewidths = None
1918	else:
1919	linewidths = [np.float(linewidths0)]
1920
1921	# Multiple point sizes?
1922	if pointsizes0.find(',') != -1:
1923	pointsizes = []
1924	Npts = len(pointsizes0.split(','))
1925	for ip in Npts:
1926	pointsizes.append(np.float(pointsizes0.split(',')[ip]))
1927	elif pointsizes0 == 'None':
1928	pointsizes = None
1929	else:
1930	pointsizes = [np.float(pointsizes0)]
1931
1932	timename = timevals.split('\|')[0]
1933	timeunit = timevals.split('\|')[1].replace('!',' ')
1934	timekind = timevals.split('\|')[2]
1935	timefmt = timevals.split('\|')[3]
1936
1937	if rangevals == 'None':
1938	valmin = 'auto'
1939	valmax = 'auto'
1940	else:
1941	valmin = rangevals.split(',')[0]
1942	valmax = rangevals.split(',')[1]
1943	if valmin != 'auto': valmin = np.float(valmin)
1944	if valmax != 'auto': valmax = np.float(valmax)
1945
1946	locleg = int(legvalues.split('\|')[0])
1947	if legvalues.split('\|')[1] == 'auto':
1948	legfontsize = 12
1949	else:
1950	legfontsize = int(legvalues.split('\|')[1])
1951
1952	# Getting values
1953	##
1954	varvalues = []
1955	dimvalues = []
1956	timvalues = []
1957	timvals0 = timvalues
1958
1959	print ' ' + fname
1960	ifn = 0
1961	mintval = 1.e20
1962	maxtval = -1.e20
1963
1964	for ifile in ncfiles:
1965	filen = ifile.split('@')[0]
1966
1967	print ' filen:',filen
1968
1969	if not os.path.isfile(filen):
1970	print errormsg
1971	print ' ' + fname + ": netCDF file '" + filen + "' does not exist !!"
1972	quit(-1)
1973
1974	objfile = NetCDFFile(filen, 'r')
1975
1976	founddvar = False
1977	for dvar in dimvname:
1978	if objfile.variables.has_key(dvar):
1979	founddvar = True
1980	vdobj = objfile.variables[dvar]
1981	if len(vdobj.shape) != 1:
1982	print errormsg
1983	print ' ' + fname + ': wrong shape:',vdobj.shape," of " + \
1984	"variable '" + dvar + "' !!"
1985	quit(-1)
1986	break
1987	if not founddvar:
1988	print errormsg
1989	print ' ' + fname + ": netCDF file '" + filen + \
1990	"' has any variable '", dimvname, "' !!"
1991	quit(-1)
1992
1993	foundvar = False
1994	for var in varname:
1995	if objfile.variables.has_key(var):
1996	foundvar = True
1997	vvobj = objfile.variables[var]
1998	if len(vvobj.shape) != 1:
1999	print errormsg
2000	print ' ' + fname + ': wrong shape:',vvobj.shape," of " + \
2001	"variable '" + var + "' !!"
2002	quit(-1)
2003
2004	break
2005	if not foundvar:
2006	print errormsg
2007	print ' ' + fname + ": netCDF file '" + filen + \
2008	"' has any variable '", varname, "' !!"
2009	quit(-1)
2010	if vdobj.units.find('month') != 1:
2011	print warnmsg
2012	print ' ' + fname + ": tranforming time units from 'months' to 'days'!!"
2013	timevals0, tunits0 = gen.CFmonthU_daysU(vdobj[:], vdobj.units)
2014	else:
2015	timevals0 = vdobj[:]
2016	tunits0 = str(vdobj.units)
2017
2018	# Getting period
2019	if ifn > 0:
2020	# Referring all times to the same reference time!
2021	reftvals = drw.coincident_CFtimes(timevals0, timeunit, tunits0)
2022	else:
2023	reftvals = timevals0
2024
2025	dimt = len(vdobj[:])
2026
2027	if period == '-1':
2028	varvalues.append(vvobj[:])
2029	dimvalues.append(reftvals)
2030	mindvals = np.min(reftvals)
2031	maxdvals = np.max(reftvals)
2032	else:
2033	ibeg=-1
2034	iend=-1
2035	tbeg = np.float(period.split(',')[0])
2036	tend = np.float(period.split(',')[1])
2037
2038	for it in range(dimt-1):
2039	if reftvals[it] <= tbeg and reftvals[it+1] > tbeg: ibeg = it
2040	if reftvals[it] <= tend and reftvals[it+1] > tend: iend = it + 1
2041	if ibeg != -1 and iend != -1: break
2042
2043	if ibeg == -1 and iend == -1:
2044	print warnmsg
2045	print ' ' + fname + ': Period:',tbeg,',',tend,'not found!!'
2046	print ' ibeg:',ibeg,'iend:',iend
2047	print ' period in file:',np.min(reftvals), np.max(reftvals)
2048	print ' getting all the period in file !!!'
2049	ibeg = 0
2050	iend = dimt
2051	elif iend == -1:
2052	iend = dimt
2053	print warnmsg
2054	print ' ' + fname + ': end of Period:',tbeg,',',tend,'not found!!'
2055	print ' getting last available time instead'
2056	print ' ibeg:',ibeg,'iend:',iend
2057	print ' period in file:',np.min(reftvals), np.max(reftvals)
2058	elif ibeg == -1:
2059	ibeg = 0
2060	print warnmsg
2061	print ' ' + fname + ': beginning of Period:',tbeg,',',tend, \
2062	'not found!!'
2063	print ' getting first available time instead'
2064	print ' ibeg:',ibeg,'iend:',iend
2065	print ' period in file:',np.min(reftvals), np.max(reftvals)
2066
2067	varvalues.append(vvobj[ibeg:iend])
2068	dimvalues.append(reftvals[ibeg:iend])
2069	mindvals = np.min(reftvals[ibeg:iend])
2070	maxdvals = np.max(reftvals[ibeg:iend])
2071
2072	dimt = iend - ibeg
2073
2074	if mindvals < mintval: mintval = mindvals
2075	if maxdvals > maxtval: maxtval = maxdvals
2076	print ' ' + fname + ": file '" + filen + "' period:", mindvals, '->', maxdvals
2077
2078	if ifn == 0:
2079	varunits = drw.units_lunits(vvobj.units)
2080
2081	objfile.close()
2082
2083	ifn = ifn + 1
2084
2085	# Times
2086
2087	dtvals = (maxtval - mintval)/dimt
2088	# dti = mintval-dtvals/2.
2089	# dte = maxtval+dtvals/2.
2090	dti = mintval
2091	dte = maxtval
2092	tvals = np.arange(dti, dte, dtvals)
2093
2094	dtiS = drw.datetimeStr_conversion(str(dti) + ',' + timeunit, 'cfTime', \
2095	'Y/m/d H-M-S')
2096	dteS = drw.datetimeStr_conversion(str(dte) + ',' + timeunit, 'cfTime', \
2097	'Y/m/d H-M-S')
2098
2099	print ' ' + fname + ': plotting from: ' + dtiS + ' to ' + dteS
2100
2101	timepos, timelabels = drw.CFtimes_plot(tvals, timeunit, timekind, timefmt)
2102
2103	# print 'Lluis min/max tval +/- dtval/2:', mintval-dtvals/2., maxtval+dtvals/2.,'dt:', len(tvals)
2104	# for it in range(len(timepos)):
2105	# print timepos[it], timelabels[it]
2106
2107	if leglabels != 'None':
2108	legvals = leglabels.split(',')
2109	else:
2110	legvals = None
2111
2112	if pointfreq0 == 'all':
2113	pointfreq = None
2114	else:
2115	pointfreq = int(pointfreq0)
2116
2117	drw.plot_lines_time(dimvalues, varvalues, valuesaxis, dimtit, legvals, vartit, \
2118	varunits, timepos, timelabels, title, locleg, legfontsize, graphk, valmin, \
2119	valmax, lines, collines, points, linewidths, pointsizes, pointfreq)
2120
2121	return
2122
2123	def draw_Neighbourghood_evol(filen, values, variable):
2124	""" Function to draw the temporal evolution of a neighbourghood around a point
2125	draw_Neighbourghood_evol(filen, values, variable)
2126	filen= netCDF file name
2127	values= [gvarname]:[dimsval]:[neigdims]:[Nneig]:[Ncol]:[timetits]:[tkinds]:
2128	[timefmts]:[gtitle]:[shadxtrms]:[cbar]:[gkind]:[ofile]
2129	[dimsval]: [dimn1]\|[val1]\|[dimv1],...,[dimnN]\|[valN]\|[dimvN] dimension names, values to get
2130	(-1, for all; no name/value pair given full length) and variable with values of the dimension
2131	NOTE: when dimsval[X,Y] == neigdims[X,Y], valX,valY --> valX,valY-Nneig/2, valX,valY+Nneig/2
2132	[neigdims]: [dimnX],[dimnY] dimensions mnames along which the neigbourghood should be defined
2133	[Nneig]: Number of grid points of the full side of the box (odd value)
2134	[Ncol]: Number of columns ('auto': square final plot)
2135	[gvarname]: name of the variable to appear in the graph
2136	[timetits]: [titX],[titY] titles of the axes ('\|' for spaces)
2137	[tkinds]: [tkindX]\|[tkindY] kinds of time to appear in the graph
2138	'Nval': according to a given number of values as 'Nval',[Nval]
2139	'exct': according to an exact time unit as 'exct',[tunit];
2140	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2141	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2142	'l': milisecond
2143	[timefmts]: [tfmtX],[tfmtY] format of the time labels
2144	[gtitle]: title of the graphic ('\|' for spaces)
2145	[shadxtrms]: Extremes for the shading
2146	[cbar]: colorbar to use
2147	[gkind]: kind of graphical output
2148	[ofile]: True/False whether the netcdf with data should be created or not
2149	variable= name of the variable
2150	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'
2151	"""
2152
2153	fname = 'draw_Neighbourghood_evol'
2154
2155	if values == 'h':
2156	print fname + '_____________________________________________________________'
2157	print draw_Neighbourghood_evol.__doc__
2158	quit()
2159
2160	expectargs = '[gvarname]:[dimsval]:[neigdims]:[Nneig]:[Ncol]:' + \
2161	'[timetits]:[tkinds]:[timefmts]:[gtitle]:[shadxtrms]:[cbar]:[gkind]:[ofile]'
2162
2163	drw.check_arguments(fname,values,expectargs,':')
2164
2165	gvarname = values.split(':')[0]
2166	dimsval = values.split(':')[1].split(',')
2167	neigdims = values.split(':')[2].split(',')
2168	Nneig = int(values.split(':')[3])
2169	Ncol0 = values.split(':')[4]
2170	timetits = values.split(':')[5].split(',')
2171	timekinds = values.split(':')[6].split('\|')
2172	timefmts = values.split(':')[7].split(',')
2173	gtitle = values.split(':')[8].replace('\|',' ')
2174	shadxtrms = values.split(':')[9].split(',')
2175	cbar = values.split(':')[10]
2176	gkind = values.split(':')[11]
2177	ofile = values.split(':')[12]
2178
2179	if Ncol0 != 'auto':
2180	Ncol = int(Ncol0)
2181	else:
2182	Ncol = Ncol0
2183
2184	timetits[0] = timetits[0].replace('\|',' ')
2185	timetits[1] = timetits[1].replace('\|',' ')
2186
2187	if np.mod(Nneig,2) == 0:
2188	print errormsg
2189	print ' ' + fname + ": an odd value for 'Nneig':",Nneig,'is required !!!'
2190	quit(-1)
2191
2192	Nneig2 = int(Nneig/2)
2193
2194	# Values to slice the variable
2195	dimvslice = {}
2196	dimvvalues = {}
2197	for dimvs in dimsval:
2198	dimn = dimvs.split('\|')[0]
2199	dimv = int(dimvs.split('\|')[1])
2200	dimnv = dimvs.split('\|')[2]
2201
2202	dimvvalues[dimn] = dimnv
2203	dimvslice[dimn] = dimv
2204
2205	ncobj = NetCDFFile(filen, 'r')
2206
2207	varobj = ncobj.variables[variable]
2208
2209	slicevar = []
2210	newdimn = []
2211	newdimsvar = {}
2212
2213	for dimn in varobj.dimensions:
2214	if not drw.searchInlist(dimvslice.keys(), dimn):
2215	dimsize = len(ncobj.dimensions[dimn])
2216	slicevar.append(slice(0, dimsize+1))
2217	newdimn.append(dimn)
2218	newdimsvar[dimn] = dimseize
2219
2220	for dimslicen in dimvslice.keys():
2221	if dimn == dimslicen:
2222	if dimvslice[dimn] != -1:
2223	if drw.searchInlist(neigdims, dimn):
2224	slicevar.append(slice(dimvslice[dimn]-Nneig2, \
2225	dimvslice[dimn]+Nneig2+1))
2226	newdimn.append(dimn)
2227	newdimsvar[dimn] = Nneig
2228	break
2229	else:
2230	slicevar.append(slice(dimvslice[dimn], dimvslice[dimn]+1))
2231	break
2232	else:
2233	dimsize = len(ncobj.dimensions[dimn])
2234	slicevar.append(slice(0, dimsize+1))
2235	newdimn.append(dimn)
2236	newdimsvar[dimn] = dimsize
2237	break
2238
2239	varv = varobj[tuple(slicevar)]
2240
2241	if len(newdimn) != 3:
2242	print errormsg
2243	print ' ' + fname + ': sliced variable with shape=', varv.shape, \
2244	' must have three dimensions',len(varv.shape),'given !!'
2245	quit(-1)
2246
2247	newdims = []
2248	for nwdims in newdimn:
2249	newdims.append(newdimsvar[nwdims])
2250
2251	# The dimension which is not in the neighbourhood dimensions must be time!
2252	for dim1 in newdimn:
2253	if not drw.searchInlist(neigdims, dim1):
2254	dimt = newdimsvar[dim1]
2255	dimtime = dim1
2256
2257	if Ncol == 'auto':
2258	dimtsqx = int(np.sqrt(dimt)) + 1
2259	dimtsqy = int(np.sqrt(dimt)) + 1
2260	else:
2261	dimtsqx = int(Ncol)
2262	dimtsqy = dimt/dimtsqx + 1
2263
2264	neighbourghood = np.ones((dimtsqyNneig,dimtsqxNneig), dtype=np.float)*fillValue
2265
2266	for it in range(dimt):
2267	ity = int(it/dimtsqx)
2268	itx = it-ity*dimtsqx
2269
2270	itty = (dimtsqy - ity - 1)*Nneig + Nneig2
2271	ittx = itx*Nneig + Nneig2
2272
2273	neighbourghood[itty-Nneig2:itty+Nneig2+1,ittx-Nneig2:ittx+Nneig2+1]= \
2274	varv[it,::-1,:]
2275
2276	variablevals = drw.variables_values(variable)
2277	if drw.searchInlist(varobj.ncattrs(), 'units'):
2278	vunits = varobj.units
2279	else:
2280	vunits = variablevals[5]
2281
2282	# Time values at the X/Y axes
2283	if ncobj.variables[dimvvalues[dimtime]].dtype == '\|S1':
2284	print ' ' + fname + ': WRF time variable!'
2285	refdate = '19491201000000'
2286	tunitsval = 'hours'
2287	dimtvalues = np.zeros((dimt), dtype=np.float)
2288	tvals = ncobj.variables[dimvvalues[dimtime]]
2289	yrref=refdate[0:4]
2290	monref=refdate[4:6]
2291	dayref=refdate[6:8]
2292	horref=refdate[8:10]
2293	minref=refdate[10:12]
2294	secref=refdate[12:14]
2295
2296	refdateS = yrref + '/' + monref + '/' + dayref + '_' + horref + ':' + \
2297	minref + ':' + secref
2298	tunits = tunitsval + ' since ' + refdateS
2299	for it in range(dimt):
2300	wrfdates = drw.datetimeStr_conversion(tvals[it,:],'WRFdatetime', 'matYmdHMS')
2301	dimtvalues[it] = drw.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
2302	else:
2303	dimtvalues = ncobj.variables[dimvvalues[dimtime]][:]
2304	tunits = ncobj.variables[newdimsvar[dimtime]].units
2305
2306	dimxv = dimtvalues[0:dimtsqx]
2307	dimyv = dimtvalues[0:dimt:dimtsqx]
2308
2309	dimn = ['time','time']
2310
2311	if ofile == 'True':
2312	ofilen = 'Neighbourghood_evol.nc'
2313	newnc = NetCDFFile(ofilen, 'w')
2314	# Dimensions
2315	newdim = newnc.createDimension('time',None)
2316	newdim = newnc.createDimension('y',dimtsqy*Nneig)
2317	newdim = newnc.createDimension('x',dimtsqx*Nneig)
2318	# Dimension values
2319	newvar = newnc.createVariable('time','f8',('time'))
2320	newvar[:] = np.arange(dimt)
2321	newattr = drw.basicvardef(newvar, 'time','time',tunits)
2322	# Neighbourhghood variable
2323	newvar = newnc.createVariable(variable + 'neigevol', 'f4', ('y','x'), \
2324	fill_value=fillValue)
2325	newvar[:] = neighbourghood
2326
2327	newnc.sync()
2328	newnc.close()
2329	print fname + ": Successfull generation of file '" + ofilen + "' !!"
2330
2331	# Time ticks
2332	timeposX, timelabelsX = drw.CFtimes_plot(dimxv, tunits, timekinds[0], timefmts[0])
2333	timeposY, timelabelsY = drw.CFtimes_plot(dimyv, tunits, timekinds[1], timefmts[1])
2334
2335	timepos = [timeposX[0:len(timeposX)], timeposY[len(timeposY):0:-1]]
2336	timelabels = [timelabelsX[0:len(timeposX)], timelabelsY[0:len(timeposY)]]
2337
2338	for i in range(2):
2339	if shadxtrms[i][0:1] != 'S':
2340	shadxtrms[i] = np.float(shadxtrms[i])
2341
2342	drw.plot_Neighbourghood_evol(neighbourghood, dimxv, dimyv, gvarname, timetits, \
2343	timepos, timelabels, cbar, Nneig, shadxtrms, vunits, gtitle, gkind, True)
2344
2345	def draw_points(filen, values):
2346	""" Function to plot a series of points
2347	[values]= [ptasciifile]:[gtit]:[mapvalues]:[kindfigure]:[pointcolor]:[pointlabels]:
2348	[locleg]:[figureko]:[figuren]
2349	[ptasciifile]:[file],[comchar],[collon],[collat],[lab]
2350	[file]: column ASCII file with the location of the points
2351	[comchar]: '\|' list of characters for commentaries
2352	[collon]: number of column with the longitude of the points
2353	[collat]: number of column with the latitude of the points
2354	[collab]: number of column with the labels of the points ('None', and will get
2355	the values from the [pointlabels] variable
2356	[gtit]: title of the figure ('\|' for spaces)
2357	[mapvalues]: drawing coastaline ([proj],[res]) or None
2358	[proj]: projection
2359	* 'cyl', cilindric
2360	* 'lcc', lambert conformal
2361	[res]: resolution:
2362	* 'c', crude
2363	* 'l', low
2364	* 'i', intermediate
2365	* 'h', high
2366	* 'f', full
2367	[kindfigure]: kind of figure
2368	'legend': only points in the map with the legend with the names
2369	'labelled',[txtsize],[txtcol]: points with the names and size, color of text
2370	[pointcolor]: color for the points ('auto' for "red")
2371	[pointlabels]: ',' of labels [only used if [collab]='None'] ('None' for no labels)
2372	[locleg]: location of the legend (0, autmoatic)
2373	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2374	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2375	9: 'upper center', 10: 'center'
2376	[figureko]: kind of the output file (pdf, png, ...)
2377	[figuren]: name of the figure
2378	[filen]= [ncfile],[lonvarn],[latvarn][,[varn],[dimvals],[vargn],[min],[max],[cbar],[varu]]
2379	[ncfile]: netCDF to use to geolocalize the points
2380	[lonvarn]: name of the variable with the longitudes
2381	[latvarn]: name of the variable with the latitudes
2382	[varn]: optional variable to add staff into the graph
2383	[dimval]: '@' list of [dimn]\|[dimval] to get the values for [varn]
2384	[dimn]: name of the dimension
2385	[dimval]: value of the dimension (no value all range)
2386	[vargn]: name of the variable in the graph
2387	[min]: minimum value for the extra variable
2388	[max]: maximum value for the extra variable
2389	[cbar]: color bar
2390	[varu]: units of the variable
2391	"""
2392	fname = 'draw_points'
2393
2394	if values == 'h':
2395	print fname + '_____________________________________________________________'
2396	print draw_points.__doc__
2397	quit()
2398
2399	expectargs = '[ptasciifile]:[gtit]:[mapvalues]:[kindfigure]:[pointcolor]:' + \
2400	'[pointlabels]:[locleg]:[figurek]:[figuren]'
2401
2402	drw.check_arguments(fname,values,expectargs,':')
2403
2404	ptasciifile = values.split(':')[0]
2405	gtit = values.split(':')[1]
2406	mapvalues = values.split(':')[2]
2407	kindfigure = values.split(':')[3]
2408	pointcolor = values.split(':')[4]
2409	pointlabels = values.split(':')[5]
2410	locleg = int(values.split(':')[6])
2411	figureko = values.split(':')[7]
2412	figuren = values.split(':')[8]
2413
2414	# Getting station locations
2415	##
2416	filev = ptasciifile.split(',')[0]
2417	comchar = ptasciifile.split(',')[1].split('\|')
2418	collon = int(ptasciifile.split(',')[2])
2419	collat = int(ptasciifile.split(',')[3])
2420	collab = ptasciifile.split(',')[4]
2421
2422	if not os.path.isfile(filev):
2423	print errormsg
2424	print ' ' + fname + ": file '" + filev + "' does not exist!!"
2425	quit(-1)
2426
2427	# Getting points position and labels
2428	oascii = open(filev, 'r')
2429	xptval = []
2430	yptval = []
2431	if collab != 'None':
2432	ptlabels = []
2433	for line in oascii:
2434	if not drw.searchInlist(comchar, line[0:1]):
2435	linevals = drw.reduce_spaces(line)
2436	xptval.append(np.float(linevals[collon].replace('\n','')))
2437	yptval.append(np.float(linevals[collat].replace('\n','')))
2438	ptlabels.append(linevals[int(collab)].replace('\n',''))
2439	else:
2440	ptlabels = None
2441	for line in oascii:
2442	if not drw.searchInlist(comchar, line[0:1]):
2443	linevals = drw.reduce_spaces(line)
2444	xptval.append(np.float(linevals[collon].replace('\n','')))
2445	yptval.append(np.float(linevals[collat].replace('\n','')))
2446
2447	oascii.close()
2448
2449	if pointlabels != 'None' and collab == 'None':
2450	ptlabels = pointlabels.split(',')
2451
2452	# Getting localization of the points
2453	##
2454	filev = filen.split(',')
2455	Nvals = len(filev)
2456
2457	ncfile = filev[0]
2458	lonvarn = filev[1]
2459	latvarn = filev[2]
2460	varn = None
2461	varextrav = None
2462	if Nvals == 10:
2463	varn = filev[3]
2464	dimvals = filev[4]
2465	varextrav = [filev[5], np.float(filev[6]), np.float(filev[7]), filev[8], \
2466	filev[9]]
2467
2468	if not os.path.isfile(ncfile):
2469	print errormsg
2470	print ' ' + fname + ": file '" + ncfile + "' does not exist!!"
2471	quit(-1)
2472
2473	onc = NetCDFFile(ncfile, 'r')
2474
2475	lonv, latv = drw.lonlat2D(onc.variables[lonvarn], onc.variables[latvarn])
2476
2477	if varn is not None:
2478	objextra = onc.variables[varn]
2479	vard = objextra.dimensions
2480	dd = {}
2481	for dn in dimvals.split('@'):
2482	ddn = dn.split('\|')[0]
2483	ddv = dn.split('\|')[1]
2484	dd[ddn] = ddv
2485	slicevar = []
2486	for dv in vard:
2487	found= False
2488	for dn in dd.keys():
2489	if dn == dv:
2490	slicevar.append(int(dd[dn]))
2491	found = True
2492	break
2493	if not found:
2494	slicevar.append(slice(0,len(onc.dimensions[dv])))
2495
2496	varextra = np.squeeze(objextra[tuple(slicevar)])
2497
2498	if mapvalues == 'None':
2499	mapV = None
2500	else:
2501	mapV = mapvalues
2502
2503	drw.plot_points(xptval, yptval, lonv, latv, varextra, varextrav, gtit, mapV, \
2504	kindfigure, pointcolor, ptlabels, locleg, figureko, figuren)
2505
2506	onc.close()
2507
2508	return
2509
2510	def draw_points_lonlat(filen, values):
2511	""" Function to plot a series of lon/lat points
2512	filen= name of the file
2513	values= [lonvarname]:[latvarname]:[gkind]:[gtit]:[ptcolor]:[pttype]:[ptsize]:[labels]:[locleg]:[figureK]
2514	[lonvarname]: name of the variable longitude
2515	[latvarname]: name of the variable latitude
2516	[gkind]: kind of graphical output
2517	[gtit]: graphic title '!' for spaces
2518	[ptcolor]: color of the points ('auto', for "red")
2519	[pttype]: type of point
2520	[ptsize]: size of point
2521	[labels]: ',' list of labels to use
2522	[locleg]: location of the legend (0, automatic)
2523	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2524	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2525	9: 'upper center', 10: 'center'
2526	[figureK]= kind of figure
2527	'legend': only points in the map with the legend with the names
2528	'labelled',[txtsize],[txtcol]: points with the names and size, color of text
2529	"""
2530	fname = 'draw_points_lonlat'
2531
2532	if values == 'h':
2533	print fname + '_____________________________________________________________'
2534	print draw_points_lonlat.__doc__
2535	quit()
2536
2537	expectargs = '[lonvarname]:[latvarname]:[gkind]:[gtit]:[ptcolor]:[pttype]:' + \
2538	'[ptsize]:[labels]:[locleg]:[figureK]'
2539
2540	drw.check_arguments(fname,values,expectargs,':')
2541
2542	lonname = values.split(':')[0]
2543	latname = values.split(':')[1]
2544	kindfigure = values.split(':')[2]
2545	gtit = values.split(':')[3].replace('!',' ')
2546	pointcolor = values.split(':')[4]
2547	pointtype = values.split(':')[5]
2548	pointsize = np.float(values.split(':')[6])
2549	labelsv = values.split(':')[7]
2550	loclegend = values.split(':')[8]
2551	figureK = values.split(':')[9]
2552
2553	fname = 'points_lonlat'
2554
2555	onc = NetCDFFile(filen, 'r')
2556	if not onc.variables.has_key(lonname):
2557	print errormsg
2558	print fname + ": file '" + filen + "' does not have longitudes '" + lonname +\
2559	"' !!"
2560	quit(-1)
2561	if not onc.variables.has_key(lonname):
2562	print errormsg
2563	print fname + ": file '" + filen + "' does not have longitudes '" + lonname +\
2564	"' !!"
2565	quit(-1)
2566
2567	olon = onc.variables[lonname]
2568	olat = onc.variables[latname]
2569
2570	Ndimlon = len(olon.shape)
2571	if Ndimlon == 1:
2572	dx = olon.shape[0]
2573	dy = olat.shape[0]
2574	if dx == dy:
2575	lonvals = olon[:]
2576	latvals = olat[:]
2577	else:
2578	lonvals0 = np.zeros((dy,dx), dtype=np.float)
2579	latvals0 = np.zeros((dy,dx), dtype=np.float)
2580	for iL in range(dy):
2581	lonvals0[iL,:] = olon[:]
2582	for il in range(dx):
2583	latvals0[:,il] = olat[:]
2584	lonvals = lonvals0.flatten()
2585	latvals = latvals0.flatten()
2586
2587	elif Ndimlon == 2:
2588	lonvals = olon[:].flatten()
2589	latvals = olat[:].flatten()
2590	elif Ndimlon == 3:
2591	lonvals = olon[1,:,:].flatten()
2592	latvals = olat[1,:,:].flatten()
2593	# Playing for Anna
2594	# lonvals = olon[:].flatten()
2595	# latvals = olat[:].flatten()
2596	elif Ndimlon == 4:
2597	lonvals = olon[1,0,:,:].flatten()
2598	latvals = olat[1,0,:,:].flatten()
2599	else:
2600	print errormsg
2601	print ' ' + fname + ': longitude size:',len(olon),' not ready!!'
2602	quit(-1)
2603
2604	if labelsv == 'None':
2605	labels = None
2606	else:
2607	labels = labelsv.split(',')
2608
2609	drw.plot_list_points(lonvals, latvals, lonname, latname, gtit, figureK, pointcolor, pointtype, \
2610	pointsize, labels, loclegend, kindfigure, fname)
2611
2612	onc.close()
2613
2614	return
2615
2616	def draw_timeSeries(filen, values, variables):
2617	""" Function to draw a time-series
2618	draw_timeSeries(filen, values, variable):
2619	filen= name of the file
2620	values= [gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[locleg]:[gkind]:[colorlines]:[pointtype]:[pointfreq]
2621	[gvarname]: name of the variable to appear in the graph
2622	[timetit]: title of the time axis (assumed x-axis, '\|' for spaces)
2623	[tkind]: kind of time to appear in the graph (assumed x-axis)
2624	'Nval': according to a given number of values as 'Nval',[Nval]
2625	'exct': according to an exact time unit as 'exct',[tunit];
2626	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2627	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2628	'l': milisecond
2629	[timefmt]: format of the time labels
2630	[title]: title of the graphic ('\|' for spaces)
2631	[locleg]: location of the legend (0, automatic)
2632	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2633	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2634	9: 'upper center', 10: 'center'
2635	[gkind]: kind of graphical output
2636	[colorlines]: ',' list of colors for the lines, None for automatic, single
2637	value all the same
2638	[pointtype]: ',' list of type of points for the lines, None for automatic, single
2639	value all the same
2640	[pointfreq]: frequency of point plotting, 'all' for all time steps
2641	variables= [varname],[timename] names of variable and variable with times
2642	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')
2643	"""
2644
2645	fname = 'draw_timeSeries'
2646
2647	if values == 'h':
2648	print fname + '_____________________________________________________________'
2649	print draw_timeSeries.__doc__
2650	quit()
2651
2652	expectargs = '[gvarname]:[timetit]:[tkind]:[timefmt]:[title]:' + \
2653	'[locleg]:[gkind]:[colorlines]:[pointtype]:[pointfreq]'
2654
2655	drw.check_arguments(fname,values,expectargs,':')
2656
2657	gvarname = values.split(':')[0]
2658	timetit = values.split(':')[1].replace('\|',' ')
2659	tkind = values.split(':')[2]
2660	timefmt = values.split(':')[3]
2661	title = values.split(':')[4].replace('\|',' ')
2662	locleg = int(values.split(':')[5])
2663	gkind = values.split(':')[6]
2664	colorlines = values.split(':')[7]
2665	pointtype = values.split(':')[8]
2666	pointfreq0 = values.split(':')[9]
2667
2668	ncobj = NetCDFFile(filen, 'r')
2669
2670	variable = variables.split(',')[0]
2671	timevar = variables.split(',')[1]
2672
2673	if not ncobj.variables.has_key(variable):
2674	print errormsg
2675	print ' ' + fname + ": file '" + filen + "' does not have variable '" + \
2676	variable + "' !!"
2677	quit(-1)
2678
2679	if not ncobj.variables.has_key(timevar):
2680	print errormsg
2681	print ' ' + fname + ": file '" + filen + "' does not have variable time '" \
2682	+ timevar + "' !!"
2683	quit(-1)
2684
2685	varobj = ncobj.variables[variable]
2686	timeobj = ncobj.variables[timevar]
2687
2688	dimt = len(timeobj[:])
2689	varvals = np.zeros((2,dimt), dtype=np.float)
2690
2691	gunits = varobj.getncattr('units')
2692	tunits = timeobj.getncattr('units')
2693
2694	varvals[0,:], valpot, newgunits, Spot = drw.pot_values(varobj[:].flatten(), gunits)
2695	varvals[1,:] = timeobj[:]
2696
2697	tseriesvals = []
2698	tseriesvals.append(varvals)
2699
2700	if colorlines == 'None':
2701	collines = None
2702	else:
2703	collines = colorlines.split(',')
2704	if pointtype == 'None':
2705	pttype = None
2706	else:
2707	pttype = pointtype.split(',')
2708
2709	if pointfreq0 == 'all':
2710	pointfreq = None
2711	else:
2712	pointfreq = int(pointfreq0)
2713
2714	drw.plot_TimeSeries(tseriesvals, Spot + drw.units_lunits(gunits), tunits, \
2715	'TimeSeries', gvarname, timetit, tkind, timefmt, title, \
2716	gvarname.replace('_','\_'), locleg, gkind, collines, pttype, pointfreq)
2717
2718	return
2719
2720	#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')
2721
2722	def draw_trajectories(trjfilens, values, observations):
2723	""" Function to draw different trajectories at the same time
2724	draw_trajectories(trjfilens, values, observations):
2725	trjfilens= [filen]@[Tint]@[map] ',' separated list of files with trajectories,
2726	time intervals and reference maps (first one will be used to plot)
2727	[filen]: name of the file to use (lines with '#', not readed) as:
2728	[t-step] [x] [y]
2729	[Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2730	[map]: [file]#[lonname]#[latname]
2731	[file]; with the [lon],[lat] matrices
2732	[lonname],[latname]; names of the longitudes and latitudes variables
2733	values=[leglabels]\|[lonlatlims]\|[title]\|[graphk]\|[mapkind]
2734	[leglabels]: ',' separated list of names for the legend
2735	[lonlatlims]: ',' list of limits of the map [lonmin, latmin, lonmax, latmax] or None
2736	[title]: title of the plot ('!' for spaces)
2737	[graphk]: kind of the graphic
2738	[mapkind]: drawing coastaline ([proj],[res]) or None
2739	[proj]: projection
2740	* 'cyl', cilindric
2741	* 'lcc', lambert conformal
2742	[res]: resolution:
2743	* 'c', crude
2744	* 'l', low
2745	* 'i', intermediate
2746	* 'h', high
2747	* 'f', full
2748	obsevations= [obsfile],[obsname],[Tint],[null]
2749	[obsfile]: name fo the File with the observations as [t-step] [lat] [lon]
2750	[obsname]: name of the observations in the graph
2751	[Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2752	[null]: null value for the observed trajectory
2753	"""
2754
2755	fname = 'draw_trajectories'
2756
2757	if values == 'h':
2758	print fname + '_____________________________________________________________'
2759	print draw_trajectories.__doc__
2760	quit()
2761
2762	expectargs = '[leglabels]\|[lonlatlims]\|[title]\|[graphk]\|[mapkind]'
2763
2764	drw.check_arguments(fname,values,expectargs,'\|')
2765
2766	trjfiles = trjfilens.split(',')
2767	leglabels = values.split('\|')[0]
2768	lonlatlims = values.split('\|')[1]
2769	title = values.split('\|')[2].replace('!',' ')
2770	graphk = values.split('\|')[3]
2771	mapkind = values.split('\|')[4]
2772
2773	Nfiles = len(trjfiles)
2774
2775	# Getting trajectotries
2776	##
2777
2778	lontrjvalues = []
2779	lattrjvalues = []
2780
2781	print ' ' + fname
2782	ifn = 0
2783	for ifile in trjfiles:
2784	filen = ifile.split('@')[0]
2785	Tint = ifile.split('@')[1]
2786
2787	print ' trajectory:',filen
2788
2789	if Tint != '-1':
2790	Tbeg = Tint
2791	Tend = ifile.split('@')[2]
2792	mapv = ifile.split('@')[3]
2793	else:
2794	mapv = ifile.split('@')[2]
2795
2796	if not os.path.isfile(filen):
2797	print errormsg
2798	print ' ' + fname + ": trajectory file '" + filen + "' does not exist !!"
2799	quit(-1)
2800
2801	# Charging longitude and latitude values
2802	##
2803	lonvals, latvals = drw.lonlat_values(mapv.split('#')[0], mapv.split('#')[1], \
2804	mapv.split('#')[2])
2805
2806	if ifn == 0: mapref = mapv
2807	ifn = ifn + 1
2808
2809	objfile = open(filen, 'r')
2810	trjtimev = []
2811	trjxv = []
2812	trjyv = []
2813
2814	for line in objfile:
2815	if line[0:1] != '#':
2816	trjtimev.append(int(line.split(' ')[0]))
2817	trjxv.append(int(line.split(' ')[1]))
2818	trjyv.append(int(line.split(' ')[2]))
2819
2820	objfile.close()
2821
2822	if Tint != '-1':
2823	lontrjvalues.append(lonvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
2824	lattrjvalues.append(latvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
2825	else:
2826	lontrjvalues.append(lonvals[trjyv[:], trjxv[:]])
2827	lattrjvalues.append(latvals[trjyv[:], trjxv[:]])
2828
2829	# lonlatlimits
2830	##
2831
2832	if lonlatlims == 'None':
2833	lonlatlimsv = None
2834	else:
2835	lonlatlimsv = np.zeros((4), dtype=np.float)
2836	lonlatlimsv[0] = np.float(lonlatlims.split(',')[0])
2837	lonlatlimsv[1] = np.float(lonlatlims.split(',')[1])
2838	lonlatlimsv[2] = np.float(lonlatlims.split(',')[2])
2839	lonlatlimsv[3] = np.float(lonlatlims.split(',')[3])
2840
2841	# lon/lat objects
2842	##
2843	objnc = NetCDFFile(mapref.split('#')[0])
2844	lonobj = objnc.variables[mapref.split('#')[1]]
2845	latobj = objnc.variables[mapref.split('#')[2]]
2846
2847	# map
2848	##
2849	if mapkind == 'None':
2850	mapkindv = None
2851	else:
2852	mapkindv = mapkind
2853
2854	if observations is None:
2855	obsname = None
2856	else:
2857	obsfile = observations.split(',')[0]
2858	obsname = observations.split(',')[1]
2859	Tint = observations.split(',')[2]
2860	null = np.float(observations.split(',')[3])
2861	print ' observational trajectory:',obsfile
2862
2863	if not os.path.isfile(obsfile):
2864	print errormsg
2865	print ' ' + fname + ": observations trajectory file '" + obsfile + \
2866	"' does not exist !!"
2867	quit(-1)
2868
2869	objfile = open(obsfile, 'r')
2870	obstrjtimev = []
2871	obstrjxv = []
2872	obstrjyv = []
2873
2874	for line in objfile:
2875	if line[0:1] != '#':
2876	lon = np.float(line.split(' ')[2])
2877	lat = np.float(line.split(' ')[1])
2878	if not lon == null and not lat == null:
2879	obstrjtimev.append(int(line.split(' ')[0]))
2880	obstrjxv.append(lon)
2881	obstrjyv.append(lat)
2882	else:
2883	obstrjtimev.append(int(line.split(' ')[0]))
2884	obstrjxv.append(None)
2885	obstrjyv.append(None)
2886
2887	objfile.close()
2888
2889	if Tint != '-1':
2890	Tint = int(observations.split(',')[2].split('@')[0])
2891	Tend = int(observations.split(',')[2].split('@')[1])
2892	lontrjvalues.append(obstrjxv[Tint:Tend+1])
2893	lattrjvalues.append(obstrjyv[Tint:Tend+1])
2894	else:
2895	lontrjvalues.append(obstrjxv[:])
2896	lattrjvalues.append(obstrjyv[:])
2897
2898	drw.plot_Trajectories(lontrjvalues, lattrjvalues, leglabels.split(','), \
2899	lonobj, latobj, lonlatlimsv, title, graphk, mapkindv, obsname)
2900
2901	objnc.close()
2902
2903	return
2904
2905	def draw_vals_trajectories(ncfile, values, variable):
2906	""" Function to draw values from the outputs from 'compute_tevolboxtraj'
2907	draw_vals_trajectories(ncfile, values, variable)
2908	ncfile= [ncfile] ',' list of files to use
2909	values= [statisticskind]:[Tint]:[labels]@[locleg]:[gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[gkind]
2910	[statisticskind]=[statistics][kind]
2911	[statistics]: which statistics to use, from: 'center', 'min', 'max', 'mean',
2912	'mean2', 'stdev'
2913	[kind]: 'box', 'circle' statistics taking the values from a box or a circle
2914	'trj': value following the trajectory
2915	[Tint]: [Tbeg]@[Tend] or None, interval of time to plot or -1 for all the times
2916	[labels]: ',' separated list of labels for the legend
2917	[locleg]: location of the legend (0, automatic)
2918	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2919	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2920	9: 'upper center', 10: 'center'
2921	[gvarname]: name of the variable to appear in the graph
2922	[timetit]: title of the time axis (assumed x-axis, '\|' for spaces)
2923	[tkind]: kind of time to appear in the graph (assumed x-axis)
2924	'Nval': according to a given number of values as 'Nval',[Nval]
2925	'exct': according to an exact time unit as 'exct',[tunit];
2926	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2927	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2928	'l': milisecond
2929	[timefmt]: format of the time labels
2930	[title]: title of the graphic ('\|' for spaces)
2931	[gkind]: kind of graphical output
2932	variable= variable to use
2933	"""
2934
2935	fname = 'draw_vals_trajectories'
2936
2937	if values == 'h':
2938	print fname + '_____________________________________________________________'
2939	print draw_vals_trajectories.__doc__
2940	quit()
2941
2942	sims = ncfile.split(',')
2943
2944	if len(values.split(':')) != 9:
2945	print errormsg
2946	print ' ' + fname + ': wrong number of values!', len(values.split(':')), \
2947	'given 9 needed!!'
2948	print ' ',values.split(':')
2949	quit(-1)
2950
2951	statisticskind = values.split(':')[0]
2952	Tint = values.split(':')[1]
2953	labels = values.split(':')[2]
2954	gvarname = values.split(':')[3]
2955	timetit = values.split(':')[4].replace('\|',' ')
2956	tkind = values.split(':')[5]
2957	timefmt = values.split(':')[6]
2958	title = values.split(':')[7].replace('\|',' ')
2959	gkind = values.split(':')[8]
2960
2961	leglabels = labels.split('@')[0].split(',')
2962	locleg = int(labels.split('@')[1])
2963
2964	Nsims = len(sims)
2965
2966	if Tint != '-1':
2967	tini = np.float(Tint.split('@')[0])
2968	tend = np.float(Tint.split('@')[1])
2969	else:
2970	tini = -1.
2971	tend = -1.
2972
2973	vartimetrjv = []
2974
2975	print ' ' + fname
2976	for trjfile in sims:
2977	print ' ' + trjfile + ' ...'
2978	if not os.path.isfile(trjfile):
2979	print errormsg
2980	print ' ' + fname + ": trajectory file: '" + trjfile + \
2981	"' does not exist !!"
2982	quit(-1)
2983
2984	trjobj = NetCDFFile(trjfile, 'r')
2985	otim = trjobj.variables['time']
2986	if not trjobj.variables.has_key(statisticskind + '_' + variable):
2987	print errormsg
2988	print ' ' + fname + ": file '" + trjfile + "' does not have variable '"+\
2989	statisticskind + '_' + variable + "' !!"
2990	quit(-1)
2991	ovar = trjobj.variables[statisticskind + '_' + variable]
2992	dimt = otim.shape[0]
2993
2994	if trjfile == sims[0]:
2995	gunits = ovar.getncattr('units')
2996	lname = ovar.getncattr('long_name')
2997	tunits = otim.getncattr('units')
2998
2999	if tini != -1:
3000	tiniid = -1
3001	tendid = -1
3002	for itv in range(dimt):
3003	if otim[itv] <= tini and otim[itv+1] >= tini: tiniid = itv
3004	if otim[itv] <= tend and otim[itv+1] >= tend: tendid = itv
3005
3006	if tiniid == -1 or tendid == -1:
3007	print errormsg
3008	print ' ' + main + ' time interval ', tini,',',tend,' not found: ', \
3009	tendid, ',', tiniid, ' !!'
3010	print ' data interval [',otim[0], otim[dimt-1],']'
3011	quit(-1)
3012	dimt = tendid - tiniid + 1
3013
3014	else:
3015	dimt = otim.shape[0]
3016
3017	valsv = np.zeros((2,dimt), dtype=np.float)
3018	# Checking for time consistency
3019	if otim.getncattr('units') != tunits:
3020	print warnmsg
3021	print ' ' + fname + ': different time units in the plot!!'
3022	newtimes = drw.coincident_CFtimes(otim[:], tunits, otim.getncattr('units'))
3023	else:
3024	newtimes = otim[:]
3025
3026	if tini == -1:
3027	valsv[1,:] = newtimes[:]
3028	valsv[0,:] = ovar[:]
3029	else:
3030	valsv[1,:] = newtimes[tiniid:tendid+1]
3031	valsv[0,:] = ovar[tiniid:tendid+1]
3032
3033	vartimetrjv.append(valsv)
3034	trjobj.close()
3035
3036	drw.plot_TimeSeries(vartimetrjv, drw.units_lunits(gunits), tunits, \
3037	'val_trajectories_' + statisticskind, gvarname, timetit, tkind, timefmt, title,\
3038	leglabels, locleg, gkind)
3039
3040	def variable_values(values):
3041	""" Function to give back values for a given variable
3042	values= [varname] name of the variable
3043	"""
3044
3045	fname = 'variable_values'
3046
3047	values = drw.variables_values(values)
3048
3049	print fname,'values:',values
3050	print fname,'variable_name:',values[0]
3051	print fname,'standard_name:',values[1]
3052	print fname,'min,max:',str(values[2]) + ',' + str(values[3])
3053	print fname,'long_name:',values[4]
3054	print fname,'units:',values[5]
3055	print fname,'shad_colors:',values[6]
3056	print fname,'all_values:',drw.numVector_String(values,',')
3057
3058	return
3059
3060	def draw_ptZvals(ncfile, values, variable):
3061	""" Function to plot a given list of points and values
3062	ncfile= netCDF file to use
3063	values= [fvname]:[XYvar]:[pointype]:[pointsize]:[graphlimits]:[nxtype]:
3064	[figuretitle]:[colorbar]:[mapvalue]:[kindfig]
3065	fvname: name of the variable in the graph
3066	XYvar: [lon],[lat] variable names
3067	ptype: type of the point
3068	ptsize: size of the point
3069	graphlimits: minLON,minLAT,maxLON,maxLAT limits of the graph 'None' for the full size
3070	nxtype: minimum and maximum type
3071	'auto': values taken from the extrems of the data
3072	[min],[max]: given minimum and maximum values
3073	figtitle: title of the figure
3074	cbar: color bar
3075	mapv: map characteristics: [proj],[res]
3076	see full documentation: http://matplotlib.org/basemap/
3077	[proj]: projection
3078	* 'cyl', cilindric
3079	* 'lcc', lambert-conformal
3080	[res]: resolution:
3081	* 'c', crude
3082	* 'l', low
3083	* 'i', intermediate
3084	* 'h', high
3085	* 'f', full
3086	kfig: kind of figure
3087	variable= name of the variable to plot
3088	"""
3089	fname = 'draw_ptZvals'
3090	import numpy.ma as ma
3091
3092	if values == 'h':
3093	print fname + '_____________________________________________________________'
3094	print draw_ptZvals.__doc__
3095	quit()
3096
3097	expectargs = '[fvname]:[XYvar]:[pointype]:[pointsize]:[graphlmits]:[nxtype]:' + \
3098	'[figuretit]:[colorbar]:[mapvalue]:[kindfig]'
3099
3100	drw.check_arguments(fname,values,expectargs,':')
3101
3102	fvname = values.split(':')[0]
3103	XYvar = values.split(':')[1]
3104	pointype = values.split(':')[2]
3105	pointsize = int(values.split(':')[3])
3106	graphlimits = values.split(':')[4]
3107	nxtype = values.split(':')[5]
3108	figuretitle = values.split(':')[6].replace('!',' ')
3109	colorbar = values.split(':')[7]
3110	mapvalue = values.split(':')[8]
3111	kindfig = values.split(':')[9]
3112
3113	onc = NetCDFFile(ncfile, 'r')
3114
3115	if not onc.variables.has_key(variable):
3116	print errormsg
3117	print ' ' + fname + ": file '" + ncfile + "' does not have variable '" + \
3118	variable + "' !!"
3119	quit(-1)
3120
3121	# points
3122	lonvarn = XYvar.split(',')[0]
3123	latvarn = XYvar.split(',')[1]
3124
3125	if not onc.variables.has_key(lonvarn):
3126	print errormsg
3127	print ' ' + fname + ": file '" + ncfile + "' does not have longitude " + \
3128	"variable '" + lonvarn + "' !!"
3129	quit(-1)
3130
3131	if not onc.variables.has_key(latvarn):
3132	print errormsg
3133	print ' ' + fname + ": file '" + ncfile + "' does not have latitude " + \
3134	"variable '" + latvarn + "' !!"
3135	quit(-1)
3136
3137	olonvar = onc.variables[lonvarn]
3138	olatvar = onc.variables[latvarn]
3139	ovarvar = onc.variables[variable]
3140
3141	Lpts = len(olonvar[:].flatten())
3142
3143	pointvalues = ma.masked_array(np.zeros((Lpts,3), dtype=np.float))
3144	pointvalues[:,0] = olonvar[:].flatten()
3145	pointvalues[:,1] = olatvar[:].flatten()
3146	pointvalues[:,2] = ovarvar[:].flatten()
3147
3148	varattrs = ovarvar.ncattrs()
3149	if drw.searchInlist(varattrs, 'units'):
3150	fvunits = ovarvar.getncattr('units')
3151	else:
3152	fvunits = drw.variables_values(variable)[5]
3153
3154	# map value
3155	if mapvalue == 'None': mapvalue = None
3156
3157	# Graph limits
3158	if graphlimits != 'None':
3159	graphlts = np.zeros((4), dtype=np.float)
3160	for il in range(4): graphlts[il] = np.float(graphlimits.split(',')[il])
3161	pointvalues[:,0] = ma.masked_outside(pointvalues[:,0], graphlts[0], \
3162	graphlts[2])
3163	pointvalues[:,1] = ma.masked_outside(pointvalues[:,1], graphlts[3], \
3164	graphlts[2])
3165
3166	# for ip in range(Lpts):
3167	# if pointvalues[ip,0] < graphlts[0] or pointvalues[ip,0] > graphlts[2] \
3168	# or pointvalues[ip,1] < graphlts[1] or pointvalues[ip,1] > graphlts[3]:
3169	# print ip,pointvalues[ip,0:2], graphlts
3170	# pointvalues[ip,2] = None
3171	else:
3172	graphlts = None
3173
3174	drw.plot_ptZvals(fvname,fvunits,pointvalues,pointype,pointsize,graphlts, nxtype, \
3175	figuretitle,colorbar,mapvalue,kindfig)
3176
3177	return
3178
3179	#draw_ptZvals('OBSnetcdf.nc', 'pracc:lon,lat:o:80:2,42,7,47,:values!of!values:Blues:cyl,l:pdf', 'pr')
3180
3181	def draw_vectors(ncfile, values, varns):
3182	""" Function to plot wind vectors
3183	values= [dimname]\|[vardimname]\|[value]:[vecvals]:[windlabs]:[mapvalues]:
3184	[gtit]:[kindfig]:[figuren]
3185	'X/Y/Z/T'\|[dimname]\|[vardimname]\|[value]: ',', list for each basic dimension '\|' separated of:
3186	[dimname]: name of the dimension in the file
3187	[vardimname]: name of the variable with the values for the dimension in the file
3188	[value]: which value of the given dimension (-1, all; [ibeg]@[iend], i-range beginning, end)
3189	No value takes all the range of the dimension
3190	[vecvals]= [frequency],[color],[length]
3191	[frequency]: [xfreq]@[yfreq] frequency of values allong each axis ('None', all grid points;
3192	'auto', computed automatically to have 20 vectors along each axis)
3193	[color]: color of the vectors
3194	'singlecol'@[colorn]: all the vectors same color ('auto': for 'red')
3195	'wind'@[colorbar]: color of the vectors according to wind speed sqrt(u^2+v^2) and given [colorbar]
3196	all vectors the same length
3197	'3rdvar'@[colorbar]@[varn]@[units]: color of the vectors according to a 3rd variable (to be added at -v) and given [colorbar]
3198	all vectors the same length
3199	[length]: length of the wind vectors ('auto', for 9)
3200	[windlabs]= [windname],[windunits]
3201	[windname]: name of the wind variable in the graph
3202	[windunits]: units of the wind variable in the graph ('None', for the value in the file)
3203	[mapvalues]= map characteristics: [proj],[res]
3204	see full documentation: http://matplotlib.org/basemap/
3205	[proj]: projection
3206	* 'cyl', cilindric
3207	* 'lcc', lambert conformal
3208	[res]: resolution:
3209	* 'c', crude
3210	* 'l', low
3211	* 'i', intermediate
3212	* 'h', high
3213	* 'f', full
3214	gtit= title of the graph ('\|', for spaces)
3215	kindfig= kind of figure
3216	figuren= name of the figure
3217	ncfile= file to use
3218	varns= [uwind],[ywind] ',' list of the name of the variables with the u-wind,y-wind component
3219	"""
3220	fname = 'draw_vectors'
3221
3222	if values == 'h':
3223	print fname + '_____________________________________________________________'
3224	print draw_vectors.__doc__
3225	quit()
3226
3227	expectargs = '[X/Y/Z/T]\|[dimname]\|[vardimname]\|[value]:[vecvals]:[windlabs]:' + \
3228	'[mapvalues]:[gtit]:[kindfig]:[figuren]'
3229
3230	drw.check_arguments(fname,values,expectargs,':')
3231
3232	dimvals = values.split(':')[0]
3233	vecvals = values.split(':')[1]
3234	windlabels = values.split(':')[2]
3235	mapvalues = values.split(':')[3]
3236	gtit = values.split(':')[4]
3237	kindfig = values.split(':')[5]
3238	figuren = values.split(':')[6]
3239
3240	of = NetCDFFile(ncfile,'r')
3241
3242	dims = {}
3243	for dimv in dimvals.split(','):
3244	dns = dimv.split('\|')
3245	dims[dns[0]] = [dns[1], dns[2], dns[3]]
3246
3247	varNs = []
3248	for dn in dims.keys():
3249	if dn == 'X':
3250	varNs.append(dims[dn][1])
3251	dimx = len(of.dimensions[dims[dn][0]])
3252	elif dn == 'Y':
3253	varNs.append(dims[dn][1])
3254	dimy = len(of.dimensions[dims[dn][0]])
3255
3256	ivar = 0
3257	for wvar in varns.split(','):
3258	if not drw.searchInlist(of.variables.keys(), wvar):
3259	print errormsg
3260	print ' ' + fname + ": file does not have variable '" + wvar + "' !!"
3261	quit(-1)
3262	if ivar == 0:
3263	varNs.append(wvar)
3264	else:
3265	varNs.append(wvar)
3266
3267	ivar = 0
3268	for varN in varNs:
3269	varslice = []
3270
3271	ovarN = of.variables[varN]
3272	vard = ovarN.dimensions
3273	for vdn in vard:
3274	found = False
3275	for dd in dims.keys():
3276	if dims[dd][0] == vdn:
3277	if dims[dd][2].find('@') != -1:
3278	rvals = dims[dd][2].split('@')
3279	varslice.append(slice(int(rvals[0]), int(rvals[1])))
3280	elif dims[dd][2] == '-1':
3281	varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3282	else:
3283	varslice.append(int(dims[dd][2]))
3284
3285	found = True
3286	break
3287	if not found:
3288	varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3289
3290	if varN == dims['X'][1]:
3291	lonvals0 = np.squeeze(ovarN[tuple(varslice)])
3292	elif varN == dims['Y'][1]:
3293	latvals0 = np.squeeze(ovarN[tuple(varslice)])
3294	elif ivar == 2:
3295	uwvals = np.squeeze(np.array(ovarN[tuple(varslice)]))
3296	elif ivar == 3:
3297	vwvals = np.squeeze(ovarN[tuple(varslice)])
3298
3299	ivar = ivar + 1
3300
3301	# print 'Final shapes:',lonvals0.shape,':',latvals0.shape,':',uwvals.shape,':',
3302	# vwvals.shape
3303
3304	if len(uwvals.shape) != 2 or len(vwvals.shape) != 2:
3305	print errormsg
3306	print ' ' + fname + ': wrong size of the wind fields! they must be ' + \
3307	'2-dimensional!'
3308	print ' u-winds shape:',uwvals.shape,'dims:',of.variables[varNs[2]]
3309	print ' v-winds shape:',vwvals.shape,'dims:',of.variables[varNs[3]]
3310	print ' provide more values for their dimensions!!'
3311	quit(-1)
3312
3313	if len(lonvals0.shape) == 1:
3314	lonvals, latvals = np.meshgrid(lonvals0, latvals0)
3315	else:
3316	lonvals = lonvals0
3317	latvals = latvals0
3318
3319	# Vector values
3320	if vecvals.split(',')[0] == 'None':
3321	freqv = None
3322	else:
3323	freqv = vecvals.split(',')[0]
3324	colorvals = vecvals.split(',')[1]
3325	coln = colorvals.split('@')[0]
3326	colv = colorvals.split('@')[1]
3327	if coln == 'singlecol':
3328	colorv = colv
3329	elif coln == 'wind':
3330	colorv = np.sqrt(uwvals2 + vwvals2)
3331	elif coln == '3rdvar':
3332	if len(varn.split(',')) != 3:
3333	print errormsg
3334	print ' ' + fname + ": color of vectors should be according to '" + \
3335	coln + "' but a third varibale is not provided !!"
3336	quit(-1)
3337	ocolvec = of.variables[varNs[4]]
3338	colorv = ocolvec[:]
3339	stdvn, lonvn, unitsvn = drw.var_3desc(ocolvec)
3340	colorvals = colorvals + '@' + stdvn + '@' + unitsvn
3341	else:
3342	print errormsg
3343	print ' ' + fname + ": color type '" + coln + "' not ready !!"
3344	quit(-1)
3345
3346	lengthv = vecvals.split(',')[2]
3347
3348	# Vector labels
3349	windname = windlabels.split(',')[0]
3350	windunits = windlabels.split(',')[1]
3351
3352	drw.plot_vector(lonvals, latvals, uwvals, vwvals, freqv, colorvals, colorv, \
3353	lengthv, windname, windunits, mapvalues, gtit, kindfig, figuren)
3354
3355	of.close()
3356
3357	return
3358
3359	def draw_basins(ncfile, values, varns):
3360	""" Function to plot river basins with their discharge vector and basins id (from 'routing.nc')
3361	values= [lonlatbox]:[mapres]:[cbarname]:[xtrmbasin]:[mapdraw]:[veclength]:[freq]:
3362	[ifreq]:[plotcountry]:[basinidn]:[gtit]:[kindfig]:[figuren]
3363	[lonlatbox]= [lonSW],[lonNE],[latSW],[latNE] coordinates of the lon/lat box
3364	[mapres]= resolution of the mapping information
3365	[cbarname]= colorbar name for the colors
3366	[xtrmbasin]= [minbasin],[maxbasin] minimum and maximum basin numbers
3367	[mapdraw]= whether to draw the map (and project the data) or not ('True/False')
3368	[veclength]= length of the vectors of discharge at each grid cell
3369	[freq]= frequency of values allong each axis (None, all grid points;
3370	'auto', computed automatically to have 20 vectors along each axis)
3371	[plotcountry]= whether country lines should be plotted or not ('True/False')
3372	[plotbasinid]= whether id of the basins should be plotted or not ('True/False')
3373	[gtit]= title of the graph ('\|', for spaces)
3374	[kindfig]= kind of figure
3375	[figuren]= name of the figure
3376	ncfile= file to use
3377	"""
3378	fname = 'draw_basins'
3379
3380	if values == 'h':
3381	print fname + '_____________________________________________________________'
3382	print draw_vectors.__doc__
3383	quit()
3384
3385	expectargs = '[lonlatbox]:[mapres]:[cbarname]:[xtrmbasin]:[mapdraw]:' + \
3386	'[veclength]:[freq]:[plotcountry]:[basinidn]:[gtit]:[kindfig]:[figuren]'
3387
3388	drw.check_arguments(fname,values,expectargs,':')
3389
3390	varn='basins'
3391	lonname = 'nav_lon'
3392	latname = 'nav_lat'
3393	flowname = 'trip'
3394
3395	lonlims =[]
3396	latlims =[]
3397
3398	lonlims.append(np.float(values.split(':')[0].split(',')[0]))
3399	lonlims.append(np.float(values.split(':')[0].split(',')[1]))
3400	latlims.append(np.float(values.split(':')[0].split(',')[2]))
3401	latlims.append(np.float(values.split(':')[0].split(',')[3]))
3402	map_res = values.split(':')[1]
3403	cbarname = values.split(':')[2]
3404	vtit = 'basins'
3405	minbasin = np.int(values.split(':')[3].split(',')[0])
3406	maxbasin = np.int(values.split(':')[3].split(',')[1])
3407	mapdraw = gen.Str_Bool(values.split(':')[4])
3408	veclength = np.float(values.split(':')[5])
3409	freq0 = values.split(':')[6]
3410	plotcountry = gen.Str_Bool(values.split(':')[7])
3411	plotbasinid = gen.Str_Bool(values.split(':')[8])
3412	gtit = values.split(':')[9].replace('\|',' ')
3413	kindfig = values.split(':')[10]
3414	figuren = values.split(':')[11]
3415
3416	if freq0 == 'None': freq = None
3417
3418	ofile = NetCDFFile(ncfile, 'r')
3419
3420	obasins = ofile.variables[varn]
3421	olon = ofile.variables[lonname]
3422	olat = ofile.variables[latname]
3423	oflow = ofile.variables[flowname]
3424
3425	lons = olon[:]
3426	lats = olat[:]
3427
3428	lon, lat = drw.lonlat2D(lons, lats)
3429
3430	nlon = lonlims[0]
3431	xlon = lonlims[1]
3432	nlat = latlims[0]
3433	xlat = latlims[1]
3434
3435	imin, imax, jmin, jmax = gen.ijlonlat(lon, lat, nlon, xlon, nlat, xlat)
3436
3437	drw.plot_basins(lon[jmin:jmax,imin:imax], lat[jmin:jmax,imin:imax], \
3438	oflow[jmin:jmax,imin:imax], freq, cbarname+'@basin@-', \
3439	obasins[jmin:jmax,imin:imax], veclength, minbasin, maxbasin, 'outflow', '-', \
3440	'cyl,'+map_res, plotcountry, plotbasinid, gtit, kindfig, figuren)
3441
3442	ofile.close()
3443
3444	return
3445
3446	def draw_basinsold(ncfile, values, varns):
3447	""" Function to plot wind basins
3448	values= [dimname]\|[vardimname]\|[value]:[vecvals]:[windlabs]:[mapvalues]:
3449	[gtit]:[kindfig]:[figuren]
3450	'X/Y/Z/T'\|[dimname]\|[vardimname]\|[value]: ',', list for each basic dimension '\|' separated of:
3451	[dimname]: name of the dimension in the file
3452	[vardimname]: name of the variable with the values for the dimension in the file
3453	[value]: which value of the given dimension (-1, all; [ibeg]@[iend], i-range beginning, end)
3454	No value takes all the range of the dimension
3455	[vecvals]= [frequency],[color],[length]
3456	[frequency]: [xfreq]@[yfreq] frequency of values allong each axis ('None', all grid points;
3457	'auto', computed automatically to have 20 vectors along each axis)
3458	[color]: [colorbar]@[varn]@[units]: color of the vectors according to a 3rd variable (to be added at -v) and given [colorbar]
3459	all vectors the same length
3460	[length]: length of the wind vectors ('auto', for 9)
3461	[windlabs]= [windname],[windunits]
3462	[windname]: name of the wind variable in the graph
3463	[windunits]: units of the wind variable in the graph ('None', for the value in the file)
3464	[mapvalues]= map characteristics: [proj],[res]
3465	see full documentation: http://matplotlib.org/basemap/
3466	[proj]: projection
3467	* 'cyl', cilindric
3468	* 'lcc', lambert conformal
3469	[res]: resolution:
3470	* 'c', crude
3471	* 'l', low
3472	* 'i', intermediate
3473	* 'h', high
3474	* 'f', full
3475	gtit= title of the graph ('\|', for spaces)
3476	kindfig= kind of figure
3477	figuren= name of the figure
3478	ncfile= file to use
3479	varns= [lon],[lat],[outflow],[basinID] ',' list of the name of the variables with the lon,lat, the outflow and the basin ID
3480	"""
3481	fname = 'draw_basins'
3482
3483	if values == 'h':
3484	print fname + '_____________________________________________________________'
3485	print draw_vectors.__doc__
3486	quit()
3487
3488	expectargs = '[X/Y/Z/T]\|[dimname]\|[vardimname]\|[value]:[vecvals]:[windlabs]:' + \
3489	'[mapvalues]:[gtit]:[kindfig]:[figuren]'
3490
3491	drw.check_arguments(fname,values,expectargs,':')
3492
3493	dimvals = values.split(':')[0]
3494	vecvals = values.split(':')[1]
3495	windlabels = values.split(':')[2]
3496	mapvalues = values.split(':')[3]
3497	gtit = values.split(':')[4]
3498	kindfig = values.split(':')[5]
3499	figuren = values.split(':')[6]
3500
3501	of = NetCDFFile(ncfile,'r')
3502
3503	dims = {}
3504	for dimv in dimvals.split(','):
3505	dns = dimv.split('\|')
3506	dims[dns[0]] = [dns[1], dns[2], dns[3]]
3507
3508	varNs = []
3509	for dn in dims.keys():
3510	if dn == 'X':
3511	varNs.append(dims[dn][1])
3512	dimx = len(of.dimensions[dims[dn][0]])
3513	elif dn == 'Y':
3514	varNs.append(dims[dn][1])
3515	dimy = len(of.dimensions[dims[dn][0]])
3516
3517	ivar = 0
3518	for wvar in varns.split(','):
3519	if not drw.searchInlist(of.variables.keys(), wvar):
3520	print errormsg
3521	print ' ' + fname + ": file does not have variable '" + wvar + "' !!"
3522	quit(-1)
3523	if ivar == 0:
3524	varNs.append(wvar)
3525	else:
3526	varNs.append(wvar)
3527
3528	ivar = 0
3529	for varN in varNs:
3530	varslice = []
3531
3532	ovarN = of.variables[varN]
3533	vard = ovarN.dimensions
3534	for vdn in vard:
3535	found = False
3536	for dd in dims.keys():
3537	if dims[dd][0] == vdn:
3538	if dims[dd][2].find('@') != -1:
3539	rvals = dims[dd][2].split('@')
3540	varslice.append(slice(int(rvals[0]), int(rvals[1])))
3541	elif dims[dd][2] == '-1':
3542	varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3543	else:
3544	varslice.append(int(dims[dd][2]))
3545
3546	found = True
3547	break
3548	if not found:
3549	varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3550
3551	if varN == dims['X'][1]:
3552	lonvals0 = np.squeeze(ovarN[tuple(varslice)])
3553	elif varN == dims['Y'][1]:
3554	latvals0 = np.squeeze(ovarN[tuple(varslice)])
3555
3556	ivar = ivar + 1
3557
3558	if len(lonvals0.shape) == 1:
3559	lonvals, latvals = np.meshgrid(lonvals0, latvals0)
3560	else:
3561	lonvals = lonvals0
3562	latvals = latvals0
3563
3564	# Vector values
3565	if vecvals.split(',')[0] == 'None':
3566	freqv = None
3567	else:
3568	freqv = vecvals.split(',')[0]
3569
3570	colorvals = vecvals.split(',')[1]
3571	if len(varn.split(',')) != 3:
3572	print errormsg
3573	print ' ' + fname + ": color of vectors should be according to '" + \
3574	coln + "' but a third varibale is not provided !!"
3575	quit(-1)
3576
3577	ocolvec = of.variables[varNs[3]]
3578	colorv = ocolvec[:]
3579	stdvn, lonvn, unitsvn = drw.var_3desc(ocolvec)
3580	colorvals = colorvals + '@' + stdvn + '@' + unitsvn
3581
3582	lengthv = vecvals.split(',')[2]
3583
3584	# Vector labels
3585	windname = windlabels.split(',')[0]
3586	windunits = windlabels.split(',')[1]
3587
3588	# Vector angles
3589	oflow = ofile.variables[varNs[2]]
3590	angle = (oflow[:] - 1)*np.pi/4
3591	xflow = np.where(oflow[:] < 9, np.float(lengthv)*np.sin(angle), 0.)
3592	yflow = np.where(oflow[:] < 9, np.float(lengthv)*np.cos(angle), 0.)
3593
3594	drw.plot_basins(lonvals, latvals, xflow, yflow, freqv, colorvals, colorv, \
3595	lengthv, windname, windunits, mapvalues, gtit, kindfig, figuren)
3596
3597	of.close()
3598
3599	return
3600
3601	def draw_river_desc(ncfile, values, riverns):
3602	""" Function to plot rivers' description from ORCHIDEE's routing scheme
3603	values= [dimname]\|[vardimname]\|[value]:[basinvals]:[upstreamvals]:[mapvalues]:
3604	[gtit]:[kindfig]:[legloc]:[figuren]
3605	'X/Y'\|[dimname]\|[vardimname]\|[value]: ',', list for each basic dimension '\|' separated of:
3606	[dimname]: name of the dimension in the file
3607	[vardimname]: name of the variable with the values for the dimension in the file
3608	[value]: which value of the given dimension (-1, all; [ibeg]@[iend], i-range beginning, end)
3609	No value takes all the range of the dimension
3610	[basinsvals]= [colorline]
3611	[basincolor]: ',' list of colors of the line to use to mark the basins contours (single value also possible)
3612	[upstreamvals]= [upstreamvarn],[colorbar]
3613	[upstreamcolor]: colorbar to use to plot the basins upstream values
3614	[mapvalues]= map characteristics: [proj],[res]
3615	see full documentation: http://matplotlib.org/basemap/
3616	[proj]: projection
3617	* 'cyl', cilindric
3618	* 'lcc', lambert conformal
3619	[res]: resolution:
3620	* 'c', crude
3621	* 'l', low
3622	* 'i', intermediate
3623	* 'h', high
3624	* 'f', full
3625	gtit= title of the graph ('\|', for spaces)
3626	kindfig= kind of figure
3627	legloc= location of the legend (0, automatic)
3628	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
3629	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
3630	9: 'upper center', 10: 'center' kfig= kind of figure
3631	figuren= name of the figure
3632	ncfile= file to use
3633	riverns= ',' list of the name of the rivers to plot
3634	"""
3635	import numpy.ma as ma
3636	fname = 'draw_river_desc'
3637
3638	if values == 'h':
3639	print fname + '_____________________________________________________________'
3640	print draw_river_desc.__doc__
3641	quit()
3642
3643	expectargs = '[X/Y/Z/T]\|[dimname]\|[vardimname]\|[value]:[basinvals]:' + \
3644	'[upstreamvals]:[mapvalues]:[gtit]:[kindfig]:[legloc]:[figuren]'
3645
3646	drw.check_arguments(fname,values,expectargs,':')
3647
3648	dimvals = values.split(':')[0]
3649	basinvals = values.split(':')[1]
3650	upstreamvals = values.split(':')[2]
3651	mapvals = values.split(':')[3]
3652	gtit = values.split(':')[4]
3653	kindfig = values.split(':')[5]
3654	legloc = int(values.split(':')[6])
3655	figuren = values.split(':')[7]
3656
3657	basincol = basinvals
3658	if basincol.find(',') != 1:
3659	basincolor = basincol.split(',')
3660	else:
3661	basincolor = [basincol]
3662
3663	upstreamcolor = upstreamvals
3664
3665	of = NetCDFFile(ncfile,'r')
3666
3667	dims = {}
3668	for dimv in dimvals.split(','):
3669	dns = dimv.split('\|')
3670	dims[dns[0]] = [dns[1], dns[2], dns[3]]
3671
3672	varNs = []
3673	for dn in dims.keys():
3674	if dn == 'X':
3675	varNs.append(dims[dn][1])
3676	dimx = len(of.dimensions[dims[dn][0]])
3677	elif dn == 'Y':
3678	varNs.append(dims[dn][1])
3679	dimy = len(of.dimensions[dims[dn][0]])
3680
3681	if riverns.find(',') != -1:
3682	riverns = riverns.split(',')
3683	else:
3684	riverns = [riverns]
3685
3686	rivers = []
3687	riversubbasins = {}
3688	riversupstream = {}
3689	riversoutflow = {}
3690	for rivern in riverns:
3691	print rivern
3692
3693	# subBasins
3694	basinvar = rivern + '_coding'
3695	if not drw.searchInlist(of.variables.keys(), basinvar):
3696	print errormsg
3697	print ' ' + fname + ": file does not have variable '" + basinvar + "' !!"
3698	quit(-1)
3699	rivers.append(rivern)
3700	obasin = of.variables[basinvar]
3701	riversubbasins[rivern] = obasin[:]
3702	if rivern == riverns[0]:
3703	finalmask = obasin[:].mask
3704	else:
3705	finalmask = finalmask * obasin[:].mask
3706
3707	# upstream
3708	upstreamvar = rivern + '_upstream'
3709	if not drw.searchInlist(of.variables.keys(), upstreamvar):
3710	print errormsg
3711	print ' ' + fname + ": file does not have variable '" + upstreamvar + "' !!"
3712	quit(-1)
3713	oupstream = of.variables[upstreamvar]
3714	riversupstream[rivern] = oupstream[:]
3715	if rivern == riverns[0]:
3716	uunits = oupstream.getncattr('units')
3717
3718	# River metadata
3719	fracvar = rivern + '_frac'
3720	if not drw.searchInlist(of.variables.keys(), fracvar):
3721	print errormsg
3722	print ' ' + fname + ": file does not have variable '" + fracvar + "' !!"
3723	quit(-1)
3724	ofrac = of.variables[fracvar]
3725	riversoutflow[rivern] = [ofrac.getncattr('Longitude_of_outflow_point'), \
3726	ofrac.getncattr('Latitude_of_outflow_point')]
3727
3728	ivar = 0
3729	for varN in varNs:
3730	varslice = []
3731
3732	ovarN = of.variables[varN]
3733	vard = ovarN.dimensions
3734	for vdn in vard:
3735	found = False
3736	for dd in dims.keys():
3737	if dims[dd][0] == vdn:
3738	if dims[dd][2].find('@') != -1:
3739	rvals = dims[dd][2].split('@')
3740	varslice.append(slice(int(rvals[0]), int(rvals[1])))
3741	elif dims[dd][2] == '-1':
3742	varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3743	else:
3744	varslice.append(int(dims[dd][2]))
3745
3746	found = True
3747	break
3748	if not found:
3749	varslice.append(slice(0,len(of.dimensions[dims[dd][0]])))
3750
3751	if varN == dims['X'][1]:
3752	lonvals0 = np.squeeze(ovarN[tuple(varslice)])
3753	elif varN == dims['Y'][1]:
3754	latvals0 = np.squeeze(ovarN[tuple(varslice)])
3755
3756	ivar = ivar + 1
3757
3758	if len(lonvals0.shape) == 1:
3759	lonvals, latvals = np.meshgrid(lonvals0, latvals0)
3760	else:
3761	lonvals = lonvals0
3762	latvals = latvals0
3763
3764	# Masking only the lon,lat with rivers
3765	malonvals = ma.masked_array(lonvals, mask=finalmask)
3766	malatvals = ma.masked_array(latvals, mask=finalmask)
3767
3768	if mapvals == 'None':
3769	mapvalues = None
3770	else:
3771	mapvalues = mapvals
3772
3773	drw.plot_river_desc(malonvals, malatvals, rivers, riversubbasins, riversupstream, riversoutflow, \
3774	basincolor, upstreamcolor, uunits, mapvalues, gtit, kindfig, legloc, figuren)
3775
3776	of.close()
3777
3778	def draw_vertical_levels(ncfile, values, varn):
3779	""" plotting vertical levels distribution
3780	draw_topo_geogrid(ncfile, values, varn)
3781	ncfile= file to use
3782	values= [zlogs]:[plogs]:[title]:[graphic_kind]:[legloc]
3783	zlogs= zlog,dzlog
3784	zlog: to use logarithmic scale on the height axis ('true/false')
3785	dzlog: to use logarithmic scale on the difference of height between levels axis ('true/false')
3786	plogs= plog,dplog
3787	plog: to use logarithmic scale on the height axis ('true/false')
3788	dplog: to use logarithmic scale on the difference of height between levels axis ('true/false')
3789	title: title of the graph ('!' for spaces)
3790	graphic_kind: kind of figure (jpg, pdf, png)
3791	legloc= location of the legend (0, automatic)
3792	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
3793	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
3794	9: 'upper center', 10: 'center' kfig= kind of figure
3795	varn= [varnheight],[varnpres]
3796	varnheight: name of the variable with the height of the vertical levels
3797	'WRFz': for WRF z-levels (computed as (PH + PHB)/g, from a PHB(0,i,j) = 0)
3798	varnpres: name of the variable with the pressure of the vertical levels ('None', for no pressure plot)
3799	'WRFp': for WRF p-levels (computed as P + PB, from a PHB(0,i,j) = 0)
3800	"""
3801	fname = 'draw_vertical_levels'
3802
3803	if values == 'h':
3804	print fname + '_____________________________________________________________'
3805	print draw_vertical_levels.__doc__
3806	quit()
3807
3808	expectargs = '[zlogs]:[plogs]:[title]:[graphic_kind]:[legloc]'
3809
3810	drw.check_arguments(fname,values,expectargs,':')
3811
3812	zlog = values.split(':')[0].split(',')[0]
3813	dzlog = values.split(':')[0].split(',')[1]
3814	plog = values.split(':')[1].split(',')[0]
3815	dplog = values.split(':')[1].split(',')[1]
3816	title = values.split(':')[2].replace('!',' ')
3817	kindfig = values.split(':')[3]
3818	legloc = int(values.split(':')[4])
3819
3820	if varn.find(',') == -1:
3821	varnheight = varn
3822	varnpres = None
3823	pvals = None
3824	print warnmsg
3825	print ' ' + fname + ': assuming no pressure variable!!'
3826	else:
3827	varnheight = varn.split(',')[0]
3828	varnpres = varn.split(',')[1]
3829	if varnpres == 'None':
3830	varnpres = None
3831	pvals = None
3832
3833	if not os.path.isfile(ncfile):
3834	print errormsg
3835	print ' ' + fname + ': file "' + ncfile + '" does not exist !!'
3836	quit(-1)
3837
3838	objf = NetCDFFile(ncfile, 'r')
3839
3840	if varnheight == 'WRFz':
3841	if not gen.searchInlist(objf.variables,'PH'):
3842	print errormsg
3843	print ' ' + fname + ": WRF file '" + ncfile + "' does not have " + \
3844	"variable 'PH' !!"
3845	quit(-1)
3846	if not gen.searchInlist(objf.variables,'PHB'):
3847	print errormsg
3848	print ' ' + fname + ": WRF file '" + ncfile + "' does not have " + \
3849	"variable 'PHB' !!"
3850	quit(-1)
3851
3852	objph = objf.variables['PH']
3853	objphb = objf.variables['PHB']
3854	geop = objph[:] + objphb[:]
3855
3856	ijz0 = gen.index_mat(geop[0,], 0.)
3857	zvals = geop[0, :, ijz0[0], ijz0[1]] / 9.8
3858	else:
3859	if not gen.searchInlist(objf.variables, varnheight):
3860	print errormsg
3861	print ' ' + fname + ": file '" + ncfile + "' does not have height " + \
3862	" variable '" + varnheight + "' !!"
3863	quit(-1)
3864	objvar = objf.variables[varn]
3865	if len(objvar.shape) == 4:
3866	print warnmsg
3867	print ' ' + fname + ": assuming that height variable '" + varnheight + \
3868	"' with shape: dt,dz,dy,dx. Tacking first time-step"
3869
3870	ijz0 = gen.index_mat(objvar[0,0,], 0.)
3871	zvals = objvar[0, :, ijz0[0], ijz0[1]]
3872	elif len(objvar.shape) == 3:
3873	print warnmsg
3874	print ' ' + fname + ": assuming that height variable '" + varnheight + \
3875	"' with shape: dz,dy,dx"
3876
3877	ijz0 = gen.index_mat(objvar[0,], 0.)
3878	zvals = objvar[:, ijz0[0], ijz0[1]]
3879
3880	elif len(objvar.shape) == 2:
3881	print warnmsg
3882	print ' ' + fname + ": assuming that height variable '" + varnheight + \
3883	"' with shape: dz,dyx"
3884
3885	ijz0 = gen.index_mat(objvar[0,], 0.)
3886	zvals = objvar[:, ijz0[0]]
3887	else:
3888	zvals = objvar[:]
3889
3890	# Pressure
3891	if varnpres is not None:
3892	if varnpres == 'WRFp':
3893	if not gen.searchInlist(objf.variables,'P'):
3894	print errormsg
3895	print ' ' + fname + ": WRF file '" + ncfile + "' does not have " + \
3896	"variable 'P' !!"
3897	quit(-1)
3898	if not gen.searchInlist(objf.variables,'PB'):
3899	print errormsg
3900	print ' ' + fname + ": WRF file '" + ncfile + "' does not have " + \
3901	"variable 'PB' !!"
3902	quit(-1)
3903
3904	objph = objf.variables['P']
3905	objphb = objf.variables['PB']
3906	pres = objph[:] + objphb[:]
3907
3908	pvals = pres[0, :, ijz0[0], ijz0[1]]
3909	else:
3910	if not gen.searchInlist(objf.variables, varnpres):
3911	print errormsg
3912	print ' ' + fname + ": file '" + ncfile + "' does not have pressure " + \
3913	" variable '" + varnpres + "' !!"
3914	quit(-1)
3915	objvar = objf.variables[varnpres]
3916	if len(objvar.shape) == 4:
3917	print warnmsg
3918	print ' ' + fname + ": assuming that pressure variable '" + varnpres + \
3919	"' with shape: dt,dz,dy,dx. Tacking first time-step"
3920
3921	pvals = objvar[0, :, ijz0[0], ijz0[1]]
3922	elif len(objvar.shape) == 3:
3923	print warnmsg
3924	print ' ' + fname + ": assuming that pressure variable '" + varnpres + \
3925	"' with shape: dz,dy,dx"
3926
3927	pvals = objvar[:, ijz0[0], ijz0[1]]
3928
3929	elif len(objvar.shape) == 2:
3930	print warnmsg
3931	print ' ' + fname + ": assuming that pressure variable '" + varnpres + \
3932	"' with shape: dz,dyx"
3933
3934	pvals = objvar[:, ijz0[0]]
3935	else:
3936	pvals = objvar[:]
3937
3938	# Logarithmic axes
3939	if zlog == 'true':
3940	zlogv = True
3941	elif zlog == 'false':
3942	zlogv = False
3943	else:
3944	print errormsg
3945	print ' ' + fname + ": wrong value for zlog: '" + zlog + "' !!"
3946	print " must be either: 'true' or 'false'"
3947	quit(-1)
3948
3949	if dzlog == 'true':
3950	dzlogv = True
3951	elif dzlog == 'false':
3952	dzlogv = False
3953	else:
3954	print errormsg
3955	print ' ' + fname + ": wrong value for dzlog: '" + dzlog + "' !!"
3956	print " must be either: 'true' or 'false'"
3957	quit(-1)
3958
3959	if pvals is not None:
3960	if plog == 'true':
3961	plogv = True
3962	elif plog == 'false':
3963	plogv = False
3964	else:
3965	print errormsg
3966	print ' ' + fname + ": wrong value for plog: '" + plog + "' !!"
3967	print " must be either: 'true' or 'false'"
3968	quit(-1)
3969	if dplog == 'true':
3970	dplogv = True
3971	elif dplog == 'false':
3972	dplogv = False
3973	else:
3974	print errormsg
3975	print ' ' + fname + ": wrong value for dplog: '" + dplog + "' !!"
3976	print " must be either: 'true' or 'false'"
3977	quit(-1)
3978
3979	drw.plot_vertical_lev(zvals, pvals, zlogv, dzlogv, plogv, dplogv, title, kindfig, legloc)
3980
3981	objf.close()
3982
3983	return
3984
3985	def draw_subbasin(ncfile, values):
3986	""" Function to plot subbasin from 'routnig.nc' ORCDHIEE
3987	ncfile= file to use produced with nc_var.py#subbasin function
3988	values= [subasiname]:[rangecolors]:[mapv]:[basinlinewidth]:[drawsubid]:[gtit]:[figkind]:[legloc]:[figurename]
3989	[subasiname]= name of the subbasin ('!' for spaces)
3990	[rcolor]= '@', list of 'r\|g\|b' 1-based colors (as much as first level sub-flow). 'None' for automatic
3991	[mapv]= map characteristics: [proj],[res]
3992	see full documentation: http://matplotlib.org/basemap/
3993	[proj]: projection
3994	* 'cyl', cilindric
3995	* 'lcc', lambert conformal
3996	[res]: resolution:
3997	* 'c', crude
3998	* 'l', low
3999	* 'i', intermediate
4000	* 'h', high
4001	* 'f', full
4002	[basinlinewidth]= with of the line to draw the basin
4003	[drawsubid]= wehther sub-flow ids should be plot or not
4004	[graphtit]= title of the graph ('\|', for spaces)
4005	[lloc]= location of the legend (0, automatic)
4006	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
4007	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
4008	9: 'upper center', 10: 'center' kfig= kind of figure
4009	[figname]= name of the figure
4010	"""
4011	fname = 'draw_subbasin'
4012
4013	if values == 'h':
4014	print fname + '_____________________________________________________________'
4015	print draw_subbasin.__doc__
4016	quit()
4017
4018	expectargs = '[subasiname]:[rangecolors]:[mapv]:[basinlinewidth]:[drawsubid]:' + \
4019	'[gtit]:[figkind]:[legloc]:[figurename]'
4020
4021	drw.check_arguments(fname,values,expectargs,':')
4022
4023	subbasiname = values.split(':')[0].replace('!',' ')
4024	rangecolors = values.split(':')[1]
4025	mapv = values.split(':')[2]
4026	basinlinewidth = np.float(values.split(':')[3])
4027	drawsubid = gen.Str_Bool(values.split(':')[4])
4028	gtit = values.split(':')[5].replace('!',' ')
4029	figkind = values.split(':')[6]
4030	legloc = int(values.split(':')[7])
4031	figurename = values.split(':')[8]
4032
4033	if not os.path.isfile(ncfile):
4034	print errormsg
4035	print ' ' + fname + ': file "' + ncfile + '" does not exist !!'
4036	quit(-1)
4037
4038	objf = NetCDFFile(ncfile, 'r')
4039
4040	searchvars = ['lon', 'lat', 'lonsubflow', 'latsubflow', 'outsubflow']
4041	for searchvar in searchvars:
4042	if not gen.searchInlist(objf.variables,searchvar):
4043	print errormsg
4044	print ' ' + fname + ": WRF file '" + ncfile + "' does not have " + \
4045	"variable '" + searchvar + "' !!"
4046	quit(-1)
4047
4048	# lon,lat
4049	olon = objf.variables['lon']
4050	olat = objf.variables['lat']
4051	lon = olon[:]
4052	lat = olat[:]
4053
4054	# sub-flow names
4055	osubnames = objf.variables['subflow']
4056	subnames = drw.get_str_nc(osubnames, osubnames.shape[1])
4057
4058	# sub-flow lat, lon
4059	latlonsub = {}
4060	outflowsub = {}
4061	osublon = objf.variables['lonsubflow']
4062	osublat = objf.variables['latsubflow']
4063	oNsubflow = objf.variables['Nsubflow']
4064	ooutsubflow = objf.variables['outsubflow']
4065	Nsubflow = oNsubflow[:]
4066	isub = 0
4067	for Ssub in subnames:
4068	sublatlon = []
4069	suboutflow = []
4070	for igrid in range(Nsubflow[isub]):
4071	sublatlon.append([osublat[isub,igrid], osublon[isub,igrid]])
4072	suboutflow.append(ooutsubflow[isub,igrid])
4073	latlonsub[Ssub] = sublatlon
4074	outflowsub[Ssub] = suboutflow
4075	isub = isub + 1
4076
4077	# colors
4078	if rangecolors == 'None':
4079	rangecols = None
4080	else:
4081	cols = rangecolors.split('@')
4082	Ncols = len(cols)
4083	rangecols = []
4084	for icol in range(Ncols):
4085	cval = cols[icol].split('\|')
4086	rangecols.append([np.float(cval[0]),np.float(cval[1]),np.float(cval[2])])
4087
4088	drw.plot_subbasin(subbasiname, lon, lat, subnames, latlonsub, outflowsub, \
4089	rangecols, mapv, basinlinewidth, drawsubid, gtit, figkind, legloc, figurename)
4090
4091	objf.close()
4092
4093	return
4094
4095	#quit()
4096
4097	####### ###### ##### #### ### ## #
4098
4099	ngraphics = "'" + drw.numVector_String(namegraphics, "', '") + "'"
4100
4101	### Options
4102	##string_operation="operation to make: " + '\n' + " out, output values -S inidim1,[inidim2,...]:enddim1,[enddim2,...]"
4103	string_operation="""operation to make:
4104	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]
4105	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]
4106	[ckind]:
4107	'cmap': as it gets from colorbar
4108	'fixc,[colname]': fixed color [colname], all stright lines
4109	'fixsignc,[colname]': fixed color [colname], >0 stright, <0 dashed line
4110	"""
4111
4112	#print string_operation
4113
4114	parser = OptionParser()
4115	parser.add_option("-f", "--netCDF_file", dest="ncfile",
4116	help="file to use", metavar="FILE")
4117	parser.add_option("-o", "--operation", type='choice', dest="operation",
4118	choices=namegraphics,
4119	help="operation to make: " + ngraphics, metavar="OPER")
4120	parser.add_option("-S", "--valueS", dest="values",
4121	help="[WHEN APPLICABLE] values to use according to the operation", metavar="VALUES")
4122	parser.add_option("-v", "--variable", dest="varname",
4123	help="[WHEN APPLICABLE] variable to check", metavar="VAR")
4124
4125	(opts, args) = parser.parse_args()
4126
4127	####### #######
4128	## MAIN
4129	#######
4130
4131	# Not checking file operation
4132	Notcheckingfile = ['draw_2D_shad_cont', 'draw_2D_shad_cont_time', \
4133	'draw_2D_shad_line', 'draw_2D_shad_line_time', 'draw_lines', 'draw_lines_time', \
4134	'draw_points', 'draw_topo_geogrid_boxes', 'draw_trajectories', \
4135	'draw_vals_trajectories', 'variable_values']
4136
4137	####### ###### ##### #### ### ## #
4138	errormsg='ERROR -- error -- ERROR -- error'
4139
4140	varn=opts.varname
4141	oper=opts.operation
4142
4143	if opts.ncfile is not None and not os.path.isfile(opts.ncfile) and \
4144	not drw.searchInlist(Notcheckingfile, oper):
4145	print errormsg
4146	print ' ' + main + ': File ' + opts.ncfile + ' does not exist !!'
4147	quit(-1)
4148
4149	if oper == 'create_movie':
4150	create_movie(opts.ncfile, opts.values, opts.varname)
4151	elif oper == 'draw_2D_shad':
4152	draw_2D_shad(opts.ncfile, opts.values, opts.varname)
4153	elif oper == 'draw_2D_shad_time':
4154	draw_2D_shad_time(opts.ncfile, opts.values, opts.varname)
4155	elif oper == 'draw_2D_shad_cont':
4156	draw_2D_shad_cont(opts.ncfile, opts.values, opts.varname)
4157	elif oper == 'draw_2D_shad_cont_time':
4158	draw_2D_shad_cont_time(opts.ncfile, opts.values, opts.varname)
4159	elif oper == 'draw_2D_shad_line':
4160	draw_2D_shad_line(opts.ncfile, opts.values, opts.varname)
4161	elif oper == 'draw_2D_shad_line_time':
4162	draw_2D_shad_line_time(opts.ncfile, opts.values, opts.varname)
4163	elif oper == 'draw_barbs':
4164	draw_barbs(opts.ncfile, opts.values, opts.varname)
4165	elif oper == 'draw_basins':
4166	draw_basins(opts.ncfile, opts.values, opts.varname)
4167	elif oper == 'draw_Neighbourghood_evol':
4168	draw_Neighbourghood_evol(opts.ncfile, opts.values, opts.varname)
4169	elif oper == 'draw_lines':
4170	draw_lines(opts.ncfile, opts.values, opts.varname)
4171	elif oper == 'draw_lines_time':
4172	draw_lines_time(opts.ncfile, opts.values, opts.varname)
4173	elif oper == 'draw_points':
4174	draw_points(opts.ncfile, opts.values)
4175	elif oper == 'draw_points_lonlat':
4176	draw_points_lonlat(opts.ncfile, opts.values)
4177	elif oper == 'draw_ptZvals':
4178	draw_ptZvals(opts.ncfile, opts.values, opts.varname)
4179	elif oper == 'draw_river_desc':
4180	draw_river_desc(opts.ncfile, opts.values, opts.varname)
4181	elif oper == 'draw_subbasin':
4182	draw_subbasin(opts.ncfile, opts.values)
4183	elif oper == 'draw_timeSeries':
4184	draw_timeSeries(opts.ncfile, opts.values, opts.varname)
4185	elif oper == 'draw_topo_geogrid':
4186	draw_topo_geogrid(opts.ncfile, opts.values)
4187	elif oper == 'draw_topo_geogrid_boxes':
4188	draw_topo_geogrid_boxes(opts.ncfile, opts.values)
4189	elif oper == 'draw_trajectories':
4190	draw_trajectories(opts.ncfile, opts.values, opts.varname)
4191	elif oper == 'draw_vals_trajectories':
4192	draw_vals_trajectories(opts.ncfile, opts.values, opts.varname)
4193	elif oper == 'draw_vectors':
4194	draw_vectors(opts.ncfile, opts.values, opts.varname)
4195	elif oper == 'draw_vertical_levels':
4196	draw_vertical_levels(opts.ncfile, opts.values, opts.varname)
4197	elif oper == 'list_graphics':
4198	# From: http://www.diveintopython.net/power_of_introspection/all_together.html
4199	import drawing as myself
4200	object = myself
4201	for opern in namegraphics:
4202	if opern != 'list_graphics':
4203	print opern + '_______ ______ _____ ____ ___ __ _'
4204	print getattr(object, opern).__doc__
4205	elif oper == 'variable_values':
4206	variable_values(opts.values)
4207	else:
4208	print errormsg
4209	print ' ' + main + ": the graphic '" + oper + "' is not ready !!"
4210	print errormsg
4211	quit()

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

Download in other formats: