Context Navigation

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

Last change on this file since 521 was 521, checked in by lfita, 9 years ago
Removing debugging prints and related quit
File size: 89.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
	16	main = 'validarion_sim.py'
	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
[339]	831	'WRFtd': dew-point temperature from WRF variables
[505]	832	'WRFwd': wind direction from WRF variables
	833	'TSwds': surface wind direction from TS variables
	834	'WRFwds': surface wind direction from WRF variables
[339]	835	'WRFws': wind speed from WRF variables
[505]	836	'TSwss': surface wind speed from TS variables
	837	'WRFwss': surface wind speed from WRF variables
[337]	838	'WRFz': height from WRF variables
	839	"""
	840	fname = 'compute_varNOcheck'
	841
	842	def __init__(self, ncobj, varn):
	843
	844	if ncobj is None:
	845	self = None
	846	self.dimensions = None
	847	self.shape = None
	848	self.__values = None
	849	else:
	850	if varn == 'WRFdens':
	851	# print ' ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \
	852	# 'DNW)/g ...'
	853	grav = 9.81
	854
	855	# Just we need in in absolute values: Size of the central grid cell
	856	## dxval = ncobj.getncattr('DX')
	857	## dyval = ncobj.getncattr('DY')
	858	## mapfac = ncobj.variables['MAPFAC_M'][:]
	859	## area = dxvaldyvalmapfac
	860	dimensions = ncobj.variables['MU'].dimensions
	861	shape = ncobj.variables['MU'].shape
	862
	863	mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
	864	dnw = ncobj.variables['DNW'][:]
	865
	866	varNOcheckv = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \
	867	dtype=np.float)
	868	levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
	869
	870	for it in range(mu.shape[0]):
	871	for iz in range(dnw.shape[1]):
	872	levval.fill(np.abs(dnw[it,iz]))
	873	varNOcheck[it,iz,:,:] = levval
	874	varNOcheck[it,iz,:,:] = mu[it,:,:]*varNOcheck[it,iz,:,:]/grav
	875
	876	elif varn == 'WRFght':
	877	# print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
	878	varNOcheckv = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
	879	dimensions = ncobj.variables['PH'].dimensions
	880	shape = ncobj.variables['PH'].shape
	881
	882	elif varn == 'WRFp':
	883	# print ' ' + fname + ': Retrieving pressure value from WRF as P + PB'
	884	varNOcheckv = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	885	dimensions = ncobj.variables['P'].dimensions
	886	shape = ncobj.variables['P'].shape
	887
	888	elif varn == 'WRFrh':
	889	# print ' ' + main + ": computing relative humidity from WRF as 'Tetens'" +\
	890	# ' equation (T,P) ...'
	891	p0=100000.
	892	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
[505]	893	tk = (ncobj.variables['T'][:])(p/p0)*(2./7.)
[337]	894	qv = ncobj.variables['QVAPOR'][:]
	895
	896	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
	897	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
	898
	899	varNOcheckv = qv/data2
	900	dimensions = ncobj.variables['P'].dimensions
	901	shape = ncobj.variables['P'].shape
	902
[505]	903	elif varn == 'TSrhs':
	904	# print ' ' + main + ": computing surface relative humidity from TSs as 'Tetens'" +\
	905	# ' equation (T,P) ...'
	906	p0=100000.
	907	p=ncobj.variables['psfc'][:]
	908	tk = (ncobj.variables['t'][:])(p/p0)*(2./7.)
	909	qv = ncobj.variables['q'][:]
	910
	911	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
	912	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
	913
	914	varNOcheckv = qv/data2
	915	dimensions = ncobj.variables['psfc'].dimensions
	916	shape = ncobj.variables['psfc'].shape
	917
	918	elif varn == 'WRFrhs':
	919	# print ' ' + main + ": computing surface relative humidity from WRF as 'Tetens'" +\
	920	# ' equation (T,P) ...'
	921	p0=100000.
	922	p=ncobj.variables['PSFC'][:]
	923	tk = (ncobj.variables['T2'][:] + 300.)(p/p0)*(2./7.)
	924	qv = ncobj.variables['Q2'][:]
	925
	926	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
	927	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
	928
	929	varNOcheckv = qv/data2
	930	dimensions = ncobj.variables['PSFC'].dimensions
	931	shape = ncobj.variables['PSFC'].shape
	932
[343]	933	elif varn == 'WRFT':
	934	# To compute CF-times from WRF kind
	935	#
	936	import datetime as dt
	937
	938	times = ncobj.variables['Times']
	939	dimt = times.shape[0]
	940	varNOcheckv = np.zeros((dimt), dtype=np.float64)
	941	self.unitsval = 'seconds since 1949-12-01 00:00:00'
	942	refdate = datetimeStr_datetime('1949-12-01_00:00:00')
	943
	944	dimensions = tuple([ncobj.variables['Times'].dimensions[0]])
	945	shape = tuple([dimt])
	946
	947	for it in range(dimt):
	948	datevalS = datetimeStr_conversion(times[it,:], 'WRFdatetime', \
	949	'YmdHMS')
	950	dateval = dt.datetime.strptime(datevalS, '%Y%m%d%H%M%S')
	951	difft = dateval - refdate
	952	varNOcheckv[it] = difft.days3600.24. + difft.seconds + \
	953	np.float(int(difft.microseconds/10.e6))
	954
[337]	955	elif varn == 'WRFt':
	956	# print ' ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
	957	p0=100000.
	958	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	959
	960	varNOcheckv = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
	961	dimensions = ncobj.variables['T'].dimensions
	962	shape = ncobj.variables['P'].shape
	963
[339]	964	elif varn == 'WRFtd':
	965	# print ' ' + main + ': computing dew-point temperature from WRF as inv_potT(T + 300) and Tetens...'
	966	# tacking from: http://en.wikipedia.org/wiki/Dew_point
	967	p0=100000.
	968	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	969
	970	temp = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
	971
	972	qv = ncobj.variables['QVAPOR'][:]
	973
	974	tk = temp - 273.15
	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['P'].shape
	985
[505]	986	elif varn == 'WRFwd':
	987	# print ' ' + main + ': computing wind direction from WRF as ATAN2PI(V,U) ...'
[339]	988	uwind = ncobj.variables['U'][:]
	989	vwind = ncobj.variables['V'][:]
	990	dx = uwind.shape[3]
	991	dy = vwind.shape[2]
[505]	992
[339]	993	# de-staggering
	994	ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
	995	va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])
	996
[505]	997	theta = np.arctan2(ua, va)
[506]	998	theta = np.where(theta < 0., theta + 2.*np.pi, theta)
	999	varNOcheckv = 360.theta/(2.np.pi)
	1000
[339]	1001	dimensions = tuple(['Time','bottom_top','south_north','west_east'])
	1002	shape = ua.shape
	1003
[505]	1004	elif varn == 'WRFwds':
	1005	# print ' ' + main + ': computing surface wind direction from WRF as ATAN2(V,U) ...'
	1006	theta = np.arctan2(ncobj.variables['V10'][:], ncobj.variables['U10'][:])
[506]	1007	theta = np.where(theta < 0., theta + 2.*np.pi, theta)
[505]	1008
	1009	varNOcheckv = 360.theta/(2.np.pi)
	1010	dimensions = ncobj.variables['V10'].dimensions
	1011	shape = ncobj.variables['V10'].shape
	1012
	1013	elif varn == 'TSwds':
	1014	# print ' ' + main + ': computing surface wind direction from TSs as ATAN2(v,u) ...'
	1015	theta = np.arctan2(ncobj.variables['v'][:], ncobj.variables['u'][:])
[506]	1016	theta = np.where(theta < 0., theta + 2.*np.pi, theta)
[505]	1017
	1018	varNOcheckv = 360.theta/(2.np.pi)
	1019	dimensions = ncobj.variables['v'].dimensions
	1020	shape = ncobj.variables['v'].shape
	1021
	1022	elif varn == 'WRFws':
	1023	# print ' ' + main + ': computing wind speed from WRF as SQRT(U2 + V2) ...'
[339]	1024	uwind = ncobj.variables['U'][:]
	1025	vwind = ncobj.variables['V'][:]
	1026	dx = uwind.shape[3]
	1027	dy = vwind.shape[2]
[505]	1028
[339]	1029	# de-staggering
	1030	ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
	1031	va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])
	1032
[505]	1033	varNOcheckv = np.sqrt(uaua + vava)
[339]	1034	dimensions = tuple(['Time','bottom_top','south_north','west_east'])
	1035	shape = ua.shape
	1036
[505]	1037	elif varn == 'TSwss':
	1038	# print ' ' + main + ': computing surface wind speed from TSs as SQRT(u2 + v2) ...'
	1039	varNOcheckv = np.sqrt(ncobj.variables['u'][:]* \
	1040	ncobj.variables['u'][:] + ncobj.variables['v'][:]* \
	1041	ncobj.variables['v'][:])
	1042	dimensions = ncobj.variables['u'].dimensions
	1043	shape = ncobj.variables['u'].shape
	1044
	1045	elif varn == 'WRFwss':
	1046	# print ' ' + main + ': computing surface wind speed from WRF as SQRT(U2 + V2) ...'
	1047	varNOcheckv = np.sqrt(ncobj.variables['U10'][:]* \
	1048	ncobj.variables['U10'][:] + ncobj.variables['V10'][:]* \
	1049	ncobj.variables['V10'][:])
	1050	dimensions = ncobj.variables['U10'].dimensions
	1051	shape = ncobj.variables['U10'].shape
	1052
[337]	1053	elif varn == 'WRFz':
	1054	grav = 9.81
	1055	# print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
	1056	varNOcheckv = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/grav
	1057	dimensions = ncobj.variables['PH'].dimensions
	1058	shape = ncobj.variables['PH'].shape
	1059
	1060	else:
[339]	1061	print errormsg
[337]	1062	print ' ' + fname + ": variable '" + varn + "' nor ready !!"
	1063	quit(-1)
	1064
	1065	self.dimensions = dimensions
	1066	self.shape = shape
	1067	self.__values = varNOcheckv
	1068
	1069	def __getitem__(self,elem):
	1070	return self.__values[elem]
	1071
[492]	1072	def adding_station_desc(onc,stdesc):
	1073	""" Function to add a station description in a netCDF file
	1074	onc= netCDF object
	1075	stdesc= station description lon, lat, height
	1076	"""
	1077	fname = 'adding_station_desc'
	1078
[500]	1079	newvar = onc.createVariable( 'station', 'c', ('StrLength'))
[492]	1080	newvar[0:len(stdesc[0])] = stdesc[0]
	1081
	1082	newdim = onc.createDimension('nst',1)
	1083
	1084	if onc.variables.has_key('lon'):
	1085	print warnmsg
	1086	print ' ' + fname + ": variable 'lon' already exist !!"
	1087	print " renaming it as 'lonst'"
	1088	lonname = 'lonst'
	1089	else:
	1090	lonname = 'lon'
	1091
[500]	1092	newvar = onc.createVariable( lonname, 'f4', ('nst'))
[492]	1093	basicvardef(newvar, lonname, 'longitude', 'degrees_West' )
	1094	newvar[:] = stdesc[1]
	1095
	1096	if onc.variables.has_key('lat'):
	1097	print warnmsg
	1098	print ' ' + fname + ": variable 'lat' already exist !!"
	1099	print " renaming it as 'latst'"
	1100	latname = 'latst'
	1101	else:
	1102	latname = 'lat'
	1103
[500]	1104	newvar = onc.createVariable( latname, 'f4', ('nst'))
[492]	1105	basicvardef(newvar, lonname, 'latitude', 'degrees_North' )
	1106	newvar[:] = stdesc[2]
	1107
	1108	if onc.variables.has_key('height'):
	1109	print warnmsg
	1110	print ' ' + fname + ": variable 'height' already exist !!"
	1111	print " renaming it as 'heightst'"
	1112	heightname = 'heightst'
	1113	else:
	1114	heightname = 'height'
	1115
[500]	1116	newvar = onc.createVariable( heightname, 'f4', ('nst'))
[492]	1117	basicvardef(newvar, heightname, 'height above sea level', 'm' )
	1118	newvar[:] = stdesc[3]
	1119
	1120	return
	1121
[517]	1122	class Quantiles(object):
	1123	""" Class to provide quantiles from a given arrayof values
	1124	"""
	1125
	1126	def __init__(self, values, Nquants):
	1127	import numpy.ma as ma
	1128
	1129	if values is None:
	1130	self.quantilesv = None
	1131
	1132	else:
	1133	self.quantilesv = []
	1134
	1135	vals0 = values.flatten()
	1136	Nvalues = len(vals0)
	1137	vals = ma.masked_equal(vals0, None)
	1138	Nvals=len(vals.compressed())
	1139
	1140	sortedvals = sorted(vals.compressed())
	1141	for iq in range(Nquants):
	1142	self.quantilesv.append(sortedvals[int((Nvals-1)*iq/Nquants)])
	1143
	1144	self.quantilesv.append(sortedvals[Nvals-1])
	1145
	1146
[516]	1147	def getting_ValidationValues(okind, dimt, ds, trjpos, ovs, ovo, tvalues, oFill, Ng, kvals):
[496]	1148	""" Function to get the values to validate accroding to the type of observation
	1149	okind= observational kind
[516]	1150	dimt= initial number of values to retrieve
[496]	1151	ds= dictionary with the names of the dimensions (sim, obs)
	1152	trjpos= positions of the multi-stations (t, Y, X) or trajectory ([Z], Y, X)
	1153	ovs= object with the values of the simulation
	1154	ovs= object with the values of the observations
	1155	tvalues= temporal values (sim. time step, obs. time step, sim t value, obs t value, t diff)
[498]	1156	oFill= Fill Value for the observations
[496]	1157	Ng= number of grid points around the observation
[514]	1158	kvals= kind of values
	1159	'instantaneous': values are taken as instantaneous values
	1160	'tbackwardSmean': simulated values are taken as time averages from back to the point
	1161	'tbackwardOmean': observed values are taken as time averages from back to the point
[496]	1162	return:
	1163	sovalues= simulated values at the observation point and time
	1164	soSvalues= values Ngrid points around the simulated point
[498]	1165	soTtvalues= inital/ending period between two consecutive obsevations (for `single-station')
[496]	1166	trjs= trajectory on the simulation space
	1167	"""
	1168	fname = 'getting_ValidationValues'
	1169
	1170	sovalues = []
	1171
[514]	1172	if kvals == 'instantaneous':
[516]	1173	dimtf = dimt
[514]	1174	elif kvals == 'tbackwardSmean':
	1175	print ' ' + fname + ':',kvals,'!!'
[515]	1176	uniqt = np.unique(tvalues[:,3])
[516]	1177	dimtf = len(uniqt)
	1178	print ' initially we got',dimt,'values which will become',dimtf
	1179	postf = {}
	1180	for it in range(dimtf):
	1181	if it == 0:
	1182	postf[uniqt[it]] = [0,0]
[517]	1183	elif it == 1:
	1184	posprev = postf[uniqt[it-1]][1]
	1185	posit = list(tvalues[:,3]).index(uniqt[it])
	1186	postf[uniqt[it]] = [posprev, posit+1]
[516]	1187	else:
	1188	posprev = postf[uniqt[it-1]][1]
[517]	1189	posit = list(tvalues[:,3]).index(uniqt[it])
	1190	postf[uniqt[it]] = [posprev+1, posit+1]
[514]	1191	elif kvals == 'tbackwardOmean':
	1192	print ' ' + fname + ':',kvals,'!!'
[515]	1193	uniqt = np.unique(tvalues[:,2])
[516]	1194	dimtf = len(uniqt)
	1195	print ' initially we got',dimt,'values which will become',dimtf
[514]	1196	else:
	1197	print errormsg
	1198	print ' ' + fname + ": kind of values '" + kvals + "' not ready!!"
	1199	quit(-1)
	1200
[496]	1201	# Simulated values spatially around
	1202	if ds.has_key('Z'):
[516]	1203	soSvalues = np.zeros((dimt, Ng2+1, Ng2+1, Ng*2+1), dtype = np.float)
[496]	1204	if okind == 'trajectory':
	1205	trjs = np.zeros((4,dimt), dtype=int)
	1206	else:
[498]	1207	trjs = None
[496]	1208	else:
	1209	soSvalues = np.zeros((dimt, Ng2+1, Ng2+1), dtype = np.float)
	1210	if okind == 'trajectory':
	1211	trjs = np.zeros((3,dimt), dtype=int)
	1212	else:
[498]	1213	trjs = None
[496]	1214
	1215	if okind == 'single-station':
	1216	soTtvalues = np.zeros((dimt,2), dtype=np.float)
	1217	else:
	1218	None
	1219
	1220	if okind == 'multi-points':
	1221	for it in range(dimt):
[516]	1222	slicev = ds['X'][0] + ':' + str(trjpos[2,it]) + '\|' + ds['Y'][0] + \
	1223	':' + str(trjpos[1,it]) + '\|' + ds['T'][0]+ ':' + str(tvalues[it][0])
[496]	1224	slicevar, dimslice = slice_variable(ovs, slicev)
	1225	sovalues.append([ slicevar, ovo[tvalues[it][1]]])
[516]	1226	slicev = ds['X'][0] + ':' + str(trjpos[2,it]-Ng) + '@' + \
	1227	str(trjpos[2,it]+Ng) + '\|' + ds['Y'][0] + ':' + \
	1228	str(trjpos[1,it]-Ng) + '@' + str(trjpos[1,it]+Ng) + '\|' + \
[496]	1229	ds['T'][0]+':'+str(tvalues[it][0])
	1230	slicevar, dimslice = slice_variable(ovs, slicev)
	1231	soSvalues[it,:,:] = slicevar
	1232
	1233	elif okind == 'single-station':
	1234	for it in range(dimt):
	1235	ito = int(tvalues[it,1])
	1236	if valdimsim.has_key('X') and valdimsim.has_key('Y'):
[516]	1237	slicev = ds['X'][0] + ':' + str(stationpos[1]) + '\|' + \
	1238	ds['Y'][0] + ':' + str(stationpos[0]) + '\|' + \
[496]	1239	ds['T'][0] + ':' + str(int(tvalues[it][0]))
	1240	else:
	1241	slicev = ds['T'][0] + ':' + str(int(tvalues[it][0]))
	1242	slicevar, dimslice = slice_variable(ovs, slicev)
	1243	if ovo[int(ito)] == oFill or ovo[int(ito)] == '--':
	1244	sovalues.append([ slicevar, fillValueF])
	1245	# elif ovo[int(ito)] != ovo[int(ito)]:
	1246	# sovalues.append([ slicevar, fillValueF])
	1247	else:
	1248	sovalues.append([ slicevar, ovo[int(ito)]])
	1249	if valdimsim.has_key('X') and valdimsim.has_key('Y'):
[516]	1250	slicev = ds['X'][0] + ':' + str(stationpos[1]-Ng) + '@' + \
	1251	str(stationpos[1]+Ng+1) + '\|' + ds['Y'][0] + ':' + \
	1252	str(stationpos[0]-Ng) + '@' + str(stationpos[0]+Ng+1) + '\|' + \
[496]	1253	ds['T'][0] + ':' + str(int(tvalues[it,0]))
	1254	else:
	1255	slicev = ds['T'][0] + ':' + str(int(tvalues[it][0]))
	1256	slicevar, dimslice = slice_variable(ovs, slicev)
	1257	soSvalues[it,:,:] = slicevar
	1258
	1259	if it == 0:
	1260	itoi = 0
	1261	itof = int(tvalues[it,1]) / 2
	1262	elif it == dimt-1:
	1263	itoi = int( (ito + int(tvalues[it-1,1])) / 2)
	1264	itof = int(tvalues[it,1])
	1265	else:
	1266	itod = int( (ito - int(tvalues[it-1,1])) / 2 )
	1267	itoi = ito - itod
	1268	itod = int( (int(tvalues[it+1,1]) - ito) / 2 )
	1269	itof = ito + itod
	1270
	1271	soTtvalues[it,0] = valdimobs['T'][itoi]
	1272	soTtvalues[it,1] = valdimobs['T'][itof]
	1273
	1274	elif okind == 'trajectory':
	1275	if ds.has_key('Z'):
	1276	for it in range(dimt):
	1277	ito = int(tvalues[it,1])
	1278	if notfound[ito] == 0:
[516]	1279	trjpos[2,ito] = index_mat(valdimsim['Z'][tvalues[it,0],:, \
[496]	1280	trjpos[1,ito],trjpos[0,ito]], valdimobs['Z'][ito])
[516]	1281	slicev = ds['X'][0]+':'+str(trjpos[0,ito]) + '\|' + \
	1282	ds['Y'][0]+':'+str(trjpos[1,ito]) + '\|' + \
	1283	ds['Z'][0]+':'+str(trjpos[2,ito]) + '\|' + \
[496]	1284	ds['T'][0]+':'+str(int(tvalues[it,0]))
	1285	slicevar, dimslice = slice_variable(ovs, slicev)
	1286	sovalues.append([ slicevar, ovo[int(ito)]])
	1287	minx = np.max([trjpos[0,ito]-Ng,0])
	1288	maxx = np.min([trjpos[0,ito]+Ng+1,ovs.shape[3]])
	1289	miny = np.max([trjpos[1,ito]-Ng,0])
	1290	maxy = np.min([trjpos[1,ito]+Ng+1,ovs.shape[2]])
	1291	minz = np.max([trjpos[2,ito]-Ng,0])
	1292	maxz = np.min([trjpos[2,ito]+Ng+1,ovs.shape[1]])
	1293
[516]	1294	slicev = ds['X'][0] + ':' + str(minx) + '@' + str(maxx) + '\|' + \
	1295	ds['Y'][0] + ':' + str(miny) + '@' + str(maxy) + '\|' + \
	1296	ds['Z'][0] + ':' + str(minz) + '@' + str(maxz) + '\|' + \
[496]	1297	ds['T'][0] + ':' + str(int(tvalues[it,0]))
	1298	slicevar, dimslice = slice_variable(ovs, slicev)
	1299
	1300	sliceS = []
	1301	sliceS.append(it)
	1302	sliceS.append(slice(0,maxz-minz))
	1303	sliceS.append(slice(0,maxy-miny))
	1304	sliceS.append(slice(0,maxx-minx))
	1305
	1306	soSvalues[tuple(sliceS)] = slicevar
	1307	if ivar == 0:
	1308	trjs[0,it] = trjpos[0,ito]
	1309	trjs[1,it] = trjpos[1,ito]
	1310	trjs[2,it] = trjpos[2,ito]
	1311	trjs[3,it] = tvalues[it,0]
	1312	else:
	1313	sovalues.append([fillValueF, fillValueF])
[516]	1314	soSvalues[it,:,:,:]= np.ones((Ng2+1,Ng2+1,Ng*2+1), \
[496]	1315	dtype = np.float)*fillValueF
	1316	# 2D trajectory
	1317	else:
	1318	for it in range(dimt):
	1319	if notfound[it] == 0:
	1320	ito = tvalues[it,1]
[516]	1321	slicev = ds['X'][0]+':'+str(trjpos[2,ito]) + '\|' + \
	1322	ds['Y'][0]+':'+str(trjpos[1,ito]) + '\|' + \
[496]	1323	ds['T'][0]+':'+str(tvalues[ito,0])
	1324	slicevar, dimslice = slice_variable(ovs, slicev)
	1325	sovalues.append([ slicevar, ovo[tvalues[it,1]]])
[516]	1326	slicev = ds['X'][0] + ':' + str(trjpos[0,it]-Ng) + '@' + \
	1327	str(trjpos[0,it]+Ng) + '\|' + ds['Y'][0] + ':' + \
	1328	str(trjpos[1,it]-Ng) + '@' + str(trjpos[1,it]+Ng) + \
[496]	1329	'\|' + ds['T'][0] + ':' + str(tvalues[it,0])
	1330	slicevar, dimslice = slice_variable(ovs, slicev)
	1331	soSvalues[it,:,:] = slicevar
	1332	else:
	1333	sovalues.append([fillValue, fillValue])
[516]	1334	soSvalues[it,:,:] = np.ones((Ng2+1,Ng2+1), \
[496]	1335	dtype = np.float)*fillValueF
	1336	print sovalues[varsimobs][:][it]
	1337	else:
	1338	print errormsg
[497]	1339	print ' ' + fname + ": observatino kind '" + okind + "' not ready!!"
[496]	1340	quit(-1)
	1341
[516]	1342	# Re-arranging final values
	1343	##
	1344	if kvals == 'instantaneous':
[520]	1345	return sovalues, soSvalues, soTtvalues, trjs
[496]	1346
[516]	1347	elif kvals == 'tbackwardSmean':
	1348	fsovalues = []
	1349	if ds.has_key('Z'):
	1350	fsoSvalues = np.zeros((dimtf, Ng2+1, Ng2+1, Ng*2+1), dtype = np.float)
	1351	if okind == 'trajectory':
	1352	ftrjs = np.zeros((4,dimtf), dtype=int)
	1353	else:
	1354	ftrjs = None
	1355	else:
	1356	fsoSvalues = np.zeros((dimtf, Ng2+1, Ng2+1), dtype = np.float)
	1357	if okind == 'trajectory':
	1358	ftrjs = np.zeros((3,dimtf), dtype=int)
	1359	else:
	1360	ftrjs = None
	1361
	1362	if okind == 'single-station':
	1363	fsoTtvalues = np.zeros((dimtf,2), dtype=np.float)
	1364	else:
	1365	None
	1366
[517]	1367	for it in range(1,dimtf):
	1368	tv = uniqt[it]
	1369	intv = postf[tv]
[516]	1370
[517]	1371	# Temporal statistics
	1372	sovs = np.array(sovalues[intv[0]:intv[1]])[:,0]
	1373	minv = np.min(sovs)
	1374	maxv = np.max(sovs)
	1375	meanv = np.mean(sovs)
	1376	stdv = np.std(sovs)
	1377
	1378	fsovalues.append([meanv, np.array(sovalues[intv[0]:intv[1]])[0,1], minv, \
	1379	maxv, stdv])
	1380	if ds.has_key('Z'):
[518]	1381	if okind == 'trajectory':
	1382	for ip in range(4):
	1383	ftrjs[ip,it] = np.mean(trjs[ip,intv[0]:intv[1]])
[517]	1384	for iz in range(2*Ng+1):
[519]	1385	for iy in range(2*Ng+1):
	1386	for ix in range(2*Ng+1):
	1387	fsoSvalues[it,iz,iy,ix] = np.mean(soSvalues[intv[0]: \
[517]	1388	intv[1],iz,iy,ix])
[518]	1389	else:
	1390	if okind == 'trajectory':
	1391	for ip in range(3):
	1392	ftrjs[ip,it] = np.mean(trjs[ip,intv[0]:intv[1]])
[519]	1393	for iy in range(2*Ng+1):
	1394	for ix in range(2*Ng+1):
	1395	fsoSvalues[it,iy,ix] = np.mean(soSvalues[intv[0]:intv[1], \
[518]	1396	iy,ix])
[519]	1397	fsoTtvalues[it,0] = soTtvalues[intv[0],0]
	1398	fsoTtvalues[it,1] = soTtvalues[intv[1],0]
[517]	1399
[520]	1400	return sovalues, soSvalues, soTtvalues, trjs
[519]	1401
[516]	1402	elif kvals == 'tbackwardOmean':
	1403	print ' ' + fname + ':',kvals,'!!'
	1404	uniqt = np.unique(tvalues[:,2])
	1405	dimtf = len(uniqt)
	1406	print ' initially we got',dimt,'values which will become',dimtf
	1407
	1408
[520]	1409	return
[496]	1410
	1411
[330]	1412	####### ###### ##### #### ### ## #
	1413
	1414	strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " + \
	1415	" and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')"
	1416
	1417	kindobs=['multi-points', 'single-station', 'trajectory']
	1418	strkObs="kind of observations; 'multi-points': multiple individual punctual obs " + \
	1419	"(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
	1420	"'trajectory': following a trajectory"
[337]	1421	simopers = ['sumc','subc','mulc','divc']
	1422	opersinf = 'sumc,[constant]: add [constant] to variables values; subc,[constant]: '+ \
	1423	'substract [constant] to variables values; mulc,[constant]: multipy by ' + \
	1424	'[constant] to variables values; divc,[constant]: divide by [constant] to ' + \
	1425	'variables values'
[505]	1426	varNOcheck = ['WRFdens', 'WRFght', 'WRFp', 'WRFrh', 'TSrhs', 'WRFrhs', 'WRFT', \
	1427	'WRFt', 'WRFtd', 'WRFwd', 'TSwds', 'WRFwds', 'WRFws', 'TSwss', 'WRFwss', 'WRFz']
	1428	varNOcheckinf = "'WRFdens': air density from WRF variables; " + \
	1429	"'WRFght': geopotentiali height from WRF variables; " + \
	1430	"'WRFp': pressure from WRF variables; " + \
	1431	"'WRFrh': relative humidty fom WRF variables; " + \
	1432	"'WRFrhs': surface relative humidity from WRF variables; " + \
	1433	"'WRFT': CF-time from WRF variables; " + \
	1434	"'WRFt': temperature from WRF variables; " + \
	1435	"'WRFtd': dew-point temperature from WRF variables; " + \
	1436	"'WRFwd': wind direction from WRF variables; " + \
	1437	"'WRFwds': surface wind direction from WRF variables; " + \
	1438	"'WRFws': wind speed from WRF variables; " + \
	1439	"'WRFwss': surface wind speed from WRF variables; " + \
	1440	"'WRFz': height from WRF variables"
[330]	1441
[337]	1442	dimshelp = "[DIM]@[simdim]@[obsdim] ',' list of couples of dimensions names from " + \
	1443	"each source ([DIM]='X','Y','Z','T'; None, no value)"
	1444	vardimshelp = "[DIM]@[simvardim]@[obsvardim] ',' list of couples of variables " + \
	1445	"names with dimensions values from each source ([DIM]='X','Y','Z','T'; None, " + \
	1446	"no value, WRFdiagnosted variables also available: " + varNOcheckinf + ")"
	1447	varshelp="[simvar]@[obsvar]@[[oper]@[val]] ',' list of couples of variables to " + \
[505]	1448	"validate and if necessary operation and value (sim. values) available " + \
	1449	"operations: " + opersinf + " (WRFdiagnosted variables also available: " + \
	1450	varNOcheckinf + ")"
[337]	1451	statsn = ['minimum', 'maximum', 'mean', 'mean2', 'standard deviation']
[344]	1452	gstatsn = ['bias', 'simobs_mean', 'sim_obsmin', 'sim_obsmax', 'sim_obsmean', 'mae', \
[503]	1453	'rmse', 'r_pearsoncorr', 'p_pearsoncorr', 'deviation_of_residuals_SDR', \
	1454	'indef_of_efficiency_IOE', 'index_of_agreement_IOA', 'fractional_mean_bias_FMB']
[346]	1455	ostatsn = ['number of points', 'minimum', 'maximum', 'mean', 'mean2', \
	1456	'standard deviation']
[337]	1457
[330]	1458	parser = OptionParser()
[337]	1459	parser.add_option("-d", "--dimensions", dest="dims", help=dimshelp, metavar="VALUES")
[330]	1460	parser.add_option("-D", "--vardimensions", dest="vardims",
[337]	1461	help=vardimshelp, metavar="VALUES")
[330]	1462	parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs,
	1463	help=strkObs, metavar="FILE")
	1464	parser.add_option("-l", "--stationLocation", dest="stloc",
[505]	1465	help="name (\| for spaces), longitude, latitude and height of the station (only for 'single-station')",
[330]	1466	metavar="FILE")
	1467	parser.add_option("-o", "--observation", dest="fobs",
	1468	help="observations file to validate", metavar="FILE")
	1469	parser.add_option("-s", "--simulation", dest="fsim",
	1470	help="simulation file to validate", metavar="FILE")
[337]	1471	parser.add_option("-t", "--trajectoryfile", dest="trajf",
	1472	help="file with grid points of the trajectory in the simulation grid ('simtrj')",
	1473	metavar="FILE")
[330]	1474	parser.add_option("-v", "--variables", dest="vars",
[337]	1475	help=varshelp, metavar="VALUES")
[330]	1476
	1477	(opts, args) = parser.parse_args()
	1478
[502]	1479	####### ###### ##### #### ### ## #
	1480	# Number of different statistics according to the temporal coincidence
	1481	# 0: Exact time
	1482	# 1: Simulation values between consecutive observed times
	1483	Nstsim = 2
	1484
	1485	stdescsim = ['E', 'B']
	1486	Lstdescsim = ['exact time', 'between observational time-steps']
	1487
[330]	1488	####### #######
	1489	## MAIN
	1490	#######
	1491
	1492	ofile='validation_sim.nc'
	1493
	1494	if opts.dims is None:
	1495	print errormsg
	1496	print ' ' + main + ': No list of dimensions are provided!!'
	1497	print ' a ',' list of values X@[dimxsim]@[dimxobs],...,T@[dimtsim]@[dimtobs]'+\
	1498	' is needed'
	1499	quit(-1)
	1500	else:
[340]	1501	simdims = {}
	1502	obsdims = {}
[330]	1503	print main +': couple of dimensions _______'
	1504	dims = {}
	1505	ds = opts.dims.split(',')
	1506	for d in ds:
	1507	dsecs = d.split('@')
	1508	if len(dsecs) != 3:
	1509	print errormsg
	1510	print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
	1511	print ' [DIM]@[dimnsim]@[dimnobs]'
	1512	quit(-1)
[500]	1513	if dsecs[1] != 'None':
	1514	dims[dsecs[0]] = [dsecs[1], dsecs[2]]
	1515	simdims[dsecs[0]] = dsecs[1]
	1516	obsdims[dsecs[0]] = dsecs[2]
[340]	1517
[500]	1518	print ' ',dsecs[0],':',dsecs[1],',',dsecs[2]
[330]	1519
	1520	if opts.vardims is None:
	1521	print errormsg
	1522	print ' ' + main + ': No list of variables with dimension values are provided!!'
	1523	print ' a ',' list of values X@[vardimxsim]@[vardimxobs],...,T@' + \
	1524	'[vardimtsim]@[vardimtobs] is needed'
	1525	quit(-1)
	1526	else:
	1527	print main +': couple of variable dimensions _______'
	1528	vardims = {}
	1529	ds = opts.vardims.split(',')
	1530	for d in ds:
	1531	dsecs = d.split('@')
	1532	if len(dsecs) != 3:
	1533	print errormsg
	1534	print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
	1535	print ' [DIM]@[vardimnsim]@[vardimnobs]'
	1536	quit(-1)
[500]	1537	if dsecs[1] != 'None':
	1538	vardims[dsecs[0]] = [dsecs[1], dsecs[2]]
	1539	print ' ',dsecs[0],':',dsecs[1],',',dsecs[2]
[330]	1540
	1541	if opts.obskind is None:
	1542	print errormsg
	1543	print ' ' + main + ': No kind of observations provided !!'
	1544	quit(-1)
	1545	else:
	1546	obskind = opts.obskind
	1547	if obskind == 'single-station':
	1548	if opts.stloc is None:
	1549	print errormsg
	1550	print ' ' + main + ': No station location provided !!'
	1551	quit(-1)
	1552	else:
[505]	1553	stationdesc = [opts.stloc.split(',')[0].replace('\|',' '), \
[492]	1554	np.float(opts.stloc.split(',')[1]), np.float(opts.stloc.split(',')[2]),\
	1555	np.float(opts.stloc.split(',')[3])]
[330]	1556
	1557	if opts.fobs is None:
	1558	print errormsg
	1559	print ' ' + main + ': No observations file is provided!!'
	1560	quit(-1)
	1561	else:
	1562	if not os.path.isfile(opts.fobs):
	1563	print errormsg
	1564	print ' ' + main + ": observations file '" + opts.fobs + "' does not exist !!"
	1565	quit(-1)
	1566
	1567	if opts.fsim is None:
	1568	print errormsg
	1569	print ' ' + main + ': No simulation file is provided!!'
	1570	quit(-1)
	1571	else:
	1572	if not os.path.isfile(opts.fsim):
	1573	print errormsg
	1574	print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
	1575	quit(-1)
	1576
	1577	if opts.vars is None:
	1578	print errormsg
	1579	print ' ' + main + ': No list of couples of variables is provided!!'
	1580	print ' a ',' list of values [varsim]@[varobs],... is needed'
	1581	quit(-1)
	1582	else:
	1583	valvars = []
[333]	1584	vs = opts.vars.split(',')
[330]	1585	for v in vs:
	1586	vsecs = v.split('@')
[333]	1587	if len(vsecs) < 2:
[330]	1588	print errormsg
	1589	print ' ' + main + ': wrong number of values in:',vsecs, \
	1590	' at least 2 are needed !!'
	1591	print ' [varsim]@[varobs]@[[oper][val]]'
	1592	quit(-1)
[333]	1593	if len(vsecs) > 2:
	1594	if not searchInlist(simopers,vsecs[2]):
	1595	print errormsg
	1596	print main + ": operation on simulation values '" + vsecs[2] + \
	1597	"' not ready !!"
	1598	quit(-1)
	1599
[330]	1600	valvars.append(vsecs)
	1601
	1602	# Openning observations trajectory
	1603	##
	1604	oobs = NetCDFFile(opts.fobs, 'r')
	1605
	1606	valdimobs = {}
	1607	for dn in dims:
	1608	print dn,':',dims[dn]
	1609	if dims[dn][1] != 'None':
	1610	if not oobs.dimensions.has_key(dims[dn][1]):
	1611	print errormsg
	1612	print ' ' + main + ": observations file does not have dimension '" + \
	1613	dims[dn][1] + "' !!"
	1614	quit(-1)
	1615	if vardims[dn][1] != 'None':
	1616	if not oobs.variables.has_key(vardims[dn][1]):
	1617	print errormsg
	1618	print ' ' + main + ": observations file does not have varibale " + \
	1619	"dimension '" + vardims[dn][1] + "' !!"
	1620	quit(-1)
	1621	valdimobs[dn] = oobs.variables[vardims[dn][1]][:]
	1622	else:
	1623	if dn == 'X':
[492]	1624	valdimobs[dn] = stationdesc[1]
[330]	1625	elif dn == 'Y':
[492]	1626	valdimobs[dn] = stationdesc[2]
[330]	1627	elif dn == 'Z':
[492]	1628	valdimobs[dn] = stationdesc[3]
[330]	1629
	1630	osim = NetCDFFile(opts.fsim, 'r')
	1631
	1632	valdimsim = {}
	1633	for dn in dims:
[464]	1634	if dims[dn][0] != 'None':
	1635	if not osim.dimensions.has_key(dims[dn][0]):
	1636	print errormsg
	1637	print ' ' + main + ": simulation file '" + opts.fsim + \
	1638	"' does not have dimension '" + dims[dn][0] + "' !!"
	1639	print ' it has: ',osim.dimensions
	1640	quit(-1)
[343]	1641
[464]	1642	if not osim.variables.has_key(vardims[dn][0]) and \
	1643	not searchInlist(varNOcheck,vardims[dn][0]):
	1644	print errormsg
	1645	print ' ' + main + ": simulation file '" + opts.fsim + \
	1646	"' does not have varibale dimension '" + vardims[dn][0] + "' !!"
	1647	print ' it has variables:',osim.variables
	1648	quit(-1)
	1649	if searchInlist(varNOcheck,vardims[dn][0]):
	1650	valdimsim[dn] = compute_varNOcheck(osim, vardims[dn][0])
[465]	1651	else:
[464]	1652	valdimsim[dn] = osim.variables[vardims[dn][0]][:]
[330]	1653
	1654	# General characteristics
[343]	1655	dimtobs = valdimobs['T'].shape[0]
	1656	dimtsim = valdimsim['T'].shape[0]
[330]	1657
	1658	print main +': observational time-steps:',dimtobs,'simulation:',dimtsim
	1659
[337]	1660	notfound = np.zeros((dimtobs), dtype=int)
	1661
[330]	1662	if obskind == 'multi-points':
	1663	trajpos = np.zeros((2,dimt),dtype=int)
[337]	1664	for it in range(dimtobs):
[330]	1665	trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'], \
	1666	[valdimobs['X'][it],valdimobss['Y'][it]])
	1667	elif obskind == 'single-station':
[498]	1668	trajpos = None
[333]	1669	stationpos = np.zeros((2), dtype=int)
[466]	1670	if valdimsim.has_key('X') and valdimsim.has_key('Y'):
	1671	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],[valdimobs['Y'], \
	1672	valdimobs['X']])
	1673	iid = 0
	1674	for idn in osim.variables[vardims['X'][0]].dimensions:
	1675	if idn == dims['X'][0]:
	1676	stationpos[1] = stsimpos[iid]
	1677	elif idn == dims['Y'][0]:
	1678	stationpos[0] = stsimpos[iid]
[333]	1679
[466]	1680	iid = iid + 1
	1681	print main + ': station point in simulation:', stationpos
	1682	print ' station position:',valdimobs['X'],',',valdimobs['Y']
	1683	print ' simulation coord.:',valdimsim['X'][tuple(stsimpos)],',', \
	1684	valdimsim['Y'][tuple(stsimpos)]
	1685	else:
	1686	print main + ': validation with two time-series !!'
[333]	1687
[330]	1688	elif obskind == 'trajectory':
[337]	1689	if opts.trajf is not None:
	1690	if not os.path.isfile(opts.fsim):
	1691	print errormsg
	1692	print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
	1693	quit(-1)
	1694	else:
	1695	otrjf = NetCDFFile(opts.fsim, 'r')
	1696	trajpos[0,:] = otrjf.variables['obssimtrj'][0]
	1697	trajpos[1,:] = otrjf.variables['obssimtrj'][1]
	1698	otrjf.close()
[330]	1699	else:
[337]	1700	if dims.has_key('Z'):
	1701	trajpos = np.zeros((3,dimtobs),dtype=int)
	1702	for it in range(dimtobs):
	1703	if np.mod(it*100./dimtobs,10.) == 0.:
	1704	print ' trajectory done: ',it*100./dimtobs,'%'
	1705	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \
	1706	[valdimobs['Y'][it],valdimobs['X'][it]])
	1707	stationpos = np.zeros((2), dtype=int)
	1708	iid = 0
	1709	for idn in osim.variables[vardims['X'][0]].dimensions:
	1710	if idn == dims['X'][0]:
	1711	stationpos[1] = stsimpos[iid]
	1712	elif idn == dims['Y'][0]:
	1713	stationpos[0] = stsimpos[iid]
	1714	iid = iid + 1
	1715	if stationpos[0] == 0 and stationpos[1] == 0: notfound[it] = 1
	1716
	1717	trajpos[0,it] = stationpos[0]
	1718	trajpos[1,it] = stationpos[1]
	1719	# In the simulation 'Z' varies with time ... non-hydrostatic model! ;)
	1720	# trajpos[2,it] = index_mat(valdimsim['Z'][it,:,stationpos[0], \
	1721	# stationpos[1]], valdimobs['Z'][it])
	1722	else:
	1723	trajpos = np.zeros((2,dimtobs),dtype=int)
	1724	for it in range(dimtobs):
	1725	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \
	1726	[valdimobs['Y'][it],valdimobss['X'][it]])
	1727	stationpos = np.zeros((2), dtype=int)
	1728	iid = 0
	1729	for idn in osim.variables[vardims['X'][0]].dimensions:
	1730	if idn == dims['X'][0]:
	1731	stationpos[1] = stsimpos[iid]
	1732	elif idn == dims['Y'][0]:
	1733	stationpos[0] = stsimpos[iid]
	1734	iid = iid + 1
[338]	1735	if stationpos[0] == 0 or stationpos[1] == 0: notfound[it] = 1
[330]	1736
[337]	1737	trajpos[0,it] = stationspos[0]
	1738	trajpos[1,it] = stationspos[1]
	1739
	1740	print main + ': not found',np.sum(notfound),'points of the trajectory'
	1741
[330]	1742	# Getting times
	1743	tobj = oobs.variables[vardims['T'][1]]
	1744	obstunits = tobj.getncattr('units')
[343]	1745	if vardims['T'][0] == 'WRFT':
	1746	tsim = valdimsim['T'][:]
	1747	simtunits = 'seconds since 1949-12-01 00:00:00'
	1748	else:
[467]	1749	tsim = osim.variables[vardims['T'][0]][:]
	1750	otsim = osim.variables[vardims['T'][0]]
	1751	simtunits = otsim.getncattr('units')
[330]	1752
[467]	1753	simobstimes = coincident_CFtimes(tsim, obstunits, simtunits)
[330]	1754
[493]	1755	#
	1756	## Looking for exact/near times
	1757	#
	1758
	1759	# Exact Coincident times
[333]	1760	##
[493]	1761	exacttvalues0 = []
	1762	for it in range(dimtsim):
	1763	ot = 0
	1764	for ito in range(ot,dimtobs-1):
[494]	1765	if valdimobs['T'][ito] == simobstimes[it]:
[493]	1766	ot = ito
	1767	exacttvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito]])
	1768
	1769	exacttvalues = np.array(exacttvalues0, dtype=np.float)
	1770	Nexactt = len(exacttvalues[:,0])
	1771
	1772	print main + ': found',Nexactt,'Temporal same values in simulation and observations'
[511]	1773
[493]	1774	# Sim Coincident times
	1775	##
[511]	1776	sumsimval = 0.
	1777	sum2simval = 0.
	1778	Nsimt = 0
[337]	1779	coindtvalues0 = []
[511]	1780	coindtvalues0st = []
	1781	tsiminit = 0
	1782	tsimend = 0
	1783
[333]	1784	for it in range(dimtsim):
	1785	ot = 0
	1786	for ito in range(ot,dimtobs-1):
	1787	if valdimobs['T'][ito] < simobstimes[it] and valdimobs['T'][ito+1] > \
	1788	simobstimes[it]:
	1789	ot = ito
	1790	tdist = simobstimes[it] - valdimobs['T'][ito]
[511]	1791	coindtvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito], \
[337]	1792	tdist])
[511]	1793	Nsimt = Nsimt + 1
	1794	if tsiminit == 0: tsiminit=simobstimes[it]
	1795	tsimend = simobstimes[it]
	1796	elif simobstimes[it] > valdimobs['T'][ito+1]:
	1797	coindtvalues0st.append([Nsimt, ito, valdimobs['T'][ito],tsimend-tsiminit])
[330]	1798
[337]	1799	coindtvalues = np.array(coindtvalues0, dtype=np.float)
[511]	1800	coindtvaluesst = np.array(coindtvalues0st, dtype=np.float)
[337]	1801
	1802	Ncoindt = len(coindtvalues[:,0])
[493]	1803	print main + ': found',Ncoindt,'Simulation time-interval (within consecutive ' + \
	1804	'observed times) coincident times between simulation and observations'
[333]	1805
[337]	1806	if Ncoindt == 0:
	1807	print warnmsg
	1808	print main + ': no coincident times found !!'
	1809	print ' stopping it'
	1810	quit(-1)
	1811
[333]	1812	# Validating
	1813	##
	1814
	1815	onewnc = NetCDFFile(ofile, 'w')
	1816
	1817	# Dimensions
	1818	newdim = onewnc.createDimension('time',None)
[511]	1819	newdim = onewnc.createDimension('betweentime',None)
[346]	1820	newdim = onewnc.createDimension('bnds',2)
[337]	1821	newdim = onewnc.createDimension('obstime',None)
[333]	1822	newdim = onewnc.createDimension('couple',2)
	1823	newdim = onewnc.createDimension('StrLength',StringLength)
[337]	1824	newdim = onewnc.createDimension('xaround',Ngrid*2+1)
	1825	newdim = onewnc.createDimension('yaround',Ngrid*2+1)
[503]	1826	newdim = onewnc.createDimension('gstats',13)
[337]	1827	newdim = onewnc.createDimension('stats',5)
[346]	1828	newdim = onewnc.createDimension('tstats',6)
[503]	1829	newdim = onewnc.createDimension('Nstsim', 2)
[333]	1830
	1831	# Variable dimensions
	1832	##
[511]	1833	newvar = onewnc.createVariable('exacttime','f8',('time'))
	1834	basicvardef(newvar, 'obstime', 'exact coincident time observations and simulation', obstunits)
[333]	1835	set_attribute(newvar, 'calendar', 'standard')
[512]	1836	print 'Lluis Nexacttvalues:',len(exacttvalues[:,3])
[511]	1837	newvar[:] = exacttvalues[:,3]
	1838
	1839	newvar = onewnc.createVariable('obstime','f8',('obstime'))
	1840	basicvardef(newvar, 'obstime', 'time observations for between values', obstunits)
	1841	set_attribute(newvar, 'calendar', 'standard')
[346]	1842	set_attribute(newvar, 'bounds', 'time_bnds')
[337]	1843	newvar[:] = coindtvalues[:,3]
[333]	1844
[511]	1845	newvar = onewnc.createVariable('betweentime','f8',('betweentime'))
	1846	basicvardef(newvar, 'obstime', 'time simulations for between values', simtunits )
	1847	set_attribute(newvar, 'calendar', 'standard')
	1848	set_attribute(newvar, 'bounds', 'time_bnds')
	1849	newvar[:] = coindtvalues[:,2]
	1850
[333]	1851	newvar = onewnc.createVariable('couple', 'c', ('couple','StrLength'))
	1852	basicvardef(newvar, 'couple', 'couples of values', '-')
	1853	writing_str_nc(newvar, ['sim','obs'], StringLength)
	1854
[337]	1855	newvar = onewnc.createVariable('statistics', 'c', ('stats','StrLength'))
	1856	basicvardef(newvar, 'statistics', 'statitics from values', '-')
	1857	writing_str_nc(newvar, statsn, StringLength)
	1858
[344]	1859	newvar = onewnc.createVariable('gstatistics', 'c', ('gstats','StrLength'))
	1860	basicvardef(newvar, 'gstatistics', 'global statitics from values', '-')
	1861	writing_str_nc(newvar, gstatsn, StringLength)
	1862
	1863	newvar = onewnc.createVariable('tstatistics', 'c', ('tstats','StrLength'))
	1864	basicvardef(newvar, 'tstatistics', 'statitics from values along time', '-')
	1865	writing_str_nc(newvar, ostatsn, StringLength)
	1866
[502]	1867	newvar = onewnc.createVariable('ksimstatistics', 'c', ('Nstsim','StrLength'))
	1868	basicvardef(newvar, 'ksimstatistics', 'kind of simulated statitics', '-')
	1869	writing_str_nc(newvar, Lstdescsim, StringLength)
	1870
	1871
[337]	1872	if obskind == 'trajectory':
	1873	if dims.has_key('Z'):
	1874	newdim = onewnc.createDimension('trj',3)
	1875	else:
	1876	newdim = onewnc.createDimension('trj',2)
	1877
	1878	newvar = onewnc.createVariable('obssimtrj','i',('obstime','trj'))
	1879	basicvardef(newvar, 'obssimtrj', 'trajectory on the simulation grid', '-')
	1880	newvar[:] = trajpos.transpose()
	1881
	1882	if dims.has_key('Z'):
	1883	newdim = onewnc.createDimension('simtrj',4)
	1884	else:
	1885	newdim = onewnc.createDimension('simtrj',3)
	1886
[333]	1887	Nvars = len(valvars)
	1888	for ivar in range(Nvars):
	1889	simobsvalues = []
	1890
	1891	varsimobs = valvars[ivar][0] + '_' + valvars[ivar][1]
[340]	1892	print ' ' + varsimobs + '... .. .'
[333]	1893
	1894	if not oobs.variables.has_key(valvars[ivar][1]):
	1895	print errormsg
	1896	print ' ' + main + ": observations file has not '" + valvars[ivar][1] + \
	1897	"' !!"
	1898	quit(-1)
[337]	1899
[333]	1900	if not osim.variables.has_key(valvars[ivar][0]):
[337]	1901	if not searchInlist(varNOcheck, valvars[ivar][0]):
	1902	print errormsg
	1903	print ' ' + main + ": simulation file has not '" + valvars[ivar][0] + \
	1904	"' !!"
	1905	quit(-1)
	1906	else:
	1907	ovsim = compute_varNOcheck(osim, valvars[ivar][0])
	1908	else:
	1909	ovsim = osim.variables[valvars[ivar][0]]
[333]	1910
[340]	1911	for idn in ovsim.dimensions:
	1912	if not searchInlist(simdims.values(),idn):
	1913	print errormsg
	1914	print ' ' + main + ": dimension '" + idn + "' of variable '" + \
	1915	valvars[ivar][0] + "' not provided as reference coordinate [X,Y,Z,T] !!"
	1916	quit(-1)
	1917
[333]	1918	ovobs = oobs.variables[valvars[ivar][1]]
[468]	1919	if searchInlist(ovobs.ncattrs(),'_FillValue'):
	1920	oFillValue = ovobs.getncattr('_FillValue')
[337]	1921	else:
[498]	1922	oFillValue = None
[333]	1923
[502]	1924	# Observed values temporally exact times
[520]	1925	Esimobsvalues, EsimobsSvalues, EsimobsTtvalues, trjsim = \
[502]	1926	getting_ValidationValues(obskind, Nexactt, dims, trajpos, ovsim, ovobs, \
[514]	1927	exacttvalues, oFillValue, Ngrid, 'instantaneous')
[501]	1928
	1929	# Observed values temporally around coincident times
[520]	1930	simobsvalues, simobsSvalues, simobsTtvalues, trjsim = \
[498]	1931	getting_ValidationValues(obskind, Ncoindt, dims, trajpos, ovsim, ovobs, \
[515]	1932	coindtvalues, oFillValue, Ngrid, 'tbackwardSmean')
[344]	1933
[498]	1934	# Re-arranging values...
[502]	1935	Earrayvals = np.array(Esimobsvalues)
[339]	1936	arrayvals = np.array(simobsvalues)
	1937	if len(valvars[ivar]) > 2:
	1938	const=np.float(valvars[ivar][3])
	1939	if valvars[ivar][2] == 'sumc':
[502]	1940	EsimobsSvalues = EsimobsSvalues + const
	1941	Earrayvals[:,0] = Earrayvals[:,0] + const
[339]	1942	simobsSvalues = simobsSvalues + const
	1943	arrayvals[:,0] = arrayvals[:,0] + const
	1944	elif valvars[ivar][2] == 'subc':
[502]	1945	EsimobsSvalues = EsimobsSvalues - const
	1946	Earrayvals[:,0] = Earrayvals[:,0] - const
[339]	1947	simobsSvalues = simobsSvalues - const
	1948	arrayvals[:,0] = arrayvals[:,0] - const
	1949	elif valvars[ivar][2] == 'mulc':
[502]	1950	EsimobsSvalues = EsimobsSvalues * const
	1951	Earrayvals[:,0] = Earrayvals[:,0] * const
[339]	1952	simobsSvalues = simobsSvalues * const
	1953	arrayvals[:,0] = arrayvals[:,0] * const
	1954	elif valvars[ivar][2] == 'divc':
[502]	1955	EsimobsSvalues = EsimobsSvalues / const
	1956	Earrayvals[:,0] = Earrayvals[:,0] / const
[339]	1957	simobsSvalues = simobsSvalues / const
	1958	arrayvals[:,0] = arrayvals[:,0] / const
	1959	else:
	1960	print errormsg
	1961	print ' ' + fname + ": operation '" + valvars[ivar][2] + "' not ready!!"
	1962	quit(-1)
	1963
[343]	1964	# statisics sim
[502]	1965	simstats = np.zeros((2,5), dtype=np.float)
	1966	simstats[0,0] = np.min(Earrayvals[:,0])
	1967	simstats[0,1] = np.max(Earrayvals[:,0])
	1968	simstats[0,2] = np.mean(Earrayvals[:,0])
	1969	simstats[0,3] = np.mean(Earrayvals[:,0]*Earrayvals[:,0])
	1970	simstats[0,4] = np.sqrt(simstats[0,3] - simstats[0,2]*simstats[0,2])
[343]	1971
[502]	1972	simstats[1,0] = np.min(arrayvals[:,0])
	1973	simstats[1,1] = np.max(arrayvals[:,0])
	1974	simstats[1,2] = np.mean(arrayvals[:,0])
	1975	simstats[1,3] = np.mean(arrayvals[:,0]*arrayvals[:,0])
	1976	simstats[1,4] = np.sqrt(simstats[1,3] - simstats[1,2]*simstats[1,2])
	1977
[343]	1978	# statisics obs
[501]	1979	# Masking 'nan'
[503]	1980	Eobsmask0 = np.where(Earrayvals[:,1] != Earrayvals[:,1], fillValueF,
	1981	Earrayvals[:,1])
[501]	1982	obsmask0 = np.where(arrayvals[:,1] != arrayvals[:,1], fillValueF, arrayvals[:,1])
[503]	1983
	1984	Eobsmask = ma.masked_equal(Eobsmask0, fillValueF)
	1985	Eobsmask2 = Eobsmask*Eobsmask
[501]	1986	obsmask = ma.masked_equal(obsmask0, fillValueF)
[343]	1987	obsmask2 = obsmask*obsmask
	1988
[503]	1989	obsstats = np.zeros((2,5), dtype=np.float)
	1990	obsstats[0,0] = obsmask.min()
	1991	obsstats[0,1] = obsmask.max()
	1992	obsstats[0,2] = obsmask.mean()
	1993	obsstats[0,3] = obsmask2.mean()
	1994	obsstats[0,4] = np.sqrt(obsstats[0,3] - obsstats[0,2]*obsstats[0,2])
	1995
	1996	obsstats[1,0] = obsmask.min()
	1997	obsstats[1,1] = obsmask.max()
	1998	obsstats[1,2] = obsmask.mean()
	1999	obsstats[1,3] = obsmask2.mean()
	2000	obsstats[1,4] = np.sqrt(obsstats[1,3] - obsstats[1,2]*obsstats[1,2])
[501]	2001
[343]	2002	# Statistics sim-obs
[502]	2003	simobsstats = np.zeros((2,13), dtype=np.float)
	2004	Ediffvals = np.zeros((Nexactt), dtype=np.float)
[343]	2005	diffvals = np.zeros((Ncoindt), dtype=np.float)
	2006
[503]	2007	Ediffvals = Earrayvals[:,0] - Eobsmask
[344]	2008	diffvals = arrayvals[:,0] - obsmask
	2009
[502]	2010	Ediff2vals = Ediffvals*Ediffvals
	2011	Esumdiff = Ediffvals.sum()
	2012	Esumdiff2 = Ediff2vals.sum()
[343]	2013
[502]	2014	diff2vals = diffvals*diffvals
	2015	sumdiff = diffvals.sum()
	2016	sumdiff2 = diff2vals.sum()
	2017
[503]	2018	simobsstats[0,0] = simstats[0,0] - obsstats[0,0]
	2019	simobsstats[0,1] = np.mean(Earrayvals[:,0]*Eobsmask)
[502]	2020	simobsstats[0,2] = Ediffvals.min()
	2021	simobsstats[0,3] = Ediffvals.max()
	2022	simobsstats[0,4] = Ediffvals.mean()
	2023	simobsstats[0,5] = np.abs(Ediffvals).mean()
	2024	simobsstats[0,6] = np.sqrt(Ediff2vals.mean())
	2025	simobsstats[0,7], simobsstats[0,8] = sts.mstats.pearsonr(Earrayvals[:,0], \
[503]	2026	Eobsmask)
[502]	2027
[503]	2028	simobsstats[1,0] = simstats[1,0] - obsstats[1,0]
[502]	2029	simobsstats[1,1] = np.mean(arrayvals[:,0]*obsmask)
	2030	simobsstats[1,2] = diffvals.min()
	2031	simobsstats[1,3] = diffvals.max()
	2032	simobsstats[1,4] = diffvals.mean()
	2033	simobsstats[1,5] = np.abs(diffvals).mean()
	2034	simobsstats[1,6] = np.sqrt(diff2vals.mean())
	2035	simobsstats[1,7], simobsstats[1,8] = sts.mstats.pearsonr(arrayvals[:,0], \
	2036	obsmask)
	2037	# From:
[503]	2038	#Willmott, C. J. 1981. 'On the validation of models. Physical Geography', 2, 184-194
	2039	#Willmott, C. J. (1984). 'On the evaluation of model performance in physical
	2040	# geography'. Spatial Statistics and Models, G. L. Gaile and C. J. Willmott, eds.,
	2041	# 443-460
	2042	#Willmott, C. J., S. G. Ackleson, R. E. Davis, J. J. Feddema, K. M. Klink, D. R.
	2043	# Legates, J. O'Donnell, and C. M. Rowe (1985), 'Statistics for the Evaluation and
	2044	# Comparison of Models', J. Geophys. Res., 90(C5), 8995-9005
	2045	#Legates, D. R., and G. J. McCabe Jr. (1999), 'Evaluating the Use of "Goodness-of-Fit"
	2046	# Measures in Hydrologic and Hydroclimatic Model Validation', Water Resour. Res.,
	2047	# 35(1), 233-241
[502]	2048	#
	2049	# Deviation of residuals (SDR)
	2050	simobsstats[0,9] = np.mean(np.sqrt(np.abs((Ediffvals-simobsstats[0,0])*(Ediffvals-\
[503]	2051	obsstats[0,0]))))
[502]	2052	simobsstats[1,9] = np.mean(np.sqrt(np.abs((diffvals-simobsstats[1,0])*(diffvals- \
[503]	2053	obsstats[1,0]))))
[502]	2054	# Index of Efficiency (IOE)
[503]	2055	Eobsbias2series = (Eobsmask - obsstats[0,0])*(Eobsmask - obsstats[0,0])
	2056	Esumobsbias2series = Eobsbias2series.sum()
	2057	obsbias2series = (obsmask - obsstats[1,0])*(obsmask - obsstats[1,0])
[502]	2058	sumobsbias2series = obsbias2series.sum()
	2059
[503]	2060	simobsstats[0,10] = 1. - Esumdiff2/(Esumobsbias2series)
[502]	2061	simobsstats[1,10] = 1. - sumdiff2/(sumobsbias2series)
	2062	# Index of Agreement (IOA)
[503]	2063	Esimbias2series = Earrayvals[:,0] - obsstats[0,0]
	2064	simbias2series = arrayvals[:,0] - obsstats[1,0]
[502]	2065
[503]	2066	Eobsbias2series = Eobsmask - obsstats[0,0]
	2067	obsbias2series = obsmask - obsstats[1,0]
[502]	2068
	2069	Eabssimbias2series = np.abs(Esimbias2series)
	2070	abssimbias2series = np.abs(simbias2series)
[503]	2071	Eabsobsbias2series = np.where(Eobsbias2series<0, -Eobsbias2series, Eobsbias2series)
	2072	absobsbias2series = np.where(obsbias2series<0, -obsbias2series, obsbias2series)
[502]	2073
[503]	2074	Eabssimobsbias2series=(Eabssimbias2series+Eabsobsbias2series)*(Eabssimbias2series+\
	2075	Eabsobsbias2series)
[502]	2076	abssimobsbias2series = (abssimbias2series+absobsbias2series)*(abssimbias2series +\
	2077	absobsbias2series)
	2078
	2079	simobsstats[0,11] = 1. - Esumdiff2/(Eabssimobsbias2series.sum())
	2080	simobsstats[1,11] = 1. - sumdiff2/(abssimobsbias2series.sum())
	2081	# Fractional Mean Bias (FMB)
[503]	2082	simobsstats[0,12]=(simstats[0,0]-obsstats[0,0])/(0.5*(simstats[0,0]+obsstats[0,0]))
	2083	simobsstats[1,12]=(simstats[1,0]-obsstats[1,0])/(0.5*(simstats[1,0]+obsstats[1,0]))
[502]	2084
[514]	2085	# Statistics around sim values exact
	2086	Earoundstats = np.zeros((5,Nexactt), dtype=np.float)
	2087	for it in range(Nexactt):
	2088	Earoundstats[0,it] = np.min(EsimobsSvalues[it,])
	2089	Earoundstats[1,it] = np.max(EsimobsSvalues[it,])
	2090	Earoundstats[2,it] = np.mean(EsimobsSvalues[it,])
	2091	Earoundstats[3,it] = np.mean(EsimobsSvalues[it,]*EsimobsSvalues[it,])
	2092	Earoundstats[4,it] = np.sqrt(Earoundstats[3,it] - Earoundstats[2,it]* \
	2093	Earoundstats[2,it])
	2094
	2095	# Statistics around sim values between
[337]	2096	aroundstats = np.zeros((5,Ncoindt), dtype=np.float)
	2097	for it in range(Ncoindt):
	2098	aroundstats[0,it] = np.min(simobsSvalues[it,])
	2099	aroundstats[1,it] = np.max(simobsSvalues[it,])
	2100	aroundstats[2,it] = np.mean(simobsSvalues[it,])
	2101	aroundstats[3,it] = np.mean(simobsSvalues[it,]*simobsSvalues[it,])
	2102	aroundstats[4,it] = np.sqrt(aroundstats[3,it] - aroundstats[2,it]* \
	2103	aroundstats[2,it])
	2104
[511]	2105	# exact sim Values to netCDF
	2106	newvar = onewnc.createVariable(valvars[ivar][0] + '_Esim', 'f', ('time'), \
[337]	2107	fill_value=fillValueF)
[511]	2108	descvar = 'exact time simulated: ' + valvars[ivar][0]
[346]	2109	basicvardef(newvar, valvars[ivar][0], descvar, ovobs.getncattr('units'))
[512]	2110	print 'Lluis NEarratvals:',len(Earrayvals[:,0])
[511]	2111	newvar[:] = Earrayvals[:,0]
	2112
	2113	# exact obs Values to netCDF
	2114	newvar = onewnc.createVariable(valvars[ivar][1] + '_Eobs', 'f', ('time'), \
	2115	fill_value=fillValueF)
	2116	descvar = 'exact time observed: ' + valvars[ivar][1]
	2117	basicvardef(newvar, valvars[ivar][1], descvar, ovobs.getncattr('units'))
	2118	newvar[:] = Earrayvals[:,1]
	2119
	2120	# between sim Values to netCDF
	2121	newvar = onewnc.createVariable(valvars[ivar][0] + '_sim', 'f', ('betweentime'), \
	2122	fill_value=fillValueF)
	2123	descvar = 'between observed time simulated: ' + valvars[ivar][0]
	2124	basicvardef(newvar, valvars[ivar][0], descvar, ovobs.getncattr('units'))
[346]	2125	newvar[:] = arrayvals[:,0]
[333]	2126
[511]	2127	# between obs Values to netCDF
	2128	newvar = onewnc.createVariable(valvars[ivar][1] + '_obs', 'f', ('obstime'), \
[346]	2129	fill_value=fillValueF)
	2130	descvar = 'observed: ' + valvars[ivar][1]
	2131	basicvardef(newvar, valvars[ivar][1], descvar, ovobs.getncattr('units'))
	2132	newvar[:] = arrayvals[:,1]
	2133
[514]	2134	# Around values exact
	2135	if not onewnc.variables.has_key(valvars[ivar][0] + 'Earound'):
	2136	if dims.has_key('Z'):
	2137	if not onewnc.dimensions.has_key('zaround'):
	2138	newdim = onewnc.createDimension('zaround',Ngrid*2+1)
	2139	newvar = onewnc.createVariable(valvars[ivar][0] + 'Earound', 'f', \
	2140	('time','zaround','yaround','xaround'), fill_value=fillValueF)
	2141	else:
	2142	newvar = onewnc.createVariable(valvars[ivar][0] + 'Earound', 'f', \
	2143	('time','yaround','xaround'), fill_value=fillValueF)
	2144
	2145	descvar = 'exact around simulated values +/- grid values: ' + valvars[ivar][0]
	2146	basicvardef(newvar, varsimobs + 'Earound', descvar, ovobs.getncattr('units'))
	2147	newvar[:] = EsimobsSvalues
	2148
	2149	# Around values between
[341]	2150	if not onewnc.variables.has_key(valvars[ivar][0] + 'around'):
	2151	if dims.has_key('Z'):
	2152	if not onewnc.dimensions.has_key('zaround'):
	2153	newdim = onewnc.createDimension('zaround',Ngrid*2+1)
	2154	newvar = onewnc.createVariable(valvars[ivar][0] + 'around', 'f', \
[514]	2155	('betweentime','zaround','yaround','xaround'),fill_value=fillValueF)
[341]	2156	else:
	2157	newvar = onewnc.createVariable(valvars[ivar][0] + 'around', 'f', \
[512]	2158	('betweentime','yaround','xaround'), fill_value=fillValueF)
[339]	2159
[341]	2160	descvar = 'around simulated values +/- grid values: ' + valvars[ivar][0]
	2161	basicvardef(newvar, varsimobs + 'around', descvar, ovobs.getncattr('units'))
	2162	newvar[:] = simobsSvalues
[337]	2163
[344]	2164	# sim Statistics
	2165	if not searchInlist(onewnc.variables,valvars[ivar][0] + 'stsim'):
[502]	2166	newvar = onewnc.createVariable(valvars[ivar][0] + 'stsim', 'f', ('Nstsim', \
	2167	'stats'), fill_value=fillValueF)
[344]	2168	descvar = 'simulated statisitcs: ' + valvars[ivar][0]
[512]	2169	basicvardef(newvar, valvars[ivar][0] + 'stsim', descvar, \
	2170	ovobs.getncattr('units'))
[344]	2171	newvar[:] = simstats
	2172
	2173	# obs Statistics
	2174	if not searchInlist(onewnc.variables,valvars[ivar][1] + 'stobs'):
[503]	2175	newvar = onewnc.createVariable(valvars[ivar][1] + 'stobs', 'f', ('Nstsim', \
	2176	'stats'), fill_value=fillValueF)
[344]	2177	descvar = 'observed statisitcs: ' + valvars[ivar][1]
	2178	basicvardef(newvar, valvars[ivar][1] + 'stobs', descvar, \
	2179	ovobs.getncattr('units'))
	2180	newvar[:] = obsstats
	2181
	2182	# sim-obs Statistics
	2183	if not searchInlist(onewnc.variables,varsimobs + 'st'):
[502]	2184	newvar = onewnc.createVariable(varsimobs + 'st', 'f', ('Nstsim', 'gstats'), \
[344]	2185	fill_value=fillValueF)
[512]	2186	descvar = 'simulated-observed tatisitcs: ' + varsimobs
[344]	2187	basicvardef(newvar, varsimobs + 'st', descvar, ovobs.getncattr('units'))
	2188	newvar[:] = simobsstats
	2189
[514]	2190	# around sim Statistics exact
	2191	if not searchInlist(onewnc.variables,valvars[ivar][0] + 'Estaround'):
	2192	newvar = onewnc.createVariable(valvars[ivar][0] + 'Estaround', 'f', \
	2193	('time','stats'), fill_value=fillValueF)
	2194	descvar = 'exact around (' + str(Ngrid) + ', ' + str(Ngrid) + \
	2195	') simulated statisitcs: ' + valvars[ivar][0]
	2196	basicvardef(newvar, valvars[ivar][0] + 'Estaround', descvar, \
	2197	ovobs.getncattr('units'))
	2198	set_attribute(newvar, 'cell_methods', 'Etime_bnds')
	2199	newvar[:] = Earoundstats.transpose()
	2200
	2201	if not searchInlist(onewnc.variables, 'Etime_bnds'):
	2202	newvar = onewnc.createVariable('Etime_bnds','f8',('time','bnds'))
	2203	basicvardef(newvar, 'Etime_bnds', 'betweentime', obstunits )
	2204	set_attribute(newvar, 'calendar', 'standard')
	2205	newvar[:] = EsimobsTtvalues
	2206
	2207	# around obs Statistics exact
	2208	if not searchInlist(onewnc.variables,valvars[ivar][1] + 'Estaround'):
	2209	newvar = onewnc.createVariable(valvars[ivar][1] + 'Estaround', 'f', \
	2210	('time','tstats'), fill_value=fillValueF)
	2211	descvar = 'exact around temporal observed statisitcs: ' + valvars[ivar][1]
	2212	basicvardef(newvar, valvars[ivar][1] + 'Estaround', descvar, \
	2213	ovobs.getncattr('units'))
	2214	set_attribute(newvar, 'cell_methods', 'Etime_bnds')
	2215
	2216	newvar[:] = aroundostats.transpose()
	2217
	2218	# around sim Statistics between
[344]	2219	if not searchInlist(onewnc.variables,valvars[ivar][0] + 'staround'):
[341]	2220	newvar = onewnc.createVariable(valvars[ivar][0] + 'staround', 'f', \
[512]	2221	('betweentime','stats'), fill_value=fillValueF)
[514]	2222	descvar = 'between around (' + str(Ngrid) + ', ' + str(Ngrid) + \
[344]	2223	') simulated statisitcs: ' + valvars[ivar][0]
	2224	basicvardef(newvar, valvars[ivar][0] + 'staround', descvar, \
	2225	ovobs.getncattr('units'))
[514]	2226	set_attribute(newvar, 'cell_methods', 'time_bnds')
[341]	2227	newvar[:] = aroundstats.transpose()
[337]	2228
[346]	2229	if not searchInlist(onewnc.variables, 'time_bnds'):
[512]	2230	newvar = onewnc.createVariable('time_bnds','f8',('betweentime','bnds'))
	2231	basicvardef(newvar, 'time_bnds', 'betweentime', obstunits )
[346]	2232	set_attribute(newvar, 'calendar', 'standard')
	2233	newvar[:] = simobsTtvalues
	2234
[514]	2235	# around obs Statistics between
[344]	2236	if not searchInlist(onewnc.variables,valvars[ivar][1] + 'staround'):
	2237	newvar = onewnc.createVariable(valvars[ivar][1] + 'staround', 'f', \
[512]	2238	('betweentime','tstats'), fill_value=fillValueF)
[514]	2239	descvar = 'betweem around temporal observed statisitcs: ' + valvars[ivar][1]
[344]	2240	basicvardef(newvar, valvars[ivar][1] + 'staround', descvar, \
	2241	ovobs.getncattr('units'))
[514]	2242	set_attribute(newvar, 'cell_methods', 'time_bnds')
[344]	2243	newvar[:] = aroundostats.transpose()
	2244
[341]	2245	onewnc.sync()
[333]	2246
[500]	2247	if trjsim is not None:
[512]	2248	newvar = onewnc.createVariable('simtrj','i',('betweentime','simtrj'))
[500]	2249	basicvardef(newvar,'simtrj','coordinates [X,Y,Z,T] of the coincident ' + \
	2250	'trajectory in sim', obstunits)
	2251	newvar[:] = trjsim.transpose()
[337]	2252
[492]	2253	# Adding three variables with the station location, longitude, latitude and height
	2254	if obskind == 'single-station':
[500]	2255	adding_station_desc(onewnc,stationdesc)
[492]	2256
[333]	2257	# Global attributes
	2258	##
	2259	set_attribute(onewnc,'author_nc','Lluis Fita')
	2260	set_attribute(onewnc,'institution_nc','Laboratoire de Meteorology Dynamique, ' + \
	2261	'LMD-Jussieu, UPMC, Paris')
	2262	set_attribute(onewnc,'country_nc','France')
	2263	set_attribute(onewnc,'script_nc',main)
	2264	set_attribute(onewnc,'version_script',version)
	2265	set_attribute(onewnc,'information', \
	2266	'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html')
	2267	set_attribute(onewnc,'simfile',opts.fsim)
	2268	set_attribute(onewnc,'obsfile',opts.fobs)
	2269
	2270	onewnc.sync()
	2271	onewnc.close()
	2272
	2273	print main + ": successfull writting of '" + ofile + "' !!"

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

Download in other formats: