Changeset 1233 in lmdz_wrf for trunk

Timestamp:

Oct 26, 2016, 3:37:25 PM (9 years ago)

Author:

lfita

Message:

Improving `compute_opersvarsfiles' with more operations and options and clarity

Location:

trunk/tools

Files:

• generic_tools.py (modified) (2 diffs)
• nc_var.py (modified) (1 diff)
• nc_var_tools.py (modified) (3 diffs)

trunk/tools/generic_tools.py

@@ -9892 +9892 @@
 
     matAshape = list(matA.shape)
-    print fname + '; Lluis matAshape:', matAshape
 
     # From: https://en.wikipedia.org/wiki/Finite_difference_coefficient
@@ -10293 +10292 @@
 
     return newmat, Lopn[opn]
-print matoperations('centerderiv,1,8,0',np.arange(81.).reshape(9,9))
-print matoperations('forwrdderiv,2,2,0',np.arange(81.).reshape(9,9))
-print matoperations('upthres,4,3',np.arange(81).reshape(9,9))
-print matoperations('lowthres@oper,37,divc,7',np.arange(81).reshape(9,9))
 
 #quit()

trunk/tools/nc_var.py

@@ -31 +31 @@
 ## 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'
+## e.g. # nc_var.py -o compute_opersvarsfiles -S 'west_east|XLONG|-1;south_north|XLAT|-1;Time|Times|3@add|wrfout_d01_2001-11-11_00:00:00|T2%west_east|XLONG|-1;south_north|XLAT|-1;Time|Times|3@subc,273.15|wrfout_d01_2001-11-11_00:00:00|None' -v 'tempC,air!temperature,C'
 
 from optparse import OptionParser

trunk/tools/nc_var_tools.py

@@ -6882 +6882 @@
 #compute_deaccum('time', 'control/sellonlatbox_PRECIP_dayaccum.nc', 'PRECIP')
 
-def compute_opersvarsfiles(values, varname):
+def compute_opersvarsfiles(values, varinfo):
     """ Function to compute opersvarfiles: operation of variables from different files
       (OPER1.FILE1_VAR1 OPER2.FILE2_VAR2), operations are going to be secuentially 
       made
-        values= [dimvons]@[operfilevars]
-          [dimnvons]: [dimo1]|[dimo2]|[...|[dimoN]] '|' seaprated list of names with the 
-            variables which contain the values of the dimensions of the compute variable
-          operfilevars: [oper1]|[file1]|[varN],[...,[operK]|[fileM]|[varN]] ',' separated list of 
-             triplets of [operation], [file], [variable name to use]
-           * [oper]: operations: add,sub,mul,div,pot
-        varname= name to the final variable
+        values= '%' separated list of dimension ranges and file,operation,variable to compute
+          [dimranges1]@[operfile1var]%[dimranges2]@[operfile2var]%[...[dimrangesM]@[operfileMvar]]
+          dimranges: ';' separated list of dimension, dimension-variable names and their ranges ('|' separated)
+            [dim1]|[dimv1]|[range];[dim2]|[dimv2]|[range];[...|[dimN]|[dimvN]|[range]]
+            with [range] as:
+                       -1, all the values
+                       -9, last values
+                      int, a single value
+            [beg,end,frq], from a beginning to an end with a given frequency
+          operfilevar: [oper]|[file]|[var] '|' separated list of triplets of [operation], [file], [variable name to use]
+           * [oper]: operations (being [prevalues] the values computed until the operation)
+             'add': adding [var] ([prevalues] + [var])
+             'addc',[modval1]: add [modval1] to [prevalues]
+             'centerderiv',[N],[ord],[dim]: un-scaled center [N]-derivative of order [ord] along dimension [dim] of [var]
+             'div': dividing by [var] ([prevalues] / [var])
+             'divc',[modval1]: [prevalues] divide by [modval1]
+             'forwrdderiv',[N],[ord],[dim]: un-scaled forward [N]-derivative of order [ord] along dimension [dim] of [var]
+             'inv': inverting [prevalues] (1/[prevalues])
+             'lowthres',[modval1],[modval2]: if [prevalues] < [modval1]; prevalues = [modval2]
+             'lowthres@oper',[modval1],[oper],[modval2]: if [prevalues] < [modval1]; prevalues = [oper] (operation as [modval2])
+             'mul': multiplying by [var] ([prevvalues] * [var])
+             'mulc',[modval1]: [prevalues] multiplied by [modval1]
+             'pot': powering with [var] ([prevalues] ** [var])
+             'potc',[modval1]: [prevalues] ** [modval1]
+             'sub': substracting [var] ([prevalues] - [var])
+             'subc',[modval1]: remove [modval1] of [prevalues]
+             'upthres',[modval1],[modval2]: if [prevalues] > [modval1]; prevalues = [modval2]
+             'upthres@oper',[modval1],[oper],[modval2]: if [prevalues] > [modval1]; prevalues = [oper] (operation with [modval2])
+           * [file]: name of the file
+           * [var]: variable to use ('None' for the case with constant operations: 'addc', 'divc', 'inv', 'lowthres', 
+              'lowthres@oper', 'mulc', 'potc', 'subc', 'upthres', 'upthres@oper'
+        varinfo= [varname],[Lvarname],[varunits]
+          varname: name of the final variable 
+          Lvarname: long name of the final variable ('!' for spaces)
+          varunits: units of the final variable ('!' for spaces)
     """
-#    import numpy.ma as ma
+    import numpy.ma as ma
     fname='compute_opersvarsfiles'
 
@@ -6902 +6930 @@
         quit()
 
-    dimovars = values.split('@')[0].split('|')
-    filevars = values.split('@')[1].split(',')
+    # List of operations which need to be started with a matrix of zeros
+    opinizero = ['add', 'sub']
+    # List of operations which need to be started with a matrix of ones
+    opinione = ['mul', 'div', 'pot']
+
+    varname = varinfo.split(',')[0]
+    Lvarname = varinfo.split(',')[1].replace('!',' ')
+    varunits = varinfo.split(',')[2].replace('!',' ')
 
     ofile = 'opersvarsfiles_' + varname + '.nc'
 
-    Nfilevars = len(filevars)
-
-    operation = 'operation: '
+    fileopers = values.split('%')
+
+    # Getting range,files,operations, ...
+    Nfileopers = len(fileopers)
+    filens = []
+    ranges = []
+    opers = []
+    varns = []
+    for fileoper in fileopers:
+        ranges.append(fileoper.split('@')[0])
+        fnop = fileoper.split('@')[1]
+        opers.append(fnop.split('|')[0])
+        filens.append(fnop.split('|')[1])
+        varns.append(fnop.split('|')[2])
+
+    print '  ' + fname + 'Files/Operations to perform _______'
+    for ifop in range(Nfileopers):
+        print '    ' + filens[ifop] + ' range: ' + ranges[ifop] + ' operation: ' +   \
+          opers[ifop] + ' variable: ' + varns[ifop]
+
     newunits = '-'
     opervars = []
-
-    for ifv in range(Nfilevars):
-        opern = filevars[ifv].split('|')[0]
-        filen = filevars[ifv].split('|')[1]
-        varn = filevars[ifv].split('|')[2]
+    # Dictionary with the sizes of the dimensions of prevalues
+    Dds = {}
+    # Dictionary with the values of the dimension-variables of the dimensions of `prevalues'
+    Ddv = {}
+    # Dictionary with the dimensions of the dimension-variables of the dimensions of `prevalues'
+    Dddv = {}
+    # Dictionary with the type of the dimension-variables of the dimensions of `prevalues'
+    Dtdv = {}
+    # Dictionary with the attributes dictionay of the dimension-variables of the dimensions of `prevalues'
+    Dadv = {}
+
+    prevalues = None
+    for ifv in range(Nfileopers):
+        opern = opers[ifv]
+        filen = filens[ifv]
+        varn = varns[ifv]
+        rng = ranges[ifv].split(';')
+
         opervars.append(varn)
 
@@ -6926 +6990 @@
         objnc = NetCDFFile(filen, 'r')
 
-        if not objnc.variables.has_key(varn):
-            print errormsg
-            print '  ' + fname + ': netCDF file "' + filen +                         \
-              '" does not have variable "' + varn + '" !!!!'
-            quit(-1)
-
-        varobj = objnc.variables[varn]
-        varvals = varobj[:]
-
-# Creation of the operated variable
-        if ifv == 0:
-            vardims = varobj.dimensions
-            vartype = varobj.dtype
-
-            if opern == 'add' or opern == 'sub':
-                newvarv = np.zeros(tuple(varobj.shape), dtype=np.float)
-            elif opern == 'mul' or opern == 'div' or opern == 'pot':
-                newvarv = np.ones(tuple(varobj.shape), dtype=np.float)
+        # Operation with a variable
+        if varn != 'None' or prevalues is None:
+            if not objnc.variables.has_key(varn):
+                print errormsg
+                print '  ' + fname + ': netCDF file "' + filen +                     \
+                  '" does not have variable "' + varn + '" !!!!'
+                quit(-1)
+
+            varobj = objnc.variables[varn]
+            # Slice values
+            dictslice = {}
+
+            Ddimvs = []
+            for rg in rng:
+                dimn = rg.split('|')[0]
+                dimv = rg.split('|')[1]
+                dimr = rg.split('|')[2]
+                if dimr.find(',') != -1:
+                    dictslice[dimn] = list(np.array(dimr.split(','), dtype=int))
+                else:
+                    dictslice[dimn] = int(dimr)
+                # Dimension-variable
+                if not objnc.variables.has_key(dimv):
+                    print errormsg
+                    print '  ' + fname + ": file '" + filen + "' does not have " +   \
+                      "dimension-variable '" + dimv + "' !!"
+                    quit(-1)
+                # Keeping the dimensions
+                Ddimvs.append([dimn, dimv, dimr])
+
+            # Getting values of the variable
+            varslice, dimvarvals = SliceVarDict(varobj, dictslice)
+            varvals = varobj[tuple(varslice)]
+
+            # Getting values of the correspondant dimensions of the variable
+            for iddd in range(len(Ddimvs)):
+                DD = Ddimvs[iddd]
+                dimn = DD[0]
+                dimv = DD[1]
+                dimr = DD[2]
+                if dimv != 'WRFt':
+                    odimv = objnc.variables[dimv]
+                    varslice, Dddv[dimn] = SliceVarDict(odimv, dictslice)
+                    Ddv[dimn] = odimv[tuple(varslice)]
+                    Dtdv[dimn] = odimv.dtype
+                    dicattrs = {}
+                    for attrn in odimv.ncattrs():
+                        attrval = odimv.getncattr(attrn)
+                        dicattrs[attrn] = attrval
+                    Dadv[dimn] = dicattrs
+                else:
+                    datesv = []
+                    wrft = objnc.variables['Times']
+                    slicewrft, dwrfd = SliceVarDict(wrft, dictslice)
+                    fdates = wrft[tuple(slicewrft)]
+                    refdate = '19491201000000'
+                    tunitsval = 'minutes'
+                    yrref=refdate[0:4]
+                    monref=refdate[4:6]
+                    dayref=refdate[6:8]
+                    horref=refdate[8:10]
+                    minref=refdate[10:12]
+                    secref=refdate[12:14]
+
+                    refdateS = yrref + '-' + monref + '-' + dayref + ' ' + horref +  \
+                      ':' + minref + ':' + secref
+
+                    dt = fdates.shape[0]
+                    cftimes = np.zeros((dt), dtype=np.float)
+                    for it in range(dt):
+                        wrfdates = gen.datetimeStr_conversion(timewrfv[it,:],        \
+                          'WRFdatetime', 'matYmdHMS')
+                        cftimes[it] = gen.realdatetime1_CFcompilant(wrfdates,        \
+                          refdate, tunitsval)
+                    tunits = tunitsval + ' since ' + refdateS
+                    Dddv[dimn] = 'Time'
+                    Ddv[dimn] = cftimes
+                    Dtdv[dimn] = np.float
+                    Dadv[dimn] = {'stadard_name': 'time', 'long_name': 'time',       \
+                      'units': tunitsval + ' since ' + refdateS}
+
+            # Getting sizes of the dimensions
+            for dvn in dimvarvals:
+                dimobj = objnc.dimensions[dvn]
+                if dimobj.isunlimited():
+                    Dds[dvn] = None
+                else:
+                    Dds[dvn] = len(dimobj)
+
+            # operation
+            if prevalues is None:
+                # Creation of an empty variable to start with
+                if gen.searchInlist(opinizero,opern):
+                    inimat = np.zeros(varvals.shape, dtype=varobj.dtype)
+                elif gen.searchInlist(opinione,opern):
+                    inimat = np.ones(varvals.shape, dtype=varobj.dtype)
+
+                prevalues, operS = gen.matoperations(opern, varvals, inimat)
             else:
-                print errormsg
-                print '  ' + fname + ': operation "' + opern + '" is not ready !!!!'
-                quit(-1)
-
-        if gen.searchInlist(varobj.ncattrs(), '_FillValue'):
-            NOvalue = varobj.getncattr('_FillValue')
-            print warnmsg
-            print '    ' + fname + ': variable with a missing value:',NOvalue
-
-        if opern == 'add':
-            if gen.searchInlist(varobj.ncattrs(), '_FillValue'):
-                if type(varvals) != type(np.arange(1)): 
-                    print type(varvals), type(np.arange(1))
-                    prevals = varvals.filled(0.)
-                else:
-                    prevals = varvals
-                newvarv = newvarv + prevals
-            else:
-                newvarv = newvarv + varvals
-            operation = operation + ' +' + varn
-            newunits = gen.variables_values(varname)[5]
-        elif opern == 'sub':
-            if gen.searchInlist(varobj.ncattrs(), '_FillValue'):
-                if type(varvals) != type(np.arange(1)): 
-                    prevals = varvals.filled(0.)
-                else:
-                    prevals = varvals
-                newvarv = newvarv - prevals
-            else:
-                newvarv = newvarv - varvals
-            operation = operation + ' -' + varn
-            newunits = gen.variables_values(varname)[5]
-        elif opern == 'mul':
-            if gen.searchInlist(varobj.ncattrs(), '_FillValue'):
-                prevals = varvals.filled(1.)
-                newvarv = newvarv * prevals
-            else:
-                newvarv = newvarv * varvals
-            operation = operation + ' *' + varn
-            newunits = gen.variables_values(varname)[5]
-        elif opern == 'div':
-            if gen.searchInlist(varobj.ncattrs(), '_FillValue'):
-                prevals = varvals.filled(1.)
-                newvarv = newvarv / prevals
-            else:
-                newvarv = newvarv / varvals
-            operation = operation + ' /' + varn
-            newunits = gen.variables_values(varname)[5]
-        elif opern == 'pot':
-            if gen.searchInlist(varobj.ncattrs(), '_FillValue'):
-                prevals = varvals.filled(1.)
-                newvarv = newvarv ** prevals
-            else:
-                newvarv = newvarv ** varvals
-            operation = operation + ' **' + varn
-            newunits = gen.variables_values(varname)[5]
+                prevalues, operS = gen.matoperations(opern, prevalues, varvals)
+            operS = ' ' + operS + '[' + varn + ']'
+        else:
+            prevalues, operS = gen.matoperations(opern, prevalues)
+            operS = ' ' + operS
+
+        if opern == opers[0]:
+            operationS = operS 
+        else:
+            operationS = operationS + operS
+
+        objnc.close()
 
 # Creation of the new file
     objofile = NetCDFFile(ofile, 'w')
 
-# Dimensions
-    for dimn in vardims:
-        dimobj = objnc.dimensions[dimn]
-        if dimobj.isunlimited():
-            print '      ' + fname + ': Unlimited dimension '
-            dimsize = None
-        else:
-            dimsize = len(dimobj)
-
-        newdim = objofile.createDimension(dimn, dimsize)
-
-# Variables for dimensions
-    for varn in dimovars:
-
-        if not gen.searchInlist(list(objnc.variables.keys()), varn):
-            print warnmsg
-            print '  ' + fname + " file does not have dimension variable '" + varn + \
-              "' skipping it !!"
-        else:
-            dvarobj = objnc.variables[varn]
-            # Adding that dimensions of the variable which might not be in the output 
-            for vdn in dvarobj.dimensions:
-                if not objofile.dimensions.has_key(vdn):
-                    addvdno = objnc.dimensions[vdn]
-                    if addvdno.isunlimited():
-                        objofile.createDimension(vdn, None)
-                    else:
-                        objofile.createDimension(vdn, len(addvdno))
-            # Adding new variable to output
-            newvar = objofile.createVariable(varn, dvarobj.dtype, dvarobj.dimensions)
-            for attrn in  dvarobj.ncattrs():
-                attrval = dvarobj.getncattr(attrn)
-                newattr = newvar.setncattr(attrn, attrval)
-            newvar[:] = dvarobj[:]
-
-    newvar = objofile.createVariable(varname, vartype, vardims)
-    newvar[:] = newvarv
-    varattrs = varobj.ncattrs()
-    if gen.searchInlist(varattrs,'long_name'): 
-        operation = varobj.getncattr('long_name') + '; ' + operation
-
-    newattr = basicvardef(newvar, varname, operation, newunits)
-    stdattrs = ['standard_name', 'long_name', 'units']
-    for attrn in varattrs:
-        if not gen.searchInlist(stdattrs, attrn) and attrn != '_FillValue':
-            attrv = varobj.getncattr(attrn)
-            newattr = set_attribute(newvar,attrn,attrv)
+    # Dimensions
+    for dimn in dimvarvals:
+        newdim = objofile.createDimension(dimn, Dds[dimn])
+    objofile.sync()
+
+    # Variables for dimensions
+    for dvn in dimvarvals:
+        for idn in Dddv[dvn]:
+            if not objofile.dimensions.has_key(idn): objofile.createDimension(idn,   \
+              Dds[idn])
+        newvar = objofile.createVariable(dvn, Dtdv[dvn], tuple(Dddv[dvn]))
+        dictatr = Dadv[dvn]
+        for attrn in dictatr.keys():
+            newattr = newvar.setncattr(attrn, dictatr[attrn])
+        newvar[:] = Ddv[dvn]
+    objofile.sync()
+
+    # Computed variable
+    vartype = prevalues.dtype
+    if type(prevalues) == type(ma.asarray([1,1])):
+        maskvalue = prevalues.fill_value
+        newvar = objofile.createVariable(varname, nctype(vartype), dimvarvals,       \
+          fill_valu=maskvalue)
+    else:
+        newvar = objofile.createVariable(varname, nctype(vartype), dimvarvals)
+    newvar[:] = prevalues[:]
+
+    newattr = basicvardef(newvar, varname, Lvarname, varunits)
 
 # Global attributes
-    for attrn in objnc.ncattrs():
-        attrval = objnc.getncattr(attrn)
-        newattr = objofile.setncattr(attrn, attrval)
-    
-    objnc.close()
+    objofile.setncattr('script',  fname)
+    objofile.setncattr('version',  '1.0')
+    objofile.setncattr('author',  'L. Fita')
+    newattr = set_attributek(objofile, 'institution', unicode('Laboratoire de M' +   \
+      unichr(233) + 't' + unichr(233) + 'orologie Dynamique'), 'U')
+    objofile.setncattr('university',  'Pierre et Marie Curie')
+    objofile.setncattr('country',  'France')
+    newattr = set_attributek(objofile, 'used_operation',  operationS, 'S')
+    newattr = set_attributek(objofile, 'used_files',  ', '.join(filens) ,'S')
+    newattr = set_attributek(objofile, 'used_ranges',  ', '.join(ranges) ,'S')
 
     objofile.sync()