Changeset 1181 in lmdz_wrf for trunk

Timestamp:

Oct 12, 2016, 11:49:16 AM (9 years ago)

Author:

lfita

Message:

Adding `field_stats_dim': Function to retrieve statistics from a field along dimensions

Location:

trunk/tools

Files:

• nc_var.py (modified) (3 diffs)
• nc_var_tools.py (modified) (11 diffs)

trunk/tools/nc_var.py

@@ -21 +21 @@
 ## e.g. # nc_var.py -o pinterp -f /media/ExtDiskC_ext4/DATA/etudes/WRF_LMDZ/WaquaL_highres/short_copies/LMDZ/AR40/histins_19790101000000-19790304000000_short.nc -S 100000.:97500.:95000.:92500.:90000.:85000.:80000.:75000.:70000.:65000.:60000.:55000.:50000.:45000.:40000.:35000.:30000.:25000.:20000.:15000.:10000.:5000.:2500.:1000.:500.:250.,1,0 -v temp,ovap
 ## e.g. # nc_var.py -o reproject -f analysis/LMDZ/AR40/hurs_histins.nc -S 'lon,lat,analysis/WRF/current/hurs_wrfout.nc,lon,lat,npp,time@time' -v hurs
+## e.g. # nc_var.py -o field_stats_dim -f /home/lluis/PY/wrfout_d01_2001-11-11_00:00:00 -S 'full,None,None,west_east,XLONG,False' -v 'T2'
 
 from optparse import OptionParser
@@ -40 +41 @@
   'dimToUnlimited', 'dimVar_creation',                                               \
   'fattradd',                                                                        \
-  'fdimadd', 'fgaddattr', 'field_stats', 'file_creation', 'file_oper_alongdims',     \
-  'filter_2dim',                                                                     \
+  'fdimadd', 'fgaddattr', 'field_stats', 'field_stats_dim', 'file_creation',         \
+  'file_oper_alongdims', 'filter_2dim',                                              \
   'flipdim', 'fvaradd', 'gaddattrk', 'gaddattr', 'get_attribute',                    \
   'get_namelist_vars', 'get_Variables',                                              \
@@ -177 +178 @@
 elif oper == 'field_stats':
     ncvar.field_stats(opts.values, opts.ncfile, opts.varname)
+elif oper == 'field_stats_dim':
+    ncvar.field_stats_dim(opts.values, opts.ncfile, opts.varname)
 elif oper == 'filter_2dim':
     ncvar.filter_2dim(opts.values, opts.ncfile, opts.varname)

trunk/tools/nc_var_tools.py

@@ -44 +44 @@
 # fdimadd: Adding dimension from another reference file
 # field_stats: Function to retrieve statistics from a field
+# field_stats_dim: Function to retrieve statistics from a field along dimensions
 # file_creation: Operation to create a file with one variable with a given set of dimensions
 # file_oper_alongdims: Function to operate a file along different dimensions of a variable
@@ -228 +229 @@
     return
 
+def add_dims(oc,nc,dnames):
+    """ Function to add dimensions from a given netCDF object to another one
+      oc: source netcdfile object
+      nc: netcdfile object to add dimensions
+      dnames: list of name of dimensions to add
+    """
+    fname = 'add_dims'
+
+    for dn in dnames:
+        print '  ' + fname + ': adding dim:',dn,' ...'
+        if oc.dimensions[dn].isunlimited():
+            newdim = nc.createDimension(dn, None)
+        else:
+            newdim = nc.createDimension(dn, len(oc.dimensions[dn]))
+
+    nc.sync()
+
+    return
+
+def add_vars(oc,nc,vnames):
+    """ Function to add variables from a given netCDF object to another one
+      oc: source netcdfile object
+      nc: netcdfile object to add dimensions
+      vnames: list of name of variables to add
+    """
+    fname = 'add_vars'
+
+    for vn in vnames:
+        if not gen.searchInlist(nc.variables,vn):
+            print '  ' + fname + ': adding var:',vn,' ...'
+
+            varo = oc.variables[vn]
+            vartype = varo.dtype
+            vardims = varo.dimensions
+            varattrs = varo.ncattrs()
+
+            for vdn in vardims:
+                if not gen.searchInlist(nc.dimensions,vdn):
+                    print warnmsg
+                    print '  ' + fname + ": adding dimension '" + vdn +              \
+                      "' from variable '" + vdn + "' which is not in file !!"
+                    add_dims(oc,nc,[vdn])
+
+            if gen.searchInlist(varattrs,'_FillValue'):
+                newvar = nc.createVariable(vn, vartype, vardims,                     \
+                  fill_value=varo.getncattr('_FillValue'))
+            else:
+                newvar = nc.createVariable(vn, vartype, vardims)
+
+            for attrn in varattrs:
+                attrv = varo.getncattr(attrn)
+                newattr = set_attribute(newvar, attrn, attrv)
+
+            newvar[:] = varo[:]
+
+        nc.sync()
+
+    return
+
+def add_globattrs(oc,nc,attrns):
+    """ Function to add global attributes from a given netCDF object to another one
+      oc: source netcdfile object
+      nc: netcdfile object to add dimensions
+      vnames: list of name of dimensions to add ('all', for all attributes)
+    """
+    fname = 'add_globattrs'
+
+    allfattrs = oc.ncattrs()
+
+    if attrns == 'all':
+        fattrns = allfattrs
+    else:
+        fattrns = attrns
+
+    for an in fattrns:
+        if not gen.searchInlist(allfattrs,an):
+            print errormsg
+            print '  ' + fname + ": file has not global attribute '" + an + "' !!"
+            quit(-1)
+
+        print '  ' + fname + ': adding attr:',an,' ...'
+        
+        av = oc.getncattr(an)
+        newattr = set_attribute(nc, an, av)
+
+    nc.sync()
+
+    return
+
 def valmod(values, ncfile, varn):
     """ Function to modify the value of a variable
@@ -5968 +6058 @@
                     varslice.append(slice(0,Ldim))
                     vardims.append(dimname)
-                print fname + ';Lluis dimname:', dimname,'dimslices:',dimslices[dimname]
-
-            print fname + '; Lluis coinc:', coinc
-            print fname + ' Lluis; shapes varvals:', varvals, 'varobj:', varobj.shape,'varslice:',varslice 
             varvals = varobj[tuple(varslice)]
 
@@ -6269 +6355 @@
             else:
                 varslice.append(slice(0,len(objfile.dimensions[dimn])))
-                print fname + '; Lluis here! dimn', dimn
 
         if varinf.FillValue is not None:
@@ -6278 +6363 @@
 
         if dimx != dimy:
-            print fname + '; Lluis; varslice:', varslice,'shape newvar:', newvar.shape
             newvar[:] = varobj[tuple(varslice)]
         else:
@@ -6311 +6395 @@
     for atrn in objfile.ncattrs():
         attrv = objfile.getncattr(atrn)
-        print fname + '; Lluis attrn, attrv, type:', atrn, attrv, type(attrv)
         newattr = set_attributek(newncobj, atrn, attrv, type(attrv))
 
@@ -10822 +10905 @@
                 convals = [int(0)]
                 convalsS = ''
-                print fname + '; Lluis here convals:', convals
             statsvariable[vn] = [minv, maxv, meanv, mean2v, varv, convals[0]]
         else:
@@ -10843 +10925 @@
 
 #field_stats('full', 'geo_em.d01.nc', 'HGT_M')
+
+def field_stats_dim(values, ncfile, varn):
+    """ Function to retrieve statistics from a field along dimensions
+    field_stats_dim(values, ncfile, varn)
+      [values]= [stats],[fillVals],[dimns],[vardims],[stdout]
+        [stats]: kind of statistics 
+          'full': all statistics given variable (min, max, mean, mean2, variability, [countVals])
+        [fillVals]: ':' list of _fillValues ('None' for any)
+        [countVals]: ':' list of Values@cond to use to count conditions ('None' for any)
+        [dims]: ':', list of names of the dimension along which statistics is required
+        [vardims]: ':', list of names of the variable with the values for the dimension along which statistics is required
+        [stdout]: Whether should be printed in standard output the results ('True'/'False')
+      [ncfile]= name of the netCDF file to use
+      [varn]= variable name to use ('all' for all variables)
+    """
+    import numpy.ma as ma
+    fname='field_stats_dim'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print field_stats_dim.__doc__
+        quit()
+
+    Lstr = 100
+    fullstats = ['min', 'max', 'mean', 'mean2', 'variability']
+
+    arguments = '[stats],[fillVals],[coutnVals],[dims],[vardims],[stdout]'
+    gen.check_arguments(fname,values,arguments,',')
+
+    ofile = fname + '.nc'
+
+    stats=values.split(',')[0]
+    fillVals=values.split(',')[1]
+    countVals=values.split(',')[2]
+    dims=gen.str_list(values.split(',')[3],':')
+    vardims=gen.str_list(values.split(',')[4],':')
+    stdout=gen.Str_Bool(values.split(',')[5])
+
+# Fill Values
+    if fillVals == 'None':
+        fillV = None
+    else:
+        if fillVals.find(':') != -1:
+            fillV = fillVals.split(':')
+            NfillV = len(fillV)
+        else:
+            fillV = [fillVals]
+            NfillV = 1
+
+    # Statistics
+    if stats == 'full':
+        statns = fullstats
+    else:
+        print errormsg
+        print '  ' + fname + ": kind of statistics '" + stats + "' not ready!!"
+        print '    valid statistics: full'
+        quit(-1)
+
+# Count Values
+    if countVals == 'None':
+        countV = None
+        NcountV = 0
+    else:
+        if countVals.find(':') != -1:
+            countV = countVals.split(':')
+            NcountV = len(countV)
+            for cV in countV:
+                statns = statns.append('count_' + str(cV))
+        else:
+            countV = [countVals]
+            NcountV = 1
+            statns = statns.append('count_' + str(countV))
+
+    Nstats = len(statns)
+
+    ncobj = NetCDFFile(ncfile, 'r')
+    for dim in dims:
+        if not gen.searchInlist(ncobj.dimensions.keys(), dim):
+            print errormsg
+            print '  ' +fname+ ": file '" + ncfile + "' does not have dimension '" + \
+              dim + "' !!"
+            quit(-1)
+    for vardim in vardims:
+        if not gen.searchInlist(ncobj.variables.keys(), vardim):
+            print errormsg
+            print '  ' + fname + ": file '" + ncfile + "' does not have variable-" + \
+              "dimension '" + vardim + "' !!"
+            quit(-1)
+
+    if varn == 'all':
+        varstats = list(ncobj.variables)
+    else:
+        varstats = [varn]
+
+    # File creation
+    onewnc = NetCDFFile(ofile, 'w')
+
+    # Dimensions
+    stdims = {}
+    rundims = []
+    runndims = []
+    dimslice = {}
+    for dim in dims:
+        add_dims(ncobj,onewnc,[dim])
+        stdims[dim] = len(ncobj.dimensions[dim])
+        runndims.append(dim)
+        rundims.append(stdims[dim])
+        dimslice[dim] = 0
+    onewnc.createDimension('stats',Nstats)
+    onewnc.createDimension('Lstring',Lstr)
+
+    # Dimension variables
+    add_vars(ncobj,onewnc,vardims)
+
+    # Statistics var
+    newvar = onewnc.createVariable('stats','c',('stats','Lstring'))
+    vardef = basicvardef(newvar, 'stats', 'statistics along running dimensions','-')
+    newvals = writing_str_nc(newvar, statns, Lstr)
+    onewnc.sync()
+
+    for vn in varstats:
+        statsvariable = {}
+        Ntotrun = 1
+        if stdout:
+            print '  ' + fname + " file: '" + ncfile + "'..."
+            print '  ' + fname + '______ ______ _____ ____ ___ __ _'
+        for vn in varstats:
+            if not ncobj.variables.has_key(vn):
+                print errormsg
+                print '  ' + fname + ": file do not have variable '" + vn + "' !!"
+                quit(-1)
+    
+        objfield = ncobj.variables[vn]
+
+        dtype = objfield.dtype
+        # Dimensions to run along
+        runvdimns = []
+        runvLdims = []
+        for vdn in objfield.dimensions:
+            if gen.searchInlist(runndims,vdn):
+                Ntotrun = Ntotrun * stdims[vdn]
+                runvdimns.append(vdn)
+                runvLdims.append(stdims[vdn])
+        statsresults = np.zeros(tuple([Nstats] + runvLdims), dtype=np.float)
+        if len(runvLdims) == 0:
+            print warnmsg
+            print '  ' + fname + ': variable without running dimension!'
+            print '    computing statistics for the whole variable'
+            Ntotrun = 1
+
+        for idd in range(Ntotrun):
+            # Getting variable slice at this iteration
+            if len(runvLdims) > 0:
+                slicedims = gen.index_flatten_mat(idd,runvLdims)
+            else:
+                slicedims = []
+            slicevar = []
+            for vdn in objfield.dimensions:
+                runfound = False
+                Srundim = ''
+                for irvdn in range(len(runvdimns)):
+                    if vdn == runvdimns[irvdn]:
+                        slicevar.append(slicedims[irvdn])
+                        runfound = True
+                        Srundim = Srundim + vdn + ': ' + str(slicedims[irvdn]) + ' '
+                        break
+                if not runfound:
+                    slicevar.append(slice(0,len(ncobj.dimensions[vdn])))
+            
+            # Getting field values
+            field = objfield[tuple(slicevar)]
+            
+            # Computing statistics
+            counts = []
+            if countV is not None:
+                for icV in range(NcountV):
+                    cV = countV[icV].split('@')[0]
+                    cC = countV[icV].split('@')[1]
+                    countval = gen.retype(cV, dtype)
+                    counts.append([gen.count_cond(field, countval, cC), cC, countval])
+
+            if fillV is not None:
+                for ifV in range(NfillV):
+                    fillval = gen.retype(fillV[ifV], dtype)
+                    field = ma.masked_equal(field, fillval)
+
+            if stdout:
+                print '   ' + vn + '[' + Srundim + ']... .. .'
+            if stats == 'full':
+                minv = np.min(field)
+                maxv = np.max(field)
+                meanv = np.mean(field)
+                mean2v = np.mean(field*field)
+                varv = np.sqrt(mean2v - meanv*meanv)
+
+                if stdout:
+                    print '    Fstats shape:',field.shape
+                    print '    Fstats min:', minv
+                    print '    Fstats max:', maxv
+                    print '    Fstats mean:', meanv
+                    print '    Fstats mean2:', mean2v
+                    print '    Fstats variability:', varv
+
+                if countV is not None:
+                    for icV in range(NcountV):
+                        convals = counts[icV]
+                        if stdout:
+                            print "    Fstats N values '" + convals[1] + "'",        \
+                              convals[2], ':', convals[0]
+                        contvalsS = 'count(', + convals[1] + ',',convals[2],')'
+                else:
+                    convals = [int(0)]
+                    convalsS = ''
+                statsvariable[vn] = [minv, maxv, meanv, mean2v, varv, convals[0]]
+                if countV is not None:
+                    statsresults[tuple([slice(0,Nstats)]+list(slicedims))] = [minv,  \
+                      maxv, meanv, mean2v, varv, convals[0]]
+                else:
+                    statsresults[tuple([slice(0,Nstats)]+list(slicedims))] = [minv,  \
+                      maxv, meanv, mean2v, varv]
+
+            else:
+                print errormsg
+                print '  ' + fname + ": kind of statistics '" + stats + "' not ready!!"
+                print '    valid statistics: full'
+                quit(-1)
+
+            Sm = 'MAT:  '
+            if stdout:
+                print '  '
+                print Sm + 'Matrix of the statisitcs _______'
+                print Sm + ' '.join('{:<15}'.format(vn) for vn in statns)
+                for key in statsvariable.keys():
+                    vals = statsvariable[key]
+                    print Sm+'{:<15}'.format(key)+' '+' '.join('{:<15g}'.format(v) for v in vals)
+
+        newvar = onewnc.createVariable(vn + '_stats', 'f4', tuple(['stats'] +        \
+          runvdimns))
+        newvar[:] = statsresults[:]
+        for attrn in objfield.ncattrs():
+            attrv = objfield.getncattr(attrn)
+            newattr = set_attribute(newvar,attrn,attrv)
+        newattr = set_attribute(newvar,'running_stats',', '.join(runvdimns))
+        onewnc.sync()
+
+    # Global attributes
+    add_globattrs(ncobj,onewnc,'all')
+    onewnc.sync()
+
+    ncobj.close()
+    onewnc.close()
+
+    print fname + ": successfull written of file '" + ofile + "' !!"
+
+    return
+
+#field_stats_dim('full,None,None,west_east,XLONG,False', '/home/lluis/PY/wrfout_d01_2001-11-11_00:00:00', 'T2')
 
 def filter_2dim(values, ncfile, varn):
@@ -11758 +12097 @@
 
 #changevartype('i', 'test.nc', 'var')
-
-def add_dims(oc,nc,dnames):
-    """ Function to add dimensions from a given netCDF object to another one
-      oc: source netcdfile object
-      nc: netcdfile object to add dimensions
-      dnames: list of name of dimensions to add
-    """
-    fname = 'add_dims'
-
-    for dn in dnames:
-        print '  ' + fname + ': adding dim:',dn,' ...'
-        if oc.dimensions[dn].isunlimited():
-            newdim = nc.createDimension(dn, None)
-        else:
-            newdim = nc.createDimension(dn, len(oc.dimensions[dn]))
-
-    nc.sync()
-
-    return
-
-def add_vars(oc,nc,vnames):
-    """ Function to add variables from a given netCDF object to another one
-      oc: source netcdfile object
-      nc: netcdfile object to add dimensions
-      vnames: list of name of variables to add
-    """
-    fname = 'add_vars'
-
-    for vn in vnames:
-        if not gen.searchInlist(nc.variables,vn):
-            print '  ' + fname + ': adding var:',vn,' ...'
-
-            varo = oc.variables[vn]
-            vartype = varo.dtype
-            vardims = varo.dimensions
-            varattrs = varo.ncattrs()
-
-            for vdn in vardims:
-                if not gen.searchInlist(nc.dimensions,vdn):
-                    print warnmsg
-                    print '  ' + fname + ": adding dimension '" + vdn +              \
-                      "' from variable '" + vdn + "' which is not in file !!"
-                    add_dims(oc,nc,[vdn])
-
-            if gen.searchInlist(varattrs,'_FillValue'):
-                newvar = nc.createVariable(vn, vartype, vardims,                     \
-                  fill_value=varo.getncattr('_FillValue'))
-            else:
-                newvar = nc.createVariable(vn, vartype, vardims)
-
-            for attrn in varattrs:
-                attrv = varo.getncattr(attrn)
-                newattr = set_attribute(newvar, attrn, attrv)
-
-            newvar[:] = varo[:]
-
-        nc.sync()
-
-    return
-
-def add_globattrs(oc,nc,attrns):
-    """ Function to add global attributes from a given netCDF object to another one
-      oc: source netcdfile object
-      nc: netcdfile object to add dimensions
-      vnames: list of name of dimensions to add ('all', for all attributes)
-    """
-    fname = 'add_globattrs'
-
-    allfattrs = oc.ncattrs()
-
-    if attrns == 'all':
-        fattrns = allfattrs
-    else:
-        fattrns = attrns
-
-    for an in fattrns:
-        if not gen.searchInlist(allfattrs,an):
-            print errormsg
-            print '  ' + fname + ": file has not global attribute '" + an + "' !!"
-            quit(-1)
-
-        print '  ' + fname + ': adding attr:',an,' ...'
-        
-        av = oc.getncattr(an)
-        newattr = set_attribute(nc, an, av)
-
-    nc.sync()
-
-    return
 
 def selvar(values, ncfile, varn):
@@ -17214 +17464 @@
 
     ncvars = ncf.variables.keys()
-    print fname + '; Lluis ncvars:', ncvars
     olon2 = ncf.variables['lon2']
     dimslon2 = olon2.dimensions
@@ -17231 +17480 @@
     olat2 = ncf.variables['lat2']
     dimslat2 = olat2.dimensions
-    print fname + '; Lluis ncvars after 1st:', ncvars, 'dimslat2:', dimslat2
     newvar = ncf.createVariable('lat','f8',dimslat2)
     newvar[:] = olat2[:]