Context Navigation

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

Last change on this file since 620 was 620, checked in by lfita, 9 years ago
Adding 'points' as paramter on 'draw_lines_time'
File size: 112.3 KB

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

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

Download in other formats: