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

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

Adding check when simulated variable has dimensions not given as reference [X,Y,Z,T]
File size: 50.7 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
13
14main = 'validarion_sim.py'
15errormsg = 'ERROR -- errror -- ERROR -- error'
16warnmsg = 'WARNING -- warning -- WARNING -- warning'
17
18# version
19version=1.0
20
21# Filling values for floats, integer and string
22fillValueF = 1.e20
23fillValueI = -99999
24fillValueS = '---'
25
26StringLength = 50
27
28# Number of grid points to take as 'environment' around the observed point
29Ngrid = 1
30
31def searchInlist(listname, nameFind):
32    """ Function to search a value within a list
33    listname = list
34    nameFind = value to find
35    >>> searInlist(['1', '2', '3', '5'], '5')
36    True
37    """
38    for x in listname:
39      if x == nameFind:
40        return True
41    return False
42
43def set_attribute(ncvar, attrname, attrvalue):
44    """ Sets a value of an attribute of a netCDF variable. Removes previous attribute value if exists
45    ncvar = object netcdf variable
46    attrname = name of the attribute
47    attrvalue = value of the attribute
48    """
49    import numpy as np
50    from netCDF4 import Dataset as NetCDFFile
51
52    attvar = ncvar.ncattrs()
53    if searchInlist(attvar, attrname):
54        attr = ncvar.delncattr(attrname)
55
56    attr = ncvar.setncattr(attrname, attrvalue)
57
58    return ncvar
59
60def basicvardef(varobj, vstname, vlname, vunits):
61    """ Function to give the basic attributes to a variable
62    varobj= netCDF variable object
63    vstname= standard name of the variable
64    vlname= long name of the variable
65    vunits= units of the variable
66    """
67    attr = varobj.setncattr('standard_name', vstname)
68    attr = varobj.setncattr('long_name', vlname)
69    attr = varobj.setncattr('units', vunits)
70
71    return
72
73def writing_str_nc(varo, values, Lchar):
74    """ Function to write string values in a netCDF variable as a chain of 1char values
75    varo= netCDF variable object
76    values = list of values to introduce
77    Lchar = length of the string in the netCDF file
78    """
79
80    Nvals = len(values)
81    for iv in range(Nvals):   
82        stringv=values[iv] 
83        charvals = np.chararray(Lchar)
84        Lstr = len(stringv)
85        charvals[Lstr:Lchar] = ''
86
87        for ich in range(Lstr):
88            charvals[ich] = stringv[ich:ich+1]
89
90        varo[iv,:] = charvals
91
92    return
93
94def index_3mat(matA,matB,matC,val):
95    """ Function to provide the coordinates of a given value inside three matrix simultaneously
96    index_mat(matA,matB,matC,val)
97      matA= matrix with one set of values
98      matB= matrix with the other set of values
99      matB= matrix with the third set of values
100      val= triplet of values to search
101    >>> 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])
102    [2 1 1]
103    """
104    fname = 'index_3mat'
105
106    matAshape = matA.shape
107    matBshape = matB.shape
108    matCshape = matC.shape
109
110    for idv in range(len(matAshape)):
111        if matAshape[idv] != matBshape[idv]:
112            print errormsg
113            print '  ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv],  \
114              'and B:',matBshape[idv],'does not coincide!!'
115            quit(-1)
116        if matAshape[idv] != matCshape[idv]:
117            print errormsg
118            print '  ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv],  \
119              'and C:',matCshape[idv],'does not coincide!!'
120            quit(-1)
121
122    minA = np.min(matA)
123    maxA = np.max(matA)
124    minB = np.min(matB)
125    maxB = np.max(matB)
126    minC = np.min(matC)
127    maxC = np.max(matC)
128
129    if val[0] < minA or val[0] > maxA:
130        print warnmsg
131        print '  ' + fname + ': first value:',val[0],'outside matA range',minA,',',  \
132          maxA,'!!'
133    if val[1] < minB or val[1] > maxB:
134        print warnmsg
135        print '  ' + fname + ': second value:',val[1],'outside matB range',minB,',',  \
136          maxB,'!!'
137    if val[2] < minC or val[2] > maxC:
138        print warnmsg
139        print '  ' + fname + ': second value:',val[2],'outside matC range',minC,',',  \
140          maxC,'!!'
141
142    dist = np.zeros(tuple(matAshape), dtype=np.float)
143    dist = np.sqrt((matA - np.float(val[0]))**2 + (matB - np.float(val[1]))**2 +     \
144      (matC - np.float(val[2]))**2)
145
146    mindist = np.min(dist)
147   
148    matlist = list(dist.flatten())
149    ifound = matlist.index(mindist)
150
151    Ndims = len(matAshape)
152    valpos = np.zeros((Ndims), dtype=int)
153    baseprevdims = np.zeros((Ndims), dtype=int)
154
155    for dimid in range(Ndims):
156        baseprevdims[dimid] = np.product(matAshape[dimid+1:Ndims])
157        if dimid == 0:
158            alreadyplaced = 0
159        else:
160            alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
161        valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
162
163    return valpos
164
165def index_2mat(matA,matB,val):
166    """ Function to provide the coordinates of a given value inside two matrix simultaneously
167    index_mat(matA,matB,val)
168      matA= matrix with one set of values
169      matB= matrix with the pother set of values
170      val= couple of values to search
171    >>> index_2mat(np.arange(27).reshape(3,3,3),np.arange(100,127).reshape(3,3,3),[22,111])
172    [2 1 1]
173    """
174    fname = 'index_2mat'
175
176    matAshape = matA.shape
177    matBshape = matB.shape
178
179    for idv in range(len(matAshape)):
180        if matAshape[idv] != matBshape[idv]:
181            print errormsg
182            print '  ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv],  \
183              'and B:',matBshape[idv],'does not coincide!!'
184            quit(-1)
185
186    minA = np.min(matA)
187    maxA = np.max(matA)
188    minB = np.min(matB)
189    maxB = np.max(matB)
190
191    Ndims = len(matAshape)
192#    valpos = np.ones((Ndims), dtype=int)*-1.
193    valpos = np.zeros((Ndims), dtype=int)
194
195    if val[0] < minA or val[0] > maxA:
196        print warnmsg
197        print '  ' + fname + ': first value:',val[0],'outside matA range',minA,',',  \
198          maxA,'!!'
199        return valpos
200    if val[1] < minB or val[1] > maxB:
201        print warnmsg
202        print '  ' + fname + ': second value:',val[1],'outside matB range',minB,',',  \
203          maxB,'!!'
204        return valpos
205
206    dist = np.zeros(tuple(matAshape), dtype=np.float)
207    dist = np.sqrt((matA - np.float(val[0]))**2 + (matB - np.float(val[1]))**2)
208
209    mindist = np.min(dist)
210   
211    if mindist != mindist:
212        print '  ' + fname + ': wrong minimal distance',mindist,'!!'
213        return valpos
214    else:
215        matlist = list(dist.flatten())
216        ifound = matlist.index(mindist)
217
218    baseprevdims = np.zeros((Ndims), dtype=int)
219    for dimid in range(Ndims):
220        baseprevdims[dimid] = np.product(matAshape[dimid+1:Ndims])
221        if dimid == 0:
222            alreadyplaced = 0
223        else:
224            alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
225        valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
226
227    return valpos
228
229def index_mat(matA,val):
230    """ Function to provide the coordinates of a given value inside a matrix
231    index_mat(matA,val)
232      matA= matrix with one set of values
233      val= couple of values to search
234    >>> index_mat(np.arange(27),22.3)
235    22
236    """
237    fname = 'index_mat'
238
239    matAshape = matA.shape
240
241    minA = np.min(matA)
242    maxA = np.max(matA)
243
244    Ndims = len(matAshape)
245#    valpos = np.ones((Ndims), dtype=int)*-1.
246    valpos = np.zeros((Ndims), dtype=int)
247
248    if val < minA or val > maxA:
249        print warnmsg
250        print '  ' + fname + ': first value:',val,'outside matA range',minA,',',     \
251          maxA,'!!'
252        return valpos
253
254    dist = np.zeros(tuple(matAshape), dtype=np.float)
255    dist = (matA - np.float(val))**2
256
257    mindist = np.min(dist)
258    if mindist != mindist:
259        print '  ' + fname + ': wrong minimal distance',mindist,'!!'
260        return valpos
261   
262    matlist = list(dist.flatten())
263    valpos = matlist.index(mindist)
264
265    return valpos
266
267def index_mat_exact(mat,val):
268    """ Function to provide the coordinates of a given exact value inside a matrix
269    index_mat(mat,val)
270      mat= matrix with values
271      val= value to search
272    >>> index_mat(np.arange(27).reshape(3,3,3),22)
273    [2 1 1]
274    """
275
276    fname = 'index_mat'
277
278    matshape = mat.shape
279
280    matlist = list(mat.flatten())
281    ifound = matlist.index(val)
282
283    Ndims = len(matshape)
284    valpos = np.zeros((Ndims), dtype=int)
285    baseprevdims = np.zeros((Ndims), dtype=int)
286
287    for dimid in range(Ndims):
288        baseprevdims[dimid] = np.product(matshape[dimid+1:Ndims])
289        if dimid == 0:
290            alreadyplaced = 0
291        else:
292            alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
293        valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
294
295    return valpos
296
297def coincident_CFtimes(tvalB, tunitA, tunitB):
298    """ Function to make coincident times for two different sets of CFtimes
299    tvalB= time values B
300    tunitA= time units times A to which we want to make coincidence
301    tunitB= time units times B
302    >>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00',
303      'hours since 1949-12-01 00:00:00')
304    [     0.   3600.   7200.  10800.  14400.  18000.  21600.  25200.  28800.  32400.]
305    >>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00',
306      'hours since 1979-12-01 00:00:00')
307    [  9.46684800e+08   9.46688400e+08   9.46692000e+08   9.46695600e+08
308       9.46699200e+08   9.46702800e+08   9.46706400e+08   9.46710000e+08
309       9.46713600e+08   9.46717200e+08]
310    """
311    import datetime as dt
312    fname = 'coincident_CFtimes'
313
314    trefA = tunitA.split(' ')[2] + ' ' + tunitA.split(' ')[3]
315    trefB = tunitB.split(' ')[2] + ' ' + tunitB.split(' ')[3]
316    tuA = tunitA.split(' ')[0]
317    tuB = tunitB.split(' ')[0]
318
319    if tuA != tuB:
320        if tuA == 'microseconds':
321            if tuB == 'microseconds':
322                tB = tvalB*1.
323            elif tuB == 'seconds':
324                tB = tvalB*10.e6
325            elif tuB == 'minutes':
326                tB = tvalB*60.*10.e6
327            elif tuB == 'hours':
328                tB = tvalB*3600.*10.e6
329            elif tuB == 'days':
330                tB = tvalB*3600.*24.*10.e6
331            else:
332                print errormsg
333                print '  ' + fname + ": combination of time untis: '" + tuA +        \
334                  "' & '" + tuB + "' not ready !!"
335                quit(-1)
336        elif tuA == 'seconds':
337            if tuB == 'microseconds':
338                tB = tvalB/10.e6
339            elif tuB == 'seconds':
340                tB = tvalB*1.
341            elif tuB == 'minutes':
342                tB = tvalB*60.
343            elif tuB == 'hours':
344                tB = tvalB*3600.
345            elif tuB == 'days':
346                tB = tvalB*3600.*24.
347            else:
348                print errormsg
349                print '  ' + fname + ": combination of time untis: '" + tuA +        \
350                  "' & '" + tuB + "' not ready !!"
351                quit(-1)
352        elif tuA == 'minutes':
353            if tuB == 'microseconds':
354                tB = tvalB/(60.*10.e6)
355            elif tuB == 'seconds':
356                tB = tvalB/60.
357            elif tuB == 'minutes':
358                tB = tvalB*1.
359            elif tuB == 'hours':
360                tB = tvalB*60.
361            elif tuB == 'days':
362                tB = tvalB*60.*24.
363            else:
364                print errormsg
365                print '  ' + fname + ": combination of time untis: '" + tuA +        \
366                  "' & '" + tuB + "' not ready !!"
367                quit(-1)
368        elif tuA == 'hours':
369            if tuB == 'microseconds':
370                tB = tvalB/(3600.*10.e6)
371            elif tuB == 'seconds':
372                tB = tvalB/3600.
373            elif tuB == 'minutes':
374                tB = tvalB/60.
375            elif tuB == 'hours':
376                tB = tvalB*1.
377            elif tuB == 'days':
378                tB = tvalB*24.
379            else:
380                print errormsg
381                print '  ' + fname + ": combination of time untis: '" + tuA +        \
382                  "' & '" + tuB + "' not ready !!"
383                quit(-1)
384        elif tuA == 'days':
385            if tuB == 'microseconds':
386                tB = tvalB/(24.*3600.*10.e6)
387            elif tuB == 'seconds':
388                tB = tvalB/(24.*3600.)
389            elif tuB == 'minutes':
390                tB = tvalB/(24.*60.)
391            elif tuB == 'hours':
392                tB = tvalB/24.
393            elif tuB == 'days':
394                tB = tvalB*1.
395            else:
396                print errormsg
397                print '  ' + fname + ": combination of time untis: '" + tuA +        \
398                  "' & '" + tuB + "' not ready !!"
399                quit(-1)
400        else:
401            print errormsg
402            print '  ' + fname + ": time untis: '" + tuA + "' not ready !!"
403            quit(-1)
404    else:
405        tB = tvalB*1.
406
407    if trefA != trefB:
408        trefTA = dt.datetime.strptime(trefA, '%Y-%m-%d %H:%M:%S')
409        trefTB = dt.datetime.strptime(trefB, '%Y-%m-%d %H:%M:%S')
410
411        difft = trefTB - trefTA
412        diffv = difft.days*24.*3600.*10.e6 + difft.seconds*10.e6 + difft.microseconds
413        print '  ' + fname + ': different reference refA:',trefTA,'refB',trefTB
414        print '    difference:',difft,':',diffv,'microseconds'
415
416        if tuA == 'microseconds':
417            tB = tB + diffv
418        elif tuA == 'seconds':
419            tB = tB + diffv/10.e6
420        elif tuA == 'minutes':
421            tB = tB + diffv/(60.*10.e6)
422        elif tuA == 'hours':
423            tB = tB + diffv/(3600.*10.e6)
424        elif tuA == 'dayss':
425            tB = tB + diffv/(24.*3600.*10.e6)
426        else:
427            print errormsg
428            print '  ' + fname + ": time untis: '" + tuA + "' not ready !!"
429            quit(-1)
430
431    return tB
432
433
434def slice_variable(varobj, dimslice):
435    """ Function to return a slice of a given variable according to values to its
436      dimensions
437    slice_variable(varobj, dimslice)
438      varobj= object wit the variable
439      dimslice= [[dimname1]:[value1]|[[dimname2]:[value2], ...] pairs of dimension
440        [value]:
441          * [integer]: which value of the dimension
442          * -1: all along the dimension
443          * -9: last value of the dimension
444          * [beg]@[end] slice from [beg] to [end]
445    """
446    fname = 'slice_variable'
447
448    if varobj == 'h':
449        print fname + '_____________________________________________________________'
450        print slice_variable.__doc__
451        quit()
452
453    vardims = varobj.dimensions
454    Ndimvar = len(vardims)
455
456    Ndimcut = len(dimslice.split('|'))
457    dimsl = dimslice.split('|')
458
459    varvalsdim = []
460    dimnslice = []
461
462    for idd in range(Ndimvar):
463        for idc in range(Ndimcut):
464            dimcutn = dimsl[idc].split(':')[0]
465            dimcutv = dimsl[idc].split(':')[1]
466            if vardims[idd] == dimcutn: 
467                posfrac = dimcutv.find('@')
468                if posfrac != -1:
469                    inifrac = int(dimcutv.split('@')[0])
470                    endfrac = int(dimcutv.split('@')[1])
471                    varvalsdim.append(slice(inifrac,endfrac))
472                    dimnslice.append(vardims[idd])
473                else:
474                    if int(dimcutv) == -1:
475                        varvalsdim.append(slice(0,varobj.shape[idd]))
476                        dimnslice.append(vardims[idd])
477                    elif int(dimcutv) == -9:
478                        varvalsdim.append(int(varobj.shape[idd])-1)
479                    else:
480                        varvalsdim.append(int(dimcutv))
481                break
482
483    varvalues = varobj[tuple(varvalsdim)]
484
485    return varvalues, dimnslice
486
487def func_compute_varNOcheck(ncobj, varn):
488    """ Function to compute variables which are not originary in the file
489      ncobj= netCDF object file
490      varn = variable to compute:
491        'WRFdens': air density from WRF variables
492        'WRFght': geopotential height from WRF variables
493        'WRFp': pressure from WRF variables
494        'WRFrh': relative humidty fom WRF variables
495        'WRFt': temperature from WRF variables
496        'WRFz': height from WRF variables
497    """
498    fname = 'compute_varNOcheck'
499
500    if varn == 'WRFdens':
501#        print '    ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \
502#          'DNW)/g ...'
503        grav = 9.81
504
505# Just we need in in absolute values: Size of the central grid cell
506##    dxval = ncobj.getncattr('DX')
507##    dyval = ncobj.getncattr('DY')
508##    mapfac = ncobj.variables['MAPFAC_M'][:]
509##    area = dxval*dyval*mapfac
510        dimensions = ncobj.variables['MU'].dimensions
511
512        mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
513        dnw = ncobj.variables['DNW'][:]
514
515        varNOcheckv = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \
516          dtype=np.float)
517        levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
518
519        for it in range(mu.shape[0]):
520            for iz in range(dnw.shape[1]):
521                levval.fill(np.abs(dnw[it,iz]))
522                varNOcheck[it,iz,:,:] = levval
523                varNOcheck[it,iz,:,:] = mu[it,:,:]*varNOcheck[it,iz,:,:]/grav
524
525    elif varn == 'WRFght':
526#        print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...'
527        varNOcheckv = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
528        dimensions = ncobj.variables['PH'].dimensions
529
530    elif varn == 'WRFp':
531#        print '  ' + fname + ': Retrieving pressure value from WRF as P + PB'
532        varNOcheckv = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
533        dimensions = ncobj.variables['P'].dimensions
534
535    elif varn == 'WRFrh':
536#        print '    ' + main + ": computing relative humidity from WRF as 'Tetens'" +\
537#         ' equation (T,P) ...'
538        p0=100000.
539        p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
540        tk = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
541        qv = ncobj.variables['QVAPOR'][:]
542
543        data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
544        data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
545
546        varNOcheckv = qv/data2
547        dimensions = ncobj.variables['P'].dimensions
548
549    elif varn == 'WRFt':
550#        print '    ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
551        p0=100000.
552        p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
553
554        varNOcheckv = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
555        dimensions = ncobj.variables['T'].dimensions
556
557    elif varn == 'WRFz':
558        grav = 9.81
559#        print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...'
560        varNOcheckv = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/grav
561        dimensions = ncobj.variables['PH'].dimensions
562
563    else:
564        print erromsg
565        print '  ' + fname + ": variable '" + varn + "' nor ready !!"
566        quit(-1)
567
568    return varNOcheck
569
570class compute_varNOcheck(object):
571    """ Class to compute variables which are not originary in the file
572      ncobj= netCDF object file
573      varn = variable to compute:
574        'WRFdens': air density from WRF variables
575        'WRFght': geopotential height from WRF variables
576        'WRFp': pressure from WRF variables
577        'WRFrh': relative humidty fom WRF variables
578        'WRFt': temperature from WRF variables
579        'WRFtd': dew-point temperature from WRF variables
580        'WRFws': wind speed from WRF variables
581        'WRFwd': wind direction from WRF variables
582        'WRFz': height from WRF variables
583    """
584    fname = 'compute_varNOcheck'
585
586    def __init__(self, ncobj, varn):
587
588        if ncobj is None:
589            self = None
590            self.dimensions = None
591            self.shape = None
592            self.__values = None
593        else:
594            if varn == 'WRFdens':
595#        print '    ' + main + ': computing air density from WRF as ((MU + MUB) * ' + \
596#          'DNW)/g ...'
597                grav = 9.81
598
599# Just we need in in absolute values: Size of the central grid cell
600##    dxval = ncobj.getncattr('DX')
601##    dyval = ncobj.getncattr('DY')
602##    mapfac = ncobj.variables['MAPFAC_M'][:]
603##    area = dxval*dyval*mapfac
604                dimensions = ncobj.variables['MU'].dimensions
605                shape = ncobj.variables['MU'].shape
606
607                mu = (ncobj.variables['MU'][:] + ncobj.variables['MUB'][:])
608                dnw = ncobj.variables['DNW'][:]
609
610                varNOcheckv = np.zeros((mu.shape[0], dnw.shape[1], mu.shape[1], mu.shape[2]), \
611                  dtype=np.float)
612                levval = np.zeros((mu.shape[1], mu.shape[2]), dtype=np.float)
613
614                for it in range(mu.shape[0]):
615                    for iz in range(dnw.shape[1]):
616                        levval.fill(np.abs(dnw[it,iz]))
617                        varNOcheck[it,iz,:,:] = levval
618                        varNOcheck[it,iz,:,:] = mu[it,:,:]*varNOcheck[it,iz,:,:]/grav
619
620            elif varn == 'WRFght':
621#        print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...'
622                varNOcheckv = ncobj.variables['PH'][:] + ncobj.variables['PHB'][:]
623                dimensions = ncobj.variables['PH'].dimensions
624                shape = ncobj.variables['PH'].shape
625
626            elif varn == 'WRFp':
627#        print '  ' + fname + ': Retrieving pressure value from WRF as P + PB'
628                varNOcheckv = ncobj.variables['P'][:] + ncobj.variables['PB'][:]
629                dimensions = ncobj.variables['P'].dimensions
630                shape = ncobj.variables['P'].shape
631
632            elif varn == 'WRFrh':
633#        print '    ' + main + ": computing relative humidity from WRF as 'Tetens'" +\
634#         ' equation (T,P) ...'
635                p0=100000.
636                p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
637                tk = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
638                qv = ncobj.variables['QVAPOR'][:]
639
640                data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
641                data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
642
643                varNOcheckv = qv/data2
644                dimensions = ncobj.variables['P'].dimensions
645                shape = ncobj.variables['P'].shape
646
647            elif varn == 'WRFt':
648#        print '    ' + main + ': computing temperature from WRF as inv_potT(T + 300) ...'
649                p0=100000.
650                p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
651
652                varNOcheckv = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
653                dimensions = ncobj.variables['T'].dimensions
654                shape = ncobj.variables['P'].shape
655
656            elif varn == 'WRFtd':
657#        print '    ' + main + ': computing dew-point temperature from WRF as inv_potT(T + 300) and Tetens...'
658# tacking from: http://en.wikipedia.org/wiki/Dew_point
659                p0=100000.
660                p=ncobj.variables['P'][:] + ncobj.variables['PB'][:]
661
662                temp = (ncobj.variables['T'][:] + 300.)*(p/p0)**(2./7.)
663
664                qv = ncobj.variables['QVAPOR'][:]
665
666                tk = temp - 273.15
667                data1 = 10.*0.6112*np.exp(17.67*(tk-273.16)/(tk-29.65))
668                data2 = 0.622*data1/(0.01*p-(1.-0.622)*data1)
669
670                rh = qv/data2
671               
672                pa = rh * data1/100.
673                varNOcheckv = 257.44*np.log(pa/6.1121)/(18.678-np.log(pa/6.1121))
674
675                dimensions = ncobj.variables['T'].dimensions
676                shape = ncobj.variables['P'].shape
677
678            elif varn == 'WRFws':
679#        print '    ' + main + ': computing wind speed from WRF as SQRT(U**2 + V**2) ...'
680                uwind = ncobj.variables['U'][:]
681                vwind = ncobj.variables['V'][:]
682                dx = uwind.shape[3]
683                dy = vwind.shape[2]
684                 
685# de-staggering
686                ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
687                va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])
688
689                varNOcheckv = np.sqrt(ua*ua + va*va)
690                dimensions = tuple(['Time','bottom_top','south_north','west_east'])
691                shape = ua.shape
692
693            elif varn == 'WRFwd':
694#        print '    ' + main + ': computing wind direction from WRF as ATAN2PI(V,U) ...'
695                uwind = ncobj.variables['U'][:]
696                vwind = ncobj.variables['V'][:]
697                dx = uwind.shape[3]
698                dy = vwind.shape[2]
699
700# de-staggering
701                ua = 0.5*(uwind[:,:,:,0:dx-1] + uwind[:,:,:,1:dx])
702                va = 0.5*(vwind[:,:,0:dy-1,:] + vwind[:,:,1:dy,:])
703
704                theta = np.arctan2(ua, va)
705                dimensions = tuple(['Time','bottom_top','south_north','west_east'])
706                shape = ua.shape
707                varNOcheckv = 360.*(1. + theta/(2.*np.pi))
708
709                dimensions = ncobj.variables['U'].dimensions
710                shape = ncobj.variables['U'].shape
711
712            elif varn == 'WRFz':
713                grav = 9.81
714#        print '    ' + main + ': computing geopotential height from WRF as PH + PHB ...'
715                varNOcheckv = (ncobj.variables['PH'][:] + ncobj.variables['PHB'][:])/grav
716                dimensions = ncobj.variables['PH'].dimensions
717                shape = ncobj.variables['PH'].shape
718
719            else:
720                print errormsg
721                print '  ' + fname + ": variable '" + varn + "' nor ready !!"
722                quit(-1)
723
724            self.dimensions = dimensions
725            self.shape = shape
726            self.__values = varNOcheckv
727
728    def __getitem__(self,elem):
729        return self.__values[elem]
730
731####### ###### ##### #### ### ## #
732
733strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " +  \
734  " and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')"
735
736kindobs=['multi-points', 'single-station', 'trajectory']
737strkObs="kind of observations; 'multi-points': multiple individual punctual obs " +  \
738  "(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
739  "'trajectory': following a trajectory"
740simopers = ['sumc','subc','mulc','divc']
741opersinf = 'sumc,[constant]: add [constant] to variables values; subc,[constant]: '+ \
742  'substract [constant] to variables values; mulc,[constant]: multipy by ' +         \
743  '[constant] to variables values; divc,[constant]: divide by [constant] to ' +      \
744  'variables values'
745varNOcheck = ['WRFdens', 'WRFght', 'WRFp', 'WRFrh', 'WRFt', 'WRFtd', 'WRFws',        \
746  'WRFwd', 'WRFz']
747varNOcheckinf = "'WRFdens': air density from WRF variables; 'WRFght': geopotential"+ \
748  " height from WRF variables; 'WRFp': pressure from WRF variables; 'WRFrh': " +     \
749  "relative humidty fom WRF variables; 'WRFt': temperature from WRF variables; " +   \
750  "'WRFtd': dew-point temperature from WRF variables; 'WRFws': wind speed from " +   \
751  "WRF variables; 'WRFwd': wind speed direction from WRF variables; 'WRFz': " +     \
752  "height from WRF variables"
753
754dimshelp = "[DIM]@[simdim]@[obsdim] ',' list of couples of dimensions names from " + \
755  "each source ([DIM]='X','Y','Z','T'; None, no value)"
756vardimshelp = "[DIM]@[simvardim]@[obsvardim] ',' list of couples of variables " +    \
757  "names with dimensions values from each source ([DIM]='X','Y','Z','T'; None, " +   \
758  "no value, WRFdiagnosted variables also available: " + varNOcheckinf + ")"
759varshelp="[simvar]@[obsvar]@[[oper]@[val]] ',' list of couples of variables to " +   \
760  "validate and if necessary operation and value operations: " + opersinf +          \
761  "(WRFdiagnosted variables also available: " + varNOcheckinf + ")"
762statsn = ['minimum', 'maximum', 'mean', 'mean2', 'standard deviation']
763
764parser = OptionParser()
765parser.add_option("-d", "--dimensions", dest="dims", help=dimshelp, metavar="VALUES")
766parser.add_option("-D", "--vardimensions", dest="vardims",
767  help=vardimshelp, metavar="VALUES")
768parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs, 
769  help=strkObs, metavar="FILE")
770parser.add_option("-l", "--stationLocation", dest="stloc", 
771  help="longitude, latitude and height of the station (only for 'single-station')", 
772  metavar="FILE")
773parser.add_option("-o", "--observation", dest="fobs",
774  help="observations file to validate", metavar="FILE")
775parser.add_option("-s", "--simulation", dest="fsim",
776  help="simulation file to validate", metavar="FILE")
777parser.add_option("-t", "--trajectoryfile", dest="trajf",
778  help="file with grid points of the trajectory in the simulation grid ('simtrj')", 
779  metavar="FILE")
780parser.add_option("-v", "--variables", dest="vars",
781  help=varshelp, metavar="VALUES")
782
783(opts, args) = parser.parse_args()
784
785#######    #######
786## MAIN
787    #######
788
789ofile='validation_sim.nc'
790
791if opts.dims is None:
792    print errormsg
793    print '  ' + main + ': No list of dimensions are provided!!'
794    print '    a ',' list of values X@[dimxsim]@[dimxobs],...,T@[dimtsim]@[dimtobs]'+\
795      ' is needed'
796    quit(-1)
797else:
798    simdims = {}
799    obsdims = {}
800    print main +': couple of dimensions _______'
801    dims = {}
802    ds = opts.dims.split(',')
803    for d in ds:
804        dsecs = d.split('@')
805        if len(dsecs) != 3:
806            print errormsg
807            print '  ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
808            print '    [DIM]@[dimnsim]@[dimnobs]'
809            quit(-1)
810        dims[dsecs[0]] = [dsecs[1], dsecs[2]]
811        simdims[dsecs[0]] = dsecs[1]
812        obsdims[dsecs[0]] = dsecs[2]
813
814        print '  ',dsecs[0],':',dsecs[1],',',dsecs[2]
815       
816if opts.vardims is None:
817    print errormsg
818    print '  ' + main + ': No list of variables with dimension values are provided!!'
819    print '    a ',' list of values X@[vardimxsim]@[vardimxobs],...,T@' +  \
820      '[vardimtsim]@[vardimtobs] is needed'
821    quit(-1)
822else:
823    print main +': couple of variable dimensions _______'
824    vardims = {}
825    ds = opts.vardims.split(',')
826    for d in ds:
827        dsecs = d.split('@')
828        if len(dsecs) != 3:
829            print errormsg
830            print '  ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
831            print '    [DIM]@[vardimnsim]@[vardimnobs]'
832            quit(-1)
833        vardims[dsecs[0]] = [dsecs[1], dsecs[2]]
834        print '  ',dsecs[0],':',dsecs[1],',',dsecs[2]
835
836if opts.obskind is None:
837    print errormsg
838    print '  ' + main + ': No kind of observations provided !!'
839    quit(-1)
840else:
841    obskind = opts.obskind
842    if obskind == 'single-station':
843        if opts.stloc is None:
844            print errormsg
845            print '  ' + main + ': No station location provided !!'
846            quit(-1)
847        else:
848            stationdesc = [np.float(opts.stloc.split(',')[0]),                       \
849              np.float(opts.stloc.split(',')[1]), np.float(opts.stloc.split(',')[2])]
850
851if opts.fobs is None:
852    print errormsg
853    print '  ' + main + ': No observations file is provided!!'
854    quit(-1)
855else:
856    if not os.path.isfile(opts.fobs):
857        print errormsg
858        print '   ' + main + ": observations file '" + opts.fobs + "' does not exist !!"
859        quit(-1)
860
861if opts.fsim is None:
862    print errormsg
863    print '  ' + main + ': No simulation file is provided!!'
864    quit(-1)
865else:
866    if not os.path.isfile(opts.fsim):
867        print errormsg
868        print '   ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
869        quit(-1)
870
871if opts.vars is None:
872    print errormsg
873    print '  ' + main + ': No list of couples of variables is provided!!'
874    print '    a ',' list of values [varsim]@[varobs],... is needed'
875    quit(-1)
876else:
877    valvars = []
878    vs = opts.vars.split(',')
879    for v in vs:
880        vsecs = v.split('@')
881        if len(vsecs) < 2:
882            print errormsg
883            print '  ' + main + ': wrong number of values in:',vsecs,                \
884              ' at least 2 are needed !!'
885            print '    [varsim]@[varobs]@[[oper][val]]'
886            quit(-1)
887        if len(vsecs) > 2:
888            if not searchInlist(simopers,vsecs[2]): 
889                print errormsg
890                print main + ": operation on simulation values '" + vsecs[2] +       \
891                  "' not ready !!"
892                quit(-1)
893
894        valvars.append(vsecs)
895
896# Openning observations trajectory
897##
898oobs = NetCDFFile(opts.fobs, 'r')
899
900valdimobs = {}
901for dn in dims:
902    print dn,':',dims[dn]
903    if dims[dn][1] != 'None':
904        if not oobs.dimensions.has_key(dims[dn][1]):
905            print errormsg
906            print '  ' + main + ": observations file does not have dimension '" +    \
907              dims[dn][1] + "' !!"
908            quit(-1)
909        if vardims[dn][1] != 'None':
910            if not oobs.variables.has_key(vardims[dn][1]):
911                print errormsg
912                print '  ' + main + ": observations file does not have varibale " +  \
913                  "dimension '" + vardims[dn][1] + "' !!"
914                quit(-1)
915            valdimobs[dn] = oobs.variables[vardims[dn][1]][:]
916    else:
917        if dn == 'X':
918            valdimobs[dn] = stationdesc[0]
919        elif dn == 'Y':
920            valdimobs[dn] = stationdesc[1]
921        elif dn == 'Z':
922            valdimobs[dn] = stationdesc[2]
923
924osim = NetCDFFile(opts.fsim, 'r')
925
926valdimsim = {}
927for dn in dims:
928    if not osim.dimensions.has_key(dims[dn][0]):
929        print errormsg
930        print '  ' + main + ": simulation file does not have dimension '" +          \
931          dims[dn][0] + "' !!"
932        quit(-1)
933    if not osim.variables.has_key(vardims[dn][0]) and not                            \
934      searchInlist(varNOcheck,vardims[dn][0]):
935        print errormsg
936        print '  ' + main + ": simulation file does not have varibale dimension '" + \
937          vardims[dn][0] + "' !!"
938        quit(-1)
939    if searchInlist(varNOcheck,vardims[dn][0]):
940        valdimsim[dn] = compute_varNOcheck(osim, vardims[dn][0])
941    else:
942        valdimsim[dn] = osim.variables[vardims[dn][0]][:]
943
944# General characteristics
945dimtobs = len(valdimobs['T'])
946dimtsim = len(valdimsim['T'])
947
948print main +': observational time-steps:',dimtobs,'simulation:',dimtsim
949
950notfound = np.zeros((dimtobs), dtype=int)
951
952if obskind == 'multi-points':
953    trajpos = np.zeros((2,dimt),dtype=int)
954    for it in range(dimtobs):
955        trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'],                    \
956          [valdimobs['X'][it],valdimobss['Y'][it]])
957elif obskind == 'single-station':
958    stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],[valdimobs['Y'],             \
959      valdimobs['X']])
960    stationpos = np.zeros((2), dtype=int)
961    iid = 0
962    for idn in osim.variables[vardims['X'][0]].dimensions:
963        if idn == dims['X'][0]:
964            stationpos[1] = stsimpos[iid]
965        elif idn == dims['Y'][0]:
966            stationpos[0] = stsimpos[iid]
967
968        iid = iid + 1
969    print main + ': station point in simulation:', stationpos
970    print '    station position:',valdimobs['X'],',',valdimobs['Y']
971    print '    simulation coord.:',valdimsim['X'][tuple(stsimpos)],',',              \
972      valdimsim['Y'][tuple(stsimpos)]
973
974elif obskind == 'trajectory':
975    if opts.trajf is not None:
976        if not os.path.isfile(opts.fsim):
977            print errormsg
978            print '   ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
979            quit(-1)
980        else:
981            otrjf = NetCDFFile(opts.fsim, 'r')
982            trajpos[0,:] = otrjf.variables['obssimtrj'][0]
983            trajpos[1,:] = otrjf.variables['obssimtrj'][1]
984            otrjf.close()
985    else:
986        if dims.has_key('Z'):
987            trajpos = np.zeros((3,dimtobs),dtype=int)
988            for it in range(dimtobs):
989                if np.mod(it*100./dimtobs,10.) == 0.:
990                    print '    trajectory done: ',it*100./dimtobs,'%'
991                stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],                 \
992                  [valdimobs['Y'][it],valdimobs['X'][it]])
993                stationpos = np.zeros((2), dtype=int)
994                iid = 0
995                for idn in osim.variables[vardims['X'][0]].dimensions:
996                    if idn == dims['X'][0]:
997                        stationpos[1] = stsimpos[iid]
998                    elif idn == dims['Y'][0]:
999                        stationpos[0] = stsimpos[iid]
1000                    iid = iid + 1
1001                if stationpos[0] == 0 and stationpos[1] == 0: notfound[it] = 1
1002             
1003                trajpos[0,it] = stationpos[0]
1004                trajpos[1,it] = stationpos[1]
1005# In the simulation 'Z' varies with time ... non-hydrostatic model! ;)
1006#                trajpos[2,it] = index_mat(valdimsim['Z'][it,:,stationpos[0],         \
1007#                  stationpos[1]], valdimobs['Z'][it])
1008        else:
1009            trajpos = np.zeros((2,dimtobs),dtype=int)
1010            for it in range(dimtobs):
1011                stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],                 \
1012                  [valdimobs['Y'][it],valdimobss['X'][it]])
1013                stationpos = np.zeros((2), dtype=int)
1014                iid = 0
1015                for idn in osim.variables[vardims['X'][0]].dimensions:
1016                    if idn == dims['X'][0]:
1017                        stationpos[1] = stsimpos[iid]
1018                    elif idn == dims['Y'][0]:
1019                        stationpos[0] = stsimpos[iid]
1020                    iid = iid + 1
1021                if stationpos[0] == 0 or stationpos[1] == 0: notfound[it] = 1
1022
1023                trajpos[0,it] = stationspos[0]
1024                trajpos[1,it] = stationspos[1]
1025
1026        print main + ': not found',np.sum(notfound),'points of the trajectory'
1027
1028# Getting times
1029tobj = oobs.variables[vardims['T'][1]]
1030obstunits = tobj.getncattr('units')
1031tobj = osim.variables[vardims['T'][0]]
1032simtunits = tobj.getncattr('units')
1033
1034simobstimes = coincident_CFtimes(valdimsim['T'], obstunits, simtunits)
1035
1036# Concident times
1037##
1038coindtvalues0 = []
1039for it in range(dimtsim):   
1040    ot = 0
1041    for ito in range(ot,dimtobs-1):
1042        if valdimobs['T'][ito] < simobstimes[it] and valdimobs['T'][ito+1] >         \
1043          simobstimes[it]:
1044            ot = ito
1045            tdist = simobstimes[it] - valdimobs['T'][ito]
1046            coindtvalues0.append([it, ito, simobstimes[it], valdimobs['T'][ito],      \
1047              tdist])
1048
1049coindtvalues = np.array(coindtvalues0, dtype=np.float)
1050
1051Ncoindt = len(coindtvalues[:,0])
1052print main + ': found',Ncoindt,'coincident times between simulation and observations'
1053
1054if Ncoindt == 0:
1055    print warnmsg
1056    print main + ': no coincident times found !!'
1057    print '  stopping it'
1058    quit(-1)
1059
1060# Validating
1061##
1062
1063onewnc = NetCDFFile(ofile, 'w')
1064
1065# Dimensions
1066newdim = onewnc.createDimension('time',None)
1067newdim = onewnc.createDimension('obstime',None)
1068newdim = onewnc.createDimension('couple',2)
1069newdim = onewnc.createDimension('StrLength',StringLength)
1070newdim = onewnc.createDimension('xaround',Ngrid*2+1)
1071newdim = onewnc.createDimension('yaround',Ngrid*2+1)
1072newdim = onewnc.createDimension('stats',5)
1073
1074# Variable dimensions
1075##
1076newvar = onewnc.createVariable('obstime','f8',('time'))
1077basicvardef(newvar, 'obstime', 'time observations', obstunits )
1078set_attribute(newvar, 'calendar', 'standard')
1079newvar[:] = coindtvalues[:,3]
1080
1081newvar = onewnc.createVariable('couple', 'c', ('couple','StrLength'))
1082basicvardef(newvar, 'couple', 'couples of values', '-')
1083writing_str_nc(newvar, ['sim','obs'], StringLength)
1084
1085newvar = onewnc.createVariable('statistics', 'c', ('stats','StrLength'))
1086basicvardef(newvar, 'statistics', 'statitics from values', '-')
1087writing_str_nc(newvar, statsn, StringLength)
1088
1089if obskind == 'trajectory':
1090    if dims.has_key('Z'):
1091        newdim = onewnc.createDimension('trj',3)
1092    else:
1093        newdim = onewnc.createDimension('trj',2)
1094
1095    newvar = onewnc.createVariable('obssimtrj','i',('obstime','trj'))
1096    basicvardef(newvar, 'obssimtrj', 'trajectory on the simulation grid', '-')
1097    newvar[:] = trajpos.transpose()
1098
1099if dims.has_key('Z'):
1100    newdim = onewnc.createDimension('simtrj',4)
1101    trjsim = np.zeros((4,Ncoindt), dtype=int)
1102    trjsimval = np.zeros((4,Ncoindt), dtype=np.float)
1103else:
1104    newdim = onewnc.createDimension('simtrj',3)
1105    trjsim = np.zeros((3,Ncoindt), dtype=int)
1106    trjsimval = np.zeros((3,Ncoindt), dtype=np.float)
1107
1108Nvars = len(valvars)
1109for ivar in range(Nvars):
1110    simobsvalues = []
1111
1112    varsimobs = valvars[ivar][0] + '_' + valvars[ivar][1]
1113    print '  ' + varsimobs + '... .. .'
1114
1115    if not oobs.variables.has_key(valvars[ivar][1]):
1116        print errormsg
1117        print '  ' + main + ": observations file has not '" + valvars[ivar][1] +     \
1118          "' !!"
1119        quit(-1)
1120
1121    if not osim.variables.has_key(valvars[ivar][0]):
1122        if not searchInlist(varNOcheck, valvars[ivar][0]):
1123            print errormsg
1124            print '  ' + main + ": simulation file has not '" + valvars[ivar][0] +   \
1125              "' !!"
1126            quit(-1)
1127        else:
1128            ovsim = compute_varNOcheck(osim, valvars[ivar][0])
1129    else:
1130        ovsim = osim.variables[valvars[ivar][0]]
1131
1132    for idn in ovsim.dimensions:
1133        if not searchInlist(simdims.values(),idn):
1134            print errormsg
1135            print '  ' + main + ": dimension '" + idn + "' of variable '" +          \
1136              valvars[ivar][0] + "' not provided as reference coordinate [X,Y,Z,T] !!"
1137            quit(-1)
1138
1139    ovobs = oobs.variables[valvars[ivar][1]]
1140    if dims.has_key('Z'):
1141        simobsSvalues = np.zeros((Ncoindt, Ngrid*2+1, Ngrid*2+1, Ngrid*2+1),         \
1142          dtype = np.float)
1143    else:
1144        simobsSvalues = np.zeros((Ncoindt, Ngrid*2+1, Ngrid*2+1), dtype = np.float)
1145
1146    if obskind == 'multi-points':
1147        for it in range(Ncoindt):
1148            slicev = dims['X'][0] + ':' + str(trajpos[2,it]) + '|' +                 \
1149              dims['Y'][0]+ ':' + str(trajpos[1,it]) + '|' +                         \
1150              dims['T'][0]+ ':' + str(coindtvalues[it][0])
1151            slicevar, dimslice = slice_variable(ovsim, slicev)
1152            simobsvalues.append([ slicevar, ovobs[coindtvalues[it][1]]])
1153            slicev = dims['X'][0] + ':' + str(trajpos[2,it]-Ngrid) + '@' +           \
1154              str(trajpos[2,it]+Ngrid) + '|' + dims['Y'][0] + ':' +                  \
1155              str(trajpos[1,it]-Ngrid) + '@' + str(trajpos[1,it]+Ngrid) + '|' +      \
1156              dims['T'][0]+':'+str(coindtvalues[it][0])
1157            slicevar, dimslice = slice_variable(ovsim, slicev)
1158            simobsSvalues[it,:,:] = slicevar
1159
1160    elif obskind == 'single-station':
1161        for it in range(Ncoindt):
1162            slicev = dims['X'][0]+':'+str(stationpos[1]) + '|' +                      \
1163              dims['Y'][0]+':'+str(stationpos[0]) + '|' +                             \
1164              dims['T'][0]+':'+str(coindtvalues[it][0])
1165            slicevar, dimslice = slice_variable(ovsim, slicev)
1166            slicev = dims['X'][0] + ':' + str(trajpos[2,it]-Ngrid) + '@' +           \
1167              str(trajpos[2,it]+Ngrid) + '|' + dims['Y'][0] + ':' +                  \
1168              str(trajpos[1,it]-Ngrid) + '@' + str(trajpos[1,it]+Ngrid) + '|' +      \
1169              dims['T'][0] + ':' + str(coindtvalues[it,0])
1170            slicevar, dimslice = slice_variable(ovsim, slicev)
1171            simobsSvalues[it,:,:] = slicevar
1172    elif obskind == 'trajectory':
1173        if dims.has_key('Z'):
1174            for it in range(Ncoindt):
1175                ito = int(coindtvalues[it,1])
1176                if notfound[ito] == 0:
1177                    trajpos[2,ito] = index_mat(valdimsim['Z'][coindtvalues[it,0],:,  \
1178                      trajpos[1,ito],trajpos[0,ito]], valdimobs['Z'][ito])
1179                    slicev = dims['X'][0]+':'+str(trajpos[0,ito]) + '|' +            \
1180                      dims['Y'][0]+':'+str(trajpos[1,ito]) + '|' +                   \
1181                      dims['Z'][0]+':'+str(trajpos[2,ito]) + '|' +                   \
1182                      dims['T'][0]+':'+str(int(coindtvalues[it,0]))
1183                    slicevar, dimslice = slice_variable(ovsim, slicev)
1184                    simobsvalues.append([ slicevar, ovobs[int(ito)]])
1185                    minx = np.max([trajpos[0,ito]-Ngrid,0])
1186                    maxx = np.min([trajpos[0,ito]+Ngrid+1,ovsim.shape[3]])
1187                    miny = np.max([trajpos[1,ito]-Ngrid,0])
1188                    maxy = np.min([trajpos[1,ito]+Ngrid+1,ovsim.shape[2]])
1189                    minz = np.max([trajpos[2,ito]-Ngrid,0])
1190                    maxz = np.min([trajpos[2,ito]+Ngrid+1,ovsim.shape[1]])
1191
1192                    slicev = dims['X'][0] + ':' + str(minx) + '@' + str(maxx) + '|' +\
1193                      dims['Y'][0] + ':' + str(miny) + '@' + str(maxy) + '|' +       \
1194                      dims['Z'][0] + ':' + str(minz) + '@' + str(maxz) + '|' +       \
1195                      dims['T'][0] + ':' + str(int(coindtvalues[it,0]))
1196                    slicevar, dimslice = slice_variable(ovsim, slicev)
1197
1198                    sliceS = []
1199                    sliceS.append(it)
1200                    sliceS.append(slice(0,maxz-minz))
1201                    sliceS.append(slice(0,maxy-miny))
1202                    sliceS.append(slice(0,maxx-minx))
1203
1204                    simobsSvalues[tuple(sliceS)] = slicevar
1205                    if ivar == 0:
1206                        trjsim[0,it] = trajpos[0,ito]
1207                        trjsim[1,it] = trajpos[1,ito]
1208                        trjsim[2,it] = trajpos[2,ito]
1209                        trjsim[3,it] = coindtvalues[it,0]
1210                else:
1211                    simobsvalues.append([fillValueF, fillValueF])
1212                    simobsSvalues[it,:,:,:]= np.ones((Ngrid*2+1,Ngrid*2+1,Ngrid*2+1),\
1213                      dtype = np.float)*fillValueF
1214        else:
1215            for it in range(Ncoindt):
1216                if notfound[it] == 0:
1217                    ito = coindtvalues[it,1]
1218                    slicev = dims['X'][0]+':'+str(trajpos[2,ito]) + '|' +            \
1219                      dims['Y'][0]+':'+str(trajpos[1,ito]) + '|' +                   \
1220                      dims['T'][0]+':'+str(coindtvalues[ito,0])
1221                    slicevar, dimslice = slice_variable(ovsim, slicev)
1222                    simobsvalues.append([ slicevar, ovobs[coindtvalues[it,1]]])
1223                    slicev = dims['X'][0] + ':' + str(trajpos[0,it]-Ngrid) + '@' +   \
1224                      str(trajpos[0,it]+Ngrid) + '|' + dims['Y'][0] + ':' +          \
1225                      str(trajpos[1,it]-Ngrid) + '@' + str(trajpos[1,it]+Ngrid) +    \
1226                      '|' + dims['T'][0] + ':' + str(coindtvalues[it,0])
1227                    slicevar, dimslice = slice_variable(ovsim, slicev)
1228                    simobsSvalues[it,:,:] = slicevar
1229                else:
1230                    simobsvalues.append([fillValue, fillValue])
1231                    simobsSvalues[it,:,:] = np.ones((Ngrid*2+1,Ngrid*2+1),           \
1232                      dtype = np.float)*fillValueF
1233                print simobsvalues[varsimobs][:][it]
1234
1235    arrayvals = np.array(simobsvalues)
1236    if len(valvars[ivar]) > 2:
1237        const=np.float(valvars[ivar][3])
1238        if valvars[ivar][2] == 'sumc':
1239            simobsSvalues = simobsSvalues + const
1240            arrayvals[:,0] = arrayvals[:,0] + const
1241        elif valvars[ivar][2] == 'subc':
1242            simobsSvalues = simobsSvalues - const
1243            arrayvals[:,0] = arrayvals[:,0] - const
1244        elif valvars[ivar][2] == 'mulc':
1245            simobsSvalues = simobsSvalues * const
1246            arrayvals[:,0] = arrayvals[:,0] * const
1247        elif valvars[ivar][2] == 'divc':
1248            simobsSvalues = simobsSvalues / const
1249            arrayvals[:,0] = arrayvals[:,0] / const
1250        else:
1251            print errormsg
1252            print '  ' + fname + ": operation '" + valvars[ivar][2] + "' not ready!!"
1253            quit(-1)
1254
1255# Statistics around values
1256    aroundstats = np.zeros((5,Ncoindt), dtype=np.float)
1257    for it in range(Ncoindt):
1258        aroundstats[0,it] = np.min(simobsSvalues[it,])
1259        aroundstats[1,it] = np.max(simobsSvalues[it,])
1260        aroundstats[2,it] = np.mean(simobsSvalues[it,])
1261        aroundstats[3,it] = np.mean(simobsSvalues[it,]*simobsSvalues[it,])
1262        aroundstats[4,it] = np.sqrt(aroundstats[3,it] - aroundstats[2,it]*           \
1263          aroundstats[2,it])
1264
1265# Values to netCDF
1266    newvar = onewnc.createVariable(varsimobs, 'f', ('time','couple'),                \
1267      fill_value=fillValueF)
1268    descvar = 'couples of simulated: ' + valvars[ivar][0] + ' and observed ' +       \
1269      valvars[ivar][1]
1270    basicvardef(newvar, varsimobs, descvar, ovobs.getncattr('units'))
1271    newvar[:] = arrayvals
1272
1273# Around values
1274    if dims.has_key('Z'):
1275        if not onewnc.dimensions.has_key('zaround'):
1276            newdim = onewnc.createDimension('zaround',Ngrid*2+1)
1277
1278        newvar = onewnc.createVariable(valvars[ivar][0] + 'around', 'f',            \
1279          ('time','zaround','yaround','xaround'), fill_value=fillValueF)
1280    else:
1281        newvar = onewnc.createVariable(valvars[ivar][0] + 'around', 'f',            \
1282          ('time','yaround','xaround'), fill_value=fillValueF)
1283
1284    descvar = 'around simulated values +/- grid values: ' + valvars[ivar][0]
1285    basicvardef(newvar, varsimobs + 'around', descvar, ovobs.getncattr('units'))
1286    newvar[:] = simobsSvalues
1287
1288# Statistics
1289    newvar = onewnc.createVariable(valvars[ivar][0] + 'staround', 'f',               \
1290      ('time','stats'), fill_value=fillValueF)
1291    descvar = 'around simulated statisitcs: ' + valvars[ivar][0]
1292    basicvardef(newvar, varsimobs + 'staround', descvar, ovobs.getncattr('units'))
1293    newvar[:] = aroundstats.transpose()
1294
1295    onewnc.sync()
1296
1297newvar = onewnc.createVariable('simtrj','i',('time','simtrj'))
1298basicvardef(newvar, 'simtrj', 'coordinates [X,Y,Z,T] of the coincident trajectory ' +\
1299  'in sim', obstunits)
1300newvar[:] = trjsim.transpose()
1301
1302# Global attributes
1303##
1304set_attribute(onewnc,'author_nc','Lluis Fita')
1305set_attribute(onewnc,'institution_nc','Laboratoire de Meteorology Dynamique, ' +    \
1306  'LMD-Jussieu, UPMC, Paris')
1307set_attribute(onewnc,'country_nc','France')
1308set_attribute(onewnc,'script_nc',main)
1309set_attribute(onewnc,'version_script',version)
1310set_attribute(onewnc,'information',                                                 \
1311  'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html')
1312set_attribute(onewnc,'simfile',opts.fsim)
1313set_attribute(onewnc,'obsfile',opts.fobs)
1314
1315onewnc.sync()
1316onewnc.close()
1317
1318print main + ": successfull writting of '" + ofile + "' !!"
Note: See TracBrowser for help on using the repository browser.