Context Navigation

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

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

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

Download in other formats: