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

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

Working and tested version for 'single-station'
File size: 27.7 KB
Line 
1
2# L. Fita, LMD-Jussieu. February 2015
3## e.g. # 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
4import numpy as np
5import os
6import re
7from optparse import OptionParser
8from netCDF4 import Dataset as NetCDFFile
9
10main = 'validarion_sim.py'
11errormsg = 'ERROR -- errror -- ERROR -- error'
12warnmsg = 'WARNING -- warning -- WARNING -- warning'
13
14# version
15version=1.0
16
17# Filling values for floats, integer and string
18fillValueF = 1.e20
19fillValueI = -99999
20fillValueS = '---'
21
22StringLength = 50
23
24def searchInlist(listname, nameFind):
25    """ Function to search a value within a list
26    listname = list
27    nameFind = value to find
28    >>> searInlist(['1', '2', '3', '5'], '5')
29    True
30    """
31    for x in listname:
32      if x == nameFind:
33        return True
34    return False
35
36def set_attribute(ncvar, attrname, attrvalue):
37    """ Sets a value of an attribute of a netCDF variable. Removes previous attribute value if exists
38    ncvar = object netcdf variable
39    attrname = name of the attribute
40    attrvalue = value of the attribute
41    """
42    import numpy as np
43    from netCDF4 import Dataset as NetCDFFile
44
45    attvar = ncvar.ncattrs()
46    if searchInlist(attvar, attrname):
47        attr = ncvar.delncattr(attrname)
48
49    attr = ncvar.setncattr(attrname, attrvalue)
50
51    return ncvar
52
53def basicvardef(varobj, vstname, vlname, vunits):
54    """ Function to give the basic attributes to a variable
55    varobj= netCDF variable object
56    vstname= standard name of the variable
57    vlname= long name of the variable
58    vunits= units of the variable
59    """
60    attr = varobj.setncattr('standard_name', vstname)
61    attr = varobj.setncattr('long_name', vlname)
62    attr = varobj.setncattr('units', vunits)
63
64    return
65
66def writing_str_nc(varo, values, Lchar):
67    """ Function to write string values in a netCDF variable as a chain of 1char values
68    varo= netCDF variable object
69    values = list of values to introduce
70    Lchar = length of the string in the netCDF file
71    """
72
73    Nvals = len(values)
74    for iv in range(Nvals):   
75        stringv=values[iv] 
76        charvals = np.chararray(Lchar)
77        Lstr = len(stringv)
78        charvals[Lstr:Lchar] = ''
79
80        for ich in range(Lstr):
81            charvals[ich] = stringv[ich:ich+1]
82
83        varo[iv,:] = charvals
84
85    return
86
87def index_3mat(matA,matB,matC,val):
88    """ Function to provide the coordinates of a given value inside three matrix simultaneously
89    index_mat(matA,matB,matC,val)
90      matA= matrix with one set of values
91      matB= matrix with the other set of values
92      matB= matrix with the third set of values
93      val= triplet of values to search
94    >>> 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])
95    [2 1 1]
96    """
97    fname = 'index_3mat'
98
99    matAshape = matA.shape
100    matBshape = matB.shape
101    matCshape = matC.shape
102
103    for idv in range(len(matAshape)):
104        if matAshape[idv] != matBshape[idv]:
105            print errormsg
106            print '  ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv],  \
107              'and B:',matBshape[idv],'does not coincide!!'
108            quit(-1)
109        if matAshape[idv] != matCshape[idv]:
110            print errormsg
111            print '  ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv],  \
112              'and C:',matCshape[idv],'does not coincide!!'
113            quit(-1)
114
115    minA = np.min(matA)
116    maxA = np.max(matA)
117    minB = np.min(matB)
118    maxB = np.max(matB)
119    minC = np.min(matC)
120    maxC = np.max(matC)
121
122    if val[0] < minA or val[0] > maxA:
123        print warnmsg
124        print '  ' + fname + ': first value:',val[0],'outside matA range',minA,',',  \
125          maxA,'!!'
126    if val[1] < minB or val[1] > maxB:
127        print warnmsg
128        print '  ' + fname + ': second value:',val[1],'outside matB range',minB,',',  \
129          maxB,'!!'
130    if val[2] < minC or val[2] > maxC:
131        print warnmsg
132        print '  ' + fname + ': second value:',val[2],'outside matC range',minC,',',  \
133          maxC,'!!'
134
135    dist = np.zeros(tuple(matAshape), dtype=np.float)
136    dist = np.sqrt((matA - np.float(val[0]))**2 + (matB - np.float(val[1]))**2 +     \
137      (matC - np.float(val[2]))**2)
138
139    mindist = np.min(dist)
140   
141    matlist = list(dist.flatten())
142    ifound = matlist.index(mindist)
143
144    Ndims = len(matAshape)
145    valpos = np.zeros((Ndims), dtype=int)
146    baseprevdims = np.zeros((Ndims), dtype=int)
147
148    for dimid in range(Ndims):
149        baseprevdims[dimid] = np.product(matAshape[dimid+1:Ndims])
150        if dimid == 0:
151            alreadyplaced = 0
152        else:
153            alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
154        valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
155
156    return valpos
157
158def index_2mat(matA,matB,val):
159    """ Function to provide the coordinates of a given value inside two matrix simultaneously
160    index_mat(matA,matB,val)
161      matA= matrix with one set of values
162      matB= matrix with the pother set of values
163      val= couple of values to search
164    >>> index_mat(np.arange(27).reshape(3,3,3),np.arange(100,127).reshape(3,3,3),[22,111])
165    [2 1 1]
166    """
167    fname = 'index_2mat'
168
169    matAshape = matA.shape
170    matBshape = matB.shape
171
172    for idv in range(len(matAshape)):
173        if matAshape[idv] != matBshape[idv]:
174            print errormsg
175            print '  ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv],  \
176              'and B:',matBshape[idv],'does not coincide!!'
177            quit(-1)
178
179    minA = np.min(matA)
180    maxA = np.max(matA)
181    minB = np.min(matB)
182    maxB = np.max(matB)
183
184    if val[0] < minA or val[0] > maxA:
185        print warnmsg
186        print '  ' + fname + ': first value:',val[0],'outside matA range',minA,',',  \
187          maxA,'!!'
188    if val[1] < minB or val[1] > maxB:
189        print warnmsg
190        print '  ' + fname + ': second value:',val[1],'outside matB range',minB,',',  \
191          maxB,'!!'
192
193    dist = np.zeros(tuple(matAshape), dtype=np.float)
194    dist = np.sqrt((matA - np.float(val[0]))**2 + (matB - np.float(val[1]))**2)
195
196    mindist = np.min(dist)
197   
198    matlist = list(dist.flatten())
199    ifound = matlist.index(mindist)
200
201    Ndims = len(matAshape)
202    valpos = np.zeros((Ndims), dtype=int)
203    baseprevdims = np.zeros((Ndims), dtype=int)
204
205    for dimid in range(Ndims):
206        baseprevdims[dimid] = np.product(matAshape[dimid+1:Ndims])
207        if dimid == 0:
208            alreadyplaced = 0
209        else:
210            alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
211        valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
212
213    return valpos
214
215def index_mat(mat,val):
216    """ Function to provide the coordinates of a given value inside a matrix
217    index_mat(mat,val)
218      mat= matrix with values
219      val= value to search
220    >>> index_mat(np.arange(27).reshape(3,3,3),22)
221    [2 1 1]
222    """
223
224    fname = 'index_mat'
225
226    matshape = mat.shape
227
228    matlist = list(mat.flatten())
229    ifound = matlist.index(val)
230
231    Ndims = len(matshape)
232    valpos = np.zeros((Ndims), dtype=int)
233    baseprevdims = np.zeros((Ndims), dtype=int)
234
235    for dimid in range(Ndims):
236        baseprevdims[dimid] = np.product(matshape[dimid+1:Ndims])
237        if dimid == 0:
238            alreadyplaced = 0
239        else:
240            alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
241        valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
242
243    return valpos
244
245def coincident_CFtimes(tvalB, tunitA, tunitB):
246    """ Function to make coincident times for two different sets of CFtimes
247    tvalB= time values B
248    tunitA= time units times A to which we want to make coincidence
249    tunitB= time units times B
250    >>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00',
251      'hours since 1949-12-01 00:00:00')
252    [     0.   3600.   7200.  10800.  14400.  18000.  21600.  25200.  28800.  32400.]
253    >>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00',
254      'hours since 1979-12-01 00:00:00')
255    [  9.46684800e+08   9.46688400e+08   9.46692000e+08   9.46695600e+08
256       9.46699200e+08   9.46702800e+08   9.46706400e+08   9.46710000e+08
257       9.46713600e+08   9.46717200e+08]
258    """
259    import datetime as dt
260    fname = 'coincident_CFtimes'
261
262    trefA = tunitA.split(' ')[2] + ' ' + tunitA.split(' ')[3]
263    trefB = tunitB.split(' ')[2] + ' ' + tunitB.split(' ')[3]
264    tuA = tunitA.split(' ')[0]
265    tuB = tunitB.split(' ')[0]
266
267    if tuA != tuB:
268        if tuA == 'microseconds':
269            if tuB == 'microseconds':
270                tB = tvalB*1.
271            elif tuB == 'seconds':
272                tB = tvalB*10.e6
273            elif tuB == 'minutes':
274                tB = tvalB*60.*10.e6
275            elif tuB == 'hours':
276                tB = tvalB*3600.*10.e6
277            elif tuB == 'days':
278                tB = tvalB*3600.*24.*10.e6
279            else:
280                print errormsg
281                print '  ' + fname + ": combination of time untis: '" + tuA +        \
282                  "' & '" + tuB + "' not ready !!"
283                quit(-1)
284        elif tuA == 'seconds':
285            if tuB == 'microseconds':
286                tB = tvalB/10.e6
287            elif tuB == 'seconds':
288                tB = tvalB*1.
289            elif tuB == 'minutes':
290                tB = tvalB*60.
291            elif tuB == 'hours':
292                tB = tvalB*3600.
293            elif tuB == 'days':
294                tB = tvalB*3600.*24.
295            else:
296                print errormsg
297                print '  ' + fname + ": combination of time untis: '" + tuA +        \
298                  "' & '" + tuB + "' not ready !!"
299                quit(-1)
300        elif tuA == 'minutes':
301            if tuB == 'microseconds':
302                tB = tvalB/(60.*10.e6)
303            elif tuB == 'seconds':
304                tB = tvalB/60.
305            elif tuB == 'minutes':
306                tB = tvalB*1.
307            elif tuB == 'hours':
308                tB = tvalB*60.
309            elif tuB == 'days':
310                tB = tvalB*60.*24.
311            else:
312                print errormsg
313                print '  ' + fname + ": combination of time untis: '" + tuA +        \
314                  "' & '" + tuB + "' not ready !!"
315                quit(-1)
316        elif tuA == 'hours':
317            if tuB == 'microseconds':
318                tB = tvalB/(3600.*10.e6)
319            elif tuB == 'seconds':
320                tB = tvalB/3600.
321            elif tuB == 'minutes':
322                tB = tvalB/60.
323            elif tuB == 'hours':
324                tB = tvalB*1.
325            elif tuB == 'days':
326                tB = tvalB*24.
327            else:
328                print errormsg
329                print '  ' + fname + ": combination of time untis: '" + tuA +        \
330                  "' & '" + tuB + "' not ready !!"
331                quit(-1)
332        elif tuA == 'days':
333            if tuB == 'microseconds':
334                tB = tvalB/(24.*3600.*10.e6)
335            elif tuB == 'seconds':
336                tB = tvalB/(24.*3600.)
337            elif tuB == 'minutes':
338                tB = tvalB/(24.*60.)
339            elif tuB == 'hours':
340                tB = tvalB/24.
341            elif tuB == 'days':
342                tB = tvalB*1.
343            else:
344                print errormsg
345                print '  ' + fname + ": combination of time untis: '" + tuA +        \
346                  "' & '" + tuB + "' not ready !!"
347                quit(-1)
348        else:
349            print errormsg
350            print '  ' + fname + ": time untis: '" + tuA + "' not ready !!"
351            quit(-1)
352    else:
353        tB = tvalB*1.
354
355    if trefA != trefB:
356        trefTA = dt.datetime.strptime(trefA, '%Y-%m-%d %H:%M:%S')
357        trefTB = dt.datetime.strptime(trefB, '%Y-%m-%d %H:%M:%S')
358
359        difft = trefTB - trefTA
360        diffv = difft.days*24.*3600.*10.e6 + difft.seconds*10.e6 + difft.microseconds
361        print '  ' + fname + ': different reference refA:',trefTA,'refB',trefTB
362        print '    difference:',difft,':',diffv,'microseconds'
363
364        if tuA == 'microseconds':
365            tB = tB + diffv
366        elif tuA == 'seconds':
367            tB = tB + diffv/10.e6
368        elif tuA == 'minutes':
369            tB = tB + diffv/(60.*10.e6)
370        elif tuA == 'hours':
371            tB = tB + diffv/(3600.*10.e6)
372        elif tuA == 'dayss':
373            tB = tB + diffv/(24.*3600.*10.e6)
374        else:
375            print errormsg
376            print '  ' + fname + ": time untis: '" + tuA + "' not ready !!"
377            quit(-1)
378
379    return tB
380
381
382def slice_variable(varobj, dimslice):
383    """ Function to return a slice of a given variable according to values to its
384      dimensions
385    slice_variable(varobj, dimslice)
386      varobj= object wit the variable
387      dimslice= [[dimname1]:[value1]|[[dimname2]:[value2], ...] pairs of dimension
388        [value]:
389          * [integer]: which value of the dimension
390          * -1: all along the dimension
391          * -9: last value of the dimension
392          * [beg]@[end] slice from [beg] to [end]
393    """
394    fname = 'slice_variable'
395
396    if varobj == 'h':
397        print fname + '_____________________________________________________________'
398        print slice_variable.__doc__
399        quit()
400
401    vardims = varobj.dimensions
402    Ndimvar = len(vardims)
403
404    Ndimcut = len(dimslice.split('|'))
405    dimsl = dimslice.split('|')
406
407    varvalsdim = []
408    dimnslice = []
409
410    for idd in range(Ndimvar):
411        for idc in range(Ndimcut):
412            dimcutn = dimsl[idc].split(':')[0]
413            dimcutv = dimsl[idc].split(':')[1]
414            if vardims[idd] == dimcutn: 
415                posfrac = dimcutv.find('@')
416                if posfrac != -1:
417                    inifrac = int(dimcutv.split('@')[0])
418                    endfrac = int(dimcutv.split('@')[1])
419                    varvalsdim.append(slice(inifrac,endfrac))
420                    dimnslice.append(vardims[idd])
421                else:
422                    if int(dimcutv) == -1:
423                        varvalsdim.append(slice(0,varobj.shape[idd]))
424                        dimnslice.append(vardims[idd])
425                    elif int(dimcutv) == -9:
426                        varvalsdim.append(int(varobj.shape[idd])-1)
427                    else:
428                        varvalsdim.append(int(dimcutv))
429                break
430
431    varvalues = varobj[tuple(varvalsdim)]
432
433    return varvalues, dimnslice
434
435
436####### ###### ##### #### ### ## #
437
438strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " +  \
439  " and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')"
440
441kindobs=['multi-points', 'single-station', 'trajectory']
442strkObs="kind of observations; 'multi-points': multiple individual punctual obs " +  \
443  "(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
444  "'trajectory': following a trajectory"
445simopers =['sumc','subc','mulc','divc']
446#sumc,[constant]: add [constant] to variables values
447#subc,[constant]: substract [constant] to variables values
448#mulc,[constant]: multipy by [constant] to variables values
449#divc,[constant]: divide by [constant] to variables values
450
451parser = OptionParser()
452parser.add_option("-d", "--dimensions", dest="dims",
453  help="[DIM]@[simdim]@[obsdim] ',' list of couples of dimensions names from each source ([DIM]='X','Y','Z','T'; None, no value)",
454  metavar="VALUES")
455parser.add_option("-D", "--vardimensions", dest="vardims",
456  help="[DIM]@[simvardim]@[obsvardim] ',' list of couples of variables names with dimensions values from each source ([DIM]='X','Y','Z','T'; None, no value)", metavar="VALUES")
457parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs, 
458  help=strkObs, metavar="FILE")
459parser.add_option("-l", "--stationLocation", dest="stloc", 
460  help="longitude, latitude and height of the station (only for 'single-station')", 
461  metavar="FILE")
462parser.add_option("-o", "--observation", dest="fobs",
463  help="observations file to validate", metavar="FILE")
464parser.add_option("-s", "--simulation", dest="fsim",
465  help="simulation file to validate", metavar="FILE")
466parser.add_option("-v", "--variables", dest="vars",
467  help="[simvar]@[obsvar]@[[oper]@[val]] ',' list of couples of variables to validate and if necessary operation and value", metavar="VALUES")
468
469(opts, args) = parser.parse_args()
470
471#######    #######
472## MAIN
473    #######
474
475ofile='validation_sim.nc'
476
477if opts.dims is None:
478    print errormsg
479    print '  ' + main + ': No list of dimensions are provided!!'
480    print '    a ',' list of values X@[dimxsim]@[dimxobs],...,T@[dimtsim]@[dimtobs]'+\
481      ' is needed'
482    quit(-1)
483else:
484    print main +': couple of dimensions _______'
485    dims = {}
486    ds = opts.dims.split(',')
487    for d in ds:
488        dsecs = d.split('@')
489        if len(dsecs) != 3:
490            print errormsg
491            print '  ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
492            print '    [DIM]@[dimnsim]@[dimnobs]'
493            quit(-1)
494        dims[dsecs[0]] = [dsecs[1], dsecs[2]]
495        print '  ',dsecs[0],':',dsecs[1],',',dsecs[2]
496       
497if opts.vardims is None:
498    print errormsg
499    print '  ' + main + ': No list of variables with dimension values are provided!!'
500    print '    a ',' list of values X@[vardimxsim]@[vardimxobs],...,T@' +  \
501      '[vardimtsim]@[vardimtobs] is needed'
502    quit(-1)
503else:
504    print main +': couple of variable dimensions _______'
505    vardims = {}
506    ds = opts.vardims.split(',')
507    for d in ds:
508        dsecs = d.split('@')
509        if len(dsecs) != 3:
510            print errormsg
511            print '  ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
512            print '    [DIM]@[vardimnsim]@[vardimnobs]'
513            quit(-1)
514        vardims[dsecs[0]] = [dsecs[1], dsecs[2]]
515        print '  ',dsecs[0],':',dsecs[1],',',dsecs[2]
516
517if opts.obskind is None:
518    print errormsg
519    print '  ' + main + ': No kind of observations provided !!'
520    quit(-1)
521else:
522    obskind = opts.obskind
523    if obskind == 'single-station':
524        if opts.stloc is None:
525            print errormsg
526            print '  ' + main + ': No station location provided !!'
527            quit(-1)
528        else:
529            stationdesc = [np.float(opts.stloc.split(',')[0]),                       \
530              np.float(opts.stloc.split(',')[1]), np.float(opts.stloc.split(',')[2])]
531
532if opts.fobs is None:
533    print errormsg
534    print '  ' + main + ': No observations file is provided!!'
535    quit(-1)
536else:
537    if not os.path.isfile(opts.fobs):
538        print errormsg
539        print '   ' + main + ": observations file '" + opts.fobs + "' does not exist !!"
540        quit(-1)
541
542if opts.fsim is None:
543    print errormsg
544    print '  ' + main + ': No simulation file is provided!!'
545    quit(-1)
546else:
547    if not os.path.isfile(opts.fsim):
548        print errormsg
549        print '   ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
550        quit(-1)
551
552if opts.vars is None:
553    print errormsg
554    print '  ' + main + ': No list of couples of variables is provided!!'
555    print '    a ',' list of values [varsim]@[varobs],... is needed'
556    quit(-1)
557else:
558    valvars = []
559    vs = opts.vars.split(',')
560    for v in vs:
561        vsecs = v.split('@')
562        if len(vsecs) < 2:
563            print errormsg
564            print '  ' + main + ': wrong number of values in:',vsecs,                \
565              ' at least 2 are needed !!'
566            print '    [varsim]@[varobs]@[[oper][val]]'
567            quit(-1)
568        if len(vsecs) > 2:
569            if not searchInlist(simopers,vsecs[2]): 
570                print errormsg
571                print main + ": operation on simulation values '" + vsecs[2] +       \
572                  "' not ready !!"
573                quit(-1)
574
575        valvars.append(vsecs)
576
577# Openning observations trajectory
578##
579oobs = NetCDFFile(opts.fobs, 'r')
580
581valdimobs = {}
582for dn in dims:
583    print dn,':',dims[dn]
584    if dims[dn][1] != 'None':
585        if not oobs.dimensions.has_key(dims[dn][1]):
586            print errormsg
587            print '  ' + main + ": observations file does not have dimension '" +    \
588              dims[dn][1] + "' !!"
589            quit(-1)
590        if vardims[dn][1] != 'None':
591            if not oobs.variables.has_key(vardims[dn][1]):
592                print errormsg
593                print '  ' + main + ": observations file does not have varibale " +  \
594                  "dimension '" + vardims[dn][1] + "' !!"
595                quit(-1)
596            valdimobs[dn] = oobs.variables[vardims[dn][1]][:]
597    else:
598        if dn == 'X':
599            valdimobs[dn] = stationdesc[0]
600        elif dn == 'Y':
601            valdimobs[dn] = stationdesc[1]
602        elif dn == 'Z':
603            valdimobs[dn] = stationdesc[2]
604
605osim = NetCDFFile(opts.fsim, 'r')
606
607valdimsim = {}
608for dn in dims:
609    if not osim.dimensions.has_key(dims[dn][0]):
610        print errormsg
611        print '  ' + main + ": simulation file does not have dimension '" +          \
612          dims[dn][0] + "' !!"
613        quit(-1)
614    if not osim.variables.has_key(vardims[dn][0]):
615        print errormsg
616        print '  ' + main + ": simulation file does not have varibale dimension '" + \
617          vardims[dn][0] + "' !!"
618        quit(-1)
619    valdimsim[dn] = osim.variables[vardims[dn][0]][:]
620
621# General characteristics
622dimtobs = len(valdimobs['T'])
623dimtsim = len(valdimsim['T'])
624
625print main +': observational time-steps:',dimtobs,'simulation:',dimtsim
626
627if obskind == 'multi-points':
628    trajpos = np.zeros((2,dimt),dtype=int)
629    for it in dimtobs:
630        trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'],                    \
631          [valdimobs['X'][it],valdimobss['Y'][it]])
632elif obskind == 'single-station':
633    stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],[valdimobs['Y'],             \
634      valdimobs['X']])
635    stationpos = np.zeros((2), dtype=int)
636    iid = 0
637    for idn in osim.variables[vardims['X'][0]].dimensions:
638        if idn == dims['X'][0]:
639            stationpos[1] = stsimpos[iid]
640        elif idn == dims['Y'][0]:
641            stationpos[0] = stsimpos[iid]
642
643        iid = iid + 1
644    print main + ': station point in simulation:', stationpos
645    print '    station position:',valdimobs['X'],',',valdimobs['Y']
646    print '    simulation coord.:',valdimsim['X'][tuple(stsimpos)],',',              \
647      valdimsim['Y'][tuple(stsimpos)]
648
649elif obskind == 'trajectory':
650    if dims.has_key('Z'):
651        trajpos = np.zeros((3,dimt),dtype=int)
652        for it in dimtobs:
653            trajpos[0:1,it] = index_2mat(valdimsim['Y'],valdimsim['X'],              \
654              [valdimobs['Y'][it],valdimobss['X'][it]])
655            trajpos[2,it] = index_mat(valdimsim['Z'],valdimobs['Z'][it])
656    else:
657        trajpos = np.zeros((2,dimt),dtype=int)
658        for it in dimtobs:
659            trajpos[:,it] = index_2mat(valdimsim['Y'],valdimsim['X'],                \
660              [valdimobs['Y'][it],valdimobss['X'][it]])
661
662# Getting times
663tobj = oobs.variables[vardims['T'][1]]
664obstunits = tobj.getncattr('units')
665tobj = osim.variables[vardims['T'][0]]
666simtunits = tobj.getncattr('units')
667
668simobstimes = coincident_CFtimes(valdimsim['T'], obstunits, simtunits)
669
670# Concident times
671##
672coindtvalues = []
673for it in range(dimtsim):   
674    ot = 0
675    for ito in range(ot,dimtobs-1):
676        if valdimobs['T'][ito] < simobstimes[it] and valdimobs['T'][ito+1] >         \
677          simobstimes[it]:
678            ot = ito
679            tdist = simobstimes[it] - valdimobs['T'][ito]
680            coindtvalues.append([it,ito,simobstimes[it],valdimobs['T'][ito],tdist])
681
682Ncoindt = len(coindtvalues)
683print main + ': found',Ncoindt,'coincident times between simulation and observations'
684
685# Validating
686##
687
688onewnc = NetCDFFile(ofile, 'w')
689
690# Dimensions
691newdim = onewnc.createDimension('time',None)
692newdim = onewnc.createDimension('couple',2)
693newdim = onewnc.createDimension('StrLength',StringLength)
694
695# Variable dimensions
696##
697newvar = onewnc.createVariable('simtime','f8',('time'))
698basicvardef(newvar, 'simtime', 'time simulation', obstunits )
699set_attribute(newvar, 'calendar', 'standard')
700newvar[:] = coindtvalues[:][2]
701
702newvar = onewnc.createVariable('obstime','f8',('time'))
703basicvardef(newvar, 'obstime', 'time observations', obstunits )
704set_attribute(newvar, 'calendar', 'standard')
705newvar[:] = coindtvalues[:][3]
706
707newvar = onewnc.createVariable('couple', 'c', ('couple','StrLength'))
708basicvardef(newvar, 'couple', 'couples of values', '-')
709writing_str_nc(newvar, ['sim','obs'], StringLength)
710
711Nvars = len(valvars)
712for ivar in range(Nvars):
713    simobsvalues = []
714
715    varsimobs = valvars[ivar][0] + '_' + valvars[ivar][1]
716
717    if not oobs.variables.has_key(valvars[ivar][1]):
718        print errormsg
719        print '  ' + main + ": observations file has not '" + valvars[ivar][1] +     \
720          "' !!"
721        quit(-1)
722    if not osim.variables.has_key(valvars[ivar][0]):
723        print errormsg
724        print '  ' + main + ": simulation file has not '" + valvars[ivar][0] + "' !!"
725        quit(-1)
726
727    ovobs = oobs.variables[valvars[ivar][1]]
728    ovsim = osim.variables[valvars[ivar][0]]
729
730    if obskind == 'multi-points':
731        for it in range(Ncoindt):
732            slicev = dims['X'][0]+':'+str(trajpos[2,it]) + '|' +                     \
733              dims['Y'][0]+':'+str(trajpos[1,it]) + '|' +                            \
734              dims['T'][0]+':'+str(coindtvalues[it][0])
735            slicevar, dimslice = slice_variable(ovsim, slicev)
736            simobsvalues.append([ ovsim[slicevar], ovobs[coindtvalues[it][1]]])
737    elif obskind == 'single-station':
738        for it in range(Ncoindt):
739            slicev = dims['X'][0]+':'+str(stationpos[1]) + '|' +                       \
740              dims['Y'][0]+':'+str(stationpos[0]) + '|' +                              \
741              dims['T'][0]+':'+str(coindtvalues[it][0])
742            slicevar, dimslice = slice_variable(ovsim, slicev)
743            simobsvalues.append([ slicevar, ovobs[coindtvalues[it][1]]])
744    elif obskind == 'trajectory':
745        if dims.has_key('Z'):
746            for it in range(Ncoindt):
747                slicev = dims['X'][0]+':'+str(trajpos[2,it]) + '|' +                 \
748                  dims['Y'][0]+':'+str(trajpos[1,it]) + '|' +                        \
749                  dims['Z'][0]+':'+str(trajpos[0,it]) + '|' +                        \
750                  dims['T'][0]+':'+str(coindtvalues[it][0])
751                slicevar, dimslice = slice_variable(ovsim, slicev)
752                simobsvalues.append([ ovsim[slicevar], ovobs[coindtvalues[it][1]]])
753                print simobsvalues[varsimobs][:][it]
754        else:
755            for it in range(Ncoindt):
756                slicev = dims['X'][0]+':'+str(trajpos[2,it]) + '|' +                 \
757                  dims['Y'][0]+':'+str(trajpos[1,it]) + '|' +                        \
758                  dims['T'][0]+':'+str(coindtvalues[it][0])
759                slicevar, dimslice = slice_variable(ovsim, slicev)
760                simobsvalues.append([ ovsim[slicevar], ovobs[coindtvalues[it][1]]])
761                print simobsvalues[varsimobs][:][it]
762
763    newvar = onewnc.createVariable(varsimobs, 'f', ('time','couple'))
764    descvar = 'couples of simulated: ' + valvars[ivar][0] + ' and observed ' +        \
765      valvars[ivar][1]
766    basicvardef(newvar, varsimobs, descvar, ovobs.getncattr('units'))
767 
768    newvar[:] = np.array(simobsvalues)
769
770    onewnc.sync()
771
772# Global attributes
773##
774set_attribute(onewnc,'author_nc','Lluis Fita')
775set_attribute(onewnc,'institution_nc','Laboratoire de Meteorology Dynamique, ' +    \
776  'LMD-Jussieu, UPMC, Paris')
777set_attribute(onewnc,'country_nc','France')
778set_attribute(onewnc,'script_nc',main)
779set_attribute(onewnc,'version_script',version)
780set_attribute(onewnc,'information',                                                 \
781  'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html')
782set_attribute(onewnc,'simfile',opts.fsim)
783set_attribute(onewnc,'obsfile',opts.fobs)
784
785onewnc.sync()
786onewnc.close()
787
788print main + ": successfull writting of '" + ofile + "' !!"
Note: See TracBrowser for help on using the repository browser.