Context Navigation

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

Last change on this file since 529 was 529, checked in by lfita, 9 years ago
Fixing minor things problem on `betweentim' 0 values at it=0
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 = 'validation_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, dt, 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	dt= 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	dtf = dt
1174	elif kvals == 'tbackwardSmean':
1175	print ' ' + fname + ':',kvals,'!!'
1176	uniqt = np.unique(tvalues[:,3])
1177	dtf = len(uniqt)
1178	print ' initially we got',dt,'values which will become',dtf
1179	postf = {}
1180	for it in range(dtf):
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	dtf = len(uniqt)
1195	print ' initially we got',dt,'values which will become',dtf
1196	print ' ==== NOT READY === '
1197	quit(-1)
1198	else:
1199	print errormsg
1200	print ' ' + fname + ": kind of values '" + kvals + "' not ready!!"
1201	quit(-1)
1202
1203	# Simulated values spatially around
1204	if ds.has_key('Z'):
1205	soSvalues = np.zeros((dt, Ng2+1, Ng2+1, Ng*2+1), dtype = np.float)
1206	if okind == 'trajectory':
1207	trjs = np.zeros((4,dt), dtype=int)
1208	else:
1209	trjs = None
1210	else:
1211	soSvalues = np.zeros((dt, Ng2+1, Ng2+1), dtype = np.float)
1212	if okind == 'trajectory':
1213	trjs = np.zeros((3,dt), dtype=int)
1214	else:
1215	trjs = None
1216
1217	if okind == 'single-station':
1218	soTtvalues = np.zeros((dt,2), dtype=np.float)
1219	else:
1220	None
1221
1222	if okind == 'multi-points':
1223	for it in range(dt):
1224	slicev = ds['X'][0] + ':' + str(trjpos[2,it]) + '\|' + ds['Y'][0] + \
1225	':' + str(trjpos[1,it]) + '\|' + ds['T'][0]+ ':' + str(tvalues[it][0])
1226	slicevar, dimslice = slice_variable(ovs, slicev)
1227	sovalues.append([ slicevar, ovo[tvalues[it][1]]])
1228	slicev = ds['X'][0] + ':' + str(trjpos[2,it]-Ng) + '@' + \
1229	str(trjpos[2,it]+Ng) + '\|' + ds['Y'][0] + ':' + \
1230	str(trjpos[1,it]-Ng) + '@' + str(trjpos[1,it]+Ng) + '\|' + \
1231	ds['T'][0]+':'+str(tvalues[it][0])
1232	slicevar, dimslice = slice_variable(ovs, slicev)
1233	soSvalues[it,:,:] = slicevar
1234
1235	elif okind == 'single-station':
1236	for it in range(dt):
1237	ito = int(tvalues[it,1])
1238	if valdimsim.has_key('X') and valdimsim.has_key('Y'):
1239	slicev = ds['X'][0] + ':' + str(stationpos[1]) + '\|' + \
1240	ds['Y'][0] + ':' + str(stationpos[0]) + '\|' + \
1241	ds['T'][0] + ':' + str(int(tvalues[it][0]))
1242	else:
1243	slicev = ds['T'][0] + ':' + str(int(tvalues[it][0]))
1244	slicevar, dimslice = slice_variable(ovs, slicev)
1245	if ovo[int(ito)] == oFill or ovo[int(ito)] == '--':
1246	sovalues.append([ slicevar, fillValueF])
1247	# elif ovo[int(ito)] != ovo[int(ito)]:
1248	# sovalues.append([ slicevar, fillValueF])
1249	else:
1250	sovalues.append([ slicevar, ovo[int(ito)]])
1251	if valdimsim.has_key('X') and valdimsim.has_key('Y'):
1252	slicev = ds['X'][0] + ':' + str(stationpos[1]-Ng) + '@' + \
1253	str(stationpos[1]+Ng+1) + '\|' + ds['Y'][0] + ':' + \
1254	str(stationpos[0]-Ng) + '@' + str(stationpos[0]+Ng+1) + '\|' + \
1255	ds['T'][0] + ':' + str(int(tvalues[it,0]))
1256	else:
1257	slicev = ds['T'][0] + ':' + str(int(tvalues[it][0]))
1258	slicevar, dimslice = slice_variable(ovs, slicev)
1259	soSvalues[it,:,:] = slicevar
1260
1261	if it == 0:
1262	itoi = 0
1263	itof = int(tvalues[it,1]) / 2
1264	elif it == dt-1:
1265	itoi = int( (ito + int(tvalues[it-1,1])) / 2)
1266	itof = int(tvalues[it,1])
1267	else:
1268	itod = int( (ito - int(tvalues[it-1,1])) / 2 )
1269	itoi = ito - itod
1270	itod = int( (int(tvalues[it+1,1]) - ito) / 2 )
1271	itof = ito + itod
1272
1273	soTtvalues[it,0] = valdimobs['T'][itoi]
1274	soTtvalues[it,1] = valdimobs['T'][itof]
1275
1276	elif okind == 'trajectory':
1277	if ds.has_key('Z'):
1278	for it in range(dt):
1279	ito = int(tvalues[it,1])
1280	if notfound[ito] == 0:
1281	trjpos[2,ito] = index_mat(valdimsim['Z'][tvalues[it,0],:, \
1282	trjpos[1,ito],trjpos[0,ito]], valdimobs['Z'][ito])
1283	slicev = ds['X'][0]+':'+str(trjpos[0,ito]) + '\|' + \
1284	ds['Y'][0]+':'+str(trjpos[1,ito]) + '\|' + \
1285	ds['Z'][0]+':'+str(trjpos[2,ito]) + '\|' + \
1286	ds['T'][0]+':'+str(int(tvalues[it,0]))
1287	slicevar, dimslice = slice_variable(ovs, slicev)
1288	sovalues.append([ slicevar, ovo[int(ito)]])
1289	minx = np.max([trjpos[0,ito]-Ng,0])
1290	maxx = np.min([trjpos[0,ito]+Ng+1,ovs.shape[3]])
1291	miny = np.max([trjpos[1,ito]-Ng,0])
1292	maxy = np.min([trjpos[1,ito]+Ng+1,ovs.shape[2]])
1293	minz = np.max([trjpos[2,ito]-Ng,0])
1294	maxz = np.min([trjpos[2,ito]+Ng+1,ovs.shape[1]])
1295
1296	slicev = ds['X'][0] + ':' + str(minx) + '@' + str(maxx) + '\|' + \
1297	ds['Y'][0] + ':' + str(miny) + '@' + str(maxy) + '\|' + \
1298	ds['Z'][0] + ':' + str(minz) + '@' + str(maxz) + '\|' + \
1299	ds['T'][0] + ':' + str(int(tvalues[it,0]))
1300	slicevar, dimslice = slice_variable(ovs, slicev)
1301
1302	sliceS = []
1303	sliceS.append(it)
1304	sliceS.append(slice(0,maxz-minz))
1305	sliceS.append(slice(0,maxy-miny))
1306	sliceS.append(slice(0,maxx-minx))
1307
1308	soSvalues[tuple(sliceS)] = slicevar
1309	if ivar == 0:
1310	trjs[0,it] = trjpos[0,ito]
1311	trjs[1,it] = trjpos[1,ito]
1312	trjs[2,it] = trjpos[2,ito]
1313	trjs[3,it] = tvalues[it,0]
1314	else:
1315	sovalues.append([fillValueF, fillValueF])
1316	soSvalues[it,:,:,:]= np.ones((Ng2+1,Ng2+1,Ng*2+1), \
1317	dtype = np.float)*fillValueF
1318	# 2D trajectory
1319	else:
1320	for it in range(dt):
1321	if notfound[it] == 0:
1322	ito = tvalues[it,1]
1323	slicev = ds['X'][0]+':'+str(trjpos[2,ito]) + '\|' + \
1324	ds['Y'][0]+':'+str(trjpos[1,ito]) + '\|' + \
1325	ds['T'][0]+':'+str(tvalues[ito,0])
1326	slicevar, dimslice = slice_variable(ovs, slicev)
1327	sovalues.append([ slicevar, ovo[tvalues[it,1]]])
1328	slicev = ds['X'][0] + ':' + str(trjpos[0,it]-Ng) + '@' + \
1329	str(trjpos[0,it]+Ng) + '\|' + ds['Y'][0] + ':' + \
1330	str(trjpos[1,it]-Ng) + '@' + str(trjpos[1,it]+Ng) + \
1331	'\|' + ds['T'][0] + ':' + str(tvalues[it,0])
1332	slicevar, dimslice = slice_variable(ovs, slicev)
1333	soSvalues[it,:,:] = slicevar
1334	else:
1335	sovalues.append([fillValue, fillValue])
1336	soSvalues[it,:,:] = np.ones((Ng2+1,Ng2+1), \
1337	dtype = np.float)*fillValueF
1338	print sovalues[varsimobs][:][it]
1339	else:
1340	print errormsg
1341	print ' ' + fname + ": observatino kind '" + okind + "' not ready!!"
1342	quit(-1)
1343
1344	# Re-arranging final values
1345	##
1346	if kvals == 'instantaneous':
1347	return sovalues, soSvalues, soTtvalues, trjs
1348
1349	elif kvals == 'tbackwardSmean':
1350	fsovalues = []
1351	if ds.has_key('Z'):
1352	fsoSvalues = np.zeros((dtf, Ng2+1, Ng2+1, Ng*2+1), dtype = np.float)
1353	if okind == 'trajectory':
1354	ftrjs = np.zeros((4,dtf), dtype=int)
1355	else:
1356	ftrjs = None
1357	else:
1358	fsoSvalues = np.zeros((dtf, Ng2+1, Ng2+1), dtype = np.float)
1359	if okind == 'trajectory':
1360	ftrjs = np.zeros((3,dtf), dtype=int)
1361	else:
1362	ftrjs = None
1363
1364	if okind == 'single-station':
1365	fsoTtvalues = np.ones((dtf,2), dtype=np.float)*fillValueF
1366	else:
1367	None
1368
1369	for it in range(1,dtf):
1370	tv = uniqt[it]
1371	intv = postf[tv]
1372
1373	# Temporal statistics
1374	sovs = np.array(sovalues[intv[0]:intv[1]])[:,0]
1375	minv = np.min(sovs)
1376	maxv = np.max(sovs)
1377	meanv = np.mean(sovs)
1378	stdv = np.std(sovs)
1379
1380	fsovalues.append([meanv, np.array(sovalues[intv[0]:intv[1]])[0,1], minv, \
1381	maxv, stdv])
1382	if ds.has_key('Z'):
1383	if okind == 'trajectory':
1384	for ip in range(4):
1385	ftrjs[ip,it] = np.mean(trjs[ip,intv[0]:intv[1]])
1386	for iz in range(2*Ng+1):
1387	for iy in range(2*Ng+1):
1388	for ix in range(2*Ng+1):
1389	fsoSvalues[it,iz,iy,ix] = np.mean(soSvalues[intv[0]: \
1390	intv[1],iz,iy,ix])
1391	else:
1392	if okind == 'trajectory':
1393	for ip in range(3):
1394	ftrjs[ip,it] = np.mean(trjs[ip,intv[0]:intv[1]])
1395	for iy in range(2*Ng+1):
1396	for ix in range(2*Ng+1):
1397	fsoSvalues[it,iy,ix] = np.mean(soSvalues[intv[0]:intv[1], \
1398	iy,ix])
1399	fsoTtvalues[it,0] = soTtvalues[intv[0],0]
1400	fsoTtvalues[it,1] = soTtvalues[intv[1],0]
1401
1402	return fsovalues, fsoSvalues, fsoTtvalues, ftrjs
1403
1404	elif kvals == 'tbackwardOmean':
1405	print ' ' + fname + ':',kvals,'!!'
1406	uniqt = np.unique(tvalues[:,2])
1407	dtf = len(uniqt)
1408	print ' initially we got',dt,'values which will become',dtf
1409
1410	return
1411
1412
1413	####### ###### ##### #### ### ## #
1414
1415	strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " + \
1416	" and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')"
1417
1418	kindobs=['multi-points', 'single-station', 'trajectory']
1419	strkObs="kind of observations; 'multi-points': multiple individual punctual obs " + \
1420	"(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
1421	"'trajectory': following a trajectory"
1422	simopers = ['sumc','subc','mulc','divc']
1423	opersinf = 'sumc,[constant]: add [constant] to variables values; subc,[constant]: '+ \
1424	'substract [constant] to variables values; mulc,[constant]: multipy by ' + \
1425	'[constant] to variables values; divc,[constant]: divide by [constant] to ' + \
1426	'variables values'
1427	varNOcheck = ['WRFdens', 'WRFght', 'WRFp', 'WRFrh', 'TSrhs', 'WRFrhs', 'WRFT', \
1428	'WRFt', 'WRFtd', 'WRFwd', 'TSwds', 'WRFwds', 'WRFws', 'TSwss', 'WRFwss', 'WRFz']
1429	varNOcheckinf = "'WRFdens': air density from WRF variables; " + \
1430	"'WRFght': geopotentiali height from WRF variables; " + \
1431	"'WRFp': pressure from WRF variables; " + \
1432	"'WRFrh': relative humidty fom WRF variables; " + \
1433	"'WRFrhs': surface relative humidity from WRF variables; " + \
1434	"'WRFT': CF-time from WRF variables; " + \
1435	"'WRFt': temperature from WRF variables; " + \
1436	"'WRFtd': dew-point temperature from WRF variables; " + \
1437	"'WRFwd': wind direction from WRF variables; " + \
1438	"'WRFwds': surface wind direction from WRF variables; " + \
1439	"'WRFws': wind speed from WRF variables; " + \
1440	"'WRFwss': surface wind speed from WRF variables; " + \
1441	"'WRFz': height from WRF variables"
1442
1443	dimshelp = "[DIM]@[simdim]@[obsdim] ',' list of couples of dimensions names from " + \
1444	"each source ([DIM]='X','Y','Z','T'; None, no value)"
1445	vardimshelp = "[DIM]@[simvardim]@[obsvardim] ',' list of couples of variables " + \
1446	"names with dimensions values from each source ([DIM]='X','Y','Z','T'; None, " + \
1447	"no value, WRFdiagnosted variables also available: " + varNOcheckinf + ")"
1448	varshelp="[simvar]@[obsvar]@[[oper]@[val]] ',' list of couples of variables to " + \
1449	"validate and if necessary operation and value (sim. values) available " + \
1450	"operations: " + opersinf + " (WRFdiagnosted variables also available: " + \
1451	varNOcheckinf + ")"
1452	statsn = ['minimum', 'maximum', 'mean', 'mean2', 'standard deviation']
1453	gstatsn = ['bias', 'simobs_mean', 'sim_obsmin', 'sim_obsmax', 'sim_obsmean', 'mae', \
1454	'rmse', 'r_pearsoncorr', 'p_pearsoncorr', 'deviation_of_residuals_SDR', \
1455	'indef_of_efficiency_IOE', 'index_of_agreement_IOA', 'fractional_mean_bias_FMB']
1456	ostatsn = ['number of points', 'minimum', 'maximum', 'mean', 'mean2', \
1457	'standard deviation']
1458
1459	parser = OptionParser()
1460	parser.add_option("-d", "--dimensions", dest="dims", help=dimshelp, metavar="VALUES")
1461	parser.add_option("-D", "--vardimensions", dest="vardims",
1462	help=vardimshelp, metavar="VALUES")
1463	parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs,
1464	help=strkObs, metavar="FILE")
1465	parser.add_option("-l", "--stationLocation", dest="stloc",
1466	help="name (\| for spaces), longitude, latitude and height of the station (only for 'single-station')",
1467	metavar="FILE")
1468	parser.add_option("-o", "--observation", dest="fobs",
1469	help="observations file to validate", metavar="FILE")
1470	parser.add_option("-s", "--simulation", dest="fsim",
1471	help="simulation file to validate", metavar="FILE")
1472	parser.add_option("-t", "--trajectoryfile", dest="trajf",
1473	help="file with grid points of the trajectory in the simulation grid ('simtrj')",
1474	metavar="FILE")
1475	parser.add_option("-v", "--variables", dest="vars",
1476	help=varshelp, metavar="VALUES")
1477
1478	(opts, args) = parser.parse_args()
1479
1480	####### ###### ##### #### ### ## #
1481	# Number of different statistics according to the temporal coincidence
1482	# 0: Exact time
1483	# 1: Simulation values closest to observed times
1484	# 2: Simulation values between consecutive observed times
1485	Nstsim = 3
1486
1487	stdescsim = ['E', 'C', 'B']
1488	prestdescsim = ['exact', 'closest', 'between']
1489	Lstdescsim = ['exact time', 'cloest time', 'between observational time-steps']
1490
1491	####### #######
1492	## MAIN
1493	#######
1494
1495	ofile='validation_sim.nc'
1496
1497	if opts.dims is None:
1498	print errormsg
1499	print ' ' + main + ': No list of dimensions are provided!!'
1500	print ' a ',' list of values X@[dimxsim]@[dimxobs],...,T@[dimtsim]@[dimtobs]'+\
1501	' is needed'
1502	quit(-1)
1503	else:
1504	simdims = {}
1505	obsdims = {}
1506	print main +': couple of dimensions _______'
1507	dims = {}
1508	ds = opts.dims.split(',')
1509	for d in ds:
1510	dsecs = d.split('@')
1511	if len(dsecs) != 3:
1512	print errormsg
1513	print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
1514	print ' [DIM]@[dimnsim]@[dimnobs]'
1515	quit(-1)
1516	if dsecs[1] != 'None':
1517	dims[dsecs[0]] = [dsecs[1], dsecs[2]]
1518	simdims[dsecs[0]] = dsecs[1]
1519	obsdims[dsecs[0]] = dsecs[2]
1520
1521	print ' ',dsecs[0],':',dsecs[1],',',dsecs[2]
1522
1523	if opts.vardims is None:
1524	print errormsg
1525	print ' ' + main + ': No list of variables with dimension values are provided!!'
1526	print ' a ',' list of values X@[vardimxsim]@[vardimxobs],...,T@' + \
1527	'[vardimtsim]@[vardimtobs] is needed'
1528	quit(-1)
1529	else:
1530	print main +': couple of variable dimensions _______'
1531	vardims = {}
1532	ds = opts.vardims.split(',')
1533	for d in ds:
1534	dsecs = d.split('@')
1535	if len(dsecs) != 3:
1536	print errormsg
1537	print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
1538	print ' [DIM]@[vardimnsim]@[vardimnobs]'
1539	quit(-1)
1540	if dsecs[1] != 'None':
1541	vardims[dsecs[0]] = [dsecs[1], dsecs[2]]
1542	print ' ',dsecs[0],':',dsecs[1],',',dsecs[2]
1543
1544	if opts.obskind is None:
1545	print errormsg
1546	print ' ' + main + ': No kind of observations provided !!'
1547	quit(-1)
1548	else:
1549	obskind = opts.obskind
1550	if obskind == 'single-station':
1551	if opts.stloc is None:
1552	print errormsg
1553	print ' ' + main + ': No station location provided !!'
1554	quit(-1)
1555	else:
1556	stationdesc = [opts.stloc.split(',')[0].replace('\|',' '), \
1557	np.float(opts.stloc.split(',')[1]), np.float(opts.stloc.split(',')[2]),\
1558	np.float(opts.stloc.split(',')[3])]
1559
1560	if opts.fobs is None:
1561	print errormsg
1562	print ' ' + main + ': No observations file is provided!!'
1563	quit(-1)
1564	else:
1565	if not os.path.isfile(opts.fobs):
1566	print errormsg
1567	print ' ' + main + ": observations file '" + opts.fobs + "' does not exist !!"
1568	quit(-1)
1569
1570	if opts.fsim is None:
1571	print errormsg
1572	print ' ' + main + ': No simulation file is provided!!'
1573	quit(-1)
1574	else:
1575	if not os.path.isfile(opts.fsim):
1576	print errormsg
1577	print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
1578	quit(-1)
1579
1580	if opts.vars is None:
1581	print errormsg
1582	print ' ' + main + ': No list of couples of variables is provided!!'
1583	print ' a ',' list of values [varsim]@[varobs],... is needed'
1584	quit(-1)
1585	else:
1586	valvars = []
1587	vs = opts.vars.split(',')
1588	for v in vs:
1589	vsecs = v.split('@')
1590	if len(vsecs) < 2:
1591	print errormsg
1592	print ' ' + main + ': wrong number of values in:',vsecs, \
1593	' at least 2 are needed !!'
1594	print ' [varsim]@[varobs]@[[oper][val]]'
1595	quit(-1)
1596	if len(vsecs) > 2:
1597	if not searchInlist(simopers,vsecs[2]):
1598	print errormsg
1599	print main + ": operation on simulation values '" + vsecs[2] + \
1600	"' not ready !!"
1601	quit(-1)
1602
1603	valvars.append(vsecs)
1604
1605	# Openning observations trajectory
1606	##
1607	oobs = NetCDFFile(opts.fobs, 'r')
1608
1609	valdimobs = {}
1610	for dn in dims:
1611	print dn,':',dims[dn]
1612	if dims[dn][1] != 'None':
1613	if not oobs.dimensions.has_key(dims[dn][1]):
1614	print errormsg
1615	print ' ' + main + ": observations file does not have dimension '" + \
1616	dims[dn][1] + "' !!"
1617	quit(-1)
1618	if vardims[dn][1] != 'None':
1619	if not oobs.variables.has_key(vardims[dn][1]):
1620	print errormsg
1621	print ' ' + main + ": observations file does not have varibale " + \
1622	"dimension '" + vardims[dn][1] + "' !!"
1623	quit(-1)
1624	valdimobs[dn] = oobs.variables[vardims[dn][1]][:]
1625	else:
1626	if dn == 'X':
1627	valdimobs[dn] = stationdesc[1]
1628	elif dn == 'Y':
1629	valdimobs[dn] = stationdesc[2]
1630	elif dn == 'Z':
1631	valdimobs[dn] = stationdesc[3]
1632
1633	osim = NetCDFFile(opts.fsim, 'r')
1634
1635	valdimsim = {}
1636	for dn in dims:
1637	if dims[dn][0] != 'None':
1638	if not osim.dimensions.has_key(dims[dn][0]):
1639	print errormsg
1640	print ' ' + main + ": simulation file '" + opts.fsim + \
1641	"' does not have dimension '" + dims[dn][0] + "' !!"
1642	print ' it has: ',osim.dimensions
1643	quit(-1)
1644
1645	if not osim.variables.has_key(vardims[dn][0]) and \
1646	not searchInlist(varNOcheck,vardims[dn][0]):
1647	print errormsg
1648	print ' ' + main + ": simulation file '" + opts.fsim + \
1649	"' does not have varibale dimension '" + vardims[dn][0] + "' !!"
1650	print ' it has variables:',osim.variables
1651	quit(-1)
1652	if searchInlist(varNOcheck,vardims[dn][0]):
1653	valdimsim[dn] = compute_varNOcheck(osim, vardims[dn][0])
1654	else:
1655	valdimsim[dn] = osim.variables[vardims[dn][0]][:]
1656
1657	# General characteristics
1658	dimtobs = valdimobs['T'].shape[0]
1659	dimtsim = valdimsim['T'].shape[0]
1660
1661	print main +': observational time-steps:',dimtobs,'simulation:',dimtsim
1662
1663	notfound = np.zeros((dimtobs), dtype=int)
1664
1665	if obskind == 'multi-points':
1666	trajpos = np.zeros((2,dimt),dtype=int)
1667	for it in range(dimtobs):
1668	trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'], \
1669	[valdimobs['X'][it],valdimobss['Y'][it]])
1670	elif obskind == 'single-station':
1671	trajpos = None
1672	stationpos = np.zeros((2), dtype=int)
1673	if valdimsim.has_key('X') and valdimsim.has_key('Y'):
1674	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],[valdimobs['Y'], \
1675	valdimobs['X']])
1676	iid = 0
1677	for idn in osim.variables[vardims['X'][0]].dimensions:
1678	if idn == dims['X'][0]:
1679	stationpos[1] = stsimpos[iid]
1680	elif idn == dims['Y'][0]:
1681	stationpos[0] = stsimpos[iid]
1682
1683	iid = iid + 1
1684	print main + ': station point in simulation:', stationpos
1685	print ' station position:',valdimobs['X'],',',valdimobs['Y']
1686	print ' simulation coord.:',valdimsim['X'][tuple(stsimpos)],',', \
1687	valdimsim['Y'][tuple(stsimpos)]
1688	else:
1689	print main + ': validation with two time-series !!'
1690
1691	elif obskind == 'trajectory':
1692	if opts.trajf is not None:
1693	if not os.path.isfile(opts.fsim):
1694	print errormsg
1695	print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
1696	quit(-1)
1697	else:
1698	otrjf = NetCDFFile(opts.fsim, 'r')
1699	trajpos[0,:] = otrjf.variables['obssimtrj'][0]
1700	trajpos[1,:] = otrjf.variables['obssimtrj'][1]
1701	otrjf.close()
1702	else:
1703	if dims.has_key('Z'):
1704	trajpos = np.zeros((3,dimtobs),dtype=int)
1705	for it in range(dimtobs):
1706	if np.mod(it*100./dimtobs,10.) == 0.:
1707	print ' trajectory done: ',it*100./dimtobs,'%'
1708	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \
1709	[valdimobs['Y'][it],valdimobs['X'][it]])
1710	stationpos = np.zeros((2), dtype=int)
1711	iid = 0
1712	for idn in osim.variables[vardims['X'][0]].dimensions:
1713	if idn == dims['X'][0]:
1714	stationpos[1] = stsimpos[iid]
1715	elif idn == dims['Y'][0]:
1716	stationpos[0] = stsimpos[iid]
1717	iid = iid + 1
1718	if stationpos[0] == 0 and stationpos[1] == 0: notfound[it] = 1
1719
1720	trajpos[0,it] = stationpos[0]
1721	trajpos[1,it] = stationpos[1]
1722	# In the simulation 'Z' varies with time ... non-hydrostatic model! ;)
1723	# trajpos[2,it] = index_mat(valdimsim['Z'][it,:,stationpos[0], \
1724	# stationpos[1]], valdimobs['Z'][it])
1725	else:
1726	trajpos = np.zeros((2,dimtobs),dtype=int)
1727	for it in range(dimtobs):
1728	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \
1729	[valdimobs['Y'][it],valdimobss['X'][it]])
1730	stationpos = np.zeros((2), dtype=int)
1731	iid = 0
1732	for idn in osim.variables[vardims['X'][0]].dimensions:
1733	if idn == dims['X'][0]:
1734	stationpos[1] = stsimpos[iid]
1735	elif idn == dims['Y'][0]:
1736	stationpos[0] = stsimpos[iid]
1737	iid = iid + 1
1738	if stationpos[0] == 0 or stationpos[1] == 0: notfound[it] = 1
1739
1740	trajpos[0,it] = stationspos[0]
1741	trajpos[1,it] = stationspos[1]
1742
1743	print main + ': not found',np.sum(notfound),'points of the trajectory'
1744
1745	# Getting times
1746	tobj = oobs.variables[vardims['T'][1]]
1747	obstunits = tobj.getncattr('units')
1748	if vardims['T'][0] == 'WRFT':
1749	tsim = valdimsim['T'][:]
1750	simtunits = 'seconds since 1949-12-01 00:00:00'
1751	else:
1752	tsim = osim.variables[vardims['T'][0]][:]
1753	otsim = osim.variables[vardims['T'][0]]
1754	simtunits = otsim.getncattr('units')
1755
1756	simobstimes = coincident_CFtimes(tsim, obstunits, simtunits)
1757
1758	#
1759	## Looking for exact/near times
1760	#
1761
1762	# Exact Coincident times
1763	##
1764	exacttvalues0 = []
1765	for it in range(dimtsim):
1766	ot = 0
1767	for ito in range(ot,dimtobs-1):
1768	if valdimobs['T'][ito] == simobstimes[it]:
1769	ot = ito
1770	exacttvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito]])
1771
1772	exacttvalues = np.array(exacttvalues0, dtype=np.float)
1773	Nexactt = len(exacttvalues[:,0])
1774
1775	print main + ': found',Nexactt,'Temporal exact values in simulation and observations'
1776
1777	# Sim Closest times
1778	##
1779	Nsimt = 0
1780	closesttvalues0 = []
1781	closesttvalues0st = []
1782	tsiminit = 0
1783	tsimend = 0
1784
1785	dtsim = simobstimes[1] - simobstimes[0]
1786
1787	for it in range(dimtsim):
1788	ot = 0
1789	for ito in range(ot,dimtobs-1):
1790	if np.abs(valdimobs['T'][ito] - simobstimes[it]) <= dtsim/2.:
1791	ot = ito
1792	tdist = simobstimes[it] - valdimobs['T'][ito]
1793	closesttvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito], \
1794	tdist])
1795	Nsimt = Nsimt + 1
1796	if tsiminit == 0: tsiminit=simobstimes[it]
1797	tsimend = simobstimes[it]
1798
1799	closesttvalues = np.array(closesttvalues0, dtype=np.float)
1800
1801	Nclosest = len(closesttvalues[:,0])
1802	print main + ': found',Nclosest,'Simulation time-values closest to observations'
1803
1804	if Nclosest == 0:
1805	print warnmsg
1806	print main + ': no cclosest times found !!'
1807	print ' stopping it'
1808	quit(-1)
1809
1810	# Sim Coincident times
1811	##
1812	Nsimt = 0
1813	coindtvalues0 = []
1814	coindtvalues0st = []
1815	tsiminit = 0
1816	tsimend = 0
1817
1818	for it in range(dimtsim):
1819	ot = 0
1820	for ito in range(ot,dimtobs-1):
1821	if valdimobs['T'][ito] < simobstimes[it] and valdimobs['T'][ito+1] > \
1822	simobstimes[it]:
1823	ot = ito
1824	tdist = simobstimes[it] - valdimobs['T'][ito]
1825	coindtvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito], \
1826	tdist])
1827	Nsimt = Nsimt + 1
1828	if tsiminit == 0: tsiminit=simobstimes[it]
1829	tsimend = simobstimes[it]
1830	elif simobstimes[it] > valdimobs['T'][ito+1]:
1831	coindtvalues0st.append([Nsimt, ito, valdimobs['T'][ito],tsimend-tsiminit])
1832
1833	coindtvalues = np.array(coindtvalues0, dtype=np.float)
1834	coindtvaluesst = np.array(coindtvalues0st, dtype=np.float)
1835
1836	Ncoindt = len(coindtvalues[:,0])
1837	print main + ': found',Ncoindt,'Simulation time-interval (within consecutive ' + \
1838	'observed times) coincident times between simulation and observations'
1839
1840	if Ncoindt == 0:
1841	print warnmsg
1842	print main + ': no coincident times found !!'
1843	print ' stopping it'
1844	quit(-1)
1845
1846	# Validating
1847	##
1848
1849	onewnc = NetCDFFile(ofile, 'w')
1850
1851	# Dimensions
1852	for kst in range(Nstsim):
1853	newdim = onewnc.createDimension(prestdescsim[kst] + 'time',None)
1854	if stdescsim[kst] != 'E':
1855	newdim = onewnc.createDimension(prestdescsim[kst] + 'obstime',None)
1856
1857	newdim = onewnc.createDimension('bnds',2)
1858	newdim = onewnc.createDimension('couple',2)
1859	newdim = onewnc.createDimension('StrLength',StringLength)
1860	newdim = onewnc.createDimension('xaround',Ngrid*2+1)
1861	newdim = onewnc.createDimension('yaround',Ngrid*2+1)
1862	newdim = onewnc.createDimension('gstats',13)
1863	newdim = onewnc.createDimension('stats',5)
1864	newdim = onewnc.createDimension('tstats',6)
1865	newdim = onewnc.createDimension('Nstsim', 3)
1866
1867	# Variable dimensions
1868	##
1869	newvar = onewnc.createVariable('couple', 'c', ('couple','StrLength'))
1870	basicvardef(newvar, 'couple', 'couples of values', '-')
1871	writing_str_nc(newvar, ['sim','obs'], StringLength)
1872
1873	newvar = onewnc.createVariable('statistics', 'c', ('stats','StrLength'))
1874	basicvardef(newvar, 'statistics', 'statitics from values', '-')
1875	writing_str_nc(newvar, statsn, StringLength)
1876
1877	newvar = onewnc.createVariable('gstatistics', 'c', ('gstats','StrLength'))
1878	basicvardef(newvar, 'gstatistics', 'global statitics from values', '-')
1879	writing_str_nc(newvar, gstatsn, StringLength)
1880
1881	newvar = onewnc.createVariable('tstatistics', 'c', ('tstats','StrLength'))
1882	basicvardef(newvar, 'tstatistics', 'statitics from values along time', '-')
1883	writing_str_nc(newvar, ostatsn, StringLength)
1884
1885	newvar = onewnc.createVariable('ksimstatistics', 'c', ('Nstsim','StrLength'))
1886	basicvardef(newvar, 'ksimstatistics', 'kind of simulated statitics', '-')
1887	writing_str_nc(newvar, Lstdescsim, StringLength)
1888
1889	if obskind == 'trajectory':
1890	if dims.has_key('Z'):
1891	newdim = onewnc.createDimension('trj',3)
1892	else:
1893	newdim = onewnc.createDimension('trj',2)
1894
1895	newvar = onewnc.createVariable('obssimtrj','i',('obstime','trj'))
1896	basicvardef(newvar, 'obssimtrj', 'trajectory on the simulation grid', '-')
1897	newvar[:] = trajpos.transpose()
1898
1899	if dims.has_key('Z'):
1900	newdim = onewnc.createDimension('simtrj',4)
1901	else:
1902	newdim = onewnc.createDimension('simtrj',3)
1903
1904	Nvars = len(valvars)
1905	for ivar in range(Nvars):
1906	simobsvalues = []
1907
1908	varsimobs = valvars[ivar][0] + '_' + valvars[ivar][1]
1909	print ' ' + varsimobs + '... .. .'
1910
1911	if not oobs.variables.has_key(valvars[ivar][1]):
1912	print errormsg
1913	print ' ' + main + ": observations file has not '" + valvars[ivar][1] + \
1914	"' !!"
1915	quit(-1)
1916
1917	if not osim.variables.has_key(valvars[ivar][0]):
1918	if not searchInlist(varNOcheck, valvars[ivar][0]):
1919	print errormsg
1920	print ' ' + main + ": simulation file has not '" + valvars[ivar][0] + \
1921	"' !!"
1922	quit(-1)
1923	else:
1924	ovsim = compute_varNOcheck(osim, valvars[ivar][0])
1925	else:
1926	ovsim = osim.variables[valvars[ivar][0]]
1927
1928	for idn in ovsim.dimensions:
1929	if not searchInlist(simdims.values(),idn):
1930	print errormsg
1931	print ' ' + main + ": dimension '" + idn + "' of variable '" + \
1932	valvars[ivar][0] + "' not provided as reference coordinate [X,Y,Z,T] !!"
1933	quit(-1)
1934
1935	ovobs = oobs.variables[valvars[ivar][1]]
1936	if searchInlist(ovobs.ncattrs(),'_FillValue'):
1937	oFillValue = ovobs.getncattr('_FillValue')
1938	else:
1939	oFillValue = None
1940
1941	for kst in range(Nstsim):
1942	simobsvalues = []
1943
1944	timedn = prestdescsim[kst] + 'time'
1945	timeobsdn = prestdescsim[kst] + 'obstime'
1946	print ' ' + prestdescsim[kst] + ' ...'
1947
1948	if stdescsim[kst] == 'E':
1949	# Observed and simualted exact times
1950	simobsvalues, simobsSvalues, simobsTtvalues, trjsim = \
1951	getting_ValidationValues(obskind, Nexactt, dims, trajpos, ovsim, \
1952	ovobs, exacttvalues, oFillValue, Ngrid, 'instantaneous')
1953
1954	if ivar == 0:
1955	vname = prestdescsim[kst] + 'time'
1956	newvar = onewnc.createVariable(vname,'f8', (timedn))
1957	basicvardef(newvar, vname, 'exact coincident time observations and '+\
1958	'simulation', obstunits)
1959	set_attribute(newvar, 'calendar', 'standard')
1960	newvar[:] = exacttvalues[:,3]
1961
1962	dimt = Nexactt
1963
1964	elif stdescsim[kst] == 'C':
1965	# Simualted closest to Observed times
1966	simobsvalues, simobsSvalues, simobsTtvalues, trjsim = \
1967	getting_ValidationValues(obskind, Nclosest, dims, trajpos, ovsim, \
1968	ovobs, closesttvalues, oFillValue, Ngrid, 'instantaneous')
1969	dimt = Nclosest
1970
1971
1972	if ivar == 0:
1973	vname = prestdescsim[kst] + 'time'
1974	newvar = onewnc.createVariable(vname,'f8', (timedn))
1975	basicvardef(newvar, vname, 'time simulations closest to observed ' + \
1976	'values', obstunits )
1977	set_attribute(newvar, 'calendar', 'standard')
1978	newvar[:] = closesttvalues[:,2]
1979
1980	vname = prestdescsim[kst] + 'obstime'
1981	newvar = onewnc.createVariable(vname,'f8', (vname))
1982	basicvardef(newvar, vname, 'closest time observations', obstunits)
1983	set_attribute(newvar, 'calendar', 'standard')
1984	newvar[:] = closesttvalues[:,3]
1985
1986	elif stdescsim[kst] == 'B':
1987	# Observed values temporally around coincident times
1988	simobsvalues, simobsSvalues, simobsTtvalues, trjsim = \
1989	getting_ValidationValues(obskind, Ncoindt, dims, trajpos, ovsim, ovobs,\
1990	coindtvalues, oFillValue, Ngrid, 'tbackwardSmean')
1991	dimt = simobsSvalues.shape[0]
1992
1993	if ivar == 0:
1994	vname = prestdescsim[kst] + 'time'
1995	newvar = onewnc.createVariable(vname,'f8', (timedn), \
1996	fill_value = fillValueF)
1997	basicvardef(newvar, vname, 'simulation time between observations', \
1998	obstunits)
1999	set_attribute(newvar, 'calendar', 'standard')
2000	set_attribute(newvar, 'bounds', 'time_bnds')
2001	newvar[:] = simobsTtvalues[:,1]
2002
2003	vname = prestdescsim[kst] + 'obstime'
2004	newvar = onewnc.createVariable(vname,'f8', (vname))
2005	basicvardef(newvar, vname, 'observed between time', obstunits )
2006	set_attribute(newvar, 'calendar', 'standard')
2007	newvar[:] = np.unique(coindtvalues[:,3])
2008
2009	# Re-arranging values...
2010	arrayvals = np.array(simobsvalues)
2011	if len(valvars[ivar]) > 2:
2012	const=np.float(valvars[ivar][3])
2013	if valvars[ivar][2] == 'sumc':
2014	simobsSvalues = simobsSvalues + const
2015	arrayvals[:,0] = arrayvals[:,0] + const
2016	elif valvars[ivar][2] == 'subc':
2017	simobsSvalues = simobsSvalues - const
2018	arrayvals[:,0] = arrayvals[:,0] - const
2019	elif valvars[ivar][2] == 'mulc':
2020	simobsSvalues = simobsSvalues * const
2021	arrayvals[:,0] = arrayvals[:,0] * const
2022	elif valvars[ivar][2] == 'divc':
2023	simobsSvalues = simobsSvalues / const
2024	arrayvals[:,0] = arrayvals[:,0] / const
2025	else:
2026	print errormsg
2027	print ' ' + fname + ": operation '"+valvars[ivar][2]+"' not ready!!"
2028	quit(-1)
2029
2030	if kst == 0:
2031	simstats = np.zeros((Nstsim,5), dtype=np.float)
2032	obsstats = np.zeros((Nstsim,5), dtype=np.float)
2033	simobsstats = np.zeros((Nstsim,13), dtype=np.float)
2034
2035	# statisics sim
2036	simstats[kst,0] = np.min(arrayvals[:,0])
2037	simstats[kst,1] = np.max(arrayvals[:,0])
2038	simstats[kst,2] = np.mean(arrayvals[:,0])
2039	simstats[kst,3] = np.mean(arrayvals[:,0]*arrayvals[:,0])
2040	simstats[kst,4] = np.sqrt(simstats[kst,3] - simstats[kst,2]*simstats[kst,2])
2041
2042	# statisics obs
2043	# Masking 'nan'
2044	obsmask0 = np.where(arrayvals[:,1] != arrayvals[:,1], fillValueF, \
2045	arrayvals[:,1])
2046
2047	obsmask = ma.masked_equal(obsmask0, fillValueF)
2048	obsmask2 = obsmask*obsmask
2049
2050	obsstats[kst,0] = obsmask.min()
2051	obsstats[kst,1] = obsmask.max()
2052	obsstats[kst,2] = obsmask.mean()
2053	obsstats[kst,3] = obsmask2.mean()
2054	obsstats[kst,4] = np.sqrt(obsstats[kst,3] - obsstats[kst,2]*obsstats[kst,2])
2055
2056	# Statistics sim-obs
2057	diffvals = np.zeros((dimt), dtype=np.float)
2058
2059	diffvals = arrayvals[:,0] - obsmask
2060
2061	diff2vals = diffvals*diffvals
2062	sumdiff = diffvals.sum()
2063	sumdiff2 = diff2vals.sum()
2064
2065	simobsstats[kst,0] = simstats[kst,0] - obsstats[kst,0]
2066	simobsstats[kst,1] = np.mean(arrayvals[:,0]*obsmask)
2067	simobsstats[kst,2] = diffvals.min()
2068	simobsstats[kst,3] = diffvals.max()
2069	simobsstats[kst,4] = diffvals.mean()
2070	simobsstats[kst,5] = np.abs(diffvals).mean()
2071	simobsstats[kst,6] = np.sqrt(diff2vals.mean())
2072	simobsstats[kst,7], simobsstats[kst,8] = sts.mstats.pearsonr(arrayvals[:,0], \
2073	obsmask)
2074	# From:
2075	#Willmott, C. J. 1981. 'On the validation of models. Physical Geography', 2, 184-194
2076	#Willmott, C. J. (1984). 'On the evaluation of model performance in physical
2077	# geography'. Spatial Statistics and Models, G. L. Gaile and C. J. Willmott, eds.,
2078	# 443-460
2079	#Willmott, C. J., S. G. Ackleson, R. E. Davis, J. J. Feddema, K. M. Klink, D. R.
2080	# Legates, J. O'Donnell, and C. M. Rowe (1985), 'Statistics for the Evaluation and
2081	# Comparison of Models', J. Geophys. Res., 90(C5), 8995-9005
2082	#Legates, D. R., and G. J. McCabe Jr. (1999), 'Evaluating the Use of "Goodness-of-Fit"
2083	# Measures in Hydrologic and Hydroclimatic Model Validation', Water Resour. Res.,
2084	# 35(1), 233-241
2085	#
2086	# Deviation of residuals (SDR)
2087	simobsstats[kst,9] = np.mean(np.sqrt(np.abs((diffvals-simobsstats[kst,0])* \
2088	(diffvals-obsstats[kst,0]))))
2089	# Index of Efficiency (IOE)
2090	obsbias2series = (obsmask - obsstats[kst,0])*(obsmask - obsstats[kst,0])
2091	sumobsbias2series = obsbias2series.sum()
2092
2093	simobsstats[kst,10] = 1. - sumdiff2/(sumobsbias2series)
2094	# Index of Agreement (IOA)
2095	simbias2series = arrayvals[:,0] - obsstats[kst,0]
2096	obsbias2series = obsmask - obsstats[kst,0]
2097
2098	abssimbias2series = np.abs(simbias2series)
2099	absobsbias2series = np.where(obsbias2series<0, -obsbias2series, \
2100	obsbias2series)
2101
2102	abssimobsbias2series = (abssimbias2series+absobsbias2series)*( \
2103	abssimbias2series + absobsbias2series)
2104
2105	simobsstats[kst,11] = 1. - sumdiff2/(abssimobsbias2series.sum())
2106	# Fractional Mean Bias (FMB)
2107	simobsstats[kst,12]=(simstats[kst,0]-obsstats[kst,0])/(0.5*(simstats[kst,0] +\
2108	obsstats[kst,0]))
2109
2110	# Statistics around sim values
2111	aroundstats = np.zeros((5,dimt), dtype=np.float)
2112	for it in range(dimt):
2113	aroundstats[0,it] = np.min(simobsSvalues[it,])
2114	aroundstats[1,it] = np.max(simobsSvalues[it,])
2115	aroundstats[2,it] = np.mean(simobsSvalues[it,])
2116	aroundstats[3,it] = np.mean(simobsSvalues[it,]*simobsSvalues[it,])
2117	aroundstats[4,it] = np.sqrt(aroundstats[3,it] - aroundstats[2,it]* \
2118	aroundstats[2,it])
2119
2120	# sim Values to netCDF
2121	newvar = onewnc.createVariable(valvars[ivar][0] + '_' + stdescsim[kst] + \
2122	'sim', 'f', (timedn), fill_value=fillValueF)
2123	descvar = prestdescsim[kst] + ' time simulated: ' + valvars[ivar][0]
2124	basicvardef(newvar, valvars[ivar][0], descvar, ovobs.getncattr('units'))
2125	newvar[:] = arrayvals[:,0]
2126
2127	# obs Values to netCDF
2128	if stdescsim[kst] != 'E':
2129	newvar = onewnc.createVariable(valvars[ivar][1] + '_' + stdescsim[kst] + \
2130	'obs', 'f', (timeobsdn), fill_value=fillValueF)
2131	else:
2132	newvar = onewnc.createVariable(valvars[ivar][1] + '_' + stdescsim[kst] + \
2133	'obs', 'f', (timedn), fill_value=fillValueF)
2134
2135	descvar = prestdescsim[kst] + ' time observed: ' + valvars[ivar][1]
2136	basicvardef(newvar, valvars[ivar][1], descvar, ovobs.getncattr('units'))
2137	newvar[:] = arrayvals[:,1]
2138
2139	# Around values
2140	if not onewnc.variables.has_key(valvars[ivar][0] + 'around'):
2141	vname = prestdescsim[kst] + valvars[ivar][0] + 'around'
2142	else:
2143	vname = prestdescsim[kst] + valvars[ivar][0] + 'Around'
2144
2145	if dims.has_key('Z'):
2146	if not onewnc.dimensions.has_key('zaround'):
2147	newdim = onewnc.createDimension('zaround',Ngrid*2+1)
2148	newvar = onewnc.createVariable(vname, 'f', (timedn,'zaround', \
2149	'yaround','xaround'), fill_value=fillValueF)
2150	else:
2151	newvar = onewnc.createVariable(vname, 'f', (timedn,'yaround','xaround'), \
2152	fill_value=fillValueF)
2153
2154	descvar = prestdescsim[kst] + 'around simulated values +/- grid values: ' + \
2155	valvars[ivar][0]
2156	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2157	newvar[:] = simobsSvalues
2158
2159
2160	# around sim Statistics
2161	if not searchInlist(onewnc.variables,prestdescsim[kst] + valvars[ivar][0] + \
2162	'staround'):
2163	vname = prestdescsim[kst] + valvars[ivar][0] + 'staround'
2164	else:
2165	vname = prestdescsim[kst] + valvars[ivar][0] + 'Staround'
2166
2167	newvar = onewnc.createVariable(vname, 'f', (timedn,'stats'), \
2168	fill_value=fillValueF)
2169	descvar = prestdescsim[kst] + ' around (' + str(Ngrid) + ', ' + str(Ngrid) +\
2170	') simulated statisitcs: ' + valvars[ivar][0]
2171	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2172	set_attribute(newvar, 'cell_methods', 'time_bnds')
2173	newvar[:] = aroundstats.transpose()
2174
2175	if stdescsim[kst] == 'B':
2176	if not searchInlist(onewnc.variables, 'time_bnds'):
2177	newvar = onewnc.createVariable('time_bnds','f8',(timedn,'bnds'))
2178	basicvardef(newvar, 'time_bnds', timedn, obstunits )
2179	set_attribute(newvar, 'calendar', 'standard')
2180	newvar[:] = simobsTtvalues
2181
2182	# sim Statistics
2183	if not searchInlist(onewnc.variables,valvars[ivar][0] + 'stsim'):
2184	vname = valvars[ivar][0] + 'stsim'
2185	else:
2186	vname = valvars[ivar][0] + 'stSim'
2187
2188	newvar = onewnc.createVariable(vname, 'f', ('Nstsim', 'stats'), \
2189	fill_value=fillValueF)
2190	descvar = 'simulated statisitcs: ' + valvars[ivar][0]
2191	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2192	newvar[:] = simstats
2193
2194	# obs Statistics
2195	if not searchInlist(onewnc.variables,valvars[ivar][1] + 'stobs'):
2196	vname = valvars[ivar][1] + 'stobs'
2197	else:
2198	vname = valvars[ivar][1] + 'stObs'
2199
2200	newvar = onewnc.createVariable(vname, 'f', ('Nstsim', 'stats'), \
2201	fill_value=fillValueF)
2202	descvar = 'observed statisitcs: ' + valvars[ivar][1]
2203	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2204	newvar[:] = obsstats
2205
2206	# sim-obs Statistics
2207	if not searchInlist(onewnc.variables,varsimobs + 'st'):
2208	vname = varsimobs + 'st'
2209	else:
2210	vname = varSimobs + 'st'
2211
2212	newvar = onewnc.createVariable(vname, 'f', ('Nstsim', 'gstats'), \
2213	fill_value=fillValueF)
2214	descvar = 'simulated-observed tatisitcs: ' + varsimobs
2215	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2216	newvar[:] = simobsstats
2217
2218	onewnc.sync()
2219
2220	if trjsim is not None:
2221	newvar = onewnc.createVariable('simtrj','i',('betweentime','simtrj'))
2222	basicvardef(newvar,'simtrj','coordinates [X,Y,Z,T] of the coincident ' + \
2223	'trajectory in sim', obstunits)
2224	newvar[:] = trjsim.transpose()
2225
2226	# Adding three variables with the station location, longitude, latitude and height
2227	if obskind == 'single-station':
2228	adding_station_desc(onewnc,stationdesc)
2229
2230	# Global attributes
2231	##
2232	set_attribute(onewnc,'author_nc','Lluis Fita')
2233	set_attribute(onewnc,'institution_nc','Laboratoire de Meteorology Dynamique, ' + \
2234	'LMD-Jussieu, UPMC, Paris')
2235	set_attribute(onewnc,'country_nc','France')
2236	set_attribute(onewnc,'script_nc',main)
2237	set_attribute(onewnc,'version_script',version)
2238	set_attribute(onewnc,'information', \
2239	'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html')
2240	set_attribute(onewnc,'simfile',opts.fsim)
2241	set_attribute(onewnc,'obsfile',opts.fobs)
2242
2243	onewnc.sync()
2244	onewnc.close()
2245
2246	print main + ": successfull writting of '" + ofile + "' !!"

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

Download in other formats: