Changeset 333 in lmdz_wrf

Timestamp:

Feb 27, 2015, 4:09:03 PM (10 years ago)

Author:

lfita

Message:

Working and tested version for 'single-station'

File:

• trunk/tools/validation_sim.py (modified) (10 diffs)

trunk/tools/validation_sim.py

@@ -1 +1 @@
 
 # L. Fita, LMD-Jussieu. February 2015
-## 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 
+## 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
 import numpy as np
 import os
@@ -19 +19 @@
 fillValueI = -99999
 fillValueS = '---'
+
+StringLength = 50
+
+def searchInlist(listname, nameFind):
+    """ Function to search a value within a list
+    listname = list
+    nameFind = value to find
+    >>> searInlist(['1', '2', '3', '5'], '5')
+    True
+    """
+    for x in listname:
+      if x == nameFind:
+        return True
+    return False
+
+def set_attribute(ncvar, attrname, attrvalue):
+    """ Sets a value of an attribute of a netCDF variable. Removes previous attribute value if exists
+    ncvar = object netcdf variable
+    attrname = name of the attribute
+    attrvalue = value of the attribute
+    """
+    import numpy as np
+    from netCDF4 import Dataset as NetCDFFile
+
+    attvar = ncvar.ncattrs()
+    if searchInlist(attvar, attrname):
+        attr = ncvar.delncattr(attrname)
+
+    attr = ncvar.setncattr(attrname, attrvalue)
+
+    return ncvar
+
+def basicvardef(varobj, vstname, vlname, vunits):
+    """ Function to give the basic attributes to a variable
+    varobj= netCDF variable object
+    vstname= standard name of the variable
+    vlname= long name of the variable
+    vunits= units of the variable
+    """
+    attr = varobj.setncattr('standard_name', vstname)
+    attr = varobj.setncattr('long_name', vlname)
+    attr = varobj.setncattr('units', vunits)
+
+    return
+
+def writing_str_nc(varo, values, Lchar):
+    """ Function to write string values in a netCDF variable as a chain of 1char values
+    varo= netCDF variable object
+    values = list of values to introduce
+    Lchar = length of the string in the netCDF file
+    """
+
+    Nvals = len(values)
+    for iv in range(Nvals):    
+        stringv=values[iv]  
+        charvals = np.chararray(Lchar)
+        Lstr = len(stringv)
+        charvals[Lstr:Lchar] = ''
+
+        for ich in range(Lstr):
+            charvals[ich] = stringv[ich:ich+1]
+
+        varo[iv,:] = charvals
+
+    return
 
 def index_3mat(matA,matB,matC,val):
@@ -314 +379 @@
     return tB
 
+
+def slice_variable(varobj, dimslice):
+    """ Function to return a slice of a given variable according to values to its 
+      dimensions
+    slice_variable(varobj, dimslice)
+      varobj= object wit the variable
+      dimslice= [[dimname1]:[value1]|[[dimname2]:[value2], ...] pairs of dimension 
+        [value]: 
+          * [integer]: which value of the dimension
+          * -1: all along the dimension
+          * -9: last value of the dimension
+          * [beg]@[end] slice from [beg] to [end]
+    """
+    fname = 'slice_variable'
+
+    if varobj == 'h':
+        print fname + '_____________________________________________________________'
+        print slice_variable.__doc__
+        quit()
+
+    vardims = varobj.dimensions
+    Ndimvar = len(vardims)
+
+    Ndimcut = len(dimslice.split('|'))
+    dimsl = dimslice.split('|')
+
+    varvalsdim = []
+    dimnslice = []
+
+    for idd in range(Ndimvar):
+        for idc in range(Ndimcut):
+            dimcutn = dimsl[idc].split(':')[0]
+            dimcutv = dimsl[idc].split(':')[1]
+            if vardims[idd] == dimcutn: 
+                posfrac = dimcutv.find('@')
+                if posfrac != -1:
+                    inifrac = int(dimcutv.split('@')[0])
+                    endfrac = int(dimcutv.split('@')[1])
+                    varvalsdim.append(slice(inifrac,endfrac))
+                    dimnslice.append(vardims[idd])
+                else:
+                    if int(dimcutv) == -1:
+                        varvalsdim.append(slice(0,varobj.shape[idd]))
+                        dimnslice.append(vardims[idd])
+                    elif int(dimcutv) == -9:
+                        varvalsdim.append(int(varobj.shape[idd])-1)
+                    else:
+                        varvalsdim.append(int(dimcutv))
+                break
+
+    varvalues = varobj[tuple(varvalsdim)]
+
+    return varvalues, dimnslice
+
+
 ####### ###### ##### #### ### ## #
 
@@ -323 +443 @@
   "(e.g., lightning strikes), 'single-station': single station on a fixed position,"+\
   "'trajectory': following a trajectory"
+simopers =['sumc','subc','mulc','divc']
 #sumc,[constant]: add [constant] to variables values
 #subc,[constant]: substract [constant] to variables values
@@ -372 +493 @@
             quit(-1)
         dims[dsecs[0]] = [dsecs[1], dsecs[2]]
-        print dsecs[0],':',dsecs[1],',',dsecs[2]
+        print '  ',dsecs[0],':',dsecs[1],',',dsecs[2]
         
-
 if opts.vardims is None:
     print errormsg
@@ -393 +513 @@
             quit(-1)
         vardims[dsecs[0]] = [dsecs[1], dsecs[2]]
-        print dsecs[0],':',dsecs[1],',',dsecs[2]
+        print '  ',dsecs[0],':',dsecs[1],',',dsecs[2]
 
 if opts.obskind is None:
@@ -437 +557 @@
 else:
     valvars = []
-    vs = opts.dims.split(',')
+    vs = opts.vars.split(',')
     for v in vs:
         vsecs = v.split('@')
-        if len(dsecs) < 2:
+        if len(vsecs) < 2:
             print errormsg
             print '  ' + main + ': wrong number of values in:',vsecs,                \
@@ -446 +566 @@
             print '    [varsim]@[varobs]@[[oper][val]]'
             quit(-1)
+        if len(vsecs) > 2:
+            if not searchInlist(simopers,vsecs[2]): 
+                print errormsg
+                print main + ": operation on simulation values '" + vsecs[2] +       \
+                  "' not ready !!"
+                quit(-1)
+
         valvars.append(vsecs)
 
@@ -503 +630 @@
         trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'],                    \
           [valdimobs['X'][it],valdimobss['Y'][it]])
-
 elif obskind == 'single-station':
     stsimpos = index_2mat(valdimsim['Y'],valdimsim['X'],[valdimobs['Y'],             \
       valdimobs['X']])
-    print main + ': station point in simulation:', stsimpos
+    stationpos = np.zeros((2), dtype=int)
+    iid = 0
+    for idn in osim.variables[vardims['X'][0]].dimensions:
+        if idn == dims['X'][0]:
+            stationpos[1] = stsimpos[iid]
+        elif idn == dims['Y'][0]:
+            stationpos[0] = stsimpos[iid]
+
+        iid = iid + 1
+    print main + ': station point in simulation:', stationpos
     print '    station position:',valdimobs['X'],',',valdimobs['Y']
     print '    simulation coord.:',valdimsim['X'][tuple(stsimpos)],',',              \
       valdimsim['Y'][tuple(stsimpos)]
+
 elif obskind == 'trajectory':
     if dims.has_key('Z'):
         trajpos = np.zeros((3,dimt),dtype=int)
         for it in dimtobs:
-            trajpos[0:1,it] = index_2mat(valdimsim['X'],valdimsim['Y'],              \
-              [valdimobs['X'][it],valdimobss['Y'][it]])
+            trajpos[0:1,it] = index_2mat(valdimsim['Y'],valdimsim['X'],              \
+              [valdimobs['Y'][it],valdimobss['X'][it]])
             trajpos[2,it] = index_mat(valdimsim['Z'],valdimobs['Z'][it])
     else:
         trajpos = np.zeros((2,dimt),dtype=int)
         for it in dimtobs:
-            trajpos[:,it] = index_2mat(valdimsim['X'],valdimsim['Y'],                \
-              [valdimobs['X'][it],valdimobss['Y'][it]])
+            trajpos[:,it] = index_2mat(valdimsim['Y'],valdimsim['X'],                \
+              [valdimobs['Y'][it],valdimobss['X'][it]])
 
 # Getting times
@@ -532 +668 @@
 simobstimes = coincident_CFtimes(valdimsim['T'], obstunits, simtunits)
 
-print 'Lluis:',simobstimes
-
+# Concident times
+##
+coindtvalues = []
+for it in range(dimtsim):    
+    ot = 0
+    for ito in range(ot,dimtobs-1):
+        if valdimobs['T'][ito] < simobstimes[it] and valdimobs['T'][ito+1] >         \
+          simobstimes[it]:
+            ot = ito
+            tdist = simobstimes[it] - valdimobs['T'][ito]
+            coindtvalues.append([it,ito,simobstimes[it],valdimobs['T'][ito],tdist])
+
+Ncoindt = len(coindtvalues)
+print main + ': found',Ncoindt,'coincident times between simulation and observations'
+
+# Validating
+##
+
+onewnc = NetCDFFile(ofile, 'w')
+
+# Dimensions
+newdim = onewnc.createDimension('time',None)
+newdim = onewnc.createDimension('couple',2)
+newdim = onewnc.createDimension('StrLength',StringLength)
+
+# Variable dimensions
+##
+newvar = onewnc.createVariable('simtime','f8',('time'))
+basicvardef(newvar, 'simtime', 'time simulation', obstunits )
+set_attribute(newvar, 'calendar', 'standard')
+newvar[:] = coindtvalues[:][2]
+
+newvar = onewnc.createVariable('obstime','f8',('time'))
+basicvardef(newvar, 'obstime', 'time observations', obstunits )
+set_attribute(newvar, 'calendar', 'standard')
+newvar[:] = coindtvalues[:][3]
+
+newvar = onewnc.createVariable('couple', 'c', ('couple','StrLength'))
+basicvardef(newvar, 'couple', 'couples of values', '-')
+writing_str_nc(newvar, ['sim','obs'], StringLength)
+
+Nvars = len(valvars)
+for ivar in range(Nvars):
+    simobsvalues = []
+
+    varsimobs = valvars[ivar][0] + '_' + valvars[ivar][1]
+
+    if not oobs.variables.has_key(valvars[ivar][1]):
+        print errormsg
+        print '  ' + main + ": observations file has not '" + valvars[ivar][1] +     \
+          "' !!"
+        quit(-1)
+    if not osim.variables.has_key(valvars[ivar][0]):
+        print errormsg
+        print '  ' + main + ": simulation file has not '" + valvars[ivar][0] + "' !!"
+        quit(-1)
+
+    ovobs = oobs.variables[valvars[ivar][1]]
+    ovsim = osim.variables[valvars[ivar][0]]
+
+    if obskind == 'multi-points':
+        for it in range(Ncoindt):
+            slicev = dims['X'][0]+':'+str(trajpos[2,it]) + '|' +                     \
+              dims['Y'][0]+':'+str(trajpos[1,it]) + '|' +                            \
+              dims['T'][0]+':'+str(coindtvalues[it][0])
+            slicevar, dimslice = slice_variable(ovsim, slicev)
+            simobsvalues.append([ ovsim[slicevar], ovobs[coindtvalues[it][1]]])
+    elif obskind == 'single-station':
+        for it in range(Ncoindt):
+            slicev = dims['X'][0]+':'+str(stationpos[1]) + '|' +                       \
+              dims['Y'][0]+':'+str(stationpos[0]) + '|' +                              \
+              dims['T'][0]+':'+str(coindtvalues[it][0])
+            slicevar, dimslice = slice_variable(ovsim, slicev)
+            simobsvalues.append([ slicevar, ovobs[coindtvalues[it][1]]])
+    elif obskind == 'trajectory':
+        if dims.has_key('Z'):
+            for it in range(Ncoindt):
+                slicev = dims['X'][0]+':'+str(trajpos[2,it]) + '|' +                 \
+                  dims['Y'][0]+':'+str(trajpos[1,it]) + '|' +                        \
+                  dims['Z'][0]+':'+str(trajpos[0,it]) + '|' +                        \
+                  dims['T'][0]+':'+str(coindtvalues[it][0])
+                slicevar, dimslice = slice_variable(ovsim, slicev)
+                simobsvalues.append([ ovsim[slicevar], ovobs[coindtvalues[it][1]]])
+                print simobsvalues[varsimobs][:][it]
+        else:
+            for it in range(Ncoindt):
+                slicev = dims['X'][0]+':'+str(trajpos[2,it]) + '|' +                 \
+                  dims['Y'][0]+':'+str(trajpos[1,it]) + '|' +                        \
+                  dims['T'][0]+':'+str(coindtvalues[it][0])
+                slicevar, dimslice = slice_variable(ovsim, slicev)
+                simobsvalues.append([ ovsim[slicevar], ovobs[coindtvalues[it][1]]])
+                print simobsvalues[varsimobs][:][it]
+
+    newvar = onewnc.createVariable(varsimobs, 'f', ('time','couple'))
+    descvar = 'couples of simulated: ' + valvars[ivar][0] + ' and observed ' +        \
+      valvars[ivar][1]
+    basicvardef(newvar, varsimobs, descvar, ovobs.getncattr('units'))
+ 
+    newvar[:] = np.array(simobsvalues)
+
+    onewnc.sync()
+
+# Global attributes
+##
+set_attribute(onewnc,'author_nc','Lluis Fita')
+set_attribute(onewnc,'institution_nc','Laboratoire de Meteorology Dynamique, ' +    \
+  'LMD-Jussieu, UPMC, Paris')
+set_attribute(onewnc,'country_nc','France')
+set_attribute(onewnc,'script_nc',main)
+set_attribute(onewnc,'version_script',version)
+set_attribute(onewnc,'information',                                                 \
+  'http://www.lmd.jussieu.fr/~lflmd/ASCIIobs_nc/index.html')
+set_attribute(onewnc,'simfile',opts.fsim)
+set_attribute(onewnc,'obsfile',opts.fobs)
+
+onewnc.sync()
+onewnc.close()
+
+print main + ": successfull writting of '" + ofile + "' !!"