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

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

Creation of the script to validate the simulations
File size: 18.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
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
22def index_3mat(matA,matB,matC,val):
23    """ Function to provide the coordinates of a given value inside three matrix simultaneously
24    index_mat(matA,matB,matC,val)
25      matA= matrix with one set of values
26      matB= matrix with the other set of values
27      matB= matrix with the third set of values
28      val= triplet of values to search
29    >>> 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])
30    [2 1 1]
31    """
32    fname = 'index_3mat'
33
34    matAshape = matA.shape
35    matBshape = matB.shape
36    matCshape = matC.shape
37
38    for idv in range(len(matAshape)):
39        if matAshape[idv] != matBshape[idv]:
40            print errormsg
41            print '  ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv],  \
42              'and B:',matBshape[idv],'does not coincide!!'
43            quit(-1)
44        if matAshape[idv] != matCshape[idv]:
45            print errormsg
46            print '  ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv],  \
47              'and C:',matCshape[idv],'does not coincide!!'
48            quit(-1)
49
50    minA = np.min(matA)
51    maxA = np.max(matA)
52    minB = np.min(matB)
53    maxB = np.max(matB)
54    minC = np.min(matC)
55    maxC = np.max(matC)
56
57    if val[0] < minA or val[0] > maxA:
58        print warnmsg
59        print '  ' + fname + ': first value:',val[0],'outside matA range',minA,',',  \
60          maxA,'!!'
61    if val[1] < minB or val[1] > maxB:
62        print warnmsg
63        print '  ' + fname + ': second value:',val[1],'outside matB range',minB,',',  \
64          maxB,'!!'
65    if val[2] < minC or val[2] > maxC:
66        print warnmsg
67        print '  ' + fname + ': second value:',val[2],'outside matC range',minC,',',  \
68          maxC,'!!'
69
70    dist = np.zeros(tuple(matAshape), dtype=np.float)
71    dist = np.sqrt((matA - np.float(val[0]))**2 + (matB - np.float(val[1]))**2 +     \
72      (matC - np.float(val[2]))**2)
73
74    mindist = np.min(dist)
75   
76    matlist = list(dist.flatten())
77    ifound = matlist.index(mindist)
78
79    Ndims = len(matAshape)
80    valpos = np.zeros((Ndims), dtype=int)
81    baseprevdims = np.zeros((Ndims), dtype=int)
82
83    for dimid in range(Ndims):
84        baseprevdims[dimid] = np.product(matAshape[dimid+1:Ndims])
85        if dimid == 0:
86            alreadyplaced = 0
87        else:
88            alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
89        valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
90
91    return valpos
92
93def index_2mat(matA,matB,val):
94    """ Function to provide the coordinates of a given value inside two matrix simultaneously
95    index_mat(matA,matB,val)
96      matA= matrix with one set of values
97      matB= matrix with the pother set of values
98      val= couple of values to search
99    >>> index_mat(np.arange(27).reshape(3,3,3),np.arange(100,127).reshape(3,3,3),[22,111])
100    [2 1 1]
101    """
102    fname = 'index_2mat'
103
104    matAshape = matA.shape
105    matBshape = matB.shape
106
107    for idv in range(len(matAshape)):
108        if matAshape[idv] != matBshape[idv]:
109            print errormsg
110            print '  ' + fname + ': Dimension',idv,'of matrices A:',matAshape[idv],  \
111              'and B:',matBshape[idv],'does not coincide!!'
112            quit(-1)
113
114    minA = np.min(matA)
115    maxA = np.max(matA)
116    minB = np.min(matB)
117    maxB = np.max(matB)
118
119    if val[0] < minA or val[0] > maxA:
120        print warnmsg
121        print '  ' + fname + ': first value:',val[0],'outside matA range',minA,',',  \
122          maxA,'!!'
123    if val[1] < minB or val[1] > maxB:
124        print warnmsg
125        print '  ' + fname + ': second value:',val[1],'outside matB range',minB,',',  \
126          maxB,'!!'
127
128    dist = np.zeros(tuple(matAshape), dtype=np.float)
129    dist = np.sqrt((matA - np.float(val[0]))**2 + (matB - np.float(val[1]))**2)
130
131    mindist = np.min(dist)
132   
133    matlist = list(dist.flatten())
134    ifound = matlist.index(mindist)
135
136    Ndims = len(matAshape)
137    valpos = np.zeros((Ndims), dtype=int)
138    baseprevdims = np.zeros((Ndims), dtype=int)
139
140    for dimid in range(Ndims):
141        baseprevdims[dimid] = np.product(matAshape[dimid+1:Ndims])
142        if dimid == 0:
143            alreadyplaced = 0
144        else:
145            alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
146        valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
147
148    return valpos
149
150def index_mat(mat,val):
151    """ Function to provide the coordinates of a given value inside a matrix
152    index_mat(mat,val)
153      mat= matrix with values
154      val= value to search
155    >>> index_mat(np.arange(27).reshape(3,3,3),22)
156    [2 1 1]
157    """
158
159    fname = 'index_mat'
160
161    matshape = mat.shape
162
163    matlist = list(mat.flatten())
164    ifound = matlist.index(val)
165
166    Ndims = len(matshape)
167    valpos = np.zeros((Ndims), dtype=int)
168    baseprevdims = np.zeros((Ndims), dtype=int)
169
170    for dimid in range(Ndims):
171        baseprevdims[dimid] = np.product(matshape[dimid+1:Ndims])
172        if dimid == 0:
173            alreadyplaced = 0
174        else:
175            alreadyplaced = np.sum(baseprevdims[0:dimid]*valpos[0:dimid])
176        valpos[dimid] = int((ifound - alreadyplaced )/ baseprevdims[dimid])
177
178    return valpos
179
180def coincident_CFtimes(tvalB, tunitA, tunitB):
181    """ Function to make coincident times for two different sets of CFtimes
182    tvalB= time values B
183    tunitA= time units times A to which we want to make coincidence
184    tunitB= time units times B
185    >>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00',
186      'hours since 1949-12-01 00:00:00')
187    [     0.   3600.   7200.  10800.  14400.  18000.  21600.  25200.  28800.  32400.]
188    >>> coincident_CFtimes(np.arange(10),'seconds since 1949-12-01 00:00:00',
189      'hours since 1979-12-01 00:00:00')
190    [  9.46684800e+08   9.46688400e+08   9.46692000e+08   9.46695600e+08
191       9.46699200e+08   9.46702800e+08   9.46706400e+08   9.46710000e+08
192       9.46713600e+08   9.46717200e+08]
193    """
194    import datetime as dt
195    fname = 'coincident_CFtimes'
196
197    trefA = tunitA.split(' ')[2] + ' ' + tunitA.split(' ')[3]
198    trefB = tunitB.split(' ')[2] + ' ' + tunitB.split(' ')[3]
199    tuA = tunitA.split(' ')[0]
200    tuB = tunitB.split(' ')[0]
201
202    if tuA != tuB:
203        if tuA == 'microseconds':
204            if tuB == 'microseconds':
205                tB = tvalB*1.
206            elif tuB == 'seconds':
207                tB = tvalB*10.e6
208            elif tuB == 'minutes':
209                tB = tvalB*60.*10.e6
210            elif tuB == 'hours':
211                tB = tvalB*3600.*10.e6
212            elif tuB == 'days':
213                tB = tvalB*3600.*24.*10.e6
214            else:
215                print errormsg
216                print '  ' + fname + ": combination of time untis: '" + tuA +        \
217                  "' & '" + tuB + "' not ready !!"
218                quit(-1)
219        elif tuA == 'seconds':
220            if tuB == 'microseconds':
221                tB = tvalB/10.e6
222            elif tuB == 'seconds':
223                tB = tvalB*1.
224            elif tuB == 'minutes':
225                tB = tvalB*60.
226            elif tuB == 'hours':
227                tB = tvalB*3600.
228            elif tuB == 'days':
229                tB = tvalB*3600.*24.
230            else:
231                print errormsg
232                print '  ' + fname + ": combination of time untis: '" + tuA +        \
233                  "' & '" + tuB + "' not ready !!"
234                quit(-1)
235        elif tuA == 'minutes':
236            if tuB == 'microseconds':
237                tB = tvalB/(60.*10.e6)
238            elif tuB == 'seconds':
239                tB = tvalB/60.
240            elif tuB == 'minutes':
241                tB = tvalB*1.
242            elif tuB == 'hours':
243                tB = tvalB*60.
244            elif tuB == 'days':
245                tB = tvalB*60.*24.
246            else:
247                print errormsg
248                print '  ' + fname + ": combination of time untis: '" + tuA +        \
249                  "' & '" + tuB + "' not ready !!"
250                quit(-1)
251        elif tuA == 'hours':
252            if tuB == 'microseconds':
253                tB = tvalB/(3600.*10.e6)
254            elif tuB == 'seconds':
255                tB = tvalB/3600.
256            elif tuB == 'minutes':
257                tB = tvalB/60.
258            elif tuB == 'hours':
259                tB = tvalB*1.
260            elif tuB == 'days':
261                tB = tvalB*24.
262            else:
263                print errormsg
264                print '  ' + fname + ": combination of time untis: '" + tuA +        \
265                  "' & '" + tuB + "' not ready !!"
266                quit(-1)
267        elif tuA == 'days':
268            if tuB == 'microseconds':
269                tB = tvalB/(24.*3600.*10.e6)
270            elif tuB == 'seconds':
271                tB = tvalB/(24.*3600.)
272            elif tuB == 'minutes':
273                tB = tvalB/(24.*60.)
274            elif tuB == 'hours':
275                tB = tvalB/24.
276            elif tuB == 'days':
277                tB = tvalB*1.
278            else:
279                print errormsg
280                print '  ' + fname + ": combination of time untis: '" + tuA +        \
281                  "' & '" + tuB + "' not ready !!"
282                quit(-1)
283        else:
284            print errormsg
285            print '  ' + fname + ": time untis: '" + tuA + "' not ready !!"
286            quit(-1)
287    else:
288        tB = tvalB*1.
289
290    if trefA != trefB:
291        trefTA = dt.datetime.strptime(trefA, '%Y-%m-%d %H:%M:%S')
292        trefTB = dt.datetime.strptime(trefB, '%Y-%m-%d %H:%M:%S')
293
294        difft = trefTB - trefTA
295        diffv = difft.days*24.*3600.*10.e6 + difft.seconds*10.e6 + difft.microseconds
296        print '  ' + fname + ': different reference refA:',trefTA,'refB',trefTB
297        print '    difference:',difft,':',diffv,'microseconds'
298
299        if tuA == 'microseconds':
300            tB = tB + diffv
301        elif tuA == 'seconds':
302            tB = tB + diffv/10.e6
303        elif tuA == 'minutes':
304            tB = tB + diffv/(60.*10.e6)
305        elif tuA == 'hours':
306            tB = tB + diffv/(3600.*10.e6)
307        elif tuA == 'dayss':
308            tB = tB + diffv/(24.*3600.*10.e6)
309        else:
310            print errormsg
311            print '  ' + fname + ": time untis: '" + tuA + "' not ready !!"
312            quit(-1)
313
314    return tB
315
316####### ###### ##### #### ### ## #
317
318strCFt="Refdate,tunits (CF reference date [YYYY][MM][DD][HH][MI][SS] format and " +  \
319  " and time units: 'weeks', 'days', 'hours', 'miuntes', 'seconds')"
320
321kindobs=['multi-points', 'single-station', 'trajectory']
322strkObs="kind of observations; 'multi-points': multiple individual punctual obs " +  \
323  "(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
324  "'trajectory': following a trajectory"
325#sumc,[constant]: add [constant] to variables values
326#subc,[constant]: substract [constant] to variables values
327#mulc,[constant]: multipy by [constant] to variables values
328#divc,[constant]: divide by [constant] to variables values
329
330parser = OptionParser()
331parser.add_option("-d", "--dimensions", dest="dims",
332  help="[DIM]@[simdim]@[obsdim] ',' list of couples of dimensions names from each source ([DIM]='X','Y','Z','T'; None, no value)",
333  metavar="VALUES")
334parser.add_option("-D", "--vardimensions", dest="vardims",
335  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")
336parser.add_option("-k", "--kindObs", dest="obskind", type='choice', choices=kindobs, 
337  help=strkObs, metavar="FILE")
338parser.add_option("-l", "--stationLocation", dest="stloc", 
339  help="longitude, latitude and height of the station (only for 'single-station')", 
340  metavar="FILE")
341parser.add_option("-o", "--observation", dest="fobs",
342  help="observations file to validate", metavar="FILE")
343parser.add_option("-s", "--simulation", dest="fsim",
344  help="simulation file to validate", metavar="FILE")
345parser.add_option("-v", "--variables", dest="vars",
346  help="[simvar]@[obsvar]@[[oper]@[val]] ',' list of couples of variables to validate and if necessary operation and value", metavar="VALUES")
347
348(opts, args) = parser.parse_args()
349
350#######    #######
351## MAIN
352    #######
353
354ofile='validation_sim.nc'
355
356if opts.dims is None:
357    print errormsg
358    print '  ' + main + ': No list of dimensions are provided!!'
359    print '    a ',' list of values X@[dimxsim]@[dimxobs],...,T@[dimtsim]@[dimtobs]'+\
360      ' is needed'
361    quit(-1)
362else:
363    print main +': couple of dimensions _______'
364    dims = {}
365    ds = opts.dims.split(',')
366    for d in ds:
367        dsecs = d.split('@')
368        if len(dsecs) != 3:
369            print errormsg
370            print '  ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
371            print '    [DIM]@[dimnsim]@[dimnobs]'
372            quit(-1)
373        dims[dsecs[0]] = [dsecs[1], dsecs[2]]
374        print dsecs[0],':',dsecs[1],',',dsecs[2]
375       
376
377if opts.vardims is None:
378    print errormsg
379    print '  ' + main + ': No list of variables with dimension values are provided!!'
380    print '    a ',' list of values X@[vardimxsim]@[vardimxobs],...,T@' +  \
381      '[vardimtsim]@[vardimtobs] is needed'
382    quit(-1)
383else:
384    print main +': couple of variable dimensions _______'
385    vardims = {}
386    ds = opts.vardims.split(',')
387    for d in ds:
388        dsecs = d.split('@')
389        if len(dsecs) != 3:
390            print errormsg
391            print '  ' + main + ': wrong number of values in:',dsecs,' 3 are needed !!'
392            print '    [DIM]@[vardimnsim]@[vardimnobs]'
393            quit(-1)
394        vardims[dsecs[0]] = [dsecs[1], dsecs[2]]
395        print dsecs[0],':',dsecs[1],',',dsecs[2]
396
397if opts.obskind is None:
398    print errormsg
399    print '  ' + main + ': No kind of observations provided !!'
400    quit(-1)
401else:
402    obskind = opts.obskind
403    if obskind == 'single-station':
404        if opts.stloc is None:
405            print errormsg
406            print '  ' + main + ': No station location provided !!'
407            quit(-1)
408        else:
409            stationdesc = [np.float(opts.stloc.split(',')[0]),                       \
410              np.float(opts.stloc.split(',')[1]), np.float(opts.stloc.split(',')[2])]
411
412if opts.fobs is None:
413    print errormsg
414    print '  ' + main + ': No observations file is provided!!'
415    quit(-1)
416else:
417    if not os.path.isfile(opts.fobs):
418        print errormsg
419        print '   ' + main + ": observations file '" + opts.fobs + "' does not exist !!"
420        quit(-1)
421
422if opts.fsim is None:
423    print errormsg
424    print '  ' + main + ': No simulation file is provided!!'
425    quit(-1)
426else:
427    if not os.path.isfile(opts.fsim):
428        print errormsg
429        print '   ' + main + ": simulation file '" + opts.fsim + "' does not exist !!"
430        quit(-1)
431
432if opts.vars is None:
433    print errormsg
434    print '  ' + main + ': No list of couples of variables is provided!!'
435    print '    a ',' list of values [varsim]@[varobs],... is needed'
436    quit(-1)
437else:
438    valvars = []
439    vs = opts.dims.split(',')
440    for v in vs:
441        vsecs = v.split('@')
442        if len(dsecs) < 2:
443            print errormsg
444            print '  ' + main + ': wrong number of values in:',vsecs,                \
445              ' at least 2 are needed !!'
446            print '    [varsim]@[varobs]@[[oper][val]]'
447            quit(-1)
448        valvars.append(vsecs)
449
450# Openning observations trajectory
451##
452oobs = NetCDFFile(opts.fobs, 'r')
453
454valdimobs = {}
455for dn in dims:
456    print dn,':',dims[dn]
457    if dims[dn][1] != 'None':
458        if not oobs.dimensions.has_key(dims[dn][1]):
459            print errormsg
460            print '  ' + main + ": observations file does not have dimension '" +    \
461              dims[dn][1] + "' !!"
462            quit(-1)
463        if vardims[dn][1] != 'None':
464            if not oobs.variables.has_key(vardims[dn][1]):
465                print errormsg
466                print '  ' + main + ": observations file does not have varibale " +  \
467                  "dimension '" + vardims[dn][1] + "' !!"
468                quit(-1)
469            valdimobs[dn] = oobs.variables[vardims[dn][1]][:]
470    else:
471        if dn == 'X':
472            valdimobs[dn] = stationdesc[0]
473        elif dn == 'Y':
474            valdimobs[dn] = stationdesc[1]
475        elif dn == 'Z':
476            valdimobs[dn] = stationdesc[2]
477
478osim = NetCDFFile(opts.fsim, 'r')
479
480valdimsim = {}
481for dn in dims:
482    if not osim.dimensions.has_key(dims[dn][0]):
483        print errormsg
484        print '  ' + main + ": simulation file does not have dimension '" +          \
485          dims[dn][0] + "' !!"
486        quit(-1)
487    if not osim.variables.has_key(vardims[dn][0]):
488        print errormsg
489        print '  ' + main + ": simulation file does not have varibale dimension '" + \
490          vardims[dn][0] + "' !!"
491        quit(-1)
492    valdimsim[dn] = osim.variables[vardims[dn][0]][:]
493
494# General characteristics
495dimtobs = len(valdimobs['T'])
496dimtsim = len(valdimsim['T'])
497
498print main +': observational time-steps:',dimtobs,'simulation:',dimtsim
499
500if obskind == 'multi-points':
501    trajpos = np.zeros((2,dimt),dtype=int)
502    for it in dimtobs:
503        trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'],                    \
504          [valdimobs['X'][it],valdimobss['Y'][it]])
505
506elif obskind == 'single-station':
507    stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],[valdimobs['Y'],             \
508      valdimobs['X']])
509    print main + ': station point in simulation:', stsimpos
510    print '    station position:',valdimobs['X'],',',valdimobs['Y']
511    print '    simulation coord.:',valdimsim['X'][tuple(stsimpos)],',',              \
512      valdimsim['Y'][tuple(stsimpos)]
513elif obskind == 'trajectory':
514    if dims.has_key('Z'):
515        trajpos = np.zeros((3,dimt),dtype=int)
516        for it in dimtobs:
517            trajpos[0:1,it] = index_2mat(valdimsim['X'],valdimsim['Y'],              \
518              [valdimobs['X'][it],valdimobss['Y'][it]])
519            trajpos[2,it] = index_mat(valdimsim['Z'],valdimobs['Z'][it])
520    else:
521        trajpos = np.zeros((2,dimt),dtype=int)
522        for it in dimtobs:
523            trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'],                \
524              [valdimobs['X'][it],valdimobss['Y'][it]])
525
526# Getting times
527tobj = oobs.variables[vardims['T'][1]]
528obstunits = tobj.getncattr('units')
529tobj = osim.variables[vardims['T'][0]]
530simtunits = tobj.getncattr('units')
531
532simobstimes = coincident_CFtimes(valdimsim['T'], obstunits, simtunits)
533
534print 'Lluis:',simobstimes
535
Note: See TracBrowser for help on using the repository browser.