Context Navigation

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

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

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

Download in other formats: