Context Navigation

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

Last change on this file since 489 was 489, checked in by lfita, 9 years ago
Adding use of `dxdy_lonlatDIMS' on 'shad_contour' Removing some residual debug printings...
File size: 96.0 KB

Line
1	import numpy as np
2	import os
3	from netCDF4 import Dataset as NetCDFFile
4	import drawing_tools as drw
5	from optparse import OptionParser
6	import sys
7	from cStringIO import StringIO
8
9	## 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
10	## 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
11	## 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
12	## 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
13	## 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'
14	## e.g. # drawing.py -f 'geo_em.d02.nc' -o draw_topo_geogrid -S '0.,3000.:None:FF_3dom d02:png:cyl,i'
15	## e.g. # drawing.py -o draw_topo_geogrid_boxes -S '0.,3000.:None:FF domains:pdf:lcc,i:d01,d02,d03' -f geo_em.d01.nc,geo_em.d02.nc,geo_em.d03.nc
16	## 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
17	## 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
18	## 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'
19	## e.g. # drawing.py -o variable_values -S PSFC
20	## 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'
21	## 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
22	## 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
23
24	main = 'drawing.py'
25
26	errormsg = 'ERROR -- error -- ERROR -- error'
27	warnmsg = 'WARNING -- waring -- WARNING -- warning'
28	fillValue=1.e20
29
30	namegraphics = ['create_movie', 'draw_2D_shad', 'draw_2D_shad_time', \
31	'draw_2D_shad_cont', 'draw_2D_shad_cont_time', 'draw_2D_shad_line', \
32	'draw_2D_shad_line_time', 'draw_timeSeries', 'draw_topo_geogrid', \
33	'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories', \
34	'draw_lines', 'draw_lines_time', 'draw_Neighbourghood_evol', 'list_graphics', \
35	'variable_values']
36
37	def draw_2D_shad(ncfile, values, varn):
38	""" plotting a fields with shading
39	draw_2D_shad(ncfile, values, varn)
40	ncfile= file to use
41	values=[vnamefs]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
42	[kindfig]:[reverse]:[mapv]:[close]
43	[vnamefs]: Name in the figure of the variable to be shaded
44	[dimvals]: ',' list of [dimname]\|[value] telling at which dimension of the
45	variable a given value is required (-1, all the length)
46	[dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
47	[colorbar]: name of the color bar
48	[smin/axv]: minimum and maximum value for the shading or:
49	'Srange': for full range
50	'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
51	'Saroundminmax@val': for minval,maxval
52	'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
53	percentile_(100-val)-median)
54	'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min@val,max@val-mean)
55	'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
56	'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
57	percentile_(100-val)-median)
58	[figt]: title of the figure ('\|' for spaces)
59	[kindfig]: kind of figure
60	[reverse]: Transformation of the values
61	* 'transpose': reverse the axes (x-->y, y-->x)
62	* 'flip'@[x/y]: flip the axis x or y
63	[mapv]: map characteristics: [proj],[res]
64	see full documentation: http://matplotlib.org/basemap/
65	[proj]: projection
66	* 'cyl', cilindric
67	* 'lcc', lamvbert conformal
68	[res]: resolution:
69	* 'c', crude
70	* 'l', low
71	* 'i', intermediate
72	* 'h', high
73	* 'f', full
74	valules= 'rh:z\|-1,x\|-1:z\|-1,x\|-1:lat:pressure:BuPu:0.,100.:rh:pdf:flip@y:None'
75	varn= [varsn] name of the variable to plot with shading
76	"""
77
78	fname = 'draw_2D_shad'
79	if values == 'h':
80	print fname + '_____________________________________________________________'
81	print draw_2D_shad.__doc__
82	quit()
83
84	expectargs = ['[vnamefs]', '[dimvals]', '[dimxvn]', '[dimyvn]', '[colorbar]', \
85	'[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', '[mapv]', '[close]']
86
87	drw.check_arguments(fname,len(expectargs),values,':',expectargs)
88
89	vnamesfig = values.split(':')[0]
90	dimvals= values.split(':')[1].replace('\|',':')
91	vdimxn = values.split(':')[2]
92	vdimyn = values.split(':')[3]
93	colbarn = values.split(':')[4]
94	shadminmax = values.split(':')[5]
95	figtitle = values.split(':')[6].replace('\|',' ')
96	figkind = values.split(':')[7]
97	revals = values.split(':')[8]
98	mapvalue = values.split(':')[9]
99	# varn = values.split(':')[10]
100
101	ncfiles = ncfile
102
103	if not os.path.isfile(ncfiles):
104	print errormsg
105	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
106	quit(-1)
107
108	objsf = NetCDFFile(ncfiles, 'r')
109
110	varns = varn.split(',')[0]
111
112	if not objsf.variables.has_key(varns):
113	print errormsg
114	print ' ' + fname + ': shading file "' + ncfiles + \
115	'" does not have variable "' + varns + '" !!'
116	quit(-1)
117
118	# Variables' values
119	objvars = objsf.variables[varns]
120
121	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
122
123	# Dimensions names
124	## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
125	## dimnamesv = []
126	## for idd in range(len(objvars.dimensions)):
127	## cutdim = False
128	## for idc in range(len(dimvals.split(','))):
129	## dimcutn = dimvals.split(',')[idc].split(':')[0]
130	## print objvars.dimensions[idd], dimcutn
131	## if objvars.dimensions[idd] == dimcutn:
132	## cutdim = True
133	## break
134	## if not cutdim: dimnamesv.append(objvars.dimensions[idd])
135	dimnamesv = [vdimyn, vdimxn]
136
137	if drw.searchInlist(objvars.ncattrs(),'units'):
138	varunits = objvars.getncattr('units')
139	else:
140	print warnmsg
141	print ' ' + fname + ": variable '" + varn + "' without units!!"
142	varunits = '-'
143
144	if not objsf.variables.has_key(vdimxn):
145	print errormsg
146	print ' ' + fname + ': shading file "' + ncfiles + \
147	'" does not have dimension variable "' + vdimxn + '" !!'
148	quit(-1)
149	if not objsf.variables.has_key(vdimyn):
150	print errormsg
151	print ' ' + fname + ': shading file "' + ncfiles + \
152	'" does not have dimension variable "' + vdimyn + '" !!'
153	quit(-1)
154
155	objdimx = objsf.variables[vdimxn]
156	objdimy = objsf.variables[vdimyn]
157	if drw.searchInlist(objdimx.ncattrs(),'units'):
158	odimxu = objdimx.getncattr('units')
159	else:
160	print warnmsg
161	print ' ' + fname + ": variable dimension '" + vdimxn + "' without units!!"
162	odimxu = '-'
163
164	if drw.searchInlist(objdimy.ncattrs(),'units'):
165	odimyu = objdimy.getncattr('units')
166	else:
167	print warnmsg
168	print ' ' + fname + ": variable dimension '" + vdimyn + "' without units!!"
169	odimyu = '-'
170
171	odimxv, odimyv = drw.dxdy_lonlatDIMS(objdimx[:], objdimy[:], objdimx.dimensions, \
172	objdimy.dimensions, dimvals.replace(':','\|').split(','))
173
174
175	# if len(objdimx.shape) <= 2:
176	## odimxv = objdimx[valshad.shape]
177	## odimyv = objdimy[valshad.shape]
178	# odimxv = objdimx[:]
179	# odimyv = objdimy[:]
180
181	# elif len(objdimx.shape) == 3:
182	## dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])]
183	## odimxv = objdimx[tuple(dimcut)]
184	## odimyv = objdimy[tuple(dimcut)]
185	# odimxv = objdimx[0,:]
186	# odimyv = objdimy[0,:]
187	# else:
188	# print errormsg
189	# print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
190	# ' not ready!!'
191	# quit(-1)
192
193	shading_nx = []
194	if shadminmax.split(',')[0][0:1] != 'S':
195	shading_nx.append(np.float(shadminmax.split(',')[0]))
196	else:
197	shading_nx.append(shadminmax.split(',')[0])
198
199	if shadminmax.split(',')[1][0:1] != 'S':
200	shading_nx.append(np.float(shadminmax.split(',')[1]))
201	else:
202	shading_nx.append(shadminmax.split(',')[1])
203
204	if mapvalue == 'None': mapvalue = None
205
206	drw.plot_2D_shadow(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu, dimnamesv,\
207	colbarn, shading_nx, varunits, figtitle, figkind, revals, mapvalue, True)
208
209	return
210
211	def draw_2D_shad_time(ncfile, values, varn):
212	""" plotting a fields with shading with time values
213	draw_2D_shad(ncfile, values, varn)
214	ncfile= file to use
215	values=[vnamefs]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
216	[kindfig]:[reverse]:[timevals]:[close]
217	[vnamefs]: Name in the figure of the variable to be shaded
218	[dimvals]: ',' list of [dimname]\|[value] telling at which dimension of the
219	variable a given value is required (-1, all the length)
220	[dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
221	[colorbar]: name of the color bar
222	[smin/axv]: minimum and maximum value for the shading or:
223	'Srange': for full range
224	'Saroundmean@val': for mean-xtrm,mean+xtrm where xtrm = np.min(mean-min@val,max@val-mean)
225	'Saroundminmax@val': for minval,maxval
226	'Saroundpercentile@val': for median-xtrm,median+xtrm where xtrm = np.min(median-percentile_(val),
227	percentile_(100-val)-median)
228	'Smean@val': for -xtrm,xtrm where xtrm = np.min(mean-min@val,max@val-mean)
229	'Smedian@val': for -xtrm,xtrm where xtrm = np.min(median-min@val,max@val-median)
230	'Spercentile@val': for -xtrm,xtrm where xtrm = np.min(median-percentile_(val),
231	percentile_(100-val)-median)
232	[figt]: title of the figure ('\|' for spaces)
233	[kindfig]: kind of figure
234	[reverse]: Transformation of the values
235	* 'transpose': reverse the axes (x-->y, y-->x)
236	* 'flip'@[x/y]: flip the axis x or y
237	[timevals]: [timen]\|[units]\|[kind]\|[tfmt]\|[label]\|[timeaxis] time labels characteristics
238	[timen]; name of the time variable
239	[units]; units string according to CF conventions ([tunits] since
240	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
241	[kind]; kind of output
242	'Nval': according to a given number of values as 'Nval',[Nval]
243	'exct': according to an exact time unit as 'exct',[tunit];
244	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
245	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
246	'l': milisecond
247	[tfmt]; desired format
248	[label]; label at the graph ('!' for spaces)
249	[close]: should figure be closed (finished)
250	values='dtcon:Time\|-1,bottom_top\|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly\|'
251	'dtcon:pdf:transpose:time\|hours!since!1949-12-01\|exct,2,d\|%d\|date!([DD])\|x:True
252	varn= [varsn] name of the variable to plot with shading
253	"""
254	fname = 'draw_2D_shad_time'
255
256	if values == 'h':
257	print fname + '_____________________________________________________________'
258	print draw_2D_shad_time.__doc__
259	quit()
260
261	farguments = ['[vnamefs]', '[dimvals]', '[dimxvn]', '[dimyvn]', '[colorbar]', \
262	'[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', '[timevals]', '[close]']
263	drw.check_arguments(fname,11,values,':',farguments)
264
265	vnamesfig = values.split(':')[0]
266	dimvals= values.split(':')[1].replace('\|',':')
267	vdimxn = values.split(':')[2]
268	vdimyn = values.split(':')[3]
269	colbarn = values.split(':')[4]
270	shadminmax = values.split(':')[5]
271	figtitle = values.split(':')[6].replace('\|',' ')
272	figkind = values.split(':')[7]
273	revals = values.split(':')[8]
274	timevals = values.split(':')[9]
275	close = values.split(':')[10]
276
277	ncfiles = ncfile
278
279	if not os.path.isfile(ncfiles):
280	print errormsg
281	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
282	quit(-1)
283
284	objsf = NetCDFFile(ncfiles, 'r')
285
286	varns = varn.split(',')[0]
287
288	if not objsf.variables.has_key(varns):
289	print errormsg
290	print ' ' + fname + ': shading file "' + ncfiles + \
291	'" does not have variable "' + varns + '" !!'
292	quit(-1)
293
294	# Variables' values
295	objvars = objsf.variables[varns]
296
297	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
298
299	dimnamesv = [vdimyn, vdimxn]
300
301	varunits = objvars.getncattr('units')
302
303	if not objsf.variables.has_key(vdimxn):
304	print errormsg
305	print ' ' + fname + ': shading file "' + ncfiles + \
306	'" does not have dimension variable "' + vdimxn + '" !!'
307	quit(-1)
308	if not objsf.variables.has_key(vdimyn):
309	print errormsg
310	print ' ' + fname + ': shading file "' + ncfiles + \
311	'" does not have dimensino variable "' + vdimyn + '" !!'
312	quit(-1)
313
314	objdimx = objsf.variables[vdimxn]
315	objdimy = objsf.variables[vdimyn]
316	odimxu = objdimx.getncattr('units')
317	odimyu = objdimy.getncattr('units')
318
319	if len(objdimx.shape) <= 2:
320	odimxv = objdimx[:]
321	odimyv = objdimy[:]
322
323	elif len(objdimx.shape) == 3:
324	odimxv = objdimx[0,:]
325	odimyv = objdimy[0,:]
326	else:
327	print errormsg
328	print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
329	' not ready!!'
330	quit(-1)
331
332	timename = timevals.split('\|')[0]
333	timeunit = timevals.split('\|')[1].replace('!',' ')
334	timekind = timevals.split('\|')[2]
335	timefmt = timevals.split('\|')[3]
336	timelabel = timevals.split('\|')[4].replace('!',' ')
337
338	if vdimxn == timename:
339	odimxv = objsf.variables[vdimxn][:]
340	odimxu = timelabel
341	timeaxis = 'x'
342	odimyv = objsf.variables[vdimyn]
343	odimyu = odimyv.getncattr('units')
344	timepos, timelabels = drw.CFtimes_plot(odimxv, timeunit, timekind, timefmt)
345	elif vdimyn == timename:
346	odimyv = objsf.variables[vdimyn][:]
347	odimyu = timelabel
348	timeaxis = 'y'
349	odimxv = objsf.variables[vdimxn]
350	odimxu = odimxv.getncattr('units')
351	timepos, timelabels = drw.CFtimes_plot(odimyv, timeunit, timekind, timefmt)
352	else:
353	print errormsg
354	print ' ' + fname + ": time variable '" + timename + "' not found!!"
355	quit(-1)
356
357	shading_nx = []
358	if shadminmax.split(',')[0][0:1] != 'S':
359	shading_nx.append(np.float(shadminmax.split(',')[0]))
360	else:
361	shading_nx.append(shadminmax.split(',')[0])
362
363	if shadminmax.split(',')[1][0:1] != 'S':
364	shading_nx.append(np.float(shadminmax.split(',')[1]))
365	else:
366	shading_nx.append(shadminmax.split(',')[1])
367
368	closeval = drw.Str_Bool(close)
369
370	drw.plot_2D_shadow_time(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu, \
371	dimnamesv, colbarn, shading_nx, varunits, figtitle, figkind, revals, timeaxis, \
372	timepos, timelabels, closeval)
373
374	return
375
376	def draw_2D_shad_cont(ncfile, values, varn):
377	""" plotting two fields, one with shading and the other with contour lines
378	draw_2D_shad_cont(ncfile, values, varn)
379	ncfile= [ncfilevars],[ncfilevarc] files to use (one value, same file)
380	values=[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:[figt]:[kindfig]:[reverse]:[mapv]
381	[vnamefs],[vnamefc]: Name in the figure of the shaded and the contour variables
382	[dimvals/c]: list of [dimname]\|[value] telling at which dimension of the
383	variable a given value is required (no dimension name, all the length)
384	[dimx/yvn]: names of the variables with the values of the dimensions for the plot
385	[colorbar]: name of the color bar
386	[ckind]: kind of contours
387	'cmap': as it gets from colorbar
388	'fixc,[colname]': fixed color [colname], all stright lines
389	'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed line
390	[clabfmt]: format of the labels in the contour (None, also possible)
391	[smin/axv]: minimum and maximum value for the shading
392	[sminc]:[smaxv]:[Nlev]: minimum, maximum and number of values for the contour
393	[figt]: title of the figure ('\|' for spaces)
394	[kindfig]: kind of figure
395	[reverse]: does the values be transposed? 'True/False',
396	[mapv]: map characteristics: [proj],[res]
397	see full documentation: http://matplotlib.org/basemap/
398	[proj]: projection
399	* 'cyl', cilindric
400	* 'lcc', lamvbert conformal
401	[res]: resolution:
402	* 'c', crude
403	* 'l', low
404	* 'i', intermediate
405	* 'h', high
406	* 'f', full
407	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'
408	varn= [varsn],[varcn] name of the variable to plot with shading variable with contour
409	"""
410
411	fname = 'draw_2D_shad_cont'
412	if values == 'h':
413	print fname + '_____________________________________________________________'
414	print draw_2D_shad_cont.__doc__
415	quit()
416
417	expectargs = '[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:'
418	expectargs = expectargs + '[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],'
419	expectargs = expectargs + '[Nlev]:[figt]:[kindfig]:[reverse]:[mapv]'
420
421	drw.check_arguments(fname,len(expectargs.split(':')),values,':',expectargs)
422
423	vnamesfig = values.split(':')[0].split(',')
424	dimvals= values.split(':')[1].replace('\|',':')
425	dimvalc= values.split(':')[2].replace('\|',':')
426	vdimxn = values.split(':')[3]
427	vdimyn = values.split(':')[4]
428	colbarn = values.split(':')[5]
429	countkind = values.split(':')[6]
430	countlabelfmt = values.split(':')[7]
431	shadminmax = values.split(':')[8]
432	contlevels = values.split(':')[9]
433	figtitle = values.split(':')[10].replace('\|',' ')
434	figkind = values.split(':')[11]
435	revals = values.split(':')[12]
436	mapvalue = values.split(':')[13]
437
438	if2filenames = ncfile.find(',')
439
440	if if2filenames != -1:
441	ncfiles = ncfile.split(',')[0]
442	ncfilec = ncfile.split(',')[1]
443	else:
444	ncfiles = ncfile
445	ncfilec = ncfile
446
447	if not os.path.isfile(ncfiles):
448	print errormsg
449	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
450	quit(-1)
451
452	if not os.path.isfile(ncfilec):
453	print errormsg
454	print ' ' + fname + ': contour file "' + ncfilec + '" does not exist !!'
455	quit(-1)
456
457	objsf = NetCDFFile(ncfiles, 'r')
458	objcf = NetCDFFile(ncfilec, 'r')
459
460	varns = varn.split(',')[0]
461	varnc = varn.split(',')[1]
462
463	if not objsf.variables.has_key(varns):
464	print errormsg
465	print ' ' + fname + ': shading file "' + ncfiles + \
466	'" does not have variable "' + varns + '" !!'
467	quit(-1)
468
469	if not objcf.variables.has_key(varnc):
470	print errormsg
471	print ' ' + fname + ': contour file "' + ncfilec + \
472	'" does not have variable "' + varnc + '" !!'
473	quit(-1)
474
475	# Variables' values
476	objvars = objsf.variables[varns]
477	objvarc = objcf.variables[varnc]
478
479	if len(objvars.shape) != len(objvarc.shape):
480	print errormsg
481	print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \
482	objvars.shape, 'different than contour variable "' + varnc + '": ', \
483	objvarc.shape,' !!!'
484	quit(-1)
485
486	for idim in range(len(objvars.shape)):
487	if objvars.shape[idim] != objvarc.shape[idim]:
488	print errormsg
489	print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \
490	objvars.shape, 'different than contour variable "' + varnc + '": ', \
491	objvarc.shape,' !!!'
492	quit(-1)
493
494	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
495	valcont, dimscont = drw.slice_variable(objvarc, dimvalc.replace(',','\|'))
496
497	# Dimensions names
498	## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
499	## dimnamesv = []
500	## for idd in range(len(objvars.dimensions)):
501	## cutdim = False
502	## for idc in range(len(dimvals.split(','))):
503	## dimcutn = dimvals.split(',')[idc].split(':')[0]
504	## print objvars.dimensions[idd], dimcutn
505	## if objvars.dimensions[idd] == dimcutn:
506	## cutdim = True
507	## break
508	## if not cutdim: dimnamesv.append(objvars.dimensions[idd])
509	dimnamesv = [vdimyn, vdimxn]
510
511	varunits = []
512	varunits.append(objvars.getncattr('units'))
513	varunits.append(objvarc.getncattr('units'))
514
515	if not objsf.variables.has_key(vdimxn):
516	print errormsg
517	print ' ' + fname + ': shading file "' + ncfiles + \
518	'" does not have dimension variable "' + vdimxn + '" !!'
519	quit(-1)
520	if not objsf.variables.has_key(vdimyn):
521	print errormsg
522	print ' ' + fname + ': shading file "' + ncfiles + \
523	'" does not have dimensino variable "' + vdimyn + '" !!'
524	quit(-1)
525
526	objdimx = objsf.variables[vdimxn]
527	objdimy = objsf.variables[vdimyn]
528	odimxu = objdimx.getncattr('units')
529	odimyu = objdimy.getncattr('units')
530
531	# Getting only that dimensions with coincident names
532	odimxv, odimyv = drw.dxdy_lonlatDIMS(objdimx[:], objdimy[:], objdimx.dimensions, \
533	objdimy.dimensions, dimvals.replace(':','\|').split(','))
534
535	# dimnvx = objdimx.dimensions
536	# cutslice = []
537	# for idimn in objdimx.dimensions:
538	# found = False
539	# for dimsn in dimsshad:
540	# if idimn == dimsn:
541	# cutslice.append(slice(0,len(objsf.dimensions[idimn])))
542	# found = True
543	# if not found: cutslice.append(0)
544	#
545	# odimxv = objdimx[tuple(cutslice)]
546	#
547	# dimnvy = objdimy.dimensions
548	# cutslice = []
549	# for idimn in objdimy.dimensions:
550	# found = False
551	# for dimsn in dimsshad:
552	# if idimn == dimsn:
553	# cutslice.append(slice(0,len(objsf.dimensions[idimn])))
554	# found = True
555	# if not found: cutslice.append(0)
556	#
557	# odimyv = objdimy[tuple(cutslice)]
558
559	# if len(objdimx.shape) <= 2:
560	# odimxv = objdimx[:]
561	# odimyv = objdimy[:]
562	# elif len(objdimx.shape) == 3:
563	# odimxv = objdimx[0,:]
564	# odimyv = objdimy[0,:]
565	# else:
566	# print errormsg
567	# print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
568	# ' not ready!!'
569	# quit(-1)
570
571	if countlabelfmt == 'None':
572	countlfmt = None
573	else:
574	countlfmt = countlabelfmt
575
576	shading_nx = np.zeros((2), dtype=np.float)
577	shading_nx[0] = np.float(shadminmax.split(',')[0])
578	shading_nx[1] = np.float(shadminmax.split(',')[1])
579
580	clevmin = np.float(contlevels.split(',')[0])
581	clevmax = np.float(contlevels.split(',')[1])
582	Nclevels = int(contlevels.split(',')[2])
583
584	levels_cont = drw.pretty_int(clevmin, clevmax, Nclevels)
585
586	if len(levels_cont) <= 1:
587	print warnmsg
588	print ' ' + fname + ': wrong contour levels:', levels_cont,' !!'
589	del(levels_cont)
590	levels_cont = np.zeros((Nclevels), dtype=np.float)
591	levels_cont = np.arange(7)*(clevmax - clevmin)/(Nclevels-1)
592	print ' generating default ones: ',levels_cont
593
594	if mapvalue == 'None': mapvalue = None
595
596	drw.plot_2D_shadow_contour(valshad, valcont, vnamesfig, odimxv, odimyv, odimxu, \
597	odimyu, dimnamesv, colbarn, countkind, countlfmt, shading_nx, levels_cont, \
598	varunits, figtitle, figkind, revals, mapvalue)
599
600	return
601
602	def draw_2D_shad_cont_time(ncfile, values, varn):
603	""" plotting two fields, one with shading and the other with contour lines being
604	one of the dimensions of time characteristics
605	draw_2D_shad_cont(ncfile, values, varn)
606	ncfile= [ncfilevars],[ncfilevarc] files to use (one value, same file)
607	values=[vnamefs]:[dimvals]:[dimvalc]:[dimxvn]:[dimyvn]:[colorbar]:[ckind]:[clabfmt]:[sminv],[smaxv]:[sminc],[smaxv],[Nlev]:[figt]:[kindfig]:[reverse]:[timevals]:[mapv]
608	[vnamefs],[vnamefc]: Name in the figure of the shaded and the contour variables
609	[dimvals/c]: list of [dimname]\|[value] telling at which dimension of the
610	variable a given value is required (no dimension name, all the length)
611	[dimx/yvn]: ',' list with the name of the variables with the values of the dimensions
612	[colorbar]: name of the color bar
613	[ckind]: kind of contours
614	'cmap': as it gets from colorbar
615	'fixc,[colname]': fixed color [colname], all stright lines
616	'fixsigc,[colname]': fixed color [colname], >0 stright, <0 dashed line
617	[clabfmt]: format of the labels in the contour (None, also possible)
618	[smin/axv]: minimum and maximum value for the shading
619	[sminc]:[smaxv]:[Nlev]: minimum, maximum and number of values for the contour
620	[figt]: title of the figure ('\|' for spaces)
621	[kindfig]: kind of figure
622	[reverse]: modification to the dimensions:
623	'transposed': transpose matrices
624	'flip',[x/y]: flip only the dimension [x] or [y]
625	[timevals]: [timen]\|[units]\|[kind]\|[tfmt]\|[label] time labels characteristics
626	[timen]; name of the time variable
627	[units]; units string according to CF conventions ([tunits] since
628	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
629	[kind]; kind of output
630	'Nval': according to a given number of values as 'Nval',[Nval]
631	'exct': according to an exact time unit as 'exct',[tunit];
632	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
633	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
634	'l': milisecond
635	[tfmt]; desired format
636	[label]; label at the graph ('!' for spaces)
637	[mapv]: map characteristics: [proj],[res]
638	see full documentation: http://matplotlib.org/basemap/
639	[proj]: projection
640	* 'cyl', cilindric
641	* 'lcc', lamvbert conformal
642	[res]: resolution:
643	* 'c', crude
644	* 'l', low
645	* 'i', intermediate
646	* 'h', high
647	* 'f', full
648	valules= 'rh,ta:z\|-1,x\|-1:z\|-1,x\|-1:lat:pressure:BuPu:fixsigc,black:%d:0.,100.:195.,305.,7:Meridonal\|average\|of\|rh\|&\|ta:pdf:flip@y:time!hours!since!1949/12/01\|exct,5d\|%d\|date!([DD]):None'
649	varn= [varsn],[varcn] name of the variable to plot with shading variable with contour
650	"""
651
652	fname = 'draw_2D_shad_cont_time'
653	if values == 'h':
654	print fname + '_____________________________________________________________'
655	print draw_2D_shad_cont_time.__doc__
656	quit()
657
658	expectargs = ['[vnamefs]', '[dimvals]', '[dimvalc]', '[dimxvn]', '[dimyvn]', \
659	'[colorbar]', '[ckind]', '[clabfmt]', '[sminv],[smaxv]', \
660	'[sminc],[smaxv],[Nlev]', '[figt]', '[kindfig]', '[reverse]', '[timevals]', \
661	'[mapv]']
662
663	drw.check_arguments(fname,len(expectargs),values,':',expectargs)
664
665	vnamesfig = values.split(':')[0].split(',')
666	dimvals= values.split(':')[1].replace('\|',':')
667	dimvalc= values.split(':')[2].replace('\|',':')
668	vdimxn = values.split(':')[3]
669	vdimyn = values.split(':')[4]
670	colbarn = values.split(':')[5]
671	countkind = values.split(':')[6]
672	countlabelfmt = values.split(':')[7]
673	shadminmax = values.split(':')[8]
674	contlevels = values.split(':')[9]
675	figtitle = values.split(':')[10].replace('\|',' ')
676	figkind = values.split(':')[11]
677	revals = values.split(':')[12]
678	timevals = values.split(':')[13]
679	mapvalue = values.split(':')[14]
680
681	if2filenames = ncfile.find(',')
682
683	if if2filenames != -1:
684	ncfiles = ncfile.split(',')[0]
685	ncfilec = ncfile.split(',')[1]
686	else:
687	ncfiles = ncfile
688	ncfilec = ncfile
689
690	if not os.path.isfile(ncfiles):
691	print errormsg
692	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
693	quit(-1)
694
695	if not os.path.isfile(ncfilec):
696	print errormsg
697	print ' ' + fname + ': contour file "' + ncfilec + '" does not exist !!'
698	quit(-1)
699
700	objsf = NetCDFFile(ncfiles, 'r')
701	objcf = NetCDFFile(ncfilec, 'r')
702
703	varns = varn.split(',')[0]
704	varnc = varn.split(',')[1]
705
706	if not objsf.variables.has_key(varns):
707	print errormsg
708	print ' ' + fname + ': shading file "' + ncfiles + \
709	'" does not have variable "' + varns + '" !!'
710	quit(-1)
711
712	if not objcf.variables.has_key(varnc):
713	print errormsg
714	print ' ' + fname + ': contour file "' + ncfilec + \
715	'" does not have variable "' + varnc + '" !!'
716	quit(-1)
717
718	# Variables' values
719	objvars = objsf.variables[varns]
720	objvarc = objcf.variables[varnc]
721
722	if len(objvars.shape) != len(objvarc.shape):
723	print errormsg
724	print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \
725	objvars.shape, 'different than contour variable "' + varnc + '": ', \
726	objvarc.shape,' !!!'
727	quit(-1)
728
729	for idim in range(len(objvars.shape)):
730	if objvars.shape[idim] != objvarc.shape[idim]:
731	print errormsg
732	print ' ' + fname + ': shading variable "' + varns + '" has a shape: ', \
733	objvars.shape, 'different than contour variable "' + varnc + '": ', \
734	objvarc.shape,' !!!'
735	quit(-1)
736
737	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
738	valcont, dimscont = drw.slice_variable(objvarc, dimvalc.replace(',','\|'))
739
740	# Dimensions names
741	## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
742	## dimnamesv = []
743	## for idd in range(len(objvars.dimensions)):
744	## cutdim = False
745	## for idc in range(len(dimvals.split(','))):
746	## dimcutn = dimvals.split(',')[idc].split(':')[0]
747	## print objvars.dimensions[idd], dimcutn
748	## if objvars.dimensions[idd] == dimcutn:
749	## cutdim = True
750	## break
751	## if not cutdim: dimnamesv.append(objvars.dimensions[idd])
752	dimnamesv = [vdimyn, vdimxn]
753
754	varunits = []
755	varunits.append(objvars.getncattr('units'))
756	varunits.append(objvarc.getncattr('units'))
757
758	if not objsf.variables.has_key(vdimxn):
759	print errormsg
760	print ' ' + fname + ': shading file "' + ncfiles + \
761	'" does not have dimension variable "' + vdimxn + '" !!'
762	quit(-1)
763	if not objsf.variables.has_key(vdimyn):
764	print errormsg
765	print ' ' + fname + ': shading file "' + ncfiles + \
766	'" does not have dimension variable "' + vdimyn + '" !!'
767	quit(-1)
768
769	timename = timevals.split('\|')[0]
770	timeunit = timevals.split('\|')[1].replace('!',' ')
771	timekind = timevals.split('\|')[2]
772	timefmt = timevals.split('\|')[3]
773	timelabel = timevals.split('\|')[4].replace('!',' ')
774
775	if vdimxn == timename:
776	timevals = objsf.variables[vdimxn][:]
777	timedims = objsf.variables[vdimxn].dimensions
778	dimt = 'x'
779	ovalaxis = objsf.variables[vdimyn]
780	ovalu = ovalaxis.getncattr('units')
781	elif vdimyn == timename:
782	timevals = objsf.variables[vdimyn][:]
783	timedims = objsf.variables[vdimyn].dimensions
784	dimt = 'y'
785	ovalaxis = objsf.variables[vdimxn]
786	ovalu = ovalaxis.getncattr('units')
787	else:
788	print errormsg
789	print ' ' + fname + ": time variable '" + timename + "' not found!!"
790	quit(-1)
791
792	timepos, timelabels = drw.CFtimes_plot(timevals, timeunit, timekind, timefmt)
793
794	# Getting only that dimensions with coincident names
795	dimnvx = ovalaxis.dimensions
796
797	cutslice = []
798	for idimn in dimsshad:
799	found = False
800	for dimsn in dimnvx:
801	if idimn == dimsn:
802	cutslice.append(slice(0,len(objsf.dimensions[idimn])))
803	found = True
804	if not found: cutslice.append(0)
805
806	ovalaxisv = ovalaxis[tuple(cutslice)]
807
808	## if len(ovalaxis.shape) <= 2:
809	## ovalaxisv = ovalaxis[:]
810
811	## elif len(ovalaxis.shape) == 3:
812	## ovalaxisv = ovalaxis[0,:]
813	## else:
814	## print errormsg
815	## print ' ' + fname + ': shape of dimension variable:', ovalaxis.shape, \
816	## ' not ready!!'
817	## quit(-1)
818
819	if countlabelfmt == 'None':
820	countlfmt = None
821	else:
822	countlfmt = countlabelfmt
823
824	shading_nx = np.zeros((2), dtype=np.float)
825	shading_nx[0] = np.float(shadminmax.split(',')[0])
826	shading_nx[1] = np.float(shadminmax.split(',')[1])
827
828	clevmin = np.float(contlevels.split(',')[0])
829	clevmax = np.float(contlevels.split(',')[1])
830	Nclevels = int(contlevels.split(',')[2])
831
832	levels_cont = drw.pretty_int(clevmin, clevmax, Nclevels)
833
834	if len(levels_cont) <= 1:
835	print warnmsg
836	print ' ' + fname + ': wrong contour levels:', levels_cont,' !!'
837	del(levels_cont)
838	levels_cont = np.zeros((Nclevels), dtype=np.float)
839	levels_cont = np.arange(7)*(clevmax - clevmin)/(Nclevels-1)
840	print ' generating default ones: ',levels_cont
841
842	if mapvalue == 'None': mapvalue = None
843
844	drw.plot_2D_shadow_contour_time(valshad, valcont, vnamesfig, ovalaxisv, \
845	timevals, timepos, timelabels, ovalu, timelabel, dimt, dimnamesv, colbarn, \
846	countkind, countlfmt, shading_nx, levels_cont, varunits, figtitle, figkind, \
847	revals, mapvalue)
848
849	return
850
851	def draw_2D_shad_line(ncfile, values, varn):
852	""" plotting a fields with shading and another with line
853	draw_2D_shad_line(ncfile, values, varn)
854	ncfile= [ncfiles],[ncfilel] file to use for the shading and for the line
855	values=[vnamefs],[vnamefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar],[colline]:[sminv],[smaxv]:[figt]:
856	[kindfig]:[reverse]:[mapv]:[close]
857	[vnamefs]: Name in the figure of the variable to be shaded
858	[vnamefl]: Name in the figure of the variable to be lined
859	[dimvals]: ',' list of [dimname]\|[value] telling at which dimension of the
860	variable a given value is required (-1, all the length)
861	[dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
862	[colorbar]: name of the color bar
863	[colline]: name of the color for the line
864	[smin/axv]: minimum and maximum value for the shading
865	[figt]: title of the figure ('\|' for spaces)
866	[kindfig]: kind of figure
867	[reverse]: Transformation of the values
868	* 'transpose': reverse the axes (x-->y, y-->x)
869	* 'flip'@[x/y]: flip the axis x or y
870	[mapv]: map characteristics: [proj],[res]
871	see full documentation: http://matplotlib.org/basemap/
872	[proj]: projection
873	* 'cyl', cilindric
874	* 'lcc', lamvbert conformal
875	[res]: resolution:
876	* 'c', crude
877	* 'l', low
878	* 'i', intermediate
879	* 'h', high
880	* 'f', full
881	valules= 'rh:z\|-1,x\|-1:z\|-1,x\|-1:lat:pressure:BuPu:0.,100.:rh:pdf:flip@y:None'
882	varn= [varsn],[varnl] name of the variable to plot with shading and with line
883	"""
884
885	fname = 'draw_2D_shad_line'
886	if values == 'h':
887	print fname + '_____________________________________________________________'
888	print draw_2D_shad_line.__doc__
889	quit()
890
891	farguments = ['[vnamefs],[vnamefl]', '[dimvals]', '[dimxvn]', '[dimyvn]', \
892	'[colorbar],[colline]', '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', \
893	'[mapv]', '[close]']
894	drw.check_arguments(fname,11,values,':',farguments)
895
896	vnamesfig = values.split(':')[0].split(',')[0]
897	dimvals= values.split(':')[1].replace('\|',':')
898	vdimxn = values.split(':')[2]
899	vdimyn = values.split(':')[3]
900	colbarn = values.split(':')[4].split(',')[0]
901	shadminmax = values.split(':')[5]
902	figtitle = values.split(':')[6].replace('\|',' ')
903	figkind = values.split(':')[7]
904	revals = values.split(':')[8]
905	mapvalue = values.split(':')[9]
906	# varn = values.split(':')[10]
907
908	ncfiles = ncfile.split(',')[0]
909
910	if not os.path.isfile(ncfiles):
911	print errormsg
912	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
913	quit(-1)
914
915	objsf = NetCDFFile(ncfiles, 'r')
916
917	varns = varn.split(',')[0]
918
919	if not objsf.variables.has_key(varns):
920	print errormsg
921	print ' ' + fname + ': shading file "' + ncfiles + \
922	'" does not have variable "' + varns + '" !!'
923	quit(-1)
924
925	# Variables' values
926	objvars = objsf.variables[varns]
927
928	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
929
930	# Dimensions names
931	## print fname + ' obj dimnames: ', objvars.dimensions, dimvals, len(dimvals.split(','))
932	## dimnamesv = []
933	## for idd in range(len(objvars.dimensions)):
934	## cutdim = False
935	## for idc in range(len(dimvals.split(','))):
936	## dimcutn = dimvals.split(',')[idc].split(':')[0]
937	## print objvars.dimensions[idd], dimcutn
938	## if objvars.dimensions[idd] == dimcutn:
939	## cutdim = True
940	## break
941	## if not cutdim: dimnamesv.append(objvars.dimensions[idd])
942	dimnamesv = [vdimyn, vdimxn]
943
944	varunits = objvars.getncattr('units')
945
946	if not objsf.variables.has_key(vdimxn):
947	print errormsg
948	print ' ' + fname + ': shading file "' + ncfiles + \
949	'" does not have dimension variable "' + vdimxn + '" !!'
950	quit(-1)
951	if not objsf.variables.has_key(vdimyn):
952	print errormsg
953	print ' ' + fname + ': shading file "' + ncfiles + \
954	'" does not have dimensino variable "' + vdimyn + '" !!'
955	quit(-1)
956
957	objdimx = objsf.variables[vdimxn]
958	objdimy = objsf.variables[vdimyn]
959	odimxu = objdimx.getncattr('units')
960	odimyu = objdimy.getncattr('units')
961
962	if len(objdimx.shape) <= 2:
963	# odimxv = objdimx[valshad.shape]
964	# odimyv = objdimy[valshad.shape]
965	odimxv = objdimx[:]
966	odimyv = objdimy[:]
967
968	elif len(objdimx.shape) == 3:
969	# dimcut = [0, slice(0,valshad.shape[0]), slice(0,valshad.shape[1])]
970	# odimxv = objdimx[tuple(dimcut)]
971	# odimyv = objdimy[tuple(dimcut)]
972	odimxv = objdimx[0,:]
973	odimyv = objdimy[0,:]
974	else:
975	print errormsg
976	print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
977	' not ready!!'
978	quit(-1)
979
980	shading_nx = np.zeros((2), dtype=np.float)
981	shading_nx[0] = np.float(shadminmax.split(',')[0])
982	shading_nx[1] = np.float(shadminmax.split(',')[1])
983
984	if mapvalue == 'None': mapvalue = None
985
986	# line plot
987	##
988	ncfilel = ncfile.split(',')[1]
989	vnamelfig = values.split(':')[0].split(',')[1]
990	varnl = varn.split(',')[1]
991	colline = values.split(':')[4].split(',')[1]
992
993	objlf = NetCDFFile(ncfilel,'r')
994	objlvar = objlf.variables[varnl]
995
996	linevals = objlvar[:]
997	varlunits = objlvar.units
998
999	drw.plot_2D_shadow_line(valshad, linevals, vnamesfig, vnamelfig, odimxv, odimyv, \
1000	odimxu, odimyu, dimnamesv, colbarn, colline, shading_nx, varunits, varlunits, \
1001	figtitle, figkind, revals, mapvalue, True)
1002
1003	objsf.close()
1004	objlf.close()
1005
1006	return
1007
1008	def draw_2D_shad_line_time(ncfile, values, varn):
1009	""" plotting a fields with shading and a line with time values
1010	draw_2D_shad_line(ncfile, values, varn)
1011	ncfile= [ncfiles],[ncfilel] files to use to draw with shading and the line
1012	values= [vnamefs],[vanemefl]:[dimvals]:[dimxvn]:[dimyvn]:[colorbar]:[sminv],[smaxv]:[figt]:
1013	[kindfig]:[reverse]:[timevals]:[close]
1014	[vnamefs]: Name in the figure of the variable to be shaded
1015	[vnamefl]: Name in the figure of the variable to be lined
1016	[dimvals]: ',' list of [dimname]\|[value] telling at which dimension of the
1017	variable a given value is required (-1, all the length)
1018	[dimx/yvn]: name of the variables with the values of the final dimensions (x,y)
1019	[colorbar]: name of the color bar
1020	[smin/axv]: minimum and maximum value for the shading
1021	[figt]: title of the figure ('\|' for spaces)
1022	[kindfig]: kind of figure
1023	[reverse]: Transformation of the values
1024	* 'transpose': reverse the axes (x-->y, y-->x)
1025	* 'flip'@[x/y]: flip the axis x or y
1026	[timevals]: [timen]\|[units]\|[kind]\|[tfmt]\|[label]\|[timeaxis] time labels characteristics
1027	[timen]; name of the time variable
1028	[units]; units string according to CF conventions ([tunits] since
1029	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
1030	[kind]; kind of output
1031	'Nval': according to a given number of values as 'Nval',[Nval]
1032	'exct': according to an exact time unit as 'exct',[tunit];
1033	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1034	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1035	'l': milisecond
1036	[tfmt]; desired format
1037	[label]; label at the graph ('!' for spaces)
1038	[close]: should figure be closed (finished)
1039	values='dtcon,prc:Time\|-1,bottom_top\|-1:presmean:time:seismic:-3.e-6,3.e-6:monthly\|'
1040	'dtcon:pdf:transpose:time\|hours!since!1949-12-01\|exct,2,d\|%d\|date!([DD])\|x:True
1041	varn= [varsn].[varln] name of the variable to plot with shading and to plot with line
1042	"""
1043	fname = 'draw_2D_shad_line_time'
1044	if values == 'h':
1045	print fname + '_____________________________________________________________'
1046	print draw_2D_shad__line_time.__doc__
1047	quit()
1048
1049	farguments = ['[vnamefs],[vanemefl]', '[dimvals]', '[dimxvn]', '[dimyvn]', \
1050	'[colorbar]', '[sminv],[smaxv]', '[figt]', '[kindfig]', '[reverse]', \
1051	'[timevals]', '[close]']
1052	drw.check_arguments(fname,11,values,':',farguments)
1053
1054	vnamesfig = values.split(':')[0].split(',')[0]
1055	dimvals= values.split(':')[1].replace('\|',':')
1056	vdimxn = values.split(':')[2]
1057	vdimyn = values.split(':')[3]
1058	colbarn = values.split(':')[4]
1059	shadminmax = values.split(':')[5]
1060	figtitle = values.split(':')[6].replace('\|',' ')
1061	figkind = values.split(':')[7]
1062	revals = values.split(':')[8]
1063	timevals = values.split(':')[9]
1064	close = values.split(':')[10]
1065
1066	ncfiles = ncfile.split(',')[0]
1067
1068	if not os.path.isfile(ncfiles):
1069	print errormsg
1070	print ' ' + fname + ': shading file "' + ncfiles + '" does not exist !!'
1071	quit(-1)
1072
1073	objsf = NetCDFFile(ncfiles, 'r')
1074
1075	varns = varn.split(',')[0]
1076
1077	if not objsf.variables.has_key(varns):
1078	print errormsg
1079	print ' ' + fname + ': shading file "' + ncfiles + \
1080	'" does not have variable "' + varns + '" !!'
1081	quit(-1)
1082
1083	# Variables' values
1084	objvars = objsf.variables[varns]
1085
1086	valshad, dimsshad = drw.slice_variable(objvars, dimvals.replace(',','\|'))
1087
1088	dimnamesv = [vdimyn, vdimxn]
1089
1090	varunits = objvars.getncattr('units')
1091
1092	if not objsf.variables.has_key(vdimxn):
1093	print errormsg
1094	print ' ' + fname + ': shading file "' + ncfiles + \
1095	'" does not have dimension variable "' + vdimxn + '" !!'
1096	quit(-1)
1097	if not objsf.variables.has_key(vdimyn):
1098	print errormsg
1099	print ' ' + fname + ': shading file "' + ncfiles + \
1100	'" does not have dimensino variable "' + vdimyn + '" !!'
1101	quit(-1)
1102
1103	objdimx = objsf.variables[vdimxn]
1104	objdimy = objsf.variables[vdimyn]
1105	odimxu = objdimx.getncattr('units')
1106	odimyu = objdimy.getncattr('units')
1107
1108	if len(objdimx.shape) <= 2:
1109	odimxv = objdimx[:]
1110	odimyv = objdimy[:]
1111
1112	elif len(objdimx.shape) == 3:
1113	odimxv = objdimx[0,:]
1114	odimyv = objdimy[0,:]
1115	else:
1116	print errormsg
1117	print ' ' + fname + ': shape of dimension variable:', objdimx.shape, \
1118	' not ready!!'
1119	quit(-1)
1120
1121	timename = timevals.split('\|')[0]
1122	timeunit = timevals.split('\|')[1].replace('!',' ')
1123	timekind = timevals.split('\|')[2]
1124	timefmt = timevals.split('\|')[3]
1125	timelabel = timevals.split('\|')[4].replace('!',' ')
1126
1127	if vdimxn == timename:
1128	odimxv = objsf.variables[vdimxn][:]
1129	odimxu = timelabel
1130	timeaxis = 'x'
1131	odimyv = objsf.variables[vdimyn]
1132	odimyu = odimyv.getncattr('units')
1133	timepos, timelabels = drw.CFtimes_plot(odimxv, timeunit, timekind, timefmt)
1134	elif vdimyn == timename:
1135	odimyv = objsf.variables[vdimyn][:]
1136	odimyu = timelabel
1137	timeaxis = 'y'
1138	odimxv = objsf.variables[vdimxn]
1139	odimxu = odimxv.getncattr('units')
1140	timepos, timelabels = drw.CFtimes_plot(odimyv, timeunit, timekind, timefmt)
1141	else:
1142	print errormsg
1143	print ' ' + fname + ": time variable '" + timename + "' not found!!"
1144	quit(-1)
1145
1146	shading_nx = np.zeros((2), dtype=np.float)
1147	shading_nx[0] = np.float(shadminmax.split(',')[0])
1148	shading_nx[1] = np.float(shadminmax.split(',')[1])
1149
1150	closeval = drw.Str_Bool(close)
1151
1152	drw.plot_2D_shadow_time(valshad, vnamesfig, odimxv, odimyv, odimxu, odimyu, \
1153	dimnamesv, colbarn, shading_nx, varunits, figtitle, figkind, revals, timeaxis, \
1154	timepos, timelabels, False)
1155
1156	# Line values
1157	##
1158	ncfilel = ncfile.split(',')[1]
1159
1160	vnamelfig = values.split(':')[0].split(',')[1]
1161	varnl = varn.split(',')[1]
1162
1163	objlf = NetCDFFile(ncfilel,'r')
1164	objlvar = objlf.variables[varnl]
1165
1166	linevals = objlvar[:]
1167	if reva0 == 'tranpose':
1168	plt.plot (linevals, odimxv, '-', color='k')
1169	else:
1170	plt.plot (odimxv, linevals, '-', color='k')
1171
1172	objsf.close()
1173	objsl.close()
1174
1175	return
1176
1177	def draw_topo_geogrid(ncfile, values):
1178	""" plotting geo_em.d[nn].nc topography from WPS files
1179	draw_topo_geogrid(ncfile, values)
1180	ncfile= geo_em.d[nn].nc file to use
1181	values= [minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]
1182	[min/max]Topo: minimum and maximum values of topography to draw
1183	lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] or None
1184	title: title of the graph
1185	graphic_kind: kind of figure (jpg, pdf, png)
1186	mapvalues: map characteristics [proj],[res]
1187	see full documentation: http://matplotlib.org/basemap/
1188	[proj]: projection
1189	* 'cyl', cilindric
1190	* 'lcc', lambert conformal
1191	[res]: resolution:
1192	* 'c', crude
1193	* 'l', low
1194	* 'i', intermediate
1195	* 'h', high
1196	* 'f', full
1197	"""
1198	fname = 'draw_topo_geogrid'
1199
1200	if values == 'h':
1201	print fname + '_____________________________________________________________'
1202	print draw_topo_geogrid.__doc__
1203	quit()
1204
1205	expectargs = ['[minTopo]','[maxTopo]', '[lonlatL]', '[title]', '[graphic_kind]', \
1206	'[mapvalues]']
1207
1208	drw.check_arguments(fname,5,values,':',expectargs)
1209
1210	mintopo = values.split(':')[0].split(',')[0]
1211	maxtopo = values.split(':')[0].split(',')[1]
1212
1213	lonlatLS = values.split(':')[1]
1214	lonlatLv = lonlatLS.split(',')[0]
1215
1216	if lonlatLv == 'None':
1217	lonlatL = None
1218	else:
1219	lonlatL = np.zeros((4), dtype=np.float)
1220	lonlatL[0] = np.float(lonlatLS.split(',')[0])
1221	lonlatL[1] = np.float(lonlatLS.split(',')[1])
1222	lonlatL[2] = np.float(lonlatLS.split(',')[2])
1223	lonlatL[3] = np.float(lonlatLS.split(',')[3])
1224
1225	grtit = values.split(':')[2]
1226	kindfig = values.split(':')[3]
1227	mapvalues = values.split(':')[4]
1228
1229	if not os.path.isfile(ncfile):
1230	print errormsg
1231	print ' ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1232	quit(-1)
1233
1234	objdomf = NetCDFFile(ncfile, 'r')
1235
1236	objhgt = objdomf.variables['HGT_M']
1237	objlon = objdomf.variables['XLONG_M']
1238	objlat = objdomf.variables['XLAT_M']
1239
1240	topography = objhgt[0,:,:]
1241
1242	drw.plot_topo_geogrid(topography, objlon, objlat, mintopo, maxtopo, lonlatL, \
1243	grtit, kindfig, mapvalues, True)
1244
1245	objdomf.close()
1246
1247	return
1248
1249	def draw_topo_geogrid_boxes(ncfiles, values):
1250	""" plotting different geo_em.d[nn].nc topography from WPS files
1251	draw_topo_geogrid_boxes(ncfile, values)
1252	ncfiles= ',' list of geo_em.d[nn].nc files to use (fisrt as topographyc reference)
1253	values= [minTopo],[maxTopo]:[lonlatL]:[title]:[graphic_kind]:[mapvalues]:[labels]
1254	[min/max]Topo: minimum and maximum values of topography to draw
1255	lonlatL: limits of longitudes and latitudes [lonmin, latmin, lonmax, latmax] or None
1256	title: title of the graph
1257	graphic_kind: kind of figure (jpg, pdf, png)
1258	mapvalues: map characteristics [proj],[res]
1259	see full documentation: http://matplotlib.org/basemap/
1260	[proj]: projection
1261	* 'cyl', cilindric
1262	* 'lcc', lambert conformal
1263	[res]: resolution:
1264	* 'c', crude
1265	* 'l', low
1266	* 'i', intermediate
1267	* 'h', high
1268	* 'f', full
1269	labels= labels to write in the graph
1270	"""
1271	# import matplotlib as mpl
1272	# mpl.use('Agg')
1273	import matplotlib.pyplot as plt
1274
1275	fname = 'draw_topo_geogrid_boxes'
1276
1277	if values == 'h':
1278	print fname + '_____________________________________________________________'
1279	print draw_topo_geogrid_boxes.__doc__
1280	quit()
1281
1282	mintopo = values.split(':')[0].split(',')[0]
1283	maxtopo = values.split(':')[0].split(',')[1]
1284
1285	lonlatLS = values.split(':')[1]
1286	lonlatLv = lonlatLS.split(',')[0]
1287
1288	if lonlatLv == 'None':
1289	lonlatL = None
1290	else:
1291	lonlatL = np.zeros((4), dtype=np.float)
1292	lonlatL[0] = np.float(lonlatLS.split(',')[0])
1293	lonlatL[1] = np.float(lonlatLS.split(',')[1])
1294	lonlatL[2] = np.float(lonlatLS.split(',')[2])
1295	lonlatL[3] = np.float(lonlatLS.split(',')[3])
1296
1297	grtit = values.split(':')[2]
1298	kindfig = values.split(':')[3]
1299	mapvalues = values.split(':')[4]
1300	labels = values.split(':')[5]
1301
1302	ncfile = ncfiles.split(',')[0]
1303	if not os.path.isfile(ncfile):
1304	print errormsg
1305	print ' ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1306	quit(-1)
1307
1308	objdomf = NetCDFFile(ncfile, 'r')
1309
1310	objhgt = objdomf.variables['HGT_M']
1311	objlon0 = objdomf.variables['XLONG_M']
1312	objlat0 = objdomf.variables['XLAT_M']
1313
1314	topography = objhgt[0,:,:]
1315
1316	Nfiles = len(ncfiles.split(','))
1317	boxlabels = labels.split(',')
1318
1319	Xboxlines = []
1320	Yboxlines = []
1321
1322	for ifile in range(Nfiles):
1323	ncfile = ncfiles.split(',')[ifile]
1324	# print ifile, ncfile
1325	if not os.path.isfile(ncfile):
1326	print errormsg
1327	print ' ' + fname + ': domain file "' + ncfile + '" does not exist !!'
1328	quit(-1)
1329
1330	objdomfi = NetCDFFile(ncfile, 'r')
1331
1332	objlon = objdomfi.variables['XLONG_M']
1333	objlat = objdomfi.variables['XLAT_M']
1334
1335	dx = objlon.shape[2]
1336	dy = objlon.shape[1]
1337
1338	Xboxlines.append(objlon[0,0,:])
1339	Yboxlines.append(objlat[0,0,:])
1340	Xboxlines.append(objlon[0,dy-1,:])
1341	Yboxlines.append(objlat[0,dy-1,:])
1342	Xboxlines.append(objlon[0,:,0])
1343	Yboxlines.append(objlat[0,:,0])
1344	Xboxlines.append(objlon[0,:,dx-1])
1345	Yboxlines.append(objlat[0,:,dx-1])
1346
1347	objdomfi.close()
1348
1349	drw.plot_topo_geogrid_boxes(topography, Xboxlines, Yboxlines, boxlabels, \
1350	objlon0, objlat0, mintopo, maxtopo, lonlatL, grtit, kindfig, mapvalues, True)
1351
1352	objdomf.close()
1353
1354	return
1355
1356	def movievalslice(origslice, dimmovien, framenum):
1357	""" Function to provide variable slice according to a geneation of a movie
1358	movievals(origslice, dimmovien, framenum)
1359	[origslice]= slice original as [dimname1]\|[val1],[...,[dimnameN]\|[valN]]
1360	([val] = -1, full length)
1361	[dimmovien]= name of the dimension to produce the movie
1362	[framenum]= value of the frame to substitue in [origslice] as
1363	[dimmovien]\|[framenum]
1364	>>> movievalslice('East_West\|-1,North_South\|-1,Time\|2','Time',0)
1365	East_West\|-1,North_South\|-1,Time\|0
1366	"""
1367
1368	fname = 'movievalslice'
1369
1370	if origslice == 'h':
1371	print fname + '_____________________________________________________________'
1372	print movievalslice.__doc__
1373	quit()
1374
1375	dims = origslice.split(',')
1376
1377	movieslice = ''
1378	idim = 0
1379
1380	for dimn in dims:
1381	dn = dimn.split('\|')[0]
1382	if dn == dimmovien:
1383	movieslice = movieslice + dn + '\|' + str(framenum)
1384	else:
1385	movieslice = movieslice + dimn
1386	if idim < len(dims)-1: movieslice = movieslice + ','
1387
1388	idim = idim + 1
1389
1390	return movieslice
1391
1392	class Capturing(list):
1393	""" Class to capture function output as a list
1394	from: http://stackoverflow.com/questions/16571150/how-to-capture-stdout-output-from-a-python-function-call
1395	"""
1396	# from cStringIO import StringIO
1397
1398	def __enter__(self):
1399	self._stdout = sys.stdout
1400	sys.stdout = self._stringio = StringIO()
1401	return self
1402	def __exit__(self, *args):
1403	self.extend(self._stringio.getvalue().splitlines())
1404	sys.stdout = self._stdout
1405
1406	def create_movie(netcdfile, values, variable):
1407	""" Function to create a movie assuming ImageMagick installed!
1408	values= [graph]#[movie_dimension]#[graph_values]
1409	[graph]: which graphic
1410	[movie_dimension]: [dimnmovie]@[dimvmovie]@[moviedelay]@[interval]
1411	[dimnmovie]; name of the dimension from which make the movie
1412	[dimvmovie]; name of the variable with the values of the dimension
1413	[moviedelay]; delay between frames
1414	[interval]; [beg]@[end]@[freq] or -1 (all)
1415	[graph_values]: values to generate the graphic
1416	netcdfile= netCDF file
1417	variable= variable to use (when applicable)
1418	"""
1419	fname = 'create_movie'
1420
1421	if values == 'h':
1422	print fname + '_____________________________________________________________'
1423	print create_movie.__doc__
1424	quit()
1425
1426	graph = values.split('#')[0]
1427	movie_dim = values.split('#')[1]
1428	graph_vals = values.split('#')[2]
1429
1430	ncobj = NetCDFFile(netcdfile, 'r')
1431
1432	# Movie dimension
1433	##
1434	dimnmovie = movie_dim.split('@')[0]
1435	dimvmovie = movie_dim.split('@')[1]
1436	moviedelay = movie_dim.split('@')[2]
1437	moviebeg = int(movie_dim.split('@')[3])
1438
1439	if not drw.searchInlist(ncobj.dimensions.keys(),dimnmovie):
1440	print errormsg
1441	print ' ' + fname + ": file '" + netcdfile + "' has not dimension '" + \
1442	dimnmovie + "' !!!"
1443	quit(-1)
1444
1445	objdmovie = ncobj.dimensions[dimnmovie]
1446	dmovie = len(objdmovie)
1447	if moviebeg != -1:
1448	moviend = int(movie_dim.split('@')[4])
1449	moviefreq = int(movie_dim.split('@')[5])
1450	else:
1451	moviebeg = 0
1452	moviend = dmovie
1453	moviefreq = 1
1454
1455	if dimvmovie == 'WRFTimes':
1456	objvdmovie = ncobj.variables['Times']
1457	vdmovieunits = ''
1458	valsdmovie = []
1459	for it in range(objvdmovie.shape[0]):
1460	valsdmovie.append(drw.datetimeStr_conversion(objvdmovie[it,:], \
1461	'WRFdatetime', 'Y/m/d H-M-S'))
1462	elif dimvmovie == 'CFtime':
1463	objvdmovie = ncobj.variables['time']
1464	vdmovieunits = ''
1465	print objvdmovie.units
1466	valsdmovie0 = drw.netCDFdatetime_realdatetime(objvdmovie.units, 'standard', \
1467	objvdmovie[:])
1468	valsdmovie = []
1469	for it in range(objvdmovie.shape[0]):
1470	valsdmovie.append(drw.datetimeStr_conversion(valsdmovie0[it,:], \
1471	'matYmdHMS', 'Y/m/d H-M-S'))
1472	else:
1473	if not drw.searchInlist(ncobj.variables.keys(),dimvmovie):
1474	print errormsg
1475	print ' ' + fname + ": file '" + netcdfile + "' has not variable '" + \
1476	dimvmovie + "' !!!"
1477	quit(-1)
1478	vdmovieunits = objvdmovie.getncattr('units')
1479	objvdmovie = ncobj.variables[dimvmovie]
1480	if len(objvdmovie.shape) == 1:
1481	vasldmovie = objvdmovie[:]
1482	else:
1483	print errormsg
1484	print ' ' + fname + ': shape', objvdmovie.shape, 'of variable with ' + \
1485	'dimension movie values not ready!!!'
1486	quit(-1)
1487
1488	ncobj.close()
1489	os.system('rm frame_*.png > /dev/null')
1490
1491	# graphic
1492	##
1493	if graph == 'draw_2D_shad':
1494	graphvals = graph_vals.split(':')
1495
1496	for iframe in range(moviebeg,moviend,moviefreq):
1497	iframeS = str(iframe).zfill(4)
1498
1499	drw.percendone((iframe-moviebeg)/moviefreq,(moviend-moviebeg)/moviefreq, \
1500	5, 'frames')
1501	titgraph = dimnmovie + '\|=\|' + str(valsdmovie[iframe]) + '\|' + \
1502	vdmovieunits
1503
1504	graphvals[1] = movievalslice(graphvals[1],dimnmovie,iframe)
1505	graphvals[6] = titgraph
1506	graphvals[7] = 'png'
1507
1508	graphv = drw.numVector_String(graphvals, ":")
1509
1510	with Capturing() as output:
1511	draw_2D_shad(netcdfile, graphv, variable)
1512
1513	os.system('mv 2Dfields_shadow.png frame_' + iframeS + '.png')
1514	else:
1515	print errormsg
1516	print ' ' + fname + ": graphic '" + graph + "' not defined !!!"
1517	quit(-1)
1518
1519	os.system('convert -delay ' + moviedelay + ' -loop 0 frame_*.png create_movie.gif')
1520
1521	print "Succesfuly creation of movie file 'create_movie.gif' !!!"
1522
1523	return
1524
1525	def draw_lines(ncfilens, values, varname):
1526	""" Function to draw different lines at the same time from different files
1527	draw_lines(ncfilens, values, varname):
1528	ncfilens= [filen] ',' separated list of netCDF files
1529	values= [dimvname]:[valuesaxis]:[dimtit]:[leglabels]:[vtit]:[title]:[locleg]:[graphk]
1530	[dimvname]: ',' list of names of the variable with he values of the common dimension
1531	[valuesaxis]: which axis will be used for the values ('x', or 'y')
1532	[dimtit]: title for the common dimension
1533	[leglabels]: ',' separated list of names for the legend
1534	[vartit]: name of the variable in the graph
1535	[title]: title of the plot ('\|' for spaces)
1536	[locleg]: location of the legend (-1, autmoatic)
1537	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
1538	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
1539	9: 'upper center', 10: 'center'
1540	[graphk]: kind of the graphic
1541	varname= variable to plot
1542	values= 'XLAT:x:latitude:32x32:$wss^{*}$:wss Taylor's turbulence term:pdf'
1543	"""
1544
1545	fname = 'draw_lines'
1546
1547	if values == 'h':
1548	print fname + '_____________________________________________________________'
1549	print draw_lines.__doc__
1550	quit()
1551
1552	expectargs = '[dimvname]:[valuesaxis]:[dimtit]:[leglabels]:[vtit]:[title]:[locleg]:[graphk]'
1553	drw.check_arguments(fname,len(expectargs.split(':')),values,':',expectargs)
1554
1555	ncfiles = ncfilens.split(',')
1556	dimvnames = values.split(':')[0]
1557	valuesaxis = values.split(':')[1]
1558	dimtit = values.split(':')[2]
1559	leglabels = values.split(':')[3].replace('_','\_')
1560	vartit = values.split(':')[4]
1561	title = values.split(':')[5].replace('\|',' ')
1562	locleg = values.split(':')[6]
1563	graphk = values.split(':')[7]
1564
1565	Nfiles = len(ncfiles)
1566
1567	# Getting trajectotries
1568	##
1569
1570	varvalues = []
1571	dimvalues = []
1572
1573	print ' ' + fname
1574	ifn = 0
1575	for ifile in ncfiles:
1576	filen = ifile.split('@')[0]
1577
1578	print ' filen:',filen
1579
1580	if not os.path.isfile(filen):
1581	print errormsg
1582	print ' ' + fname + ": netCDF file '" + filen + "' does not exist !!"
1583	quit(-1)
1584
1585	objfile = NetCDFFile(filen, 'r')
1586
1587	if dimvnames.find(',') != -1:
1588	dimvname = dimvnames.split(',')
1589	else:
1590	dimvname = [dimvnames]
1591
1592	found = False
1593	for dvn in dimvname:
1594	if objfile.variables.has_key(dvn):
1595	found = True
1596	break
1597
1598	if not found:
1599	print errormsg
1600	print ' ' + fname + ": netCDF file '" + filen + \
1601	"' does not have variables '" + dimvnames + "' !!"
1602	quit(-1)
1603
1604	if not objfile.variables.has_key(varname):
1605	print errormsg
1606	print ' ' + fname + ": netCDF file '" + filen + \
1607	"' does not have variable '" + varname + "' !!"
1608	quit(-1)
1609
1610	vvobj = objfile.variables[varname]
1611	if len(vvobj.shape) != 1:
1612	print errormsg
1613	print ' ' + fname + ': wrong shape:',vvobj.shape," of variable '" + \
1614	varname + "' !!"
1615	quit(-1)
1616
1617	for dimvn in dimvname:
1618	if drw.searchInlist(objfile.variables, dimvn):
1619	vdobj = objfile.variables[dimvn]
1620	if len(vdobj.shape) != 1:
1621	print errormsg
1622	print ' ' + fname + ': wrong shape:',vdobj.shape, \
1623	" of variable '" + dimvn + "' !!"
1624	quit(-1)
1625	break
1626
1627	varvalues.append(vvobj[:])
1628	dimvalues.append(vdobj[:])
1629
1630	if ifn == 0:
1631	varunits = vvobj.units
1632
1633	objfile.close()
1634
1635	ifn = ifn + 1
1636
1637	drw.plot_lines(dimvalues, varvalues, valuesaxis, dimtit, leglabels.split(','), \
1638	vartit, varunits, title, locleg, graphk)
1639
1640	return
1641
1642	def draw_lines_time(ncfilens, values, varname):
1643	""" Function to draw different lines at the same time from different files with times
1644	draw_lines_time(ncfilens, values, varname):
1645	ncfilens= [filen] ',' separated list of netCDF files
1646	values= [dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];[timevals];[locleg];[graphk]
1647	[dimvname]: name of the variable with he values of the common dimension
1648	[valuesaxis]: which axis will be used for the values ('x', or 'y')
1649	[dimtit]: title for the common dimension
1650	[leglabels]: ',' separated list of names for the legend
1651	[vartit]: name of the variable in the graph
1652	[title]: title of the plot ('\|' for spaces)
1653	[timevals]: [timen]\|[units]\|[kind]\|[tfmt] time labels characteristics
1654	[timen]; name of the time variable
1655	[units]; units string according to CF conventions ([tunits] since
1656	[YYYY]-[MM]-[DD] [[HH]:[MI]:[SS]], '!' for spaces)
1657	[kind]; kind of output
1658	'Nval': according to a given number of values as 'Nval',[Nval]
1659	'exct': according to an exact time unit as 'exct',[tunit];
1660	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1661	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1662	'l': milisecond
1663	[tfmt]; desired format
1664	[locleg]: location of the legend (-1, autmoatic)
1665	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
1666	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
1667	9: 'upper center', 10: 'center'
1668	[graphk]: kind of the graphic
1669	varname= variable to plot
1670	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'
1671	"""
1672
1673	fname = 'draw_lines_time'
1674
1675	if values == 'h':
1676	print fname + '_____________________________________________________________'
1677	print draw_lines_time.__doc__
1678	quit()
1679
1680	expectargs = '[dimvname];[valuesaxis];[dimtit];[leglabels];[vtit];[title];'
1681	expectargs = expectargs + '[timevals];[locleg];[graphk]'
1682	drw.check_arguments(fname,len(expectargs.split(';')),values,';',expectargs)
1683
1684	ncfiles = ncfilens.split(',')
1685	dimvname = values.split(';')[0]
1686	valuesaxis = values.split(';')[1]
1687	dimtit = values.split(';')[2]
1688	leglabels = values.split(';')[3].replace('_','\_')
1689	vartit = values.split(';')[4]
1690	title = values.split(';')[5].replace('\|',' ')
1691	timevals = values.split(';')[6]
1692	locleg = int(values.split(';')[7])
1693	graphk = values.split(';')[8]
1694
1695	Nfiles = len(ncfiles)
1696
1697	# Getting values
1698	##
1699	varvalues = []
1700	dimvalues = []
1701	timvalues = []
1702	timvals0 = timvalues
1703
1704	print ' ' + fname
1705	ifn = 0
1706	mintval = 1.e20
1707	maxtval = -1.e20
1708
1709	for ifile in ncfiles:
1710	filen = ifile.split('@')[0]
1711
1712	print ' filen:',filen
1713
1714	if not os.path.isfile(filen):
1715	print errormsg
1716	print ' ' + fname + ": netCDF file '" + filen + "' does not exist !!"
1717	quit(-1)
1718
1719	objfile = NetCDFFile(filen, 'r')
1720
1721	if not objfile.variables.has_key(dimvname):
1722	print errormsg
1723	print ' ' + fname + ": netCDF file '" + filen + \
1724	"' does not have variable '" + dimvname + "' !!"
1725	quit(-1)
1726
1727	if not objfile.variables.has_key(varname):
1728	print errormsg
1729	print ' ' + fname + ": netCDF file '" + filen + \
1730	"' does not have variable '" + varname + "' !!"
1731	quit(-1)
1732
1733	vvobj = objfile.variables[varname]
1734	if len(vvobj.shape) != 1:
1735	print errormsg
1736	print ' ' + fname + ': wrong shape:',vvobj.shape," of variable '" + \
1737	varname + "' !!"
1738	quit(-1)
1739
1740	vdobj = objfile.variables[dimvname]
1741	if len(vdobj.shape) != 1:
1742	print errormsg
1743	print ' ' + fname + ': wrong shape:',vdobj.shape," of variable '" + \
1744	dimvname + "' !!"
1745	quit(-1)
1746
1747	varvalues.append(vvobj[:])
1748	dimvalues.append(vdobj[:])
1749
1750	mindvals = np.min(vdobj[:])
1751	maxdvals = np.max(vdobj[:])
1752
1753	if mindvals < mintval: mintval = mindvals
1754	if maxdvals > maxtval: maxtval = maxdvals
1755
1756	if ifn == 0:
1757	varunits = drw.units_lunits(vvobj.units)
1758
1759	objfile.close()
1760
1761	ifn = ifn + 1
1762
1763	# Times
1764	timename = timevals.split('\|')[0]
1765	timeunit = timevals.split('\|')[1].replace('!',' ')
1766	timekind = timevals.split('\|')[2]
1767	timefmt = timevals.split('\|')[3]
1768
1769	dtvals = (maxtval - mintval)/5.
1770	tvals = np.arange(mintval, maxtval, dtvals/2.)
1771
1772	timepos, timelabels = drw.CFtimes_plot(tvals, timeunit, timekind, timefmt)
1773
1774	drw.plot_lines_time(dimvalues, varvalues, valuesaxis, dimtit, leglabels.split(','), \
1775	vartit, varunits, timepos, timelabels, title, locleg, graphk)
1776
1777	return
1778
1779	def draw_Neighbourghood_evol(filen, values, variable):
1780	""" Function to draw the temporal evolution of a neighbourghood around a point
1781	draw_Neighbourghood_evol(filen, values, variable)
1782	filen= netCDF file name
1783	values= [gvarname]:[dimsval]:[neigdims]:[Nneig]:[Ncol]:[timetits]:[tkinds]:
1784	[timefmts]:[gtitle]:[shadxtrms]:[cbar]:[gkind]:[ofile]
1785	[dimsval]: [dimn1]\|[val1]\|[dimv1],...,[dimnN]\|[valN]\|[dimvN] dimension names, values to get
1786	(-1, for all; no name/value pair given full length) and variable with values of the dimension
1787	NOTE: when dimsval[X,Y] == neigdims[X,Y], valX,valY --> valX,valY-Nneig/2, valX,valY+Nneig/2
1788	[neigdims]: [dimnX],[dimnY] dimensions mnames along which the neigbourghood should be defined
1789	[Nneig]: Number of grid points of the full side of the box (odd value)
1790	[Ncol]: Number of columns ('auto': square final plot)
1791	[gvarname]: name of the variable to appear in the graph
1792	[timetits]: [titX],[titY] titles of the axes ('\|' for spaces)
1793	[tkinds]: [tkindX]\|[tkindY] kinds of time to appear in the graph
1794	'Nval': according to a given number of values as 'Nval',[Nval]
1795	'exct': according to an exact time unit as 'exct',[tunit];
1796	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
1797	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
1798	'l': milisecond
1799	[timefmts]: [tfmtX],[tfmtY] format of the time labels
1800	[gtitle]: title of the graphic ('\|' for spaces)
1801	[shadxtrms]: Extremes for the shading
1802	[cbar]: colorbar to use
1803	[gkind]: kind of graphical output
1804	[ofile]: True/False whether the netcdf with data should be created or not
1805	variable= name of the variable
1806	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'
1807	"""
1808
1809	fname = 'draw_Neighbourghood_evol'
1810
1811	if values == 'h':
1812	print fname + '_____________________________________________________________'
1813	print draw_Neighbourghood_evol.__doc__
1814	quit()
1815
1816	expectargs = ['[gvarname]', '[dimsval]', '[neigdims]', '[Nneig]', '[Ncol]', \
1817	'[timetits]', '[tkinds]', '[timefmts]', '[gtitle]', '[shadxtrms]', '[cbar]', \
1818	'[gkind]', '[ofile]']
1819
1820	drw.check_arguments(fname,len(expectargs),values,':',expectargs)
1821
1822	gvarname = values.split(':')[0]
1823	dimsval = values.split(':')[1].split(',')
1824	neigdims = values.split(':')[2].split(',')
1825	Nneig = int(values.split(':')[3])
1826	Ncol0 = values.split(':')[4]
1827	timetits = values.split(':')[5].split(',')
1828	timekinds = values.split(':')[6].split('\|')
1829	timefmts = values.split(':')[7].split(',')
1830	gtitle = values.split(':')[8].replace('\|',' ')
1831	shadxtrms = values.split(':')[9].split(',')
1832	cbar = values.split(':')[10]
1833	gkind = values.split(':')[11]
1834	ofile = values.split(':')[12]
1835
1836	if Ncol0 != 'auto':
1837	Ncol = int(Ncol0)
1838	else:
1839	Ncol = Ncol0
1840
1841	timetits[0] = timetits[0].replace('\|',' ')
1842	timetits[1] = timetits[1].replace('\|',' ')
1843
1844	if np.mod(Nneig,2) == 0:
1845	print errormsg
1846	print ' ' + fname + ": an odd value for 'Nneig':",Nneig,'is required !!!'
1847	quit(-1)
1848
1849	Nneig2 = int(Nneig/2)
1850
1851	# Values to slice the variable
1852	dimvslice = {}
1853	dimvvalues = {}
1854	for dimvs in dimsval:
1855	dimn = dimvs.split('\|')[0]
1856	dimv = int(dimvs.split('\|')[1])
1857	dimnv = dimvs.split('\|')[2]
1858
1859	dimvvalues[dimn] = dimnv
1860	dimvslice[dimn] = dimv
1861
1862	ncobj = NetCDFFile(filen, 'r')
1863
1864	varobj = ncobj.variables[variable]
1865
1866	slicevar = []
1867	newdimn = []
1868	newdimsvar = {}
1869
1870	for dimn in varobj.dimensions:
1871	if not drw.searchInlist(dimvslice.keys(), dimn):
1872	dimsize = len(ncobj.dimensions[dimn])
1873	slicevar.append(slice(0, dimsize+1))
1874	newdimn.append(dimn)
1875	newdimsvar[dimn] = dimseize
1876
1877	for dimslicen in dimvslice.keys():
1878	if dimn == dimslicen:
1879	if dimvslice[dimn] != -1:
1880	if drw.searchInlist(neigdims, dimn):
1881	slicevar.append(slice(dimvslice[dimn]-Nneig2, \
1882	dimvslice[dimn]+Nneig2+1))
1883	newdimn.append(dimn)
1884	newdimsvar[dimn] = Nneig
1885	break
1886	else:
1887	slicevar.append(slice(dimvslice[dimn], dimvslice[dimn]+1))
1888	break
1889	else:
1890	dimsize = len(ncobj.dimensions[dimn])
1891	slicevar.append(slice(0, dimsize+1))
1892	newdimn.append(dimn)
1893	newdimsvar[dimn] = dimsize
1894	break
1895
1896	varv = varobj[tuple(slicevar)]
1897
1898	if len(newdimn) != 3:
1899	print errormsg
1900	print ' ' + fname + ': sliced variable with shape=', varv.shape, \
1901	' must have three dimensions',len(varv.shape),'given !!'
1902	quit(-1)
1903
1904	newdims = []
1905	for nwdims in newdimn:
1906	newdims.append(newdimsvar[nwdims])
1907
1908	# The dimension which is not in the neighbourhood dimensions must be time!
1909	for dim1 in newdimn:
1910	if not drw.searchInlist(neigdims, dim1):
1911	dimt = newdimsvar[dim1]
1912	dimtime = dim1
1913
1914	if Ncol == 'auto':
1915	dimtsqx = int(np.sqrt(dimt)) + 1
1916	dimtsqy = int(np.sqrt(dimt)) + 1
1917	else:
1918	dimtsqx = int(Ncol)
1919	dimtsqy = dimt/dimtsqx + 1
1920
1921	neighbourghood = np.ones((dimtsqyNneig,dimtsqxNneig), dtype=np.float)*fillValue
1922
1923	for it in range(dimt):
1924	ity = int(it/dimtsqx)
1925	itx = it-ity*dimtsqx
1926
1927	itty = (dimtsqy - ity - 1)*Nneig + Nneig2
1928	ittx = itx*Nneig + Nneig2
1929
1930	neighbourghood[itty-Nneig2:itty+Nneig2+1,ittx-Nneig2:ittx+Nneig2+1]= \
1931	varv[it,::-1,:]
1932
1933	variablevals = drw.variables_values(variable)
1934	if drw.searchInlist(varobj.ncattrs(), 'units'):
1935	vunits = varobj.units
1936	else:
1937	vunits = variablevals[5]
1938
1939	# Time values at the X/Y axes
1940	if ncobj.variables[dimvvalues[dimtime]].dtype == '\|S1':
1941	print ' ' + fname + ': WRF time variable!'
1942	refdate = '19491201000000'
1943	tunitsval = 'hours'
1944	dimtvalues = np.zeros((dimt), dtype=np.float)
1945	tvals = ncobj.variables[dimvvalues[dimtime]]
1946	yrref=refdate[0:4]
1947	monref=refdate[4:6]
1948	dayref=refdate[6:8]
1949	horref=refdate[8:10]
1950	minref=refdate[10:12]
1951	secref=refdate[12:14]
1952
1953	refdateS = yrref + '/' + monref + '/' + dayref + '_' + horref + ':' + \
1954	minref + ':' + secref
1955	tunits = tunitsval + ' since ' + refdateS
1956	for it in range(dimt):
1957	wrfdates = drw.datetimeStr_conversion(tvals[it,:],'WRFdatetime', 'matYmdHMS')
1958	dimtvalues[it] = drw.realdatetime1_CFcompilant(wrfdates, refdate, tunitsval)
1959	else:
1960	dimtvalues = ncobj.variables[dimvvalues[dimtime]][:]
1961	tunits = ncobj.variables[newdimsvar[dimtime]].units
1962
1963	dimxv = dimtvalues[0:dimtsqx]
1964	dimyv = dimtvalues[0:dimt:dimtsqx]
1965
1966	dimn = ['time','time']
1967
1968	if ofile == 'True':
1969	ofilen = 'Neighbourghood_evol.nc'
1970	newnc = NetCDFFile(ofilen, 'w')
1971	# Dimensions
1972	newdim = newnc.createDimension('time',None)
1973	newdim = newnc.createDimension('y',dimtsqy*Nneig)
1974	newdim = newnc.createDimension('x',dimtsqx*Nneig)
1975	# Dimension values
1976	newvar = newnc.createVariable('time','f8',('time'))
1977	newvar[:] = np.arange(dimt)
1978	newattr = drw.basicvardef(newvar, 'time','time',tunits)
1979	# Neighbourhghood variable
1980	newvar = newnc.createVariable(variable + 'neigevol', 'f4', ('y','x'), \
1981	fill_value=fillValue)
1982	newvar[:] = neighbourghood
1983
1984	newnc.sync()
1985	newnc.close()
1986	print fname + ": Successfull generation of file '" + ofilen + "' !!"
1987
1988	# Time ticks
1989	timeposX, timelabelsX = drw.CFtimes_plot(dimxv, tunits, timekinds[0], timefmts[0])
1990	timeposY, timelabelsY = drw.CFtimes_plot(dimyv, tunits, timekinds[1], timefmts[1])
1991
1992	timepos = [timeposX[0:len(timeposX)], timeposY[len(timeposY):0:-1]]
1993	timelabels = [timelabelsX[0:len(timeposX)], timelabelsY[0:len(timeposY)]]
1994
1995	for i in range(2):
1996	if shadxtrms[i][0:1] != 'S':
1997	shadxtrms[i] = np.float(shadxtrms[i])
1998
1999	drw.plot_Neighbourghood_evol(neighbourghood, dimxv, dimyv, gvarname, timetits, \
2000	timepos, timelabels, cbar, Nneig, shadxtrms, vunits, gtitle, gkind, True)
2001
2002	def draw_timeSeries(filen, values, variables):
2003	""" Function to draw a time-series
2004	draw_timeSeries(filen, values, variable):
2005	filen= name of the file
2006	values= [gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[locleg]:[gkind]
2007	[gvarname]: name of the variable to appear in the graph
2008	[timetit]: title of the time axis (assumed x-axis, '\|' for spaces)
2009	[tkind]: kind of time to appear in the graph (assumed x-axis)
2010	'Nval': according to a given number of values as 'Nval',[Nval]
2011	'exct': according to an exact time unit as 'exct',[tunit];
2012	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2013	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2014	'l': milisecond
2015	[timefmt]: format of the time labels
2016	[title]: title of the graphic ('\|' for spaces)
2017	[locleg]: location of the legend (-1, autmoatic)
2018	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2019	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2020	9: 'upper center', 10: 'center'
2021	[gkind]: kind of graphical output
2022	variables= [varname],[timename] names of variable and variable with times
2023	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')
2024	"""
2025
2026	fname = 'draw_timeSeries'
2027
2028	if values == 'h':
2029	print fname + '_____________________________________________________________'
2030	print draw_timeSeries.__doc__
2031	quit()
2032
2033	expectargs = ['[gvarname]', '[timetit]', '[tkind]', '[timefmt]', '[title]', \
2034	'[locleg]', '[gkind]']
2035
2036	drw.check_arguments(fname,len(expectargs),values,':',expectargs)
2037
2038	gvarname = values.split(':')[0]
2039	timetit = values.split(':')[1].replace('\|',' ')
2040	tkind = values.split(':')[2]
2041	timefmt = values.split(':')[3]
2042	title = values.split(':')[4].replace('\|',' ')
2043	locleg = int(values.split(':')[5])
2044	gkind = values.split(':')[6]
2045
2046	ncobj = NetCDFFile(filen, 'r')
2047
2048	variable = variables.split(',')[0]
2049	timevar = variables.split(',')[1]
2050
2051	if not ncobj.variables.has_key(variable):
2052	print errormsg
2053	print ' ' + fname + ": file '" + filen + "' does not have variable '" + \
2054	variable + "' !!"
2055	quit(-1)
2056
2057	if not ncobj.variables.has_key(timevar):
2058	print errormsg
2059	print ' ' + fname + ": file '" + filen + "' does not have variable time '" \
2060	+ timevar + "' !!"
2061	quit(-1)
2062
2063	varobj = ncobj.variables[variable]
2064	timeobj = ncobj.variables[timevar]
2065
2066	dimt = len(timeobj[:])
2067	varvals = np.zeros((2,dimt), dtype=np.float)
2068
2069	gunits = varobj.getncattr('units')
2070	tunits = timeobj.getncattr('units')
2071
2072	varvals[0,:], valpot, newgunits, Spot = drw.pot_values(varobj[:].flatten(), gunits)
2073	varvals[1,:] = timeobj[:]
2074
2075	tseriesvals = []
2076	tseriesvals.append(varvals)
2077
2078	drw.plot_TimeSeries(tseriesvals, Spot + drw.units_lunits(gunits), tunits, \
2079	'TimeSeries', gvarname, timetit, tkind, timefmt, title, \
2080	gvarname.replace('_','\_'), locleg, gkind)
2081
2082	return
2083
2084	#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')
2085
2086	def draw_trajectories(trjfilens, values, observations):
2087	""" Function to draw different trajectories at the same time
2088	draw_trajectories(trjfilens, values, observations):
2089	trjfilens= [filen]@[Tint]@[map] ',' separated list of files with trajectories,
2090	time intervals and reference maps (first one will be used to plot)
2091	[filen]: name of the file to use (lines with '#', not readed) as:
2092	[t-step] [x] [y]
2093	[Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2094	[map]: [file]#[lonname]#[latname]
2095	[file]; with the [lon],[lat] matrices
2096	[lonname],[latname]; names of the longitudes and latitudes variables
2097	values=[leglabels]\|[lonlatlims]\|[title]\|[graphk]\|[mapkind]
2098	[leglabels]: ',' separated list of names for the legend
2099	[lonlatlims]: ',' list of limits of the map [lonmin, latmin, lonmax, latmax] or None
2100	[title]: title of the plot ('!' for spaces)
2101	[graphk]: kind of the graphic
2102	[mapkind]: drawing coastaline ([proj],[res]) or None
2103	[proj]: projection
2104	* 'cyl', cilindric
2105	* 'lcc', lambert conformal
2106	[res]: resolution:
2107	* 'c', crude
2108	* 'l', low
2109	* 'i', intermediate
2110	* 'h', high
2111	* 'f', full
2112	obsevations= [obsfile],[obsname],[Tint],[null]
2113	[obsfile]: name fo the File with the observations as [t-step] [lat] [lon]
2114	[obsname]: name of the observations in the graph
2115	[Tint]: interval of time as [Tbeg]@[Tend] or -1 for all the interval
2116	[null]: null value for the observed trajectory
2117	"""
2118
2119	fname = 'draw_trajectories'
2120
2121	if values == 'h':
2122	print fname + '_____________________________________________________________'
2123	print draw_trajectories.__doc__
2124	quit()
2125
2126	expectargs = '[leglabels]\|[lonlatlims]\|[title]\|[graphk]\|[mapkind]'
2127
2128	drw.check_arguments(fname,len(expectargs.split('\|')),values,'\|',expectargs)
2129
2130	trjfiles = trjfilens.split(',')
2131	leglabels = values.split('\|')[0]
2132	lonlatlims = values.split('\|')[1]
2133	title = values.split('\|')[2].replace('!',' ')
2134	graphk = values.split('\|')[3]
2135	mapkind = values.split('\|')[4]
2136
2137	Nfiles = len(trjfiles)
2138
2139	# Getting trajectotries
2140	##
2141
2142	lontrjvalues = []
2143	lattrjvalues = []
2144
2145	print ' ' + fname
2146	ifn = 0
2147	for ifile in trjfiles:
2148	filen = ifile.split('@')[0]
2149	Tint = ifile.split('@')[1]
2150
2151	print ' trajectory:',filen
2152
2153	if Tint != '-1':
2154	Tbeg = Tint
2155	Tend = ifile.split('@')[2]
2156	mapv = ifile.split('@')[3]
2157	else:
2158	mapv = ifile.split('@')[2]
2159
2160	if not os.path.isfile(filen):
2161	print errormsg
2162	print ' ' + fname + ": trajectory file '" + filen + "' does not exist !!"
2163	quit(-1)
2164
2165	# Charging longitude and latitude values
2166	##
2167	lonvals, latvals = drw.lonlat_values(mapv.split('#')[0], mapv.split('#')[1], \
2168	mapv.split('#')[2])
2169
2170	if ifn == 0: mapref = mapv
2171	ifn = ifn + 1
2172
2173	objfile = open(filen, 'r')
2174	trjtimev = []
2175	trjxv = []
2176	trjyv = []
2177
2178	for line in objfile:
2179	if line[0:1] != '#':
2180	trjtimev.append(int(line.split(' ')[0]))
2181	trjxv.append(int(line.split(' ')[1]))
2182	trjyv.append(int(line.split(' ')[2]))
2183
2184	objfile.close()
2185
2186	if Tint != '-1':
2187	lontrjvalues.append(lonvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
2188	lattrjvalues.append(latvals[trjyv[Tint:Tend+1], trjxv[Tint:Tend+1]])
2189	else:
2190	lontrjvalues.append(lonvals[trjyv[:], trjxv[:]])
2191	lattrjvalues.append(latvals[trjyv[:], trjxv[:]])
2192
2193	# lonlatlimits
2194	##
2195
2196	if lonlatlims == 'None':
2197	lonlatlimsv = None
2198	else:
2199	lonlatlimsv = np.zeros((4), dtype=np.float)
2200	lonlatlimsv[0] = np.float(lonlatlims.split(',')[0])
2201	lonlatlimsv[1] = np.float(lonlatlims.split(',')[1])
2202	lonlatlimsv[2] = np.float(lonlatlims.split(',')[2])
2203	lonlatlimsv[3] = np.float(lonlatlims.split(',')[3])
2204
2205	# lon/lat objects
2206	##
2207	objnc = NetCDFFile(mapref.split('#')[0])
2208	lonobj = objnc.variables[mapref.split('#')[1]]
2209	latobj = objnc.variables[mapref.split('#')[2]]
2210
2211	# map
2212	##
2213	if mapkind == 'None':
2214	mapkindv = None
2215	else:
2216	mapkindv = mapkind
2217
2218	if observations is None:
2219	obsname = None
2220	else:
2221	obsfile = observations.split(',')[0]
2222	obsname = observations.split(',')[1]
2223	Tint = observations.split(',')[2]
2224	null = np.float(observations.split(',')[3])
2225	print ' observational trajectory:',obsfile
2226
2227	if not os.path.isfile(obsfile):
2228	print errormsg
2229	print ' ' + fname + ": observations trajectory file '" + obsfile + \
2230	"' does not exist !!"
2231	quit(-1)
2232
2233	objfile = open(obsfile, 'r')
2234	obstrjtimev = []
2235	obstrjxv = []
2236	obstrjyv = []
2237
2238	for line in objfile:
2239	if line[0:1] != '#':
2240	lon = np.float(line.split(' ')[2])
2241	lat = np.float(line.split(' ')[1])
2242	if not lon == null and not lat == null:
2243	obstrjtimev.append(int(line.split(' ')[0]))
2244	obstrjxv.append(lon)
2245	obstrjyv.append(lat)
2246	else:
2247	obstrjtimev.append(int(line.split(' ')[0]))
2248	obstrjxv.append(None)
2249	obstrjyv.append(None)
2250
2251	objfile.close()
2252
2253	if Tint != '-1':
2254	Tint = int(observations.split(',')[2].split('@')[0])
2255	Tend = int(observations.split(',')[2].split('@')[1])
2256	lontrjvalues.append(obstrjxv[Tint:Tend+1])
2257	lattrjvalues.append(obstrjyv[Tint:Tend+1])
2258	else:
2259	lontrjvalues.append(obstrjxv[:])
2260	lattrjvalues.append(obstrjyv[:])
2261
2262	drw.plot_Trajectories(lontrjvalues, lattrjvalues, leglabels.split(','), \
2263	lonobj, latobj, lonlatlimsv, title, graphk, mapkindv, obsname)
2264
2265	objnc.close()
2266
2267	return
2268
2269	def draw_vals_trajectories(ncfile, values, variable):
2270	""" Function to draw values from the outputs from 'compute_tevolboxtraj'
2271	draw_vals_trajectories(ncfile, values, variable)
2272	ncfile= [ncfile] ',' list of files to use
2273	values= [statisticskind]:[Tint]:[labels]@[locleg]:[gvarname]:[timetit]:[tkind]:[timefmt]:[title]:[gkind]
2274	[statisticskind]=[statistics][kind]
2275	[statistics]: which statistics to use, from: 'center', 'min', 'max', 'mean',
2276	'mean2', 'stdev'
2277	[kind]: 'box', 'circle' statistics taking the values from a box or a circle
2278	'trj': value following the trajectory
2279	[Tint]: [Tbeg]@[Tend] or None, interval of time to plot or -1 for all the times
2280	[labels]: ',' separated list of labels for the legend
2281	[locleg]: location of the legend (-1, autmoatic)
2282	1: 'upper right', 2: 'upper left', 3: 'lower left', 4: 'lower right',
2283	5: 'right', 6: 'center left', 7: 'center right', 8: 'lower center',
2284	9: 'upper center', 10: 'center'
2285	[gvarname]: name of the variable to appear in the graph
2286	[timetit]: title of the time axis (assumed x-axis, '\|' for spaces)
2287	[tkind]: kind of time to appear in the graph (assumed x-axis)
2288	'Nval': according to a given number of values as 'Nval',[Nval]
2289	'exct': according to an exact time unit as 'exct',[tunit];
2290	tunit= [Nunits],[tu]; [tu]= 'c': centuries, 'y': year, 'm': month,
2291	'w': week, 'd': day, 'h': hour, 'i': minute, 's': second,
2292	'l': milisecond
2293	[timefmt]: format of the time labels
2294	[title]: title of the graphic ('\|' for spaces)
2295	[gkind]: kind of graphical output
2296	variable= variable to use
2297	"""
2298
2299	fname = 'draw_vals_trajectories'
2300
2301	if values == 'h':
2302	print fname + '_____________________________________________________________'
2303	print draw_vals_trajectories.__doc__
2304	quit()
2305
2306	sims = ncfile.split(',')
2307
2308	if len(values.split(':')) != 9:
2309	print errormsg
2310	print ' ' + fname + ': wrong number of values!', len(values.split(':')), \
2311	'given 9 needed!!'
2312	print ' ',values.split(':')
2313	quit(-1)
2314
2315	statisticskind = values.split(':')[0]
2316	Tint = values.split(':')[1]
2317	labels = values.split(':')[2]
2318	gvarname = values.split(':')[3]
2319	timetit = values.split(':')[4].replace('\|',' ')
2320	tkind = values.split(':')[5]
2321	timefmt = values.split(':')[6]
2322	title = values.split(':')[7].replace('\|',' ')
2323	gkind = values.split(':')[8]
2324
2325	leglabels = labels.split('@')[0].split(',')
2326	locleg = int(labels.split('@')[1])
2327
2328	Nsims = len(sims)
2329
2330	if Tint != '-1':
2331	tini = np.float(Tint.split('@')[0])
2332	tend = np.float(Tint.split('@')[1])
2333	else:
2334	tini = -1.
2335	tend = -1.
2336
2337	vartimetrjv = []
2338
2339	print ' ' + fname
2340	for trjfile in sims:
2341	print ' ' + trjfile + ' ...'
2342	if not os.path.isfile(trjfile):
2343	print errormsg
2344	print ' ' + fname + ": trajectory file: '" + trjfile + \
2345	"' does not exist !!"
2346	quit(-1)
2347
2348	trjobj = NetCDFFile(trjfile, 'r')
2349	otim = trjobj.variables['time']
2350	if not trjobj.variables.has_key(statisticskind + '_' + variable):
2351	print errormsg
2352	print ' ' + fname + ": file '" + trjfile + "' does not have variable '"+\
2353	statisticskind + '_' + variable + "' !!"
2354	quit(-1)
2355	ovar = trjobj.variables[statisticskind + '_' + variable]
2356	dimt = otim.shape[0]
2357
2358	if trjfile == sims[0]:
2359	gunits = ovar.getncattr('units')
2360	lname = ovar.getncattr('long_name')
2361	tunits = otim.getncattr('units')
2362
2363	if tini != -1:
2364	tiniid = -1
2365	tendid = -1
2366	for itv in range(dimt):
2367	if otim[itv] <= tini and otim[itv+1] >= tini: tiniid = itv
2368	if otim[itv] <= tend and otim[itv+1] >= tend: tendid = itv
2369
2370	if tiniid == -1 or tendid == -1:
2371	print errormsg
2372	print ' ' + main + ' time interval ', tini,',',tend,' not found: ', \
2373	tendid, ',', tiniid, ' !!'
2374	print ' data interval [',otim[0], otim[dimt-1],']'
2375	quit(-1)
2376	dimt = tendid - tiniid + 1
2377
2378	else:
2379	dimt = otim.shape[0]
2380
2381	valsv = np.zeros((2,dimt), dtype=np.float)
2382	# Checking for time consistency
2383	if otim.getncattr('units') != tunits:
2384	print warnmsg
2385	print ' ' + fname + ': different time units in the plot!!'
2386	newtimes = drw.coincident_CFtimes(otim[:], tunits, otim.getncattr('units'))
2387	else:
2388	newtimes = otim[:]
2389
2390	if tini == -1:
2391	valsv[1,:] = newtimes[:]
2392	valsv[0,:] = ovar[:]
2393	else:
2394	valsv[1,:] = newtimes[tiniid:tendid+1]
2395	valsv[0,:] = ovar[tiniid:tendid+1]
2396
2397	vartimetrjv.append(valsv)
2398	trjobj.close()
2399
2400	drw.plot_TimeSeries(vartimetrjv, drw.units_lunits(gunits), tunits, \
2401	'val_trajectories_' + statisticskind, gvarname, timetit, tkind, timefmt, title,\
2402	leglabels, locleg, gkind)
2403
2404	def variable_values(values):
2405	""" Function to give back values for a given variable
2406	values= [varname] name of the variable
2407	"""
2408
2409	fname = 'variable_values'
2410
2411	values = drw.variables_values(values)
2412
2413	print fname,'values:',values
2414	print fname,'variable_name:',values[0]
2415	print fname,'standard_name:',values[1]
2416	print fname,'min,max:',str(values[2]) + ',' + str(values[3])
2417	print fname,'long_name:',values[4]
2418	print fname,'units:',values[5]
2419	print fname,'shad_colors:',values[6]
2420	print fname,'all_values:',drw.numVector_String(values,',')
2421
2422	return
2423
2424	####### ###### ##### #### ### ## #
2425
2426	ngraphics = "'" + drw.numVector_String(namegraphics, "', '") + "'"
2427
2428	### Options
2429	##string_operation="operation to make: " + '\n' + " out, output values -S inidim1,[inidim2,...]:enddim1,[enddim2,...]"
2430	string_operation="""operation to make:
2431	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]
2432	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]
2433	[ckind]:
2434	'cmap': as it gets from colorbar
2435	'fixc,[colname]': fixed color [colname], all stright lines
2436	'fixsignc,[colname]': fixed color [colname], >0 stright, <0 dashed line
2437	"""
2438
2439	#print string_operation
2440
2441	parser = OptionParser()
2442	parser.add_option("-f", "--netCDF_file", dest="ncfile",
2443	help="file to use", metavar="FILE")
2444	parser.add_option("-o", "--operation", type='choice', dest="operation",
2445	choices=namegraphics,
2446	help="operation to make: " + ngraphics, metavar="OPER")
2447	parser.add_option("-S", "--valueS", dest="values",
2448	help="[WHEN APPLICABLE] values to use according to the operation", metavar="VALUES")
2449	parser.add_option("-v", "--variable", dest="varname",
2450	help="[WHEN APPLICABLE] variable to check", metavar="VAR")
2451
2452	(opts, args) = parser.parse_args()
2453
2454	####### #######
2455	## MAIN
2456	#######
2457
2458	# Not checking file operation
2459	Notcheckingfile = ['draw_2D_shad_cont', 'draw_2D_shad_cont_time', \
2460	'draw_2D_shad_line', 'draw_2D_shad_line_time', 'draw_lines', 'draw_lines_time', \
2461	'draw_topo_geogrid_boxes', 'draw_trajectories', 'draw_vals_trajectories', \
2462	'variable_values']
2463
2464	####### ###### ##### #### ### ## #
2465	errormsg='ERROR -- error -- ERROR -- error'
2466
2467	varn=opts.varname
2468	oper=opts.operation
2469
2470	if opts.ncfile is not None and not os.path.isfile(opts.ncfile) and \
2471	not drw.searchInlist(Notcheckingfile, oper):
2472	print errormsg
2473	print ' ' + main + ': File ' + opts.ncfile + ' does not exist !!'
2474	quit(-1)
2475
2476	if oper == 'create_movie':
2477	create_movie(opts.ncfile, opts.values, opts.varname)
2478	elif oper == 'draw_2D_shad':
2479	draw_2D_shad(opts.ncfile, opts.values, opts.varname)
2480	elif oper == 'draw_2D_shad_time':
2481	draw_2D_shad_time(opts.ncfile, opts.values, opts.varname)
2482	elif oper == 'draw_2D_shad_cont':
2483	draw_2D_shad_cont(opts.ncfile, opts.values, opts.varname)
2484	elif oper == 'draw_2D_shad_cont_time':
2485	draw_2D_shad_cont_time(opts.ncfile, opts.values, opts.varname)
2486	elif oper == 'draw_2D_shad_line':
2487	draw_2D_shad_line(opts.ncfile, opts.values, opts.varname)
2488	elif oper == 'draw_2D_shad_line_time':
2489	draw_2D_shad_line_time(opts.ncfile, opts.values, opts.varname)
2490	elif oper == 'draw_Neighbourghood_evol':
2491	draw_Neighbourghood_evol(opts.ncfile, opts.values, opts.varname)
2492	elif oper == 'draw_lines':
2493	draw_lines(opts.ncfile, opts.values, opts.varname)
2494	elif oper == 'draw_lines_time':
2495	draw_lines_time(opts.ncfile, opts.values, opts.varname)
2496	elif oper == 'draw_timeSeries':
2497	draw_timeSeries(opts.ncfile, opts.values, opts.varname)
2498	elif oper == 'draw_topo_geogrid':
2499	draw_topo_geogrid(opts.ncfile, opts.values)
2500	elif oper == 'draw_topo_geogrid_boxes':
2501	draw_topo_geogrid_boxes(opts.ncfile, opts.values)
2502	elif oper == 'draw_trajectories':
2503	draw_trajectories(opts.ncfile, opts.values, opts.varname)
2504	elif oper == 'draw_vals_trajectories':
2505	draw_vals_trajectories(opts.ncfile, opts.values, opts.varname)
2506	elif oper == 'list_graphics':
2507	# From: http://www.diveintopython.net/power_of_introspection/all_together.html
2508	import drawing as myself
2509	object = myself
2510	for opern in namegraphics:
2511	if opern != 'list_graphics':
2512	print opern + '_______ ______ _____ ____ ___ __ _'
2513	print getattr(object, opern).__doc__
2514	elif oper == 'variable_values':
2515	variable_values(opts.values)
2516	else:
2517	print errormsg
2518	print ' ' + main + ": the graphic '" + oper + "' is not ready !!"
2519	print errormsg
2520	quit()

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

Download in other formats: