Changeset 662 in lmdz_wrf for trunk

Timestamp:

Jan 8, 2016, 6:08:34 PM (9 years ago)

Author:

lfita

Message:

Adding `SpatialWeightedMean?' to compute spatial averaged means

Location:

trunk/tools

Files:

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

trunk/tools/nc_var.py

@@ -36 +36 @@
   'remapnn',                                                                         \
   'seasmean', 'sellonlatbox', 'sellonlatlevbox', 'selvar', 'setvar_asciivalues',     \
-  'sorttimesmat', 'spacemean', 'statcompare_files', 'submns', 'subyrs', 'TimeInf',   \
+  'sorttimesmat', 'spacemean', 'SpatialWeightedMean', 'statcompare_files', 'submns', \
+  'subyrs', 'TimeInf',                                                               \
   'TimeSplitmat', 'timemean', 'valmod', 'valmod_dim','varaddattrk', 'varaddattr',    \
   'varaddref',                                                                       \
@@ -220 +221 @@
 elif oper == 'spacemean':
     ncvar.spacemean(opts.ncfile, opts.varname)
+elif oper == 'SpatialWeightedMean':
+    ncvar.SpatialWeightedMean(opts.values, opts.ncfile, opts.varname)
 elif oper == 'statcompare_files':
     ncvar.statcompare_files(opts.values)

trunk/tools/nc_var_tools.py

@@ -17110 +17110 @@
     return
 
-
+def DimsLoop(ovar,seldims):
+    """ Function to provide the shape of the dimensions which are not selected
+      ova= variable object
+      seldims= list of dimensions to select
+    """
+    fname = 'DimsLoop'
+
+    newshape = []
+    noloopdims = []
+
+    idim = 0
+    for dim in list(ovar.dimensions):
+        if not searchInlist(seldims,dim):
+            newshape.append(ovar.shape[idim])
+            noloopdims.append(dim)
+        idim = idim + 1
+
+    return newshape, noloopdims
+
+def SliceVar(ovar,dloop,dloopindex):
+    """ Function to slice a given variable throughout a given list of dimensions
+      ovar= variable object to slice
+      dloop= list of dimensions to get a single value
+      dloopindex= value of the loop dimenion to get
+    """
+    fname = 'SliceVar'
+
+    varshape = ovar.shape
+
+    slicevals = []
+    idim = 0
+    for dim in ovar.dimensions:
+        if searchInlist(dloop,dim):
+            slicevals.append(dloopindex[idim])
+        else:
+            slicevals.append(slice(0,varshape[idim]))
+        idim = idim + 1
+
+    return slicevals
+
+def SpatialWeightedMean(values, filen, varn):
+    """ Function to compute the spatial mean using weights from a netCDF file
+      values= [weightskind],[xdimname],[ydimname]
+        [weightskind]: type of weights:
+          * 'varnames',[varname],[oper]: using a variable [varname] with an operation [oper] as area size
+            [oper]: available operations
+              'inv': inverse of the values 1/[varname]
+              'nothing': direct values
+          * 'reglonlat',[lonname],[latname]: regular lon/lat projection to compute the area size
+        [xdimname]: name of the dimension for the x-axis
+        [ydimname]: name of the dimension for the y-axis
+      filen= name of the netCDF file
+      varn= name of the variable
+    """
+    fname = 'SpatialWeightedMean'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print SpatialWeightedMean.__doc__
+        quit()
+
+    operations = ['inv','nothing']
+
+    weightk = values.split(',')[0]
+    xdimname = values.split(',')[3]
+    ydimname = values.split(',')[4]
+
+    ofile = 'SpatialWeightedMean_' + weightk + '.nc'
+
+    if weightk == 'varnames':
+        arguments = '[weightk],[varname],[oper],[xdimname],[ydimname]'
+        check_arguments(fname, len(values.split(',')), arguments, ',',               \
+          len(arguments.split(',')))
+        varname = values.split(',')[1]
+        oper = values.split(',')[2]     
+        if not searchInlist(operations, oper):
+            print errormsg
+            print '  ' + fname + ": operation '" + oper + "' not ready!!"
+            print '  ready operations:',operations
+            quit(-1)
+    elif weightk == 'reglonlat':
+        arguments = '[weightk],[lonname],[latname],[xdimname],[ydimname]'
+        check_arguments(fname, len(values.split(',')), arguments, ',',               \
+          len(arguments.split(',')))
+        lonname = values.split(',')[1]
+        latname = values.split(',')[2]
+    else:
+        print errormsg
+        print '  ' + fname + ": wieght type '" + weightk + "' not ready!!"
+        quit(-1)
+
+# Operations
+##
+    onc = NetCDFFile(filen, 'r')
+    if not onc.variables.has_key(varn):
+        print errormsg
+        print '  ' + fname + ": file '" + filen + "' does not have variable '" +     \
+          varn + "' !!"
+        quit(-1)
+    iovar = onc.variables[varn]
+    if not searchInlist(list(iovar.dimensions),xdimname) and not                     \
+      searchInlist(list(iovar.dimensions),ydimname):
+        print errormsg
+        print '  ' + fname + ": variable '" + varn + "' does not have dimension '" + \
+          xdimname + "' neither '" + ydimname + "' !!"
+        quit(-1)
+
+# Getting shape of output variables for future average loop
+##
+    loopshape, dimsloop = DimsLoop(iovar,[xdimname, ydimname])
+
+    if len(loopshape) > 3:
+        print erromsg
+        print '  ' + fname + ': output size not ready!!'
+        quit(-1)
+
+    ivarv = iovar[:]
+
+    if weightk == 'varnames':
+        if not onc.variables.has_key(varname):
+            print errormsg
+            print '  ' + fname + ": file '" + filen + "' does not have variable '" + \
+              varname + "' !!"
+            quit(-1)
+        iovarwgt = onc.variables[varname]
+
+        ivarwgtv = iovarwgt[:]
+        if oper == 'inv':
+            longvarname = 'using inverse of variable ' + varname + ' as space weights'
+            if len(loopshape) == 1:
+                newvals = np.ones((loopshape[0]), dtype=np.float)*fillValueF
+                for id1 in range(loopshape[0]):
+                    slicevalues = SliceVar(iovar,dimsloop,[id1])
+                    newvals[id1] = np.mean(ivarv[tuple(slicevalues)] /               \
+                      ivarwgtv[tuple(slicevalues)])
+            elif len(loopshape) == 2:
+                newvals = np.ones((loopshape[0],loopshape[1]),dtype=np.float)*       \
+                  fillValueF
+                for id1 in range(loopshape[0]):
+                    for id2 in range(loopshape[1]):
+                        slicevalues = SliceVar(iovar,dimsloop,[id1,id2])
+                        newvals[id1,id2] = np.mean(ivarv[tuple(slicevalues)] /       \
+                          ivarwgtv[tuple(slicevalues)])
+            elif len(loopshape) == 3:
+                newvals = np.ones((loopshape[0],loopshape[1],loopshape[2]),          \
+                  dtype=np.float)*fillValueF
+                for id1 in range(loopshape[0]):
+                    for id2 in range(loopshape[1]):
+                        for id3 in range(loopshape[1]):
+                            slicevalues = SliceVar(iovar,dimsloop,[id1,id2,id3])
+                            newvals[id1,id2,id3]= np.mean(ivarv[tuple(slicevalues)]/ \
+                              ivarwgtv[tuple(slicevalues)])
+        elif oper == 'nothing':
+            longvarname = 'using variable ' + varname + ' as space weights'
+            if len(loopshape) == 1:
+                newvals = np.ones((loopshape[0]), dtype=np.float)*fillValueF
+                for id1 in range(loopshape[0]):
+                    slicevalues = SliceVar(iovar,dimsloop,[id1])
+                    newvals[id1] = np.mean(ivarv[tuple(slicevalues)]*                \
+                      ivarwgtv[tuple(slicevalues)])
+            elif len(loopshape) == 2:
+                newvals = np.ones((loopshape[0],loopshape[1]), dtype=np.float)*      \
+                  fillValueF
+                for id1 in range(loopshape[0]):
+                    for id2 in range(loopshape[1]):
+                        slicevalues = SliceVar(iovar,dimsloop,[id1,id2])
+                        newvals[id1,id2] = np.mean(ivarv[tuple(slicevalues)]*        \
+                          ivarwgtv[tuple(slicevalues)])
+            elif len(loopshape) == 3:
+                newvals = np.ones((loopshape[0],loopshape[1],loopshape[2]),          \
+                  dtype=np.float)*fillValueF
+                for id1 in range(loopshape[0]):
+                    for id2 in range(loopshape[1]):
+                        for id3 in range(loopshape[2]):
+                            slicevalues = SliceVar(iovar,dimsloop,[id1,id2.id3])
+                            newvals[id1,id2,id3]=np.mean(ivarv[tuple(slicevalues)]*  \
+                              ivarwgtv[tuple(slicevalues)])
+    elif weightk == 'reglonlat':
+        longvarname = 'using cos(latitude) variable ' + latname + ' as space weights'
+        if not onc.variables.has_key(lonname):
+            print errormsg
+            print '  ' + fname + ": file '" + filen + "' does not have variable '" + \
+              lonname + "' !!"
+            quit(-1)
+        if not onc.variables.has_key(latname):
+            print errormsg
+            print '  ' + fname + ": file '" + filen + "' does not have variable '" + \
+              latname + "' !!"
+            quit(-1)
+
+        iolon = onc.variables[lonname]
+        iolat = onc.variables[latname]
+        ilonv = iolon[:]
+        ilatv = iolat[:]
+
+        lonv, latv = lonlat2D(ilonv, ilatv)
+
+        ivarwgtv = np.abs(np.cos(latv*np.pi/180.))
+        if len(loopshape) == 1:
+            newvals = np.ones((loopshape[0]), dtype=np.float)*fillValueF
+            for id1 in range(loopshape[0]):
+                slicevalues = SliceVar(iovar,dimsloop,[id1])
+                newvals[id1] = np.mean(ivarv[tuple(slicevalues)] * ivarwgtv)
+        elif len(loopshape) == 2:
+            newvals = np.ones((loopshape[0],loopshape[1]), dtype=np.float)*fillValueF
+            for id1 in range(loopshape[0]):
+                for id2 in range(loopshape[1]):
+                    slicevalues = SliceVar(iovar,dimsloop,[id1,id2])
+                    newvals[id1,id2] = np.mean(ivarv[tuple(slicevalues)] * ivarwgtv)
+        elif len(loopshape) == 3:
+            newvals = np.ones((loopshape[0],loopshape[1],loopshape[2]),              \
+              dtype=np.float)*fillValueF
+            for id1 in range(loopshape[0]):
+                for id2 in range(loopshape[1]):
+                    for id3 in range(loopshape[2]):
+                        slicevalues = SliceVar(iovar,dimsloop,[id1,id2,id3])
+                        newvals[id1,id2,id3] = np.mean(ivarv[tuple(slicevalues)] *   \
+                          ivarwgtv)
+
+# Writting output file
+##
+    onewnc = NetCDFFile(ofile, 'w')
+
+# Dimensions creation
+##
+    for dim in dimsloop:
+        newdim = onewnc.createDimension(dim, len(onc.dimensions[dim]))
+
+# Output variable
+##
+    newvar = onewnc.createVariable(varn + 'spaceweightmean', 'f4', tuple(dimsloop),  \
+      fill_value=fillValueF)
+    basicvardef(newvar, varn + 'spaceweightmean', 'space weighted ' + varn + ' ' +   \
+      longvarname, iovar.getncattr('units'))
+    newvar[:] = newvals
+
+# Global attributes
+##
+    onewnc.setncattr('script',  fname)
+    onewnc.setncattr('version',  '1.0')
+    onewnc.setncattr('author',  'L. Fita')
+    onewnc.setncattr('institution',  'Laboratoire de Meteorology Dynamique')
+    onewnc.setncattr('university',  'Pierre et Marie Curie')
+    onewnc.setncattr('country',  'France')
+    for attrs in onc.ncattrs():
+        attrv = onc.getncattr(attrs)
+        attr = set_attribute(onewnc, attrs, attrv)
+
+    onc.close()
+
+    print fname + ": Successfull written of file: '" + ofile + "' !!"
+
+    return
+
+#SpatialWeightedMean('varnames,MAPFAC_M,inv,west_east,south_north','wrfout_d01_2001-11-11_00:00:00','HGT')
+#SpatialWeightedMean('reglonlat,XLONG,XLAT,west_east,south_north','wrfout_d01_2001-11-11_00:00:00','HGT')
 #quit()