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

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

Working version with the function and the addition of the station description
File size: 70.8 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
987def adding_station_desc(onc,stdesc):
988    """ Function to add a station description in a netCDF file
989      onc= netCDF object
990      stdesc= station description lon, lat, height
991    """
992    fname = 'adding_station_desc'
993
994    newvar = onc.createVariable( 'station', 'c', ('StrLength'))
995    newvar[0:len(stdesc[0])] = stdesc[0]
996
997    newdim = onc.createDimension('nst',1)
998
999    if onc.variables.has_key('lon'):
1000        print warnmsg
1001        print '  ' + fname + ": variable 'lon' already exist !!"
1002        print "    renaming it as 'lonst'"
1003        lonname = 'lonst'
1004    else:
1005        lonname = 'lon'
1006
1007    newvar = onc.createVariable( lonname, 'f4', ('nst'))
1008    basicvardef(newvar, lonname, 'longitude', 'degrees_West' )
1009    newvar[:] = stdesc[1]
1010
1011    if onc.variables.has_key('lat'):
1012        print warnmsg
1013        print '  ' + fname + ": variable 'lat' already exist !!"
1014        print "    renaming it as 'latst'"
1015        latname = 'latst'
1016    else:
1017        latname = 'lat'
1018
1019    newvar = onc.createVariable( latname, 'f4', ('nst'))
1020    basicvardef(newvar, lonname, 'latitude', 'degrees_North' )
1021    newvar[:] = stdesc[2]
1022
1023    if onc.variables.has_key('height'):
1024        print warnmsg
1025        print '  ' + fname + ": variable 'height' already exist !!"
1026        print "    renaming it as 'heightst'"
1027        heightname = 'heightst'
1028    else:
1029        heightname = 'height'
1030
1031    newvar = onc.createVariable( heightname, 'f4', ('nst'))
1032    basicvardef(newvar, heightname, 'height above sea level', 'm' )
1033    newvar[:] = stdesc[3]
1034
1035    return
1036
1037def getting_ValidationValues(okind, dimt, ds, trjpos, ovs, ovo, tvalues, oFill, Ng):
1038    """ Function to get the values to validate accroding to the type of observation
1039      okind= observational kind
1040      dimt= number of values to retrieve
1041      ds= dictionary with the names of the dimensions (sim, obs)
1042      trjpos= positions of the multi-stations (t, Y, X) or trajectory ([Z], Y, X)
1043      ovs= object with the values of the simulation
1044      ovs= object with the values of the observations
1045      tvalues= temporal values (sim. time step, obs. time step, sim t value, obs t value, t diff)
1046      oFill= Fill Value for the observations
1047      Ng= number of grid points around the observation
1048    return:
1049      sovalues= simulated values at the observation point and time
1050      soSvalues= values Ngrid points around the simulated point
1051      soTvalues= values around observed times (for `single-station')
1052      soTtvalues= inital/ending period between two consecutive obsevations (for `single-station')
1053      trjs= trajectory on the simulation space
1054    """
1055    fname = 'getting_ValidationValues'
1056
1057    sovalues = []
1058
1059# Simulated values spatially around
1060    if ds.has_key('Z'):
1061        soSvalues = np.zeros((dimt, Ng*2+1, Ng*2+1, Ng*2+1),            \
1062          dtype = np.float)
1063        if okind == 'trajectory':
1064            trjs = np.zeros((4,dimt), dtype=int)
1065        else:
1066            trjs = None
1067    else:
1068        soSvalues = np.zeros((dimt, Ng*2+1, Ng*2+1), dtype = np.float)
1069        if okind == 'trajectory':
1070            trjs = np.zeros((3,dimt), dtype=int)
1071        else:
1072            trjs = None
1073
1074    if okind == 'single-station':
1075        soTvalues = {}
1076        soTtvalues = np.zeros((dimt,2), dtype=np.float)
1077    else:
1078        None
1079
1080    if okind == 'multi-points':
1081        for it in range(dimt):
1082            slicev = ds['X'][0] + ':' + str(trjpos[2,it]) + '|' +                 \
1083              ds['Y'][0]+ ':' + str(trjpos[1,it]) + '|' +                         \
1084              ds['T'][0]+ ':' + str(tvalues[it][0])
1085            slicevar, dimslice = slice_variable(ovs, slicev)
1086            sovalues.append([ slicevar, ovo[tvalues[it][1]]])
1087            slicev = ds['X'][0] + ':' + str(trjpos[2,it]-Ng) + '@' +           \
1088              str(trjpos[2,it]+Ng) + '|' + ds['Y'][0] + ':' +                  \
1089              str(trjpos[1,it]-Ng) + '@' + str(trjpos[1,it]+Ng) + '|' +      \
1090              ds['T'][0]+':'+str(tvalues[it][0])
1091            slicevar, dimslice = slice_variable(ovs, slicev)
1092            soSvalues[it,:,:] = slicevar
1093
1094    elif okind == 'single-station':
1095        for it in range(dimt):
1096            ito = int(tvalues[it,1])
1097            if valdimsim.has_key('X') and valdimsim.has_key('Y'):
1098                slicev = ds['X'][0] + ':' + str(stationpos[1]) + '|' +             \
1099                  ds['Y'][0] + ':' + str(stationpos[0]) + '|' +                    \
1100                  ds['T'][0] + ':' + str(int(tvalues[it][0]))
1101            else:
1102                slicev = ds['T'][0] + ':' + str(int(tvalues[it][0]))
1103            slicevar, dimslice = slice_variable(ovs, slicev)
1104            if ovo[int(ito)] == oFill or ovo[int(ito)] == '--':
1105                sovalues.append([ slicevar, fillValueF])
1106#            elif ovo[int(ito)] != ovo[int(ito)]:
1107#                sovalues.append([ slicevar, fillValueF])
1108            else:
1109                sovalues.append([ slicevar, ovo[int(ito)]])
1110            if valdimsim.has_key('X') and valdimsim.has_key('Y'):
1111                slicev = ds['X'][0] + ':' + str(stationpos[1]-Ng) + '@' +       \
1112                  str(stationpos[1]+Ng+1) + '|' + ds['Y'][0] + ':' +            \
1113                  str(stationpos[0]-Ng) + '@' + str(stationpos[0]+Ng+1) + '|' +\
1114                  ds['T'][0] + ':' + str(int(tvalues[it,0]))
1115            else:
1116                slicev = ds['T'][0] + ':' + str(int(tvalues[it][0]))
1117            slicevar, dimslice = slice_variable(ovs, slicev)
1118            soSvalues[it,:,:] = slicevar
1119
1120            if it == 0:
1121                itoi = 0
1122                itof = int(tvalues[it,1]) / 2
1123            elif it == dimt-1:
1124                itoi = int( (ito + int(tvalues[it-1,1])) / 2)
1125                itof = int(tvalues[it,1])
1126            else:
1127                itod = int( (ito - int(tvalues[it-1,1])) / 2 ) 
1128                itoi = ito - itod
1129                itod = int( (int(tvalues[it+1,1]) - ito) / 2 )
1130                itof = ito + itod
1131
1132            slicev = ds['T'][1] + ':' + str(itoi) + '@' + str(itof + 1)
1133
1134            slicevar, dimslice = slice_variable(ovo, slicev)
1135            soTvalues[str(it)] = slicevar
1136
1137            soTtvalues[it,0] = valdimobs['T'][itoi]
1138            soTtvalues[it,1] = valdimobs['T'][itof]
1139
1140    elif okind == 'trajectory':
1141        if ds.has_key('Z'):
1142            for it in range(dimt):
1143                ito = int(tvalues[it,1])
1144                if notfound[ito] == 0:
1145                    trjpos[2,ito] = index_mat(valdimsim['Z'][tvalues[it,0],:,  \
1146                      trjpos[1,ito],trjpos[0,ito]], valdimobs['Z'][ito])
1147                    slicev = ds['X'][0]+':'+str(trjpos[0,ito]) + '|' +            \
1148                      ds['Y'][0]+':'+str(trjpos[1,ito]) + '|' +                   \
1149                      ds['Z'][0]+':'+str(trjpos[2,ito]) + '|' +                   \
1150                      ds['T'][0]+':'+str(int(tvalues[it,0]))
1151                    slicevar, dimslice = slice_variable(ovs, slicev)
1152                    sovalues.append([ slicevar, ovo[int(ito)]])
1153                    minx = np.max([trjpos[0,ito]-Ng,0])
1154                    maxx = np.min([trjpos[0,ito]+Ng+1,ovs.shape[3]])
1155                    miny = np.max([trjpos[1,ito]-Ng,0])
1156                    maxy = np.min([trjpos[1,ito]+Ng+1,ovs.shape[2]])
1157                    minz = np.max([trjpos[2,ito]-Ng,0])
1158                    maxz = np.min([trjpos[2,ito]+Ng+1,ovs.shape[1]])
1159
1160                    slicev = ds['X'][0] + ':' + str(minx) + '@' + str(maxx) + '|' +\
1161                      ds['Y'][0] + ':' + str(miny) + '@' + str(maxy) + '|' +       \
1162                      ds['Z'][0] + ':' + str(minz) + '@' + str(maxz) + '|' +       \
1163                      ds['T'][0] + ':' + str(int(tvalues[it,0]))
1164                    slicevar, dimslice = slice_variable(ovs, slicev)
1165
1166                    sliceS = []
1167                    sliceS.append(it)
1168                    sliceS.append(slice(0,maxz-minz))
1169                    sliceS.append(slice(0,maxy-miny))
1170                    sliceS.append(slice(0,maxx-minx))
1171
1172                    soSvalues[tuple(sliceS)] = slicevar
1173                    if ivar == 0:
1174                        trjs[0,it] = trjpos[0,ito]
1175                        trjs[1,it] = trjpos[1,ito]
1176                        trjs[2,it] = trjpos[2,ito]
1177                        trjs[3,it] = tvalues[it,0]
1178                else:
1179                    sovalues.append([fillValueF, fillValueF])
1180                    soSvalues[it,:,:,:]= np.ones((Ng*2+1,Ng*2+1,Ng*2+1),         \
1181                      dtype = np.float)*fillValueF
1182# 2D trajectory
1183        else:
1184            for it in range(dimt):
1185                if notfound[it] == 0:
1186                    ito = tvalues[it,1]
1187                    slicev = ds['X'][0]+':'+str(trjpos[2,ito]) + '|' +            \
1188                      ds['Y'][0]+':'+str(trjpos[1,ito]) + '|' +                   \
1189                      ds['T'][0]+':'+str(tvalues[ito,0])
1190                    slicevar, dimslice = slice_variable(ovs, slicev)
1191                    sovalues.append([ slicevar, ovo[tvalues[it,1]]])
1192                    slicev = ds['X'][0] + ':' + str(trjpos[0,it]-Ng) + '@' +   \
1193                      str(trjpos[0,it]+Ng) + '|' + ds['Y'][0] + ':' +          \
1194                      str(trjpos[1,it]-Ng) + '@' + str(trjpos[1,it]+Ng) +    \
1195                      '|' + ds['T'][0] + ':' + str(tvalues[it,0])
1196                    slicevar, dimslice = slice_variable(ovs, slicev)
1197                    soSvalues[it,:,:] = slicevar
1198                else:
1199                    sovalues.append([fillValue, fillValue])
1200                    soSvalues[it,:,:] = np.ones((Ng*2+1,Ng*2+1),           \
1201                      dtype = np.float)*fillValueF
1202                print sovalues[varsimobs][:][it]
1203    else:
1204        print errormsg
1205        print '  ' + fname + ": observatino kind '" + okind + "' not ready!!"
1206        quit(-1)
1207
1208
1209    return sovalues, soSvalues, soTvalues, soTtvalues, trjs
1210
1211
1212####### ###### ##### #### ### ## #
1213
1214strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " +  \
1215  " and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')"
1216
1217kindobs=['multi-points', 'single-station', 'trajectory']
1218strkObs="kind of observations; 'multi-points': multiple individual punctual obs " +  \
1219  "(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
1220  "'trajectory': following a trajectory"
1221simopers = ['sumc','subc','mulc','divc']
1222opersinf = 'sumc,[constant]: add [constant] to variables values; subc,[constant]: '+ \
1223  'substract [constant] to variables values; mulc,[constant]: multipy by ' +         \
1224  '[constant] to variables values; divc,[constant]: divide by [constant] to ' +      \
1225  'variables values'
1226varNOcheck = ['WRFdens', 'WRFght', 'WRFp', 'WRFrh', 'WRFT', 'WRFt', 'WRFtd', 'WRFws',\
1227  'WRFwd', 'WRFz']
1228varNOcheckinf = "'WRFdens': air density from WRF variables; 'WRFght': geopotential"+ \
1229  " height from WRF variables; 'WRFp': pressure from WRF variables; 'WRFrh': " +     \
1230  "relative humidty fom WRF variables; 'WRFT': CF-time from WRF variables'WRFt'; " + \
1231  " temperature from WRF variables; 'WRFtd': dew-point temperature from WRF " +      \
1232  "variables; 'WRFws': wind speed from WRF variables; 'WRFwd': wind speed " +        \
1233  "direction from WRF variables; 'WRFz': height from WRF variables"
1234
1235dimshelp = "[DIM]@[simdim]@[obsdim] ',' list of couples of dimensions names from " + \
1236  "each source ([DIM]='X','Y','Z','T'; None, no value)"
1237vardimshelp = "[DIM]@[simvardim]@[obsvardim] ',' list of couples of variables " +    \
1238  "names with dimensions values from each source ([DIM]='X','Y','Z','T'; None, " +   \
1239  "no value, WRFdiagnosted variables also available: " + varNOcheckinf + ")"
1240varshelp="[simvar]@[obsvar]@[[oper]@[val]] ',' list of couples of variables to " +   \
1241  "validate and if necessary operation and value operations: " + opersinf +          \
1242  "(WRFdiagnosted variables also available: " + varNOcheckinf + ")"
1243statsn = ['minimum', 'maximum', 'mean', 'mean2', 'standard deviation']
1244gstatsn = ['bias', 'simobs_mean', 'sim_obsmin', 'sim_obsmax', 'sim_obsmean', 'mae',  \
1245  'rmse', 'r_pearsoncorr', 'p_pearsoncorr']
1246ostatsn = ['number of points', 'minimum', 'maximum', 'mean', 'mean2',                \
1247  'standard deviation']
1248
1249parser = OptionParser()
1250parser.add_option("-d", "--dimensions", dest="dims", help=dimshelp, metavar="VALUES")
1251parser.add_option("-D", "--vardimensions", dest="vardims",
1252  help=vardimshelp, metavar="VALUES")
1253parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs, 
1254  help=strkObs, metavar="FILE")
1255parser.add_option("-l", "--stationLocation", dest="stloc", 
1256  help="name, longitude, latitude and height of the station (only for 'single-station')", 
1257  metavar="FILE")
1258parser.add_option("-o", "--observation", dest="fobs",
1259  help="observations file to validate", metavar="FILE")
1260parser.add_option("-s", "--simulation", dest="fsim",
1261  help="simulation file to validate", metavar="FILE")
1262parser.add_option("-t", "--trajectoryfile", dest="trajf",
1263  help="file with grid points of the trajectory in the simulation grid ('simtrj')", 
1264  metavar="FILE")
1265parser.add_option("-v", "--variables", dest="vars",
1266  help=varshelp, metavar="VALUES")
1267
1268(opts, args) = parser.parse_args()
1269
1270#######    #######
1271## MAIN
1272    #######
1273
1274ofile='validation_sim.nc'
1275
1276if opts.dims is None:
1277    print errormsg
1278    print '  ' + main + ': No list of dimensions are provided!!'
1279    print '    a ',' list of values X@[dimxsim]@[dimxobs],...,T@[dimtsim]@[dimtobs]'+\
1280      ' is needed'
1281    quit(-1)
1282else:
1283    simdims = {}
1284    obsdims = {}
1285    print main +': couple of dimensions _______'
1286    dims = {}
1287    ds = opts.dims.split(',')
1288    for d in ds:
1289        dsecs = d.split('@')
1290        if len(dsecs) != 3:
1291            print errormsg
1292            print '  ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
1293            print '    [DIM]@[dimnsim]@[dimnobs]'
1294            quit(-1)
1295        if dsecs[1] != 'None':
1296            dims[dsecs[0]] = [dsecs[1], dsecs[2]]
1297            simdims[dsecs[0]] = dsecs[1]
1298            obsdims[dsecs[0]] = dsecs[2]
1299
1300            print '  ',dsecs[0],':',dsecs[1],',',dsecs[2]
1301       
1302if opts.vardims is None:
1303    print errormsg
1304    print '  ' + main + ': No list of variables with dimension values are provided!!'
1305    print '    a ',' list of values X@[vardimxsim]@[vardimxobs],...,T@' +  \
1306      '[vardimtsim]@[vardimtobs] is needed'
1307    quit(-1)
1308else:
1309    print main +': couple of variable dimensions _______'
1310    vardims = {}
1311    ds = opts.vardims.split(',')
1312    for d in ds:
1313        dsecs = d.split('@')
1314        if len(dsecs) != 3:
1315            print errormsg
1316            print '  ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
1317            print '    [DIM]@[vardimnsim]@[vardimnobs]'
1318            quit(-1)
1319        if dsecs[1] != 'None':
1320            vardims[dsecs[0]] = [dsecs[1], dsecs[2]]
1321            print '  ',dsecs[0],':',dsecs[1],',',dsecs[2]
1322
1323if opts.obskind is None:
1324    print errormsg
1325    print '  ' + main + ': No kind of observations provided !!'
1326    quit(-1)
1327else:
1328    obskind = opts.obskind
1329    if obskind == 'single-station':
1330        if opts.stloc is None:
1331            print errormsg
1332            print '  ' + main + ': No station location provided !!'
1333            quit(-1)
1334        else:
1335            stationdesc = [opts.stloc.split(',')[0],                                 \
1336              np.float(opts.stloc.split(',')[1]), np.float(opts.stloc.split(',')[2]),\
1337              np.float(opts.stloc.split(',')[3])]
1338
1339if opts.fobs is None:
1340    print errormsg
1341    print '  ' + main + ': No observations file is provided!!'
1342    quit(-1)
1343else:
1344    if not os.path.isfile(opts.fobs):
1345        print errormsg
1346        print '   ' + main + ": observations file '" + opts.fobs + "' does not exist !!"
1347        quit(-1)
1348
1349if opts.fsim is None:
1350    print errormsg
1351    print '  ' + main + ': No simulation file is provided!!'
1352    quit(-1)
1353else:
1354    if not os.path.isfile(opts.fsim):
1355        print errormsg
1356        print '   ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
1357        quit(-1)
1358
1359if opts.vars is None:
1360    print errormsg
1361    print '  ' + main + ': No list of couples of variables is provided!!'
1362    print '    a ',' list of values [varsim]@[varobs],... is needed'
1363    quit(-1)
1364else:
1365    valvars = []
1366    vs = opts.vars.split(',')
1367    for v in vs:
1368        vsecs = v.split('@')
1369        if len(vsecs) < 2:
1370            print errormsg
1371            print '  ' + main + ': wrong number of values in:',vsecs,                \
1372              ' at least 2 are needed !!'
1373            print '    [varsim]@[varobs]@[[oper][val]]'
1374            quit(-1)
1375        if len(vsecs) > 2:
1376            if not searchInlist(simopers,vsecs[2]): 
1377                print errormsg
1378                print main + ": operation on simulation values '" + vsecs[2] +       \
1379                  "' not ready !!"
1380                quit(-1)
1381
1382        valvars.append(vsecs)
1383
1384# Openning observations trajectory
1385##
1386oobs = NetCDFFile(opts.fobs, 'r')
1387
1388valdimobs = {}
1389for dn in dims:
1390    print dn,':',dims[dn]
1391    if dims[dn][1] != 'None':
1392        if not oobs.dimensions.has_key(dims[dn][1]):
1393            print errormsg
1394            print '  ' + main + ": observations file does not have dimension '" +    \
1395              dims[dn][1] + "' !!"
1396            quit(-1)
1397        if vardims[dn][1] != 'None':
1398            if not oobs.variables.has_key(vardims[dn][1]):
1399                print errormsg
1400                print '  ' + main + ": observations file does not have varibale " +  \
1401                  "dimension '" + vardims[dn][1] + "' !!"
1402                quit(-1)
1403            valdimobs[dn] = oobs.variables[vardims[dn][1]][:]
1404    else:
1405        if dn == 'X':
1406            valdimobs[dn] = stationdesc[1]
1407        elif dn == 'Y':
1408            valdimobs[dn] = stationdesc[2]
1409        elif dn == 'Z':
1410            valdimobs[dn] = stationdesc[3]
1411
1412osim = NetCDFFile(opts.fsim, 'r')
1413
1414valdimsim = {}
1415for dn in dims:
1416    if dims[dn][0] != 'None':
1417        if not osim.dimensions.has_key(dims[dn][0]):
1418            print errormsg
1419            print '  ' + main + ": simulation file '" + opts.fsim +                  \
1420              "' does not have dimension '" + dims[dn][0] + "' !!"
1421            print '    it has: ',osim.dimensions
1422            quit(-1)
1423
1424        if not osim.variables.has_key(vardims[dn][0]) and                            \
1425          not searchInlist(varNOcheck,vardims[dn][0]):
1426            print errormsg
1427            print '  ' + main + ": simulation file '" + opts.fsim +                  \
1428              "' does not have varibale dimension '" + vardims[dn][0] + "' !!"
1429            print '    it has variables:',osim.variables
1430            quit(-1)
1431        if searchInlist(varNOcheck,vardims[dn][0]):
1432            valdimsim[dn] = compute_varNOcheck(osim, vardims[dn][0])
1433        else:
1434            valdimsim[dn] = osim.variables[vardims[dn][0]][:]
1435
1436# General characteristics
1437dimtobs = valdimobs['T'].shape[0]
1438dimtsim = valdimsim['T'].shape[0]
1439
1440print main +': observational time-steps:',dimtobs,'simulation:',dimtsim
1441
1442notfound = np.zeros((dimtobs), dtype=int)
1443
1444if obskind == 'multi-points':
1445    trajpos = np.zeros((2,dimt),dtype=int)
1446    for it in range(dimtobs):
1447        trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'],                    \
1448          [valdimobs['X'][it],valdimobss['Y'][it]])
1449elif obskind == 'single-station':
1450    trajpos = None
1451    stationpos = np.zeros((2), dtype=int)
1452    if valdimsim.has_key('X') and valdimsim.has_key('Y'):
1453        stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],[valdimobs['Y'],         \
1454          valdimobs['X']])
1455        iid = 0
1456        for idn in osim.variables[vardims['X'][0]].dimensions:
1457            if idn == dims['X'][0]:
1458                stationpos[1] = stsimpos[iid]
1459            elif idn == dims['Y'][0]:
1460                stationpos[0] = stsimpos[iid]
1461
1462            iid = iid + 1
1463        print main + ': station point in simulation:', stationpos
1464        print '    station position:',valdimobs['X'],',',valdimobs['Y']
1465        print '    simulation coord.:',valdimsim['X'][tuple(stsimpos)],',',          \
1466          valdimsim['Y'][tuple(stsimpos)]
1467    else:
1468        print main + ': validation with two time-series !!'
1469
1470elif obskind == 'trajectory':
1471    if opts.trajf is not None:
1472        if not os.path.isfile(opts.fsim):
1473            print errormsg
1474            print '   ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
1475            quit(-1)
1476        else:
1477            otrjf = NetCDFFile(opts.fsim, 'r')
1478            trajpos[0,:] = otrjf.variables['obssimtrj'][0]
1479            trajpos[1,:] = otrjf.variables['obssimtrj'][1]
1480            otrjf.close()
1481    else:
1482        if dims.has_key('Z'):
1483            trajpos = np.zeros((3,dimtobs),dtype=int)
1484            for it in range(dimtobs):
1485                if np.mod(it*100./dimtobs,10.) == 0.:
1486                    print '    trajectory done: ',it*100./dimtobs,'%'
1487                stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],                 \
1488                  [valdimobs['Y'][it],valdimobs['X'][it]])
1489                stationpos = np.zeros((2), dtype=int)
1490                iid = 0
1491                for idn in osim.variables[vardims['X'][0]].dimensions:
1492                    if idn == dims['X'][0]:
1493                        stationpos[1] = stsimpos[iid]
1494                    elif idn == dims['Y'][0]:
1495                        stationpos[0] = stsimpos[iid]
1496                    iid = iid + 1
1497                if stationpos[0] == 0 and stationpos[1] == 0: notfound[it] = 1
1498             
1499                trajpos[0,it] = stationpos[0]
1500                trajpos[1,it] = stationpos[1]
1501# In the simulation 'Z' varies with time ... non-hydrostatic model! ;)
1502#                trajpos[2,it] = index_mat(valdimsim['Z'][it,:,stationpos[0],         \
1503#                  stationpos[1]], valdimobs['Z'][it])
1504        else:
1505            trajpos = np.zeros((2,dimtobs),dtype=int)
1506            for it in range(dimtobs):
1507                stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],                 \
1508                  [valdimobs['Y'][it],valdimobss['X'][it]])
1509                stationpos = np.zeros((2), dtype=int)
1510                iid = 0
1511                for idn in osim.variables[vardims['X'][0]].dimensions:
1512                    if idn == dims['X'][0]:
1513                        stationpos[1] = stsimpos[iid]
1514                    elif idn == dims['Y'][0]:
1515                        stationpos[0] = stsimpos[iid]
1516                    iid = iid + 1
1517                if stationpos[0] == 0 or stationpos[1] == 0: notfound[it] = 1
1518
1519                trajpos[0,it] = stationspos[0]
1520                trajpos[1,it] = stationspos[1]
1521
1522        print main + ': not found',np.sum(notfound),'points of the trajectory'
1523
1524# Getting times
1525tobj = oobs.variables[vardims['T'][1]]
1526obstunits = tobj.getncattr('units')
1527if vardims['T'][0] == 'WRFT':
1528    tsim = valdimsim['T'][:]
1529    simtunits = 'seconds since 1949-12-01 00:00:00'
1530else:
1531    tsim = osim.variables[vardims['T'][0]][:]
1532    otsim = osim.variables[vardims['T'][0]]
1533    simtunits = otsim.getncattr('units')
1534
1535simobstimes = coincident_CFtimes(tsim, obstunits, simtunits)
1536
1537#
1538## Looking for exact/near times
1539#
1540
1541# Exact Coincident times
1542##
1543exacttvalues0 = []
1544for it in range(dimtsim):   
1545    ot = 0
1546    for ito in range(ot,dimtobs-1):
1547        if valdimobs['T'][ito] == simobstimes[it]:
1548            ot = ito
1549            exacttvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito]])
1550
1551exacttvalues = np.array(exacttvalues0, dtype=np.float)
1552Nexactt = len(exacttvalues[:,0])
1553
1554print main + ': found',Nexactt,'Temporal same values in simulation and observations'
1555# Sim Coincident times
1556##
1557coindtvalues0 = []
1558for it in range(dimtsim):   
1559    ot = 0
1560    for ito in range(ot,dimtobs-1):
1561        if valdimobs['T'][ito] < simobstimes[it] and valdimobs['T'][ito+1] >         \
1562          simobstimes[it]:
1563            ot = ito
1564            tdist = simobstimes[it] - valdimobs['T'][ito]
1565            coindtvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito],      \
1566              tdist])
1567
1568coindtvalues = np.array(coindtvalues0, dtype=np.float)
1569
1570Ncoindt = len(coindtvalues[:,0])
1571print main + ': found',Ncoindt,'Simulation time-interval (within consecutive ' +     \
1572  'observed times) coincident times between simulation and observations'
1573
1574if Ncoindt == 0:
1575    print warnmsg
1576    print main + ': no coincident times found !!'
1577    print '  stopping it'
1578    quit(-1)
1579
1580# Validating
1581##
1582
1583onewnc = NetCDFFile(ofile, 'w')
1584
1585# Dimensions
1586newdim = onewnc.createDimension('time',None)
1587newdim = onewnc.createDimension('bnds',2)
1588newdim = onewnc.createDimension('obstime',None)
1589newdim = onewnc.createDimension('couple',2)
1590newdim = onewnc.createDimension('StrLength',StringLength)
1591newdim = onewnc.createDimension('xaround',Ngrid*2+1)
1592newdim = onewnc.createDimension('yaround',Ngrid*2+1)
1593newdim = onewnc.createDimension('gstats',9)
1594newdim = onewnc.createDimension('stats',5)
1595newdim = onewnc.createDimension('tstats',6)
1596
1597# Variable dimensions
1598##
1599newvar = onewnc.createVariable('obstime','f8',('time'))
1600basicvardef(newvar, 'obstime', 'time observations', obstunits )
1601set_attribute(newvar, 'calendar', 'standard')
1602set_attribute(newvar, 'bounds', 'time_bnds')
1603newvar[:] = coindtvalues[:,3]
1604
1605newvar = onewnc.createVariable('couple', 'c', ('couple','StrLength'))
1606basicvardef(newvar, 'couple', 'couples of values', '-')
1607writing_str_nc(newvar, ['sim','obs'], StringLength)
1608
1609newvar = onewnc.createVariable('statistics', 'c', ('stats','StrLength'))
1610basicvardef(newvar, 'statistics', 'statitics from values', '-')
1611writing_str_nc(newvar, statsn, StringLength)
1612
1613newvar = onewnc.createVariable('gstatistics', 'c', ('gstats','StrLength'))
1614basicvardef(newvar, 'gstatistics', 'global statitics from values', '-')
1615writing_str_nc(newvar, gstatsn, StringLength)
1616
1617newvar = onewnc.createVariable('tstatistics', 'c', ('tstats','StrLength'))
1618basicvardef(newvar, 'tstatistics', 'statitics from values along time', '-')
1619writing_str_nc(newvar, ostatsn, StringLength)
1620
1621if obskind == 'trajectory':
1622    if dims.has_key('Z'):
1623        newdim = onewnc.createDimension('trj',3)
1624    else:
1625        newdim = onewnc.createDimension('trj',2)
1626
1627    newvar = onewnc.createVariable('obssimtrj','i',('obstime','trj'))
1628    basicvardef(newvar, 'obssimtrj', 'trajectory on the simulation grid', '-')
1629    newvar[:] = trajpos.transpose()
1630
1631if dims.has_key('Z'):
1632    newdim = onewnc.createDimension('simtrj',4)
1633else:
1634    newdim = onewnc.createDimension('simtrj',3)
1635
1636Nvars = len(valvars)
1637for ivar in range(Nvars):
1638    simobsvalues = []
1639
1640    varsimobs = valvars[ivar][0] + '_' + valvars[ivar][1]
1641    print '  ' + varsimobs + '... .. .'
1642
1643    if not oobs.variables.has_key(valvars[ivar][1]):
1644        print errormsg
1645        print '  ' + main + ": observations file has not '" + valvars[ivar][1] +     \
1646          "' !!"
1647        quit(-1)
1648
1649    if not osim.variables.has_key(valvars[ivar][0]):
1650        if not searchInlist(varNOcheck, valvars[ivar][0]):
1651            print errormsg
1652            print '  ' + main + ": simulation file has not '" + valvars[ivar][0] +   \
1653              "' !!"
1654            quit(-1)
1655        else:
1656            ovsim = compute_varNOcheck(osim, valvars[ivar][0])
1657    else:
1658        ovsim = osim.variables[valvars[ivar][0]]
1659
1660    for idn in ovsim.dimensions:
1661        if not searchInlist(simdims.values(),idn):
1662            print errormsg
1663            print '  ' + main + ": dimension '" + idn + "' of variable '" +          \
1664              valvars[ivar][0] + "' not provided as reference coordinate [X,Y,Z,T] !!"
1665            quit(-1)
1666
1667    ovobs = oobs.variables[valvars[ivar][1]]
1668    if searchInlist(ovobs.ncattrs(),'_FillValue'): 
1669        oFillValue = ovobs.getncattr('_FillValue')
1670    else:
1671        oFillValue = None
1672
1673# Observed values temporally around coincident times
1674    simobsvalues, simobsSvalues, simobsTvalues, simobsTtvalues, trjsim =             \
1675        getting_ValidationValues(obskind, Ncoindt, dims, trajpos, ovsim, ovobs,      \
1676        coindtvalues, oFillValue, Ngrid)
1677
1678# Re-arranging values...
1679    arrayvals = np.array(simobsvalues)
1680    if len(valvars[ivar]) > 2:
1681        const=np.float(valvars[ivar][3])
1682        if valvars[ivar][2] == 'sumc':
1683            simobsSvalues = simobsSvalues + const
1684            arrayvals[:,0] = arrayvals[:,0] + const
1685        elif valvars[ivar][2] == 'subc':
1686            simobsSvalues = simobsSvalues - const
1687            arrayvals[:,0] = arrayvals[:,0] - const
1688        elif valvars[ivar][2] == 'mulc':
1689            simobsSvalues = simobsSvalues * const
1690            arrayvals[:,0] = arrayvals[:,0] * const
1691        elif valvars[ivar][2] == 'divc':
1692            simobsSvalues = simobsSvalues / const
1693            arrayvals[:,0] = arrayvals[:,0] / const
1694        else:
1695            print errormsg
1696            print '  ' + fname + ": operation '" + valvars[ivar][2] + "' not ready!!"
1697            quit(-1)
1698
1699# statisics sim
1700    simstats = np.zeros((5), dtype=np.float)
1701    simstats[0] = np.min(arrayvals[:,0])
1702    simstats[1] = np.max(arrayvals[:,0])
1703    simstats[2] = np.mean(arrayvals[:,0])
1704    simstats[3] = np.mean(arrayvals[:,0]*arrayvals[:,0])
1705    simstats[4] = np.sqrt(simstats[3] - simstats[2]*simstats[2])
1706
1707# statisics obs
1708    obsmask = ma.masked_equal(arrayvals[:,1], fillValueF)
1709    obsmask2 = obsmask*obsmask
1710
1711    obsstats = np.zeros((5), dtype=np.float)
1712    obsstats[0] = obsmask.min()
1713    obsstats[1] = obsmask.max()
1714    obsstats[2] = obsmask.mean()
1715    obsstats[3] = obsmask2.mean()
1716    obsstats[4] = np.sqrt(obsstats[3] - obsstats[2]*obsstats[2])
1717
1718# Statistics sim-obs
1719    simobsstats = np.zeros((9), dtype=np.float)
1720    diffvals = np.zeros((Ncoindt), dtype=np.float)
1721
1722    diffvals = arrayvals[:,0] - obsmask
1723
1724    simobsstats[0] = simstats[0] - obsstats[0]
1725    simobsstats[1] = np.mean(arrayvals[:,0]*obsmask)
1726    simobsstats[2] = np.min(diffvals)
1727    simobsstats[3] = np.max(diffvals)
1728    simobsstats[4] = np.mean(diffvals)
1729    simobsstats[5] = np.mean(np.abs(diffvals))
1730    simobsstats[6] = np.sqrt(np.mean(diffvals*diffvals))
1731    simobsstats[7], simobsstats[8] = sts.pearsonr(arrayvals[:,0], arrayvals[:,1])
1732
1733# Statistics around sim values
1734    aroundstats = np.zeros((5,Ncoindt), dtype=np.float)
1735    for it in range(Ncoindt):
1736        aroundstats[0,it] = np.min(simobsSvalues[it,])
1737        aroundstats[1,it] = np.max(simobsSvalues[it,])
1738        aroundstats[2,it] = np.mean(simobsSvalues[it,])
1739        aroundstats[3,it] = np.mean(simobsSvalues[it,]*simobsSvalues[it,])
1740        aroundstats[4,it] = np.sqrt(aroundstats[3,it] - aroundstats[2,it]*           \
1741          aroundstats[2,it])
1742
1743# Statistics around obs values
1744    aroundostats = np.zeros((6,Ncoindt), dtype=np.float)
1745
1746    for it in range(Ncoindt):
1747        obsmask = ma.masked_equal(simobsTvalues[str(it)], fillValueF)
1748        obsmask2 = obsmask*obsmask
1749
1750        aroundostats[0,it] = len(obsmask.flatten())
1751        aroundostats[1,it] = obsmask.min()
1752        aroundostats[2,it] = obsmask.max()
1753        aroundostats[3,it] = obsmask.mean()
1754        aroundostats[4,it] = obsmask2.mean()
1755        aroundostats[5,it] = np.sqrt(aroundostats[4,it] - aroundostats[3,it]*        \
1756          aroundostats[3,it])
1757
1758# sim Values to netCDF
1759    newvar = onewnc.createVariable(valvars[ivar][0] + '_sim', 'f', ('time'),         \
1760      fill_value=fillValueF)
1761    descvar = 'simulated: ' + valvars[ivar][0]
1762    basicvardef(newvar, valvars[ivar][0], descvar, ovobs.getncattr('units'))
1763    newvar[:] = arrayvals[:,0]
1764
1765# obs Values to netCDF
1766    newvar = onewnc.createVariable(valvars[ivar][1] + '_obs', 'f', ('time'),         \
1767      fill_value=fillValueF)
1768    descvar = 'observed: ' + valvars[ivar][1]
1769    basicvardef(newvar, valvars[ivar][1], descvar, ovobs.getncattr('units'))
1770    newvar[:] = arrayvals[:,1]
1771
1772# Around values
1773    if not onewnc.variables.has_key(valvars[ivar][0] + 'around'):
1774        if dims.has_key('Z'):
1775            if not onewnc.dimensions.has_key('zaround'):
1776                newdim = onewnc.createDimension('zaround',Ngrid*2+1)
1777                newvar = onewnc.createVariable(valvars[ivar][0] + 'around', 'f',     \
1778                  ('time','zaround','yaround','xaround'), fill_value=fillValueF)
1779        else:
1780            newvar = onewnc.createVariable(valvars[ivar][0] + 'around', 'f',         \
1781              ('time','yaround','xaround'), fill_value=fillValueF)
1782
1783        descvar = 'around simulated values +/- grid values: ' + valvars[ivar][0]
1784        basicvardef(newvar, varsimobs + 'around', descvar, ovobs.getncattr('units'))
1785        newvar[:] = simobsSvalues
1786
1787# sim Statistics
1788    if not searchInlist(onewnc.variables,valvars[ivar][0] + 'stsim'):
1789        newvar = onewnc.createVariable(valvars[ivar][0] + 'stsim', 'f', ('stats'),   \
1790          fill_value=fillValueF)
1791        descvar = 'simulated statisitcs: ' + valvars[ivar][0]
1792        basicvardef(newvar, valvars[ivar][0] + 'stsim', descvar, ovobs.getncattr('units'))
1793        newvar[:] = simstats
1794
1795# obs Statistics
1796    if not searchInlist(onewnc.variables,valvars[ivar][1] + 'stobs'):
1797        newvar = onewnc.createVariable(valvars[ivar][1] + 'stobs', 'f', ('stats'),   \
1798          fill_value=fillValueF)
1799        descvar = 'observed statisitcs: ' + valvars[ivar][1]
1800        basicvardef(newvar, valvars[ivar][1] + 'stobs', descvar,                     \
1801          ovobs.getncattr('units'))
1802        newvar[:] = obsstats
1803
1804# sim-obs Statistics
1805    if not searchInlist(onewnc.variables,varsimobs + 'st'):
1806        newvar = onewnc.createVariable(varsimobs + 'st', 'f', ('gstats'),            \
1807          fill_value=fillValueF)
1808        descvar = 'simulated-observed statisitcs: ' + varsimobs
1809        basicvardef(newvar, varsimobs + 'st', descvar, ovobs.getncattr('units'))
1810        newvar[:] = simobsstats
1811
1812# around sim Statistics
1813    if not searchInlist(onewnc.variables,valvars[ivar][0] + 'staround'):
1814        newvar = onewnc.createVariable(valvars[ivar][0] + 'staround', 'f',           \
1815          ('time','stats'), fill_value=fillValueF)
1816        descvar = 'around (' +  str(Ngrid) + ', ' + str(Ngrid) +                     \
1817          ') simulated statisitcs: ' + valvars[ivar][0]
1818        basicvardef(newvar, valvars[ivar][0] + 'staround', descvar,                  \
1819          ovobs.getncattr('units'))
1820        newvar[:] = aroundstats.transpose()
1821
1822    if not searchInlist(onewnc.variables, 'time_bnds'):
1823        newvar = onewnc.createVariable('time_bnds','f8',('time','bnds'))
1824        basicvardef(newvar, 'time_bnds', 'time', obstunits )
1825        set_attribute(newvar, 'calendar', 'standard')
1826        newvar[:] = simobsTtvalues
1827
1828# around obs Statistics
1829    if not searchInlist(onewnc.variables,valvars[ivar][1] + 'staround'):
1830        newvar = onewnc.createVariable(valvars[ivar][1] + 'staround', 'f',           \
1831          ('time','tstats'), fill_value=fillValueF)
1832        descvar = 'around temporal observed statisitcs: ' + valvars[ivar][1]
1833        basicvardef(newvar, valvars[ivar][1] + 'staround', descvar,                  \
1834          ovobs.getncattr('units'))
1835        set_attribute(newvar, 'cell_methods', 'statistics')
1836
1837        newvar[:] = aroundostats.transpose()
1838
1839        onewnc.sync()
1840
1841if trjsim is not None:
1842    newvar = onewnc.createVariable('simtrj','i',('time','simtrj'))
1843    basicvardef(newvar,'simtrj','coordinates [X,Y,Z,T] of the coincident ' +         \
1844      'trajectory in sim', obstunits)
1845    newvar[:] = trjsim.transpose()
1846
1847# Adding three variables with the station location, longitude, latitude and height
1848if obskind == 'single-station':
1849    adding_station_desc(onewnc,stationdesc)
1850
1851# Global attributes
1852##
1853set_attribute(onewnc,'author_nc','Lluis Fita')
1854set_attribute(onewnc,'institution_nc','Laboratoire de Meteorology Dynamique, ' +    \
1855  'LMD-Jussieu, UPMC, Paris')
1856set_attribute(onewnc,'country_nc','France')
1857set_attribute(onewnc,'script_nc',main)
1858set_attribute(onewnc,'version_script',version)
1859set_attribute(onewnc,'information',                                                 \
1860  'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html')
1861set_attribute(onewnc,'simfile',opts.fsim)
1862set_attribute(onewnc,'obsfile',opts.fobs)
1863
1864onewnc.sync()
1865onewnc.close()
1866
1867print main + ": successfull writting of '" + ofile + "' !!"
Note: See TracBrowser for help on using the repository browser.