Context Navigation

validation_sim.py @ 1444

Last change on this file since 1444 was 610, checked in by lfita, 10 years ago
Fixing wrong calculation of statistics Adding statistics of the values `between'
File size: 97.0 KB

Rev	Line
[330]	1
	2	# L. Fita, LMD-Jussieu. February 2015
[337]	3	## e.g. sfcEneAvigon # validation_sim.py -d X@west_east@None,Y@south_north@None,T@Time@time -D X@XLONG@longitude,Y@XLAT@latitude,T@time@time -k single-station -l 4.878773,43.915876,12. -o /home/lluis/DATA/obs/HyMeX/IOP15/sfcEnergyBalance_Avignon/OBSnetcdf.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v HFX@H,LH@LE,GRDFLX@G
[339]	4	## e.g. AIREP # validation_sim.py -d X@west_east@lon2D,Y@south_north@lat2D,Z@bottom_top@z2D,T@Time@time -D X@XLONG@longitude,Y@XLAT@latitude,Z@WRFz@alti,T@time@time -k trajectory -o /home/lluis/DATA/obs/HyMeX/IOP15/AIREP/2012/10/AIREP_121018.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v WRFt@t,WRFtd@td,WRFws@u,WRFwd@dd
[340]	5	## e.g. ATRCore # validation_sim.py -d X@west_east@lon2D,Y@south_north@lat2D,Z@bottom_top@z2D,T@Time@CFtime -D X@XLONG@longitude,Y@XLAT@latitude,Z@WRFz@altitude,T@time@time -k trajectory -o /home/lluis/DATA/obs/HyMeX/IOP15/ATRCore/V3/ATR_1Hz-HYMEXBDD-SOP1-v3_20121018_as120051.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v WRFt@air_temperature@subc@273.15,WRFp@air_pressure,WRFrh@relative_humidity,WRFrh@relative_humidity_Rosemount,WRFwd@wind_from_direction,WRFws@wind_speed
	6	## e.g. BAMED # validation_sim.py -d X@west_east@lon2D,Y@south_north@lat2D,Z@bottom_top@z2D,T@Time@CFtime -D X@XLONG@longitude,Y@XLAT@latitude,Z@WRFz@altitude,T@time@time -k trajectory -o /home/lluis/DATA/obs/HyMeX/IOP15/BAMED/BAMED_SOP1_B12_TOT5.nc -s /home/lluis/PY/wrfout_d01_2012-10-18_00:00:00.tests -v WRFt@tas_north,WRFp@pressure,WRFrh@hus,U@uas,V@vas
[337]	7
[330]	8	import numpy as np
	9	import os
	10	import re
	11	from optparse import OptionParser
	12	from netCDF4 import Dataset as NetCDFFile
[343]	13	from scipy import stats as sts
	14	import numpy.ma as ma
[330]	15
[528]	16	main = 'validation_sim.py'
[330]	17	errormsg = 'ERROR -- errror -- ERROR -- error'
	18	warnmsg = 'WARNING -- warning -- WARNING -- warning'
	19
	20	# version
	21	version=1.0
	22
	23	# Filling values for floats, integer and string
	24	fillValueF = 1.e20
	25	fillValueI = -99999
	26	fillValueS = '---'
	27
[333]	28	StringLength = 50
	29
[337]	30	# Number of grid points to take as 'environment' around the observed point
	31	Ngrid = 1
	32
[333]	33	def searchInlist(listname, nameFind):
	34	""" Function to search a value within a list
	35	listname = list
	36	nameFind = value to find
	37	>>> searInlist(['1', '2', '3', '5'], '5')
	38	True
	39	"""
	40	for x in listname:
	41	if x == nameFind:
	42	return True
	43	return False
	44
	45	def set_attribute(ncvar, attrname, attrvalue):
	46	""" Sets a value of an attribute of a netCDF variable. Removes previous attribute value if exists
	47	ncvar = object netcdf variable
	48	attrname = name of the attribute
	49	attrvalue = value of the attribute
	50	"""
	51	import numpy as np
	52	from netCDF4 import Dataset as NetCDFFile
	53
	54	attvar = ncvar.ncattrs()
	55	if searchInlist(attvar, attrname):
	56	attr = ncvar.delncattr(attrname)
	57
	58	attr = ncvar.setncattr(attrname, attrvalue)
	59
	60	return ncvar
	61
	62	def basicvardef(varobj, vstname, vlname, vunits):
	63	""" Function to give the basic attributes to a variable
	64	varobj= netCDF variable object
	65	vstname= standard name of the variable
	66	vlname= long name of the variable
	67	vunits= units of the variable
	68	"""
	69	attr = varobj.setncattr('standard_name', vstname)
	70	attr = varobj.setncattr('long_name', vlname)
	71	attr = varobj.setncattr('units', vunits)
	72
	73	return
	74
	75	def writing_str_nc(varo, values, Lchar):
	76	""" Function to write string values in a netCDF variable as a chain of 1char values
	77	varo= netCDF variable object
	78	values = list of values to introduce
	79	Lchar = length of the string in the netCDF file
	80	"""
	81
	82	Nvals = len(values)
	83	for iv in range(Nvals):
	84	stringv=values[iv]
	85	charvals = np.chararray(Lchar)
	86	Lstr = len(stringv)
	87	charvals[Lstr:Lchar] = ''
	88
	89	for ich in range(Lstr):
	90	charvals[ich] = stringv[ich:ich+1]
	91
	92	varo[iv,:] = charvals
	93
	94	return
	95
[601]	96	def CFtimes_datetime_NOfile(times, units, calendar):
	97	""" Provide date/time array from velues of netCDF CF-compilant time variable
	98	times= time values
	99	units= CF units of the variable time
	100	output:
	101	array(dimt, 0) = year
	102	array(dimt, 1) = month
	103	array(dimt, 2) = day
	104	array(dimt, 3) = hour
	105	array(dimt, 4) = minute
	106	array(dimt, 5) = second
	107	"""
	108	import datetime as dt
	109	fname = 'CFtimes_datetime_NOfile'
	110
	111	txtunits = units.split(' ')
	112	tunits = txtunits[0]
	113	Srefdate = txtunits[len(txtunits) - 1]
	114	# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
	115	##
	116	timeval = Srefdate.find(':')
	117
	118	if not timeval == -1:
	119	# print ' refdate with time!'
	120	refdate = datetimeStr_datetime(txtunits[len(txtunits) - 2] + '_' + Srefdate)
	121	else:
	122	refdate = dateStr_date(Srefdate)
	123
	124	dimt = len(times)
	125	realdates = np.zeros((dimt, 6), dtype=int)
	126
	127	secsDay=3600*24.
	128
	129	# Checking calendar!
	130	##
	131	y360 = False
	132	daycal360 = ['earth_360d', '360d', '360days', '360_day']
	133	if searchInlist(daycal360,calendar):
	134	print warnmsg
	135	print ' ' + fname + ': calendar of 12 months of 30 days !!'
	136	y360 = True
	137
	138	## Not in timedelta
	139	# if tunits == 'years':
	140	# for it in range(dimt):
	141	# realdate = refdate + dt.timedelta(years=float(times[it]))
	142	# realdates[it] = int(realdate.year)
	143	# elif tunits == 'months':
	144	# for it in range(dimt):
	145	# realdate = refdate + dt.timedelta(months=float(times[it]))
	146	# realdates[it] = int(realdate.year)
	147	if y360:
	148	if tunits == 'weeks':
	149	for it in range(dimt):
	150	deltat = dt.timedelta(weeks=float(times[it]))
	151	Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
	152	realdates[it,:] = date_juliandate(refdate.year,Tsecs)
	153	elif tunits == 'days':
	154	for it in range(dimt):
	155	deltat = dt.timedelta(days=float(times[it]))
	156	Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
	157	realdates[it,:] = date_juliandate(refdate.year,Tsecs)
	158	elif tunits == 'hours':
	159	for it in range(dimt):
	160	realdate = dt.timedelta(hours=float(times[it]))
	161	Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
	162	realdates[it,:] = date_juliandate(refdate.year,Tsecs)
	163	elif tunits == 'minutes':
	164	for it in range(dimt):
	165	realdate = dt.timedelta(minutes=float(times[it]))
	166	Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
	167	realdates[it,:] = date_juliandate(refdate.year,Tsecs)
	168	elif tunits == 'seconds':
	169	for it in range(dimt):
	170	realdate = dt.timedelta(seconds=float(times[it]))
	171	Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
	172	realdates[it,:] = date_juliandate(refdate.year,Tsecs)
	173	elif tunits == 'miliseconds':
	174	for it in range(dimt):
	175	realdate = dt.timedelta(miliseconds=float(times[it]))
	176	Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
	177	realdates[it,:] = date_juliandate(refdate.year,Tsecs)
	178	else:
	179	print errormsg
	180	print ' CFtimes_datetime: time units "' + tunits + '" not ready!!!!'
	181	quit(-1)
	182	else:
	183	if tunits == 'weeks':
	184	for it in range(dimt):
	185	realdate = refdate + dt.timedelta(weeks=float(times[it]))
	186	realdates[it,0] = int(realdate.year)
	187	realdates[it,1] = int(realdate.month)
	188	realdates[it,2] = int(realdate.day)
	189	realdates[it,3] = int(realdate.hour)
	190	realdates[it,4] = int(realdate.second)
	191	realdates[it,5] = int(realdate.minute)
	192	elif tunits == 'days':
	193	for it in range(dimt):
	194	realdate = refdate + dt.timedelta(days=float(times[it]))
	195	realdates[it,0] = int(realdate.year)
	196	realdates[it,1] = int(realdate.month)
	197	realdates[it,2] = int(realdate.day)
	198	realdates[it,3] = int(realdate.hour)
	199	realdates[it,4] = int(realdate.second)
	200	realdates[it,5] = int(realdate.minute)
	201	elif tunits == 'hours':
	202	for it in range(dimt):
	203	realdate = refdate + dt.timedelta(hours=float(times[it]))
	204	realdates[it,0] = int(realdate.year)
	205	realdates[it,1] = int(realdate.month)
	206	realdates[it,2] = int(realdate.day)
	207	realdates[it,3] = int(realdate.hour)
	208	realdates[it,4] = int(realdate.second)
	209	realdates[it,5] = int(realdate.minute)
	210	elif tunits == 'minutes':
	211	for it in range(dimt):
	212	realdate = refdate + dt.timedelta(minutes=float(times[it]))
	213	realdates[it,0] = int(realdate.year)
	214	realdates[it,1] = int(realdate.month)
	215	realdates[it,2] = int(realdate.day)
	216	realdates[it,3] = int(realdate.hour)
	217	realdates[it,4] = int(realdate.second)
	218	realdates[it,5] = int(realdate.minute)
	219	elif tunits == 'seconds':
	220	for it in range(dimt):
	221	realdate = refdate + dt.timedelta(seconds=float(times[it]))
	222	realdates[it,0] = int(realdate.year)
	223	realdates[it,1] = int(realdate.month)
	224	realdates[it,2] = int(realdate.day)
	225	realdates[it,3] = int(realdate.hour)
	226	realdates[it,4] = int(realdate.second)
	227	realdates[it,5] = int(realdate.minute)
	228	elif tunits == 'milliseconds':
	229	for it in range(dimt):
	230	realdate = refdate + dt.timedelta(milliseconds=float(times[it]))
	231	realdates[it,0] = int(realdate.year)
	232	realdates[it,1] = int(realdate.month)
	233	realdates[it,2] = int(realdate.day)
	234	realdates[it,3] = int(realdate.hour)
	235	realdates[it,4] = int(realdate.second)
	236	realdates[it,5] = int(realdate.minute)
	237	else:
	238	print errormsg
	239	print ' CFtimes_datetime: time units "' + tunits + '" not ready!!!!'
	240	quit(-1)
	241
	242	return realdates
	243
[330]	244	def index_3mat(matA,matB,matC,val):
	245	""" Function to provide the coordinates of a given value inside three matrix simultaneously
	246	index_mat(matA,matB,matC,val)
	247	matA= matrix with one set of values
	248	matB= matrix with the other set of values
	249	matB= matrix with the third set of values
	250	val= triplet of values to search
	251	>>> index_mat(np.arange(27).reshape(3,3,3),np.arange(100,127).reshape(3,3,3),np.arange(200,227).reshape(3,3,3),[22,122,222])
	252	[2 1 1]
	253	"""
	254	fname = 'index_3mat'
	255
	256	matAshape = matA.shape
	257	matBshape = matB.shape
	258	matCshape = matC.shape
	259
	260	for idv in range(len(matAshape)):
	261	if matAshape[idv] != matBshape[idv]:
	262	print errormsg
	263	print ' ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv], \
	264	'and B:',matBshape[idv],'does not coincide!!'
	265	quit(-1)
	266	if matAshape[idv] != matCshape[idv]:
	267	print errormsg
	268	print ' ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv], \
	269	'and C:',matCshape[idv],'does not coincide!!'
	270	quit(-1)
	271
	272	minA = np.min(matA)
	273	maxA = np.max(matA)
	274	minB = np.min(matB)
	275	maxB = np.max(matB)
	276	minC = np.min(matC)
	277	maxC = np.max(matC)
	278
	279	if val[0] < minA or val[0] > maxA:
	280	print warnmsg
	281	print ' ' + fname + ': first value:',val[0],'outside matA range',minA,',', \
	282	maxA,'!!'
	283	if val[1] < minB or val[1] > maxB:
	284	print warnmsg
	285	print ' ' + fname + ': second value:',val[1],'outside matB range',minB,',', \
	286	maxB,'!!'
	287	if val[2] < minC or val[2] > maxC:
	288	print warnmsg
	289	print ' ' + fname + ': second value:',val[2],'outside matC range',minC,',', \
	290	maxC,'!!'
	291
	292	dist = np.zeros(tuple(matAshape), dtype=np.float)
	293	dist = np.sqrt((matA - np.float(val[0]))2 + (matB - np.float(val[1]))2 + \
	294	(matC - np.float(val[2]))**2)
	295
	296	mindist = np.min(dist)
	297
	298	matlist = list(dist.flatten())
	299	ifound = matlist.index(mindist)
	300
	301	Ndims = len(matAshape)
	302	valpos = np.zeros((Ndims), dtype=int)
	303	baseprevdims = np.zeros((Ndims), dtype=int)
	304
	305	for dimid in range(Ndims):
	306	baseprevdims[dimid] = np.product(matAshape[dimid+1:Ndims])
	307	if dimid == 0:
	308	alreadyplaced = 0
	309	else:
	310	alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
	311	valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
	312
	313	return valpos
	314
	315	def index_2mat(matA,matB,val):
	316	""" Function to provide the coordinates of a given value inside two matrix simultaneously
	317	index_mat(matA,matB,val)
	318	matA= matrix with one set of values
	319	matB= matrix with the pother set of values
	320	val= couple of values to search
[337]	321	>>> index_2mat(np.arange(27).reshape(3,3,3),np.arange(100,127).reshape(3,3,3),[22,111])
[330]	322	[2 1 1]
	323	"""
	324	fname = 'index_2mat'
	325
	326	matAshape = matA.shape
	327	matBshape = matB.shape
	328
	329	for idv in range(len(matAshape)):
	330	if matAshape[idv] != matBshape[idv]:
	331	print errormsg
	332	print ' ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv], \
	333	'and B:',matBshape[idv],'does not coincide!!'
	334	quit(-1)
	335
	336	minA = np.min(matA)
	337	maxA = np.max(matA)
	338	minB = np.min(matB)
	339	maxB = np.max(matB)
	340
[337]	341	Ndims = len(matAshape)
	342	# valpos = np.ones((Ndims), dtype=int)*-1.
	343	valpos = np.zeros((Ndims), dtype=int)
	344
[330]	345	if val[0] < minA or val[0] > maxA:
	346	print warnmsg
	347	print ' ' + fname + ': first value:',val[0],'outside matA range',minA,',', \
	348	maxA,'!!'
[337]	349	return valpos
[330]	350	if val[1] < minB or val[1] > maxB:
	351	print warnmsg
	352	print ' ' + fname + ': second value:',val[1],'outside matB range',minB,',', \
	353	maxB,'!!'
[337]	354	return valpos
[330]	355
	356	dist = np.zeros(tuple(matAshape), dtype=np.float)
	357	dist = np.sqrt((matA - np.float(val[0]))2 + (matB - np.float(val[1]))2)
	358
	359	mindist = np.min(dist)
	360
[337]	361	if mindist != mindist:
	362	print ' ' + fname + ': wrong minimal distance',mindist,'!!'
	363	return valpos
	364	else:
	365	matlist = list(dist.flatten())
	366	ifound = matlist.index(mindist)
[330]	367
	368	baseprevdims = np.zeros((Ndims), dtype=int)
	369	for dimid in range(Ndims):
	370	baseprevdims[dimid] = np.product(matAshape[dimid+1:Ndims])
	371	if dimid == 0:
	372	alreadyplaced = 0
	373	else:
	374	alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
	375	valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
	376
	377	return valpos
	378
[337]	379	def index_mat(matA,val):
[330]	380	""" Function to provide the coordinates of a given value inside a matrix
[337]	381	index_mat(matA,val)
	382	matA= matrix with one set of values
	383	val= couple of values to search
	384	>>> index_mat(np.arange(27),22.3)
	385	22
	386	"""
	387	fname = 'index_mat'
	388
	389	matAshape = matA.shape
	390
	391	minA = np.min(matA)
	392	maxA = np.max(matA)
	393
	394	Ndims = len(matAshape)
	395	# valpos = np.ones((Ndims), dtype=int)*-1.
	396	valpos = np.zeros((Ndims), dtype=int)
	397
	398	if val < minA or val > maxA:
	399	print warnmsg
	400	print ' ' + fname + ': first value:',val,'outside matA range',minA,',', \
	401	maxA,'!!'
	402	return valpos
	403
	404	dist = np.zeros(tuple(matAshape), dtype=np.float)
	405	dist = (matA - np.float(val))**2
	406
	407	mindist = np.min(dist)
	408	if mindist != mindist:
	409	print ' ' + fname + ': wrong minimal distance',mindist,'!!'
	410	return valpos
	411
	412	matlist = list(dist.flatten())
	413	valpos = matlist.index(mindist)
	414
	415	return valpos
	416
	417	def index_mat_exact(mat,val):
	418	""" Function to provide the coordinates of a given exact value inside a matrix
[330]	419	index_mat(mat,val)
	420	mat= matrix with values
	421	val= value to search
	422	>>> index_mat(np.arange(27).reshape(3,3,3),22)
	423	[2 1 1]
	424	"""
	425
	426	fname = 'index_mat'
	427
	428	matshape = mat.shape
	429
	430	matlist = list(mat.flatten())
	431	ifound = matlist.index(val)
	432
	433	Ndims = len(matshape)
	434	valpos = np.zeros((Ndims), dtype=int)
	435	baseprevdims = np.zeros((Ndims), dtype=int)
	436
	437	for dimid in range(Ndims):
	438	baseprevdims[dimid] = np.product(matshape[dimid+1:Ndims])
	439	if dimid == 0:
	440	alreadyplaced = 0
	441	else:
	442	alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
	443	valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
	444
	445	return valpos
	446
[343]	447	def datetimeStr_datetime(StringDT):
	448	""" Function to transform a string date ([YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format) to a date object
	449	>>> datetimeStr_datetime('1976-02-17_00:00:00')
	450	1976-02-17 00:00:00
	451	"""
	452	import datetime as dt
	453
	454	fname = 'datetimeStr_datetime'
	455
	456	dateD = np.zeros((3), dtype=int)
	457	timeT = np.zeros((3), dtype=int)
	458
	459	dateD[0] = int(StringDT[0:4])
	460	dateD[1] = int(StringDT[5:7])
	461	dateD[2] = int(StringDT[8:10])
	462
	463	trefT = StringDT.find(':')
	464	if not trefT == -1:
	465	# print ' ' + fname + ': refdate with time!'
	466	timeT[0] = int(StringDT[11:13])
	467	timeT[1] = int(StringDT[14:16])
	468	timeT[2] = int(StringDT[17:19])
	469
	470	if int(dateD[0]) == 0:
	471	print warnmsg
	472	print ' ' + fname + ': 0 reference year!! changing to 1'
	473	dateD[0] = 1
	474
	475	newdatetime = dt.datetime(dateD[0], dateD[1], dateD[2], timeT[0], timeT[1], timeT[2])
	476
	477	return newdatetime
	478
	479	def datetimeStr_conversion(StringDT,typeSi,typeSo):
	480	""" Function to transform a string date to an another date object
	481	StringDT= string with the date and time
	482	typeSi= type of datetime string input
	483	typeSo= type of datetime string output
	484	[typeSi/o]
	485	'cfTime': [time],[units]; ]time in CF-convention format [units] = [tunits] since [refdate]
	486	'matYmdHMS': numerical vector with [[YYYY], [MM], [DD], [HH], [MI], [SS]]
	487	'YmdHMS': [YYYY][MM][DD][HH][MI][SS] format
	488	'Y-m-d_H:M:S': [YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format
	489	'Y-m-d H:M:S': [YYYY]-[MM]-[DD] [HH]:[MI]:[SS] format
	490	'Y/m/d H-M-S': [YYYY]/[MM]/[DD] [HH]-[MI]-[SS] format
	491	'WRFdatetime': [Y], [Y], [Y], [Y], '-', [M], [M], '-', [D], [D], '_', [H],
	492	[H], ':', [M], [M], ':', [S], [S]
	493	>>> datetimeStr_conversion('1976-02-17_08:32:05','Y-m-d_H:M:S','matYmdHMS')
	494	[1976 2 17 8 32 5]
	495	>>> datetimeStr_conversion(str(137880)+',minutes since 1979-12-01_00:00:00','cfTime','Y/m/d H-M-S')
	496	1980/03/05 18-00-00
	497	"""
	498	import datetime as dt
	499
	500	fname = 'datetimeStr_conversion'
	501
	502	if StringDT[0:1] == 'h':
	503	print fname + '_____________________________________________________________'
	504	print datetimeStr_conversion.__doc__
	505	quit()
	506
	507	if typeSi == 'cfTime':
	508	timeval = np.float(StringDT.split(',')[0])
	509	tunits = StringDT.split(',')[1].split(' ')[0]
	510	Srefdate = StringDT.split(',')[1].split(' ')[2]
	511
	512	# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
	513	##
	514	yrref=Srefdate[0:4]
	515	monref=Srefdate[5:7]
	516	dayref=Srefdate[8:10]
	517
	518	trefT = Srefdate.find(':')
	519	if not trefT == -1:
	520	# print ' ' + fname + ': refdate with time!'
	521	horref=Srefdate[11:13]
	522	minref=Srefdate[14:16]
	523	secref=Srefdate[17:19]
	524	refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref + \
	525	'_' + horref + ':' + minref + ':' + secref)
	526	else:
	527	refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref + \
	528	+ '_00:00:00')
	529
	530	if tunits == 'weeks':
	531	newdate = refdate + dt.timedelta(weeks=float(timeval))
	532	elif tunits == 'days':
	533	newdate = refdate + dt.timedelta(days=float(timeval))
	534	elif tunits == 'hours':
	535	newdate = refdate + dt.timedelta(hours=float(timeval))
	536	elif tunits == 'minutes':
	537	newdate = refdate + dt.timedelta(minutes=float(timeval))
	538	elif tunits == 'seconds':
	539	newdate = refdate + dt.timedelta(seconds=float(timeval))
	540	elif tunits == 'milliseconds':
	541	newdate = refdate + dt.timedelta(milliseconds=float(timeval))
	542	else:
	543	print errormsg
	544	print ' timeref_datetime: time units "' + tunits + '" not ready!!!!'
	545	quit(-1)
	546
	547	yr = newdate.year
	548	mo = newdate.month
	549	da = newdate.day
	550	ho = newdate.hour
	551	mi = newdate.minute
	552	se = newdate.second
	553	elif typeSi == 'matYmdHMS':
	554	yr = StringDT[0]
	555	mo = StringDT[1]
	556	da = StringDT[2]
	557	ho = StringDT[3]
	558	mi = StringDT[4]
	559	se = StringDT[5]
	560	elif typeSi == 'YmdHMS':
	561	yr = int(StringDT[0:4])
	562	mo = int(StringDT[4:6])
	563	da = int(StringDT[6:8])
	564	ho = int(StringDT[8:10])
	565	mi = int(StringDT[10:12])
	566	se = int(StringDT[12:14])
	567	elif typeSi == 'Y-m-d_H:M:S':
	568	dateDT = StringDT.split('_')
	569	dateD = dateDT[0].split('-')
	570	timeT = dateDT[1].split(':')
	571	yr = int(dateD[0])
	572	mo = int(dateD[1])
	573	da = int(dateD[2])
	574	ho = int(timeT[0])
	575	mi = int(timeT[1])
	576	se = int(timeT[2])
	577	elif typeSi == 'Y-m-d H:M:S':
	578	dateDT = StringDT.split(' ')
	579	dateD = dateDT[0].split('-')
	580	timeT = dateDT[1].split(':')
	581	yr = int(dateD[0])
	582	mo = int(dateD[1])
	583	da = int(dateD[2])
	584	ho = int(timeT[0])
	585	mi = int(timeT[1])
	586	se = int(timeT[2])
	587	elif typeSi == 'Y/m/d H-M-S':
	588	dateDT = StringDT.split(' ')
	589	dateD = dateDT[0].split('/')
	590	timeT = dateDT[1].split('-')
	591	yr = int(dateD[0])
	592	mo = int(dateD[1])
	593	da = int(dateD[2])
	594	ho = int(timeT[0])
	595	mi = int(timeT[1])
	596	se = int(timeT[2])
	597	elif typeSi == 'WRFdatetime':
	598	yr = int(StringDT[0])1000 + int(StringDT[1])100 + int(StringDT[2])*10 + \
	599	int(StringDT[3])
	600	mo = int(StringDT[5])*10 + int(StringDT[6])
	601	da = int(StringDT[8])*10 + int(StringDT[9])
	602	ho = int(StringDT[11])*10 + int(StringDT[12])
	603	mi = int(StringDT[14])*10 + int(StringDT[15])
	604	se = int(StringDT[17])*10 + int(StringDT[18])
	605	else:
	606	print errormsg
	607	print ' ' + fname + ': type of String input date "' + typeSi + \
	608	'" not ready !!!!'
	609	quit(-1)
	610
	611	if typeSo == 'matYmdHMS':
	612	dateYmdHMS = np.zeros((6), dtype=int)
	613	dateYmdHMS[0] = yr
	614	dateYmdHMS[1] = mo
	615	dateYmdHMS[2] = da
	616	dateYmdHMS[3] = ho
	617	dateYmdHMS[4] = mi
	618	dateYmdHMS[5] = se
	619	elif typeSo == 'YmdHMS':
	620	dateYmdHMS = str(yr).zfill(4) + str(mo).zfill(2) + str(da).zfill(2) + \
	621	str(ho).zfill(2) + str(mi).zfill(2) + str(se).zfill(2)
	622	elif typeSo == 'Y-m-d_H:M:S':
	623	dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' + \
	624	str(da).zfill(2) + '_' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \
	625	str(se).zfill(2)
	626	elif typeSo == 'Y-m-d H:M:S':
	627	dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' + \
	628	str(da).zfill(2) + ' ' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \
	629	str(se).zfill(2)
	630	elif typeSo == 'Y/m/d H-M-S':
	631	dateYmdHMS = str(yr).zfill(4) + '/' + str(mo).zfill(2) + '/' + \
	632	str(da).zfill(2) + ' ' + str(ho).zfill(2) + '-' + str(mi).zfill(2) + '-' + \
	633	str(se).zfill(2)
	634	elif typeSo == 'WRFdatetime':
	635	dateYmdHMS = []
	636	yM = yr/1000
	637	yC = (yr-yM*1000)/100
	638	yD = (yr-yM1000-yC100)/10
	639	yU = yr-yM1000-yC100-yD*10
	640
	641	mD = mo/10
	642	mU = mo-mD*10
	643
	644	dD = da/10
	645	dU = da-dD*10
	646
	647	hD = ho/10
	648	hU = ho-hD*10
	649
	650	miD = mi/10
	651	miU = mi-miD*10
	652
	653	sD = se/10
	654	sU = se-sD*10
	655
	656	dateYmdHMS.append(str(yM))
	657	dateYmdHMS.append(str(yC))
	658	dateYmdHMS.append(str(yD))
	659	dateYmdHMS.append(str(yU))
	660	dateYmdHMS.append('-')
	661	dateYmdHMS.append(str(mD))
	662	dateYmdHMS.append(str(mU))
	663	dateYmdHMS.append('-')
	664	dateYmdHMS.append(str(dD))
	665	dateYmdHMS.append(str(dU))
	666	dateYmdHMS.append('_')
	667	dateYmdHMS.append(str(hD))
	668	dateYmdHMS.append(str(hU))
	669	dateYmdHMS.append(':')
	670	dateYmdHMS.append(str(miD))
	671	dateYmdHMS.append(str(miU))
	672	dateYmdHMS.append(':')
	673	dateYmdHMS.append(str(sD))
	674	dateYmdHMS.append(str(sU))
	675	else:
	676	print errormsg
	677	print ' ' + fname + ': type of output date "' + typeSo + '" not ready !!!!'
	678	quit(-1)
	679
	680	return dateYmdHMS
	681
[330]	682	def coincident_CFtimes(tvalB, tunitA, tunitB):
	683	""" Function to make coincident times for two different sets of CFtimes
	684	tvalB= time values B
	685	tunitA= time units times A to which we want to make coincidence
	686	tunitB= time units times B
	687	>>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00',
	688	'hours since 1949-12-01 00:00:00')
	689	[ 0. 3600. 7200. 10800. 14400. 18000. 21600. 25200. 28800. 32400.]
	690	>>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00',
	691	'hours since 1979-12-01 00:00:00')
	692	[ 9.46684800e+08 9.46688400e+08 9.46692000e+08 9.46695600e+08
	693	9.46699200e+08 9.46702800e+08 9.46706400e+08 9.46710000e+08
	694	9.46713600e+08 9.46717200e+08]
	695	"""
	696	import datetime as dt
	697	fname = 'coincident_CFtimes'
	698
	699	trefA = tunitA.split(' ')[2] + ' ' + tunitA.split(' ')[3]
	700	trefB = tunitB.split(' ')[2] + ' ' + tunitB.split(' ')[3]
	701	tuA = tunitA.split(' ')[0]
	702	tuB = tunitB.split(' ')[0]
	703
	704	if tuA != tuB:
	705	if tuA == 'microseconds':
	706	if tuB == 'microseconds':
	707	tB = tvalB*1.
	708	elif tuB == 'seconds':
	709	tB = tvalB*10.e6
	710	elif tuB == 'minutes':
	711	tB = tvalB60.10.e6
	712	elif tuB == 'hours':
	713	tB = tvalB3600.10.e6
	714	elif tuB == 'days':
	715	tB = tvalB3600.24.*10.e6
	716	else:
	717	print errormsg
	718	print ' ' + fname + ": combination of time untis: '" + tuA + \
	719	"' & '" + tuB + "' not ready !!"
	720	quit(-1)
	721	elif tuA == 'seconds':
	722	if tuB == 'microseconds':
	723	tB = tvalB/10.e6
	724	elif tuB == 'seconds':
	725	tB = tvalB*1.
	726	elif tuB == 'minutes':
	727	tB = tvalB*60.
	728	elif tuB == 'hours':
	729	tB = tvalB*3600.
	730	elif tuB == 'days':
	731	tB = tvalB3600.24.
	732	else:
	733	print errormsg
	734	print ' ' + fname + ": combination of time untis: '" + tuA + \
	735	"' & '" + tuB + "' not ready !!"
	736	quit(-1)
	737	elif tuA == 'minutes':
	738	if tuB == 'microseconds':
	739	tB = tvalB/(60.*10.e6)
	740	elif tuB == 'seconds':
	741	tB = tvalB/60.
	742	elif tuB == 'minutes':
	743	tB = tvalB*1.
	744	elif tuB == 'hours':
	745	tB = tvalB*60.
	746	elif tuB == 'days':
	747	tB = tvalB60.24.
	748	else:
	749	print errormsg
	750	print ' ' + fname + ": combination of time untis: '" + tuA + \
	751	"' & '" + tuB + "' not ready !!"
	752	quit(-1)
	753	elif tuA == 'hours':
	754	if tuB == 'microseconds':
	755	tB = tvalB/(3600.*10.e6)
	756	elif tuB == 'seconds':
	757	tB = tvalB/3600.
	758	elif tuB == 'minutes':
	759	tB = tvalB/60.
	760	elif tuB == 'hours':
	761	tB = tvalB*1.
	762	elif tuB == 'days':
	763	tB = tvalB*24.
	764	else:
	765	print errormsg
	766	print ' ' + fname + ": combination of time untis: '" + tuA + \
	767	"' & '" + tuB + "' not ready !!"
	768	quit(-1)
	769	elif tuA == 'days':
	770	if tuB == 'microseconds':
	771	tB = tvalB/(24.3600.10.e6)
	772	elif tuB == 'seconds':
	773	tB = tvalB/(24.*3600.)
	774	elif tuB == 'minutes':
	775	tB = tvalB/(24.*60.)
	776	elif tuB == 'hours':
	777	tB = tvalB/24.
	778	elif tuB == 'days':
	779	tB = tvalB*1.
	780	else:
	781	print errormsg
	782	print ' ' + fname + ": combination of time untis: '" + tuA + \
	783	"' & '" + tuB + "' not ready !!"
	784	quit(-1)
	785	else:
	786	print errormsg
	787	print ' ' + fname + ": time untis: '" + tuA + "' not ready !!"
	788	quit(-1)
	789	else:
	790	tB = tvalB*1.
	791
	792	if trefA != trefB:
	793	trefTA = dt.datetime.strptime(trefA, '%Y-%m-%d %H:%M:%S')
	794	trefTB = dt.datetime.strptime(trefB, '%Y-%m-%d %H:%M:%S')
	795
	796	difft = trefTB - trefTA
[602]	797	diffv = difft.days24.3600. + difft.seconds + difft.microseconds/1.e6
[330]	798	print ' ' + fname + ': different reference refA:',trefTA,'refB',trefTB
[602]	799	print ' difference:',difft,':',diffv,'(in seconds)'
[330]	800
	801	if tuA == 'microseconds':
[602]	802	tB = tB + diffv*1.e6
	803	elif tuA == 'seconds':
[330]	804	tB = tB + diffv
	805	elif tuA == 'minutes':
[602]	806	tB = tB + diffv/(60.)
[330]	807	elif tuA == 'hours':
[602]	808	tB = tB + diffv/(3600.)
	809	elif tuA == 'days':
	810	tB = tB + diffv/(24.*3600.)
[330]	811	else:
	812	print errormsg
	813	print ' ' + fname + ": time untis: '" + tuA + "' not ready !!"
	814	quit(-1)
	815
	816	return tB
	817
[333]	818	def slice_variable(varobj, dimslice):
	819	""" Function to return a slice of a given variable according to values to its
	820	dimensions
	821	slice_variable(varobj, dimslice)
	822	varobj= object wit the variable
	823	dimslice= [[dimname1]:[value1]\|[[dimname2]:[value2], ...] pairs of dimension
	824	[value]:
	825	* [integer]: which value of the dimension
	826	* -1: all along the dimension
	827	* -9: last value of the dimension
	828	* [beg]@[end] slice from [beg] to [end]
	829	"""
	830	fname = 'slice_variable'
	831
	832	if varobj == 'h':
	833	print fname + '_____________________________________________________________'
	834	print slice_variable.__doc__
	835	quit()
	836
	837	vardims = varobj.dimensions
	838	Ndimvar = len(vardims)
	839
	840	Ndimcut = len(dimslice.split('\|'))
	841	dimsl = dimslice.split('\|')
	842
	843	varvalsdim = []
	844	dimnslice = []
	845
	846	for idd in range(Ndimvar):
	847	for idc in range(Ndimcut):
	848	dimcutn = dimsl[idc].split(':')[0]
	849	dimcutv = dimsl[idc].split(':')[1]
	850	if vardims[idd] == dimcutn:
	851	posfrac = dimcutv.find('@')
	852	if posfrac != -1:
	853	inifrac = int(dimcutv.split('@')[0])
	854	endfrac = int(dimcutv.split('@')[1])
	855	varvalsdim.append(slice(inifrac,endfrac))
	856	dimnslice.append(vardims[idd])
	857	else:
	858	if int(dimcutv) == -1:
	859	varvalsdim.append(slice(0,varobj.shape[idd]))
	860	dimnslice.append(vardims[idd])
	861	elif int(dimcutv) == -9:
	862	varvalsdim.append(int(varobj.shape[idd])-1)
	863	else:
	864	varvalsdim.append(int(dimcutv))
	865	break
	866
	867	varvalues = varobj[tuple(varvalsdim)]
	868
	869	return varvalues, dimnslice
	870
[337]	871	def func_compute_varNOcheck(ncobj, varn):
	872	""" Function to compute variables which are not originary in the file
	873	ncobj= netCDF object file
	874	varn = variable to compute:
	875	'WRFdens': air density from WRF variables
	876	'WRFght': geopotential height from WRF variables
	877	'WRFp': pressure from WRF variables
	878	'WRFrh': relative humidty fom WRF variables
	879	'WRFt': temperature from WRF variables
[505]	880	'WRFwds': surface wind direction from WRF variables
	881	'WRFwss': surface wind speed from WRF variables
[337]	882	'WRFz': height from WRF variables
	883	"""
	884	fname = 'compute_varNOcheck'
[333]	885
[337]	886	if varn == 'WRFdens':
	887	# print ' ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \
	888	# 'DNW)/g ...'
	889	grav = 9.81
	890
	891	# Just we need in in absolute values: Size of the central grid cell
	892	## dxval = ncobj.getncattr('DX')
	893	## dyval = ncobj.getncattr('DY')
	894	## mapfac = ncobj.variables['MAPFAC_M'][:]
	895	## area = dxvaldyvalmapfac
	896	dimensions = ncobj.variables['MU'].dimensions
	897
	898	mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
	899	dnw = ncobj.variables['DNW'][:]
	900
	901	varNOcheckv = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \
	902	dtype=np.float)
	903	levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
	904
	905	for it in range(mu.shape[0]):
	906	for iz in range(dnw.shape[1]):
	907	levval.fill(np.abs(dnw[it,iz]))
	908	varNOcheck[it,iz,:,:] = levval
	909	varNOcheck[it,iz,:,:] = mu[it,:,:]*varNOcheck[it,iz,:,:]/grav
	910
	911	elif varn == 'WRFght':
	912	# print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
	913	varNOcheckv = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
	914	dimensions = ncobj.variables['PH'].dimensions
	915
	916	elif varn == 'WRFp':
	917	# print ' ' + fname + ': Retrieving pressure value from WRF as P + PB'
	918	varNOcheckv = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	919	dimensions = ncobj.variables['P'].dimensions
	920
	921	elif varn == 'WRFrh':
	922	# print ' ' + main + ": computing relative humidity from WRF as 'Tetens'" +\
	923	# ' equation (T,P) ...'
	924	p0=100000.
	925	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	926	tk = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
	927	qv = ncobj.variables['QVAPOR'][:]
	928
	929	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
	930	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
	931
	932	varNOcheckv = qv/data2
	933	dimensions = ncobj.variables['P'].dimensions
	934
	935	elif varn == 'WRFt':
	936	# print ' ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
	937	p0=100000.
	938	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	939
	940	varNOcheckv = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
	941	dimensions = ncobj.variables['T'].dimensions
	942
[505]	943	elif varn == 'WRFwds':
	944	# print ' ' + main + ': computing surface wind direction from WRF as ATAN2(V,U) ...'
	945	varNOcheckv = np.arctan2(ncobj.variables['V10'][:], ncobj.variables['U10'][:])
	946	dimensions = ncobj.variables['V10'].dimensions
	947
	948	elif varn == 'WRFwss':
	949	# print ' ' + main + ': computing surface wind speed from WRF as SQRT(U2 + V2) ...'
	950	varNOcheckv = np.sqrt(ncobj.variables['U10'][:]*ncobj.variables['U10'][:] + \
	951	ncobj.variables['V10'][:]*ncobj.variables['V10'][:])
	952	dimensions = ncobj.variables['U10'].dimensions
	953
[337]	954	elif varn == 'WRFz':
	955	grav = 9.81
	956	# print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
	957	varNOcheckv = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/grav
	958	dimensions = ncobj.variables['PH'].dimensions
	959
	960	else:
	961	print erromsg
	962	print ' ' + fname + ": variable '" + varn + "' nor ready !!"
	963	quit(-1)
	964
	965	return varNOcheck
	966
	967	class compute_varNOcheck(object):
	968	""" Class to compute variables which are not originary in the file
	969	ncobj= netCDF object file
	970	varn = variable to compute:
	971	'WRFdens': air density from WRF variables
	972	'WRFght': geopotential height from WRF variables
	973	'WRFp': pressure from WRF variables
	974	'WRFrh': relative humidty fom WRF variables
[505]	975	'TSrhs': surface relative humidty fom TS variables
	976	'WRFrhs': surface relative humidty fom WRF variables
[343]	977	'WRFT': CF-time from WRF variables
[337]	978	'WRFt': temperature from WRF variables
[564]	979	'TStd': dew-point temperature from TS variables
[339]	980	'WRFtd': dew-point temperature from WRF variables
[505]	981	'WRFwd': wind direction from WRF variables
	982	'TSwds': surface wind direction from TS variables
	983	'WRFwds': surface wind direction from WRF variables
[339]	984	'WRFws': wind speed from WRF variables
[505]	985	'TSwss': surface wind speed from TS variables
	986	'WRFwss': surface wind speed from WRF variables
[337]	987	'WRFz': height from WRF variables
	988	"""
	989	fname = 'compute_varNOcheck'
	990
	991	def __init__(self, ncobj, varn):
	992
	993	if ncobj is None:
	994	self = None
	995	self.dimensions = None
	996	self.shape = None
	997	self.__values = None
	998	else:
	999	if varn == 'WRFdens':
	1000	# print ' ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \
	1001	# 'DNW)/g ...'
	1002	grav = 9.81
	1003
	1004	# Just we need in in absolute values: Size of the central grid cell
	1005	## dxval = ncobj.getncattr('DX')
	1006	## dyval = ncobj.getncattr('DY')
	1007	## mapfac = ncobj.variables['MAPFAC_M'][:]
	1008	## area = dxvaldyvalmapfac
	1009	dimensions = ncobj.variables['MU'].dimensions
	1010	shape = ncobj.variables['MU'].shape
	1011
	1012	mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
	1013	dnw = ncobj.variables['DNW'][:]
	1014
	1015	varNOcheckv = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \
	1016	dtype=np.float)
	1017	levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
	1018
	1019	for it in range(mu.shape[0]):
	1020	for iz in range(dnw.shape[1]):
	1021	levval.fill(np.abs(dnw[it,iz]))
	1022	varNOcheck[it,iz,:,:] = levval
	1023	varNOcheck[it,iz,:,:] = mu[it,:,:]*varNOcheck[it,iz,:,:]/grav
	1024
	1025	elif varn == 'WRFght':
	1026	# print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
	1027	varNOcheckv = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
	1028	dimensions = ncobj.variables['PH'].dimensions
	1029	shape = ncobj.variables['PH'].shape
	1030
	1031	elif varn == 'WRFp':
	1032	# print ' ' + fname + ': Retrieving pressure value from WRF as P + PB'
	1033	varNOcheckv = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	1034	dimensions = ncobj.variables['P'].dimensions
	1035	shape = ncobj.variables['P'].shape
	1036
	1037	elif varn == 'WRFrh':
	1038	# print ' ' + main + ": computing relative humidity from WRF as 'Tetens'" +\
	1039	# ' equation (T,P) ...'
	1040	p0=100000.
	1041	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
[505]	1042	tk = (ncobj.variables['T'][:])(p/p0)*(2./7.)
[337]	1043	qv = ncobj.variables['QVAPOR'][:]
	1044
	1045	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
	1046	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
	1047
	1048	varNOcheckv = qv/data2
	1049	dimensions = ncobj.variables['P'].dimensions
	1050	shape = ncobj.variables['P'].shape
	1051
[505]	1052	elif varn == 'TSrhs':
	1053	# print ' ' + main + ": computing surface relative humidity from TSs as 'Tetens'" +\
	1054	# ' equation (T,P) ...'
	1055	p0=100000.
	1056	p=ncobj.variables['psfc'][:]
	1057	tk = (ncobj.variables['t'][:])(p/p0)*(2./7.)
	1058	qv = ncobj.variables['q'][:]
	1059
	1060	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
	1061	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
	1062
	1063	varNOcheckv = qv/data2
	1064	dimensions = ncobj.variables['psfc'].dimensions
	1065	shape = ncobj.variables['psfc'].shape
	1066
	1067	elif varn == 'WRFrhs':
	1068	# print ' ' + main + ": computing surface relative humidity from WRF as 'Tetens'" +\
	1069	# ' equation (T,P) ...'
	1070	p0=100000.
	1071	p=ncobj.variables['PSFC'][:]
	1072	tk = (ncobj.variables['T2'][:] + 300.)(p/p0)*(2./7.)
	1073	qv = ncobj.variables['Q2'][:]
	1074
	1075	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
	1076	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
	1077
	1078	varNOcheckv = qv/data2
	1079	dimensions = ncobj.variables['PSFC'].dimensions
	1080	shape = ncobj.variables['PSFC'].shape
	1081
[343]	1082	elif varn == 'WRFT':
	1083	# To compute CF-times from WRF kind
	1084	#
	1085	import datetime as dt
	1086
	1087	times = ncobj.variables['Times']
	1088	dimt = times.shape[0]
	1089	varNOcheckv = np.zeros((dimt), dtype=np.float64)
	1090	self.unitsval = 'seconds since 1949-12-01 00:00:00'
	1091	refdate = datetimeStr_datetime('1949-12-01_00:00:00')
	1092
	1093	dimensions = tuple([ncobj.variables['Times'].dimensions[0]])
	1094	shape = tuple([dimt])
	1095
	1096	for it in range(dimt):
	1097	datevalS = datetimeStr_conversion(times[it,:], 'WRFdatetime', \
	1098	'YmdHMS')
	1099	dateval = dt.datetime.strptime(datevalS, '%Y%m%d%H%M%S')
	1100	difft = dateval - refdate
	1101	varNOcheckv[it] = difft.days3600.24. + difft.seconds + \
	1102	np.float(int(difft.microseconds/10.e6))
	1103
[337]	1104	elif varn == 'WRFt':
	1105	# print ' ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
	1106	p0=100000.
	1107	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	1108
	1109	varNOcheckv = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
	1110	dimensions = ncobj.variables['T'].dimensions
	1111	shape = ncobj.variables['P'].shape
	1112
[564]	1113	elif varn == 'TStd':
	1114	# print ' ' + main + ': computing dew-point temperature from TS as t and Tetens...'
	1115	# tacking from: http://en.wikipedia.org/wiki/Dew_point
[573]	1116	p=ncobj.variables['psfc'][:]
[564]	1117
	1118	temp = ncobj.variables['t'][:]
	1119
	1120	qv = ncobj.variables['q'][:]
	1121
	1122	tk = temp
	1123	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
	1124	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
	1125
	1126	rh = qv/data2
	1127
[609]	1128	pa = rh * data1
[564]	1129	varNOcheckv = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
	1130
	1131	dimensions = ncobj.variables['t'].dimensions
	1132	shape = ncobj.variables['t'].shape
	1133
[339]	1134	elif varn == 'WRFtd':
	1135	# print ' ' + main + ': computing dew-point temperature from WRF as inv_potT(T + 300) and Tetens...'
	1136	# tacking from: http://en.wikipedia.org/wiki/Dew_point
	1137	p0=100000.
	1138	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	1139
	1140	temp = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
	1141
	1142	qv = ncobj.variables['QVAPOR'][:]
	1143
	1144	tk = temp - 273.15
	1145	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
	1146	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
	1147
	1148	rh = qv/data2
	1149
[609]	1150	pa = rh * data1
[339]	1151	varNOcheckv = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
	1152
	1153	dimensions = ncobj.variables['T'].dimensions
	1154	shape = ncobj.variables['P'].shape
	1155
[505]	1156	elif varn == 'WRFwd':
	1157	# print ' ' + main + ': computing wind direction from WRF as ATAN2PI(V,U) ...'
[339]	1158	uwind = ncobj.variables['U'][:]
	1159	vwind = ncobj.variables['V'][:]
	1160	dx = uwind.shape[3]
	1161	dy = vwind.shape[2]
[505]	1162
[339]	1163	# de-staggering
	1164	ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
	1165	va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])
	1166
[505]	1167	theta = np.arctan2(ua, va)
[506]	1168	theta = np.where(theta < 0., theta + 2.*np.pi, theta)
	1169	varNOcheckv = 360.theta/(2.np.pi)
	1170
[339]	1171	dimensions = tuple(['Time','bottom_top','south_north','west_east'])
	1172	shape = ua.shape
	1173
[505]	1174	elif varn == 'WRFwds':
	1175	# print ' ' + main + ': computing surface wind direction from WRF as ATAN2(V,U) ...'
	1176	theta = np.arctan2(ncobj.variables['V10'][:], ncobj.variables['U10'][:])
[506]	1177	theta = np.where(theta < 0., theta + 2.*np.pi, theta)
[505]	1178
	1179	varNOcheckv = 360.theta/(2.np.pi)
	1180	dimensions = ncobj.variables['V10'].dimensions
	1181	shape = ncobj.variables['V10'].shape
	1182
	1183	elif varn == 'TSwds':
	1184	# print ' ' + main + ': computing surface wind direction from TSs as ATAN2(v,u) ...'
	1185	theta = np.arctan2(ncobj.variables['v'][:], ncobj.variables['u'][:])
[506]	1186	theta = np.where(theta < 0., theta + 2.*np.pi, theta)
[505]	1187
	1188	varNOcheckv = 360.theta/(2.np.pi)
	1189	dimensions = ncobj.variables['v'].dimensions
	1190	shape = ncobj.variables['v'].shape
	1191
	1192	elif varn == 'WRFws':
	1193	# print ' ' + main + ': computing wind speed from WRF as SQRT(U2 + V2) ...'
[339]	1194	uwind = ncobj.variables['U'][:]
	1195	vwind = ncobj.variables['V'][:]
	1196	dx = uwind.shape[3]
	1197	dy = vwind.shape[2]
[505]	1198
[339]	1199	# de-staggering
	1200	ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
	1201	va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])
	1202
[505]	1203	varNOcheckv = np.sqrt(uaua + vava)
[339]	1204	dimensions = tuple(['Time','bottom_top','south_north','west_east'])
	1205	shape = ua.shape
	1206
[505]	1207	elif varn == 'TSwss':
	1208	# print ' ' + main + ': computing surface wind speed from TSs as SQRT(u2 + v2) ...'
	1209	varNOcheckv = np.sqrt(ncobj.variables['u'][:]* \
	1210	ncobj.variables['u'][:] + ncobj.variables['v'][:]* \
	1211	ncobj.variables['v'][:])
	1212	dimensions = ncobj.variables['u'].dimensions
	1213	shape = ncobj.variables['u'].shape
	1214
	1215	elif varn == 'WRFwss':
	1216	# print ' ' + main + ': computing surface wind speed from WRF as SQRT(U2 + V2) ...'
	1217	varNOcheckv = np.sqrt(ncobj.variables['U10'][:]* \
	1218	ncobj.variables['U10'][:] + ncobj.variables['V10'][:]* \
	1219	ncobj.variables['V10'][:])
	1220	dimensions = ncobj.variables['U10'].dimensions
	1221	shape = ncobj.variables['U10'].shape
	1222
[337]	1223	elif varn == 'WRFz':
	1224	grav = 9.81
	1225	# print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
	1226	varNOcheckv = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/grav
	1227	dimensions = ncobj.variables['PH'].dimensions
	1228	shape = ncobj.variables['PH'].shape
	1229
	1230	else:
[339]	1231	print errormsg
[337]	1232	print ' ' + fname + ": variable '" + varn + "' nor ready !!"
	1233	quit(-1)
	1234
	1235	self.dimensions = dimensions
	1236	self.shape = shape
	1237	self.__values = varNOcheckv
	1238
	1239	def __getitem__(self,elem):
	1240	return self.__values[elem]
	1241
[492]	1242	def adding_station_desc(onc,stdesc):
	1243	""" Function to add a station description in a netCDF file
	1244	onc= netCDF object
	1245	stdesc= station description lon, lat, height
	1246	"""
	1247	fname = 'adding_station_desc'
	1248
[500]	1249	newvar = onc.createVariable( 'station', 'c', ('StrLength'))
[492]	1250	newvar[0:len(stdesc[0])] = stdesc[0]
	1251
	1252	newdim = onc.createDimension('nst',1)
	1253
	1254	if onc.variables.has_key('lon'):
	1255	print warnmsg
	1256	print ' ' + fname + ": variable 'lon' already exist !!"
	1257	print " renaming it as 'lonst'"
	1258	lonname = 'lonst'
	1259	else:
	1260	lonname = 'lon'
	1261
[500]	1262	newvar = onc.createVariable( lonname, 'f4', ('nst'))
[492]	1263	basicvardef(newvar, lonname, 'longitude', 'degrees_West' )
	1264	newvar[:] = stdesc[1]
	1265
	1266	if onc.variables.has_key('lat'):
	1267	print warnmsg
	1268	print ' ' + fname + ": variable 'lat' already exist !!"
	1269	print " renaming it as 'latst'"
	1270	latname = 'latst'
	1271	else:
	1272	latname = 'lat'
	1273
[500]	1274	newvar = onc.createVariable( latname, 'f4', ('nst'))
[492]	1275	basicvardef(newvar, lonname, 'latitude', 'degrees_North' )
	1276	newvar[:] = stdesc[2]
	1277
	1278	if onc.variables.has_key('height'):
	1279	print warnmsg
	1280	print ' ' + fname + ": variable 'height' already exist !!"
	1281	print " renaming it as 'heightst'"
	1282	heightname = 'heightst'
	1283	else:
	1284	heightname = 'height'
	1285
[500]	1286	newvar = onc.createVariable( heightname, 'f4', ('nst'))
[492]	1287	basicvardef(newvar, heightname, 'height above sea level', 'm' )
	1288	newvar[:] = stdesc[3]
	1289
	1290	return
	1291
[517]	1292	class Quantiles(object):
	1293	""" Class to provide quantiles from a given arrayof values
	1294	"""
	1295
	1296	def __init__(self, values, Nquants):
	1297	import numpy.ma as ma
	1298
	1299	if values is None:
	1300	self.quantilesv = None
	1301
	1302	else:
	1303	self.quantilesv = []
	1304
	1305	vals0 = values.flatten()
	1306	Nvalues = len(vals0)
	1307	vals = ma.masked_equal(vals0, None)
	1308	Nvals=len(vals.compressed())
	1309
	1310	sortedvals = sorted(vals.compressed())
	1311	for iq in range(Nquants):
	1312	self.quantilesv.append(sortedvals[int((Nvals-1)*iq/Nquants)])
	1313
	1314	self.quantilesv.append(sortedvals[Nvals-1])
	1315
	1316
[529]	1317	def getting_ValidationValues(okind, dt, ds, trjpos, ovs, ovo, tvalues, oFill, Ng, kvals):
[496]	1318	""" Function to get the values to validate accroding to the type of observation
	1319	okind= observational kind
[529]	1320	dt= initial number of values to retrieve
[496]	1321	ds= dictionary with the names of the dimensions (sim, obs)
	1322	trjpos= positions of the multi-stations (t, Y, X) or trajectory ([Z], Y, X)
	1323	ovs= object with the values of the simulation
	1324	ovs= object with the values of the observations
	1325	tvalues= temporal values (sim. time step, obs. time step, sim t value, obs t value, t diff)
[498]	1326	oFill= Fill Value for the observations
[496]	1327	Ng= number of grid points around the observation
[514]	1328	kvals= kind of values
	1329	'instantaneous': values are taken as instantaneous values
	1330	'tbackwardSmean': simulated values are taken as time averages from back to the point
	1331	'tbackwardOmean': observed values are taken as time averages from back to the point
[496]	1332	return:
	1333	sovalues= simulated values at the observation point and time
	1334	soSvalues= values Ngrid points around the simulated point
[498]	1335	soTtvalues= inital/ending period between two consecutive obsevations (for `single-station')
[496]	1336	trjs= trajectory on the simulation space
	1337	"""
	1338	fname = 'getting_ValidationValues'
	1339
	1340	sovalues = []
	1341
[514]	1342	if kvals == 'instantaneous':
[529]	1343	dtf = dt
[514]	1344	elif kvals == 'tbackwardSmean':
	1345	print ' ' + fname + ':',kvals,'!!'
[515]	1346	uniqt = np.unique(tvalues[:,3])
[529]	1347	dtf = len(uniqt)
	1348	print ' initially we got',dt,'values which will become',dtf
[516]	1349	postf = {}
[529]	1350	for it in range(dtf):
[516]	1351	if it == 0:
	1352	postf[uniqt[it]] = [0,0]
[517]	1353	elif it == 1:
	1354	posprev = postf[uniqt[it-1]][1]
	1355	posit = list(tvalues[:,3]).index(uniqt[it])
	1356	postf[uniqt[it]] = [posprev, posit+1]
[516]	1357	else:
	1358	posprev = postf[uniqt[it-1]][1]
[517]	1359	posit = list(tvalues[:,3]).index(uniqt[it])
[607]	1360	postf[uniqt[it]] = [posprev+1, np.min([posit+1,dt-1])]
[514]	1361	elif kvals == 'tbackwardOmean':
	1362	print ' ' + fname + ':',kvals,'!!'
[515]	1363	uniqt = np.unique(tvalues[:,2])
[529]	1364	dtf = len(uniqt)
	1365	print ' initially we got',dt,'values which will become',dtf
[527]	1366	print ' ==== NOT READY === '
	1367	quit(-1)
[514]	1368	else:
	1369	print errormsg
	1370	print ' ' + fname + ": kind of values '" + kvals + "' not ready!!"
	1371	quit(-1)
	1372
[496]	1373	# Simulated values spatially around
	1374	if ds.has_key('Z'):
[529]	1375	soSvalues = np.zeros((dt, Ng2+1, Ng2+1, Ng*2+1), dtype = np.float)
[496]	1376	if okind == 'trajectory':
[529]	1377	trjs = np.zeros((4,dt), dtype=int)
[496]	1378	else:
[498]	1379	trjs = None
[496]	1380	else:
[529]	1381	soSvalues = np.zeros((dt, Ng2+1, Ng2+1), dtype = np.float)
[496]	1382	if okind == 'trajectory':
[529]	1383	trjs = np.zeros((3,dt), dtype=int)
[496]	1384	else:
[498]	1385	trjs = None
[496]	1386
	1387	if okind == 'single-station':
[529]	1388	soTtvalues = np.zeros((dt,2), dtype=np.float)
[496]	1389	else:
	1390	None
	1391
	1392	if okind == 'multi-points':
[529]	1393	for it in range(dt):
[516]	1394	slicev = ds['X'][0] + ':' + str(trjpos[2,it]) + '\|' + ds['Y'][0] + \
	1395	':' + str(trjpos[1,it]) + '\|' + ds['T'][0]+ ':' + str(tvalues[it][0])
[496]	1396	slicevar, dimslice = slice_variable(ovs, slicev)
	1397	sovalues.append([ slicevar, ovo[tvalues[it][1]]])
[516]	1398	slicev = ds['X'][0] + ':' + str(trjpos[2,it]-Ng) + '@' + \
	1399	str(trjpos[2,it]+Ng) + '\|' + ds['Y'][0] + ':' + \
	1400	str(trjpos[1,it]-Ng) + '@' + str(trjpos[1,it]+Ng) + '\|' + \
[496]	1401	ds['T'][0]+':'+str(tvalues[it][0])
	1402	slicevar, dimslice = slice_variable(ovs, slicev)
	1403	soSvalues[it,:,:] = slicevar
	1404
	1405	elif okind == 'single-station':
[529]	1406	for it in range(dt):
[496]	1407	ito = int(tvalues[it,1])
	1408	if valdimsim.has_key('X') and valdimsim.has_key('Y'):
[516]	1409	slicev = ds['X'][0] + ':' + str(stationpos[1]) + '\|' + \
	1410	ds['Y'][0] + ':' + str(stationpos[0]) + '\|' + \
[496]	1411	ds['T'][0] + ':' + str(int(tvalues[it][0]))
	1412	else:
	1413	slicev = ds['T'][0] + ':' + str(int(tvalues[it][0]))
	1414	slicevar, dimslice = slice_variable(ovs, slicev)
	1415	if ovo[int(ito)] == oFill or ovo[int(ito)] == '--':
	1416	sovalues.append([ slicevar, fillValueF])
	1417	# elif ovo[int(ito)] != ovo[int(ito)]:
	1418	# sovalues.append([ slicevar, fillValueF])
	1419	else:
	1420	sovalues.append([ slicevar, ovo[int(ito)]])
	1421	if valdimsim.has_key('X') and valdimsim.has_key('Y'):
[516]	1422	slicev = ds['X'][0] + ':' + str(stationpos[1]-Ng) + '@' + \
	1423	str(stationpos[1]+Ng+1) + '\|' + ds['Y'][0] + ':' + \
	1424	str(stationpos[0]-Ng) + '@' + str(stationpos[0]+Ng+1) + '\|' + \
[496]	1425	ds['T'][0] + ':' + str(int(tvalues[it,0]))
	1426	else:
	1427	slicev = ds['T'][0] + ':' + str(int(tvalues[it][0]))
	1428	slicevar, dimslice = slice_variable(ovs, slicev)
	1429	soSvalues[it,:,:] = slicevar
	1430
	1431	if it == 0:
	1432	itoi = 0
	1433	itof = int(tvalues[it,1]) / 2
[529]	1434	elif it == dt-1:
[496]	1435	itoi = int( (ito + int(tvalues[it-1,1])) / 2)
	1436	itof = int(tvalues[it,1])
	1437	else:
	1438	itod = int( (ito - int(tvalues[it-1,1])) / 2 )
	1439	itoi = ito - itod
	1440	itod = int( (int(tvalues[it+1,1]) - ito) / 2 )
	1441	itof = ito + itod
	1442
	1443	soTtvalues[it,0] = valdimobs['T'][itoi]
	1444	soTtvalues[it,1] = valdimobs['T'][itof]
	1445
	1446	elif okind == 'trajectory':
	1447	if ds.has_key('Z'):
[529]	1448	for it in range(dt):
[496]	1449	ito = int(tvalues[it,1])
	1450	if notfound[ito] == 0:
[516]	1451	trjpos[2,ito] = index_mat(valdimsim['Z'][tvalues[it,0],:, \
[496]	1452	trjpos[1,ito],trjpos[0,ito]], valdimobs['Z'][ito])
[516]	1453	slicev = ds['X'][0]+':'+str(trjpos[0,ito]) + '\|' + \
	1454	ds['Y'][0]+':'+str(trjpos[1,ito]) + '\|' + \
	1455	ds['Z'][0]+':'+str(trjpos[2,ito]) + '\|' + \
[496]	1456	ds['T'][0]+':'+str(int(tvalues[it,0]))
	1457	slicevar, dimslice = slice_variable(ovs, slicev)
	1458	sovalues.append([ slicevar, ovo[int(ito)]])
	1459	minx = np.max([trjpos[0,ito]-Ng,0])
	1460	maxx = np.min([trjpos[0,ito]+Ng+1,ovs.shape[3]])
	1461	miny = np.max([trjpos[1,ito]-Ng,0])
	1462	maxy = np.min([trjpos[1,ito]+Ng+1,ovs.shape[2]])
	1463	minz = np.max([trjpos[2,ito]-Ng,0])
	1464	maxz = np.min([trjpos[2,ito]+Ng+1,ovs.shape[1]])
	1465
[516]	1466	slicev = ds['X'][0] + ':' + str(minx) + '@' + str(maxx) + '\|' + \
	1467	ds['Y'][0] + ':' + str(miny) + '@' + str(maxy) + '\|' + \
	1468	ds['Z'][0] + ':' + str(minz) + '@' + str(maxz) + '\|' + \
[496]	1469	ds['T'][0] + ':' + str(int(tvalues[it,0]))
	1470	slicevar, dimslice = slice_variable(ovs, slicev)
	1471
	1472	sliceS = []
	1473	sliceS.append(it)
	1474	sliceS.append(slice(0,maxz-minz))
	1475	sliceS.append(slice(0,maxy-miny))
	1476	sliceS.append(slice(0,maxx-minx))
	1477
	1478	soSvalues[tuple(sliceS)] = slicevar
	1479	if ivar == 0:
	1480	trjs[0,it] = trjpos[0,ito]
	1481	trjs[1,it] = trjpos[1,ito]
	1482	trjs[2,it] = trjpos[2,ito]
	1483	trjs[3,it] = tvalues[it,0]
	1484	else:
	1485	sovalues.append([fillValueF, fillValueF])
[516]	1486	soSvalues[it,:,:,:]= np.ones((Ng2+1,Ng2+1,Ng*2+1), \
[496]	1487	dtype = np.float)*fillValueF
	1488	# 2D trajectory
	1489	else:
[529]	1490	for it in range(dt):
[496]	1491	if notfound[it] == 0:
	1492	ito = tvalues[it,1]
[516]	1493	slicev = ds['X'][0]+':'+str(trjpos[2,ito]) + '\|' + \
	1494	ds['Y'][0]+':'+str(trjpos[1,ito]) + '\|' + \
[496]	1495	ds['T'][0]+':'+str(tvalues[ito,0])
	1496	slicevar, dimslice = slice_variable(ovs, slicev)
	1497	sovalues.append([ slicevar, ovo[tvalues[it,1]]])
[516]	1498	slicev = ds['X'][0] + ':' + str(trjpos[0,it]-Ng) + '@' + \
	1499	str(trjpos[0,it]+Ng) + '\|' + ds['Y'][0] + ':' + \
	1500	str(trjpos[1,it]-Ng) + '@' + str(trjpos[1,it]+Ng) + \
[496]	1501	'\|' + ds['T'][0] + ':' + str(tvalues[it,0])
	1502	slicevar, dimslice = slice_variable(ovs, slicev)
	1503	soSvalues[it,:,:] = slicevar
	1504	else:
	1505	sovalues.append([fillValue, fillValue])
[516]	1506	soSvalues[it,:,:] = np.ones((Ng2+1,Ng2+1), \
[496]	1507	dtype = np.float)*fillValueF
	1508	print sovalues[varsimobs][:][it]
	1509	else:
	1510	print errormsg
[497]	1511	print ' ' + fname + ": observatino kind '" + okind + "' not ready!!"
[496]	1512	quit(-1)
	1513
[516]	1514	# Re-arranging final values
	1515	##
	1516	if kvals == 'instantaneous':
[520]	1517	return sovalues, soSvalues, soTtvalues, trjs
[496]	1518
[516]	1519	elif kvals == 'tbackwardSmean':
	1520	fsovalues = []
	1521	if ds.has_key('Z'):
[529]	1522	fsoSvalues = np.zeros((dtf, Ng2+1, Ng2+1, Ng*2+1), dtype = np.float)
[516]	1523	if okind == 'trajectory':
[529]	1524	ftrjs = np.zeros((4,dtf), dtype=int)
[516]	1525	else:
	1526	ftrjs = None
	1527	else:
[529]	1528	fsoSvalues = np.zeros((dtf, Ng2+1, Ng2+1), dtype = np.float)
[516]	1529	if okind == 'trajectory':
[529]	1530	ftrjs = np.zeros((3,dtf), dtype=int)
[516]	1531	else:
	1532	ftrjs = None
	1533
	1534	if okind == 'single-station':
[529]	1535	fsoTtvalues = np.ones((dtf,2), dtype=np.float)*fillValueF
[516]	1536	else:
	1537	None
	1538
[529]	1539	for it in range(1,dtf):
[517]	1540	tv = uniqt[it]
	1541	intv = postf[tv]
[516]	1542
[517]	1543	# Temporal statistics
[583]	1544	if len(np.array(sovalues[intv[0]:intv[1]]).shape) != 1:
	1545	sovs = np.array(sovalues[intv[0]:intv[1]])[:,0]
	1546	minv = np.min(sovs)
	1547	maxv = np.max(sovs)
	1548	meanv = np.mean(sovs)
	1549	stdv = np.std(sovs)
[517]	1550
[584]	1551	fsovalues.append([meanv, np.array(sovalues[intv[0]:intv[1]])[0,1], \
[583]	1552	minv, maxv, stdv])
	1553	else:
[584]	1554	fsovalues.append([fillValueF, fillValueF, fillValueF, fillValueF, \
	1555	fillValueF, fillValueF])
[517]	1556	if ds.has_key('Z'):
[518]	1557	if okind == 'trajectory':
	1558	for ip in range(4):
	1559	ftrjs[ip,it] = np.mean(trjs[ip,intv[0]:intv[1]])
[517]	1560	for iz in range(2*Ng+1):
[519]	1561	for iy in range(2*Ng+1):
	1562	for ix in range(2*Ng+1):
	1563	fsoSvalues[it,iz,iy,ix] = np.mean(soSvalues[intv[0]: \
[517]	1564	intv[1],iz,iy,ix])
[518]	1565	else:
	1566	if okind == 'trajectory':
	1567	for ip in range(3):
	1568	ftrjs[ip,it] = np.mean(trjs[ip,intv[0]:intv[1]])
[519]	1569	for iy in range(2*Ng+1):
	1570	for ix in range(2*Ng+1):
	1571	fsoSvalues[it,iy,ix] = np.mean(soSvalues[intv[0]:intv[1], \
[518]	1572	iy,ix])
[519]	1573	fsoTtvalues[it,0] = soTtvalues[intv[0],0]
	1574	fsoTtvalues[it,1] = soTtvalues[intv[1],0]
[517]	1575
[527]	1576	return fsovalues, fsoSvalues, fsoTtvalues, ftrjs
[519]	1577
[516]	1578	elif kvals == 'tbackwardOmean':
	1579	print ' ' + fname + ':',kvals,'!!'
	1580	uniqt = np.unique(tvalues[:,2])
[529]	1581	dtf = len(uniqt)
	1582	print ' initially we got',dt,'values which will become',dtf
[516]	1583
[520]	1584	return
[496]	1585
	1586
[330]	1587	####### ###### ##### #### ### ## #
	1588
	1589	strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " + \
	1590	" and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')"
	1591
	1592	kindobs=['multi-points', 'single-station', 'trajectory']
	1593	strkObs="kind of observations; 'multi-points': multiple individual punctual obs " + \
	1594	"(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
	1595	"'trajectory': following a trajectory"
[337]	1596	simopers = ['sumc','subc','mulc','divc']
	1597	opersinf = 'sumc,[constant]: add [constant] to variables values; subc,[constant]: '+ \
	1598	'substract [constant] to variables values; mulc,[constant]: multipy by ' + \
	1599	'[constant] to variables values; divc,[constant]: divide by [constant] to ' + \
	1600	'variables values'
[505]	1601	varNOcheck = ['WRFdens', 'WRFght', 'WRFp', 'WRFrh', 'TSrhs', 'WRFrhs', 'WRFT', \
[572]	1602	'WRFt', 'TStd', 'WRFtd', 'WRFwd', 'TSwds', 'WRFwds', 'WRFws', 'TSwss', 'WRFwss', \
	1603	'WRFz']
[505]	1604	varNOcheckinf = "'WRFdens': air density from WRF variables; " + \
	1605	"'WRFght': geopotentiali height from WRF variables; " + \
	1606	"'WRFp': pressure from WRF variables; " + \
	1607	"'WRFrh': relative humidty fom WRF variables; " + \
[564]	1608	"'TSrhs': surface relative humidity from TS variables; " + \
[505]	1609	"'WRFrhs': surface relative humidity from WRF variables; " + \
	1610	"'WRFT': CF-time from WRF variables; " + \
	1611	"'WRFt': temperature from WRF variables; " + \
[564]	1612	"'TStd': dew-point temperature from TS variables; " + \
[505]	1613	"'WRFtd': dew-point temperature from WRF variables; " + \
	1614	"'WRFwd': wind direction from WRF variables; " + \
[564]	1615	"'TSwds': surface wind direction from TS variables; " + \
[505]	1616	"'WRFwds': surface wind direction from WRF variables; " + \
	1617	"'WRFws': wind speed from WRF variables; " + \
[564]	1618	"'TSwss': surface wind speed from TS variables; " + \
[505]	1619	"'WRFwss': surface wind speed from WRF variables; " + \
	1620	"'WRFz': height from WRF variables"
[330]	1621
[337]	1622	dimshelp = "[DIM]@[simdim]@[obsdim] ',' list of couples of dimensions names from " + \
	1623	"each source ([DIM]='X','Y','Z','T'; None, no value)"
	1624	vardimshelp = "[DIM]@[simvardim]@[obsvardim] ',' list of couples of variables " + \
	1625	"names with dimensions values from each source ([DIM]='X','Y','Z','T'; None, " + \
	1626	"no value, WRFdiagnosted variables also available: " + varNOcheckinf + ")"
	1627	varshelp="[simvar]@[obsvar]@[[oper]@[val]] ',' list of couples of variables to " + \
[505]	1628	"validate and if necessary operation and value (sim. values) available " + \
	1629	"operations: " + opersinf + " (WRFdiagnosted variables also available: " + \
	1630	varNOcheckinf + ")"
[337]	1631	statsn = ['minimum', 'maximum', 'mean', 'mean2', 'standard deviation']
[344]	1632	gstatsn = ['bias', 'simobs_mean', 'sim_obsmin', 'sim_obsmax', 'sim_obsmean', 'mae', \
[503]	1633	'rmse', 'r_pearsoncorr', 'p_pearsoncorr', 'deviation_of_residuals_SDR', \
	1634	'indef_of_efficiency_IOE', 'index_of_agreement_IOA', 'fractional_mean_bias_FMB']
[346]	1635	ostatsn = ['number of points', 'minimum', 'maximum', 'mean', 'mean2', \
	1636	'standard deviation']
[337]	1637
[330]	1638	parser = OptionParser()
[337]	1639	parser.add_option("-d", "--dimensions", dest="dims", help=dimshelp, metavar="VALUES")
[330]	1640	parser.add_option("-D", "--vardimensions", dest="vardims",
[337]	1641	help=vardimshelp, metavar="VALUES")
[330]	1642	parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs,
	1643	help=strkObs, metavar="FILE")
	1644	parser.add_option("-l", "--stationLocation", dest="stloc",
[505]	1645	help="name (\| for spaces), longitude, latitude and height of the station (only for 'single-station')",
[330]	1646	metavar="FILE")
	1647	parser.add_option("-o", "--observation", dest="fobs",
	1648	help="observations file to validate", metavar="FILE")
	1649	parser.add_option("-s", "--simulation", dest="fsim",
	1650	help="simulation file to validate", metavar="FILE")
[337]	1651	parser.add_option("-t", "--trajectoryfile", dest="trajf",
	1652	help="file with grid points of the trajectory in the simulation grid ('simtrj')",
	1653	metavar="FILE")
[330]	1654	parser.add_option("-v", "--variables", dest="vars",
[337]	1655	help=varshelp, metavar="VALUES")
[330]	1656
	1657	(opts, args) = parser.parse_args()
	1658
[502]	1659	####### ###### ##### #### ### ## #
	1660	# Number of different statistics according to the temporal coincidence
	1661	# 0: Exact time
[523]	1662	# 1: Simulation values closest to observed times
	1663	# 2: Simulation values between consecutive observed times
	1664	Nstsim = 3
[502]	1665
[523]	1666	stdescsim = ['E', 'C', 'B']
[524]	1667	prestdescsim = ['exact', 'closest', 'between']
[523]	1668	Lstdescsim = ['exact time', 'cloest time', 'between observational time-steps']
[502]	1669
[330]	1670	####### #######
	1671	## MAIN
	1672	#######
	1673
	1674	ofile='validation_sim.nc'
	1675
	1676	if opts.dims is None:
	1677	print errormsg
	1678	print ' ' + main + ': No list of dimensions are provided!!'
	1679	print ' a ',' list of values X@[dimxsim]@[dimxobs],...,T@[dimtsim]@[dimtobs]'+\
	1680	' is needed'
	1681	quit(-1)
	1682	else:
[340]	1683	simdims = {}
	1684	obsdims = {}
[330]	1685	print main +': couple of dimensions _______'
	1686	dims = {}
	1687	ds = opts.dims.split(',')
	1688	for d in ds:
	1689	dsecs = d.split('@')
	1690	if len(dsecs) != 3:
	1691	print errormsg
	1692	print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
	1693	print ' [DIM]@[dimnsim]@[dimnobs]'
	1694	quit(-1)
[500]	1695	if dsecs[1] != 'None':
	1696	dims[dsecs[0]] = [dsecs[1], dsecs[2]]
	1697	simdims[dsecs[0]] = dsecs[1]
	1698	obsdims[dsecs[0]] = dsecs[2]
[340]	1699
[500]	1700	print ' ',dsecs[0],':',dsecs[1],',',dsecs[2]
[330]	1701
	1702	if opts.vardims is None:
	1703	print errormsg
	1704	print ' ' + main + ': No list of variables with dimension values are provided!!'
	1705	print ' a ',' list of values X@[vardimxsim]@[vardimxobs],...,T@' + \
	1706	'[vardimtsim]@[vardimtobs] is needed'
	1707	quit(-1)
	1708	else:
	1709	print main +': couple of variable dimensions _______'
	1710	vardims = {}
	1711	ds = opts.vardims.split(',')
	1712	for d in ds:
	1713	dsecs = d.split('@')
	1714	if len(dsecs) != 3:
	1715	print errormsg
	1716	print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
	1717	print ' [DIM]@[vardimnsim]@[vardimnobs]'
	1718	quit(-1)
[500]	1719	if dsecs[1] != 'None':
	1720	vardims[dsecs[0]] = [dsecs[1], dsecs[2]]
	1721	print ' ',dsecs[0],':',dsecs[1],',',dsecs[2]
[330]	1722
	1723	if opts.obskind is None:
	1724	print errormsg
	1725	print ' ' + main + ': No kind of observations provided !!'
	1726	quit(-1)
	1727	else:
	1728	obskind = opts.obskind
	1729	if obskind == 'single-station':
	1730	if opts.stloc is None:
	1731	print errormsg
	1732	print ' ' + main + ': No station location provided !!'
	1733	quit(-1)
	1734	else:
[505]	1735	stationdesc = [opts.stloc.split(',')[0].replace('\|',' '), \
[492]	1736	np.float(opts.stloc.split(',')[1]), np.float(opts.stloc.split(',')[2]),\
	1737	np.float(opts.stloc.split(',')[3])]
[330]	1738
	1739	if opts.fobs is None:
	1740	print errormsg
	1741	print ' ' + main + ': No observations file is provided!!'
	1742	quit(-1)
	1743	else:
	1744	if not os.path.isfile(opts.fobs):
	1745	print errormsg
	1746	print ' ' + main + ": observations file '" + opts.fobs + "' does not exist !!"
	1747	quit(-1)
	1748
	1749	if opts.fsim is None:
	1750	print errormsg
	1751	print ' ' + main + ': No simulation file is provided!!'
	1752	quit(-1)
	1753	else:
	1754	if not os.path.isfile(opts.fsim):
	1755	print errormsg
	1756	print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
	1757	quit(-1)
	1758
	1759	if opts.vars is None:
	1760	print errormsg
	1761	print ' ' + main + ': No list of couples of variables is provided!!'
	1762	print ' a ',' list of values [varsim]@[varobs],... is needed'
	1763	quit(-1)
	1764	else:
	1765	valvars = []
[333]	1766	vs = opts.vars.split(',')
[330]	1767	for v in vs:
	1768	vsecs = v.split('@')
[333]	1769	if len(vsecs) < 2:
[330]	1770	print errormsg
	1771	print ' ' + main + ': wrong number of values in:',vsecs, \
	1772	' at least 2 are needed !!'
	1773	print ' [varsim]@[varobs]@[[oper][val]]'
	1774	quit(-1)
[333]	1775	if len(vsecs) > 2:
	1776	if not searchInlist(simopers,vsecs[2]):
	1777	print errormsg
	1778	print main + ": operation on simulation values '" + vsecs[2] + \
	1779	"' not ready !!"
	1780	quit(-1)
	1781
[330]	1782	valvars.append(vsecs)
	1783
	1784	# Openning observations trajectory
	1785	##
	1786	oobs = NetCDFFile(opts.fobs, 'r')
	1787
	1788	valdimobs = {}
	1789	for dn in dims:
	1790	print dn,':',dims[dn]
	1791	if dims[dn][1] != 'None':
	1792	if not oobs.dimensions.has_key(dims[dn][1]):
	1793	print errormsg
	1794	print ' ' + main + ": observations file does not have dimension '" + \
	1795	dims[dn][1] + "' !!"
	1796	quit(-1)
	1797	if vardims[dn][1] != 'None':
	1798	if not oobs.variables.has_key(vardims[dn][1]):
	1799	print errormsg
	1800	print ' ' + main + ": observations file does not have varibale " + \
	1801	"dimension '" + vardims[dn][1] + "' !!"
	1802	quit(-1)
	1803	valdimobs[dn] = oobs.variables[vardims[dn][1]][:]
	1804	else:
	1805	if dn == 'X':
[492]	1806	valdimobs[dn] = stationdesc[1]
[330]	1807	elif dn == 'Y':
[492]	1808	valdimobs[dn] = stationdesc[2]
[330]	1809	elif dn == 'Z':
[492]	1810	valdimobs[dn] = stationdesc[3]
[330]	1811
	1812	osim = NetCDFFile(opts.fsim, 'r')
	1813
	1814	valdimsim = {}
	1815	for dn in dims:
[464]	1816	if dims[dn][0] != 'None':
	1817	if not osim.dimensions.has_key(dims[dn][0]):
	1818	print errormsg
	1819	print ' ' + main + ": simulation file '" + opts.fsim + \
	1820	"' does not have dimension '" + dims[dn][0] + "' !!"
	1821	print ' it has: ',osim.dimensions
	1822	quit(-1)
[343]	1823
[464]	1824	if not osim.variables.has_key(vardims[dn][0]) and \
	1825	not searchInlist(varNOcheck,vardims[dn][0]):
	1826	print errormsg
	1827	print ' ' + main + ": simulation file '" + opts.fsim + \
	1828	"' does not have varibale dimension '" + vardims[dn][0] + "' !!"
	1829	print ' it has variables:',osim.variables
	1830	quit(-1)
	1831	if searchInlist(varNOcheck,vardims[dn][0]):
	1832	valdimsim[dn] = compute_varNOcheck(osim, vardims[dn][0])
[465]	1833	else:
[464]	1834	valdimsim[dn] = osim.variables[vardims[dn][0]][:]
[330]	1835
	1836	# General characteristics
[343]	1837	dimtobs = valdimobs['T'].shape[0]
	1838	dimtsim = valdimsim['T'].shape[0]
[330]	1839
	1840	print main +': observational time-steps:',dimtobs,'simulation:',dimtsim
	1841
[337]	1842	notfound = np.zeros((dimtobs), dtype=int)
	1843
[330]	1844	if obskind == 'multi-points':
	1845	trajpos = np.zeros((2,dimt),dtype=int)
[337]	1846	for it in range(dimtobs):
[330]	1847	trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'], \
	1848	[valdimobs['X'][it],valdimobss['Y'][it]])
	1849	elif obskind == 'single-station':
[498]	1850	trajpos = None
[333]	1851	stationpos = np.zeros((2), dtype=int)
[466]	1852	if valdimsim.has_key('X') and valdimsim.has_key('Y'):
	1853	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],[valdimobs['Y'], \
	1854	valdimobs['X']])
	1855	iid = 0
	1856	for idn in osim.variables[vardims['X'][0]].dimensions:
	1857	if idn == dims['X'][0]:
	1858	stationpos[1] = stsimpos[iid]
	1859	elif idn == dims['Y'][0]:
	1860	stationpos[0] = stsimpos[iid]
[333]	1861
[466]	1862	iid = iid + 1
	1863	print main + ': station point in simulation:', stationpos
	1864	print ' station position:',valdimobs['X'],',',valdimobs['Y']
	1865	print ' simulation coord.:',valdimsim['X'][tuple(stsimpos)],',', \
	1866	valdimsim['Y'][tuple(stsimpos)]
	1867	else:
	1868	print main + ': validation with two time-series !!'
[333]	1869
[330]	1870	elif obskind == 'trajectory':
[337]	1871	if opts.trajf is not None:
	1872	if not os.path.isfile(opts.fsim):
	1873	print errormsg
	1874	print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
	1875	quit(-1)
	1876	else:
	1877	otrjf = NetCDFFile(opts.fsim, 'r')
	1878	trajpos[0,:] = otrjf.variables['obssimtrj'][0]
	1879	trajpos[1,:] = otrjf.variables['obssimtrj'][1]
	1880	otrjf.close()
[330]	1881	else:
[337]	1882	if dims.has_key('Z'):
	1883	trajpos = np.zeros((3,dimtobs),dtype=int)
	1884	for it in range(dimtobs):
	1885	if np.mod(it*100./dimtobs,10.) == 0.:
	1886	print ' trajectory done: ',it*100./dimtobs,'%'
	1887	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \
	1888	[valdimobs['Y'][it],valdimobs['X'][it]])
	1889	stationpos = np.zeros((2), dtype=int)
	1890	iid = 0
	1891	for idn in osim.variables[vardims['X'][0]].dimensions:
	1892	if idn == dims['X'][0]:
	1893	stationpos[1] = stsimpos[iid]
	1894	elif idn == dims['Y'][0]:
	1895	stationpos[0] = stsimpos[iid]
	1896	iid = iid + 1
	1897	if stationpos[0] == 0 and stationpos[1] == 0: notfound[it] = 1
	1898
	1899	trajpos[0,it] = stationpos[0]
	1900	trajpos[1,it] = stationpos[1]
	1901	# In the simulation 'Z' varies with time ... non-hydrostatic model! ;)
	1902	# trajpos[2,it] = index_mat(valdimsim['Z'][it,:,stationpos[0], \
	1903	# stationpos[1]], valdimobs['Z'][it])
	1904	else:
	1905	trajpos = np.zeros((2,dimtobs),dtype=int)
	1906	for it in range(dimtobs):
	1907	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \
	1908	[valdimobs['Y'][it],valdimobss['X'][it]])
	1909	stationpos = np.zeros((2), dtype=int)
	1910	iid = 0
	1911	for idn in osim.variables[vardims['X'][0]].dimensions:
	1912	if idn == dims['X'][0]:
	1913	stationpos[1] = stsimpos[iid]
	1914	elif idn == dims['Y'][0]:
	1915	stationpos[0] = stsimpos[iid]
	1916	iid = iid + 1
[338]	1917	if stationpos[0] == 0 or stationpos[1] == 0: notfound[it] = 1
[330]	1918
[337]	1919	trajpos[0,it] = stationspos[0]
	1920	trajpos[1,it] = stationspos[1]
	1921
	1922	print main + ': not found',np.sum(notfound),'points of the trajectory'
	1923
[330]	1924	# Getting times
	1925	tobj = oobs.variables[vardims['T'][1]]
	1926	obstunits = tobj.getncattr('units')
[343]	1927	if vardims['T'][0] == 'WRFT':
	1928	tsim = valdimsim['T'][:]
	1929	simtunits = 'seconds since 1949-12-01 00:00:00'
	1930	else:
[467]	1931	tsim = osim.variables[vardims['T'][0]][:]
	1932	otsim = osim.variables[vardims['T'][0]]
	1933	simtunits = otsim.getncattr('units')
[330]	1934
[467]	1935	simobstimes = coincident_CFtimes(tsim, obstunits, simtunits)
[330]	1936
[601]	1937	iobst = CFtimes_datetime_NOfile([valdimobs['T'][0]], obstunits, 'standard')
	1938	fobst = CFtimes_datetime_NOfile([valdimobs['T'][dimtobs-1]], obstunits, 'standard')
	1939	isimt = CFtimes_datetime_NOfile([simobstimes[0]], obstunits, 'standard')
	1940	fsimt = CFtimes_datetime_NOfile([simobstimes[dimtsim-1]], obstunits, 'standard')
	1941	print 'check of times _______'
	1942	print ' * Observations'
	1943	print ' - first time:', iobst
	1944	print ' - last time:', fobst
	1945	print ' * Simulations'
	1946	print ' - first time:', isimt
	1947	print ' - last time:', fsimt
	1948
[493]	1949	#
	1950	## Looking for exact/near times
	1951	#
	1952
	1953	# Exact Coincident times
[333]	1954	##
[493]	1955	exacttvalues0 = []
	1956	for it in range(dimtsim):
	1957	ot = 0
	1958	for ito in range(ot,dimtobs-1):
[494]	1959	if valdimobs['T'][ito] == simobstimes[it]:
[493]	1960	ot = ito
	1961	exacttvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito]])
	1962
	1963	exacttvalues = np.array(exacttvalues0, dtype=np.float)
	1964
[579]	1965	if len(exacttvalues) == 0:
	1966	print warnmsg
	1967	print ' ' + main + ': no exact values found!'
	1968	Nexactt = 0
[601]	1969	# quit(-1)
[579]	1970	else:
	1971	Nexactt = len(exacttvalues[:,0])
	1972
[528]	1973	print main + ': found',Nexactt,'Temporal exact values in simulation and observations'
[511]	1974
[523]	1975	# Sim Closest times
	1976	##
	1977	Nsimt = 0
	1978	closesttvalues0 = []
	1979	closesttvalues0st = []
	1980	tsiminit = 0
	1981	tsimend = 0
	1982
	1983	dtsim = simobstimes[1] - simobstimes[0]
	1984
	1985	for it in range(dimtsim):
	1986	ot = 0
	1987	for ito in range(ot,dimtobs-1):
	1988	if np.abs(valdimobs['T'][ito] - simobstimes[it]) <= dtsim/2.:
	1989	ot = ito
	1990	tdist = simobstimes[it] - valdimobs['T'][ito]
	1991	closesttvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito], \
	1992	tdist])
	1993	Nsimt = Nsimt + 1
	1994	if tsiminit == 0: tsiminit=simobstimes[it]
	1995	tsimend = simobstimes[it]
	1996
	1997	closesttvalues = np.array(closesttvalues0, dtype=np.float)
	1998
[601]	1999	print 'closesttvales shape:',closesttvalues.shape
[523]	2000	Nclosest = len(closesttvalues[:,0])
	2001	print main + ': found',Nclosest,'Simulation time-values closest to observations'
	2002
	2003	if Nclosest == 0:
	2004	print warnmsg
	2005	print main + ': no cclosest times found !!'
	2006	print ' stopping it'
	2007	quit(-1)
	2008
[493]	2009	# Sim Coincident times
	2010	##
[511]	2011	Nsimt = 0
[337]	2012	coindtvalues0 = []
[511]	2013	coindtvalues0st = []
	2014	tsiminit = 0
	2015	tsimend = 0
	2016
[333]	2017	for it in range(dimtsim):
	2018	ot = 0
	2019	for ito in range(ot,dimtobs-1):
	2020	if valdimobs['T'][ito] < simobstimes[it] and valdimobs['T'][ito+1] > \
	2021	simobstimes[it]:
	2022	ot = ito
	2023	tdist = simobstimes[it] - valdimobs['T'][ito]
[511]	2024	coindtvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito], \
[337]	2025	tdist])
[511]	2026	Nsimt = Nsimt + 1
	2027	if tsiminit == 0: tsiminit=simobstimes[it]
	2028	tsimend = simobstimes[it]
	2029	elif simobstimes[it] > valdimobs['T'][ito+1]:
	2030	coindtvalues0st.append([Nsimt, ito, valdimobs['T'][ito],tsimend-tsiminit])
[330]	2031
[337]	2032	coindtvalues = np.array(coindtvalues0, dtype=np.float)
[511]	2033	coindtvaluesst = np.array(coindtvalues0st, dtype=np.float)
[337]	2034
	2035	Ncoindt = len(coindtvalues[:,0])
[493]	2036	print main + ': found',Ncoindt,'Simulation time-interval (within consecutive ' + \
	2037	'observed times) coincident times between simulation and observations'
[333]	2038
[337]	2039	if Ncoindt == 0:
	2040	print warnmsg
	2041	print main + ': no coincident times found !!'
	2042	print ' stopping it'
	2043	quit(-1)
	2044
[333]	2045	# Validating
	2046	##
	2047
	2048	onewnc = NetCDFFile(ofile, 'w')
	2049
	2050	# Dimensions
[524]	2051	for kst in range(Nstsim):
	2052	newdim = onewnc.createDimension(prestdescsim[kst] + 'time',None)
	2053	if stdescsim[kst] != 'E':
	2054	newdim = onewnc.createDimension(prestdescsim[kst] + 'obstime',None)
	2055
[346]	2056	newdim = onewnc.createDimension('bnds',2)
[333]	2057	newdim = onewnc.createDimension('couple',2)
	2058	newdim = onewnc.createDimension('StrLength',StringLength)
[337]	2059	newdim = onewnc.createDimension('xaround',Ngrid*2+1)
	2060	newdim = onewnc.createDimension('yaround',Ngrid*2+1)
[503]	2061	newdim = onewnc.createDimension('gstats',13)
[337]	2062	newdim = onewnc.createDimension('stats',5)
[346]	2063	newdim = onewnc.createDimension('tstats',6)
[523]	2064	newdim = onewnc.createDimension('Nstsim', 3)
[333]	2065
	2066	# Variable dimensions
	2067	##
	2068	newvar = onewnc.createVariable('couple', 'c', ('couple','StrLength'))
	2069	basicvardef(newvar, 'couple', 'couples of values', '-')
	2070	writing_str_nc(newvar, ['sim','obs'], StringLength)
	2071
[337]	2072	newvar = onewnc.createVariable('statistics', 'c', ('stats','StrLength'))
	2073	basicvardef(newvar, 'statistics', 'statitics from values', '-')
	2074	writing_str_nc(newvar, statsn, StringLength)
	2075
[344]	2076	newvar = onewnc.createVariable('gstatistics', 'c', ('gstats','StrLength'))
	2077	basicvardef(newvar, 'gstatistics', 'global statitics from values', '-')
	2078	writing_str_nc(newvar, gstatsn, StringLength)
	2079
	2080	newvar = onewnc.createVariable('tstatistics', 'c', ('tstats','StrLength'))
	2081	basicvardef(newvar, 'tstatistics', 'statitics from values along time', '-')
	2082	writing_str_nc(newvar, ostatsn, StringLength)
	2083
[502]	2084	newvar = onewnc.createVariable('ksimstatistics', 'c', ('Nstsim','StrLength'))
	2085	basicvardef(newvar, 'ksimstatistics', 'kind of simulated statitics', '-')
	2086	writing_str_nc(newvar, Lstdescsim, StringLength)
	2087
[337]	2088	if obskind == 'trajectory':
	2089	if dims.has_key('Z'):
	2090	newdim = onewnc.createDimension('trj',3)
	2091	else:
	2092	newdim = onewnc.createDimension('trj',2)
	2093
	2094	newvar = onewnc.createVariable('obssimtrj','i',('obstime','trj'))
	2095	basicvardef(newvar, 'obssimtrj', 'trajectory on the simulation grid', '-')
	2096	newvar[:] = trajpos.transpose()
	2097
	2098	if dims.has_key('Z'):
	2099	newdim = onewnc.createDimension('simtrj',4)
	2100	else:
	2101	newdim = onewnc.createDimension('simtrj',3)
	2102
[333]	2103	Nvars = len(valvars)
	2104	for ivar in range(Nvars):
	2105	simobsvalues = []
	2106
	2107	varsimobs = valvars[ivar][0] + '_' + valvars[ivar][1]
[340]	2108	print ' ' + varsimobs + '... .. .'
[333]	2109
	2110	if not oobs.variables.has_key(valvars[ivar][1]):
	2111	print errormsg
	2112	print ' ' + main + ": observations file has not '" + valvars[ivar][1] + \
	2113	"' !!"
	2114	quit(-1)
[337]	2115
[333]	2116	if not osim.variables.has_key(valvars[ivar][0]):
[337]	2117	if not searchInlist(varNOcheck, valvars[ivar][0]):
	2118	print errormsg
	2119	print ' ' + main + ": simulation file has not '" + valvars[ivar][0] + \
	2120	"' !!"
	2121	quit(-1)
	2122	else:
	2123	ovsim = compute_varNOcheck(osim, valvars[ivar][0])
	2124	else:
	2125	ovsim = osim.variables[valvars[ivar][0]]
[333]	2126
[340]	2127	for idn in ovsim.dimensions:
	2128	if not searchInlist(simdims.values(),idn):
	2129	print errormsg
	2130	print ' ' + main + ": dimension '" + idn + "' of variable '" + \
	2131	valvars[ivar][0] + "' not provided as reference coordinate [X,Y,Z,T] !!"
	2132	quit(-1)
	2133
[333]	2134	ovobs = oobs.variables[valvars[ivar][1]]
[468]	2135	if searchInlist(ovobs.ncattrs(),'_FillValue'):
	2136	oFillValue = ovobs.getncattr('_FillValue')
[337]	2137	else:
[498]	2138	oFillValue = None
[333]	2139
[523]	2140	for kst in range(Nstsim):
[529]	2141	simobsvalues = []
	2142
[524]	2143	timedn = prestdescsim[kst] + 'time'
[529]	2144	timeobsdn = prestdescsim[kst] + 'obstime'
[526]	2145	print ' ' + prestdescsim[kst] + ' ...'
[501]	2146
[581]	2147	if stdescsim[kst] == 'E':
[523]	2148	# Observed and simualted exact times
	2149	simobsvalues, simobsSvalues, simobsTtvalues, trjsim = \
	2150	getting_ValidationValues(obskind, Nexactt, dims, trajpos, ovsim, \
	2151	ovobs, exacttvalues, oFillValue, Ngrid, 'instantaneous')
	2152
[526]	2153	if ivar == 0:
	2154	vname = prestdescsim[kst] + 'time'
	2155	newvar = onewnc.createVariable(vname,'f8', (timedn))
	2156	basicvardef(newvar, vname, 'exact coincident time observations and '+\
	2157	'simulation', obstunits)
	2158	set_attribute(newvar, 'calendar', 'standard')
[581]	2159	if Nexactt == 0:
	2160	newvar[:] = np.float64(0.)
	2161	else:
	2162	newvar[:] = exacttvalues[:,3]
[583]	2163	if Nexactt == 0:
	2164	simobsSvalues = np.zeros((1,2), dtype=np.float)
	2165	simobsvalues = np.zeros((1,2), dtype=np.float)
[524]	2166
[526]	2167	dimt = Nexactt
	2168
	2169	elif stdescsim[kst] == 'C':
[523]	2170	# Simualted closest to Observed times
	2171	simobsvalues, simobsSvalues, simobsTtvalues, trjsim = \
	2172	getting_ValidationValues(obskind, Nclosest, dims, trajpos, ovsim, \
	2173	ovobs, closesttvalues, oFillValue, Ngrid, 'instantaneous')
[526]	2174	dimt = Nclosest
[523]	2175
[526]	2176	if ivar == 0:
	2177	vname = prestdescsim[kst] + 'time'
	2178	newvar = onewnc.createVariable(vname,'f8', (timedn))
	2179	basicvardef(newvar, vname, 'time simulations closest to observed ' + \
	2180	'values', obstunits )
	2181	set_attribute(newvar, 'calendar', 'standard')
	2182	newvar[:] = closesttvalues[:,2]
[523]	2183
[526]	2184	vname = prestdescsim[kst] + 'obstime'
	2185	newvar = onewnc.createVariable(vname,'f8', (vname))
	2186	basicvardef(newvar, vname, 'closest time observations', obstunits)
	2187	set_attribute(newvar, 'calendar', 'standard')
	2188	newvar[:] = closesttvalues[:,3]
[524]	2189
[526]	2190	elif stdescsim[kst] == 'B':
[501]	2191	# Observed values temporally around coincident times
[523]	2192	simobsvalues, simobsSvalues, simobsTtvalues, trjsim = \
	2193	getting_ValidationValues(obskind, Ncoindt, dims, trajpos, ovsim, ovobs,\
	2194	coindtvalues, oFillValue, Ngrid, 'tbackwardSmean')
[528]	2195	dimt = simobsSvalues.shape[0]
[344]	2196
[526]	2197	if ivar == 0:
	2198	vname = prestdescsim[kst] + 'time'
[529]	2199	newvar = onewnc.createVariable(vname,'f8', (timedn), \
	2200	fill_value = fillValueF)
[526]	2201	basicvardef(newvar, vname, 'simulation time between observations', \
[529]	2202	obstunits)
[526]	2203	set_attribute(newvar, 'calendar', 'standard')
	2204	set_attribute(newvar, 'bounds', 'time_bnds')
	2205	newvar[:] = simobsTtvalues[:,1]
[523]	2206
[526]	2207	vname = prestdescsim[kst] + 'obstime'
	2208	newvar = onewnc.createVariable(vname,'f8', (vname))
	2209	basicvardef(newvar, vname, 'observed between time', obstunits )
	2210	set_attribute(newvar, 'calendar', 'standard')
[528]	2211	newvar[:] = np.unique(coindtvalues[:,3])
[524]	2212
[498]	2213	# Re-arranging values...
[584]	2214	# For an incomprensible reason it is not working?
	2215	# arrayvals = np.array(simobsvalues)
	2216	arrayvals = np.zeros((len(simobsvalues),2), dtype=np.float)
	2217	for it in range(len(simobsvalues)):
[610]	2218	arrayvals[it,:] = simobsvalues[it][0:2]
[584]	2219
[524]	2220	if len(valvars[ivar]) > 2:
	2221	const=np.float(valvars[ivar][3])
	2222	if valvars[ivar][2] == 'sumc':
	2223	simobsSvalues = simobsSvalues + const
	2224	arrayvals[:,0] = arrayvals[:,0] + const
	2225	elif valvars[ivar][2] == 'subc':
	2226	simobsSvalues = simobsSvalues - const
	2227	arrayvals[:,0] = arrayvals[:,0] - const
	2228	elif valvars[ivar][2] == 'mulc':
	2229	simobsSvalues = simobsSvalues * const
	2230	arrayvals[:,0] = arrayvals[:,0] * const
	2231	elif valvars[ivar][2] == 'divc':
	2232	simobsSvalues = simobsSvalues / const
	2233	arrayvals[:,0] = arrayvals[:,0] / const
	2234	else:
	2235	print errormsg
	2236	print ' ' + fname + ": operation '"+valvars[ivar][2]+"' not ready!!"
	2237	quit(-1)
[339]	2238
[610]	2239	# for it in range(len(arrayvals[:,0])):
	2240	# print it,arrayvals[it,:],':',simobsvalues[it]
	2241	# quit()
	2242
[528]	2243	if kst == 0:
	2244	simstats = np.zeros((Nstsim,5), dtype=np.float)
	2245	obsstats = np.zeros((Nstsim,5), dtype=np.float)
	2246	simobsstats = np.zeros((Nstsim,13), dtype=np.float)
	2247
[343]	2248	# statisics sim
[524]	2249	simstats[kst,0] = np.min(arrayvals[:,0])
	2250	simstats[kst,1] = np.max(arrayvals[:,0])
	2251	simstats[kst,2] = np.mean(arrayvals[:,0])
	2252	simstats[kst,3] = np.mean(arrayvals[:,0]*arrayvals[:,0])
	2253	simstats[kst,4] = np.sqrt(simstats[kst,3] - simstats[kst,2]*simstats[kst,2])
[343]	2254
	2255	# statisics obs
[501]	2256	# Masking 'nan'
[524]	2257	obsmask0 = np.where(arrayvals[:,1] != arrayvals[:,1], fillValueF, \
	2258	arrayvals[:,1])
[503]	2259
[524]	2260	obsmask = ma.masked_equal(obsmask0, fillValueF)
	2261	obsmask2 = obsmask*obsmask
[343]	2262
[524]	2263	obsstats[kst,0] = obsmask.min()
	2264	obsstats[kst,1] = obsmask.max()
	2265	obsstats[kst,2] = obsmask.mean()
	2266	obsstats[kst,3] = obsmask2.mean()
	2267	obsstats[kst,4] = np.sqrt(obsstats[kst,3] - obsstats[kst,2]*obsstats[kst,2])
[503]	2268
[343]	2269	# Statistics sim-obs
[526]	2270	diffvals = np.zeros((dimt), dtype=np.float)
[343]	2271
[524]	2272	diffvals = arrayvals[:,0] - obsmask
[344]	2273
[524]	2274	diff2vals = diffvals*diffvals
	2275	sumdiff = diffvals.sum()
	2276	sumdiff2 = diff2vals.sum()
[343]	2277
[524]	2278	simobsstats[kst,0] = simstats[kst,0] - obsstats[kst,0]
	2279	simobsstats[kst,1] = np.mean(arrayvals[:,0]*obsmask)
	2280	simobsstats[kst,2] = diffvals.min()
	2281	simobsstats[kst,3] = diffvals.max()
	2282	simobsstats[kst,4] = diffvals.mean()
	2283	simobsstats[kst,5] = np.abs(diffvals).mean()
	2284	simobsstats[kst,6] = np.sqrt(diff2vals.mean())
	2285	simobsstats[kst,7], simobsstats[kst,8] = sts.mstats.pearsonr(arrayvals[:,0], \
	2286	obsmask)
[502]	2287	# From:
[503]	2288	#Willmott, C. J. 1981. 'On the validation of models. Physical Geography', 2, 184-194
	2289	#Willmott, C. J. (1984). 'On the evaluation of model performance in physical
	2290	# geography'. Spatial Statistics and Models, G. L. Gaile and C. J. Willmott, eds.,
	2291	# 443-460
	2292	#Willmott, C. J., S. G. Ackleson, R. E. Davis, J. J. Feddema, K. M. Klink, D. R.
	2293	# Legates, J. O'Donnell, and C. M. Rowe (1985), 'Statistics for the Evaluation and
	2294	# Comparison of Models', J. Geophys. Res., 90(C5), 8995-9005
	2295	#Legates, D. R., and G. J. McCabe Jr. (1999), 'Evaluating the Use of "Goodness-of-Fit"
	2296	# Measures in Hydrologic and Hydroclimatic Model Validation', Water Resour. Res.,
	2297	# 35(1), 233-241
[502]	2298	#
	2299	# Deviation of residuals (SDR)
[524]	2300	simobsstats[kst,9] = np.mean(np.sqrt(np.abs((diffvals-simobsstats[kst,0])* \
	2301	(diffvals-obsstats[kst,0]))))
[502]	2302	# Index of Efficiency (IOE)
[524]	2303	obsbias2series = (obsmask - obsstats[kst,0])*(obsmask - obsstats[kst,0])
	2304	sumobsbias2series = obsbias2series.sum()
[502]	2305
[524]	2306	simobsstats[kst,10] = 1. - sumdiff2/(sumobsbias2series)
[502]	2307	# Index of Agreement (IOA)
[524]	2308	simbias2series = arrayvals[:,0] - obsstats[kst,0]
	2309	obsbias2series = obsmask - obsstats[kst,0]
[502]	2310
[524]	2311	abssimbias2series = np.abs(simbias2series)
	2312	absobsbias2series = np.where(obsbias2series<0, -obsbias2series, \
	2313	obsbias2series)
[502]	2314
[524]	2315	abssimobsbias2series = (abssimbias2series+absobsbias2series)*( \
	2316	abssimbias2series + absobsbias2series)
[502]	2317
[524]	2318	simobsstats[kst,11] = 1. - sumdiff2/(abssimobsbias2series.sum())
[502]	2319	# Fractional Mean Bias (FMB)
[526]	2320	simobsstats[kst,12]=(simstats[kst,0]-obsstats[kst,0])/(0.5*(simstats[kst,0] +\
[524]	2321	obsstats[kst,0]))
[502]	2322
[524]	2323	# Statistics around sim values
[528]	2324	aroundstats = np.zeros((5,dimt), dtype=np.float)
[526]	2325	for it in range(dimt):
[524]	2326	aroundstats[0,it] = np.min(simobsSvalues[it,])
	2327	aroundstats[1,it] = np.max(simobsSvalues[it,])
	2328	aroundstats[2,it] = np.mean(simobsSvalues[it,])
	2329	aroundstats[3,it] = np.mean(simobsSvalues[it,]*simobsSvalues[it,])
	2330	aroundstats[4,it] = np.sqrt(aroundstats[3,it] - aroundstats[2,it]* \
	2331	aroundstats[2,it])
[514]	2332
[524]	2333	# sim Values to netCDF
	2334	newvar = onewnc.createVariable(valvars[ivar][0] + '_' + stdescsim[kst] + \
	2335	'sim', 'f', (timedn), fill_value=fillValueF)
	2336	descvar = prestdescsim[kst] + ' time simulated: ' + valvars[ivar][0]
	2337	basicvardef(newvar, valvars[ivar][0], descvar, ovobs.getncattr('units'))
[581]	2338	if stdescsim[kst] == 'E' and Nexactt == 0:
	2339	newvar[:] = fillValueF
	2340	else:
	2341	newvar[:] = arrayvals[:,0]
[337]	2342
[524]	2343	# obs Values to netCDF
[529]	2344	if stdescsim[kst] != 'E':
	2345	newvar = onewnc.createVariable(valvars[ivar][1] + '_' + stdescsim[kst] + \
	2346	'obs', 'f', (timeobsdn), fill_value=fillValueF)
	2347	else:
	2348	newvar = onewnc.createVariable(valvars[ivar][1] + '_' + stdescsim[kst] + \
[524]	2349	'obs', 'f', (timedn), fill_value=fillValueF)
[529]	2350
	2351	descvar = prestdescsim[kst] + ' time observed: ' + valvars[ivar][1]
[524]	2352	basicvardef(newvar, valvars[ivar][1], descvar, ovobs.getncattr('units'))
[511]	2353
[581]	2354	if stdescsim[kst] == 'E' and Nexactt == 0:
	2355	newvar[:] = fillValueF
	2356	else:
	2357	newvar[:] = arrayvals[:,1]
	2358
[610]	2359	# Including statistics of the between simulated values
	2360	if stdescsim[kst] == 'B':
	2361	if not searchInlist(onewnc.dimensions, 'Bstats'):
	2362	newdim = onewnc.createDimension('Bstats',3)
	2363
	2364	newvar = onewnc.createVariable('Bstats','c',('Bstats','StrLength'))
	2365	descvar = prestdescsim[kst] + ' time simulated: ' + valvars[ivar][0]
	2366	basicvardef(newvar, 'Bstats', 'Between values statistics', \
	2367	ovobs.getncattr('units'))
	2368	Bstsvals=['minimum','maximum','stdandard deviation']
	2369	writing_str_nc(newvar,Bstsvals,StringLength)
	2370
	2371	newvar = onewnc.createVariable(valvars[ivar][0] + '_Bstats_sim', 'f', \
	2372	(timedn,'Bstats'), fill_value=fillValueF)
	2373	descvar = prestdescsim[kst] + ' time simulated: ' + valvars[ivar][0] + \
	2374	' statistics'
	2375	basicvardef(newvar, valvars[ivar][0], descvar, ovobs.getncattr('units'))
	2376	newvar[0,:] = [fillValueF, fillValueF, fillValueF]
	2377	for it in range(dimt-1):
	2378	newvar[it+1,:] = simobsvalues[it][2:5]
	2379
[524]	2380	# Around values
	2381	if not onewnc.variables.has_key(valvars[ivar][0] + 'around'):
	2382	vname = prestdescsim[kst] + valvars[ivar][0] + 'around'
[514]	2383	else:
[524]	2384	vname = prestdescsim[kst] + valvars[ivar][0] + 'Around'
[514]	2385
[341]	2386	if dims.has_key('Z'):
	2387	if not onewnc.dimensions.has_key('zaround'):
	2388	newdim = onewnc.createDimension('zaround',Ngrid*2+1)
[524]	2389	newvar = onewnc.createVariable(vname, 'f', (timedn,'zaround', \
	2390	'yaround','xaround'), fill_value=fillValueF)
[341]	2391	else:
[524]	2392	newvar = onewnc.createVariable(vname, 'f', (timedn,'yaround','xaround'), \
	2393	fill_value=fillValueF)
[339]	2394
[524]	2395	descvar = prestdescsim[kst] + 'around simulated values +/- grid values: ' + \
	2396	valvars[ivar][0]
	2397	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
[337]	2398
[581]	2399	if stdescsim[kst] == 'E' and Nexactt == 0:
	2400	newvar[:] = np.ones((0,Ngrid2+1,Ngrid2+1))*fillValueF
	2401	else:
	2402	newvar[:] = simobsSvalues
[344]	2403
[581]	2404
[529]	2405	# around sim Statistics
[524]	2406	if not searchInlist(onewnc.variables,prestdescsim[kst] + valvars[ivar][0] + \
	2407	'staround'):
	2408	vname = prestdescsim[kst] + valvars[ivar][0] + 'staround'
	2409	else:
	2410	vname = prestdescsim[kst] + valvars[ivar][0] + 'Staround'
[344]	2411
[524]	2412	newvar = onewnc.createVariable(vname, 'f', (timedn,'stats'), \
[344]	2413	fill_value=fillValueF)
[524]	2414	descvar = prestdescsim[kst] + ' around (' + str(Ngrid) + ', ' + str(Ngrid) +\
[514]	2415	') simulated statisitcs: ' + valvars[ivar][0]
[524]	2416	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
[514]	2417	set_attribute(newvar, 'cell_methods', 'time_bnds')
[581]	2418	if stdescsim[kst] == 'E' and Nexactt == 0:
	2419	newvar[:] = np.ones((0,5))*fillValueF
	2420	else:
	2421	newvar[:] = aroundstats.transpose()
[337]	2422
[524]	2423	if stdescsim[kst] == 'B':
	2424	if not searchInlist(onewnc.variables, 'time_bnds'):
[526]	2425	newvar = onewnc.createVariable('time_bnds','f8',(timedn,'bnds'))
[524]	2426	basicvardef(newvar, 'time_bnds', timedn, obstunits )
	2427	set_attribute(newvar, 'calendar', 'standard')
	2428	newvar[:] = simobsTtvalues
[346]	2429
[524]	2430	# sim Statistics
	2431	if not searchInlist(onewnc.variables,valvars[ivar][0] + 'stsim'):
	2432	vname = valvars[ivar][0] + 'stsim'
	2433	else:
	2434	vname = valvars[ivar][0] + 'stSim'
[333]	2435
[524]	2436	newvar = onewnc.createVariable(vname, 'f', ('Nstsim', 'stats'), \
	2437	fill_value=fillValueF)
	2438	descvar = 'simulated statisitcs: ' + valvars[ivar][0]
	2439	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
	2440	newvar[:] = simstats
	2441
	2442	# obs Statistics
	2443	if not searchInlist(onewnc.variables,valvars[ivar][1] + 'stobs'):
	2444	vname = valvars[ivar][1] + 'stobs'
	2445	else:
	2446	vname = valvars[ivar][1] + 'stObs'
	2447
	2448	newvar = onewnc.createVariable(vname, 'f', ('Nstsim', 'stats'), \
	2449	fill_value=fillValueF)
	2450	descvar = 'observed statisitcs: ' + valvars[ivar][1]
	2451	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
	2452	newvar[:] = obsstats
	2453
	2454	# sim-obs Statistics
	2455	if not searchInlist(onewnc.variables,varsimobs + 'st'):
	2456	vname = varsimobs + 'st'
	2457	else:
	2458	vname = varSimobs + 'st'
	2459
	2460	newvar = onewnc.createVariable(vname, 'f', ('Nstsim', 'gstats'), \
	2461	fill_value=fillValueF)
	2462	descvar = 'simulated-observed tatisitcs: ' + varsimobs
	2463	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
	2464	newvar[:] = simobsstats
	2465
	2466	onewnc.sync()
	2467
[500]	2468	if trjsim is not None:
[512]	2469	newvar = onewnc.createVariable('simtrj','i',('betweentime','simtrj'))
[500]	2470	basicvardef(newvar,'simtrj','coordinates [X,Y,Z,T] of the coincident ' + \
	2471	'trajectory in sim', obstunits)
	2472	newvar[:] = trjsim.transpose()
[337]	2473
[492]	2474	# Adding three variables with the station location, longitude, latitude and height
	2475	if obskind == 'single-station':
[500]	2476	adding_station_desc(onewnc,stationdesc)
[492]	2477
[333]	2478	# Global attributes
	2479	##
	2480	set_attribute(onewnc,'author_nc','Lluis Fita')
	2481	set_attribute(onewnc,'institution_nc','Laboratoire de Meteorology Dynamique, ' + \
	2482	'LMD-Jussieu, UPMC, Paris')
	2483	set_attribute(onewnc,'country_nc','France')
	2484	set_attribute(onewnc,'script_nc',main)
	2485	set_attribute(onewnc,'version_script',version)
	2486	set_attribute(onewnc,'information', \
	2487	'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html')
	2488	set_attribute(onewnc,'simfile',opts.fsim)
	2489	set_attribute(onewnc,'obsfile',opts.fobs)
	2490
	2491	onewnc.sync()
	2492	onewnc.close()
	2493
	2494	print main + ": successfull writting of '" + ofile + "' !!"

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

Context Navigation

source: lmdz_wrf/trunk/tools/validation_sim.py @ 1444

Download in other formats: