Context Navigation

validation_sim.py @ 1784

Last change on this file since 1784 was 610, checked in by lfita, 9 years ago
Fixing wrong calculation of statistics Adding statistics of the values `between'
File size: 97.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 = '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 CFtimes_datetime_NOfile(times, units, calendar):
97	""" Provide date/time array from velues of netCDF CF-compilant time variable
98	times= time values
99	units= CF units of the variable time
100	output:
101	array(dimt, 0) = year
102	array(dimt, 1) = month
103	array(dimt, 2) = day
104	array(dimt, 3) = hour
105	array(dimt, 4) = minute
106	array(dimt, 5) = second
107	"""
108	import datetime as dt
109	fname = 'CFtimes_datetime_NOfile'
110
111	txtunits = units.split(' ')
112	tunits = txtunits[0]
113	Srefdate = txtunits[len(txtunits) - 1]
114	# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
115	##
116	timeval = Srefdate.find(':')
117
118	if not timeval == -1:
119	# print ' refdate with time!'
120	refdate = datetimeStr_datetime(txtunits[len(txtunits) - 2] + '_' + Srefdate)
121	else:
122	refdate = dateStr_date(Srefdate)
123
124	dimt = len(times)
125	realdates = np.zeros((dimt, 6), dtype=int)
126
127	secsDay=3600*24.
128
129	# Checking calendar!
130	##
131	y360 = False
132	daycal360 = ['earth_360d', '360d', '360days', '360_day']
133	if searchInlist(daycal360,calendar):
134	print warnmsg
135	print ' ' + fname + ': calendar of 12 months of 30 days !!'
136	y360 = True
137
138	## Not in timedelta
139	# if tunits == 'years':
140	# for it in range(dimt):
141	# realdate = refdate + dt.timedelta(years=float(times[it]))
142	# realdates[it] = int(realdate.year)
143	# elif tunits == 'months':
144	# for it in range(dimt):
145	# realdate = refdate + dt.timedelta(months=float(times[it]))
146	# realdates[it] = int(realdate.year)
147	if y360:
148	if tunits == 'weeks':
149	for it in range(dimt):
150	deltat = dt.timedelta(weeks=float(times[it]))
151	Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
152	realdates[it,:] = date_juliandate(refdate.year,Tsecs)
153	elif tunits == 'days':
154	for it in range(dimt):
155	deltat = dt.timedelta(days=float(times[it]))
156	Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
157	realdates[it,:] = date_juliandate(refdate.year,Tsecs)
158	elif tunits == 'hours':
159	for it in range(dimt):
160	realdate = dt.timedelta(hours=float(times[it]))
161	Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
162	realdates[it,:] = date_juliandate(refdate.year,Tsecs)
163	elif tunits == 'minutes':
164	for it in range(dimt):
165	realdate = dt.timedelta(minutes=float(times[it]))
166	Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
167	realdates[it,:] = date_juliandate(refdate.year,Tsecs)
168	elif tunits == 'seconds':
169	for it in range(dimt):
170	realdate = dt.timedelta(seconds=float(times[it]))
171	Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
172	realdates[it,:] = date_juliandate(refdate.year,Tsecs)
173	elif tunits == 'miliseconds':
174	for it in range(dimt):
175	realdate = dt.timedelta(miliseconds=float(times[it]))
176	Tsecs = deltat.days*secsDay + deltat.seconds + deltat.microseconds/1000.
177	realdates[it,:] = date_juliandate(refdate.year,Tsecs)
178	else:
179	print errormsg
180	print ' CFtimes_datetime: time units "' + tunits + '" not ready!!!!'
181	quit(-1)
182	else:
183	if tunits == 'weeks':
184	for it in range(dimt):
185	realdate = refdate + dt.timedelta(weeks=float(times[it]))
186	realdates[it,0] = int(realdate.year)
187	realdates[it,1] = int(realdate.month)
188	realdates[it,2] = int(realdate.day)
189	realdates[it,3] = int(realdate.hour)
190	realdates[it,4] = int(realdate.second)
191	realdates[it,5] = int(realdate.minute)
192	elif tunits == 'days':
193	for it in range(dimt):
194	realdate = refdate + dt.timedelta(days=float(times[it]))
195	realdates[it,0] = int(realdate.year)
196	realdates[it,1] = int(realdate.month)
197	realdates[it,2] = int(realdate.day)
198	realdates[it,3] = int(realdate.hour)
199	realdates[it,4] = int(realdate.second)
200	realdates[it,5] = int(realdate.minute)
201	elif tunits == 'hours':
202	for it in range(dimt):
203	realdate = refdate + dt.timedelta(hours=float(times[it]))
204	realdates[it,0] = int(realdate.year)
205	realdates[it,1] = int(realdate.month)
206	realdates[it,2] = int(realdate.day)
207	realdates[it,3] = int(realdate.hour)
208	realdates[it,4] = int(realdate.second)
209	realdates[it,5] = int(realdate.minute)
210	elif tunits == 'minutes':
211	for it in range(dimt):
212	realdate = refdate + dt.timedelta(minutes=float(times[it]))
213	realdates[it,0] = int(realdate.year)
214	realdates[it,1] = int(realdate.month)
215	realdates[it,2] = int(realdate.day)
216	realdates[it,3] = int(realdate.hour)
217	realdates[it,4] = int(realdate.second)
218	realdates[it,5] = int(realdate.minute)
219	elif tunits == 'seconds':
220	for it in range(dimt):
221	realdate = refdate + dt.timedelta(seconds=float(times[it]))
222	realdates[it,0] = int(realdate.year)
223	realdates[it,1] = int(realdate.month)
224	realdates[it,2] = int(realdate.day)
225	realdates[it,3] = int(realdate.hour)
226	realdates[it,4] = int(realdate.second)
227	realdates[it,5] = int(realdate.minute)
228	elif tunits == 'milliseconds':
229	for it in range(dimt):
230	realdate = refdate + dt.timedelta(milliseconds=float(times[it]))
231	realdates[it,0] = int(realdate.year)
232	realdates[it,1] = int(realdate.month)
233	realdates[it,2] = int(realdate.day)
234	realdates[it,3] = int(realdate.hour)
235	realdates[it,4] = int(realdate.second)
236	realdates[it,5] = int(realdate.minute)
237	else:
238	print errormsg
239	print ' CFtimes_datetime: time units "' + tunits + '" not ready!!!!'
240	quit(-1)
241
242	return realdates
243
244	def index_3mat(matA,matB,matC,val):
245	""" Function to provide the coordinates of a given value inside three matrix simultaneously
246	index_mat(matA,matB,matC,val)
247	matA= matrix with one set of values
248	matB= matrix with the other set of values
249	matB= matrix with the third set of values
250	val= triplet of values to search
251	>>> 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])
252	[2 1 1]
253	"""
254	fname = 'index_3mat'
255
256	matAshape = matA.shape
257	matBshape = matB.shape
258	matCshape = matC.shape
259
260	for idv in range(len(matAshape)):
261	if matAshape[idv] != matBshape[idv]:
262	print errormsg
263	print ' ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv], \
264	'and B:',matBshape[idv],'does not coincide!!'
265	quit(-1)
266	if matAshape[idv] != matCshape[idv]:
267	print errormsg
268	print ' ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv], \
269	'and C:',matCshape[idv],'does not coincide!!'
270	quit(-1)
271
272	minA = np.min(matA)
273	maxA = np.max(matA)
274	minB = np.min(matB)
275	maxB = np.max(matB)
276	minC = np.min(matC)
277	maxC = np.max(matC)
278
279	if val[0] < minA or val[0] > maxA:
280	print warnmsg
281	print ' ' + fname + ': first value:',val[0],'outside matA range',minA,',', \
282	maxA,'!!'
283	if val[1] < minB or val[1] > maxB:
284	print warnmsg
285	print ' ' + fname + ': second value:',val[1],'outside matB range',minB,',', \
286	maxB,'!!'
287	if val[2] < minC or val[2] > maxC:
288	print warnmsg
289	print ' ' + fname + ': second value:',val[2],'outside matC range',minC,',', \
290	maxC,'!!'
291
292	dist = np.zeros(tuple(matAshape), dtype=np.float)
293	dist = np.sqrt((matA - np.float(val[0]))2 + (matB - np.float(val[1]))2 + \
294	(matC - np.float(val[2]))**2)
295
296	mindist = np.min(dist)
297
298	matlist = list(dist.flatten())
299	ifound = matlist.index(mindist)
300
301	Ndims = len(matAshape)
302	valpos = np.zeros((Ndims), dtype=int)
303	baseprevdims = np.zeros((Ndims), dtype=int)
304
305	for dimid in range(Ndims):
306	baseprevdims[dimid] = np.product(matAshape[dimid+1:Ndims])
307	if dimid == 0:
308	alreadyplaced = 0
309	else:
310	alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
311	valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
312
313	return valpos
314
315	def index_2mat(matA,matB,val):
316	""" Function to provide the coordinates of a given value inside two matrix simultaneously
317	index_mat(matA,matB,val)
318	matA= matrix with one set of values
319	matB= matrix with the pother set of values
320	val= couple of values to search
321	>>> index_2mat(np.arange(27).reshape(3,3,3),np.arange(100,127).reshape(3,3,3),[22,111])
322	[2 1 1]
323	"""
324	fname = 'index_2mat'
325
326	matAshape = matA.shape
327	matBshape = matB.shape
328
329	for idv in range(len(matAshape)):
330	if matAshape[idv] != matBshape[idv]:
331	print errormsg
332	print ' ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv], \
333	'and B:',matBshape[idv],'does not coincide!!'
334	quit(-1)
335
336	minA = np.min(matA)
337	maxA = np.max(matA)
338	minB = np.min(matB)
339	maxB = np.max(matB)
340
341	Ndims = len(matAshape)
342	# valpos = np.ones((Ndims), dtype=int)*-1.
343	valpos = np.zeros((Ndims), dtype=int)
344
345	if val[0] < minA or val[0] > maxA:
346	print warnmsg
347	print ' ' + fname + ': first value:',val[0],'outside matA range',minA,',', \
348	maxA,'!!'
349	return valpos
350	if val[1] < minB or val[1] > maxB:
351	print warnmsg
352	print ' ' + fname + ': second value:',val[1],'outside matB range',minB,',', \
353	maxB,'!!'
354	return valpos
355
356	dist = np.zeros(tuple(matAshape), dtype=np.float)
357	dist = np.sqrt((matA - np.float(val[0]))2 + (matB - np.float(val[1]))2)
358
359	mindist = np.min(dist)
360
361	if mindist != mindist:
362	print ' ' + fname + ': wrong minimal distance',mindist,'!!'
363	return valpos
364	else:
365	matlist = list(dist.flatten())
366	ifound = matlist.index(mindist)
367
368	baseprevdims = np.zeros((Ndims), dtype=int)
369	for dimid in range(Ndims):
370	baseprevdims[dimid] = np.product(matAshape[dimid+1:Ndims])
371	if dimid == 0:
372	alreadyplaced = 0
373	else:
374	alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
375	valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
376
377	return valpos
378
379	def index_mat(matA,val):
380	""" Function to provide the coordinates of a given value inside a matrix
381	index_mat(matA,val)
382	matA= matrix with one set of values
383	val= couple of values to search
384	>>> index_mat(np.arange(27),22.3)
385	22
386	"""
387	fname = 'index_mat'
388
389	matAshape = matA.shape
390
391	minA = np.min(matA)
392	maxA = np.max(matA)
393
394	Ndims = len(matAshape)
395	# valpos = np.ones((Ndims), dtype=int)*-1.
396	valpos = np.zeros((Ndims), dtype=int)
397
398	if val < minA or val > maxA:
399	print warnmsg
400	print ' ' + fname + ': first value:',val,'outside matA range',minA,',', \
401	maxA,'!!'
402	return valpos
403
404	dist = np.zeros(tuple(matAshape), dtype=np.float)
405	dist = (matA - np.float(val))**2
406
407	mindist = np.min(dist)
408	if mindist != mindist:
409	print ' ' + fname + ': wrong minimal distance',mindist,'!!'
410	return valpos
411
412	matlist = list(dist.flatten())
413	valpos = matlist.index(mindist)
414
415	return valpos
416
417	def index_mat_exact(mat,val):
418	""" Function to provide the coordinates of a given exact value inside a matrix
419	index_mat(mat,val)
420	mat= matrix with values
421	val= value to search
422	>>> index_mat(np.arange(27).reshape(3,3,3),22)
423	[2 1 1]
424	"""
425
426	fname = 'index_mat'
427
428	matshape = mat.shape
429
430	matlist = list(mat.flatten())
431	ifound = matlist.index(val)
432
433	Ndims = len(matshape)
434	valpos = np.zeros((Ndims), dtype=int)
435	baseprevdims = np.zeros((Ndims), dtype=int)
436
437	for dimid in range(Ndims):
438	baseprevdims[dimid] = np.product(matshape[dimid+1:Ndims])
439	if dimid == 0:
440	alreadyplaced = 0
441	else:
442	alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
443	valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
444
445	return valpos
446
447	def datetimeStr_datetime(StringDT):
448	""" Function to transform a string date ([YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format) to a date object
449	>>> datetimeStr_datetime('1976-02-17_00:00:00')
450	1976-02-17 00:00:00
451	"""
452	import datetime as dt
453
454	fname = 'datetimeStr_datetime'
455
456	dateD = np.zeros((3), dtype=int)
457	timeT = np.zeros((3), dtype=int)
458
459	dateD[0] = int(StringDT[0:4])
460	dateD[1] = int(StringDT[5:7])
461	dateD[2] = int(StringDT[8:10])
462
463	trefT = StringDT.find(':')
464	if not trefT == -1:
465	# print ' ' + fname + ': refdate with time!'
466	timeT[0] = int(StringDT[11:13])
467	timeT[1] = int(StringDT[14:16])
468	timeT[2] = int(StringDT[17:19])
469
470	if int(dateD[0]) == 0:
471	print warnmsg
472	print ' ' + fname + ': 0 reference year!! changing to 1'
473	dateD[0] = 1
474
475	newdatetime = dt.datetime(dateD[0], dateD[1], dateD[2], timeT[0], timeT[1], timeT[2])
476
477	return newdatetime
478
479	def datetimeStr_conversion(StringDT,typeSi,typeSo):
480	""" Function to transform a string date to an another date object
481	StringDT= string with the date and time
482	typeSi= type of datetime string input
483	typeSo= type of datetime string output
484	[typeSi/o]
485	'cfTime': [time],[units]; ]time in CF-convention format [units] = [tunits] since [refdate]
486	'matYmdHMS': numerical vector with [[YYYY], [MM], [DD], [HH], [MI], [SS]]
487	'YmdHMS': [YYYY][MM][DD][HH][MI][SS] format
488	'Y-m-d_H:M:S': [YYYY]-[MM]-[DD]_[HH]:[MI]:[SS] format
489	'Y-m-d H:M:S': [YYYY]-[MM]-[DD] [HH]:[MI]:[SS] format
490	'Y/m/d H-M-S': [YYYY]/[MM]/[DD] [HH]-[MI]-[SS] format
491	'WRFdatetime': [Y], [Y], [Y], [Y], '-', [M], [M], '-', [D], [D], '_', [H],
492	[H], ':', [M], [M], ':', [S], [S]
493	>>> datetimeStr_conversion('1976-02-17_08:32:05','Y-m-d_H:M:S','matYmdHMS')
494	[1976 2 17 8 32 5]
495	>>> datetimeStr_conversion(str(137880)+',minutes since 1979-12-01_00:00:00','cfTime','Y/m/d H-M-S')
496	1980/03/05 18-00-00
497	"""
498	import datetime as dt
499
500	fname = 'datetimeStr_conversion'
501
502	if StringDT[0:1] == 'h':
503	print fname + '_____________________________________________________________'
504	print datetimeStr_conversion.__doc__
505	quit()
506
507	if typeSi == 'cfTime':
508	timeval = np.float(StringDT.split(',')[0])
509	tunits = StringDT.split(',')[1].split(' ')[0]
510	Srefdate = StringDT.split(',')[1].split(' ')[2]
511
512	# Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
513	##
514	yrref=Srefdate[0:4]
515	monref=Srefdate[5:7]
516	dayref=Srefdate[8:10]
517
518	trefT = Srefdate.find(':')
519	if not trefT == -1:
520	# print ' ' + fname + ': refdate with time!'
521	horref=Srefdate[11:13]
522	minref=Srefdate[14:16]
523	secref=Srefdate[17:19]
524	refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref + \
525	'_' + horref + ':' + minref + ':' + secref)
526	else:
527	refdate = datetimeStr_datetime( yrref + '-' + monref + '-' + dayref + \
528	+ '_00:00:00')
529
530	if tunits == 'weeks':
531	newdate = refdate + dt.timedelta(weeks=float(timeval))
532	elif tunits == 'days':
533	newdate = refdate + dt.timedelta(days=float(timeval))
534	elif tunits == 'hours':
535	newdate = refdate + dt.timedelta(hours=float(timeval))
536	elif tunits == 'minutes':
537	newdate = refdate + dt.timedelta(minutes=float(timeval))
538	elif tunits == 'seconds':
539	newdate = refdate + dt.timedelta(seconds=float(timeval))
540	elif tunits == 'milliseconds':
541	newdate = refdate + dt.timedelta(milliseconds=float(timeval))
542	else:
543	print errormsg
544	print ' timeref_datetime: time units "' + tunits + '" not ready!!!!'
545	quit(-1)
546
547	yr = newdate.year
548	mo = newdate.month
549	da = newdate.day
550	ho = newdate.hour
551	mi = newdate.minute
552	se = newdate.second
553	elif typeSi == 'matYmdHMS':
554	yr = StringDT[0]
555	mo = StringDT[1]
556	da = StringDT[2]
557	ho = StringDT[3]
558	mi = StringDT[4]
559	se = StringDT[5]
560	elif typeSi == 'YmdHMS':
561	yr = int(StringDT[0:4])
562	mo = int(StringDT[4:6])
563	da = int(StringDT[6:8])
564	ho = int(StringDT[8:10])
565	mi = int(StringDT[10:12])
566	se = int(StringDT[12:14])
567	elif typeSi == 'Y-m-d_H:M:S':
568	dateDT = StringDT.split('_')
569	dateD = dateDT[0].split('-')
570	timeT = dateDT[1].split(':')
571	yr = int(dateD[0])
572	mo = int(dateD[1])
573	da = int(dateD[2])
574	ho = int(timeT[0])
575	mi = int(timeT[1])
576	se = int(timeT[2])
577	elif typeSi == 'Y-m-d H:M:S':
578	dateDT = StringDT.split(' ')
579	dateD = dateDT[0].split('-')
580	timeT = dateDT[1].split(':')
581	yr = int(dateD[0])
582	mo = int(dateD[1])
583	da = int(dateD[2])
584	ho = int(timeT[0])
585	mi = int(timeT[1])
586	se = int(timeT[2])
587	elif typeSi == 'Y/m/d H-M-S':
588	dateDT = StringDT.split(' ')
589	dateD = dateDT[0].split('/')
590	timeT = dateDT[1].split('-')
591	yr = int(dateD[0])
592	mo = int(dateD[1])
593	da = int(dateD[2])
594	ho = int(timeT[0])
595	mi = int(timeT[1])
596	se = int(timeT[2])
597	elif typeSi == 'WRFdatetime':
598	yr = int(StringDT[0])1000 + int(StringDT[1])100 + int(StringDT[2])*10 + \
599	int(StringDT[3])
600	mo = int(StringDT[5])*10 + int(StringDT[6])
601	da = int(StringDT[8])*10 + int(StringDT[9])
602	ho = int(StringDT[11])*10 + int(StringDT[12])
603	mi = int(StringDT[14])*10 + int(StringDT[15])
604	se = int(StringDT[17])*10 + int(StringDT[18])
605	else:
606	print errormsg
607	print ' ' + fname + ': type of String input date "' + typeSi + \
608	'" not ready !!!!'
609	quit(-1)
610
611	if typeSo == 'matYmdHMS':
612	dateYmdHMS = np.zeros((6), dtype=int)
613	dateYmdHMS[0] = yr
614	dateYmdHMS[1] = mo
615	dateYmdHMS[2] = da
616	dateYmdHMS[3] = ho
617	dateYmdHMS[4] = mi
618	dateYmdHMS[5] = se
619	elif typeSo == 'YmdHMS':
620	dateYmdHMS = str(yr).zfill(4) + str(mo).zfill(2) + str(da).zfill(2) + \
621	str(ho).zfill(2) + str(mi).zfill(2) + str(se).zfill(2)
622	elif typeSo == 'Y-m-d_H:M:S':
623	dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' + \
624	str(da).zfill(2) + '_' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \
625	str(se).zfill(2)
626	elif typeSo == 'Y-m-d H:M:S':
627	dateYmdHMS = str(yr).zfill(4) + '-' + str(mo).zfill(2) + '-' + \
628	str(da).zfill(2) + ' ' + str(ho).zfill(2) + ':' + str(mi).zfill(2) + ':' + \
629	str(se).zfill(2)
630	elif typeSo == 'Y/m/d H-M-S':
631	dateYmdHMS = str(yr).zfill(4) + '/' + str(mo).zfill(2) + '/' + \
632	str(da).zfill(2) + ' ' + str(ho).zfill(2) + '-' + str(mi).zfill(2) + '-' + \
633	str(se).zfill(2)
634	elif typeSo == 'WRFdatetime':
635	dateYmdHMS = []
636	yM = yr/1000
637	yC = (yr-yM*1000)/100
638	yD = (yr-yM1000-yC100)/10
639	yU = yr-yM1000-yC100-yD*10
640
641	mD = mo/10
642	mU = mo-mD*10
643
644	dD = da/10
645	dU = da-dD*10
646
647	hD = ho/10
648	hU = ho-hD*10
649
650	miD = mi/10
651	miU = mi-miD*10
652
653	sD = se/10
654	sU = se-sD*10
655
656	dateYmdHMS.append(str(yM))
657	dateYmdHMS.append(str(yC))
658	dateYmdHMS.append(str(yD))
659	dateYmdHMS.append(str(yU))
660	dateYmdHMS.append('-')
661	dateYmdHMS.append(str(mD))
662	dateYmdHMS.append(str(mU))
663	dateYmdHMS.append('-')
664	dateYmdHMS.append(str(dD))
665	dateYmdHMS.append(str(dU))
666	dateYmdHMS.append('_')
667	dateYmdHMS.append(str(hD))
668	dateYmdHMS.append(str(hU))
669	dateYmdHMS.append(':')
670	dateYmdHMS.append(str(miD))
671	dateYmdHMS.append(str(miU))
672	dateYmdHMS.append(':')
673	dateYmdHMS.append(str(sD))
674	dateYmdHMS.append(str(sU))
675	else:
676	print errormsg
677	print ' ' + fname + ': type of output date "' + typeSo + '" not ready !!!!'
678	quit(-1)
679
680	return dateYmdHMS
681
682	def coincident_CFtimes(tvalB, tunitA, tunitB):
683	""" Function to make coincident times for two different sets of CFtimes
684	tvalB= time values B
685	tunitA= time units times A to which we want to make coincidence
686	tunitB= time units times B
687	>>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00',
688	'hours since 1949-12-01 00:00:00')
689	[ 0. 3600. 7200. 10800. 14400. 18000. 21600. 25200. 28800. 32400.]
690	>>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00',
691	'hours since 1979-12-01 00:00:00')
692	[ 9.46684800e+08 9.46688400e+08 9.46692000e+08 9.46695600e+08
693	9.46699200e+08 9.46702800e+08 9.46706400e+08 9.46710000e+08
694	9.46713600e+08 9.46717200e+08]
695	"""
696	import datetime as dt
697	fname = 'coincident_CFtimes'
698
699	trefA = tunitA.split(' ')[2] + ' ' + tunitA.split(' ')[3]
700	trefB = tunitB.split(' ')[2] + ' ' + tunitB.split(' ')[3]
701	tuA = tunitA.split(' ')[0]
702	tuB = tunitB.split(' ')[0]
703
704	if tuA != tuB:
705	if tuA == 'microseconds':
706	if tuB == 'microseconds':
707	tB = tvalB*1.
708	elif tuB == 'seconds':
709	tB = tvalB*10.e6
710	elif tuB == 'minutes':
711	tB = tvalB60.10.e6
712	elif tuB == 'hours':
713	tB = tvalB3600.10.e6
714	elif tuB == 'days':
715	tB = tvalB3600.24.*10.e6
716	else:
717	print errormsg
718	print ' ' + fname + ": combination of time untis: '" + tuA + \
719	"' & '" + tuB + "' not ready !!"
720	quit(-1)
721	elif tuA == 'seconds':
722	if tuB == 'microseconds':
723	tB = tvalB/10.e6
724	elif tuB == 'seconds':
725	tB = tvalB*1.
726	elif tuB == 'minutes':
727	tB = tvalB*60.
728	elif tuB == 'hours':
729	tB = tvalB*3600.
730	elif tuB == 'days':
731	tB = tvalB3600.24.
732	else:
733	print errormsg
734	print ' ' + fname + ": combination of time untis: '" + tuA + \
735	"' & '" + tuB + "' not ready !!"
736	quit(-1)
737	elif tuA == 'minutes':
738	if tuB == 'microseconds':
739	tB = tvalB/(60.*10.e6)
740	elif tuB == 'seconds':
741	tB = tvalB/60.
742	elif tuB == 'minutes':
743	tB = tvalB*1.
744	elif tuB == 'hours':
745	tB = tvalB*60.
746	elif tuB == 'days':
747	tB = tvalB60.24.
748	else:
749	print errormsg
750	print ' ' + fname + ": combination of time untis: '" + tuA + \
751	"' & '" + tuB + "' not ready !!"
752	quit(-1)
753	elif tuA == 'hours':
754	if tuB == 'microseconds':
755	tB = tvalB/(3600.*10.e6)
756	elif tuB == 'seconds':
757	tB = tvalB/3600.
758	elif tuB == 'minutes':
759	tB = tvalB/60.
760	elif tuB == 'hours':
761	tB = tvalB*1.
762	elif tuB == 'days':
763	tB = tvalB*24.
764	else:
765	print errormsg
766	print ' ' + fname + ": combination of time untis: '" + tuA + \
767	"' & '" + tuB + "' not ready !!"
768	quit(-1)
769	elif tuA == 'days':
770	if tuB == 'microseconds':
771	tB = tvalB/(24.3600.10.e6)
772	elif tuB == 'seconds':
773	tB = tvalB/(24.*3600.)
774	elif tuB == 'minutes':
775	tB = tvalB/(24.*60.)
776	elif tuB == 'hours':
777	tB = tvalB/24.
778	elif tuB == 'days':
779	tB = tvalB*1.
780	else:
781	print errormsg
782	print ' ' + fname + ": combination of time untis: '" + tuA + \
783	"' & '" + tuB + "' not ready !!"
784	quit(-1)
785	else:
786	print errormsg
787	print ' ' + fname + ": time untis: '" + tuA + "' not ready !!"
788	quit(-1)
789	else:
790	tB = tvalB*1.
791
792	if trefA != trefB:
793	trefTA = dt.datetime.strptime(trefA, '%Y-%m-%d %H:%M:%S')
794	trefTB = dt.datetime.strptime(trefB, '%Y-%m-%d %H:%M:%S')
795
796	difft = trefTB - trefTA
797	diffv = difft.days24.3600. + difft.seconds + difft.microseconds/1.e6
798	print ' ' + fname + ': different reference refA:',trefTA,'refB',trefTB
799	print ' difference:',difft,':',diffv,'(in seconds)'
800
801	if tuA == 'microseconds':
802	tB = tB + diffv*1.e6
803	elif tuA == 'seconds':
804	tB = tB + diffv
805	elif tuA == 'minutes':
806	tB = tB + diffv/(60.)
807	elif tuA == 'hours':
808	tB = tB + diffv/(3600.)
809	elif tuA == 'days':
810	tB = tB + diffv/(24.*3600.)
811	else:
812	print errormsg
813	print ' ' + fname + ": time untis: '" + tuA + "' not ready !!"
814	quit(-1)
815
816	return tB
817
818	def slice_variable(varobj, dimslice):
819	""" Function to return a slice of a given variable according to values to its
820	dimensions
821	slice_variable(varobj, dimslice)
822	varobj= object wit the variable
823	dimslice= [[dimname1]:[value1]\|[[dimname2]:[value2], ...] pairs of dimension
824	[value]:
825	* [integer]: which value of the dimension
826	* -1: all along the dimension
827	* -9: last value of the dimension
828	* [beg]@[end] slice from [beg] to [end]
829	"""
830	fname = 'slice_variable'
831
832	if varobj == 'h':
833	print fname + '_____________________________________________________________'
834	print slice_variable.__doc__
835	quit()
836
837	vardims = varobj.dimensions
838	Ndimvar = len(vardims)
839
840	Ndimcut = len(dimslice.split('\|'))
841	dimsl = dimslice.split('\|')
842
843	varvalsdim = []
844	dimnslice = []
845
846	for idd in range(Ndimvar):
847	for idc in range(Ndimcut):
848	dimcutn = dimsl[idc].split(':')[0]
849	dimcutv = dimsl[idc].split(':')[1]
850	if vardims[idd] == dimcutn:
851	posfrac = dimcutv.find('@')
852	if posfrac != -1:
853	inifrac = int(dimcutv.split('@')[0])
854	endfrac = int(dimcutv.split('@')[1])
855	varvalsdim.append(slice(inifrac,endfrac))
856	dimnslice.append(vardims[idd])
857	else:
858	if int(dimcutv) == -1:
859	varvalsdim.append(slice(0,varobj.shape[idd]))
860	dimnslice.append(vardims[idd])
861	elif int(dimcutv) == -9:
862	varvalsdim.append(int(varobj.shape[idd])-1)
863	else:
864	varvalsdim.append(int(dimcutv))
865	break
866
867	varvalues = varobj[tuple(varvalsdim)]
868
869	return varvalues, dimnslice
870
871	def func_compute_varNOcheck(ncobj, varn):
872	""" Function to compute variables which are not originary in the file
873	ncobj= netCDF object file
874	varn = variable to compute:
875	'WRFdens': air density from WRF variables
876	'WRFght': geopotential height from WRF variables
877	'WRFp': pressure from WRF variables
878	'WRFrh': relative humidty fom WRF variables
879	'WRFt': temperature from WRF variables
880	'WRFwds': surface wind direction from WRF variables
881	'WRFwss': surface wind speed from WRF variables
882	'WRFz': height from WRF variables
883	"""
884	fname = 'compute_varNOcheck'
885
886	if varn == 'WRFdens':
887	# print ' ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \
888	# 'DNW)/g ...'
889	grav = 9.81
890
891	# Just we need in in absolute values: Size of the central grid cell
892	## dxval = ncobj.getncattr('DX')
893	## dyval = ncobj.getncattr('DY')
894	## mapfac = ncobj.variables['MAPFAC_M'][:]
895	## area = dxvaldyvalmapfac
896	dimensions = ncobj.variables['MU'].dimensions
897
898	mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
899	dnw = ncobj.variables['DNW'][:]
900
901	varNOcheckv = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \
902	dtype=np.float)
903	levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
904
905	for it in range(mu.shape[0]):
906	for iz in range(dnw.shape[1]):
907	levval.fill(np.abs(dnw[it,iz]))
908	varNOcheck[it,iz,:,:] = levval
909	varNOcheck[it,iz,:,:] = mu[it,:,:]*varNOcheck[it,iz,:,:]/grav
910
911	elif varn == 'WRFght':
912	# print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
913	varNOcheckv = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
914	dimensions = ncobj.variables['PH'].dimensions
915
916	elif varn == 'WRFp':
917	# print ' ' + fname + ': Retrieving pressure value from WRF as P + PB'
918	varNOcheckv = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
919	dimensions = ncobj.variables['P'].dimensions
920
921	elif varn == 'WRFrh':
922	# print ' ' + main + ": computing relative humidity from WRF as 'Tetens'" +\
923	# ' equation (T,P) ...'
924	p0=100000.
925	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
926	tk = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
927	qv = ncobj.variables['QVAPOR'][:]
928
929	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
930	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
931
932	varNOcheckv = qv/data2
933	dimensions = ncobj.variables['P'].dimensions
934
935	elif varn == 'WRFt':
936	# print ' ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
937	p0=100000.
938	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
939
940	varNOcheckv = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
941	dimensions = ncobj.variables['T'].dimensions
942
943	elif varn == 'WRFwds':
944	# print ' ' + main + ': computing surface wind direction from WRF as ATAN2(V,U) ...'
945	varNOcheckv = np.arctan2(ncobj.variables['V10'][:], ncobj.variables['U10'][:])
946	dimensions = ncobj.variables['V10'].dimensions
947
948	elif varn == 'WRFwss':
949	# print ' ' + main + ': computing surface wind speed from WRF as SQRT(U2 + V2) ...'
950	varNOcheckv = np.sqrt(ncobj.variables['U10'][:]*ncobj.variables['U10'][:] + \
951	ncobj.variables['V10'][:]*ncobj.variables['V10'][:])
952	dimensions = ncobj.variables['U10'].dimensions
953
954	elif varn == 'WRFz':
955	grav = 9.81
956	# print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
957	varNOcheckv = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/grav
958	dimensions = ncobj.variables['PH'].dimensions
959
960	else:
961	print erromsg
962	print ' ' + fname + ": variable '" + varn + "' nor ready !!"
963	quit(-1)
964
965	return varNOcheck
966
967	class compute_varNOcheck(object):
968	""" Class to compute variables which are not originary in the file
969	ncobj= netCDF object file
970	varn = variable to compute:
971	'WRFdens': air density from WRF variables
972	'WRFght': geopotential height from WRF variables
973	'WRFp': pressure from WRF variables
974	'WRFrh': relative humidty fom WRF variables
975	'TSrhs': surface relative humidty fom TS variables
976	'WRFrhs': surface relative humidty fom WRF variables
977	'WRFT': CF-time from WRF variables
978	'WRFt': temperature from WRF variables
979	'TStd': dew-point temperature from TS variables
980	'WRFtd': dew-point temperature from WRF variables
981	'WRFwd': wind direction from WRF variables
982	'TSwds': surface wind direction from TS variables
983	'WRFwds': surface wind direction from WRF variables
984	'WRFws': wind speed from WRF variables
985	'TSwss': surface wind speed from TS variables
986	'WRFwss': surface wind speed from WRF variables
987	'WRFz': height from WRF variables
988	"""
989	fname = 'compute_varNOcheck'
990
991	def __init__(self, ncobj, varn):
992
993	if ncobj is None:
994	self = None
995	self.dimensions = None
996	self.shape = None
997	self.__values = None
998	else:
999	if varn == 'WRFdens':
1000	# print ' ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \
1001	# 'DNW)/g ...'
1002	grav = 9.81
1003
1004	# Just we need in in absolute values: Size of the central grid cell
1005	## dxval = ncobj.getncattr('DX')
1006	## dyval = ncobj.getncattr('DY')
1007	## mapfac = ncobj.variables['MAPFAC_M'][:]
1008	## area = dxvaldyvalmapfac
1009	dimensions = ncobj.variables['MU'].dimensions
1010	shape = ncobj.variables['MU'].shape
1011
1012	mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
1013	dnw = ncobj.variables['DNW'][:]
1014
1015	varNOcheckv = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \
1016	dtype=np.float)
1017	levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
1018
1019	for it in range(mu.shape[0]):
1020	for iz in range(dnw.shape[1]):
1021	levval.fill(np.abs(dnw[it,iz]))
1022	varNOcheck[it,iz,:,:] = levval
1023	varNOcheck[it,iz,:,:] = mu[it,:,:]*varNOcheck[it,iz,:,:]/grav
1024
1025	elif varn == 'WRFght':
1026	# print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
1027	varNOcheckv = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
1028	dimensions = ncobj.variables['PH'].dimensions
1029	shape = ncobj.variables['PH'].shape
1030
1031	elif varn == 'WRFp':
1032	# print ' ' + fname + ': Retrieving pressure value from WRF as P + PB'
1033	varNOcheckv = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
1034	dimensions = ncobj.variables['P'].dimensions
1035	shape = ncobj.variables['P'].shape
1036
1037	elif varn == 'WRFrh':
1038	# print ' ' + main + ": computing relative humidity from WRF as 'Tetens'" +\
1039	# ' equation (T,P) ...'
1040	p0=100000.
1041	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
1042	tk = (ncobj.variables['T'][:])(p/p0)*(2./7.)
1043	qv = ncobj.variables['QVAPOR'][:]
1044
1045	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
1046	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
1047
1048	varNOcheckv = qv/data2
1049	dimensions = ncobj.variables['P'].dimensions
1050	shape = ncobj.variables['P'].shape
1051
1052	elif varn == 'TSrhs':
1053	# print ' ' + main + ": computing surface relative humidity from TSs as 'Tetens'" +\
1054	# ' equation (T,P) ...'
1055	p0=100000.
1056	p=ncobj.variables['psfc'][:]
1057	tk = (ncobj.variables['t'][:])(p/p0)*(2./7.)
1058	qv = ncobj.variables['q'][:]
1059
1060	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
1061	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
1062
1063	varNOcheckv = qv/data2
1064	dimensions = ncobj.variables['psfc'].dimensions
1065	shape = ncobj.variables['psfc'].shape
1066
1067	elif varn == 'WRFrhs':
1068	# print ' ' + main + ": computing surface relative humidity from WRF as 'Tetens'" +\
1069	# ' equation (T,P) ...'
1070	p0=100000.
1071	p=ncobj.variables['PSFC'][:]
1072	tk = (ncobj.variables['T2'][:] + 300.)(p/p0)*(2./7.)
1073	qv = ncobj.variables['Q2'][:]
1074
1075	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
1076	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
1077
1078	varNOcheckv = qv/data2
1079	dimensions = ncobj.variables['PSFC'].dimensions
1080	shape = ncobj.variables['PSFC'].shape
1081
1082	elif varn == 'WRFT':
1083	# To compute CF-times from WRF kind
1084	#
1085	import datetime as dt
1086
1087	times = ncobj.variables['Times']
1088	dimt = times.shape[0]
1089	varNOcheckv = np.zeros((dimt), dtype=np.float64)
1090	self.unitsval = 'seconds since 1949-12-01 00:00:00'
1091	refdate = datetimeStr_datetime('1949-12-01_00:00:00')
1092
1093	dimensions = tuple([ncobj.variables['Times'].dimensions[0]])
1094	shape = tuple([dimt])
1095
1096	for it in range(dimt):
1097	datevalS = datetimeStr_conversion(times[it,:], 'WRFdatetime', \
1098	'YmdHMS')
1099	dateval = dt.datetime.strptime(datevalS, '%Y%m%d%H%M%S')
1100	difft = dateval - refdate
1101	varNOcheckv[it] = difft.days3600.24. + difft.seconds + \
1102	np.float(int(difft.microseconds/10.e6))
1103
1104	elif varn == 'WRFt':
1105	# print ' ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
1106	p0=100000.
1107	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
1108
1109	varNOcheckv = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
1110	dimensions = ncobj.variables['T'].dimensions
1111	shape = ncobj.variables['P'].shape
1112
1113	elif varn == 'TStd':
1114	# print ' ' + main + ': computing dew-point temperature from TS as t and Tetens...'
1115	# tacking from: http://en.wikipedia.org/wiki/Dew_point
1116	p=ncobj.variables['psfc'][:]
1117
1118	temp = ncobj.variables['t'][:]
1119
1120	qv = ncobj.variables['q'][:]
1121
1122	tk = temp
1123	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
1124	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
1125
1126	rh = qv/data2
1127
1128	pa = rh * data1
1129	varNOcheckv = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
1130
1131	dimensions = ncobj.variables['t'].dimensions
1132	shape = ncobj.variables['t'].shape
1133
1134	elif varn == 'WRFtd':
1135	# print ' ' + main + ': computing dew-point temperature from WRF as inv_potT(T + 300) and Tetens...'
1136	# tacking from: http://en.wikipedia.org/wiki/Dew_point
1137	p0=100000.
1138	p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
1139
1140	temp = (ncobj.variables['T'][:] + 300.)(p/p0)*(2./7.)
1141
1142	qv = ncobj.variables['QVAPOR'][:]
1143
1144	tk = temp - 273.15
1145	data1 = 10.0.6112np.exp(17.67*(tk-273.16)/(tk-29.65))
1146	data2 = 0.622data1/(0.01p-(1.-0.622)*data1)
1147
1148	rh = qv/data2
1149
1150	pa = rh * data1
1151	varNOcheckv = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
1152
1153	dimensions = ncobj.variables['T'].dimensions
1154	shape = ncobj.variables['P'].shape
1155
1156	elif varn == 'WRFwd':
1157	# print ' ' + main + ': computing wind direction from WRF as ATAN2PI(V,U) ...'
1158	uwind = ncobj.variables['U'][:]
1159	vwind = ncobj.variables['V'][:]
1160	dx = uwind.shape[3]
1161	dy = vwind.shape[2]
1162
1163	# de-staggering
1164	ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
1165	va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])
1166
1167	theta = np.arctan2(ua, va)
1168	theta = np.where(theta < 0., theta + 2.*np.pi, theta)
1169	varNOcheckv = 360.theta/(2.np.pi)
1170
1171	dimensions = tuple(['Time','bottom_top','south_north','west_east'])
1172	shape = ua.shape
1173
1174	elif varn == 'WRFwds':
1175	# print ' ' + main + ': computing surface wind direction from WRF as ATAN2(V,U) ...'
1176	theta = np.arctan2(ncobj.variables['V10'][:], ncobj.variables['U10'][:])
1177	theta = np.where(theta < 0., theta + 2.*np.pi, theta)
1178
1179	varNOcheckv = 360.theta/(2.np.pi)
1180	dimensions = ncobj.variables['V10'].dimensions
1181	shape = ncobj.variables['V10'].shape
1182
1183	elif varn == 'TSwds':
1184	# print ' ' + main + ': computing surface wind direction from TSs as ATAN2(v,u) ...'
1185	theta = np.arctan2(ncobj.variables['v'][:], ncobj.variables['u'][:])
1186	theta = np.where(theta < 0., theta + 2.*np.pi, theta)
1187
1188	varNOcheckv = 360.theta/(2.np.pi)
1189	dimensions = ncobj.variables['v'].dimensions
1190	shape = ncobj.variables['v'].shape
1191
1192	elif varn == 'WRFws':
1193	# print ' ' + main + ': computing wind speed from WRF as SQRT(U2 + V2) ...'
1194	uwind = ncobj.variables['U'][:]
1195	vwind = ncobj.variables['V'][:]
1196	dx = uwind.shape[3]
1197	dy = vwind.shape[2]
1198
1199	# de-staggering
1200	ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
1201	va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])
1202
1203	varNOcheckv = np.sqrt(uaua + vava)
1204	dimensions = tuple(['Time','bottom_top','south_north','west_east'])
1205	shape = ua.shape
1206
1207	elif varn == 'TSwss':
1208	# print ' ' + main + ': computing surface wind speed from TSs as SQRT(u2 + v2) ...'
1209	varNOcheckv = np.sqrt(ncobj.variables['u'][:]* \
1210	ncobj.variables['u'][:] + ncobj.variables['v'][:]* \
1211	ncobj.variables['v'][:])
1212	dimensions = ncobj.variables['u'].dimensions
1213	shape = ncobj.variables['u'].shape
1214
1215	elif varn == 'WRFwss':
1216	# print ' ' + main + ': computing surface wind speed from WRF as SQRT(U2 + V2) ...'
1217	varNOcheckv = np.sqrt(ncobj.variables['U10'][:]* \
1218	ncobj.variables['U10'][:] + ncobj.variables['V10'][:]* \
1219	ncobj.variables['V10'][:])
1220	dimensions = ncobj.variables['U10'].dimensions
1221	shape = ncobj.variables['U10'].shape
1222
1223	elif varn == 'WRFz':
1224	grav = 9.81
1225	# print ' ' + main + ': computing geopotential height from WRF as PH + PHB ...'
1226	varNOcheckv = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/grav
1227	dimensions = ncobj.variables['PH'].dimensions
1228	shape = ncobj.variables['PH'].shape
1229
1230	else:
1231	print errormsg
1232	print ' ' + fname + ": variable '" + varn + "' nor ready !!"
1233	quit(-1)
1234
1235	self.dimensions = dimensions
1236	self.shape = shape
1237	self.__values = varNOcheckv
1238
1239	def __getitem__(self,elem):
1240	return self.__values[elem]
1241
1242	def adding_station_desc(onc,stdesc):
1243	""" Function to add a station description in a netCDF file
1244	onc= netCDF object
1245	stdesc= station description lon, lat, height
1246	"""
1247	fname = 'adding_station_desc'
1248
1249	newvar = onc.createVariable( 'station', 'c', ('StrLength'))
1250	newvar[0:len(stdesc[0])] = stdesc[0]
1251
1252	newdim = onc.createDimension('nst',1)
1253
1254	if onc.variables.has_key('lon'):
1255	print warnmsg
1256	print ' ' + fname + ": variable 'lon' already exist !!"
1257	print " renaming it as 'lonst'"
1258	lonname = 'lonst'
1259	else:
1260	lonname = 'lon'
1261
1262	newvar = onc.createVariable( lonname, 'f4', ('nst'))
1263	basicvardef(newvar, lonname, 'longitude', 'degrees_West' )
1264	newvar[:] = stdesc[1]
1265
1266	if onc.variables.has_key('lat'):
1267	print warnmsg
1268	print ' ' + fname + ": variable 'lat' already exist !!"
1269	print " renaming it as 'latst'"
1270	latname = 'latst'
1271	else:
1272	latname = 'lat'
1273
1274	newvar = onc.createVariable( latname, 'f4', ('nst'))
1275	basicvardef(newvar, lonname, 'latitude', 'degrees_North' )
1276	newvar[:] = stdesc[2]
1277
1278	if onc.variables.has_key('height'):
1279	print warnmsg
1280	print ' ' + fname + ": variable 'height' already exist !!"
1281	print " renaming it as 'heightst'"
1282	heightname = 'heightst'
1283	else:
1284	heightname = 'height'
1285
1286	newvar = onc.createVariable( heightname, 'f4', ('nst'))
1287	basicvardef(newvar, heightname, 'height above sea level', 'm' )
1288	newvar[:] = stdesc[3]
1289
1290	return
1291
1292	class Quantiles(object):
1293	""" Class to provide quantiles from a given arrayof values
1294	"""
1295
1296	def __init__(self, values, Nquants):
1297	import numpy.ma as ma
1298
1299	if values is None:
1300	self.quantilesv = None
1301
1302	else:
1303	self.quantilesv = []
1304
1305	vals0 = values.flatten()
1306	Nvalues = len(vals0)
1307	vals = ma.masked_equal(vals0, None)
1308	Nvals=len(vals.compressed())
1309
1310	sortedvals = sorted(vals.compressed())
1311	for iq in range(Nquants):
1312	self.quantilesv.append(sortedvals[int((Nvals-1)*iq/Nquants)])
1313
1314	self.quantilesv.append(sortedvals[Nvals-1])
1315
1316
1317	def getting_ValidationValues(okind, dt, ds, trjpos, ovs, ovo, tvalues, oFill, Ng, kvals):
1318	""" Function to get the values to validate accroding to the type of observation
1319	okind= observational kind
1320	dt= initial number of values to retrieve
1321	ds= dictionary with the names of the dimensions (sim, obs)
1322	trjpos= positions of the multi-stations (t, Y, X) or trajectory ([Z], Y, X)
1323	ovs= object with the values of the simulation
1324	ovs= object with the values of the observations
1325	tvalues= temporal values (sim. time step, obs. time step, sim t value, obs t value, t diff)
1326	oFill= Fill Value for the observations
1327	Ng= number of grid points around the observation
1328	kvals= kind of values
1329	'instantaneous': values are taken as instantaneous values
1330	'tbackwardSmean': simulated values are taken as time averages from back to the point
1331	'tbackwardOmean': observed values are taken as time averages from back to the point
1332	return:
1333	sovalues= simulated values at the observation point and time
1334	soSvalues= values Ngrid points around the simulated point
1335	soTtvalues= inital/ending period between two consecutive obsevations (for `single-station')
1336	trjs= trajectory on the simulation space
1337	"""
1338	fname = 'getting_ValidationValues'
1339
1340	sovalues = []
1341
1342	if kvals == 'instantaneous':
1343	dtf = dt
1344	elif kvals == 'tbackwardSmean':
1345	print ' ' + fname + ':',kvals,'!!'
1346	uniqt = np.unique(tvalues[:,3])
1347	dtf = len(uniqt)
1348	print ' initially we got',dt,'values which will become',dtf
1349	postf = {}
1350	for it in range(dtf):
1351	if it == 0:
1352	postf[uniqt[it]] = [0,0]
1353	elif it == 1:
1354	posprev = postf[uniqt[it-1]][1]
1355	posit = list(tvalues[:,3]).index(uniqt[it])
1356	postf[uniqt[it]] = [posprev, posit+1]
1357	else:
1358	posprev = postf[uniqt[it-1]][1]
1359	posit = list(tvalues[:,3]).index(uniqt[it])
1360	postf[uniqt[it]] = [posprev+1, np.min([posit+1,dt-1])]
1361	elif kvals == 'tbackwardOmean':
1362	print ' ' + fname + ':',kvals,'!!'
1363	uniqt = np.unique(tvalues[:,2])
1364	dtf = len(uniqt)
1365	print ' initially we got',dt,'values which will become',dtf
1366	print ' ==== NOT READY === '
1367	quit(-1)
1368	else:
1369	print errormsg
1370	print ' ' + fname + ": kind of values '" + kvals + "' not ready!!"
1371	quit(-1)
1372
1373	# Simulated values spatially around
1374	if ds.has_key('Z'):
1375	soSvalues = np.zeros((dt, Ng2+1, Ng2+1, Ng*2+1), dtype = np.float)
1376	if okind == 'trajectory':
1377	trjs = np.zeros((4,dt), dtype=int)
1378	else:
1379	trjs = None
1380	else:
1381	soSvalues = np.zeros((dt, Ng2+1, Ng2+1), dtype = np.float)
1382	if okind == 'trajectory':
1383	trjs = np.zeros((3,dt), dtype=int)
1384	else:
1385	trjs = None
1386
1387	if okind == 'single-station':
1388	soTtvalues = np.zeros((dt,2), dtype=np.float)
1389	else:
1390	None
1391
1392	if okind == 'multi-points':
1393	for it in range(dt):
1394	slicev = ds['X'][0] + ':' + str(trjpos[2,it]) + '\|' + ds['Y'][0] + \
1395	':' + str(trjpos[1,it]) + '\|' + ds['T'][0]+ ':' + str(tvalues[it][0])
1396	slicevar, dimslice = slice_variable(ovs, slicev)
1397	sovalues.append([ slicevar, ovo[tvalues[it][1]]])
1398	slicev = ds['X'][0] + ':' + str(trjpos[2,it]-Ng) + '@' + \
1399	str(trjpos[2,it]+Ng) + '\|' + ds['Y'][0] + ':' + \
1400	str(trjpos[1,it]-Ng) + '@' + str(trjpos[1,it]+Ng) + '\|' + \
1401	ds['T'][0]+':'+str(tvalues[it][0])
1402	slicevar, dimslice = slice_variable(ovs, slicev)
1403	soSvalues[it,:,:] = slicevar
1404
1405	elif okind == 'single-station':
1406	for it in range(dt):
1407	ito = int(tvalues[it,1])
1408	if valdimsim.has_key('X') and valdimsim.has_key('Y'):
1409	slicev = ds['X'][0] + ':' + str(stationpos[1]) + '\|' + \
1410	ds['Y'][0] + ':' + str(stationpos[0]) + '\|' + \
1411	ds['T'][0] + ':' + str(int(tvalues[it][0]))
1412	else:
1413	slicev = ds['T'][0] + ':' + str(int(tvalues[it][0]))
1414	slicevar, dimslice = slice_variable(ovs, slicev)
1415	if ovo[int(ito)] == oFill or ovo[int(ito)] == '--':
1416	sovalues.append([ slicevar, fillValueF])
1417	# elif ovo[int(ito)] != ovo[int(ito)]:
1418	# sovalues.append([ slicevar, fillValueF])
1419	else:
1420	sovalues.append([ slicevar, ovo[int(ito)]])
1421	if valdimsim.has_key('X') and valdimsim.has_key('Y'):
1422	slicev = ds['X'][0] + ':' + str(stationpos[1]-Ng) + '@' + \
1423	str(stationpos[1]+Ng+1) + '\|' + ds['Y'][0] + ':' + \
1424	str(stationpos[0]-Ng) + '@' + str(stationpos[0]+Ng+1) + '\|' + \
1425	ds['T'][0] + ':' + str(int(tvalues[it,0]))
1426	else:
1427	slicev = ds['T'][0] + ':' + str(int(tvalues[it][0]))
1428	slicevar, dimslice = slice_variable(ovs, slicev)
1429	soSvalues[it,:,:] = slicevar
1430
1431	if it == 0:
1432	itoi = 0
1433	itof = int(tvalues[it,1]) / 2
1434	elif it == dt-1:
1435	itoi = int( (ito + int(tvalues[it-1,1])) / 2)
1436	itof = int(tvalues[it,1])
1437	else:
1438	itod = int( (ito - int(tvalues[it-1,1])) / 2 )
1439	itoi = ito - itod
1440	itod = int( (int(tvalues[it+1,1]) - ito) / 2 )
1441	itof = ito + itod
1442
1443	soTtvalues[it,0] = valdimobs['T'][itoi]
1444	soTtvalues[it,1] = valdimobs['T'][itof]
1445
1446	elif okind == 'trajectory':
1447	if ds.has_key('Z'):
1448	for it in range(dt):
1449	ito = int(tvalues[it,1])
1450	if notfound[ito] == 0:
1451	trjpos[2,ito] = index_mat(valdimsim['Z'][tvalues[it,0],:, \
1452	trjpos[1,ito],trjpos[0,ito]], valdimobs['Z'][ito])
1453	slicev = ds['X'][0]+':'+str(trjpos[0,ito]) + '\|' + \
1454	ds['Y'][0]+':'+str(trjpos[1,ito]) + '\|' + \
1455	ds['Z'][0]+':'+str(trjpos[2,ito]) + '\|' + \
1456	ds['T'][0]+':'+str(int(tvalues[it,0]))
1457	slicevar, dimslice = slice_variable(ovs, slicev)
1458	sovalues.append([ slicevar, ovo[int(ito)]])
1459	minx = np.max([trjpos[0,ito]-Ng,0])
1460	maxx = np.min([trjpos[0,ito]+Ng+1,ovs.shape[3]])
1461	miny = np.max([trjpos[1,ito]-Ng,0])
1462	maxy = np.min([trjpos[1,ito]+Ng+1,ovs.shape[2]])
1463	minz = np.max([trjpos[2,ito]-Ng,0])
1464	maxz = np.min([trjpos[2,ito]+Ng+1,ovs.shape[1]])
1465
1466	slicev = ds['X'][0] + ':' + str(minx) + '@' + str(maxx) + '\|' + \
1467	ds['Y'][0] + ':' + str(miny) + '@' + str(maxy) + '\|' + \
1468	ds['Z'][0] + ':' + str(minz) + '@' + str(maxz) + '\|' + \
1469	ds['T'][0] + ':' + str(int(tvalues[it,0]))
1470	slicevar, dimslice = slice_variable(ovs, slicev)
1471
1472	sliceS = []
1473	sliceS.append(it)
1474	sliceS.append(slice(0,maxz-minz))
1475	sliceS.append(slice(0,maxy-miny))
1476	sliceS.append(slice(0,maxx-minx))
1477
1478	soSvalues[tuple(sliceS)] = slicevar
1479	if ivar == 0:
1480	trjs[0,it] = trjpos[0,ito]
1481	trjs[1,it] = trjpos[1,ito]
1482	trjs[2,it] = trjpos[2,ito]
1483	trjs[3,it] = tvalues[it,0]
1484	else:
1485	sovalues.append([fillValueF, fillValueF])
1486	soSvalues[it,:,:,:]= np.ones((Ng2+1,Ng2+1,Ng*2+1), \
1487	dtype = np.float)*fillValueF
1488	# 2D trajectory
1489	else:
1490	for it in range(dt):
1491	if notfound[it] == 0:
1492	ito = tvalues[it,1]
1493	slicev = ds['X'][0]+':'+str(trjpos[2,ito]) + '\|' + \
1494	ds['Y'][0]+':'+str(trjpos[1,ito]) + '\|' + \
1495	ds['T'][0]+':'+str(tvalues[ito,0])
1496	slicevar, dimslice = slice_variable(ovs, slicev)
1497	sovalues.append([ slicevar, ovo[tvalues[it,1]]])
1498	slicev = ds['X'][0] + ':' + str(trjpos[0,it]-Ng) + '@' + \
1499	str(trjpos[0,it]+Ng) + '\|' + ds['Y'][0] + ':' + \
1500	str(trjpos[1,it]-Ng) + '@' + str(trjpos[1,it]+Ng) + \
1501	'\|' + ds['T'][0] + ':' + str(tvalues[it,0])
1502	slicevar, dimslice = slice_variable(ovs, slicev)
1503	soSvalues[it,:,:] = slicevar
1504	else:
1505	sovalues.append([fillValue, fillValue])
1506	soSvalues[it,:,:] = np.ones((Ng2+1,Ng2+1), \
1507	dtype = np.float)*fillValueF
1508	print sovalues[varsimobs][:][it]
1509	else:
1510	print errormsg
1511	print ' ' + fname + ": observatino kind '" + okind + "' not ready!!"
1512	quit(-1)
1513
1514	# Re-arranging final values
1515	##
1516	if kvals == 'instantaneous':
1517	return sovalues, soSvalues, soTtvalues, trjs
1518
1519	elif kvals == 'tbackwardSmean':
1520	fsovalues = []
1521	if ds.has_key('Z'):
1522	fsoSvalues = np.zeros((dtf, Ng2+1, Ng2+1, Ng*2+1), dtype = np.float)
1523	if okind == 'trajectory':
1524	ftrjs = np.zeros((4,dtf), dtype=int)
1525	else:
1526	ftrjs = None
1527	else:
1528	fsoSvalues = np.zeros((dtf, Ng2+1, Ng2+1), dtype = np.float)
1529	if okind == 'trajectory':
1530	ftrjs = np.zeros((3,dtf), dtype=int)
1531	else:
1532	ftrjs = None
1533
1534	if okind == 'single-station':
1535	fsoTtvalues = np.ones((dtf,2), dtype=np.float)*fillValueF
1536	else:
1537	None
1538
1539	for it in range(1,dtf):
1540	tv = uniqt[it]
1541	intv = postf[tv]
1542
1543	# Temporal statistics
1544	if len(np.array(sovalues[intv[0]:intv[1]]).shape) != 1:
1545	sovs = np.array(sovalues[intv[0]:intv[1]])[:,0]
1546	minv = np.min(sovs)
1547	maxv = np.max(sovs)
1548	meanv = np.mean(sovs)
1549	stdv = np.std(sovs)
1550
1551	fsovalues.append([meanv, np.array(sovalues[intv[0]:intv[1]])[0,1], \
1552	minv, maxv, stdv])
1553	else:
1554	fsovalues.append([fillValueF, fillValueF, fillValueF, fillValueF, \
1555	fillValueF, fillValueF])
1556	if ds.has_key('Z'):
1557	if okind == 'trajectory':
1558	for ip in range(4):
1559	ftrjs[ip,it] = np.mean(trjs[ip,intv[0]:intv[1]])
1560	for iz in range(2*Ng+1):
1561	for iy in range(2*Ng+1):
1562	for ix in range(2*Ng+1):
1563	fsoSvalues[it,iz,iy,ix] = np.mean(soSvalues[intv[0]: \
1564	intv[1],iz,iy,ix])
1565	else:
1566	if okind == 'trajectory':
1567	for ip in range(3):
1568	ftrjs[ip,it] = np.mean(trjs[ip,intv[0]:intv[1]])
1569	for iy in range(2*Ng+1):
1570	for ix in range(2*Ng+1):
1571	fsoSvalues[it,iy,ix] = np.mean(soSvalues[intv[0]:intv[1], \
1572	iy,ix])
1573	fsoTtvalues[it,0] = soTtvalues[intv[0],0]
1574	fsoTtvalues[it,1] = soTtvalues[intv[1],0]
1575
1576	return fsovalues, fsoSvalues, fsoTtvalues, ftrjs
1577
1578	elif kvals == 'tbackwardOmean':
1579	print ' ' + fname + ':',kvals,'!!'
1580	uniqt = np.unique(tvalues[:,2])
1581	dtf = len(uniqt)
1582	print ' initially we got',dt,'values which will become',dtf
1583
1584	return
1585
1586
1587	####### ###### ##### #### ### ## #
1588
1589	strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " + \
1590	" and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')"
1591
1592	kindobs=['multi-points', 'single-station', 'trajectory']
1593	strkObs="kind of observations; 'multi-points': multiple individual punctual obs " + \
1594	"(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
1595	"'trajectory': following a trajectory"
1596	simopers = ['sumc','subc','mulc','divc']
1597	opersinf = 'sumc,[constant]: add [constant] to variables values; subc,[constant]: '+ \
1598	'substract [constant] to variables values; mulc,[constant]: multipy by ' + \
1599	'[constant] to variables values; divc,[constant]: divide by [constant] to ' + \
1600	'variables values'
1601	varNOcheck = ['WRFdens', 'WRFght', 'WRFp', 'WRFrh', 'TSrhs', 'WRFrhs', 'WRFT', \
1602	'WRFt', 'TStd', 'WRFtd', 'WRFwd', 'TSwds', 'WRFwds', 'WRFws', 'TSwss', 'WRFwss', \
1603	'WRFz']
1604	varNOcheckinf = "'WRFdens': air density from WRF variables; " + \
1605	"'WRFght': geopotentiali height from WRF variables; " + \
1606	"'WRFp': pressure from WRF variables; " + \
1607	"'WRFrh': relative humidty fom WRF variables; " + \
1608	"'TSrhs': surface relative humidity from TS variables; " + \
1609	"'WRFrhs': surface relative humidity from WRF variables; " + \
1610	"'WRFT': CF-time from WRF variables; " + \
1611	"'WRFt': temperature from WRF variables; " + \
1612	"'TStd': dew-point temperature from TS variables; " + \
1613	"'WRFtd': dew-point temperature from WRF variables; " + \
1614	"'WRFwd': wind direction from WRF variables; " + \
1615	"'TSwds': surface wind direction from TS variables; " + \
1616	"'WRFwds': surface wind direction from WRF variables; " + \
1617	"'WRFws': wind speed from WRF variables; " + \
1618	"'TSwss': surface wind speed from TS variables; " + \
1619	"'WRFwss': surface wind speed from WRF variables; " + \
1620	"'WRFz': height from WRF variables"
1621
1622	dimshelp = "[DIM]@[simdim]@[obsdim] ',' list of couples of dimensions names from " + \
1623	"each source ([DIM]='X','Y','Z','T'; None, no value)"
1624	vardimshelp = "[DIM]@[simvardim]@[obsvardim] ',' list of couples of variables " + \
1625	"names with dimensions values from each source ([DIM]='X','Y','Z','T'; None, " + \
1626	"no value, WRFdiagnosted variables also available: " + varNOcheckinf + ")"
1627	varshelp="[simvar]@[obsvar]@[[oper]@[val]] ',' list of couples of variables to " + \
1628	"validate and if necessary operation and value (sim. values) available " + \
1629	"operations: " + opersinf + " (WRFdiagnosted variables also available: " + \
1630	varNOcheckinf + ")"
1631	statsn = ['minimum', 'maximum', 'mean', 'mean2', 'standard deviation']
1632	gstatsn = ['bias', 'simobs_mean', 'sim_obsmin', 'sim_obsmax', 'sim_obsmean', 'mae', \
1633	'rmse', 'r_pearsoncorr', 'p_pearsoncorr', 'deviation_of_residuals_SDR', \
1634	'indef_of_efficiency_IOE', 'index_of_agreement_IOA', 'fractional_mean_bias_FMB']
1635	ostatsn = ['number of points', 'minimum', 'maximum', 'mean', 'mean2', \
1636	'standard deviation']
1637
1638	parser = OptionParser()
1639	parser.add_option("-d", "--dimensions", dest="dims", help=dimshelp, metavar="VALUES")
1640	parser.add_option("-D", "--vardimensions", dest="vardims",
1641	help=vardimshelp, metavar="VALUES")
1642	parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs,
1643	help=strkObs, metavar="FILE")
1644	parser.add_option("-l", "--stationLocation", dest="stloc",
1645	help="name (\| for spaces), longitude, latitude and height of the station (only for 'single-station')",
1646	metavar="FILE")
1647	parser.add_option("-o", "--observation", dest="fobs",
1648	help="observations file to validate", metavar="FILE")
1649	parser.add_option("-s", "--simulation", dest="fsim",
1650	help="simulation file to validate", metavar="FILE")
1651	parser.add_option("-t", "--trajectoryfile", dest="trajf",
1652	help="file with grid points of the trajectory in the simulation grid ('simtrj')",
1653	metavar="FILE")
1654	parser.add_option("-v", "--variables", dest="vars",
1655	help=varshelp, metavar="VALUES")
1656
1657	(opts, args) = parser.parse_args()
1658
1659	####### ###### ##### #### ### ## #
1660	# Number of different statistics according to the temporal coincidence
1661	# 0: Exact time
1662	# 1: Simulation values closest to observed times
1663	# 2: Simulation values between consecutive observed times
1664	Nstsim = 3
1665
1666	stdescsim = ['E', 'C', 'B']
1667	prestdescsim = ['exact', 'closest', 'between']
1668	Lstdescsim = ['exact time', 'cloest time', 'between observational time-steps']
1669
1670	####### #######
1671	## MAIN
1672	#######
1673
1674	ofile='validation_sim.nc'
1675
1676	if opts.dims is None:
1677	print errormsg
1678	print ' ' + main + ': No list of dimensions are provided!!'
1679	print ' a ',' list of values X@[dimxsim]@[dimxobs],...,T@[dimtsim]@[dimtobs]'+\
1680	' is needed'
1681	quit(-1)
1682	else:
1683	simdims = {}
1684	obsdims = {}
1685	print main +': couple of dimensions _______'
1686	dims = {}
1687	ds = opts.dims.split(',')
1688	for d in ds:
1689	dsecs = d.split('@')
1690	if len(dsecs) != 3:
1691	print errormsg
1692	print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
1693	print ' [DIM]@[dimnsim]@[dimnobs]'
1694	quit(-1)
1695	if dsecs[1] != 'None':
1696	dims[dsecs[0]] = [dsecs[1], dsecs[2]]
1697	simdims[dsecs[0]] = dsecs[1]
1698	obsdims[dsecs[0]] = dsecs[2]
1699
1700	print ' ',dsecs[0],':',dsecs[1],',',dsecs[2]
1701
1702	if opts.vardims is None:
1703	print errormsg
1704	print ' ' + main + ': No list of variables with dimension values are provided!!'
1705	print ' a ',' list of values X@[vardimxsim]@[vardimxobs],...,T@' + \
1706	'[vardimtsim]@[vardimtobs] is needed'
1707	quit(-1)
1708	else:
1709	print main +': couple of variable dimensions _______'
1710	vardims = {}
1711	ds = opts.vardims.split(',')
1712	for d in ds:
1713	dsecs = d.split('@')
1714	if len(dsecs) != 3:
1715	print errormsg
1716	print ' ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
1717	print ' [DIM]@[vardimnsim]@[vardimnobs]'
1718	quit(-1)
1719	if dsecs[1] != 'None':
1720	vardims[dsecs[0]] = [dsecs[1], dsecs[2]]
1721	print ' ',dsecs[0],':',dsecs[1],',',dsecs[2]
1722
1723	if opts.obskind is None:
1724	print errormsg
1725	print ' ' + main + ': No kind of observations provided !!'
1726	quit(-1)
1727	else:
1728	obskind = opts.obskind
1729	if obskind == 'single-station':
1730	if opts.stloc is None:
1731	print errormsg
1732	print ' ' + main + ': No station location provided !!'
1733	quit(-1)
1734	else:
1735	stationdesc = [opts.stloc.split(',')[0].replace('\|',' '), \
1736	np.float(opts.stloc.split(',')[1]), np.float(opts.stloc.split(',')[2]),\
1737	np.float(opts.stloc.split(',')[3])]
1738
1739	if opts.fobs is None:
1740	print errormsg
1741	print ' ' + main + ': No observations file is provided!!'
1742	quit(-1)
1743	else:
1744	if not os.path.isfile(opts.fobs):
1745	print errormsg
1746	print ' ' + main + ": observations file '" + opts.fobs + "' does not exist !!"
1747	quit(-1)
1748
1749	if opts.fsim is None:
1750	print errormsg
1751	print ' ' + main + ': No simulation file is provided!!'
1752	quit(-1)
1753	else:
1754	if not os.path.isfile(opts.fsim):
1755	print errormsg
1756	print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
1757	quit(-1)
1758
1759	if opts.vars is None:
1760	print errormsg
1761	print ' ' + main + ': No list of couples of variables is provided!!'
1762	print ' a ',' list of values [varsim]@[varobs],... is needed'
1763	quit(-1)
1764	else:
1765	valvars = []
1766	vs = opts.vars.split(',')
1767	for v in vs:
1768	vsecs = v.split('@')
1769	if len(vsecs) < 2:
1770	print errormsg
1771	print ' ' + main + ': wrong number of values in:',vsecs, \
1772	' at least 2 are needed !!'
1773	print ' [varsim]@[varobs]@[[oper][val]]'
1774	quit(-1)
1775	if len(vsecs) > 2:
1776	if not searchInlist(simopers,vsecs[2]):
1777	print errormsg
1778	print main + ": operation on simulation values '" + vsecs[2] + \
1779	"' not ready !!"
1780	quit(-1)
1781
1782	valvars.append(vsecs)
1783
1784	# Openning observations trajectory
1785	##
1786	oobs = NetCDFFile(opts.fobs, 'r')
1787
1788	valdimobs = {}
1789	for dn in dims:
1790	print dn,':',dims[dn]
1791	if dims[dn][1] != 'None':
1792	if not oobs.dimensions.has_key(dims[dn][1]):
1793	print errormsg
1794	print ' ' + main + ": observations file does not have dimension '" + \
1795	dims[dn][1] + "' !!"
1796	quit(-1)
1797	if vardims[dn][1] != 'None':
1798	if not oobs.variables.has_key(vardims[dn][1]):
1799	print errormsg
1800	print ' ' + main + ": observations file does not have varibale " + \
1801	"dimension '" + vardims[dn][1] + "' !!"
1802	quit(-1)
1803	valdimobs[dn] = oobs.variables[vardims[dn][1]][:]
1804	else:
1805	if dn == 'X':
1806	valdimobs[dn] = stationdesc[1]
1807	elif dn == 'Y':
1808	valdimobs[dn] = stationdesc[2]
1809	elif dn == 'Z':
1810	valdimobs[dn] = stationdesc[3]
1811
1812	osim = NetCDFFile(opts.fsim, 'r')
1813
1814	valdimsim = {}
1815	for dn in dims:
1816	if dims[dn][0] != 'None':
1817	if not osim.dimensions.has_key(dims[dn][0]):
1818	print errormsg
1819	print ' ' + main + ": simulation file '" + opts.fsim + \
1820	"' does not have dimension '" + dims[dn][0] + "' !!"
1821	print ' it has: ',osim.dimensions
1822	quit(-1)
1823
1824	if not osim.variables.has_key(vardims[dn][0]) and \
1825	not searchInlist(varNOcheck,vardims[dn][0]):
1826	print errormsg
1827	print ' ' + main + ": simulation file '" + opts.fsim + \
1828	"' does not have varibale dimension '" + vardims[dn][0] + "' !!"
1829	print ' it has variables:',osim.variables
1830	quit(-1)
1831	if searchInlist(varNOcheck,vardims[dn][0]):
1832	valdimsim[dn] = compute_varNOcheck(osim, vardims[dn][0])
1833	else:
1834	valdimsim[dn] = osim.variables[vardims[dn][0]][:]
1835
1836	# General characteristics
1837	dimtobs = valdimobs['T'].shape[0]
1838	dimtsim = valdimsim['T'].shape[0]
1839
1840	print main +': observational time-steps:',dimtobs,'simulation:',dimtsim
1841
1842	notfound = np.zeros((dimtobs), dtype=int)
1843
1844	if obskind == 'multi-points':
1845	trajpos = np.zeros((2,dimt),dtype=int)
1846	for it in range(dimtobs):
1847	trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'], \
1848	[valdimobs['X'][it],valdimobss['Y'][it]])
1849	elif obskind == 'single-station':
1850	trajpos = None
1851	stationpos = np.zeros((2), dtype=int)
1852	if valdimsim.has_key('X') and valdimsim.has_key('Y'):
1853	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],[valdimobs['Y'], \
1854	valdimobs['X']])
1855	iid = 0
1856	for idn in osim.variables[vardims['X'][0]].dimensions:
1857	if idn == dims['X'][0]:
1858	stationpos[1] = stsimpos[iid]
1859	elif idn == dims['Y'][0]:
1860	stationpos[0] = stsimpos[iid]
1861
1862	iid = iid + 1
1863	print main + ': station point in simulation:', stationpos
1864	print ' station position:',valdimobs['X'],',',valdimobs['Y']
1865	print ' simulation coord.:',valdimsim['X'][tuple(stsimpos)],',', \
1866	valdimsim['Y'][tuple(stsimpos)]
1867	else:
1868	print main + ': validation with two time-series !!'
1869
1870	elif obskind == 'trajectory':
1871	if opts.trajf is not None:
1872	if not os.path.isfile(opts.fsim):
1873	print errormsg
1874	print ' ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
1875	quit(-1)
1876	else:
1877	otrjf = NetCDFFile(opts.fsim, 'r')
1878	trajpos[0,:] = otrjf.variables['obssimtrj'][0]
1879	trajpos[1,:] = otrjf.variables['obssimtrj'][1]
1880	otrjf.close()
1881	else:
1882	if dims.has_key('Z'):
1883	trajpos = np.zeros((3,dimtobs),dtype=int)
1884	for it in range(dimtobs):
1885	if np.mod(it*100./dimtobs,10.) == 0.:
1886	print ' trajectory done: ',it*100./dimtobs,'%'
1887	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \
1888	[valdimobs['Y'][it],valdimobs['X'][it]])
1889	stationpos = np.zeros((2), dtype=int)
1890	iid = 0
1891	for idn in osim.variables[vardims['X'][0]].dimensions:
1892	if idn == dims['X'][0]:
1893	stationpos[1] = stsimpos[iid]
1894	elif idn == dims['Y'][0]:
1895	stationpos[0] = stsimpos[iid]
1896	iid = iid + 1
1897	if stationpos[0] == 0 and stationpos[1] == 0: notfound[it] = 1
1898
1899	trajpos[0,it] = stationpos[0]
1900	trajpos[1,it] = stationpos[1]
1901	# In the simulation 'Z' varies with time ... non-hydrostatic model! ;)
1902	# trajpos[2,it] = index_mat(valdimsim['Z'][it,:,stationpos[0], \
1903	# stationpos[1]], valdimobs['Z'][it])
1904	else:
1905	trajpos = np.zeros((2,dimtobs),dtype=int)
1906	for it in range(dimtobs):
1907	stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'], \
1908	[valdimobs['Y'][it],valdimobss['X'][it]])
1909	stationpos = np.zeros((2), dtype=int)
1910	iid = 0
1911	for idn in osim.variables[vardims['X'][0]].dimensions:
1912	if idn == dims['X'][0]:
1913	stationpos[1] = stsimpos[iid]
1914	elif idn == dims['Y'][0]:
1915	stationpos[0] = stsimpos[iid]
1916	iid = iid + 1
1917	if stationpos[0] == 0 or stationpos[1] == 0: notfound[it] = 1
1918
1919	trajpos[0,it] = stationspos[0]
1920	trajpos[1,it] = stationspos[1]
1921
1922	print main + ': not found',np.sum(notfound),'points of the trajectory'
1923
1924	# Getting times
1925	tobj = oobs.variables[vardims['T'][1]]
1926	obstunits = tobj.getncattr('units')
1927	if vardims['T'][0] == 'WRFT':
1928	tsim = valdimsim['T'][:]
1929	simtunits = 'seconds since 1949-12-01 00:00:00'
1930	else:
1931	tsim = osim.variables[vardims['T'][0]][:]
1932	otsim = osim.variables[vardims['T'][0]]
1933	simtunits = otsim.getncattr('units')
1934
1935	simobstimes = coincident_CFtimes(tsim, obstunits, simtunits)
1936
1937	iobst = CFtimes_datetime_NOfile([valdimobs['T'][0]], obstunits, 'standard')
1938	fobst = CFtimes_datetime_NOfile([valdimobs['T'][dimtobs-1]], obstunits, 'standard')
1939	isimt = CFtimes_datetime_NOfile([simobstimes[0]], obstunits, 'standard')
1940	fsimt = CFtimes_datetime_NOfile([simobstimes[dimtsim-1]], obstunits, 'standard')
1941	print 'check of times _______'
1942	print ' * Observations'
1943	print ' - first time:', iobst
1944	print ' - last time:', fobst
1945	print ' * Simulations'
1946	print ' - first time:', isimt
1947	print ' - last time:', fsimt
1948
1949	#
1950	## Looking for exact/near times
1951	#
1952
1953	# Exact Coincident times
1954	##
1955	exacttvalues0 = []
1956	for it in range(dimtsim):
1957	ot = 0
1958	for ito in range(ot,dimtobs-1):
1959	if valdimobs['T'][ito] == simobstimes[it]:
1960	ot = ito
1961	exacttvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito]])
1962
1963	exacttvalues = np.array(exacttvalues0, dtype=np.float)
1964
1965	if len(exacttvalues) == 0:
1966	print warnmsg
1967	print ' ' + main + ': no exact values found!'
1968	Nexactt = 0
1969	# quit(-1)
1970	else:
1971	Nexactt = len(exacttvalues[:,0])
1972
1973	print main + ': found',Nexactt,'Temporal exact values in simulation and observations'
1974
1975	# Sim Closest times
1976	##
1977	Nsimt = 0
1978	closesttvalues0 = []
1979	closesttvalues0st = []
1980	tsiminit = 0
1981	tsimend = 0
1982
1983	dtsim = simobstimes[1] - simobstimes[0]
1984
1985	for it in range(dimtsim):
1986	ot = 0
1987	for ito in range(ot,dimtobs-1):
1988	if np.abs(valdimobs['T'][ito] - simobstimes[it]) <= dtsim/2.:
1989	ot = ito
1990	tdist = simobstimes[it] - valdimobs['T'][ito]
1991	closesttvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito], \
1992	tdist])
1993	Nsimt = Nsimt + 1
1994	if tsiminit == 0: tsiminit=simobstimes[it]
1995	tsimend = simobstimes[it]
1996
1997	closesttvalues = np.array(closesttvalues0, dtype=np.float)
1998
1999	print 'closesttvales shape:',closesttvalues.shape
2000	Nclosest = len(closesttvalues[:,0])
2001	print main + ': found',Nclosest,'Simulation time-values closest to observations'
2002
2003	if Nclosest == 0:
2004	print warnmsg
2005	print main + ': no cclosest times found !!'
2006	print ' stopping it'
2007	quit(-1)
2008
2009	# Sim Coincident times
2010	##
2011	Nsimt = 0
2012	coindtvalues0 = []
2013	coindtvalues0st = []
2014	tsiminit = 0
2015	tsimend = 0
2016
2017	for it in range(dimtsim):
2018	ot = 0
2019	for ito in range(ot,dimtobs-1):
2020	if valdimobs['T'][ito] < simobstimes[it] and valdimobs['T'][ito+1] > \
2021	simobstimes[it]:
2022	ot = ito
2023	tdist = simobstimes[it] - valdimobs['T'][ito]
2024	coindtvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito], \
2025	tdist])
2026	Nsimt = Nsimt + 1
2027	if tsiminit == 0: tsiminit=simobstimes[it]
2028	tsimend = simobstimes[it]
2029	elif simobstimes[it] > valdimobs['T'][ito+1]:
2030	coindtvalues0st.append([Nsimt, ito, valdimobs['T'][ito],tsimend-tsiminit])
2031
2032	coindtvalues = np.array(coindtvalues0, dtype=np.float)
2033	coindtvaluesst = np.array(coindtvalues0st, dtype=np.float)
2034
2035	Ncoindt = len(coindtvalues[:,0])
2036	print main + ': found',Ncoindt,'Simulation time-interval (within consecutive ' + \
2037	'observed times) coincident times between simulation and observations'
2038
2039	if Ncoindt == 0:
2040	print warnmsg
2041	print main + ': no coincident times found !!'
2042	print ' stopping it'
2043	quit(-1)
2044
2045	# Validating
2046	##
2047
2048	onewnc = NetCDFFile(ofile, 'w')
2049
2050	# Dimensions
2051	for kst in range(Nstsim):
2052	newdim = onewnc.createDimension(prestdescsim[kst] + 'time',None)
2053	if stdescsim[kst] != 'E':
2054	newdim = onewnc.createDimension(prestdescsim[kst] + 'obstime',None)
2055
2056	newdim = onewnc.createDimension('bnds',2)
2057	newdim = onewnc.createDimension('couple',2)
2058	newdim = onewnc.createDimension('StrLength',StringLength)
2059	newdim = onewnc.createDimension('xaround',Ngrid*2+1)
2060	newdim = onewnc.createDimension('yaround',Ngrid*2+1)
2061	newdim = onewnc.createDimension('gstats',13)
2062	newdim = onewnc.createDimension('stats',5)
2063	newdim = onewnc.createDimension('tstats',6)
2064	newdim = onewnc.createDimension('Nstsim', 3)
2065
2066	# Variable dimensions
2067	##
2068	newvar = onewnc.createVariable('couple', 'c', ('couple','StrLength'))
2069	basicvardef(newvar, 'couple', 'couples of values', '-')
2070	writing_str_nc(newvar, ['sim','obs'], StringLength)
2071
2072	newvar = onewnc.createVariable('statistics', 'c', ('stats','StrLength'))
2073	basicvardef(newvar, 'statistics', 'statitics from values', '-')
2074	writing_str_nc(newvar, statsn, StringLength)
2075
2076	newvar = onewnc.createVariable('gstatistics', 'c', ('gstats','StrLength'))
2077	basicvardef(newvar, 'gstatistics', 'global statitics from values', '-')
2078	writing_str_nc(newvar, gstatsn, StringLength)
2079
2080	newvar = onewnc.createVariable('tstatistics', 'c', ('tstats','StrLength'))
2081	basicvardef(newvar, 'tstatistics', 'statitics from values along time', '-')
2082	writing_str_nc(newvar, ostatsn, StringLength)
2083
2084	newvar = onewnc.createVariable('ksimstatistics', 'c', ('Nstsim','StrLength'))
2085	basicvardef(newvar, 'ksimstatistics', 'kind of simulated statitics', '-')
2086	writing_str_nc(newvar, Lstdescsim, StringLength)
2087
2088	if obskind == 'trajectory':
2089	if dims.has_key('Z'):
2090	newdim = onewnc.createDimension('trj',3)
2091	else:
2092	newdim = onewnc.createDimension('trj',2)
2093
2094	newvar = onewnc.createVariable('obssimtrj','i',('obstime','trj'))
2095	basicvardef(newvar, 'obssimtrj', 'trajectory on the simulation grid', '-')
2096	newvar[:] = trajpos.transpose()
2097
2098	if dims.has_key('Z'):
2099	newdim = onewnc.createDimension('simtrj',4)
2100	else:
2101	newdim = onewnc.createDimension('simtrj',3)
2102
2103	Nvars = len(valvars)
2104	for ivar in range(Nvars):
2105	simobsvalues = []
2106
2107	varsimobs = valvars[ivar][0] + '_' + valvars[ivar][1]
2108	print ' ' + varsimobs + '... .. .'
2109
2110	if not oobs.variables.has_key(valvars[ivar][1]):
2111	print errormsg
2112	print ' ' + main + ": observations file has not '" + valvars[ivar][1] + \
2113	"' !!"
2114	quit(-1)
2115
2116	if not osim.variables.has_key(valvars[ivar][0]):
2117	if not searchInlist(varNOcheck, valvars[ivar][0]):
2118	print errormsg
2119	print ' ' + main + ": simulation file has not '" + valvars[ivar][0] + \
2120	"' !!"
2121	quit(-1)
2122	else:
2123	ovsim = compute_varNOcheck(osim, valvars[ivar][0])
2124	else:
2125	ovsim = osim.variables[valvars[ivar][0]]
2126
2127	for idn in ovsim.dimensions:
2128	if not searchInlist(simdims.values(),idn):
2129	print errormsg
2130	print ' ' + main + ": dimension '" + idn + "' of variable '" + \
2131	valvars[ivar][0] + "' not provided as reference coordinate [X,Y,Z,T] !!"
2132	quit(-1)
2133
2134	ovobs = oobs.variables[valvars[ivar][1]]
2135	if searchInlist(ovobs.ncattrs(),'_FillValue'):
2136	oFillValue = ovobs.getncattr('_FillValue')
2137	else:
2138	oFillValue = None
2139
2140	for kst in range(Nstsim):
2141	simobsvalues = []
2142
2143	timedn = prestdescsim[kst] + 'time'
2144	timeobsdn = prestdescsim[kst] + 'obstime'
2145	print ' ' + prestdescsim[kst] + ' ...'
2146
2147	if stdescsim[kst] == 'E':
2148	# Observed and simualted exact times
2149	simobsvalues, simobsSvalues, simobsTtvalues, trjsim = \
2150	getting_ValidationValues(obskind, Nexactt, dims, trajpos, ovsim, \
2151	ovobs, exacttvalues, oFillValue, Ngrid, 'instantaneous')
2152
2153	if ivar == 0:
2154	vname = prestdescsim[kst] + 'time'
2155	newvar = onewnc.createVariable(vname,'f8', (timedn))
2156	basicvardef(newvar, vname, 'exact coincident time observations and '+\
2157	'simulation', obstunits)
2158	set_attribute(newvar, 'calendar', 'standard')
2159	if Nexactt == 0:
2160	newvar[:] = np.float64(0.)
2161	else:
2162	newvar[:] = exacttvalues[:,3]
2163	if Nexactt == 0:
2164	simobsSvalues = np.zeros((1,2), dtype=np.float)
2165	simobsvalues = np.zeros((1,2), dtype=np.float)
2166
2167	dimt = Nexactt
2168
2169	elif stdescsim[kst] == 'C':
2170	# Simualted closest to Observed times
2171	simobsvalues, simobsSvalues, simobsTtvalues, trjsim = \
2172	getting_ValidationValues(obskind, Nclosest, dims, trajpos, ovsim, \
2173	ovobs, closesttvalues, oFillValue, Ngrid, 'instantaneous')
2174	dimt = Nclosest
2175
2176	if ivar == 0:
2177	vname = prestdescsim[kst] + 'time'
2178	newvar = onewnc.createVariable(vname,'f8', (timedn))
2179	basicvardef(newvar, vname, 'time simulations closest to observed ' + \
2180	'values', obstunits )
2181	set_attribute(newvar, 'calendar', 'standard')
2182	newvar[:] = closesttvalues[:,2]
2183
2184	vname = prestdescsim[kst] + 'obstime'
2185	newvar = onewnc.createVariable(vname,'f8', (vname))
2186	basicvardef(newvar, vname, 'closest time observations', obstunits)
2187	set_attribute(newvar, 'calendar', 'standard')
2188	newvar[:] = closesttvalues[:,3]
2189
2190	elif stdescsim[kst] == 'B':
2191	# Observed values temporally around coincident times
2192	simobsvalues, simobsSvalues, simobsTtvalues, trjsim = \
2193	getting_ValidationValues(obskind, Ncoindt, dims, trajpos, ovsim, ovobs,\
2194	coindtvalues, oFillValue, Ngrid, 'tbackwardSmean')
2195	dimt = simobsSvalues.shape[0]
2196
2197	if ivar == 0:
2198	vname = prestdescsim[kst] + 'time'
2199	newvar = onewnc.createVariable(vname,'f8', (timedn), \
2200	fill_value = fillValueF)
2201	basicvardef(newvar, vname, 'simulation time between observations', \
2202	obstunits)
2203	set_attribute(newvar, 'calendar', 'standard')
2204	set_attribute(newvar, 'bounds', 'time_bnds')
2205	newvar[:] = simobsTtvalues[:,1]
2206
2207	vname = prestdescsim[kst] + 'obstime'
2208	newvar = onewnc.createVariable(vname,'f8', (vname))
2209	basicvardef(newvar, vname, 'observed between time', obstunits )
2210	set_attribute(newvar, 'calendar', 'standard')
2211	newvar[:] = np.unique(coindtvalues[:,3])
2212
2213	# Re-arranging values...
2214	# For an incomprensible reason it is not working?
2215	# arrayvals = np.array(simobsvalues)
2216	arrayvals = np.zeros((len(simobsvalues),2), dtype=np.float)
2217	for it in range(len(simobsvalues)):
2218	arrayvals[it,:] = simobsvalues[it][0:2]
2219
2220	if len(valvars[ivar]) > 2:
2221	const=np.float(valvars[ivar][3])
2222	if valvars[ivar][2] == 'sumc':
2223	simobsSvalues = simobsSvalues + const
2224	arrayvals[:,0] = arrayvals[:,0] + const
2225	elif valvars[ivar][2] == 'subc':
2226	simobsSvalues = simobsSvalues - const
2227	arrayvals[:,0] = arrayvals[:,0] - const
2228	elif valvars[ivar][2] == 'mulc':
2229	simobsSvalues = simobsSvalues * const
2230	arrayvals[:,0] = arrayvals[:,0] * const
2231	elif valvars[ivar][2] == 'divc':
2232	simobsSvalues = simobsSvalues / const
2233	arrayvals[:,0] = arrayvals[:,0] / const
2234	else:
2235	print errormsg
2236	print ' ' + fname + ": operation '"+valvars[ivar][2]+"' not ready!!"
2237	quit(-1)
2238
2239	# for it in range(len(arrayvals[:,0])):
2240	# print it,arrayvals[it,:],':',simobsvalues[it]
2241	# quit()
2242
2243	if kst == 0:
2244	simstats = np.zeros((Nstsim,5), dtype=np.float)
2245	obsstats = np.zeros((Nstsim,5), dtype=np.float)
2246	simobsstats = np.zeros((Nstsim,13), dtype=np.float)
2247
2248	# statisics sim
2249	simstats[kst,0] = np.min(arrayvals[:,0])
2250	simstats[kst,1] = np.max(arrayvals[:,0])
2251	simstats[kst,2] = np.mean(arrayvals[:,0])
2252	simstats[kst,3] = np.mean(arrayvals[:,0]*arrayvals[:,0])
2253	simstats[kst,4] = np.sqrt(simstats[kst,3] - simstats[kst,2]*simstats[kst,2])
2254
2255	# statisics obs
2256	# Masking 'nan'
2257	obsmask0 = np.where(arrayvals[:,1] != arrayvals[:,1], fillValueF, \
2258	arrayvals[:,1])
2259
2260	obsmask = ma.masked_equal(obsmask0, fillValueF)
2261	obsmask2 = obsmask*obsmask
2262
2263	obsstats[kst,0] = obsmask.min()
2264	obsstats[kst,1] = obsmask.max()
2265	obsstats[kst,2] = obsmask.mean()
2266	obsstats[kst,3] = obsmask2.mean()
2267	obsstats[kst,4] = np.sqrt(obsstats[kst,3] - obsstats[kst,2]*obsstats[kst,2])
2268
2269	# Statistics sim-obs
2270	diffvals = np.zeros((dimt), dtype=np.float)
2271
2272	diffvals = arrayvals[:,0] - obsmask
2273
2274	diff2vals = diffvals*diffvals
2275	sumdiff = diffvals.sum()
2276	sumdiff2 = diff2vals.sum()
2277
2278	simobsstats[kst,0] = simstats[kst,0] - obsstats[kst,0]
2279	simobsstats[kst,1] = np.mean(arrayvals[:,0]*obsmask)
2280	simobsstats[kst,2] = diffvals.min()
2281	simobsstats[kst,3] = diffvals.max()
2282	simobsstats[kst,4] = diffvals.mean()
2283	simobsstats[kst,5] = np.abs(diffvals).mean()
2284	simobsstats[kst,6] = np.sqrt(diff2vals.mean())
2285	simobsstats[kst,7], simobsstats[kst,8] = sts.mstats.pearsonr(arrayvals[:,0], \
2286	obsmask)
2287	# From:
2288	#Willmott, C. J. 1981. 'On the validation of models. Physical Geography', 2, 184-194
2289	#Willmott, C. J. (1984). 'On the evaluation of model performance in physical
2290	# geography'. Spatial Statistics and Models, G. L. Gaile and C. J. Willmott, eds.,
2291	# 443-460
2292	#Willmott, C. J., S. G. Ackleson, R. E. Davis, J. J. Feddema, K. M. Klink, D. R.
2293	# Legates, J. O'Donnell, and C. M. Rowe (1985), 'Statistics for the Evaluation and
2294	# Comparison of Models', J. Geophys. Res., 90(C5), 8995-9005
2295	#Legates, D. R., and G. J. McCabe Jr. (1999), 'Evaluating the Use of "Goodness-of-Fit"
2296	# Measures in Hydrologic and Hydroclimatic Model Validation', Water Resour. Res.,
2297	# 35(1), 233-241
2298	#
2299	# Deviation of residuals (SDR)
2300	simobsstats[kst,9] = np.mean(np.sqrt(np.abs((diffvals-simobsstats[kst,0])* \
2301	(diffvals-obsstats[kst,0]))))
2302	# Index of Efficiency (IOE)
2303	obsbias2series = (obsmask - obsstats[kst,0])*(obsmask - obsstats[kst,0])
2304	sumobsbias2series = obsbias2series.sum()
2305
2306	simobsstats[kst,10] = 1. - sumdiff2/(sumobsbias2series)
2307	# Index of Agreement (IOA)
2308	simbias2series = arrayvals[:,0] - obsstats[kst,0]
2309	obsbias2series = obsmask - obsstats[kst,0]
2310
2311	abssimbias2series = np.abs(simbias2series)
2312	absobsbias2series = np.where(obsbias2series<0, -obsbias2series, \
2313	obsbias2series)
2314
2315	abssimobsbias2series = (abssimbias2series+absobsbias2series)*( \
2316	abssimbias2series + absobsbias2series)
2317
2318	simobsstats[kst,11] = 1. - sumdiff2/(abssimobsbias2series.sum())
2319	# Fractional Mean Bias (FMB)
2320	simobsstats[kst,12]=(simstats[kst,0]-obsstats[kst,0])/(0.5*(simstats[kst,0] +\
2321	obsstats[kst,0]))
2322
2323	# Statistics around sim values
2324	aroundstats = np.zeros((5,dimt), dtype=np.float)
2325	for it in range(dimt):
2326	aroundstats[0,it] = np.min(simobsSvalues[it,])
2327	aroundstats[1,it] = np.max(simobsSvalues[it,])
2328	aroundstats[2,it] = np.mean(simobsSvalues[it,])
2329	aroundstats[3,it] = np.mean(simobsSvalues[it,]*simobsSvalues[it,])
2330	aroundstats[4,it] = np.sqrt(aroundstats[3,it] - aroundstats[2,it]* \
2331	aroundstats[2,it])
2332
2333	# sim Values to netCDF
2334	newvar = onewnc.createVariable(valvars[ivar][0] + '_' + stdescsim[kst] + \
2335	'sim', 'f', (timedn), fill_value=fillValueF)
2336	descvar = prestdescsim[kst] + ' time simulated: ' + valvars[ivar][0]
2337	basicvardef(newvar, valvars[ivar][0], descvar, ovobs.getncattr('units'))
2338	if stdescsim[kst] == 'E' and Nexactt == 0:
2339	newvar[:] = fillValueF
2340	else:
2341	newvar[:] = arrayvals[:,0]
2342
2343	# obs Values to netCDF
2344	if stdescsim[kst] != 'E':
2345	newvar = onewnc.createVariable(valvars[ivar][1] + '_' + stdescsim[kst] + \
2346	'obs', 'f', (timeobsdn), fill_value=fillValueF)
2347	else:
2348	newvar = onewnc.createVariable(valvars[ivar][1] + '_' + stdescsim[kst] + \
2349	'obs', 'f', (timedn), fill_value=fillValueF)
2350
2351	descvar = prestdescsim[kst] + ' time observed: ' + valvars[ivar][1]
2352	basicvardef(newvar, valvars[ivar][1], descvar, ovobs.getncattr('units'))
2353
2354	if stdescsim[kst] == 'E' and Nexactt == 0:
2355	newvar[:] = fillValueF
2356	else:
2357	newvar[:] = arrayvals[:,1]
2358
2359	# Including statistics of the between simulated values
2360	if stdescsim[kst] == 'B':
2361	if not searchInlist(onewnc.dimensions, 'Bstats'):
2362	newdim = onewnc.createDimension('Bstats',3)
2363
2364	newvar = onewnc.createVariable('Bstats','c',('Bstats','StrLength'))
2365	descvar = prestdescsim[kst] + ' time simulated: ' + valvars[ivar][0]
2366	basicvardef(newvar, 'Bstats', 'Between values statistics', \
2367	ovobs.getncattr('units'))
2368	Bstsvals=['minimum','maximum','stdandard deviation']
2369	writing_str_nc(newvar,Bstsvals,StringLength)
2370
2371	newvar = onewnc.createVariable(valvars[ivar][0] + '_Bstats_sim', 'f', \
2372	(timedn,'Bstats'), fill_value=fillValueF)
2373	descvar = prestdescsim[kst] + ' time simulated: ' + valvars[ivar][0] + \
2374	' statistics'
2375	basicvardef(newvar, valvars[ivar][0], descvar, ovobs.getncattr('units'))
2376	newvar[0,:] = [fillValueF, fillValueF, fillValueF]
2377	for it in range(dimt-1):
2378	newvar[it+1,:] = simobsvalues[it][2:5]
2379
2380	# Around values
2381	if not onewnc.variables.has_key(valvars[ivar][0] + 'around'):
2382	vname = prestdescsim[kst] + valvars[ivar][0] + 'around'
2383	else:
2384	vname = prestdescsim[kst] + valvars[ivar][0] + 'Around'
2385
2386	if dims.has_key('Z'):
2387	if not onewnc.dimensions.has_key('zaround'):
2388	newdim = onewnc.createDimension('zaround',Ngrid*2+1)
2389	newvar = onewnc.createVariable(vname, 'f', (timedn,'zaround', \
2390	'yaround','xaround'), fill_value=fillValueF)
2391	else:
2392	newvar = onewnc.createVariable(vname, 'f', (timedn,'yaround','xaround'), \
2393	fill_value=fillValueF)
2394
2395	descvar = prestdescsim[kst] + 'around simulated values +/- grid values: ' + \
2396	valvars[ivar][0]
2397	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2398
2399	if stdescsim[kst] == 'E' and Nexactt == 0:
2400	newvar[:] = np.ones((0,Ngrid2+1,Ngrid2+1))*fillValueF
2401	else:
2402	newvar[:] = simobsSvalues
2403
2404
2405	# around sim Statistics
2406	if not searchInlist(onewnc.variables,prestdescsim[kst] + valvars[ivar][0] + \
2407	'staround'):
2408	vname = prestdescsim[kst] + valvars[ivar][0] + 'staround'
2409	else:
2410	vname = prestdescsim[kst] + valvars[ivar][0] + 'Staround'
2411
2412	newvar = onewnc.createVariable(vname, 'f', (timedn,'stats'), \
2413	fill_value=fillValueF)
2414	descvar = prestdescsim[kst] + ' around (' + str(Ngrid) + ', ' + str(Ngrid) +\
2415	') simulated statisitcs: ' + valvars[ivar][0]
2416	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2417	set_attribute(newvar, 'cell_methods', 'time_bnds')
2418	if stdescsim[kst] == 'E' and Nexactt == 0:
2419	newvar[:] = np.ones((0,5))*fillValueF
2420	else:
2421	newvar[:] = aroundstats.transpose()
2422
2423	if stdescsim[kst] == 'B':
2424	if not searchInlist(onewnc.variables, 'time_bnds'):
2425	newvar = onewnc.createVariable('time_bnds','f8',(timedn,'bnds'))
2426	basicvardef(newvar, 'time_bnds', timedn, obstunits )
2427	set_attribute(newvar, 'calendar', 'standard')
2428	newvar[:] = simobsTtvalues
2429
2430	# sim Statistics
2431	if not searchInlist(onewnc.variables,valvars[ivar][0] + 'stsim'):
2432	vname = valvars[ivar][0] + 'stsim'
2433	else:
2434	vname = valvars[ivar][0] + 'stSim'
2435
2436	newvar = onewnc.createVariable(vname, 'f', ('Nstsim', 'stats'), \
2437	fill_value=fillValueF)
2438	descvar = 'simulated statisitcs: ' + valvars[ivar][0]
2439	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2440	newvar[:] = simstats
2441
2442	# obs Statistics
2443	if not searchInlist(onewnc.variables,valvars[ivar][1] + 'stobs'):
2444	vname = valvars[ivar][1] + 'stobs'
2445	else:
2446	vname = valvars[ivar][1] + 'stObs'
2447
2448	newvar = onewnc.createVariable(vname, 'f', ('Nstsim', 'stats'), \
2449	fill_value=fillValueF)
2450	descvar = 'observed statisitcs: ' + valvars[ivar][1]
2451	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2452	newvar[:] = obsstats
2453
2454	# sim-obs Statistics
2455	if not searchInlist(onewnc.variables,varsimobs + 'st'):
2456	vname = varsimobs + 'st'
2457	else:
2458	vname = varSimobs + 'st'
2459
2460	newvar = onewnc.createVariable(vname, 'f', ('Nstsim', 'gstats'), \
2461	fill_value=fillValueF)
2462	descvar = 'simulated-observed tatisitcs: ' + varsimobs
2463	basicvardef(newvar, vname, descvar, ovobs.getncattr('units'))
2464	newvar[:] = simobsstats
2465
2466	onewnc.sync()
2467
2468	if trjsim is not None:
2469	newvar = onewnc.createVariable('simtrj','i',('betweentime','simtrj'))
2470	basicvardef(newvar,'simtrj','coordinates [X,Y,Z,T] of the coincident ' + \
2471	'trajectory in sim', obstunits)
2472	newvar[:] = trjsim.transpose()
2473
2474	# Adding three variables with the station location, longitude, latitude and height
2475	if obskind == 'single-station':
2476	adding_station_desc(onewnc,stationdesc)
2477
2478	# Global attributes
2479	##
2480	set_attribute(onewnc,'author_nc','Lluis Fita')
2481	set_attribute(onewnc,'institution_nc','Laboratoire de Meteorology Dynamique, ' + \
2482	'LMD-Jussieu, UPMC, Paris')
2483	set_attribute(onewnc,'country_nc','France')
2484	set_attribute(onewnc,'script_nc',main)
2485	set_attribute(onewnc,'version_script',version)
2486	set_attribute(onewnc,'information', \
2487	'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html')
2488	set_attribute(onewnc,'simfile',opts.fsim)
2489	set_attribute(onewnc,'obsfile',opts.fobs)
2490
2491	onewnc.sync()
2492	onewnc.close()
2493
2494	print main + ": successfull writting of '" + ofile + "' !!"

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

Context Navigation

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

Download in other formats: