Context Navigation

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

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

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

Download in other formats: