Context Navigation

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

Last change on this file since 406 was 346, checked in by lfita, 10 years ago
Splitting simulated and observed variables Adding time_bnds
File size: 66.3 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
	671	def slice_variable(varobj, dimslice):
	672	""" Function to return a slice of a given variable according to values to its
	673	dimensions
	674	slice_variable(varobj, dimslice)
	675	varobj= object wit the variable
	676	dimslice= [[dimname1]:[value1]\|[[dimname2]:[value2], ...] pairs of dimension
	677	[value]:
	678	* [integer]: which value of the dimension
	679	* -1: all along the dimension
	680	* -9: last value of the dimension
	681	* [beg]@[end] slice from [beg] to [end]
	682	"""
	683	fname = 'slice_variable'
	684
	685	if varobj == 'h':
	686	print fname + '_____________________________________________________________'
	687	print slice_variable.__doc__
	688	quit()
	689
	690	vardims = varobj.dimensions
	691	Ndimvar = len(vardims)
	692
	693	Ndimcut = len(dimslice.split('\|'))
	694	dimsl = dimslice.split('\|')
	695
	696	varvalsdim = []
	697	dimnslice = []
	698
	699	for idd in range(Ndimvar):
	700	for idc in range(Ndimcut):
	701	dimcutn = dimsl[idc].split(':')[0]
	702	dimcutv = dimsl[idc].split(':')[1]
	703	if vardims[idd] == dimcutn:
	704	posfrac = dimcutv.find('@')
	705	if posfrac != -1:
	706	inifrac = int(dimcutv.split('@')[0])
	707	endfrac = int(dimcutv.split('@')[1])
	708	varvalsdim.append(slice(inifrac,endfrac))
	709	dimnslice.append(vardims[idd])
	710	else:
	711	if int(dimcutv) == -1:
	712	varvalsdim.append(slice(0,varobj.shape[idd]))
	713	dimnslice.append(vardims[idd])
	714	elif int(dimcutv) == -9:
	715	varvalsdim.append(int(varobj.shape[idd])-1)
	716	else:
	717	varvalsdim.append(int(dimcutv))
	718	break
	719
	720	varvalues = varobj[tuple(varvalsdim)]
	721
	722	return varvalues, dimnslice
	723
[337]	724	def func_compute_varNOcheck(ncobj, varn):
	725	""" Function to compute variables which are not originary in the file
	726	ncobj= netCDF object file
	727	varn = variable to compute:
	728	'WRFdens': air density from WRF variables
	729	'WRFght': geopotential height from WRF variables
	730	'WRFp': pressure from WRF variables
	731	'WRFrh': relative humidty fom WRF variables
	732	'WRFt': temperature from WRF variables
	733	'WRFz': height from WRF variables
	734	"""
	735	fname = 'compute_varNOcheck'
[333]	736
[337]	737	if varn == 'WRFdens':
	738	# print ' ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \
	739	# 'DNW)/g ...'
	740	grav = 9.81
	741
	742	# Just we need in in absolute values: Size of the central grid cell
	743	## dxval = ncobj.getncattr('DX')
	744	## dyval = ncobj.getncattr('DY')
	745	## mapfac = ncobj.variables['MAPFAC_M'][:]
	746	## area = dxvaldyvalmapfac
	747	dimensions = ncobj.variables['MU'].dimensions
	748
	749	mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
	750	dnw = ncobj.variables['DNW'][:]
	751
	752	varNOcheckv = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \
	753	dtype=np.float)
	754	levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
	755
	756	for it in range(mu.shape[0]):
	757	for iz in range(dnw.shape[1]):
	758	levval.fill(np.abs(dnw[it,iz]))
	759	varNOcheck[it,iz,:,:] = levval
	760	varNOcheck[it,iz,:,:] = mu[it,:,:]*varNOcheck[it,iz,:,:]/grav
	761
	762	elif varn == 'WRFght':
	763	# print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
	764	varNOcheckv = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
	765	dimensions = ncobj.variables['PH'].dimensions
	766
	767	elif varn == 'WRFp':
	768	# print ' ' + fname + ': Retrieving pressure value from WRF as P + PB'
	769	varNOcheckv = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	770	dimensions = ncobj.variables['P'].dimensions
	771
	772	elif varn == 'WRFrh':
	773	# print ' ' + main + ": computing relative humidity from WRF as 'Tetens'" +\
	774	# ' equation (T,P) ...'
	775	p0=100000.
	776	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	777	tk = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
	778	qv = ncobj.variables['QVAPOR'][:]
	779
	780	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
	781	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
	782
	783	varNOcheckv = qv/data2
	784	dimensions = ncobj.variables['P'].dimensions
	785
	786	elif varn == 'WRFt':
	787	# print ' ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
	788	p0=100000.
	789	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	790
	791	varNOcheckv = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
	792	dimensions = ncobj.variables['T'].dimensions
	793
	794	elif varn == 'WRFz':
	795	grav = 9.81
	796	# print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
	797	varNOcheckv = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/grav
	798	dimensions = ncobj.variables['PH'].dimensions
	799
	800	else:
	801	print erromsg
	802	print ' ' + fname + ": variable '" + varn + "' nor ready !!"
	803	quit(-1)
	804
	805	return varNOcheck
	806
	807	class compute_varNOcheck(object):
	808	""" Class to compute variables which are not originary in the file
	809	ncobj= netCDF object file
	810	varn = variable to compute:
	811	'WRFdens': air density from WRF variables
	812	'WRFght': geopotential height from WRF variables
	813	'WRFp': pressure from WRF variables
	814	'WRFrh': relative humidty fom WRF variables
[343]	815	'WRFT': CF-time from WRF variables
[337]	816	'WRFt': temperature from WRF variables
[339]	817	'WRFtd': dew-point temperature from WRF variables
	818	'WRFws': wind speed from WRF variables
	819	'WRFwd': wind direction from WRF variables
[337]	820	'WRFz': height from WRF variables
	821	"""
	822	fname = 'compute_varNOcheck'
	823
	824	def __init__(self, ncobj, varn):
	825
	826	if ncobj is None:
	827	self = None
	828	self.dimensions = None
	829	self.shape = None
	830	self.__values = None
	831	else:
	832	if varn == 'WRFdens':
	833	# print ' ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \
	834	# 'DNW)/g ...'
	835	grav = 9.81
	836
	837	# Just we need in in absolute values: Size of the central grid cell
	838	## dxval = ncobj.getncattr('DX')
	839	## dyval = ncobj.getncattr('DY')
	840	## mapfac = ncobj.variables['MAPFAC_M'][:]
	841	## area = dxvaldyvalmapfac
	842	dimensions = ncobj.variables['MU'].dimensions
	843	shape = ncobj.variables['MU'].shape
	844
	845	mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
	846	dnw = ncobj.variables['DNW'][:]
	847
	848	varNOcheckv = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \
	849	dtype=np.float)
	850	levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
	851
	852	for it in range(mu.shape[0]):
	853	for iz in range(dnw.shape[1]):
	854	levval.fill(np.abs(dnw[it,iz]))
	855	varNOcheck[it,iz,:,:] = levval
	856	varNOcheck[it,iz,:,:] = mu[it,:,:]*varNOcheck[it,iz,:,:]/grav
	857
	858	elif varn == 'WRFght':
	859	# print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
	860	varNOcheckv = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
	861	dimensions = ncobj.variables['PH'].dimensions
	862	shape = ncobj.variables['PH'].shape
	863
	864	elif varn == 'WRFp':
	865	# print ' ' + fname + ': Retrieving pressure value from WRF as P + PB'
	866	varNOcheckv = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	867	dimensions = ncobj.variables['P'].dimensions
	868	shape = ncobj.variables['P'].shape
	869
	870	elif varn == 'WRFrh':
	871	# print ' ' + main + ": computing relative humidity from WRF as 'Tetens'" +\
	872	# ' equation (T,P) ...'
	873	p0=100000.
	874	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	875	tk = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
	876	qv = ncobj.variables['QVAPOR'][:]
	877
	878	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
	879	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
	880
	881	varNOcheckv = qv/data2
	882	dimensions = ncobj.variables['P'].dimensions
	883	shape = ncobj.variables['P'].shape
	884
[343]	885	elif varn == 'WRFT':
	886	# To compute CF-times from WRF kind
	887	#
	888	import datetime as dt
	889
	890	times = ncobj.variables['Times']
	891	dimt = times.shape[0]
	892	varNOcheckv = np.zeros((dimt), dtype=np.float64)
	893	self.unitsval = 'seconds since 1949-12-01 00:00:00'
	894	refdate = datetimeStr_datetime('1949-12-01_00:00:00')
	895
	896	dimensions = tuple([ncobj.variables['Times'].dimensions[0]])
	897	shape = tuple([dimt])
	898
	899	for it in range(dimt):
	900	datevalS = datetimeStr_conversion(times[it,:], 'WRFdatetime', \
	901	'YmdHMS')
	902	dateval = dt.datetime.strptime(datevalS, '%Y%m%d%H%M%S')
	903	difft = dateval - refdate
	904	varNOcheckv[it] = difft.days3600.24. + difft.seconds + \
	905	np.float(int(difft.microseconds/10.e6))
	906
[337]	907	elif varn == 'WRFt':
	908	# print ' ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
	909	p0=100000.
	910	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	911
	912	varNOcheckv = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
	913	dimensions = ncobj.variables['T'].dimensions
	914	shape = ncobj.variables['P'].shape
	915
[339]	916	elif varn == 'WRFtd':
	917	# print ' ' + main + ': computing dew-point temperature from WRF as inv_potT(T + 300) and Tetens...'
	918	# tacking from: http://en.wikipedia.org/wiki/Dew_point
	919	p0=100000.
	920	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
	921
	922	temp = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
	923
	924	qv = ncobj.variables['QVAPOR'][:]
	925
	926	tk = temp - 273.15
	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	rh = qv/data2
	931
	932	pa = rh * data1/100.
	933	varNOcheckv = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
	934
	935	dimensions = ncobj.variables['T'].dimensions
	936	shape = ncobj.variables['P'].shape
	937
	938	elif varn == 'WRFws':
	939	# print ' ' + main + ': computing wind speed from WRF as SQRT(U2 + V2) ...'
	940	uwind = ncobj.variables['U'][:]
	941	vwind = ncobj.variables['V'][:]
	942	dx = uwind.shape[3]
	943	dy = vwind.shape[2]
	944
	945	# de-staggering
	946	ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
	947	va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])
	948
	949	varNOcheckv = np.sqrt(uaua + vava)
	950	dimensions = tuple(['Time','bottom_top','south_north','west_east'])
	951	shape = ua.shape
	952
	953	elif varn == 'WRFwd':
	954	# print ' ' + main + ': computing wind direction from WRF as ATAN2PI(V,U) ...'
	955	uwind = ncobj.variables['U'][:]
	956	vwind = ncobj.variables['V'][:]
	957	dx = uwind.shape[3]
	958	dy = vwind.shape[2]
	959
	960	# de-staggering
	961	ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
	962	va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])
	963
	964	theta = np.arctan2(ua, va)
	965	dimensions = tuple(['Time','bottom_top','south_north','west_east'])
	966	shape = ua.shape
	967	varNOcheckv = 360.(1. + theta/(2.np.pi))
	968
[337]	969	elif varn == 'WRFz':
	970	grav = 9.81
	971	# print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
	972	varNOcheckv = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/grav
	973	dimensions = ncobj.variables['PH'].dimensions
	974	shape = ncobj.variables['PH'].shape
	975
	976	else:
[339]	977	print errormsg
[337]	978	print ' ' + fname + ": variable '" + varn + "' nor ready !!"
	979	quit(-1)
	980
	981	self.dimensions = dimensions
	982	self.shape = shape
	983	self.__values = varNOcheckv
	984
	985	def __getitem__(self,elem):
	986	return self.__values[elem]
	987
[330]	988	####### ###### ##### #### ### ## #
	989
	990	strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " + \
	991	" and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')"
	992
	993	kindobs=['multi-points', 'single-station', 'trajectory']
	994	strkObs="kind of observations; 'multi-points': multiple individual punctual obs " + \
	995	"(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
	996	"'trajectory': following a trajectory"
[337]	997	simopers = ['sumc','subc','mulc','divc']
	998	opersinf = 'sumc,[constant]: add [constant] to variables values; subc,[constant]: '+ \
	999	'substract [constant] to variables values; mulc,[constant]: multipy by ' + \
	1000	'[constant] to variables values; divc,[constant]: divide by [constant] to ' + \
	1001	'variables values'
[343]	1002	varNOcheck = ['WRFdens', 'WRFght', 'WRFp', 'WRFrh', 'WRFT', 'WRFt', 'WRFtd', 'WRFws',\
[339]	1003	'WRFwd', 'WRFz']
[337]	1004	varNOcheckinf = "'WRFdens': air density from WRF variables; 'WRFght': geopotential"+ \
	1005	" height from WRF variables; 'WRFp': pressure from WRF variables; 'WRFrh': " + \
[343]	1006	"relative humidty fom WRF variables; 'WRFT': CF-time from WRF variables'WRFt'; " + \
	1007	" temperature from WRF variables; 'WRFtd': dew-point temperature from WRF " + \
	1008	"variables; 'WRFws': wind speed from WRF variables; 'WRFwd': wind speed " + \
	1009	"direction from WRF variables; 'WRFz': height from WRF variables"
[330]	1010
[337]	1011	dimshelp = "[DIM]@[simdim]@[obsdim] ',' list of couples of dimensions names from " + \
	1012	"each source ([DIM]='X','Y','Z','T'; None, no value)"
	1013	vardimshelp = "[DIM]@[simvardim]@[obsvardim] ',' list of couples of variables " + \
	1014	"names with dimensions values from each source ([DIM]='X','Y','Z','T'; None, " + \
	1015	"no value, WRFdiagnosted variables also available: " + varNOcheckinf + ")"
	1016	varshelp="[simvar]@[obsvar]@[[oper]@[val]] ',' list of couples of variables to " + \
	1017	"validate and if necessary operation and value operations: " + opersinf + \
	1018	"(WRFdiagnosted variables also available: " + varNOcheckinf + ")"
	1019	statsn = ['minimum', 'maximum', 'mean', 'mean2', 'standard deviation']
[344]	1020	gstatsn = ['bias', 'simobs_mean', 'sim_obsmin', 'sim_obsmax', 'sim_obsmean', 'mae', \
	1021	'rmse', 'r_pearsoncorr', 'p_pearsoncorr']
[346]	1022	ostatsn = ['number of points', 'minimum', 'maximum', 'mean', 'mean2', \
	1023	'standard deviation']
[337]	1024
[330]	1025	parser = OptionParser()
[337]	1026	parser.add_option("-d", "--dimensions", dest="dims", help=dimshelp, metavar="VALUES")
[330]	1027	parser.add_option("-D", "--vardimensions", dest="vardims",
[337]	1028	help=vardimshelp, metavar="VALUES")
[330]	1029	parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs,
	1030	help=strkObs, metavar="FILE")
	1031	parser.add_option("-l", "--stationLocation", dest="stloc",
	1032	help="longitude, latitude and height of the station (only for 'single-station')",
	1033	metavar="FILE")
	1034	parser.add_option("-o", "--observation", dest="fobs",
	1035	help="observations file to validate", metavar="FILE")
	1036	parser.add_option("-s", "--simulation", dest="fsim",
	1037	help="simulation file to validate", metavar="FILE")
[337]	1038	parser.add_option("-t", "--trajectoryfile", dest="trajf",
	1039	help="file with grid points of the trajectory in the simulation grid ('simtrj')",
	1040	metavar="FILE")
[330]	1041	parser.add_option("-v", "--variables", dest="vars",
[337]	1042	help=varshelp, metavar="VALUES")
[330]	1043
	1044	(opts, args) = parser.parse_args()
	1045
	1046	####### #######
	1047	## MAIN
	1048	#######
	1049
	1050	ofile='validation_sim.nc'
	1051
	1052	if opts.dims is None:
	1053	print errormsg
	1054	print ' ' + main + ': No list of dimensions are provided!!'
	1055	print ' a ',' list of values X@[dimxsim]@[dimxobs],...,T@[dimtsim]@[dimtobs]'+\
	1056	' is needed'
	1057	quit(-1)
	1058	else:
[340]	1059	simdims = {}
	1060	obsdims = {}
[330]	1061	print main +': couple of dimensions _______'
	1062	dims = {}
	1063	ds = opts.dims.split(',')
	1064	for d in ds:
	1065	dsecs = d.split('@')
	1066	if len(dsecs) != 3:
	1067	print errormsg
	1068	print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
	1069	print ' [DIM]@[dimnsim]@[dimnobs]'
	1070	quit(-1)
	1071	dims[dsecs[0]] = [dsecs[1], dsecs[2]]
[340]	1072	simdims[dsecs[0]] = dsecs[1]
	1073	obsdims[dsecs[0]] = dsecs[2]
	1074
[333]	1075	print ' ',dsecs[0],':',dsecs[1],',',dsecs[2]
[330]	1076
	1077	if opts.vardims is None:
	1078	print errormsg
	1079	print ' ' + main + ': No list of variables with dimension values are provided!!'
	1080	print ' a ',' list of values X@[vardimxsim]@[vardimxobs],...,T@' + \
	1081	'[vardimtsim]@[vardimtobs] is needed'
	1082	quit(-1)
	1083	else:
	1084	print main +': couple of variable dimensions _______'
	1085	vardims = {}
	1086	ds = opts.vardims.split(',')
	1087	for d in ds:
	1088	dsecs = d.split('@')
	1089	if len(dsecs) != 3:
	1090	print errormsg
	1091	print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
	1092	print ' [DIM]@[vardimnsim]@[vardimnobs]'
	1093	quit(-1)
	1094	vardims[dsecs[0]] = [dsecs[1], dsecs[2]]
[333]	1095	print ' ',dsecs[0],':',dsecs[1],',',dsecs[2]
[330]	1096
	1097	if opts.obskind is None:
	1098	print errormsg
	1099	print ' ' + main + ': No kind of observations provided !!'
	1100	quit(-1)
	1101	else:
	1102	obskind = opts.obskind
	1103	if obskind == 'single-station':
	1104	if opts.stloc is None:
	1105	print errormsg
	1106	print ' ' + main + ': No station location provided !!'
	1107	quit(-1)
	1108	else:
	1109	stationdesc = [np.float(opts.stloc.split(',')[0]), \
	1110	np.float(opts.stloc.split(',')[1]), np.float(opts.stloc.split(',')[2])]
	1111
	1112	if opts.fobs is None:
	1113	print errormsg
	1114	print ' ' + main + ': No observations file is provided!!'
	1115	quit(-1)
	1116	else:
	1117	if not os.path.isfile(opts.fobs):
	1118	print errormsg
	1119	print ' ' + main + ": observations file '" + opts.fobs + "' does not exist !!"
	1120	quit(-1)
	1121
	1122	if opts.fsim is None:
	1123	print errormsg
	1124	print ' ' + main + ': No simulation file is provided!!'
	1125	quit(-1)
	1126	else:
	1127	if not os.path.isfile(opts.fsim):
	1128	print errormsg
	1129	print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
	1130	quit(-1)
	1131
	1132	if opts.vars is None:
	1133	print errormsg
	1134	print ' ' + main + ': No list of couples of variables is provided!!'
	1135	print ' a ',' list of values [varsim]@[varobs],... is needed'
	1136	quit(-1)
	1137	else:
	1138	valvars = []
[333]	1139	vs = opts.vars.split(',')
[330]	1140	for v in vs:
	1141	vsecs = v.split('@')
[333]	1142	if len(vsecs) < 2:
[330]	1143	print errormsg
	1144	print ' ' + main + ': wrong number of values in:',vsecs, \
	1145	' at least 2 are needed !!'
	1146	print ' [varsim]@[varobs]@[[oper][val]]'
	1147	quit(-1)
[333]	1148	if len(vsecs) > 2:
	1149	if not searchInlist(simopers,vsecs[2]):
	1150	print errormsg
	1151	print main + ": operation on simulation values '" + vsecs[2] + \
	1152	"' not ready !!"
	1153	quit(-1)
	1154
[330]	1155	valvars.append(vsecs)
	1156
	1157	# Openning observations trajectory
	1158	##
	1159	oobs = NetCDFFile(opts.fobs, 'r')
	1160
	1161	valdimobs = {}
	1162	for dn in dims:
	1163	print dn,':',dims[dn]
	1164	if dims[dn][1] != 'None':
	1165	if not oobs.dimensions.has_key(dims[dn][1]):
	1166	print errormsg
	1167	print ' ' + main + ": observations file does not have dimension '" + \
	1168	dims[dn][1] + "' !!"
	1169	quit(-1)
	1170	if vardims[dn][1] != 'None':
	1171	if not oobs.variables.has_key(vardims[dn][1]):
	1172	print errormsg
	1173	print ' ' + main + ": observations file does not have varibale " + \
	1174	"dimension '" + vardims[dn][1] + "' !!"
	1175	quit(-1)
	1176	valdimobs[dn] = oobs.variables[vardims[dn][1]][:]
	1177	else:
	1178	if dn == 'X':
	1179	valdimobs[dn] = stationdesc[0]
	1180	elif dn == 'Y':
	1181	valdimobs[dn] = stationdesc[1]
	1182	elif dn == 'Z':
	1183	valdimobs[dn] = stationdesc[2]
	1184
	1185	osim = NetCDFFile(opts.fsim, 'r')
	1186
	1187	valdimsim = {}
	1188	for dn in dims:
	1189	if not osim.dimensions.has_key(dims[dn][0]):
	1190	print errormsg
[343]	1191	print ' ' + main + ": simulation file '" + opts.fsim + "' does not have " + \
	1192	"dimension '" + dims[dn][0] + "' !!"
	1193	print ' it has: ',osim.dimensions
[330]	1194	quit(-1)
[343]	1195
[337]	1196	if not osim.variables.has_key(vardims[dn][0]) and not \
	1197	searchInlist(varNOcheck,vardims[dn][0]):
[330]	1198	print errormsg
[343]	1199	print ' ' + main + ": simulation file '" + opts.fsim + "' does not have " + \
	1200	"varibale dimension '" + vardims[dn][0] + "' !!"
	1201	print ' it has variables:',osim.variables
[330]	1202	quit(-1)
[337]	1203	if searchInlist(varNOcheck,vardims[dn][0]):
	1204	valdimsim[dn] = compute_varNOcheck(osim, vardims[dn][0])
	1205	else:
	1206	valdimsim[dn] = osim.variables[vardims[dn][0]][:]
[330]	1207
	1208	# General characteristics
[343]	1209	dimtobs = valdimobs['T'].shape[0]
	1210	dimtsim = valdimsim['T'].shape[0]
[330]	1211
	1212	print main +': observational time-steps:',dimtobs,'simulation:',dimtsim
	1213
[337]	1214	notfound = np.zeros((dimtobs), dtype=int)
	1215
[330]	1216	if obskind == 'multi-points':
	1217	trajpos = np.zeros((2,dimt),dtype=int)
[337]	1218	for it in range(dimtobs):
[330]	1219	trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'], \
	1220	[valdimobs['X'][it],valdimobss['Y'][it]])
	1221	elif obskind == 'single-station':
	1222	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],[valdimobs['Y'], \
	1223	valdimobs['X']])
[333]	1224	stationpos = np.zeros((2), dtype=int)
	1225	iid = 0
	1226	for idn in osim.variables[vardims['X'][0]].dimensions:
	1227	if idn == dims['X'][0]:
	1228	stationpos[1] = stsimpos[iid]
	1229	elif idn == dims['Y'][0]:
	1230	stationpos[0] = stsimpos[iid]
	1231
	1232	iid = iid + 1
	1233	print main + ': station point in simulation:', stationpos
[330]	1234	print ' station position:',valdimobs['X'],',',valdimobs['Y']
	1235	print ' simulation coord.:',valdimsim['X'][tuple(stsimpos)],',', \
	1236	valdimsim['Y'][tuple(stsimpos)]
[333]	1237
[330]	1238	elif obskind == 'trajectory':
[337]	1239	if opts.trajf is not None:
	1240	if not os.path.isfile(opts.fsim):
	1241	print errormsg
	1242	print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
	1243	quit(-1)
	1244	else:
	1245	otrjf = NetCDFFile(opts.fsim, 'r')
	1246	trajpos[0,:] = otrjf.variables['obssimtrj'][0]
	1247	trajpos[1,:] = otrjf.variables['obssimtrj'][1]
	1248	otrjf.close()
[330]	1249	else:
[337]	1250	if dims.has_key('Z'):
	1251	trajpos = np.zeros((3,dimtobs),dtype=int)
	1252	for it in range(dimtobs):
	1253	if np.mod(it*100./dimtobs,10.) == 0.:
	1254	print ' trajectory done: ',it*100./dimtobs,'%'
	1255	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \
	1256	[valdimobs['Y'][it],valdimobs['X'][it]])
	1257	stationpos = np.zeros((2), dtype=int)
	1258	iid = 0
	1259	for idn in osim.variables[vardims['X'][0]].dimensions:
	1260	if idn == dims['X'][0]:
	1261	stationpos[1] = stsimpos[iid]
	1262	elif idn == dims['Y'][0]:
	1263	stationpos[0] = stsimpos[iid]
	1264	iid = iid + 1
	1265	if stationpos[0] == 0 and stationpos[1] == 0: notfound[it] = 1
	1266
	1267	trajpos[0,it] = stationpos[0]
	1268	trajpos[1,it] = stationpos[1]
	1269	# In the simulation 'Z' varies with time ... non-hydrostatic model! ;)
	1270	# trajpos[2,it] = index_mat(valdimsim['Z'][it,:,stationpos[0], \
	1271	# stationpos[1]], valdimobs['Z'][it])
	1272	else:
	1273	trajpos = np.zeros((2,dimtobs),dtype=int)
	1274	for it in range(dimtobs):
	1275	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \
	1276	[valdimobs['Y'][it],valdimobss['X'][it]])
	1277	stationpos = np.zeros((2), dtype=int)
	1278	iid = 0
	1279	for idn in osim.variables[vardims['X'][0]].dimensions:
	1280	if idn == dims['X'][0]:
	1281	stationpos[1] = stsimpos[iid]
	1282	elif idn == dims['Y'][0]:
	1283	stationpos[0] = stsimpos[iid]
	1284	iid = iid + 1
[338]	1285	if stationpos[0] == 0 or stationpos[1] == 0: notfound[it] = 1
[330]	1286
[337]	1287	trajpos[0,it] = stationspos[0]
	1288	trajpos[1,it] = stationspos[1]
	1289
	1290	print main + ': not found',np.sum(notfound),'points of the trajectory'
	1291
[330]	1292	# Getting times
	1293	tobj = oobs.variables[vardims['T'][1]]
	1294	obstunits = tobj.getncattr('units')
[343]	1295	if vardims['T'][0] == 'WRFT':
	1296	tsim = valdimsim['T'][:]
	1297	simtunits = 'seconds since 1949-12-01 00:00:00'
	1298	else:
	1299	tsim = osim.variables[vardims['T'][0]]
	1300	simtunits = tobj.getncattr('units')
[330]	1301
[343]	1302	simobstimes = coincident_CFtimes(valdimsim['T'][:], obstunits, simtunits)
[330]	1303
[333]	1304	# Concident times
	1305	##
[337]	1306	coindtvalues0 = []
[333]	1307	for it in range(dimtsim):
	1308	ot = 0
	1309	for ito in range(ot,dimtobs-1):
	1310	if valdimobs['T'][ito] < simobstimes[it] and valdimobs['T'][ito+1] > \
	1311	simobstimes[it]:
	1312	ot = ito
	1313	tdist = simobstimes[it] - valdimobs['T'][ito]
[337]	1314	coindtvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito], \
	1315	tdist])
[330]	1316
[337]	1317	coindtvalues = np.array(coindtvalues0, dtype=np.float)
	1318
	1319	Ncoindt = len(coindtvalues[:,0])
[333]	1320	print main + ': found',Ncoindt,'coincident times between simulation and observations'
	1321
[337]	1322	if Ncoindt == 0:
	1323	print warnmsg
	1324	print main + ': no coincident times found !!'
	1325	print ' stopping it'
	1326	quit(-1)
	1327
[333]	1328	# Validating
	1329	##
	1330
	1331	onewnc = NetCDFFile(ofile, 'w')
	1332
	1333	# Dimensions
	1334	newdim = onewnc.createDimension('time',None)
[346]	1335	newdim = onewnc.createDimension('bnds',2)
[337]	1336	newdim = onewnc.createDimension('obstime',None)
[333]	1337	newdim = onewnc.createDimension('couple',2)
	1338	newdim = onewnc.createDimension('StrLength',StringLength)
[337]	1339	newdim = onewnc.createDimension('xaround',Ngrid*2+1)
	1340	newdim = onewnc.createDimension('yaround',Ngrid*2+1)
[344]	1341	newdim = onewnc.createDimension('gstats',9)
[337]	1342	newdim = onewnc.createDimension('stats',5)
[346]	1343	newdim = onewnc.createDimension('tstats',6)
[333]	1344
	1345	# Variable dimensions
	1346	##
	1347	newvar = onewnc.createVariable('obstime','f8',('time'))
	1348	basicvardef(newvar, 'obstime', 'time observations', obstunits )
	1349	set_attribute(newvar, 'calendar', 'standard')
[346]	1350	set_attribute(newvar, 'bounds', 'time_bnds')
[337]	1351	newvar[:] = coindtvalues[:,3]
[333]	1352
	1353	newvar = onewnc.createVariable('couple', 'c', ('couple','StrLength'))
	1354	basicvardef(newvar, 'couple', 'couples of values', '-')
	1355	writing_str_nc(newvar, ['sim','obs'], StringLength)
	1356
[337]	1357	newvar = onewnc.createVariable('statistics', 'c', ('stats','StrLength'))
	1358	basicvardef(newvar, 'statistics', 'statitics from values', '-')
	1359	writing_str_nc(newvar, statsn, StringLength)
	1360
[344]	1361	newvar = onewnc.createVariable('gstatistics', 'c', ('gstats','StrLength'))
	1362	basicvardef(newvar, 'gstatistics', 'global statitics from values', '-')
	1363	writing_str_nc(newvar, gstatsn, StringLength)
	1364
	1365	newvar = onewnc.createVariable('tstatistics', 'c', ('tstats','StrLength'))
	1366	basicvardef(newvar, 'tstatistics', 'statitics from values along time', '-')
	1367	writing_str_nc(newvar, ostatsn, StringLength)
	1368
[337]	1369	if obskind == 'trajectory':
	1370	if dims.has_key('Z'):
	1371	newdim = onewnc.createDimension('trj',3)
	1372	else:
	1373	newdim = onewnc.createDimension('trj',2)
	1374
	1375	newvar = onewnc.createVariable('obssimtrj','i',('obstime','trj'))
	1376	basicvardef(newvar, 'obssimtrj', 'trajectory on the simulation grid', '-')
	1377	newvar[:] = trajpos.transpose()
	1378
	1379	if dims.has_key('Z'):
	1380	newdim = onewnc.createDimension('simtrj',4)
	1381	trjsim = np.zeros((4,Ncoindt), dtype=int)
	1382	trjsimval = np.zeros((4,Ncoindt), dtype=np.float)
	1383	else:
	1384	newdim = onewnc.createDimension('simtrj',3)
	1385	trjsim = np.zeros((3,Ncoindt), dtype=int)
	1386	trjsimval = np.zeros((3,Ncoindt), dtype=np.float)
	1387
[333]	1388	Nvars = len(valvars)
	1389	for ivar in range(Nvars):
	1390	simobsvalues = []
	1391
	1392	varsimobs = valvars[ivar][0] + '_' + valvars[ivar][1]
[340]	1393	print ' ' + varsimobs + '... .. .'
[333]	1394
	1395	if not oobs.variables.has_key(valvars[ivar][1]):
	1396	print errormsg
	1397	print ' ' + main + ": observations file has not '" + valvars[ivar][1] + \
	1398	"' !!"
	1399	quit(-1)
[337]	1400
[333]	1401	if not osim.variables.has_key(valvars[ivar][0]):
[337]	1402	if not searchInlist(varNOcheck, valvars[ivar][0]):
	1403	print errormsg
	1404	print ' ' + main + ": simulation file has not '" + valvars[ivar][0] + \
	1405	"' !!"
	1406	quit(-1)
	1407	else:
	1408	ovsim = compute_varNOcheck(osim, valvars[ivar][0])
	1409	else:
	1410	ovsim = osim.variables[valvars[ivar][0]]
[333]	1411
[340]	1412	for idn in ovsim.dimensions:
	1413	if not searchInlist(simdims.values(),idn):
	1414	print errormsg
	1415	print ' ' + main + ": dimension '" + idn + "' of variable '" + \
	1416	valvars[ivar][0] + "' not provided as reference coordinate [X,Y,Z,T] !!"
	1417	quit(-1)
	1418
[333]	1419	ovobs = oobs.variables[valvars[ivar][1]]
[344]	1420
	1421	# Simulated values spatially around coincident times
[337]	1422	if dims.has_key('Z'):
	1423	simobsSvalues = np.zeros((Ncoindt, Ngrid2+1, Ngrid2+1, Ngrid*2+1), \
	1424	dtype = np.float)
	1425	else:
	1426	simobsSvalues = np.zeros((Ncoindt, Ngrid2+1, Ngrid2+1), dtype = np.float)
[333]	1427
[344]	1428	# Observed values temporally around coincident times
	1429	simobsTvalues = {}
	1430	simobsTtvalues = np.zeros((Ncoindt,2), dtype=np.float)
	1431
[333]	1432	if obskind == 'multi-points':
	1433	for it in range(Ncoindt):
[337]	1434	slicev = dims['X'][0] + ':' + str(trajpos[2,it]) + '\|' + \
	1435	dims['Y'][0]+ ':' + str(trajpos[1,it]) + '\|' + \
	1436	dims['T'][0]+ ':' + str(coindtvalues[it][0])
	1437	slicevar, dimslice = slice_variable(ovsim, slicev)
	1438	simobsvalues.append([ slicevar, ovobs[coindtvalues[it][1]]])
	1439	slicev = dims['X'][0] + ':' + str(trajpos[2,it]-Ngrid) + '@' + \
	1440	str(trajpos[2,it]+Ngrid) + '\|' + dims['Y'][0] + ':' + \
	1441	str(trajpos[1,it]-Ngrid) + '@' + str(trajpos[1,it]+Ngrid) + '\|' + \
[333]	1442	dims['T'][0]+':'+str(coindtvalues[it][0])
	1443	slicevar, dimslice = slice_variable(ovsim, slicev)
[337]	1444	simobsSvalues[it,:,:] = slicevar
	1445
[333]	1446	elif obskind == 'single-station':
	1447	for it in range(Ncoindt):
[343]	1448	ito = int(coindtvalues[it,1])
	1449	slicev = dims['X'][0] + ':' + str(stationpos[1]) + '\|' + \
	1450	dims['Y'][0] + ':' + str(stationpos[0]) + '\|' + \
	1451	dims['T'][0] + ':' + str(int(coindtvalues[it][0]))
[333]	1452	slicevar, dimslice = slice_variable(ovsim, slicev)
[343]	1453	if ovobs[int(ito)] != ovobs[int(ito)]:
	1454	simobsvalues.append([ slicevar, fillValueF])
	1455	else:
	1456	simobsvalues.append([ slicevar, ovobs[int(ito)]])
	1457	slicev = dims['X'][0] + ':' + str(stationpos[1]-Ngrid) + '@' + \
	1458	str(stationpos[1]+Ngrid+1) + '\|' + dims['Y'][0] + ':' + \
	1459	str(stationpos[0]-Ngrid) + '@' + str(stationpos[0]+Ngrid+1) + '\|' + \
	1460	dims['T'][0] + ':' + str(int(coindtvalues[it,0]))
[337]	1461	slicevar, dimslice = slice_variable(ovsim, slicev)
	1462	simobsSvalues[it,:,:] = slicevar
[343]	1463
[344]	1464	if it == 0:
	1465	itoi = 0
	1466	itof = int(coindtvalues[it,1]) / 2
	1467	elif it == Ncoindt-1:
	1468	itoi = int( (ito + int(coindtvalues[it-1,1])) / 2)
	1469	itof = int(coindtvalues[it,1])
	1470	else:
	1471	itod = int( (ito - int(coindtvalues[it-1,1])) / 2 )
	1472	itoi = ito - itod
	1473	itod = int( (int(coindtvalues[it+1,1]) - ito) / 2 )
	1474	itof = ito + itod
	1475
	1476	slicev = dims['T'][1] + ':' + str(itoi) + '@' + str(itof + 1)
	1477
	1478	slicevar, dimslice = slice_variable(ovobs, slicev)
	1479	simobsTvalues[str(it)] = slicevar
	1480
	1481	simobsTtvalues[it,0] = valdimobs['T'][itoi]
	1482	simobsTtvalues[it,1] = valdimobs['T'][itof]
	1483
[333]	1484	elif obskind == 'trajectory':
	1485	if dims.has_key('Z'):
	1486	for it in range(Ncoindt):
[338]	1487	ito = int(coindtvalues[it,1])
	1488	if notfound[ito] == 0:
	1489	trajpos[2,ito] = index_mat(valdimsim['Z'][coindtvalues[it,0],:, \
	1490	trajpos[1,ito],trajpos[0,ito]], valdimobs['Z'][ito])
[337]	1491	slicev = dims['X'][0]+':'+str(trajpos[0,ito]) + '\|' + \
	1492	dims['Y'][0]+':'+str(trajpos[1,ito]) + '\|' + \
	1493	dims['Z'][0]+':'+str(trajpos[2,ito]) + '\|' + \
	1494	dims['T'][0]+':'+str(int(coindtvalues[it,0]))
	1495	slicevar, dimslice = slice_variable(ovsim, slicev)
	1496	simobsvalues.append([ slicevar, ovobs[int(ito)]])
[338]	1497	minx = np.max([trajpos[0,ito]-Ngrid,0])
	1498	maxx = np.min([trajpos[0,ito]+Ngrid+1,ovsim.shape[3]])
	1499	miny = np.max([trajpos[1,ito]-Ngrid,0])
	1500	maxy = np.min([trajpos[1,ito]+Ngrid+1,ovsim.shape[2]])
	1501	minz = np.max([trajpos[2,ito]-Ngrid,0])
	1502	maxz = np.min([trajpos[2,ito]+Ngrid+1,ovsim.shape[1]])
	1503
	1504	slicev = dims['X'][0] + ':' + str(minx) + '@' + str(maxx) + '\|' +\
	1505	dims['Y'][0] + ':' + str(miny) + '@' + str(maxy) + '\|' + \
	1506	dims['Z'][0] + ':' + str(minz) + '@' + str(maxz) + '\|' + \
	1507	dims['T'][0] + ':' + str(int(coindtvalues[it,0]))
[337]	1508	slicevar, dimslice = slice_variable(ovsim, slicev)
[338]	1509
	1510	sliceS = []
	1511	sliceS.append(it)
	1512	sliceS.append(slice(0,maxz-minz))
	1513	sliceS.append(slice(0,maxy-miny))
	1514	sliceS.append(slice(0,maxx-minx))
[340]	1515
[338]	1516	simobsSvalues[tuple(sliceS)] = slicevar
[337]	1517	if ivar == 0:
	1518	trjsim[0,it] = trajpos[0,ito]
	1519	trjsim[1,it] = trajpos[1,ito]
	1520	trjsim[2,it] = trajpos[2,ito]
	1521	trjsim[3,it] = coindtvalues[it,0]
	1522	else:
	1523	simobsvalues.append([fillValueF, fillValueF])
	1524	simobsSvalues[it,:,:,:]= np.ones((Ngrid2+1,Ngrid2+1,Ngrid*2+1),\
	1525	dtype = np.float)*fillValueF
[333]	1526	else:
	1527	for it in range(Ncoindt):
[337]	1528	if notfound[it] == 0:
	1529	ito = coindtvalues[it,1]
	1530	slicev = dims['X'][0]+':'+str(trajpos[2,ito]) + '\|' + \
	1531	dims['Y'][0]+':'+str(trajpos[1,ito]) + '\|' + \
	1532	dims['T'][0]+':'+str(coindtvalues[ito,0])
	1533	slicevar, dimslice = slice_variable(ovsim, slicev)
	1534	simobsvalues.append([ slicevar, ovobs[coindtvalues[it,1]]])
	1535	slicev = dims['X'][0] + ':' + str(trajpos[0,it]-Ngrid) + '@' + \
	1536	str(trajpos[0,it]+Ngrid) + '\|' + dims['Y'][0] + ':' + \
	1537	str(trajpos[1,it]-Ngrid) + '@' + str(trajpos[1,it]+Ngrid) + \
	1538	'\|' + dims['T'][0] + ':' + str(coindtvalues[it,0])
	1539	slicevar, dimslice = slice_variable(ovsim, slicev)
	1540	simobsSvalues[it,:,:] = slicevar
	1541	else:
	1542	simobsvalues.append([fillValue, fillValue])
	1543	simobsSvalues[it,:,:] = np.ones((Ngrid2+1,Ngrid2+1), \
	1544	dtype = np.float)*fillValueF
[333]	1545	print simobsvalues[varsimobs][:][it]
	1546
[339]	1547	arrayvals = np.array(simobsvalues)
	1548	if len(valvars[ivar]) > 2:
	1549	const=np.float(valvars[ivar][3])
	1550	if valvars[ivar][2] == 'sumc':
	1551	simobsSvalues = simobsSvalues + const
	1552	arrayvals[:,0] = arrayvals[:,0] + const
	1553	elif valvars[ivar][2] == 'subc':
	1554	simobsSvalues = simobsSvalues - const
	1555	arrayvals[:,0] = arrayvals[:,0] - const
	1556	elif valvars[ivar][2] == 'mulc':
	1557	simobsSvalues = simobsSvalues * const
	1558	arrayvals[:,0] = arrayvals[:,0] * const
	1559	elif valvars[ivar][2] == 'divc':
	1560	simobsSvalues = simobsSvalues / const
	1561	arrayvals[:,0] = arrayvals[:,0] / const
	1562	else:
	1563	print errormsg
	1564	print ' ' + fname + ": operation '" + valvars[ivar][2] + "' not ready!!"
	1565	quit(-1)
	1566
[343]	1567	# statisics sim
	1568	simstats = np.zeros((5), dtype=np.float)
	1569	simstats[0] = np.min(arrayvals[:,0])
	1570	simstats[1] = np.max(arrayvals[:,0])
	1571	simstats[2] = np.mean(arrayvals[:,0])
	1572	simstats[3] = np.mean(arrayvals[:,0]*arrayvals[:,0])
	1573	simstats[4] = np.sqrt(simstats[3] - simstats[2]*simstats[2])
	1574
	1575	# statisics obs
[344]	1576	obsmask = ma.masked_equal(arrayvals[:,1], fillValueF)
[343]	1577	obsmask2 = obsmask*obsmask
	1578
	1579	obsstats = np.zeros((5), dtype=np.float)
	1580	obsstats[0] = obsmask.min()
	1581	obsstats[1] = obsmask.max()
	1582	obsstats[2] = obsmask.mean()
	1583	obsstats[3] = obsmask2.mean()
	1584	obsstats[4] = np.sqrt(obsstats[3] - obsstats[2]*obsstats[2])
	1585
	1586	# Statistics sim-obs
[344]	1587	simobsstats = np.zeros((9), dtype=np.float)
[343]	1588	diffvals = np.zeros((Ncoindt), dtype=np.float)
	1589
[344]	1590	diffvals = arrayvals[:,0] - obsmask
	1591
[343]	1592	simobsstats[0] = simstats[0] - obsstats[0]
[344]	1593	simobsstats[1] = np.mean(arrayvals[:,0]*obsmask)
	1594	simobsstats[2] = np.min(diffvals)
	1595	simobsstats[3] = np.max(diffvals)
	1596	simobsstats[4] = np.mean(diffvals)
	1597	simobsstats[5] = np.mean(np.abs(diffvals))
	1598	simobsstats[6] = np.sqrt(np.mean(diffvals*diffvals))
	1599	simobsstats[7], simobsstats[8] = sts.pearsonr(arrayvals[:,0], arrayvals[:,1])
[343]	1600
[344]	1601	# Statistics around sim values
[337]	1602	aroundstats = np.zeros((5,Ncoindt), dtype=np.float)
	1603	for it in range(Ncoindt):
	1604	aroundstats[0,it] = np.min(simobsSvalues[it,])
	1605	aroundstats[1,it] = np.max(simobsSvalues[it,])
	1606	aroundstats[2,it] = np.mean(simobsSvalues[it,])
	1607	aroundstats[3,it] = np.mean(simobsSvalues[it,]*simobsSvalues[it,])
	1608	aroundstats[4,it] = np.sqrt(aroundstats[3,it] - aroundstats[2,it]* \
	1609	aroundstats[2,it])
	1610
[344]	1611	# Statistics around obs values
[346]	1612	aroundostats = np.zeros((6,Ncoindt), dtype=np.float)
[344]	1613
	1614	for it in range(Ncoindt):
	1615	obsmask = ma.masked_equal(simobsTvalues[str(it)], fillValueF)
	1616	obsmask2 = obsmask*obsmask
	1617
	1618	aroundostats[0,it] = len(obsmask.flatten())
[346]	1619	aroundostats[1,it] = obsmask.min()
	1620	aroundostats[2,it] = obsmask.max()
	1621	aroundostats[3,it] = obsmask.mean()
	1622	aroundostats[4,it] = obsmask2.mean()
	1623	aroundostats[5,it] = np.sqrt(aroundostats[4,it] - aroundostats[3,it]* \
	1624	aroundostats[3,it])
[344]	1625
[346]	1626	# sim Values to netCDF
	1627	newvar = onewnc.createVariable(valvars[ivar][0], 'f', ('time'), \
[337]	1628	fill_value=fillValueF)
[346]	1629	descvar = 'simulated: ' + valvars[ivar][0]
	1630	basicvardef(newvar, valvars[ivar][0], descvar, ovobs.getncattr('units'))
	1631	newvar[:] = arrayvals[:,0]
[333]	1632
[346]	1633	# obs Values to netCDF
	1634	newvar = onewnc.createVariable(valvars[ivar][1], 'f', ('time'), \
	1635	fill_value=fillValueF)
	1636	descvar = 'observed: ' + valvars[ivar][1]
	1637	basicvardef(newvar, valvars[ivar][1], descvar, ovobs.getncattr('units'))
	1638	newvar[:] = arrayvals[:,1]
	1639
[337]	1640	# Around values
[341]	1641	if not onewnc.variables.has_key(valvars[ivar][0] + 'around'):
	1642	if dims.has_key('Z'):
	1643	if not onewnc.dimensions.has_key('zaround'):
	1644	newdim = onewnc.createDimension('zaround',Ngrid*2+1)
	1645	newvar = onewnc.createVariable(valvars[ivar][0] + 'around', 'f', \
	1646	('time','zaround','yaround','xaround'), fill_value=fillValueF)
	1647	else:
	1648	newvar = onewnc.createVariable(valvars[ivar][0] + 'around', 'f', \
	1649	('time','yaround','xaround'), fill_value=fillValueF)
[339]	1650
[341]	1651	descvar = 'around simulated values +/- grid values: ' + valvars[ivar][0]
	1652	basicvardef(newvar, varsimobs + 'around', descvar, ovobs.getncattr('units'))
	1653	newvar[:] = simobsSvalues
[337]	1654
[344]	1655	# sim Statistics
	1656	if not searchInlist(onewnc.variables,valvars[ivar][0] + 'stsim'):
	1657	newvar = onewnc.createVariable(valvars[ivar][0] + 'stsim', 'f', ('stats'), \
	1658	fill_value=fillValueF)
	1659	descvar = 'simulated statisitcs: ' + valvars[ivar][0]
	1660	basicvardef(newvar, valvars[ivar][0] + 'stsim', descvar, ovobs.getncattr('units'))
	1661	newvar[:] = simstats
	1662
	1663	# obs Statistics
	1664	if not searchInlist(onewnc.variables,valvars[ivar][1] + 'stobs'):
	1665	newvar = onewnc.createVariable(valvars[ivar][1] + 'stobs', 'f', ('stats'), \
	1666	fill_value=fillValueF)
	1667	descvar = 'observed statisitcs: ' + valvars[ivar][1]
	1668	basicvardef(newvar, valvars[ivar][1] + 'stobs', descvar, \
	1669	ovobs.getncattr('units'))
	1670	newvar[:] = obsstats
	1671
	1672	# sim-obs Statistics
	1673	if not searchInlist(onewnc.variables,varsimobs + 'st'):
	1674	newvar = onewnc.createVariable(varsimobs + 'st', 'f', ('gstats'), \
	1675	fill_value=fillValueF)
	1676	descvar = 'simulated-observed statisitcs: ' + varsimobs
	1677	basicvardef(newvar, varsimobs + 'st', descvar, ovobs.getncattr('units'))
	1678	newvar[:] = simobsstats
	1679
	1680	# around sim Statistics
	1681	if not searchInlist(onewnc.variables,valvars[ivar][0] + 'staround'):
[341]	1682	newvar = onewnc.createVariable(valvars[ivar][0] + 'staround', 'f', \
	1683	('time','stats'), fill_value=fillValueF)
[344]	1684	descvar = 'around (' + str(Ngrid) + ', ' + str(Ngrid) + \
	1685	') simulated statisitcs: ' + valvars[ivar][0]
	1686	basicvardef(newvar, valvars[ivar][0] + 'staround', descvar, \
	1687	ovobs.getncattr('units'))
[341]	1688	newvar[:] = aroundstats.transpose()
[337]	1689
[346]	1690	if not searchInlist(onewnc.variables, 'time_bnds'):
	1691	newvar = onewnc.createVariable('time_bnds','f8',('time','bnds'))
	1692	basicvardef(newvar, 'time_bnds', 'time', obstunits )
	1693	set_attribute(newvar, 'calendar', 'standard')
	1694	newvar[:] = simobsTtvalues
	1695
[344]	1696	# around obs Statistics
	1697	if not searchInlist(onewnc.variables,valvars[ivar][1] + 'staround'):
	1698	newvar = onewnc.createVariable(valvars[ivar][1] + 'staround', 'f', \
	1699	('time','tstats'), fill_value=fillValueF)
	1700	descvar = 'around temporal observed statisitcs: ' + valvars[ivar][1]
	1701	basicvardef(newvar, valvars[ivar][1] + 'staround', descvar, \
	1702	ovobs.getncattr('units'))
[346]	1703	set_attribute(newvar, 'cell_methods', 'statistics')
	1704
[344]	1705	newvar[:] = aroundostats.transpose()
	1706
[341]	1707	onewnc.sync()
[333]	1708
[337]	1709	newvar = onewnc.createVariable('simtrj','i',('time','simtrj'))
	1710	basicvardef(newvar, 'simtrj', 'coordinates [X,Y,Z,T] of the coincident trajectory ' +\
	1711	'in sim', obstunits)
	1712	newvar[:] = trjsim.transpose()
	1713
[333]	1714	# Global attributes
	1715	##
	1716	set_attribute(onewnc,'author_nc','Lluis Fita')
	1717	set_attribute(onewnc,'institution_nc','Laboratoire de Meteorology Dynamique, ' + \
	1718	'LMD-Jussieu, UPMC, Paris')
	1719	set_attribute(onewnc,'country_nc','France')
	1720	set_attribute(onewnc,'script_nc',main)
	1721	set_attribute(onewnc,'version_script',version)
	1722	set_attribute(onewnc,'information', \
	1723	'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html')
	1724	set_attribute(onewnc,'simfile',opts.fsim)
	1725	set_attribute(onewnc,'obsfile',opts.fobs)
	1726
	1727	onewnc.sync()
	1728	onewnc.close()
	1729
	1730	print main + ": successfull writting of '" + ofile + "' !!"

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

Download in other formats: