Context Navigation

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

Last change on this file since 525 was 525, checked in by lfita, 9 years ago
Fixing type in 'if'
File size: 86.5 KB

Line
1
2	# L. Fita, LMD-Jussieu. February 2015
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
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
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
7
8	import numpy as np
9	import os
10	import re
11	from optparse import OptionParser
12	from netCDF4 import Dataset as NetCDFFile
13	from scipy import stats as sts
14	import numpy.ma as ma
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
28	StringLength = 50
29
30	# Number of grid points to take as 'environment' around the observed point
31	Ngrid = 1
32
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
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
173	>>> index_2mat(np.arange(27).reshape(3,3,3),np.arange(100,127).reshape(3,3,3),[22,111])
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
193	Ndims = len(matAshape)
194	# valpos = np.ones((Ndims), dtype=int)*-1.
195	valpos = np.zeros((Ndims), dtype=int)
196
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,'!!'
201	return valpos
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,'!!'
206	return valpos
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
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)
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
231	def index_mat(matA,val):
232	""" Function to provide the coordinates of a given value inside a matrix
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
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
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
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
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
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
732	'WRFwds': surface wind direction from WRF variables
733	'WRFwss': surface wind speed from WRF variables
734	'WRFz': height from WRF variables
735	"""
736	fname = 'compute_varNOcheck'
737
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
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
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
827	'TSrhs': surface relative humidty fom TS variables
828	'WRFrhs': surface relative humidty fom WRF variables
829	'WRFT': CF-time from WRF variables
830	'WRFt': temperature from WRF variables
831	'WRFtd': dew-point temperature from WRF variables
832	'WRFwd': wind direction from WRF variables
833	'TSwds': surface wind direction from TS variables
834	'WRFwds': surface wind direction from WRF variables
835	'WRFws': wind speed from WRF variables
836	'TSwss': surface wind speed from TS variables
837	'WRFwss': surface wind speed from WRF variables
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'][:]
893	tk = (ncobj.variables['T'][:])(p/p0)*(2./7.)
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
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
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
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
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
986	elif varn == 'WRFwd':
987	# print ' ' + main + ': computing wind direction from WRF as ATAN2PI(V,U) ...'
988	uwind = ncobj.variables['U'][:]
989	vwind = ncobj.variables['V'][:]
990	dx = uwind.shape[3]
991	dy = vwind.shape[2]
992
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
997	theta = np.arctan2(ua, va)
998	theta = np.where(theta < 0., theta + 2.*np.pi, theta)
999	varNOcheckv = 360.theta/(2.np.pi)
1000
1001	dimensions = tuple(['Time','bottom_top','south_north','west_east'])
1002	shape = ua.shape
1003
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'][:])
1007	theta = np.where(theta < 0., theta + 2.*np.pi, theta)
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'][:])
1016	theta = np.where(theta < 0., theta + 2.*np.pi, theta)
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) ...'
1024	uwind = ncobj.variables['U'][:]
1025	vwind = ncobj.variables['V'][:]
1026	dx = uwind.shape[3]
1027	dy = vwind.shape[2]
1028
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
1033	varNOcheckv = np.sqrt(uaua + vava)
1034	dimensions = tuple(['Time','bottom_top','south_north','west_east'])
1035	shape = ua.shape
1036
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
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:
1061	print errormsg
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
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
1079	newvar = onc.createVariable( 'station', 'c', ('StrLength'))
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
1092	newvar = onc.createVariable( lonname, 'f4', ('nst'))
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
1104	newvar = onc.createVariable( latname, 'f4', ('nst'))
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
1116	newvar = onc.createVariable( heightname, 'f4', ('nst'))
1117	basicvardef(newvar, heightname, 'height above sea level', 'm' )
1118	newvar[:] = stdesc[3]
1119
1120	return
1121
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
1147	def getting_ValidationValues(okind, dimt, ds, trjpos, ovs, ovo, tvalues, oFill, Ng, kvals):
1148	""" Function to get the values to validate accroding to the type of observation
1149	okind= observational kind
1150	dimt= initial number of values to retrieve
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)
1156	oFill= Fill Value for the observations
1157	Ng= number of grid points around the observation
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
1162	return:
1163	sovalues= simulated values at the observation point and time
1164	soSvalues= values Ngrid points around the simulated point
1165	soTtvalues= inital/ending period between two consecutive obsevations (for `single-station')
1166	trjs= trajectory on the simulation space
1167	"""
1168	fname = 'getting_ValidationValues'
1169
1170	sovalues = []
1171
1172	if kvals == 'instantaneous':
1173	dimtf = dimt
1174	elif kvals == 'tbackwardSmean':
1175	print ' ' + fname + ':',kvals,'!!'
1176	uniqt = np.unique(tvalues[:,3])
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]
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]
1187	else:
1188	posprev = postf[uniqt[it-1]][1]
1189	posit = list(tvalues[:,3]).index(uniqt[it])
1190	postf[uniqt[it]] = [posprev+1, posit+1]
1191	elif kvals == 'tbackwardOmean':
1192	print ' ' + fname + ':',kvals,'!!'
1193	uniqt = np.unique(tvalues[:,2])
1194	dimtf = len(uniqt)
1195	print ' initially we got',dimt,'values which will become',dimtf
1196	else:
1197	print errormsg
1198	print ' ' + fname + ": kind of values '" + kvals + "' not ready!!"
1199	quit(-1)
1200
1201	# Simulated values spatially around
1202	if ds.has_key('Z'):
1203	soSvalues = np.zeros((dimt, Ng2+1, Ng2+1, Ng*2+1), dtype = np.float)
1204	if okind == 'trajectory':
1205	trjs = np.zeros((4,dimt), dtype=int)
1206	else:
1207	trjs = None
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:
1213	trjs = None
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):
1222	slicev = ds['X'][0] + ':' + str(trjpos[2,it]) + '\|' + ds['Y'][0] + \
1223	':' + str(trjpos[1,it]) + '\|' + ds['T'][0]+ ':' + str(tvalues[it][0])
1224	slicevar, dimslice = slice_variable(ovs, slicev)
1225	sovalues.append([ slicevar, ovo[tvalues[it][1]]])
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) + '\|' + \
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'):
1237	slicev = ds['X'][0] + ':' + str(stationpos[1]) + '\|' + \
1238	ds['Y'][0] + ':' + str(stationpos[0]) + '\|' + \
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'):
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) + '\|' + \
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:
1279	trjpos[2,ito] = index_mat(valdimsim['Z'][tvalues[it,0],:, \
1280	trjpos[1,ito],trjpos[0,ito]], valdimobs['Z'][ito])
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]) + '\|' + \
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
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) + '\|' + \
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])
1314	soSvalues[it,:,:,:]= np.ones((Ng2+1,Ng2+1,Ng*2+1), \
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]
1321	slicev = ds['X'][0]+':'+str(trjpos[2,ito]) + '\|' + \
1322	ds['Y'][0]+':'+str(trjpos[1,ito]) + '\|' + \
1323	ds['T'][0]+':'+str(tvalues[ito,0])
1324	slicevar, dimslice = slice_variable(ovs, slicev)
1325	sovalues.append([ slicevar, ovo[tvalues[it,1]]])
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) + \
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])
1334	soSvalues[it,:,:] = np.ones((Ng2+1,Ng2+1), \
1335	dtype = np.float)*fillValueF
1336	print sovalues[varsimobs][:][it]
1337	else:
1338	print errormsg
1339	print ' ' + fname + ": observatino kind '" + okind + "' not ready!!"
1340	quit(-1)
1341
1342	# Re-arranging final values
1343	##
1344	if kvals == 'instantaneous':
1345	return sovalues, soSvalues, soTtvalues, trjs
1346
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
1367	for it in range(1,dimtf):
1368	tv = uniqt[it]
1369	intv = postf[tv]
1370
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'):
1381	if okind == 'trajectory':
1382	for ip in range(4):
1383	ftrjs[ip,it] = np.mean(trjs[ip,intv[0]:intv[1]])
1384	for iz in range(2*Ng+1):
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]: \
1388	intv[1],iz,iy,ix])
1389	else:
1390	if okind == 'trajectory':
1391	for ip in range(3):
1392	ftrjs[ip,it] = np.mean(trjs[ip,intv[0]:intv[1]])
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], \
1396	iy,ix])
1397	fsoTtvalues[it,0] = soTtvalues[intv[0],0]
1398	fsoTtvalues[it,1] = soTtvalues[intv[1],0]
1399
1400	return sovalues, soSvalues, soTtvalues, trjs
1401
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
1409	return
1410
1411
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"
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'
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"
1441
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 " + \
1448	"validate and if necessary operation and value (sim. values) available " + \
1449	"operations: " + opersinf + " (WRFdiagnosted variables also available: " + \
1450	varNOcheckinf + ")"
1451	statsn = ['minimum', 'maximum', 'mean', 'mean2', 'standard deviation']
1452	gstatsn = ['bias', 'simobs_mean', 'sim_obsmin', 'sim_obsmax', 'sim_obsmean', 'mae', \
1453	'rmse', 'r_pearsoncorr', 'p_pearsoncorr', 'deviation_of_residuals_SDR', \
1454	'indef_of_efficiency_IOE', 'index_of_agreement_IOA', 'fractional_mean_bias_FMB']
1455	ostatsn = ['number of points', 'minimum', 'maximum', 'mean', 'mean2', \
1456	'standard deviation']
1457
1458	parser = OptionParser()
1459	parser.add_option("-d", "--dimensions", dest="dims", help=dimshelp, metavar="VALUES")
1460	parser.add_option("-D", "--vardimensions", dest="vardims",
1461	help=vardimshelp, metavar="VALUES")
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",
1465	help="name (\| for spaces), longitude, latitude and height of the station (only for 'single-station')",
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")
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")
1474	parser.add_option("-v", "--variables", dest="vars",
1475	help=varshelp, metavar="VALUES")
1476
1477	(opts, args) = parser.parse_args()
1478
1479	####### ###### ##### #### ### ## #
1480	# Number of different statistics according to the temporal coincidence
1481	# 0: Exact time
1482	# 1: Simulation values closest to observed times
1483	# 2: Simulation values between consecutive observed times
1484	Nstsim = 3
1485
1486	stdescsim = ['E', 'C', 'B']
1487	prestdescsim = ['exact', 'closest', 'between']
1488	Lstdescsim = ['exact time', 'cloest time', 'between observational time-steps']
1489
1490	####### #######
1491	## MAIN
1492	#######
1493
1494	ofile='validation_sim.nc'
1495
1496	if opts.dims is None:
1497	print errormsg
1498	print ' ' + main + ': No list of dimensions are provided!!'
1499	print ' a ',' list of values X@[dimxsim]@[dimxobs],...,T@[dimtsim]@[dimtobs]'+\
1500	' is needed'
1501	quit(-1)
1502	else:
1503	simdims = {}
1504	obsdims = {}
1505	print main +': couple of dimensions _______'
1506	dims = {}
1507	ds = opts.dims.split(',')
1508	for d in ds:
1509	dsecs = d.split('@')
1510	if len(dsecs) != 3:
1511	print errormsg
1512	print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
1513	print ' [DIM]@[dimnsim]@[dimnobs]'
1514	quit(-1)
1515	if dsecs[1] != 'None':
1516	dims[dsecs[0]] = [dsecs[1], dsecs[2]]
1517	simdims[dsecs[0]] = dsecs[1]
1518	obsdims[dsecs[0]] = dsecs[2]
1519
1520	print ' ',dsecs[0],':',dsecs[1],',',dsecs[2]
1521
1522	if opts.vardims is None:
1523	print errormsg
1524	print ' ' + main + ': No list of variables with dimension values are provided!!'
1525	print ' a ',' list of values X@[vardimxsim]@[vardimxobs],...,T@' + \
1526	'[vardimtsim]@[vardimtobs] is needed'
1527	quit(-1)
1528	else:
1529	print main +': couple of variable dimensions _______'
1530	vardims = {}
1531	ds = opts.vardims.split(',')
1532	for d in ds:
1533	dsecs = d.split('@')
1534	if len(dsecs) != 3:
1535	print errormsg
1536	print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
1537	print ' [DIM]@[vardimnsim]@[vardimnobs]'
1538	quit(-1)
1539	if dsecs[1] != 'None':
1540	vardims[dsecs[0]] = [dsecs[1], dsecs[2]]
1541	print ' ',dsecs[0],':',dsecs[1],',',dsecs[2]
1542
1543	if opts.obskind is None:
1544	print errormsg
1545	print ' ' + main + ': No kind of observations provided !!'
1546	quit(-1)
1547	else:
1548	obskind = opts.obskind
1549	if obskind == 'single-station':
1550	if opts.stloc is None:
1551	print errormsg
1552	print ' ' + main + ': No station location provided !!'
1553	quit(-1)
1554	else:
1555	stationdesc = [opts.stloc.split(',')[0].replace('\|',' '), \
1556	np.float(opts.stloc.split(',')[1]), np.float(opts.stloc.split(',')[2]),\
1557	np.float(opts.stloc.split(',')[3])]
1558
1559	if opts.fobs is None:
1560	print errormsg
1561	print ' ' + main + ': No observations file is provided!!'
1562	quit(-1)
1563	else:
1564	if not os.path.isfile(opts.fobs):
1565	print errormsg
1566	print ' ' + main + ": observations file '" + opts.fobs + "' does not exist !!"
1567	quit(-1)
1568
1569	if opts.fsim is None:
1570	print errormsg
1571	print ' ' + main + ': No simulation file is provided!!'
1572	quit(-1)
1573	else:
1574	if not os.path.isfile(opts.fsim):
1575	print errormsg
1576	print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
1577	quit(-1)
1578
1579	if opts.vars is None:
1580	print errormsg
1581	print ' ' + main + ': No list of couples of variables is provided!!'
1582	print ' a ',' list of values [varsim]@[varobs],... is needed'
1583	quit(-1)
1584	else:
1585	valvars = []
1586	vs = opts.vars.split(',')
1587	for v in vs:
1588	vsecs = v.split('@')
1589	if len(vsecs) < 2:
1590	print errormsg
1591	print ' ' + main + ': wrong number of values in:',vsecs, \
1592	' at least 2 are needed !!'
1593	print ' [varsim]@[varobs]@[[oper][val]]'
1594	quit(-1)
1595	if len(vsecs) > 2:
1596	if not searchInlist(simopers,vsecs[2]):
1597	print errormsg
1598	print main + ": operation on simulation values '" + vsecs[2] + \
1599	"' not ready !!"
1600	quit(-1)
1601
1602	valvars.append(vsecs)
1603
1604	# Openning observations trajectory
1605	##
1606	oobs = NetCDFFile(opts.fobs, 'r')
1607
1608	valdimobs = {}
1609	for dn in dims:
1610	print dn,':',dims[dn]
1611	if dims[dn][1] != 'None':
1612	if not oobs.dimensions.has_key(dims[dn][1]):
1613	print errormsg
1614	print ' ' + main + ": observations file does not have dimension '" + \
1615	dims[dn][1] + "' !!"
1616	quit(-1)
1617	if vardims[dn][1] != 'None':
1618	if not oobs.variables.has_key(vardims[dn][1]):
1619	print errormsg
1620	print ' ' + main + ": observations file does not have varibale " + \
1621	"dimension '" + vardims[dn][1] + "' !!"
1622	quit(-1)
1623	valdimobs[dn] = oobs.variables[vardims[dn][1]][:]
1624	else:
1625	if dn == 'X':
1626	valdimobs[dn] = stationdesc[1]
1627	elif dn == 'Y':
1628	valdimobs[dn] = stationdesc[2]
1629	elif dn == 'Z':
1630	valdimobs[dn] = stationdesc[3]
1631
1632	osim = NetCDFFile(opts.fsim, 'r')
1633
1634	valdimsim = {}
1635	for dn in dims:
1636	if dims[dn][0] != 'None':
1637	if not osim.dimensions.has_key(dims[dn][0]):
1638	print errormsg
1639	print ' ' + main + ": simulation file '" + opts.fsim + \
1640	"' does not have dimension '" + dims[dn][0] + "' !!"
1641	print ' it has: ',osim.dimensions
1642	quit(-1)
1643
1644	if not osim.variables.has_key(vardims[dn][0]) and \
1645	not searchInlist(varNOcheck,vardims[dn][0]):
1646	print errormsg
1647	print ' ' + main + ": simulation file '" + opts.fsim + \
1648	"' does not have varibale dimension '" + vardims[dn][0] + "' !!"
1649	print ' it has variables:',osim.variables
1650	quit(-1)
1651	if searchInlist(varNOcheck,vardims[dn][0]):
1652	valdimsim[dn] = compute_varNOcheck(osim, vardims[dn][0])
1653	else:
1654	valdimsim[dn] = osim.variables[vardims[dn][0]][:]
1655
1656	# General characteristics
1657	dimtobs = valdimobs['T'].shape[0]
1658	dimtsim = valdimsim['T'].shape[0]
1659
1660	print main +': observational time-steps:',dimtobs,'simulation:',dimtsim
1661
1662	notfound = np.zeros((dimtobs), dtype=int)
1663
1664	if obskind == 'multi-points':
1665	trajpos = np.zeros((2,dimt),dtype=int)
1666	for it in range(dimtobs):
1667	trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'], \
1668	[valdimobs['X'][it],valdimobss['Y'][it]])
1669	elif obskind == 'single-station':
1670	trajpos = None
1671	stationpos = np.zeros((2), dtype=int)
1672	if valdimsim.has_key('X') and valdimsim.has_key('Y'):
1673	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],[valdimobs['Y'], \
1674	valdimobs['X']])
1675	iid = 0
1676	for idn in osim.variables[vardims['X'][0]].dimensions:
1677	if idn == dims['X'][0]:
1678	stationpos[1] = stsimpos[iid]
1679	elif idn == dims['Y'][0]:
1680	stationpos[0] = stsimpos[iid]
1681
1682	iid = iid + 1
1683	print main + ': station point in simulation:', stationpos
1684	print ' station position:',valdimobs['X'],',',valdimobs['Y']
1685	print ' simulation coord.:',valdimsim['X'][tuple(stsimpos)],',', \
1686	valdimsim['Y'][tuple(stsimpos)]
1687	else:
1688	print main + ': validation with two time-series !!'
1689
1690	elif obskind == 'trajectory':
1691	if opts.trajf is not None:
1692	if not os.path.isfile(opts.fsim):
1693	print errormsg
1694	print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
1695	quit(-1)
1696	else:
1697	otrjf = NetCDFFile(opts.fsim, 'r')
1698	trajpos[0,:] = otrjf.variables['obssimtrj'][0]
1699	trajpos[1,:] = otrjf.variables['obssimtrj'][1]
1700	otrjf.close()
1701	else:
1702	if dims.has_key('Z'):
1703	trajpos = np.zeros((3,dimtobs),dtype=int)
1704	for it in range(dimtobs):
1705	if np.mod(it*100./dimtobs,10.) == 0.:
1706	print ' trajectory done: ',it*100./dimtobs,'%'
1707	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \
1708	[valdimobs['Y'][it],valdimobs['X'][it]])
1709	stationpos = np.zeros((2), dtype=int)
1710	iid = 0
1711	for idn in osim.variables[vardims['X'][0]].dimensions:
1712	if idn == dims['X'][0]:
1713	stationpos[1] = stsimpos[iid]
1714	elif idn == dims['Y'][0]:
1715	stationpos[0] = stsimpos[iid]
1716	iid = iid + 1
1717	if stationpos[0] == 0 and stationpos[1] == 0: notfound[it] = 1
1718
1719	trajpos[0,it] = stationpos[0]
1720	trajpos[1,it] = stationpos[1]
1721	# In the simulation 'Z' varies with time ... non-hydrostatic model! ;)
1722	# trajpos[2,it] = index_mat(valdimsim['Z'][it,:,stationpos[0], \
1723	# stationpos[1]], valdimobs['Z'][it])
1724	else:
1725	trajpos = np.zeros((2,dimtobs),dtype=int)
1726	for it in range(dimtobs):
1727	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \
1728	[valdimobs['Y'][it],valdimobss['X'][it]])
1729	stationpos = np.zeros((2), dtype=int)
1730	iid = 0
1731	for idn in osim.variables[vardims['X'][0]].dimensions:
1732	if idn == dims['X'][0]:
1733	stationpos[1] = stsimpos[iid]
1734	elif idn == dims['Y'][0]:
1735	stationpos[0] = stsimpos[iid]
1736	iid = iid + 1
1737	if stationpos[0] == 0 or stationpos[1] == 0: notfound[it] = 1
1738
1739	trajpos[0,it] = stationspos[0]
1740	trajpos[1,it] = stationspos[1]
1741
1742	print main + ': not found',np.sum(notfound),'points of the trajectory'
1743
1744	# Getting times
1745	tobj = oobs.variables[vardims['T'][1]]
1746	obstunits = tobj.getncattr('units')
1747	if vardims['T'][0] == 'WRFT':
1748	tsim = valdimsim['T'][:]
1749	simtunits = 'seconds since 1949-12-01 00:00:00'
1750	else:
1751	tsim = osim.variables[vardims['T'][0]][:]
1752	otsim = osim.variables[vardims['T'][0]]
1753	simtunits = otsim.getncattr('units')
1754
1755	simobstimes = coincident_CFtimes(tsim, obstunits, simtunits)
1756
1757	#
1758	## Looking for exact/near times
1759	#
1760
1761	# Exact Coincident times
1762	##
1763	exacttvalues0 = []
1764	for it in range(dimtsim):
1765	ot = 0
1766	for ito in range(ot,dimtobs-1):
1767	if valdimobs['T'][ito] == simobstimes[it]:
1768	ot = ito
1769	exacttvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito]])
1770
1771	exacttvalues = np.array(exacttvalues0, dtype=np.float)
1772	Nexactt = len(exacttvalues[:,0])
1773
1774	print main + ': found',Nexactt,'Temporal same values in simulation and observations'
1775
1776	# Sim Closest times
1777	##
1778	Nsimt = 0
1779	closesttvalues0 = []
1780	closesttvalues0st = []
1781	tsiminit = 0
1782	tsimend = 0
1783
1784	dtsim = simobstimes[1] - simobstimes[0]
1785
1786	for it in range(dimtsim):
1787	ot = 0
1788	for ito in range(ot,dimtobs-1):
1789	if np.abs(valdimobs['T'][ito] - simobstimes[it]) <= dtsim/2.:
1790	ot = ito
1791	tdist = simobstimes[it] - valdimobs['T'][ito]
1792	closesttvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito], \
1793	tdist])
1794	Nsimt = Nsimt + 1
1795	if tsiminit == 0: tsiminit=simobstimes[it]
1796	tsimend = simobstimes[it]
1797
1798	closesttvalues = np.array(closesttvalues0, dtype=np.float)
1799
1800	Nclosest = len(closesttvalues[:,0])
1801	print main + ': found',Nclosest,'Simulation time-values closest to observations'
1802
1803	if Nclosest == 0:
1804	print warnmsg
1805	print main + ': no cclosest times found !!'
1806	print ' stopping it'
1807	quit(-1)
1808
1809	# Sim Coincident times
1810	##
1811	Nsimt = 0
1812	coindtvalues0 = []
1813	coindtvalues0st = []
1814	tsiminit = 0
1815	tsimend = 0
1816
1817	for it in range(dimtsim):
1818	ot = 0
1819	for ito in range(ot,dimtobs-1):
1820	if valdimobs['T'][ito] < simobstimes[it] and valdimobs['T'][ito+1] > \
1821	simobstimes[it]:
1822	ot = ito
1823	tdist = simobstimes[it] - valdimobs['T'][ito]
1824	coindtvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito], \
1825	tdist])
1826	Nsimt = Nsimt + 1
1827	if tsiminit == 0: tsiminit=simobstimes[it]
1828	tsimend = simobstimes[it]
1829	elif simobstimes[it] > valdimobs['T'][ito+1]:
1830	coindtvalues0st.append([Nsimt, ito, valdimobs['T'][ito],tsimend-tsiminit])
1831
1832	coindtvalues = np.array(coindtvalues0, dtype=np.float)
1833	coindtvaluesst = np.array(coindtvalues0st, dtype=np.float)
1834
1835	Ncoindt = len(coindtvalues[:,0])
1836	print main + ': found',Ncoindt,'Simulation time-interval (within consecutive ' + \
1837	'observed times) coincident times between simulation and observations'
1838
1839	if Ncoindt == 0:
1840	print warnmsg
1841	print main + ': no coincident times found !!'
1842	print ' stopping it'
1843	quit(-1)
1844
1845	# Validating
1846	##
1847
1848	onewnc = NetCDFFile(ofile, 'w')
1849
1850	# Dimensions
1851	for kst in range(Nstsim):
1852	newdim = onewnc.createDimension(prestdescsim[kst] + 'time',None)
1853	if stdescsim[kst] != 'E':
1854	newdim = onewnc.createDimension(prestdescsim[kst] + 'obstime',None)
1855
1856	newdim = onewnc.createDimension('bnds',2)
1857	newdim = onewnc.createDimension('couple',2)
1858	newdim = onewnc.createDimension('StrLength',StringLength)
1859	newdim = onewnc.createDimension('xaround',Ngrid*2+1)
1860	newdim = onewnc.createDimension('yaround',Ngrid*2+1)
1861	newdim = onewnc.createDimension('gstats',13)
1862	newdim = onewnc.createDimension('stats',5)
1863	newdim = onewnc.createDimension('tstats',6)
1864	newdim = onewnc.createDimension('Nstsim', 3)
1865
1866	# Variable dimensions
1867	##
1868	newvar = onewnc.createVariable('couple', 'c', ('couple','StrLength'))
1869	basicvardef(newvar, 'couple', 'couples of values', '-')
1870	writing_str_nc(newvar, ['sim','obs'], StringLength)
1871
1872	newvar = onewnc.createVariable('statistics', 'c', ('stats','StrLength'))
1873	basicvardef(newvar, 'statistics', 'statitics from values', '-')
1874	writing_str_nc(newvar, statsn, StringLength)
1875
1876	newvar = onewnc.createVariable('gstatistics', 'c', ('gstats','StrLength'))
1877	basicvardef(newvar, 'gstatistics', 'global statitics from values', '-')
1878	writing_str_nc(newvar, gstatsn, StringLength)
1879
1880	newvar = onewnc.createVariable('tstatistics', 'c', ('tstats','StrLength'))
1881	basicvardef(newvar, 'tstatistics', 'statitics from values along time', '-')
1882	writing_str_nc(newvar, ostatsn, StringLength)
1883
1884	newvar = onewnc.createVariable('ksimstatistics', 'c', ('Nstsim','StrLength'))
1885	basicvardef(newvar, 'ksimstatistics', 'kind of simulated statitics', '-')
1886	writing_str_nc(newvar, Lstdescsim, StringLength)
1887
1888	if obskind == 'trajectory':
1889	if dims.has_key('Z'):
1890	newdim = onewnc.createDimension('trj',3)
1891	else:
1892	newdim = onewnc.createDimension('trj',2)
1893
1894	newvar = onewnc.createVariable('obssimtrj','i',('obstime','trj'))
1895	basicvardef(newvar, 'obssimtrj', 'trajectory on the simulation grid', '-')
1896	newvar[:] = trajpos.transpose()
1897
1898	if dims.has_key('Z'):
1899	newdim = onewnc.createDimension('simtrj',4)
1900	else:
1901	newdim = onewnc.createDimension('simtrj',3)
1902
1903	Nvars = len(valvars)
1904	for ivar in range(Nvars):
1905	simobsvalues = []
1906
1907	varsimobs = valvars[ivar][0] + '_' + valvars[ivar][1]
1908	print ' ' + varsimobs + '... .. .'
1909
1910	if not oobs.variables.has_key(valvars[ivar][1]):
1911	print errormsg
1912	print ' ' + main + ": observations file has not '" + valvars[ivar][1] + \
1913	"' !!"
1914	quit(-1)
1915
1916	if not osim.variables.has_key(valvars[ivar][0]):
1917	if not searchInlist(varNOcheck, valvars[ivar][0]):
1918	print errormsg
1919	print ' ' + main + ": simulation file has not '" + valvars[ivar][0] + \
1920	"' !!"
1921	quit(-1)
1922	else:
1923	ovsim = compute_varNOcheck(osim, valvars[ivar][0])
1924	else:
1925	ovsim = osim.variables[valvars[ivar][0]]
1926
1927	for idn in ovsim.dimensions:
1928	if not searchInlist(simdims.values(),idn):
1929	print errormsg
1930	print ' ' + main + ": dimension '" + idn + "' of variable '" + \
1931	valvars[ivar][0] + "' not provided as reference coordinate [X,Y,Z,T] !!"
1932	quit(-1)
1933
1934	ovobs = oobs.variables[valvars[ivar][1]]
1935	if searchInlist(ovobs.ncattrs(),'_FillValue'):
1936	oFillValue = ovobs.getncattr('_FillValue')
1937	else:
1938	oFillValue = None
1939
1940	for kst in range(Nstsim):
1941	timedn = prestdescsim[kst] + 'time'
1942
1943	if stdescsim[kst] == 'E':
1944	# Observed and simualted exact times
1945	simobsvalues, simobsSvalues, simobsTtvalues, trjsim = \
1946	getting_ValidationValues(obskind, Nexactt, dims, trajpos, ovsim, \
1947	ovobs, exacttvalues, oFillValue, Ngrid, 'instantaneous')
1948
1949	vname = prestdescsim[kst] + 'time'
1950
1951	newvar = onewnc.createVariable(vname,'f8', (timedn))
1952	basicvardef(newvar, vname, 'exact coincident time observations and ' + \
1953	'simulation', obstunits)
1954	set_attribute(newvar, 'calendar', 'standard')
1955	newvar[:] = exacttvalues[:,3]
1956
1957	elif == stdescsim[kst] == 'C':
1958	# Simualted closest to Observed times
1959	simobsvalues, simobsSvalues, simobsTtvalues, trjsim = \
1960	getting_ValidationValues(obskind, Nclosest, dims, trajpos, ovsim, \
1961	ovobs, closesttvalues, oFillValue, Ngrid, 'instantaneous')
1962
1963	vname = prestdescsim[kst] + 'time'
1964
1965	newvar = onewnc.createVariable(vname,'f8', (timedn))
1966	basicvardef(newvar, vname, 'time simulations closest to observed values',\
1967	obstunits )
1968	set_attribue(newvar, 'calendar', 'standard')
1969	newvar[:] = closesttvalues[:,2]
1970
1971	vname = (prestdescsim[kst] + 'obstime'
1972	newvar = onewnc.createVariable(vname,'f8', (vname))
1973	basicvardef(newvar, vname, 'closest time observations', obstunits)
1974	set_attribute(newvar, 'calendar', 'standard')
1975	newvar[:] = closesttvalues[:,3]
1976
1977	elif == stdescsim[kst] == 'B':
1978	# Observed values temporally around coincident times
1979	simobsvalues, simobsSvalues, simobsTtvalues, trjsim = \
1980	getting_ValidationValues(obskind, Ncoindt, dims, trajpos, ovsim, ovobs,\
1981	coindtvalues, oFillValue, Ngrid, 'tbackwardSmean')
1982
1983	vname = prestdescsim[kst] + 'time'
1984	newvar = onewnc.createVariable(vname,'f8', (timedn))
1985	basicvardef(newvar, vname, 'simulation time between observations', \
1986	obstunits )
1987	set_attribute(newvar, 'calendar', 'standard')
1988	set_attribute(newvar, 'bounds', 'time_bnds')
1989	newvar[:] = simobsTtvalues[:,1]
1990
1991	vname = prestdescsim[kst] + 'obstime'
1992	newvar = onewnc.createVariable(vname,'f8', (vname))
1993	basicvardef(newvar, vname, 'observed between time', obstunits )
1994	set_attribute(newvar, 'calendar', 'standard')
1995	newvar[:] = np.unique(coindtvalues[:,2])
1996
1997	# Re-arranging values...
1998	arrayvals = np.array(simobsvalues)
1999	if len(valvars[ivar]) > 2:
2000	const=np.float(valvars[ivar][3])
2001	if valvars[ivar][2] == 'sumc':
2002	simobsSvalues = simobsSvalues + const
2003	arrayvals[:,0] = arrayvals[:,0] + const
2004	elif valvars[ivar][2] == 'subc':
2005	simobsSvalues = simobsSvalues - const
2006	arrayvals[:,0] = arrayvals[:,0] - const
2007	elif valvars[ivar][2] == 'mulc':
2008	simobsSvalues = simobsSvalues * const
2009	arrayvals[:,0] = arrayvals[:,0] * const
2010	elif valvars[ivar][2] == 'divc':
2011	simobsSvalues = simobsSvalues / const
2012	arrayvals[:,0] = arrayvals[:,0] / const
2013	else:
2014	print errormsg
2015	print ' ' + fname + ": operation '"+valvars[ivar][2]+"' not ready!!"
2016	quit(-1)
2017
2018	# statisics sim
2019	simstats = np.zeros((Nstsim,5), dtype=np.float)
2020	simstats[kst,0] = np.min(arrayvals[:,0])
2021	simstats[kst,1] = np.max(arrayvals[:,0])
2022	simstats[kst,2] = np.mean(arrayvals[:,0])
2023	simstats[kst,3] = np.mean(arrayvals[:,0]*arrayvals[:,0])
2024	simstats[kst,4] = np.sqrt(simstats[kst,3] - simstats[kst,2]*simstats[kst,2])
2025
2026	# statisics obs
2027	# Masking 'nan'
2028	obsmask0 = np.where(arrayvals[:,1] != arrayvals[:,1], fillValueF, \
2029	arrayvals[:,1])
2030
2031	obsmask = ma.masked_equal(obsmask0, fillValueF)
2032	obsmask2 = obsmask*obsmask
2033
2034	obsstats = np.zeros((Nstsim,5), dtype=np.float)
2035	obsstats[kst,0] = obsmask.min()
2036	obsstats[kst,1] = obsmask.max()
2037	obsstats[kst,2] = obsmask.mean()
2038	obsstats[kst,3] = obsmask2.mean()
2039	obsstats[kst,4] = np.sqrt(obsstats[kst,3] - obsstats[kst,2]*obsstats[kst,2])
2040
2041	# Statistics sim-obs
2042	simobsstats = np.zeros((Nstsim,13), dtype=np.float)
2043	diffvals = np.zeros((Ncoindt), dtype=np.float)
2044
2045	diffvals = arrayvals[:,0] - obsmask
2046
2047	diff2vals = diffvals*diffvals
2048	sumdiff = diffvals.sum()
2049	sumdiff2 = diff2vals.sum()
2050
2051	simobsstats[kst,0] = simstats[kst,0] - obsstats[kst,0]
2052	simobsstats[kst,1] = np.mean(arrayvals[:,0]*obsmask)
2053	simobsstats[kst,2] = diffvals.min()
2054	simobsstats[kst,3] = diffvals.max()
2055	simobsstats[kst,4] = diffvals.mean()
2056	simobsstats[kst,5] = np.abs(diffvals).mean()
2057	simobsstats[kst,6] = np.sqrt(diff2vals.mean())
2058	simobsstats[kst,7], simobsstats[kst,8] = sts.mstats.pearsonr(arrayvals[:,0], \
2059	obsmask)
2060	# From:
2061	#Willmott, C. J. 1981. 'On the validation of models. Physical Geography', 2, 184-194
2062	#Willmott, C. J. (1984). 'On the evaluation of model performance in physical
2063	# geography'. Spatial Statistics and Models, G. L. Gaile and C. J. Willmott, eds.,
2064	# 443-460
2065	#Willmott, C. J., S. G. Ackleson, R. E. Davis, J. J. Feddema, K. M. Klink, D. R.
2066	# Legates, J. O'Donnell, and C. M. Rowe (1985), 'Statistics for the Evaluation and
2067	# Comparison of Models', J. Geophys. Res., 90(C5), 8995-9005
2068	#Legates, D. R., and G. J. McCabe Jr. (1999), 'Evaluating the Use of "Goodness-of-Fit"
2069	# Measures in Hydrologic and Hydroclimatic Model Validation', Water Resour. Res.,
2070	# 35(1), 233-241
2071	#
2072	# Deviation of residuals (SDR)
2073	simobsstats[kst,9] = np.mean(np.sqrt(np.abs((diffvals-simobsstats[kst,0])* \
2074	(diffvals-obsstats[kst,0]))))
2075	# Index of Efficiency (IOE)
2076	obsbias2series = (obsmask - obsstats[kst,0])*(obsmask - obsstats[kst,0])
2077	sumobsbias2series = obsbias2series.sum()
2078
2079	simobsstats[kst,10] = 1. - sumdiff2/(sumobsbias2series)
2080	# Index of Agreement (IOA)
2081	simbias2series = arrayvals[:,0] - obsstats[kst,0]
2082	obsbias2series = obsmask - obsstats[kst,0]
2083
2084	abssimbias2series = np.abs(simbias2series)
2085	absobsbias2series = np.where(obsbias2series<0, -obsbias2series, \
2086	obsbias2series)
2087
2088	abssimobsbias2series = (abssimbias2series+absobsbias2series)*( \
2089	abssimbias2series + absobsbias2series)
2090
2091	simobsstats[kst,11] = 1. - sumdiff2/(abssimobsbias2series.sum())
2092	# Fractional Mean Bias (FMB)
2093	simobsstats[ksy,12]=(simstats[ksy,0]-obsstats[kst,0])/(0.5*(simstats[kst,0] +\
2094	obsstats[kst,0]))
2095
2096	# Statistics around sim values
2097	aroundstats = np.zeros((5,Ncoindt), dtype=np.float)
2098	for it in range(Ncoindt):
2099	aroundstats[0,it] = np.min(simobsSvalues[it,])
2100	aroundstats[1,it] = np.max(simobsSvalues[it,])
2101	aroundstats[2,it] = np.mean(simobsSvalues[it,])
2102	aroundstats[3,it] = np.mean(simobsSvalues[it,]*simobsSvalues[it,])
2103	aroundstats[4,it] = np.sqrt(aroundstats[3,it] - aroundstats[2,it]* \
2104	aroundstats[2,it])
2105
2106	# sim Values to netCDF
2107	newvar = onewnc.createVariable(valvars[ivar][0] + '_' + stdescsim[kst] + \
2108	'sim', 'f', (timedn), fill_value=fillValueF)
2109	descvar = prestdescsim[kst] + ' time simulated: ' + valvars[ivar][0]
2110	basicvardef(newvar, valvars[ivar][0], descvar, ovobs.getncattr('units'))
2111	newvar[:] = arrayvals[:,0]
2112
2113	# obs Values to netCDF
2114	newvar = onewnc.createVariable(valvars[ivar][1] + '_' + stdescsim[kst] + \
2115	'obs', 'f', (timedn), fill_value=fillValueF)
2116	descvar = prestdescsim[kst] + ' time simulated: ' + valvars[ivar][1]
2117	basicvardef(newvar, valvars[ivar][1], descvar, ovobs.getncattr('units'))
2118	newvar[:] = arrayvals[:,1]
2119
2120	# Around values
2121	if not onewnc.variables.has_key(valvars[ivar][0] + 'around'):
2122	vname = prestdescsim[kst] + valvars[ivar][0] + 'around'
2123	else:
2124	vname = prestdescsim[kst] + valvars[ivar][0] + 'Around'
2125
2126	if dims.has_key('Z'):
2127	if not onewnc.dimensions.has_key('zaround'):
2128	newdim = onewnc.createDimension('zaround',Ngrid*2+1)
2129	newvar = onewnc.createVariable(vname, 'f', (timedn,'zaround', \
2130	'yaround','xaround'), fill_value=fillValueF)
2131	else:
2132	newvar = onewnc.createVariable(vname, 'f', (timedn,'yaround','xaround'), \
2133	fill_value=fillValueF)
2134
2135	descvar = prestdescsim[kst] + 'around simulated values +/- grid values: ' + \
2136	valvars[ivar][0]
2137	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2138	newvar[:] = simobsSvalues
2139
2140
2141	# around sim Statistics between
2142	if not searchInlist(onewnc.variables,prestdescsim[kst] + valvars[ivar][0] + \
2143	'staround'):
2144	vname = prestdescsim[kst] + valvars[ivar][0] + 'staround'
2145	else:
2146	vname = prestdescsim[kst] + valvars[ivar][0] + 'Staround'
2147
2148	newvar = onewnc.createVariable(vname, 'f', (timedn,'stats'), \
2149	fill_value=fillValueF)
2150	descvar = prestdescsim[kst] + ' around (' + str(Ngrid) + ', ' + str(Ngrid) +\
2151	') simulated statisitcs: ' + valvars[ivar][0]
2152	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2153	set_attribute(newvar, 'cell_methods', 'time_bnds')
2154	newvar[:] = aroundstats.transpose()
2155
2156	if stdescsim[kst] == 'B':
2157	if not searchInlist(onewnc.variables, 'time_bnds'):
2158	newvar = onewnc.createVariable('time_bnds','f8',(imedn,'bnds'))
2159	basicvardef(newvar, 'time_bnds', timedn, obstunits )
2160	set_attribute(newvar, 'calendar', 'standard')
2161	newvar[:] = simobsTtvalues
2162
2163	# around obs Statistics between
2164	if not searchInlist(onewnc.variables, prestdescsim[kst] + valvars[ivar][1] + \
2165	'staround'):
2166	vname = prestdescsim[kst] + valvars[ivar][1] + 'staround'
2167	else:
2168	vname = prestdescsim[kst] + valvars[ivar][1] + 'Staround'
2169
2170	newvar = onewnc.createVariable(vname, 'f', (timedn,'tstats'), \
2171	fill_value=fillValueF)
2172	descvar = prestdescsim[kst] + ' around temporal observed statisitcs: ' + \
2173	valvars[ivar][1]
2174	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2175	set_attribute(newvar, 'cell_methods', 'time_bnds')
2176	newvar[:] = aroundostats.transpose()
2177
2178	# sim Statistics
2179	if not searchInlist(onewnc.variables,valvars[ivar][0] + 'stsim'):
2180	vname = valvars[ivar][0] + 'stsim'
2181	else:
2182	vname = valvars[ivar][0] + 'stSim'
2183
2184	newvar = onewnc.createVariable(vname, 'f', ('Nstsim', 'stats'), \
2185	fill_value=fillValueF)
2186	descvar = 'simulated statisitcs: ' + valvars[ivar][0]
2187	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2188	newvar[:] = simstats
2189
2190	# obs Statistics
2191	if not searchInlist(onewnc.variables,valvars[ivar][1] + 'stobs'):
2192	vname = valvars[ivar][1] + 'stobs'
2193	else:
2194	vname = valvars[ivar][1] + 'stObs'
2195
2196	newvar = onewnc.createVariable(vname, 'f', ('Nstsim', 'stats'), \
2197	fill_value=fillValueF)
2198	descvar = 'observed statisitcs: ' + valvars[ivar][1]
2199	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2200	newvar[:] = obsstats
2201
2202	# sim-obs Statistics
2203	if not searchInlist(onewnc.variables,varsimobs + 'st'):
2204	vname = varsimobs + 'st'
2205	else:
2206	vname = varSimobs + 'st'
2207
2208	newvar = onewnc.createVariable(vname, 'f', ('Nstsim', 'gstats'), \
2209	fill_value=fillValueF)
2210	descvar = 'simulated-observed tatisitcs: ' + varsimobs
2211	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2212	newvar[:] = simobsstats
2213
2214	onewnc.sync()
2215
2216	if trjsim is not None:
2217	newvar = onewnc.createVariable('simtrj','i',('betweentime','simtrj'))
2218	basicvardef(newvar,'simtrj','coordinates [X,Y,Z,T] of the coincident ' + \
2219	'trajectory in sim', obstunits)
2220	newvar[:] = trjsim.transpose()
2221
2222	# Adding three variables with the station location, longitude, latitude and height
2223	if obskind == 'single-station':
2224	adding_station_desc(onewnc,stationdesc)
2225
2226	# Global attributes
2227	##
2228	set_attribute(onewnc,'author_nc','Lluis Fita')
2229	set_attribute(onewnc,'institution_nc','Laboratoire de Meteorology Dynamique, ' + \
2230	'LMD-Jussieu, UPMC, Paris')
2231	set_attribute(onewnc,'country_nc','France')
2232	set_attribute(onewnc,'script_nc',main)
2233	set_attribute(onewnc,'version_script',version)
2234	set_attribute(onewnc,'information', \
2235	'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html')
2236	set_attribute(onewnc,'simfile',opts.fsim)
2237	set_attribute(onewnc,'obsfile',opts.fobs)
2238
2239	onewnc.sync()
2240	onewnc.close()
2241
2242	print main + ": successfull writting of '" + ofile + "' !!"

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

Download in other formats: