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

Last change on this file since 485 was 468, checked in by lfita, 10 years ago

'Final'? correction for the `None' vardim for the simulation
Solved issue with the masked values from the observations
File size: 67.1 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
8import numpy as np
9import os
10import re
11from optparse import OptionParser
12from netCDF4 import Dataset as NetCDFFile
13from scipy import stats as sts
14import numpy.ma as ma
15
16main = 'validarion_sim.py'
17errormsg = 'ERROR -- errror -- ERROR -- error'
18warnmsg = 'WARNING -- warning -- WARNING -- warning'
19
20# version
21version=1.0
22
23# Filling values for floats, integer and string
24fillValueF = 1.e20
25fillValueI = -99999
26fillValueS = '---'
27
28StringLength = 50
29
30# Number of grid points to take as 'environment' around the observed point
31Ngrid = 1
32
33def 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
45def 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
62def 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
75def 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
96def 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
167def 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
231def 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
269def 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
299def 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
331def 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-yM*1000-yC*100)/10
491        yU = yr-yM*1000-yC*100-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
534def 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 = tvalB*60.*10.e6
564            elif tuB == 'hours':
565                tB = tvalB*3600.*10.e6
566            elif tuB == 'days':
567                tB = tvalB*3600.*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 = tvalB*3600.*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 = tvalB*60.*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.days*24.*3600.*10.e6 + difft.seconds*10.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
670def 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
723def 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        'WRFz': height from WRF variables
733    """
734    fname = 'compute_varNOcheck'
735
736    if varn == 'WRFdens':
737#        print '    ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \
738#          'DNW)/g ...'
739        grav = 9.81
740
741# Just we need in in absolute values: Size of the central grid cell
742##    dxval = ncobj.getncattr('DX')
743##    dyval = ncobj.getncattr('DY')
744##    mapfac = ncobj.variables['MAPFAC_M'][:]
745##    area = dxval*dyval*mapfac
746        dimensions = ncobj.variables['MU'].dimensions
747
748        mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
749        dnw = ncobj.variables['DNW'][:]
750
751        varNOcheckv = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \
752          dtype=np.float)
753        levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
754
755        for it in range(mu.shape[0]):
756            for iz in range(dnw.shape[1]):
757                levval.fill(np.abs(dnw[it,iz]))
758                varNOcheck[it,iz,:,:] = levval
759                varNOcheck[it,iz,:,:] = mu[it,:,:]*varNOcheck[it,iz,:,:]/grav
760
761    elif varn == 'WRFght':
762#        print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...'
763        varNOcheckv = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
764        dimensions = ncobj.variables['PH'].dimensions
765
766    elif varn == 'WRFp':
767#        print '  ' + fname + ': Retrieving pressure value from WRF as P + PB'
768        varNOcheckv = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
769        dimensions = ncobj.variables['P'].dimensions
770
771    elif varn == 'WRFrh':
772#        print '    ' + main + ": computing relative humidity from WRF as 'Tetens'" +\
773#         ' equation (T,P) ...'
774        p0=100000.
775        p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
776        tk = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
777        qv = ncobj.variables['QVAPOR'][:]
778
779        data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
780        data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
781
782        varNOcheckv = qv/data2
783        dimensions = ncobj.variables['P'].dimensions
784
785    elif varn == 'WRFt':
786#        print '    ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
787        p0=100000.
788        p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
789
790        varNOcheckv = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
791        dimensions = ncobj.variables['T'].dimensions
792
793    elif varn == 'WRFz':
794        grav = 9.81
795#        print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...'
796        varNOcheckv = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/grav
797        dimensions = ncobj.variables['PH'].dimensions
798
799    else:
800        print erromsg
801        print '  ' + fname + ": variable '" + varn + "' nor ready !!"
802        quit(-1)
803
804    return varNOcheck
805
806class compute_varNOcheck(object):
807    """ Class to compute variables which are not originary in the file
808      ncobj= netCDF object file
809      varn = variable to compute:
810        'WRFdens': air density from WRF variables
811        'WRFght': geopotential height from WRF variables
812        'WRFp': pressure from WRF variables
813        'WRFrh': relative humidty fom WRF variables
814        'WRFT': CF-time from WRF variables
815        'WRFt': temperature from WRF variables
816        'WRFtd': dew-point temperature from WRF variables
817        'WRFws': wind speed from WRF variables
818        'WRFwd': wind direction from WRF variables
819        'WRFz': height from WRF variables
820    """
821    fname = 'compute_varNOcheck'
822
823    def __init__(self, ncobj, varn):
824
825        if ncobj is None:
826            self = None
827            self.dimensions = None
828            self.shape = None
829            self.__values = None
830        else:
831            if varn == 'WRFdens':
832#        print '    ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \
833#          'DNW)/g ...'
834                grav = 9.81
835
836# Just we need in in absolute values: Size of the central grid cell
837##    dxval = ncobj.getncattr('DX')
838##    dyval = ncobj.getncattr('DY')
839##    mapfac = ncobj.variables['MAPFAC_M'][:]
840##    area = dxval*dyval*mapfac
841                dimensions = ncobj.variables['MU'].dimensions
842                shape = ncobj.variables['MU'].shape
843
844                mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
845                dnw = ncobj.variables['DNW'][:]
846
847                varNOcheckv = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \
848                  dtype=np.float)
849                levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
850
851                for it in range(mu.shape[0]):
852                    for iz in range(dnw.shape[1]):
853                        levval.fill(np.abs(dnw[it,iz]))
854                        varNOcheck[it,iz,:,:] = levval
855                        varNOcheck[it,iz,:,:] = mu[it,:,:]*varNOcheck[it,iz,:,:]/grav
856
857            elif varn == 'WRFght':
858#        print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...'
859                varNOcheckv = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
860                dimensions = ncobj.variables['PH'].dimensions
861                shape = ncobj.variables['PH'].shape
862
863            elif varn == 'WRFp':
864#        print '  ' + fname + ': Retrieving pressure value from WRF as P + PB'
865                varNOcheckv = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
866                dimensions = ncobj.variables['P'].dimensions
867                shape = ncobj.variables['P'].shape
868
869            elif varn == 'WRFrh':
870#        print '    ' + main + ": computing relative humidity from WRF as 'Tetens'" +\
871#         ' equation (T,P) ...'
872                p0=100000.
873                p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
874                tk = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
875                qv = ncobj.variables['QVAPOR'][:]
876
877                data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
878                data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
879
880                varNOcheckv = qv/data2
881                dimensions = ncobj.variables['P'].dimensions
882                shape = ncobj.variables['P'].shape
883
884            elif varn == 'WRFT':
885# To compute CF-times from WRF kind
886#
887                import datetime as dt
888
889                times = ncobj.variables['Times']
890                dimt = times.shape[0]
891                varNOcheckv = np.zeros((dimt), dtype=np.float64)
892                self.unitsval = 'seconds since 1949-12-01 00:00:00'
893                refdate = datetimeStr_datetime('1949-12-01_00:00:00')
894
895                dimensions = tuple([ncobj.variables['Times'].dimensions[0]])
896                shape = tuple([dimt])
897
898                for it in range(dimt):
899                    datevalS = datetimeStr_conversion(times[it,:], 'WRFdatetime',    \
900                      'YmdHMS')
901                    dateval = dt.datetime.strptime(datevalS, '%Y%m%d%H%M%S')
902                    difft = dateval - refdate
903                    varNOcheckv[it] = difft.days*3600.*24. + difft.seconds +        \
904                          np.float(int(difft.microseconds/10.e6))
905
906            elif varn == 'WRFt':
907#        print '    ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
908                p0=100000.
909                p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
910
911                varNOcheckv = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
912                dimensions = ncobj.variables['T'].dimensions
913                shape = ncobj.variables['P'].shape
914
915            elif varn == 'WRFtd':
916#        print '    ' + main + ': computing dew-point temperature from WRF as inv_potT(T + 300) and Tetens...'
917# tacking from: http://en.wikipedia.org/wiki/Dew_point
918                p0=100000.
919                p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
920
921                temp = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
922
923                qv = ncobj.variables['QVAPOR'][:]
924
925                tk = temp - 273.15
926                data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
927                data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
928
929                rh = qv/data2
930               
931                pa = rh * data1/100.
932                varNOcheckv = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
933
934                dimensions = ncobj.variables['T'].dimensions
935                shape = ncobj.variables['P'].shape
936
937            elif varn == 'WRFws':
938#        print '    ' + main + ': computing wind speed from WRF as SQRT(U**2 + V**2) ...'
939                uwind = ncobj.variables['U'][:]
940                vwind = ncobj.variables['V'][:]
941                dx = uwind.shape[3]
942                dy = vwind.shape[2]
943                 
944# de-staggering
945                ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
946                va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])
947
948                varNOcheckv = np.sqrt(ua*ua + va*va)
949                dimensions = tuple(['Time','bottom_top','south_north','west_east'])
950                shape = ua.shape
951
952            elif varn == 'WRFwd':
953#        print '    ' + main + ': computing wind direction from WRF as ATAN2PI(V,U) ...'
954                uwind = ncobj.variables['U'][:]
955                vwind = ncobj.variables['V'][:]
956                dx = uwind.shape[3]
957                dy = vwind.shape[2]
958
959# de-staggering
960                ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
961                va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])
962
963                theta = np.arctan2(ua, va)
964                dimensions = tuple(['Time','bottom_top','south_north','west_east'])
965                shape = ua.shape
966                varNOcheckv = 360.*(1. + theta/(2.*np.pi))
967
968            elif varn == 'WRFz':
969                grav = 9.81
970#        print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...'
971                varNOcheckv = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/grav
972                dimensions = ncobj.variables['PH'].dimensions
973                shape = ncobj.variables['PH'].shape
974
975            else:
976                print errormsg
977                print '  ' + fname + ": variable '" + varn + "' nor ready !!"
978                quit(-1)
979
980            self.dimensions = dimensions
981            self.shape = shape
982            self.__values = varNOcheckv
983
984    def __getitem__(self,elem):
985        return self.__values[elem]
986
987####### ###### ##### #### ### ## #
988
989strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " +  \
990  " and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')"
991
992kindobs=['multi-points', 'single-station', 'trajectory']
993strkObs="kind of observations; 'multi-points': multiple individual punctual obs " +  \
994  "(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
995  "'trajectory': following a trajectory"
996simopers = ['sumc','subc','mulc','divc']
997opersinf = 'sumc,[constant]: add [constant] to variables values; subc,[constant]: '+ \
998  'substract [constant] to variables values; mulc,[constant]: multipy by ' +         \
999  '[constant] to variables values; divc,[constant]: divide by [constant] to ' +      \
1000  'variables values'
1001varNOcheck = ['WRFdens', 'WRFght', 'WRFp', 'WRFrh', 'WRFT', 'WRFt', 'WRFtd', 'WRFws',\
1002  'WRFwd', 'WRFz']
1003varNOcheckinf = "'WRFdens': air density from WRF variables; 'WRFght': geopotential"+ \
1004  " height from WRF variables; 'WRFp': pressure from WRF variables; 'WRFrh': " +     \
1005  "relative humidty fom WRF variables; 'WRFT': CF-time from WRF variables'WRFt'; " + \
1006  " temperature from WRF variables; 'WRFtd': dew-point temperature from WRF " +      \
1007  "variables; 'WRFws': wind speed from WRF variables; 'WRFwd': wind speed " +        \
1008  "direction from WRF variables; 'WRFz': height from WRF variables"
1009
1010dimshelp = "[DIM]@[simdim]@[obsdim] ',' list of couples of dimensions names from " + \
1011  "each source ([DIM]='X','Y','Z','T'; None, no value)"
1012vardimshelp = "[DIM]@[simvardim]@[obsvardim] ',' list of couples of variables " +    \
1013  "names with dimensions values from each source ([DIM]='X','Y','Z','T'; None, " +   \
1014  "no value, WRFdiagnosted variables also available: " + varNOcheckinf + ")"
1015varshelp="[simvar]@[obsvar]@[[oper]@[val]] ',' list of couples of variables to " +   \
1016  "validate and if necessary operation and value operations: " + opersinf +          \
1017  "(WRFdiagnosted variables also available: " + varNOcheckinf + ")"
1018statsn = ['minimum', 'maximum', 'mean', 'mean2', 'standard deviation']
1019gstatsn = ['bias', 'simobs_mean', 'sim_obsmin', 'sim_obsmax', 'sim_obsmean', 'mae',  \
1020  'rmse', 'r_pearsoncorr', 'p_pearsoncorr']
1021ostatsn = ['number of points', 'minimum', 'maximum', 'mean', 'mean2',                \
1022  'standard deviation']
1023
1024parser = OptionParser()
1025parser.add_option("-d", "--dimensions", dest="dims", help=dimshelp, metavar="VALUES")
1026parser.add_option("-D", "--vardimensions", dest="vardims",
1027  help=vardimshelp, metavar="VALUES")
1028parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs, 
1029  help=strkObs, metavar="FILE")
1030parser.add_option("-l", "--stationLocation", dest="stloc", 
1031  help="longitude, latitude and height of the station (only for 'single-station')", 
1032  metavar="FILE")
1033parser.add_option("-o", "--observation", dest="fobs",
1034  help="observations file to validate", metavar="FILE")
1035parser.add_option("-s", "--simulation", dest="fsim",
1036  help="simulation file to validate", metavar="FILE")
1037parser.add_option("-t", "--trajectoryfile", dest="trajf",
1038  help="file with grid points of the trajectory in the simulation grid ('simtrj')", 
1039  metavar="FILE")
1040parser.add_option("-v", "--variables", dest="vars",
1041  help=varshelp, metavar="VALUES")
1042
1043(opts, args) = parser.parse_args()
1044
1045#######    #######
1046## MAIN
1047    #######
1048
1049ofile='validation_sim.nc'
1050
1051if opts.dims is None:
1052    print errormsg
1053    print '  ' + main + ': No list of dimensions are provided!!'
1054    print '    a ',' list of values X@[dimxsim]@[dimxobs],...,T@[dimtsim]@[dimtobs]'+\
1055      ' is needed'
1056    quit(-1)
1057else:
1058    simdims = {}
1059    obsdims = {}
1060    print main +': couple of dimensions _______'
1061    dims = {}
1062    ds = opts.dims.split(',')
1063    for d in ds:
1064        dsecs = d.split('@')
1065        if len(dsecs) != 3:
1066            print errormsg
1067            print '  ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
1068            print '    [DIM]@[dimnsim]@[dimnobs]'
1069            quit(-1)
1070        dims[dsecs[0]] = [dsecs[1], dsecs[2]]
1071        simdims[dsecs[0]] = dsecs[1]
1072        obsdims[dsecs[0]] = dsecs[2]
1073
1074        print '  ',dsecs[0],':',dsecs[1],',',dsecs[2]
1075       
1076if opts.vardims is None:
1077    print errormsg
1078    print '  ' + main + ': No list of variables with dimension values are provided!!'
1079    print '    a ',' list of values X@[vardimxsim]@[vardimxobs],...,T@' +  \
1080      '[vardimtsim]@[vardimtobs] is needed'
1081    quit(-1)
1082else:
1083    print main +': couple of variable dimensions _______'
1084    vardims = {}
1085    ds = opts.vardims.split(',')
1086    for d in ds:
1087        dsecs = d.split('@')
1088        if len(dsecs) != 3:
1089            print errormsg
1090            print '  ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
1091            print '    [DIM]@[vardimnsim]@[vardimnobs]'
1092            quit(-1)
1093        vardims[dsecs[0]] = [dsecs[1], dsecs[2]]
1094        print '  ',dsecs[0],':',dsecs[1],',',dsecs[2]
1095
1096if opts.obskind is None:
1097    print errormsg
1098    print '  ' + main + ': No kind of observations provided !!'
1099    quit(-1)
1100else:
1101    obskind = opts.obskind
1102    if obskind == 'single-station':
1103        if opts.stloc is None:
1104            print errormsg
1105            print '  ' + main + ': No station location provided !!'
1106            quit(-1)
1107        else:
1108            stationdesc = [np.float(opts.stloc.split(',')[0]),                       \
1109              np.float(opts.stloc.split(',')[1]), np.float(opts.stloc.split(',')[2])]
1110
1111if opts.fobs is None:
1112    print errormsg
1113    print '  ' + main + ': No observations file is provided!!'
1114    quit(-1)
1115else:
1116    if not os.path.isfile(opts.fobs):
1117        print errormsg
1118        print '   ' + main + ": observations file '" + opts.fobs + "' does not exist !!"
1119        quit(-1)
1120
1121if opts.fsim is None:
1122    print errormsg
1123    print '  ' + main + ': No simulation file is provided!!'
1124    quit(-1)
1125else:
1126    if not os.path.isfile(opts.fsim):
1127        print errormsg
1128        print '   ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
1129        quit(-1)
1130
1131if opts.vars is None:
1132    print errormsg
1133    print '  ' + main + ': No list of couples of variables is provided!!'
1134    print '    a ',' list of values [varsim]@[varobs],... is needed'
1135    quit(-1)
1136else:
1137    valvars = []
1138    vs = opts.vars.split(',')
1139    for v in vs:
1140        vsecs = v.split('@')
1141        if len(vsecs) < 2:
1142            print errormsg
1143            print '  ' + main + ': wrong number of values in:',vsecs,                \
1144              ' at least 2 are needed !!'
1145            print '    [varsim]@[varobs]@[[oper][val]]'
1146            quit(-1)
1147        if len(vsecs) > 2:
1148            if not searchInlist(simopers,vsecs[2]): 
1149                print errormsg
1150                print main + ": operation on simulation values '" + vsecs[2] +       \
1151                  "' not ready !!"
1152                quit(-1)
1153
1154        valvars.append(vsecs)
1155
1156# Openning observations trajectory
1157##
1158oobs = NetCDFFile(opts.fobs, 'r')
1159
1160valdimobs = {}
1161for dn in dims:
1162    print dn,':',dims[dn]
1163    if dims[dn][1] != 'None':
1164        if not oobs.dimensions.has_key(dims[dn][1]):
1165            print errormsg
1166            print '  ' + main + ": observations file does not have dimension '" +    \
1167              dims[dn][1] + "' !!"
1168            quit(-1)
1169        if vardims[dn][1] != 'None':
1170            if not oobs.variables.has_key(vardims[dn][1]):
1171                print errormsg
1172                print '  ' + main + ": observations file does not have varibale " +  \
1173                  "dimension '" + vardims[dn][1] + "' !!"
1174                quit(-1)
1175            valdimobs[dn] = oobs.variables[vardims[dn][1]][:]
1176    else:
1177        if dn == 'X':
1178            valdimobs[dn] = stationdesc[0]
1179        elif dn == 'Y':
1180            valdimobs[dn] = stationdesc[1]
1181        elif dn == 'Z':
1182            valdimobs[dn] = stationdesc[2]
1183
1184osim = NetCDFFile(opts.fsim, 'r')
1185
1186valdimsim = {}
1187for dn in dims:
1188    if dims[dn][0] != 'None':
1189        if not osim.dimensions.has_key(dims[dn][0]):
1190            print errormsg
1191            print '  ' + main + ": simulation file '" + opts.fsim +                  \
1192              "' does not have dimension '" + dims[dn][0] + "' !!"
1193            print '    it has: ',osim.dimensions
1194            quit(-1)
1195
1196        if not osim.variables.has_key(vardims[dn][0]) and                            \
1197          not searchInlist(varNOcheck,vardims[dn][0]):
1198            print errormsg
1199            print '  ' + main + ": simulation file '" + opts.fsim +                  \
1200              "' does not have varibale dimension '" + vardims[dn][0] + "' !!"
1201            print '    it has variables:',osim.variables
1202            quit(-1)
1203        if searchInlist(varNOcheck,vardims[dn][0]):
1204            valdimsim[dn] = compute_varNOcheck(osim, vardims[dn][0])
1205        else:
1206            valdimsim[dn] = osim.variables[vardims[dn][0]][:]
1207
1208# General characteristics
1209dimtobs = valdimobs['T'].shape[0]
1210dimtsim = valdimsim['T'].shape[0]
1211
1212print main +': observational time-steps:',dimtobs,'simulation:',dimtsim
1213
1214notfound = np.zeros((dimtobs), dtype=int)
1215
1216if obskind == 'multi-points':
1217    trajpos = np.zeros((2,dimt),dtype=int)
1218    for it in range(dimtobs):
1219        trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'],                    \
1220          [valdimobs['X'][it],valdimobss['Y'][it]])
1221elif obskind == 'single-station':
1222    stationpos = np.zeros((2), dtype=int)
1223    if valdimsim.has_key('X') and valdimsim.has_key('Y'):
1224        stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],[valdimobs['Y'],         \
1225          valdimobs['X']])
1226        iid = 0
1227        for idn in osim.variables[vardims['X'][0]].dimensions:
1228            if idn == dims['X'][0]:
1229                stationpos[1] = stsimpos[iid]
1230            elif idn == dims['Y'][0]:
1231                stationpos[0] = stsimpos[iid]
1232
1233            iid = iid + 1
1234        print main + ': station point in simulation:', stationpos
1235        print '    station position:',valdimobs['X'],',',valdimobs['Y']
1236        print '    simulation coord.:',valdimsim['X'][tuple(stsimpos)],',',          \
1237          valdimsim['Y'][tuple(stsimpos)]
1238    else:
1239        print main + ': validation with two time-series !!'
1240
1241elif obskind == 'trajectory':
1242    if opts.trajf is not None:
1243        if not os.path.isfile(opts.fsim):
1244            print errormsg
1245            print '   ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
1246            quit(-1)
1247        else:
1248            otrjf = NetCDFFile(opts.fsim, 'r')
1249            trajpos[0,:] = otrjf.variables['obssimtrj'][0]
1250            trajpos[1,:] = otrjf.variables['obssimtrj'][1]
1251            otrjf.close()
1252    else:
1253        if dims.has_key('Z'):
1254            trajpos = np.zeros((3,dimtobs),dtype=int)
1255            for it in range(dimtobs):
1256                if np.mod(it*100./dimtobs,10.) == 0.:
1257                    print '    trajectory done: ',it*100./dimtobs,'%'
1258                stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],                 \
1259                  [valdimobs['Y'][it],valdimobs['X'][it]])
1260                stationpos = np.zeros((2), dtype=int)
1261                iid = 0
1262                for idn in osim.variables[vardims['X'][0]].dimensions:
1263                    if idn == dims['X'][0]:
1264                        stationpos[1] = stsimpos[iid]
1265                    elif idn == dims['Y'][0]:
1266                        stationpos[0] = stsimpos[iid]
1267                    iid = iid + 1
1268                if stationpos[0] == 0 and stationpos[1] == 0: notfound[it] = 1
1269             
1270                trajpos[0,it] = stationpos[0]
1271                trajpos[1,it] = stationpos[1]
1272# In the simulation 'Z' varies with time ... non-hydrostatic model! ;)
1273#                trajpos[2,it] = index_mat(valdimsim['Z'][it,:,stationpos[0],         \
1274#                  stationpos[1]], valdimobs['Z'][it])
1275        else:
1276            trajpos = np.zeros((2,dimtobs),dtype=int)
1277            for it in range(dimtobs):
1278                stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],                 \
1279                  [valdimobs['Y'][it],valdimobss['X'][it]])
1280                stationpos = np.zeros((2), dtype=int)
1281                iid = 0
1282                for idn in osim.variables[vardims['X'][0]].dimensions:
1283                    if idn == dims['X'][0]:
1284                        stationpos[1] = stsimpos[iid]
1285                    elif idn == dims['Y'][0]:
1286                        stationpos[0] = stsimpos[iid]
1287                    iid = iid + 1
1288                if stationpos[0] == 0 or stationpos[1] == 0: notfound[it] = 1
1289
1290                trajpos[0,it] = stationspos[0]
1291                trajpos[1,it] = stationspos[1]
1292
1293        print main + ': not found',np.sum(notfound),'points of the trajectory'
1294
1295# Getting times
1296tobj = oobs.variables[vardims['T'][1]]
1297obstunits = tobj.getncattr('units')
1298if vardims['T'][0] == 'WRFT':
1299    tsim = valdimsim['T'][:]
1300    simtunits = 'seconds since 1949-12-01 00:00:00'
1301else:
1302    tsim = osim.variables[vardims['T'][0]][:]
1303    otsim = osim.variables[vardims['T'][0]]
1304    simtunits = otsim.getncattr('units')
1305
1306simobstimes = coincident_CFtimes(tsim, obstunits, simtunits)
1307
1308# Concident times
1309##
1310coindtvalues0 = []
1311for it in range(dimtsim):   
1312    ot = 0
1313    for ito in range(ot,dimtobs-1):
1314        if valdimobs['T'][ito] < simobstimes[it] and valdimobs['T'][ito+1] >         \
1315          simobstimes[it]:
1316            ot = ito
1317            tdist = simobstimes[it] - valdimobs['T'][ito]
1318            coindtvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito],      \
1319              tdist])
1320
1321coindtvalues = np.array(coindtvalues0, dtype=np.float)
1322
1323Ncoindt = len(coindtvalues[:,0])
1324print main + ': found',Ncoindt,'coincident times between simulation and observations'
1325
1326if Ncoindt == 0:
1327    print warnmsg
1328    print main + ': no coincident times found !!'
1329    print '  stopping it'
1330    quit(-1)
1331
1332# Validating
1333##
1334
1335onewnc = NetCDFFile(ofile, 'w')
1336
1337# Dimensions
1338newdim = onewnc.createDimension('time',None)
1339newdim = onewnc.createDimension('bnds',2)
1340newdim = onewnc.createDimension('obstime',None)
1341newdim = onewnc.createDimension('couple',2)
1342newdim = onewnc.createDimension('StrLength',StringLength)
1343newdim = onewnc.createDimension('xaround',Ngrid*2+1)
1344newdim = onewnc.createDimension('yaround',Ngrid*2+1)
1345newdim = onewnc.createDimension('gstats',9)
1346newdim = onewnc.createDimension('stats',5)
1347newdim = onewnc.createDimension('tstats',6)
1348
1349# Variable dimensions
1350##
1351newvar = onewnc.createVariable('obstime','f8',('time'))
1352basicvardef(newvar, 'obstime', 'time observations', obstunits )
1353set_attribute(newvar, 'calendar', 'standard')
1354set_attribute(newvar, 'bounds', 'time_bnds')
1355newvar[:] = coindtvalues[:,3]
1356
1357newvar = onewnc.createVariable('couple', 'c', ('couple','StrLength'))
1358basicvardef(newvar, 'couple', 'couples of values', '-')
1359writing_str_nc(newvar, ['sim','obs'], StringLength)
1360
1361newvar = onewnc.createVariable('statistics', 'c', ('stats','StrLength'))
1362basicvardef(newvar, 'statistics', 'statitics from values', '-')
1363writing_str_nc(newvar, statsn, StringLength)
1364
1365newvar = onewnc.createVariable('gstatistics', 'c', ('gstats','StrLength'))
1366basicvardef(newvar, 'gstatistics', 'global statitics from values', '-')
1367writing_str_nc(newvar, gstatsn, StringLength)
1368
1369newvar = onewnc.createVariable('tstatistics', 'c', ('tstats','StrLength'))
1370basicvardef(newvar, 'tstatistics', 'statitics from values along time', '-')
1371writing_str_nc(newvar, ostatsn, StringLength)
1372
1373if obskind == 'trajectory':
1374    if dims.has_key('Z'):
1375        newdim = onewnc.createDimension('trj',3)
1376    else:
1377        newdim = onewnc.createDimension('trj',2)
1378
1379    newvar = onewnc.createVariable('obssimtrj','i',('obstime','trj'))
1380    basicvardef(newvar, 'obssimtrj', 'trajectory on the simulation grid', '-')
1381    newvar[:] = trajpos.transpose()
1382
1383if dims.has_key('Z'):
1384    newdim = onewnc.createDimension('simtrj',4)
1385    trjsim = np.zeros((4,Ncoindt), dtype=int)
1386    trjsimval = np.zeros((4,Ncoindt), dtype=np.float)
1387else:
1388    newdim = onewnc.createDimension('simtrj',3)
1389    trjsim = np.zeros((3,Ncoindt), dtype=int)
1390    trjsimval = np.zeros((3,Ncoindt), dtype=np.float)
1391
1392Nvars = len(valvars)
1393for ivar in range(Nvars):
1394    simobsvalues = []
1395
1396    varsimobs = valvars[ivar][0] + '_' + valvars[ivar][1]
1397    print '  ' + varsimobs + '... .. .'
1398
1399    if not oobs.variables.has_key(valvars[ivar][1]):
1400        print errormsg
1401        print '  ' + main + ": observations file has not '" + valvars[ivar][1] +     \
1402          "' !!"
1403        quit(-1)
1404
1405    if not osim.variables.has_key(valvars[ivar][0]):
1406        if not searchInlist(varNOcheck, valvars[ivar][0]):
1407            print errormsg
1408            print '  ' + main + ": simulation file has not '" + valvars[ivar][0] +   \
1409              "' !!"
1410            quit(-1)
1411        else:
1412            ovsim = compute_varNOcheck(osim, valvars[ivar][0])
1413    else:
1414        ovsim = osim.variables[valvars[ivar][0]]
1415
1416    for idn in ovsim.dimensions:
1417        if not searchInlist(simdims.values(),idn):
1418            print errormsg
1419            print '  ' + main + ": dimension '" + idn + "' of variable '" +          \
1420              valvars[ivar][0] + "' not provided as reference coordinate [X,Y,Z,T] !!"
1421            quit(-1)
1422
1423    ovobs = oobs.variables[valvars[ivar][1]]
1424    if searchInlist(ovobs.ncattrs(),'_FillValue'): 
1425        oFillValue = ovobs.getncattr('_FillValue')
1426
1427# Simulated values spatially around coincident times
1428    if dims.has_key('Z'):
1429        simobsSvalues = np.zeros((Ncoindt, Ngrid*2+1, Ngrid*2+1, Ngrid*2+1),         \
1430          dtype = np.float)
1431    else:
1432        simobsSvalues = np.zeros((Ncoindt, Ngrid*2+1, Ngrid*2+1), dtype = np.float)
1433
1434# Observed values temporally around coincident times
1435    simobsTvalues = {}
1436    simobsTtvalues = np.zeros((Ncoindt,2), dtype=np.float)
1437
1438    if obskind == 'multi-points':
1439        for it in range(Ncoindt):
1440            slicev = dims['X'][0] + ':' + str(trajpos[2,it]) + '|' +                 \
1441              dims['Y'][0]+ ':' + str(trajpos[1,it]) + '|' +                         \
1442              dims['T'][0]+ ':' + str(coindtvalues[it][0])
1443            slicevar, dimslice = slice_variable(ovsim, slicev)
1444            simobsvalues.append([ slicevar, ovobs[coindtvalues[it][1]]])
1445            slicev = dims['X'][0] + ':' + str(trajpos[2,it]-Ngrid) + '@' +           \
1446              str(trajpos[2,it]+Ngrid) + '|' + dims['Y'][0] + ':' +                  \
1447              str(trajpos[1,it]-Ngrid) + '@' + str(trajpos[1,it]+Ngrid) + '|' +      \
1448              dims['T'][0]+':'+str(coindtvalues[it][0])
1449            slicevar, dimslice = slice_variable(ovsim, slicev)
1450            simobsSvalues[it,:,:] = slicevar
1451
1452    elif obskind == 'single-station':
1453        for it in range(Ncoindt):
1454            ito = int(coindtvalues[it,1])
1455            if valdimsim.has_key('X') and valdimsim.has_key('Y'):
1456                slicev = dims['X'][0] + ':' + str(stationpos[1]) + '|' +             \
1457                  dims['Y'][0] + ':' + str(stationpos[0]) + '|' +                    \
1458                  dims['T'][0] + ':' + str(int(coindtvalues[it][0]))
1459            else:
1460                slicev = dims['T'][0] + ':' + str(int(coindtvalues[it][0]))
1461            slicevar, dimslice = slice_variable(ovsim, slicev)
1462            if ovobs[int(ito)] == oFillValue or ovobs[int(ito)] == '--':
1463                simobsvalues.append([ slicevar, fillValueF])
1464#            elif ovobs[int(ito)] != ovobs[int(ito)]:
1465#                simobsvalues.append([ slicevar, fillValueF])
1466            else:
1467                simobsvalues.append([ slicevar, ovobs[int(ito)]])
1468            if valdimsim.has_key('X') and valdimsim.has_key('Y'):
1469                slicev = dims['X'][0] + ':' + str(stationpos[1]-Ngrid) + '@' +       \
1470                  str(stationpos[1]+Ngrid+1) + '|' + dims['Y'][0] + ':' +            \
1471                  str(stationpos[0]-Ngrid) + '@' + str(stationpos[0]+Ngrid+1) + '|' +\
1472                  dims['T'][0] + ':' + str(int(coindtvalues[it,0]))
1473            else:
1474                slicev = dims['T'][0] + ':' + str(int(coindtvalues[it][0]))
1475            slicevar, dimslice = slice_variable(ovsim, slicev)
1476            simobsSvalues[it,:,:] = slicevar
1477
1478            if it == 0:
1479                itoi = 0
1480                itof = int(coindtvalues[it,1]) / 2
1481            elif it == Ncoindt-1:
1482                itoi = int( (ito + int(coindtvalues[it-1,1])) / 2)
1483                itof = int(coindtvalues[it,1])
1484            else:
1485                itod = int( (ito - int(coindtvalues[it-1,1])) / 2 ) 
1486                itoi = ito - itod
1487                itod = int( (int(coindtvalues[it+1,1]) - ito) / 2 )
1488                itof = ito + itod
1489
1490            slicev = dims['T'][1] + ':' + str(itoi) + '@' + str(itof + 1)
1491
1492            slicevar, dimslice = slice_variable(ovobs, slicev)
1493            simobsTvalues[str(it)] = slicevar
1494
1495            simobsTtvalues[it,0] = valdimobs['T'][itoi]
1496            simobsTtvalues[it,1] = valdimobs['T'][itof]
1497
1498    elif obskind == 'trajectory':
1499        if dims.has_key('Z'):
1500            for it in range(Ncoindt):
1501                ito = int(coindtvalues[it,1])
1502                if notfound[ito] == 0:
1503                    trajpos[2,ito] = index_mat(valdimsim['Z'][coindtvalues[it,0],:,  \
1504                      trajpos[1,ito],trajpos[0,ito]], valdimobs['Z'][ito])
1505                    slicev = dims['X'][0]+':'+str(trajpos[0,ito]) + '|' +            \
1506                      dims['Y'][0]+':'+str(trajpos[1,ito]) + '|' +                   \
1507                      dims['Z'][0]+':'+str(trajpos[2,ito]) + '|' +                   \
1508                      dims['T'][0]+':'+str(int(coindtvalues[it,0]))
1509                    slicevar, dimslice = slice_variable(ovsim, slicev)
1510                    simobsvalues.append([ slicevar, ovobs[int(ito)]])
1511                    minx = np.max([trajpos[0,ito]-Ngrid,0])
1512                    maxx = np.min([trajpos[0,ito]+Ngrid+1,ovsim.shape[3]])
1513                    miny = np.max([trajpos[1,ito]-Ngrid,0])
1514                    maxy = np.min([trajpos[1,ito]+Ngrid+1,ovsim.shape[2]])
1515                    minz = np.max([trajpos[2,ito]-Ngrid,0])
1516                    maxz = np.min([trajpos[2,ito]+Ngrid+1,ovsim.shape[1]])
1517
1518                    slicev = dims['X'][0] + ':' + str(minx) + '@' + str(maxx) + '|' +\
1519                      dims['Y'][0] + ':' + str(miny) + '@' + str(maxy) + '|' +       \
1520                      dims['Z'][0] + ':' + str(minz) + '@' + str(maxz) + '|' +       \
1521                      dims['T'][0] + ':' + str(int(coindtvalues[it,0]))
1522                    slicevar, dimslice = slice_variable(ovsim, slicev)
1523
1524                    sliceS = []
1525                    sliceS.append(it)
1526                    sliceS.append(slice(0,maxz-minz))
1527                    sliceS.append(slice(0,maxy-miny))
1528                    sliceS.append(slice(0,maxx-minx))
1529
1530                    simobsSvalues[tuple(sliceS)] = slicevar
1531                    if ivar == 0:
1532                        trjsim[0,it] = trajpos[0,ito]
1533                        trjsim[1,it] = trajpos[1,ito]
1534                        trjsim[2,it] = trajpos[2,ito]
1535                        trjsim[3,it] = coindtvalues[it,0]
1536                else:
1537                    simobsvalues.append([fillValueF, fillValueF])
1538                    simobsSvalues[it,:,:,:]= np.ones((Ngrid*2+1,Ngrid*2+1,Ngrid*2+1),\
1539                      dtype = np.float)*fillValueF
1540        else:
1541            for it in range(Ncoindt):
1542                if notfound[it] == 0:
1543                    ito = coindtvalues[it,1]
1544                    slicev = dims['X'][0]+':'+str(trajpos[2,ito]) + '|' +            \
1545                      dims['Y'][0]+':'+str(trajpos[1,ito]) + '|' +                   \
1546                      dims['T'][0]+':'+str(coindtvalues[ito,0])
1547                    slicevar, dimslice = slice_variable(ovsim, slicev)
1548                    simobsvalues.append([ slicevar, ovobs[coindtvalues[it,1]]])
1549                    slicev = dims['X'][0] + ':' + str(trajpos[0,it]-Ngrid) + '@' +   \
1550                      str(trajpos[0,it]+Ngrid) + '|' + dims['Y'][0] + ':' +          \
1551                      str(trajpos[1,it]-Ngrid) + '@' + str(trajpos[1,it]+Ngrid) +    \
1552                      '|' + dims['T'][0] + ':' + str(coindtvalues[it,0])
1553                    slicevar, dimslice = slice_variable(ovsim, slicev)
1554                    simobsSvalues[it,:,:] = slicevar
1555                else:
1556                    simobsvalues.append([fillValue, fillValue])
1557                    simobsSvalues[it,:,:] = np.ones((Ngrid*2+1,Ngrid*2+1),           \
1558                      dtype = np.float)*fillValueF
1559                print simobsvalues[varsimobs][:][it]
1560
1561    arrayvals = np.array(simobsvalues)
1562    if len(valvars[ivar]) > 2:
1563        const=np.float(valvars[ivar][3])
1564        if valvars[ivar][2] == 'sumc':
1565            simobsSvalues = simobsSvalues + const
1566            arrayvals[:,0] = arrayvals[:,0] + const
1567        elif valvars[ivar][2] == 'subc':
1568            simobsSvalues = simobsSvalues - const
1569            arrayvals[:,0] = arrayvals[:,0] - const
1570        elif valvars[ivar][2] == 'mulc':
1571            simobsSvalues = simobsSvalues * const
1572            arrayvals[:,0] = arrayvals[:,0] * const
1573        elif valvars[ivar][2] == 'divc':
1574            simobsSvalues = simobsSvalues / const
1575            arrayvals[:,0] = arrayvals[:,0] / const
1576        else:
1577            print errormsg
1578            print '  ' + fname + ": operation '" + valvars[ivar][2] + "' not ready!!"
1579            quit(-1)
1580
1581# statisics sim
1582    simstats = np.zeros((5), dtype=np.float)
1583    simstats[0] = np.min(arrayvals[:,0])
1584    simstats[1] = np.max(arrayvals[:,0])
1585    simstats[2] = np.mean(arrayvals[:,0])
1586    simstats[3] = np.mean(arrayvals[:,0]*arrayvals[:,0])
1587    simstats[4] = np.sqrt(simstats[3] - simstats[2]*simstats[2])
1588
1589# statisics obs
1590    obsmask = ma.masked_equal(arrayvals[:,1], fillValueF)
1591    obsmask2 = obsmask*obsmask
1592
1593    obsstats = np.zeros((5), dtype=np.float)
1594    obsstats[0] = obsmask.min()
1595    obsstats[1] = obsmask.max()
1596    obsstats[2] = obsmask.mean()
1597    obsstats[3] = obsmask2.mean()
1598    obsstats[4] = np.sqrt(obsstats[3] - obsstats[2]*obsstats[2])
1599
1600# Statistics sim-obs
1601    simobsstats = np.zeros((9), dtype=np.float)
1602    diffvals = np.zeros((Ncoindt), dtype=np.float)
1603
1604    diffvals = arrayvals[:,0] - obsmask
1605
1606    simobsstats[0] = simstats[0] - obsstats[0]
1607    simobsstats[1] = np.mean(arrayvals[:,0]*obsmask)
1608    simobsstats[2] = np.min(diffvals)
1609    simobsstats[3] = np.max(diffvals)
1610    simobsstats[4] = np.mean(diffvals)
1611    simobsstats[5] = np.mean(np.abs(diffvals))
1612    simobsstats[6] = np.sqrt(np.mean(diffvals*diffvals))
1613    simobsstats[7], simobsstats[8] = sts.pearsonr(arrayvals[:,0], arrayvals[:,1])
1614
1615# Statistics around sim values
1616    aroundstats = np.zeros((5,Ncoindt), dtype=np.float)
1617    for it in range(Ncoindt):
1618        aroundstats[0,it] = np.min(simobsSvalues[it,])
1619        aroundstats[1,it] = np.max(simobsSvalues[it,])
1620        aroundstats[2,it] = np.mean(simobsSvalues[it,])
1621        aroundstats[3,it] = np.mean(simobsSvalues[it,]*simobsSvalues[it,])
1622        aroundstats[4,it] = np.sqrt(aroundstats[3,it] - aroundstats[2,it]*           \
1623          aroundstats[2,it])
1624
1625# Statistics around obs values
1626    aroundostats = np.zeros((6,Ncoindt), dtype=np.float)
1627
1628    for it in range(Ncoindt):
1629        obsmask = ma.masked_equal(simobsTvalues[str(it)], fillValueF)
1630        obsmask2 = obsmask*obsmask
1631
1632        aroundostats[0,it] = len(obsmask.flatten())
1633        aroundostats[1,it] = obsmask.min()
1634        aroundostats[2,it] = obsmask.max()
1635        aroundostats[3,it] = obsmask.mean()
1636        aroundostats[4,it] = obsmask2.mean()
1637        aroundostats[5,it] = np.sqrt(aroundostats[4,it] - aroundostats[3,it]*        \
1638          aroundostats[3,it])
1639
1640# sim Values to netCDF
1641    newvar = onewnc.createVariable(valvars[ivar][0] + '_sim', 'f', ('time'),         \
1642      fill_value=fillValueF)
1643    descvar = 'simulated: ' + valvars[ivar][0]
1644    basicvardef(newvar, valvars[ivar][0], descvar, ovobs.getncattr('units'))
1645    newvar[:] = arrayvals[:,0]
1646
1647# obs Values to netCDF
1648    newvar = onewnc.createVariable(valvars[ivar][1] + '_obs', 'f', ('time'),         \
1649      fill_value=fillValueF)
1650    descvar = 'observed: ' + valvars[ivar][1]
1651    basicvardef(newvar, valvars[ivar][1], descvar, ovobs.getncattr('units'))
1652    newvar[:] = arrayvals[:,1]
1653
1654# Around values
1655    if not onewnc.variables.has_key(valvars[ivar][0] + 'around'):
1656        if dims.has_key('Z'):
1657            if not onewnc.dimensions.has_key('zaround'):
1658                newdim = onewnc.createDimension('zaround',Ngrid*2+1)
1659                newvar = onewnc.createVariable(valvars[ivar][0] + 'around', 'f',     \
1660                  ('time','zaround','yaround','xaround'), fill_value=fillValueF)
1661        else:
1662            newvar = onewnc.createVariable(valvars[ivar][0] + 'around', 'f',         \
1663              ('time','yaround','xaround'), fill_value=fillValueF)
1664
1665        descvar = 'around simulated values +/- grid values: ' + valvars[ivar][0]
1666        basicvardef(newvar, varsimobs + 'around', descvar, ovobs.getncattr('units'))
1667        newvar[:] = simobsSvalues
1668
1669# sim Statistics
1670    if not searchInlist(onewnc.variables,valvars[ivar][0] + 'stsim'):
1671        newvar = onewnc.createVariable(valvars[ivar][0] + 'stsim', 'f', ('stats'),   \
1672          fill_value=fillValueF)
1673        descvar = 'simulated statisitcs: ' + valvars[ivar][0]
1674        basicvardef(newvar, valvars[ivar][0] + 'stsim', descvar, ovobs.getncattr('units'))
1675        newvar[:] = simstats
1676
1677# obs Statistics
1678    if not searchInlist(onewnc.variables,valvars[ivar][1] + 'stobs'):
1679        newvar = onewnc.createVariable(valvars[ivar][1] + 'stobs', 'f', ('stats'),   \
1680          fill_value=fillValueF)
1681        descvar = 'observed statisitcs: ' + valvars[ivar][1]
1682        basicvardef(newvar, valvars[ivar][1] + 'stobs', descvar,                     \
1683          ovobs.getncattr('units'))
1684        newvar[:] = obsstats
1685
1686# sim-obs Statistics
1687    if not searchInlist(onewnc.variables,varsimobs + 'st'):
1688        newvar = onewnc.createVariable(varsimobs + 'st', 'f', ('gstats'),            \
1689          fill_value=fillValueF)
1690        descvar = 'simulated-observed statisitcs: ' + varsimobs
1691        basicvardef(newvar, varsimobs + 'st', descvar, ovobs.getncattr('units'))
1692        newvar[:] = simobsstats
1693
1694# around sim Statistics
1695    if not searchInlist(onewnc.variables,valvars[ivar][0] + 'staround'):
1696        newvar = onewnc.createVariable(valvars[ivar][0] + 'staround', 'f',           \
1697          ('time','stats'), fill_value=fillValueF)
1698        descvar = 'around (' +  str(Ngrid) + ', ' + str(Ngrid) +                     \
1699          ') simulated statisitcs: ' + valvars[ivar][0]
1700        basicvardef(newvar, valvars[ivar][0] + 'staround', descvar,                  \
1701          ovobs.getncattr('units'))
1702        newvar[:] = aroundstats.transpose()
1703
1704    if not searchInlist(onewnc.variables, 'time_bnds'):
1705        newvar = onewnc.createVariable('time_bnds','f8',('time','bnds'))
1706        basicvardef(newvar, 'time_bnds', 'time', obstunits )
1707        set_attribute(newvar, 'calendar', 'standard')
1708        newvar[:] = simobsTtvalues
1709
1710# around obs Statistics
1711    if not searchInlist(onewnc.variables,valvars[ivar][1] + 'staround'):
1712        newvar = onewnc.createVariable(valvars[ivar][1] + 'staround', 'f',           \
1713          ('time','tstats'), fill_value=fillValueF)
1714        descvar = 'around temporal observed statisitcs: ' + valvars[ivar][1]
1715        basicvardef(newvar, valvars[ivar][1] + 'staround', descvar,                  \
1716          ovobs.getncattr('units'))
1717        set_attribute(newvar, 'cell_methods', 'statistics')
1718
1719        newvar[:] = aroundostats.transpose()
1720
1721        onewnc.sync()
1722
1723newvar = onewnc.createVariable('simtrj','i',('time','simtrj'))
1724basicvardef(newvar, 'simtrj', 'coordinates [X,Y,Z,T] of the coincident trajectory ' +\
1725  'in sim', obstunits)
1726newvar[:] = trjsim.transpose()
1727
1728# Global attributes
1729##
1730set_attribute(onewnc,'author_nc','Lluis Fita')
1731set_attribute(onewnc,'institution_nc','Laboratoire de Meteorology Dynamique, ' +    \
1732  'LMD-Jussieu, UPMC, Paris')
1733set_attribute(onewnc,'country_nc','France')
1734set_attribute(onewnc,'script_nc',main)
1735set_attribute(onewnc,'version_script',version)
1736set_attribute(onewnc,'information',                                                 \
1737  'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html')
1738set_attribute(onewnc,'simfile',opts.fsim)
1739set_attribute(onewnc,'obsfile',opts.fobs)
1740
1741onewnc.sync()
1742onewnc.close()
1743
1744print main + ": successfull writting of '" + ofile + "' !!"
Note: See TracBrowser for help on using the repository browser.