Context Navigation

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

Last change on this file since 521 was 521, checked in by lfita, 9 years ago
Removing debugging prints and related quit
File size: 89.0 KB

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 between consecutive observed times
1483	Nstsim = 2
1484
1485	stdescsim = ['E', 'B']
1486	Lstdescsim = ['exact time', 'between observational time-steps']
1487
1488	####### #######
1489	## MAIN
1490	#######
1491
1492	ofile='validation_sim.nc'
1493
1494	if opts.dims is None:
1495	print errormsg
1496	print ' ' + main + ': No list of dimensions are provided!!'
1497	print ' a ',' list of values X@[dimxsim]@[dimxobs],...,T@[dimtsim]@[dimtobs]'+\
1498	' is needed'
1499	quit(-1)
1500	else:
1501	simdims = {}
1502	obsdims = {}
1503	print main +': couple of dimensions _______'
1504	dims = {}
1505	ds = opts.dims.split(',')
1506	for d in ds:
1507	dsecs = d.split('@')
1508	if len(dsecs) != 3:
1509	print errormsg
1510	print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
1511	print ' [DIM]@[dimnsim]@[dimnobs]'
1512	quit(-1)
1513	if dsecs[1] != 'None':
1514	dims[dsecs[0]] = [dsecs[1], dsecs[2]]
1515	simdims[dsecs[0]] = dsecs[1]
1516	obsdims[dsecs[0]] = dsecs[2]
1517
1518	print ' ',dsecs[0],':',dsecs[1],',',dsecs[2]
1519
1520	if opts.vardims is None:
1521	print errormsg
1522	print ' ' + main + ': No list of variables with dimension values are provided!!'
1523	print ' a ',' list of values X@[vardimxsim]@[vardimxobs],...,T@' + \
1524	'[vardimtsim]@[vardimtobs] is needed'
1525	quit(-1)
1526	else:
1527	print main +': couple of variable dimensions _______'
1528	vardims = {}
1529	ds = opts.vardims.split(',')
1530	for d in ds:
1531	dsecs = d.split('@')
1532	if len(dsecs) != 3:
1533	print errormsg
1534	print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
1535	print ' [DIM]@[vardimnsim]@[vardimnobs]'
1536	quit(-1)
1537	if dsecs[1] != 'None':
1538	vardims[dsecs[0]] = [dsecs[1], dsecs[2]]
1539	print ' ',dsecs[0],':',dsecs[1],',',dsecs[2]
1540
1541	if opts.obskind is None:
1542	print errormsg
1543	print ' ' + main + ': No kind of observations provided !!'
1544	quit(-1)
1545	else:
1546	obskind = opts.obskind
1547	if obskind == 'single-station':
1548	if opts.stloc is None:
1549	print errormsg
1550	print ' ' + main + ': No station location provided !!'
1551	quit(-1)
1552	else:
1553	stationdesc = [opts.stloc.split(',')[0].replace('\|',' '), \
1554	np.float(opts.stloc.split(',')[1]), np.float(opts.stloc.split(',')[2]),\
1555	np.float(opts.stloc.split(',')[3])]
1556
1557	if opts.fobs is None:
1558	print errormsg
1559	print ' ' + main + ': No observations file is provided!!'
1560	quit(-1)
1561	else:
1562	if not os.path.isfile(opts.fobs):
1563	print errormsg
1564	print ' ' + main + ": observations file '" + opts.fobs + "' does not exist !!"
1565	quit(-1)
1566
1567	if opts.fsim is None:
1568	print errormsg
1569	print ' ' + main + ': No simulation file is provided!!'
1570	quit(-1)
1571	else:
1572	if not os.path.isfile(opts.fsim):
1573	print errormsg
1574	print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
1575	quit(-1)
1576
1577	if opts.vars is None:
1578	print errormsg
1579	print ' ' + main + ': No list of couples of variables is provided!!'
1580	print ' a ',' list of values [varsim]@[varobs],... is needed'
1581	quit(-1)
1582	else:
1583	valvars = []
1584	vs = opts.vars.split(',')
1585	for v in vs:
1586	vsecs = v.split('@')
1587	if len(vsecs) < 2:
1588	print errormsg
1589	print ' ' + main + ': wrong number of values in:',vsecs, \
1590	' at least 2 are needed !!'
1591	print ' [varsim]@[varobs]@[[oper][val]]'
1592	quit(-1)
1593	if len(vsecs) > 2:
1594	if not searchInlist(simopers,vsecs[2]):
1595	print errormsg
1596	print main + ": operation on simulation values '" + vsecs[2] + \
1597	"' not ready !!"
1598	quit(-1)
1599
1600	valvars.append(vsecs)
1601
1602	# Openning observations trajectory
1603	##
1604	oobs = NetCDFFile(opts.fobs, 'r')
1605
1606	valdimobs = {}
1607	for dn in dims:
1608	print dn,':',dims[dn]
1609	if dims[dn][1] != 'None':
1610	if not oobs.dimensions.has_key(dims[dn][1]):
1611	print errormsg
1612	print ' ' + main + ": observations file does not have dimension '" + \
1613	dims[dn][1] + "' !!"
1614	quit(-1)
1615	if vardims[dn][1] != 'None':
1616	if not oobs.variables.has_key(vardims[dn][1]):
1617	print errormsg
1618	print ' ' + main + ": observations file does not have varibale " + \
1619	"dimension '" + vardims[dn][1] + "' !!"
1620	quit(-1)
1621	valdimobs[dn] = oobs.variables[vardims[dn][1]][:]
1622	else:
1623	if dn == 'X':
1624	valdimobs[dn] = stationdesc[1]
1625	elif dn == 'Y':
1626	valdimobs[dn] = stationdesc[2]
1627	elif dn == 'Z':
1628	valdimobs[dn] = stationdesc[3]
1629
1630	osim = NetCDFFile(opts.fsim, 'r')
1631
1632	valdimsim = {}
1633	for dn in dims:
1634	if dims[dn][0] != 'None':
1635	if not osim.dimensions.has_key(dims[dn][0]):
1636	print errormsg
1637	print ' ' + main + ": simulation file '" + opts.fsim + \
1638	"' does not have dimension '" + dims[dn][0] + "' !!"
1639	print ' it has: ',osim.dimensions
1640	quit(-1)
1641
1642	if not osim.variables.has_key(vardims[dn][0]) and \
1643	not searchInlist(varNOcheck,vardims[dn][0]):
1644	print errormsg
1645	print ' ' + main + ": simulation file '" + opts.fsim + \
1646	"' does not have varibale dimension '" + vardims[dn][0] + "' !!"
1647	print ' it has variables:',osim.variables
1648	quit(-1)
1649	if searchInlist(varNOcheck,vardims[dn][0]):
1650	valdimsim[dn] = compute_varNOcheck(osim, vardims[dn][0])
1651	else:
1652	valdimsim[dn] = osim.variables[vardims[dn][0]][:]
1653
1654	# General characteristics
1655	dimtobs = valdimobs['T'].shape[0]
1656	dimtsim = valdimsim['T'].shape[0]
1657
1658	print main +': observational time-steps:',dimtobs,'simulation:',dimtsim
1659
1660	notfound = np.zeros((dimtobs), dtype=int)
1661
1662	if obskind == 'multi-points':
1663	trajpos = np.zeros((2,dimt),dtype=int)
1664	for it in range(dimtobs):
1665	trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'], \
1666	[valdimobs['X'][it],valdimobss['Y'][it]])
1667	elif obskind == 'single-station':
1668	trajpos = None
1669	stationpos = np.zeros((2), dtype=int)
1670	if valdimsim.has_key('X') and valdimsim.has_key('Y'):
1671	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],[valdimobs['Y'], \
1672	valdimobs['X']])
1673	iid = 0
1674	for idn in osim.variables[vardims['X'][0]].dimensions:
1675	if idn == dims['X'][0]:
1676	stationpos[1] = stsimpos[iid]
1677	elif idn == dims['Y'][0]:
1678	stationpos[0] = stsimpos[iid]
1679
1680	iid = iid + 1
1681	print main + ': station point in simulation:', stationpos
1682	print ' station position:',valdimobs['X'],',',valdimobs['Y']
1683	print ' simulation coord.:',valdimsim['X'][tuple(stsimpos)],',', \
1684	valdimsim['Y'][tuple(stsimpos)]
1685	else:
1686	print main + ': validation with two time-series !!'
1687
1688	elif obskind == 'trajectory':
1689	if opts.trajf is not None:
1690	if not os.path.isfile(opts.fsim):
1691	print errormsg
1692	print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
1693	quit(-1)
1694	else:
1695	otrjf = NetCDFFile(opts.fsim, 'r')
1696	trajpos[0,:] = otrjf.variables['obssimtrj'][0]
1697	trajpos[1,:] = otrjf.variables['obssimtrj'][1]
1698	otrjf.close()
1699	else:
1700	if dims.has_key('Z'):
1701	trajpos = np.zeros((3,dimtobs),dtype=int)
1702	for it in range(dimtobs):
1703	if np.mod(it*100./dimtobs,10.) == 0.:
1704	print ' trajectory done: ',it*100./dimtobs,'%'
1705	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \
1706	[valdimobs['Y'][it],valdimobs['X'][it]])
1707	stationpos = np.zeros((2), dtype=int)
1708	iid = 0
1709	for idn in osim.variables[vardims['X'][0]].dimensions:
1710	if idn == dims['X'][0]:
1711	stationpos[1] = stsimpos[iid]
1712	elif idn == dims['Y'][0]:
1713	stationpos[0] = stsimpos[iid]
1714	iid = iid + 1
1715	if stationpos[0] == 0 and stationpos[1] == 0: notfound[it] = 1
1716
1717	trajpos[0,it] = stationpos[0]
1718	trajpos[1,it] = stationpos[1]
1719	# In the simulation 'Z' varies with time ... non-hydrostatic model! ;)
1720	# trajpos[2,it] = index_mat(valdimsim['Z'][it,:,stationpos[0], \
1721	# stationpos[1]], valdimobs['Z'][it])
1722	else:
1723	trajpos = np.zeros((2,dimtobs),dtype=int)
1724	for it in range(dimtobs):
1725	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \
1726	[valdimobs['Y'][it],valdimobss['X'][it]])
1727	stationpos = np.zeros((2), dtype=int)
1728	iid = 0
1729	for idn in osim.variables[vardims['X'][0]].dimensions:
1730	if idn == dims['X'][0]:
1731	stationpos[1] = stsimpos[iid]
1732	elif idn == dims['Y'][0]:
1733	stationpos[0] = stsimpos[iid]
1734	iid = iid + 1
1735	if stationpos[0] == 0 or stationpos[1] == 0: notfound[it] = 1
1736
1737	trajpos[0,it] = stationspos[0]
1738	trajpos[1,it] = stationspos[1]
1739
1740	print main + ': not found',np.sum(notfound),'points of the trajectory'
1741
1742	# Getting times
1743	tobj = oobs.variables[vardims['T'][1]]
1744	obstunits = tobj.getncattr('units')
1745	if vardims['T'][0] == 'WRFT':
1746	tsim = valdimsim['T'][:]
1747	simtunits = 'seconds since 1949-12-01 00:00:00'
1748	else:
1749	tsim = osim.variables[vardims['T'][0]][:]
1750	otsim = osim.variables[vardims['T'][0]]
1751	simtunits = otsim.getncattr('units')
1752
1753	simobstimes = coincident_CFtimes(tsim, obstunits, simtunits)
1754
1755	#
1756	## Looking for exact/near times
1757	#
1758
1759	# Exact Coincident times
1760	##
1761	exacttvalues0 = []
1762	for it in range(dimtsim):
1763	ot = 0
1764	for ito in range(ot,dimtobs-1):
1765	if valdimobs['T'][ito] == simobstimes[it]:
1766	ot = ito
1767	exacttvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito]])
1768
1769	exacttvalues = np.array(exacttvalues0, dtype=np.float)
1770	Nexactt = len(exacttvalues[:,0])
1771
1772	print main + ': found',Nexactt,'Temporal same values in simulation and observations'
1773
1774	# Sim Coincident times
1775	##
1776	sumsimval = 0.
1777	sum2simval = 0.
1778	Nsimt = 0
1779	coindtvalues0 = []
1780	coindtvalues0st = []
1781	tsiminit = 0
1782	tsimend = 0
1783
1784	for it in range(dimtsim):
1785	ot = 0
1786	for ito in range(ot,dimtobs-1):
1787	if valdimobs['T'][ito] < simobstimes[it] and valdimobs['T'][ito+1] > \
1788	simobstimes[it]:
1789	ot = ito
1790	tdist = simobstimes[it] - valdimobs['T'][ito]
1791	coindtvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito], \
1792	tdist])
1793	Nsimt = Nsimt + 1
1794	if tsiminit == 0: tsiminit=simobstimes[it]
1795	tsimend = simobstimes[it]
1796	elif simobstimes[it] > valdimobs['T'][ito+1]:
1797	coindtvalues0st.append([Nsimt, ito, valdimobs['T'][ito],tsimend-tsiminit])
1798
1799	coindtvalues = np.array(coindtvalues0, dtype=np.float)
1800	coindtvaluesst = np.array(coindtvalues0st, dtype=np.float)
1801
1802	Ncoindt = len(coindtvalues[:,0])
1803	print main + ': found',Ncoindt,'Simulation time-interval (within consecutive ' + \
1804	'observed times) coincident times between simulation and observations'
1805
1806	if Ncoindt == 0:
1807	print warnmsg
1808	print main + ': no coincident times found !!'
1809	print ' stopping it'
1810	quit(-1)
1811
1812	# Validating
1813	##
1814
1815	onewnc = NetCDFFile(ofile, 'w')
1816
1817	# Dimensions
1818	newdim = onewnc.createDimension('time',None)
1819	newdim = onewnc.createDimension('betweentime',None)
1820	newdim = onewnc.createDimension('bnds',2)
1821	newdim = onewnc.createDimension('obstime',None)
1822	newdim = onewnc.createDimension('couple',2)
1823	newdim = onewnc.createDimension('StrLength',StringLength)
1824	newdim = onewnc.createDimension('xaround',Ngrid*2+1)
1825	newdim = onewnc.createDimension('yaround',Ngrid*2+1)
1826	newdim = onewnc.createDimension('gstats',13)
1827	newdim = onewnc.createDimension('stats',5)
1828	newdim = onewnc.createDimension('tstats',6)
1829	newdim = onewnc.createDimension('Nstsim', 2)
1830
1831	# Variable dimensions
1832	##
1833	newvar = onewnc.createVariable('exacttime','f8',('time'))
1834	basicvardef(newvar, 'obstime', 'exact coincident time observations and simulation', obstunits)
1835	set_attribute(newvar, 'calendar', 'standard')
1836	print 'Lluis Nexacttvalues:',len(exacttvalues[:,3])
1837	newvar[:] = exacttvalues[:,3]
1838
1839	newvar = onewnc.createVariable('obstime','f8',('obstime'))
1840	basicvardef(newvar, 'obstime', 'time observations for between values', obstunits)
1841	set_attribute(newvar, 'calendar', 'standard')
1842	set_attribute(newvar, 'bounds', 'time_bnds')
1843	newvar[:] = coindtvalues[:,3]
1844
1845	newvar = onewnc.createVariable('betweentime','f8',('betweentime'))
1846	basicvardef(newvar, 'obstime', 'time simulations for between values', simtunits )
1847	set_attribute(newvar, 'calendar', 'standard')
1848	set_attribute(newvar, 'bounds', 'time_bnds')
1849	newvar[:] = coindtvalues[:,2]
1850
1851	newvar = onewnc.createVariable('couple', 'c', ('couple','StrLength'))
1852	basicvardef(newvar, 'couple', 'couples of values', '-')
1853	writing_str_nc(newvar, ['sim','obs'], StringLength)
1854
1855	newvar = onewnc.createVariable('statistics', 'c', ('stats','StrLength'))
1856	basicvardef(newvar, 'statistics', 'statitics from values', '-')
1857	writing_str_nc(newvar, statsn, StringLength)
1858
1859	newvar = onewnc.createVariable('gstatistics', 'c', ('gstats','StrLength'))
1860	basicvardef(newvar, 'gstatistics', 'global statitics from values', '-')
1861	writing_str_nc(newvar, gstatsn, StringLength)
1862
1863	newvar = onewnc.createVariable('tstatistics', 'c', ('tstats','StrLength'))
1864	basicvardef(newvar, 'tstatistics', 'statitics from values along time', '-')
1865	writing_str_nc(newvar, ostatsn, StringLength)
1866
1867	newvar = onewnc.createVariable('ksimstatistics', 'c', ('Nstsim','StrLength'))
1868	basicvardef(newvar, 'ksimstatistics', 'kind of simulated statitics', '-')
1869	writing_str_nc(newvar, Lstdescsim, StringLength)
1870
1871
1872	if obskind == 'trajectory':
1873	if dims.has_key('Z'):
1874	newdim = onewnc.createDimension('trj',3)
1875	else:
1876	newdim = onewnc.createDimension('trj',2)
1877
1878	newvar = onewnc.createVariable('obssimtrj','i',('obstime','trj'))
1879	basicvardef(newvar, 'obssimtrj', 'trajectory on the simulation grid', '-')
1880	newvar[:] = trajpos.transpose()
1881
1882	if dims.has_key('Z'):
1883	newdim = onewnc.createDimension('simtrj',4)
1884	else:
1885	newdim = onewnc.createDimension('simtrj',3)
1886
1887	Nvars = len(valvars)
1888	for ivar in range(Nvars):
1889	simobsvalues = []
1890
1891	varsimobs = valvars[ivar][0] + '_' + valvars[ivar][1]
1892	print ' ' + varsimobs + '... .. .'
1893
1894	if not oobs.variables.has_key(valvars[ivar][1]):
1895	print errormsg
1896	print ' ' + main + ": observations file has not '" + valvars[ivar][1] + \
1897	"' !!"
1898	quit(-1)
1899
1900	if not osim.variables.has_key(valvars[ivar][0]):
1901	if not searchInlist(varNOcheck, valvars[ivar][0]):
1902	print errormsg
1903	print ' ' + main + ": simulation file has not '" + valvars[ivar][0] + \
1904	"' !!"
1905	quit(-1)
1906	else:
1907	ovsim = compute_varNOcheck(osim, valvars[ivar][0])
1908	else:
1909	ovsim = osim.variables[valvars[ivar][0]]
1910
1911	for idn in ovsim.dimensions:
1912	if not searchInlist(simdims.values(),idn):
1913	print errormsg
1914	print ' ' + main + ": dimension '" + idn + "' of variable '" + \
1915	valvars[ivar][0] + "' not provided as reference coordinate [X,Y,Z,T] !!"
1916	quit(-1)
1917
1918	ovobs = oobs.variables[valvars[ivar][1]]
1919	if searchInlist(ovobs.ncattrs(),'_FillValue'):
1920	oFillValue = ovobs.getncattr('_FillValue')
1921	else:
1922	oFillValue = None
1923
1924	# Observed values temporally exact times
1925	Esimobsvalues, EsimobsSvalues, EsimobsTtvalues, trjsim = \
1926	getting_ValidationValues(obskind, Nexactt, dims, trajpos, ovsim, ovobs, \
1927	exacttvalues, oFillValue, Ngrid, 'instantaneous')
1928
1929	# Observed values temporally around coincident times
1930	simobsvalues, simobsSvalues, simobsTtvalues, trjsim = \
1931	getting_ValidationValues(obskind, Ncoindt, dims, trajpos, ovsim, ovobs, \
1932	coindtvalues, oFillValue, Ngrid, 'tbackwardSmean')
1933
1934	# Re-arranging values...
1935	Earrayvals = np.array(Esimobsvalues)
1936	arrayvals = np.array(simobsvalues)
1937	if len(valvars[ivar]) > 2:
1938	const=np.float(valvars[ivar][3])
1939	if valvars[ivar][2] == 'sumc':
1940	EsimobsSvalues = EsimobsSvalues + const
1941	Earrayvals[:,0] = Earrayvals[:,0] + const
1942	simobsSvalues = simobsSvalues + const
1943	arrayvals[:,0] = arrayvals[:,0] + const
1944	elif valvars[ivar][2] == 'subc':
1945	EsimobsSvalues = EsimobsSvalues - const
1946	Earrayvals[:,0] = Earrayvals[:,0] - const
1947	simobsSvalues = simobsSvalues - const
1948	arrayvals[:,0] = arrayvals[:,0] - const
1949	elif valvars[ivar][2] == 'mulc':
1950	EsimobsSvalues = EsimobsSvalues * const
1951	Earrayvals[:,0] = Earrayvals[:,0] * const
1952	simobsSvalues = simobsSvalues * const
1953	arrayvals[:,0] = arrayvals[:,0] * const
1954	elif valvars[ivar][2] == 'divc':
1955	EsimobsSvalues = EsimobsSvalues / const
1956	Earrayvals[:,0] = Earrayvals[:,0] / const
1957	simobsSvalues = simobsSvalues / const
1958	arrayvals[:,0] = arrayvals[:,0] / const
1959	else:
1960	print errormsg
1961	print ' ' + fname + ": operation '" + valvars[ivar][2] + "' not ready!!"
1962	quit(-1)
1963
1964	# statisics sim
1965	simstats = np.zeros((2,5), dtype=np.float)
1966	simstats[0,0] = np.min(Earrayvals[:,0])
1967	simstats[0,1] = np.max(Earrayvals[:,0])
1968	simstats[0,2] = np.mean(Earrayvals[:,0])
1969	simstats[0,3] = np.mean(Earrayvals[:,0]*Earrayvals[:,0])
1970	simstats[0,4] = np.sqrt(simstats[0,3] - simstats[0,2]*simstats[0,2])
1971
1972	simstats[1,0] = np.min(arrayvals[:,0])
1973	simstats[1,1] = np.max(arrayvals[:,0])
1974	simstats[1,2] = np.mean(arrayvals[:,0])
1975	simstats[1,3] = np.mean(arrayvals[:,0]*arrayvals[:,0])
1976	simstats[1,4] = np.sqrt(simstats[1,3] - simstats[1,2]*simstats[1,2])
1977
1978	# statisics obs
1979	# Masking 'nan'
1980	Eobsmask0 = np.where(Earrayvals[:,1] != Earrayvals[:,1], fillValueF,
1981	Earrayvals[:,1])
1982	obsmask0 = np.where(arrayvals[:,1] != arrayvals[:,1], fillValueF, arrayvals[:,1])
1983
1984	Eobsmask = ma.masked_equal(Eobsmask0, fillValueF)
1985	Eobsmask2 = Eobsmask*Eobsmask
1986	obsmask = ma.masked_equal(obsmask0, fillValueF)
1987	obsmask2 = obsmask*obsmask
1988
1989	obsstats = np.zeros((2,5), dtype=np.float)
1990	obsstats[0,0] = obsmask.min()
1991	obsstats[0,1] = obsmask.max()
1992	obsstats[0,2] = obsmask.mean()
1993	obsstats[0,3] = obsmask2.mean()
1994	obsstats[0,4] = np.sqrt(obsstats[0,3] - obsstats[0,2]*obsstats[0,2])
1995
1996	obsstats[1,0] = obsmask.min()
1997	obsstats[1,1] = obsmask.max()
1998	obsstats[1,2] = obsmask.mean()
1999	obsstats[1,3] = obsmask2.mean()
2000	obsstats[1,4] = np.sqrt(obsstats[1,3] - obsstats[1,2]*obsstats[1,2])
2001
2002	# Statistics sim-obs
2003	simobsstats = np.zeros((2,13), dtype=np.float)
2004	Ediffvals = np.zeros((Nexactt), dtype=np.float)
2005	diffvals = np.zeros((Ncoindt), dtype=np.float)
2006
2007	Ediffvals = Earrayvals[:,0] - Eobsmask
2008	diffvals = arrayvals[:,0] - obsmask
2009
2010	Ediff2vals = Ediffvals*Ediffvals
2011	Esumdiff = Ediffvals.sum()
2012	Esumdiff2 = Ediff2vals.sum()
2013
2014	diff2vals = diffvals*diffvals
2015	sumdiff = diffvals.sum()
2016	sumdiff2 = diff2vals.sum()
2017
2018	simobsstats[0,0] = simstats[0,0] - obsstats[0,0]
2019	simobsstats[0,1] = np.mean(Earrayvals[:,0]*Eobsmask)
2020	simobsstats[0,2] = Ediffvals.min()
2021	simobsstats[0,3] = Ediffvals.max()
2022	simobsstats[0,4] = Ediffvals.mean()
2023	simobsstats[0,5] = np.abs(Ediffvals).mean()
2024	simobsstats[0,6] = np.sqrt(Ediff2vals.mean())
2025	simobsstats[0,7], simobsstats[0,8] = sts.mstats.pearsonr(Earrayvals[:,0], \
2026	Eobsmask)
2027
2028	simobsstats[1,0] = simstats[1,0] - obsstats[1,0]
2029	simobsstats[1,1] = np.mean(arrayvals[:,0]*obsmask)
2030	simobsstats[1,2] = diffvals.min()
2031	simobsstats[1,3] = diffvals.max()
2032	simobsstats[1,4] = diffvals.mean()
2033	simobsstats[1,5] = np.abs(diffvals).mean()
2034	simobsstats[1,6] = np.sqrt(diff2vals.mean())
2035	simobsstats[1,7], simobsstats[1,8] = sts.mstats.pearsonr(arrayvals[:,0], \
2036	obsmask)
2037	# From:
2038	#Willmott, C. J. 1981. 'On the validation of models. Physical Geography', 2, 184-194
2039	#Willmott, C. J. (1984). 'On the evaluation of model performance in physical
2040	# geography'. Spatial Statistics and Models, G. L. Gaile and C. J. Willmott, eds.,
2041	# 443-460
2042	#Willmott, C. J., S. G. Ackleson, R. E. Davis, J. J. Feddema, K. M. Klink, D. R.
2043	# Legates, J. O'Donnell, and C. M. Rowe (1985), 'Statistics for the Evaluation and
2044	# Comparison of Models', J. Geophys. Res., 90(C5), 8995-9005
2045	#Legates, D. R., and G. J. McCabe Jr. (1999), 'Evaluating the Use of "Goodness-of-Fit"
2046	# Measures in Hydrologic and Hydroclimatic Model Validation', Water Resour. Res.,
2047	# 35(1), 233-241
2048	#
2049	# Deviation of residuals (SDR)
2050	simobsstats[0,9] = np.mean(np.sqrt(np.abs((Ediffvals-simobsstats[0,0])*(Ediffvals-\
2051	obsstats[0,0]))))
2052	simobsstats[1,9] = np.mean(np.sqrt(np.abs((diffvals-simobsstats[1,0])*(diffvals- \
2053	obsstats[1,0]))))
2054	# Index of Efficiency (IOE)
2055	Eobsbias2series = (Eobsmask - obsstats[0,0])*(Eobsmask - obsstats[0,0])
2056	Esumobsbias2series = Eobsbias2series.sum()
2057	obsbias2series = (obsmask - obsstats[1,0])*(obsmask - obsstats[1,0])
2058	sumobsbias2series = obsbias2series.sum()
2059
2060	simobsstats[0,10] = 1. - Esumdiff2/(Esumobsbias2series)
2061	simobsstats[1,10] = 1. - sumdiff2/(sumobsbias2series)
2062	# Index of Agreement (IOA)
2063	Esimbias2series = Earrayvals[:,0] - obsstats[0,0]
2064	simbias2series = arrayvals[:,0] - obsstats[1,0]
2065
2066	Eobsbias2series = Eobsmask - obsstats[0,0]
2067	obsbias2series = obsmask - obsstats[1,0]
2068
2069	Eabssimbias2series = np.abs(Esimbias2series)
2070	abssimbias2series = np.abs(simbias2series)
2071	Eabsobsbias2series = np.where(Eobsbias2series<0, -Eobsbias2series, Eobsbias2series)
2072	absobsbias2series = np.where(obsbias2series<0, -obsbias2series, obsbias2series)
2073
2074	Eabssimobsbias2series=(Eabssimbias2series+Eabsobsbias2series)*(Eabssimbias2series+\
2075	Eabsobsbias2series)
2076	abssimobsbias2series = (abssimbias2series+absobsbias2series)*(abssimbias2series +\
2077	absobsbias2series)
2078
2079	simobsstats[0,11] = 1. - Esumdiff2/(Eabssimobsbias2series.sum())
2080	simobsstats[1,11] = 1. - sumdiff2/(abssimobsbias2series.sum())
2081	# Fractional Mean Bias (FMB)
2082	simobsstats[0,12]=(simstats[0,0]-obsstats[0,0])/(0.5*(simstats[0,0]+obsstats[0,0]))
2083	simobsstats[1,12]=(simstats[1,0]-obsstats[1,0])/(0.5*(simstats[1,0]+obsstats[1,0]))
2084
2085	# Statistics around sim values exact
2086	Earoundstats = np.zeros((5,Nexactt), dtype=np.float)
2087	for it in range(Nexactt):
2088	Earoundstats[0,it] = np.min(EsimobsSvalues[it,])
2089	Earoundstats[1,it] = np.max(EsimobsSvalues[it,])
2090	Earoundstats[2,it] = np.mean(EsimobsSvalues[it,])
2091	Earoundstats[3,it] = np.mean(EsimobsSvalues[it,]*EsimobsSvalues[it,])
2092	Earoundstats[4,it] = np.sqrt(Earoundstats[3,it] - Earoundstats[2,it]* \
2093	Earoundstats[2,it])
2094
2095	# Statistics around sim values between
2096	aroundstats = np.zeros((5,Ncoindt), dtype=np.float)
2097	for it in range(Ncoindt):
2098	aroundstats[0,it] = np.min(simobsSvalues[it,])
2099	aroundstats[1,it] = np.max(simobsSvalues[it,])
2100	aroundstats[2,it] = np.mean(simobsSvalues[it,])
2101	aroundstats[3,it] = np.mean(simobsSvalues[it,]*simobsSvalues[it,])
2102	aroundstats[4,it] = np.sqrt(aroundstats[3,it] - aroundstats[2,it]* \
2103	aroundstats[2,it])
2104
2105	# exact sim Values to netCDF
2106	newvar = onewnc.createVariable(valvars[ivar][0] + '_Esim', 'f', ('time'), \
2107	fill_value=fillValueF)
2108	descvar = 'exact time simulated: ' + valvars[ivar][0]
2109	basicvardef(newvar, valvars[ivar][0], descvar, ovobs.getncattr('units'))
2110	print 'Lluis NEarratvals:',len(Earrayvals[:,0])
2111	newvar[:] = Earrayvals[:,0]
2112
2113	# exact obs Values to netCDF
2114	newvar = onewnc.createVariable(valvars[ivar][1] + '_Eobs', 'f', ('time'), \
2115	fill_value=fillValueF)
2116	descvar = 'exact time observed: ' + valvars[ivar][1]
2117	basicvardef(newvar, valvars[ivar][1], descvar, ovobs.getncattr('units'))
2118	newvar[:] = Earrayvals[:,1]
2119
2120	# between sim Values to netCDF
2121	newvar = onewnc.createVariable(valvars[ivar][0] + '_sim', 'f', ('betweentime'), \
2122	fill_value=fillValueF)
2123	descvar = 'between observed time simulated: ' + valvars[ivar][0]
2124	basicvardef(newvar, valvars[ivar][0], descvar, ovobs.getncattr('units'))
2125	newvar[:] = arrayvals[:,0]
2126
2127	# between obs Values to netCDF
2128	newvar = onewnc.createVariable(valvars[ivar][1] + '_obs', 'f', ('obstime'), \
2129	fill_value=fillValueF)
2130	descvar = 'observed: ' + valvars[ivar][1]
2131	basicvardef(newvar, valvars[ivar][1], descvar, ovobs.getncattr('units'))
2132	newvar[:] = arrayvals[:,1]
2133
2134	# Around values exact
2135	if not onewnc.variables.has_key(valvars[ivar][0] + 'Earound'):
2136	if dims.has_key('Z'):
2137	if not onewnc.dimensions.has_key('zaround'):
2138	newdim = onewnc.createDimension('zaround',Ngrid*2+1)
2139	newvar = onewnc.createVariable(valvars[ivar][0] + 'Earound', 'f', \
2140	('time','zaround','yaround','xaround'), fill_value=fillValueF)
2141	else:
2142	newvar = onewnc.createVariable(valvars[ivar][0] + 'Earound', 'f', \
2143	('time','yaround','xaround'), fill_value=fillValueF)
2144
2145	descvar = 'exact around simulated values +/- grid values: ' + valvars[ivar][0]
2146	basicvardef(newvar, varsimobs + 'Earound', descvar, ovobs.getncattr('units'))
2147	newvar[:] = EsimobsSvalues
2148
2149	# Around values between
2150	if not onewnc.variables.has_key(valvars[ivar][0] + 'around'):
2151	if dims.has_key('Z'):
2152	if not onewnc.dimensions.has_key('zaround'):
2153	newdim = onewnc.createDimension('zaround',Ngrid*2+1)
2154	newvar = onewnc.createVariable(valvars[ivar][0] + 'around', 'f', \
2155	('betweentime','zaround','yaround','xaround'),fill_value=fillValueF)
2156	else:
2157	newvar = onewnc.createVariable(valvars[ivar][0] + 'around', 'f', \
2158	('betweentime','yaround','xaround'), fill_value=fillValueF)
2159
2160	descvar = 'around simulated values +/- grid values: ' + valvars[ivar][0]
2161	basicvardef(newvar, varsimobs + 'around', descvar, ovobs.getncattr('units'))
2162	newvar[:] = simobsSvalues
2163
2164	# sim Statistics
2165	if not searchInlist(onewnc.variables,valvars[ivar][0] + 'stsim'):
2166	newvar = onewnc.createVariable(valvars[ivar][0] + 'stsim', 'f', ('Nstsim', \
2167	'stats'), fill_value=fillValueF)
2168	descvar = 'simulated statisitcs: ' + valvars[ivar][0]
2169	basicvardef(newvar, valvars[ivar][0] + 'stsim', descvar, \
2170	ovobs.getncattr('units'))
2171	newvar[:] = simstats
2172
2173	# obs Statistics
2174	if not searchInlist(onewnc.variables,valvars[ivar][1] + 'stobs'):
2175	newvar = onewnc.createVariable(valvars[ivar][1] + 'stobs', 'f', ('Nstsim', \
2176	'stats'), fill_value=fillValueF)
2177	descvar = 'observed statisitcs: ' + valvars[ivar][1]
2178	basicvardef(newvar, valvars[ivar][1] + 'stobs', descvar, \
2179	ovobs.getncattr('units'))
2180	newvar[:] = obsstats
2181
2182	# sim-obs Statistics
2183	if not searchInlist(onewnc.variables,varsimobs + 'st'):
2184	newvar = onewnc.createVariable(varsimobs + 'st', 'f', ('Nstsim', 'gstats'), \
2185	fill_value=fillValueF)
2186	descvar = 'simulated-observed tatisitcs: ' + varsimobs
2187	basicvardef(newvar, varsimobs + 'st', descvar, ovobs.getncattr('units'))
2188	newvar[:] = simobsstats
2189
2190	# around sim Statistics exact
2191	if not searchInlist(onewnc.variables,valvars[ivar][0] + 'Estaround'):
2192	newvar = onewnc.createVariable(valvars[ivar][0] + 'Estaround', 'f', \
2193	('time','stats'), fill_value=fillValueF)
2194	descvar = 'exact around (' + str(Ngrid) + ', ' + str(Ngrid) + \
2195	') simulated statisitcs: ' + valvars[ivar][0]
2196	basicvardef(newvar, valvars[ivar][0] + 'Estaround', descvar, \
2197	ovobs.getncattr('units'))
2198	set_attribute(newvar, 'cell_methods', 'Etime_bnds')
2199	newvar[:] = Earoundstats.transpose()
2200
2201	if not searchInlist(onewnc.variables, 'Etime_bnds'):
2202	newvar = onewnc.createVariable('Etime_bnds','f8',('time','bnds'))
2203	basicvardef(newvar, 'Etime_bnds', 'betweentime', obstunits )
2204	set_attribute(newvar, 'calendar', 'standard')
2205	newvar[:] = EsimobsTtvalues
2206
2207	# around obs Statistics exact
2208	if not searchInlist(onewnc.variables,valvars[ivar][1] + 'Estaround'):
2209	newvar = onewnc.createVariable(valvars[ivar][1] + 'Estaround', 'f', \
2210	('time','tstats'), fill_value=fillValueF)
2211	descvar = 'exact around temporal observed statisitcs: ' + valvars[ivar][1]
2212	basicvardef(newvar, valvars[ivar][1] + 'Estaround', descvar, \
2213	ovobs.getncattr('units'))
2214	set_attribute(newvar, 'cell_methods', 'Etime_bnds')
2215
2216	newvar[:] = aroundostats.transpose()
2217
2218	# around sim Statistics between
2219	if not searchInlist(onewnc.variables,valvars[ivar][0] + 'staround'):
2220	newvar = onewnc.createVariable(valvars[ivar][0] + 'staround', 'f', \
2221	('betweentime','stats'), fill_value=fillValueF)
2222	descvar = 'between around (' + str(Ngrid) + ', ' + str(Ngrid) + \
2223	') simulated statisitcs: ' + valvars[ivar][0]
2224	basicvardef(newvar, valvars[ivar][0] + 'staround', descvar, \
2225	ovobs.getncattr('units'))
2226	set_attribute(newvar, 'cell_methods', 'time_bnds')
2227	newvar[:] = aroundstats.transpose()
2228
2229	if not searchInlist(onewnc.variables, 'time_bnds'):
2230	newvar = onewnc.createVariable('time_bnds','f8',('betweentime','bnds'))
2231	basicvardef(newvar, 'time_bnds', 'betweentime', obstunits )
2232	set_attribute(newvar, 'calendar', 'standard')
2233	newvar[:] = simobsTtvalues
2234
2235	# around obs Statistics between
2236	if not searchInlist(onewnc.variables,valvars[ivar][1] + 'staround'):
2237	newvar = onewnc.createVariable(valvars[ivar][1] + 'staround', 'f', \
2238	('betweentime','tstats'), fill_value=fillValueF)
2239	descvar = 'betweem around temporal observed statisitcs: ' + valvars[ivar][1]
2240	basicvardef(newvar, valvars[ivar][1] + 'staround', descvar, \
2241	ovobs.getncattr('units'))
2242	set_attribute(newvar, 'cell_methods', 'time_bnds')
2243	newvar[:] = aroundostats.transpose()
2244
2245	onewnc.sync()
2246
2247	if trjsim is not None:
2248	newvar = onewnc.createVariable('simtrj','i',('betweentime','simtrj'))
2249	basicvardef(newvar,'simtrj','coordinates [X,Y,Z,T] of the coincident ' + \
2250	'trajectory in sim', obstunits)
2251	newvar[:] = trjsim.transpose()
2252
2253	# Adding three variables with the station location, longitude, latitude and height
2254	if obskind == 'single-station':
2255	adding_station_desc(onewnc,stationdesc)
2256
2257	# Global attributes
2258	##
2259	set_attribute(onewnc,'author_nc','Lluis Fita')
2260	set_attribute(onewnc,'institution_nc','Laboratoire de Meteorology Dynamique, ' + \
2261	'LMD-Jussieu, UPMC, Paris')
2262	set_attribute(onewnc,'country_nc','France')
2263	set_attribute(onewnc,'script_nc',main)
2264	set_attribute(onewnc,'version_script',version)
2265	set_attribute(onewnc,'information', \
2266	'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html')
2267	set_attribute(onewnc,'simfile',opts.fsim)
2268	set_attribute(onewnc,'obsfile',opts.fobs)
2269
2270	onewnc.sync()
2271	onewnc.close()
2272
2273	print main + ": successfull writting of '" + ofile + "' !!"

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

Download in other formats: