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

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

Fixing non-comprensible error!?

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