Context Navigation

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

Last change on this file since 947 was 947, checked in by lfita, 9 years ago
Fixing `months' time units in` draw_lines_time'
File size: 157.7 KB

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

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

Download in other formats: