Changeset 2280 in lmdz_wrf for trunk/tools

Timestamp:

Jan 22, 2019, 7:49:47 PM (6 years ago)

Author:

lfita

Message:

Working on:

• `compute_slices_stats_areaweighted': Function to compute stats of variables of a file splitting variables by slices along sets of ranges for a series of variables weighting by the area covered by each grid (defined as polygon) within the slice

File:

• trunk/tools/nc_var_tools.py (modified) (1 diff)

trunk/tools/nc_var_tools.py

@@ -27690 +27690 @@
 #area_weighted('yes:min,max,mean,stddev',fvals,'ct_values')
 
+def compute_slices_stats_areaweighted(values, ncfile, variable):
+    """ Function to compute stats of variables of a file splitting variables by 
+        slices along sets of ranges for a series of variables weighting by the area
+        covered by each grid (defined as polygon) within the slice
+      values=[slicevarsinf]@[dimvars]@[slicestatsdim]
+        [slicevarsinf] ';' separated list of [varn1],[minvar1],[maxvar1],
+          [slcevar1]; ...;[varnN],[minvarN],[maxvarN],[slcevarN]
+          varn[i]: name of the variable i
+          minvar[i]: minimum value to start the slices for variable i
+          maxvar[i]: maximum value to end the slices for variable i
+          slcevar[i]: length of slices in variable i
+            NOTE: slices will be:
+               [minvar[i]-slicevar[i]/2., minvar[i]+slicevar[i]/2.), ..., 
+               [maxvar[i]-slicevar[i]/2., maxvar[i]+slicevar[i]/2.)
+        [dimvars]: ':' separated list of [dimn]|[vardimn] dimension name [dimn] and 
+          its associated dimension-variable to get its values to provide values for 
+          the final file. NOTE: vardimn must have 'bounds' attribute in order to 
+          be able to compute the area-weights or it will be assumed that dimension
+          has no extension
+        [slicestatsdim]: ',' list of name of dimensions to do not use for computing 
+          statistics of variables at each slice ('any' for not leaving any 
+          dimension out, resulting on scalar statistics)
+      ncfile= netCDF file to use
+      variable: ',' list of variables ('all' for all variables)
+    """
+    fname = 'compute_slices_stats_areaweighted'
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print compute_slices_stats.__doc__
+        quit()
+
+    expectargs = '[slicevarsinf]@[dimvars]@[slicestatsdim]'
+    gen.check_arguments(fname, values, expectargs, '@')
+
+    varvalues0 = values.split('@')[0]
+    dimvars0 = values.split('@')[1]
+    slicestatsdim = values.split('@')[2]
+
+    varvalues = varvalues0.split(';')
+
+    varvals = {}
+    slcvarns = []
+    for varvls in varvalues:
+        expectargs = '[varni],[minvari],[maxvari],[slcevari]'
+        gen.check_arguments(fname, varvls, expectargs, ',')
+
+        varni = varvls.split(',')[0]
+        minvari = np.float(varvls.split(',')[1])
+        maxvari = np.float(varvls.split(',')[2])
+        slicevari = np.float(varvls.split(',')[3])
+        varvals[varni] = [minvari, maxvari, slicevari]
+        slcvarns.append(varni)
+    Nslcvars = len(slcvarns)
+
+    dimvars = gen.stringS_dictvar(dimvars0, Dc=':', DVs='|')
+    
+    onc = NetCDFFile(ncfile, 'r')
+
+    # Varibales to slice
+    if variable == 'all':
+        varns = onc.variables.keys()
+    else:
+        varns = gen.str_list(variable, ',')
+
+    # Getting dimensions
+    for dn in dimvars.keys():
+        if not gen.searchInlist(onc.dimensions,dn):
+            print errormsg
+            print '  ' +fname+ ": file '" + ncfile + "' does not have dimension: '"+ \
+              dn + "' !!"
+            dimns = list(onc.dimensions)
+            dimns.sort()
+            print '    available ones:', dimns
+            quit(-1)
+
+    slcvarnS = ''
+    slicesinf = {}
+    # Dictionary with the dimension-variables with bounds
+    vardimbndsvars = {}
+    # Dictionary with the result of slicing a dimension-variable with bounds
+    vardimbndslice = {}
+
+    # Preparing slices and output file
+    for varn in slcvarns:
+        if not onc.variables.has_key(varn):
+            print errormsg
+            print '  ' + fname + ": file '" + ncfile + "' does not have variable " + \
+              "i: '" + varn + "' !!"
+            varns = onc.variables.keys()
+            varns.sort()
+            print '    available ones:', varns
+            quit(-1)
+
+        ovar = onc.variables[varn]
+        vu = ovar.units
+        dimnv = list(ovar.dimensions)
+          
+        vvls = varvals[varn]
+        minvar = vvls[0]
+        maxvar = vvls[1]
+        slicevar = vvls[2]
+
+        dvar = (maxvar-minvar+slicevar)/slicevar
+        slices = np.arange(minvar, maxvar+slicevar, slicevar)
+        Nslices = slices.shape[0]
+
+        print '    ' + fname + ': slices ____'
+        print '      var:', varn, ':', Nslices, '(', minvar, ',', maxvar+slicevar,   \
+          ',', slicevar, ')'
+
+        # Slices var
+        sliceshape = [Nslices] + list(ovar.shape)
+        slcvar = np.zeros(tuple(sliceshape), dtype=bool)
+        slcvalsc = np.zeros((Nslices), dtype=np.float)
+        slcvals = np.zeros((Nslices,2), dtype=np.float)
+        varv = ovar[:]
+        for islc in range(Nslices):
+            slcvalsc[islc] = slices[islc]
+            slcvals[islc,0] = slices[islc]-slicevar/2.
+            slcvals[islc,1] = slices[islc]+slicevar/2.
+
+            #slcvar[islc,] = ma.masked_inside(varv, slcvals[islc,0], slcvals[islc,1]).mask
+            slcvar[islc,] = ma.masked_outside(varv, slcvals[islc,0], slcvals[islc,1]).mask
+
+        slicesinf[varn] = [Nslices, dimnv, slcvar, slcvalsc, slcvals]
+
+        if varn == slcvarns[0]:
+            onewnc = NetCDFFile(fname + '.nc', 'w')
+            newdim = onewnc.createDimension('slice_bnds', 2)
+            slcvarnS = varn
+        else:
+            if varn == slcvarns[Nslcvars-1]: slcvarnS = slcvarnS + ' & ' + varn
+            else: slcvarnS = slcvarnS + ', ' + varn
+
+        # dimensions
+        newdim = onewnc.createDimension('slice_'+varn, Nslices)
+
+        # variable dimensions
+        for dn in dimnv:
+            if not gen.searchInlist(onewnc.dimensions,dn): add_dims(onc, onewnc, [dn])
+
+        newvar = onewnc.createVariable('slice_'+varn, 'f', ('slice_'+varn))
+        newvar[:] = slcvalsc[:]
+        basicvardef(newvar, 'slice_'+varn, 'slices for variable ' + varn, vu)
+
+        newvar = onewnc.createVariable('slice_'+varn+'_bnds', 'f', ('slice_'+varn,  \
+          'slice_bnds'))
+        newvar[:] = slcvals[:]
+        basicvardef(newvar, 'slice_'+varn+'_bnds', 'boundaries of slices for ' +     \
+          'variable ' + varn, vu)
+
+        slcdims = ['slice_'+varn] + dimnv
+        slcshape = [Nslices] + list(ovar.shape)
+        slcv = np.zeros(tuple(slcshape), dtype=int)
+        for islc in range(Nslices):
+            vvslc = ~slcvar[islc,]
+            slcv[islc,] = np.where(vvslc, 1, 0)
+        newvar = onewnc.createVariable(varn+'sliced', 'i', tuple(slcdims))
+        newvar[:] = slcv[:]
+        basicvardef(newvar, varn+'sliced', 'sliced variable ' + varn, vu)
+
+        onewnc.sync()
+
+        # Looking for boundaries
+        varbnds = []
+        dn = varn
+        vdarattrs = ovar.ncattrs()
+        dv = ovar[:]
+
+        if gen.searchInlist(vdarattrs,'bounds'):
+            boundsv = ovar.getncattr('bounds')
+            print '    ' + infmsg
+            print '    ' +fname+ ": slicing variable '" +dn+ "' with bounds !!"
+            print '      bounds found:', boundsv
+            print '      getting them to retrieve the slices'
+            bndvarns = boundsv
+            varbnds.append(dn)
+            ovarbnds = onc.variables[boundsv]
+            varbnds = ovarbnds[:]
+
+            # Slicing following boundaries, dimensions are 1D thus expand to 2D
+            #    filling as NW, NE, SE, SW
+            if gen.searchInlist(slcvarns, dn) and not vardimbndslice.has_key(dn):
+                print infmsg
+                print '    ' + fname + ": slicing dimension '" + dn + "' ... "
+                varslcv = slicesinf[dn]
+                Nslices = varslcv[0]
+                ovardims = varslcv[1]
+                slcvar = varslcv[2]
+
+                if len(dv.shape) == 1:
+                    # Shapping 2D slices
+                    reflon = np.zeros((1,Nslices), dtype=np.float)
+                    reflon[0,:] = slcvalsc
+                    reflat = np.zeros((1,Nslices), dtype=np.float)
+                    xslice2D = np.zeros((4,1,Nslices), dtype=np.float)
+                    xslice2D[0,0,:] = slcvals[:,1]
+                    xslice2D[1,0,:] = slcvals[:,1]
+                    xslice2D[2,0,:] = slcvals[:,0]
+                    xslice2D[3,0,:] = slcvals[:,0]
+                    yslice2D = np.zeros((4,1,Nslices), dtype=np.float)
+                    yslice2D[0,0,:] = -1.
+                    yslice2D[1,0,:] = 1.
+                    yslice2D[2,0,:] = 1.
+                    yslice2D[3,0,:] = -1.
+    
+                    dd = len(onc.dimensions[dn])
+                    getlon = np.zeros((1,dd), dtype=np.float)
+                    getlon[0,:] = dv
+                    getlat = np.zeros((1,dd), dtype=np.float)
+                    xdimvarbnds2D = np.zeros((4,1,dd), dtype=np.float)
+                    xdimvarbnds2D[0,0,:] = ovarbnds[:,1]
+                    xdimvarbnds2D[1,0,:] = ovarbnds[:,1]
+                    xdimvarbnds2D[2,0,:] = ovarbnds[:,0]
+                    xdimvarbnds2D[3,0,:] = ovarbnds[:,0]
+                    ydimvarbnds2D = np.zeros((4,1,dd), dtype=np.float)
+                    ydimvarbnds2D[0,0,:] = -1.
+                    ydimvarbnds2D[1,0,:] = 1.
+                    ydimvarbnds2D[2,0,:] = 1.
+                    ydimvarbnds2D[3,0,:] = -1.
+
+                    dxref = Nslices
+                    dyref = 1
+                    dxget = dd
+                    dyget = 1
+
+                else:
+                    print errormsg
+                    print '  '+fname+": rank ", dv.shape, " of slicing varibale '" + \
+                      dn + "' not ready !!"
+                    print '    function only works with 1D slicing variables'
+                    quit(-1)
+
+                reflont = reflon.transpose()
+                reflatt = reflat.transpose()
+                refvarxbndst = xslice2D.transpose()
+                refvarybndst = yslice2D.transpose()
+                getlont = getlon.transpose()
+                getlatt = getlat.transpose()
+                getvarxbndst = xdimvarbnds2D.transpose()
+                getvarybndst = ydimvarbnds2D.transpose()
+
+                Ngridsint, gridsint, percenst =                                      \
+                  fsci.module_scientific.spacepercen(xcavals=reflont,                \
+                  ycavals=reflatt, xbavals=refvarxbndst, ybavals=refvarybndst,       \
+                  xcbvals=getlont, ycbvals=getlatt, xbbvals=getvarxbndst,            \
+                  ybbvals=getvarybndst, strict=False, dxa=dxref, dya=dyref,          \
+                  navertexmax=4, dxb=dxget, dyb=dyget, dxyb=dxget*dyget,             \
+                  nbvertexmax=4)
+
+                Ngridsin = Ngridsint.transpose()
+                gridsin = gridsint.transpose()
+                percens = percenst.transpose()
+
+                Ngridsinmax = np.max(Ngridsin)
+        else:
+            print infmsg
+            print '  ' + fname + ": slicing variable '" + dn + "' without bounds !!"
+            print '    directly using grid point values'
+
+            Ngridsin = np.zeros((1,Nslices), dtype=int)           
+            for islc in range(Nslices):
+                Ngridsin[0,islc] = np.sum(~slcvar[islc,])
+
+            Ngridsinmax = np.max(Ngridsin)
+            gridsin = np.zeros((2,Ngridsinmax,1,Nslices), dtype=int)
+            percens = np.zeros((Ngridsinmax,1,Nslices), dtype=np.float)
+            for islc in range(Nslices):
+                unmasked = slcvar[islc,]
+                indices = gen.multi_index_mat(unmasked,False)
+                ainds = np.array(indices)
+                if Ngridsin[0,islc] > 0:
+                    gridsin[0,0:Ngridsin[0,islc],0,islc] = ainds[:,0]
+                    gridsin[1,0:Ngridsin[0,islc],0,islc] = ainds[:,1]
+                    percens[0:Ngridsin[0,islc],0,islc] = 1.
+
+        # Let's avoid memory issues...
+        lvalues = [Ngridsin, gridsin, percens]
+        vardimbndslice[dn] = lvalues
+
+        # Writting it to the file variables space-weight
+        if not gen.searchInlist(onewnc.dimensions, dn+'slice'):
+            newdim = onewnc.createDimension(dn+'slice',Nslices)
+            newdim = onewnc.createDimension(dn+'bnds',4)
+            newdim = onewnc.createDimension(dn+'gridin',Ngridsinmax)
+        if not gen.searchInlist(onewnc.dimensions, 'coord'):
+            newdim = onewnc.createDimension('coord',2)
+
+        newvar = onewnc.createVariable(dn+'Ngrid','i',(dn+'slice'))
+
+        newvar[:] = Ngridsin[:]
+        basicvardef(newvar, dn+'Ngrid', "number of grids cells from " +     \
+          ".get. laying within .ref.", '-')
+        newvar.setncattr('coordinates',dn)
+
+        innewvar = onewnc.createVariable(dn+'gridin', 'i', ('coord',         \
+          dn+'gridin',dn+'slice'))
+        basicvardef(innewvar, dn+'gridin', "coordinates of the grids " +     \
+          "cells from .get. laying within .ref.",'-')
+        innewvar.setncattr('coordinates',dn+'slice')
+        
+        pnewvar = onewnc.createVariable(dn+'gridpercen','f',(dn+'gridin',    \
+          dn+'slice'))
+        basicvardef(pnewvar ,dn+'gridpercen', "percentages of the " +        \
+         "grids cells from .get. laying within .ref.", '1')
+        pnewvar.setncattr('coordinates',dn)
+        for j in range(1):
+           for i in range(Nslices):
+                innewvar[:,0:Ngridsin[j,i],i] = gridsin[:,0:Ngridsin[j,i],j,i]
+                pnewvar[0:Ngridsin[j,i],i]= percens[0:Ngridsin[j,i],j,i]
+        onewnc.sync()
+
+    for varn in varns:
+        print '    ' + varn + ' ...'
+        ovar = onc.variables[varn]
+        varv = ovar[:]
+        vdimns = list(ovar.dimensions)
+        vshape = list(ovar.shape)
+        vu = ovar.units
+
+        varbnds = []
+        for dn in vdimns:
+            if not gen.searchInlist(onewnc.dimensions,dn): add_dims(onc, onewnc, [dn])
+
+            odn = variables[dn]
+            dv = odn[:]
+
+            # Looking for boundaries
+            vdarattrs = odn.ncattrs()
+            if vdarattrs.has_key('bounds'):
+                print '    ' + infmsg
+                print '    ' +fname+ ": dimension variable '" +dn+ "' with bounds !!"
+                print '      bounds found:', varattrs['bounds']
+                print '      getting them to retrieve the slices'
+                bndvarns = varattrs['bounds']
+                varbnds.append(dn)
+                ovarbnds = onc.variables[varbnds]
+
+                # Slicing following boundaries, dimensions are 1D thus expand to 2D
+                #    filling as NW, NE, SE, SW
+                if gen.searchInlist(slcvarns, dn) and not vardimbndslice.has_key(dn):
+                    print infmsg
+                    print '    ' + fname + ": slicing dimension '" + dn + "' ... "
+                    varslcv = slicesinf[dn]
+                    Nslices = varslcv[0]
+                    ovardims = varslcv[1]
+                    slcvar = varslcv[2]
+
+                    # Shapping 2D slices
+                    reflat = np.zeros((Nslices), dtype=np.float)
+                    xslice2D = np.zeros((Nslices,4), dtype=np.float)
+                    xslice2D[:,0] = slcvar[1,:]
+                    xslice2D[:,1] = slcvar[1,:]
+                    xslice2D[:,2] = slcvar[0,:]
+                    xslice2D[:,3] = slcvar[0,:]
+                    yslice2D = np.zeros((Nslices,4), dtype=np.float)
+                    yslice2D[:,0] = -1.
+                    yslice2D[:,1] = 1.
+                    yslice2D[:,2] = 1.
+                    yslice2D[:,3] = -1.
+
+                    dd = len(onc.dimensions[dn])
+                    getlat = np.zeros((dd), dtype=np.float)
+                    xdimvarbnds2D = np.zeros((dd,4), dtype=np.float)
+                    xdimvarbnds2D[:,0] = ovarbnds[1,:]
+                    xdimvarbnds2D[:,1] = ovarbnds[1,:]
+                    xdimvarbnds2D[:,2] = ovarbnds[0,:]
+                    xdimvarbnds2D[:,3] = ovarbnds[0,:]
+                    ydimvarbnds2D = np.zeros((dd,4), dtype=np.float)
+                    ydimvarbnds2D[:,0] = -1.
+                    ydimvarbnds2D[:,1] = 1.
+                    ydimvarbnds2D[:,2] = 1.
+                    ydimvarbnds2D[:,3] = -1.
+
+                    reflont = slcvar.transpose()
+                    reflatt = reflat.transpose()
+                    refvarxbndst = xslice2D.transpose()
+                    refvarybndst = yslice2D.transpose()
+                    getlont = dv.transpose()
+                    getlatt = getlat.transpose()
+                    getvarxbndst = xdimvarbnds2D.transpose()
+                    getvarybndst = ydimvarbnds2D.transpose()
+
+                    Ngridsint, gridsint, percenst =                                  \
+                      fsci.module_scientific.spacepercen(xcavals=reflont,            \
+                      ycavals=reflatt, xbavals=refvarxbndst, ybavals=refvarybndst,   \
+                      xcbvals=getlont, ycbvals=getlatt, xbbvals=getvarxbndst,        \
+                      ybbvals=getvarybndst, strict=strict, dxa=Nslices, dya=1,       \
+                      navertexmax=4, dxb=dd, dyb=1, dxyb=dd, nbvertexmax=4)
+
+                    Ngridsin = Ngridsint.transpose()
+                    gridsin = gridsint.transpose()
+                    percens = percenst.transpose()
+
+                    Ngridsinmax = np.max(Ngridsin)
+
+                    # Let's avoid memory issues...
+                    lvalues = [Ngridsin, gridsin, percens]
+                    vardimbndslice[dn] = lvalues
+
+                    # Writting it to the file variables space-weight
+                    if not gen.searchInlist(onewnc.dimensions, dn+'bnds'):
+                        newdim = onewnc.createDimension(dn+'bnds',4)
+                        newdim = onewnc.createDimension(dn+'gridin',Ngridsinmax)
+                    if not gen.searchInlist(onewnc.dimensions, 'coord'):
+                        newdim = onewnc.createDimension('coord',2)
+
+                    newvar = onewnc.createVariable(dn + 'Ngrid','i',(dn))
+                    newvar[:] = Ngridsin[:]
+                    basicvardef(newvar, dn+'Ngrid', "number of grids cells from " +  \
+                      ".get. laying within .ref.", '-')
+                    newvar.setncattr('coordinates',dn)
+
+                    innewvar = onewnc.createVariable(dn+'gridin', 'i', ('coord',     \
+                      'gridin',dn))
+                    basicvardef(innewvar, dn+'gridin', "coordinates of the grids " + \
+                      "cells from .get. laying within .ref.",'-')
+                    innewvar.setncattr('coordinates',dn)
+
+                    pnewvar = onewnc.createVariable(dn+'gridpercen','f',('gridin',dn))
+                    basicvardef(pnewvar ,dn+'gridpercen', "percentages of the " +    \
+                      "grids cells from .get. laying within .ref.", '1')
+                    pnewvar.setncattr('coordinates',dn)
+
+                    for j in range(1):
+                        for i in range(dd):
+                            innewvar[:,0:Ngridsin[j,i],j,i] =                        \
+                              gridsin[:,0:Ngridsin[j,i],j,i]
+                            pnewvar[0:Ngridsin[j,i],j,i]= percens[0:Ngridsin[j,i],j,i]
+                    onewnc.sync()
+
+        # mask in var slice
+        maskvarslcs = {}
+        allNslices = []
+        slicedimn = []
+        for vslcn in slcvarns:
+            varslcv = slicesinf[vslcn]
+            Nslices = varslcv[0]
+            ovardims = varslcv[1]
+            slcvar = varslcv[2]
+
+            NOTcoinc = list(set(ovardims) - set(vdimns))
+            if len(NOTcoinc) > 0:
+                notcoincdim = {}
+                for dn in NOTcoinc: notcoincdim[dn] = 0
+            else:
+                notcoincdim = None
+            maskvarslice = []
+            for islc in range(Nslices):
+                varmask = gen.mat_dimreshape(list(ovardims), slcvar[islc,], vdimns,  \
+                  vshape, notcoincdim)
+                maskvarslice.append(varmask)
+
+            maskvarslcs[vslcn] = maskvarslice
+            allNslices.append(Nslices)
+            slicedimn.append('slice_'+vslcn)
+
+        newvarshape = allNslices + vshape
+        print '    Shape final sliced variable:', newvarshape
+
+        # variables
+  
+        # This can get quite huge !!
+        #newvard = slicedimn + vdimns
+        #newvar = onewnc.createVariable(varn+'sliced', 'f', tuple(newvard),           \
+        #  fill_value=gen.fillValueF)
+        #basicvardef(newvar, varn+'sliced', varn + 'sliced by '+ slcvarnS, vu)
+        #add_varattrs(onc, onewnc, [varn], [varn+'sliced'])
+
+        masknodims = []
+        vrank = np.arange(len(vdimns))
+        if slicestatsdim == 'any':
+            slcmin = np.ones(tuple(allNslices), dtype=ovar.dtype)*gen.fillValueF
+            slcmax = np.ones(tuple(allNslices), dtype=ovar.dtype)*gen.fillValueF
+            slcmean = np.ones(tuple(allNslices), dtype=ovar.dtype)*gen.fillValueF
+            slcstd= np.ones(tuple(allNslices), dtype=ovar.dtype)*gen.fillValueF
+            newvard = slicedimn
+        else:
+            nodims = gen.str_list(slicestatsdim, ',')
+            stvardims = vdimns + []
+            stvdn = []
+            stvid = []
+            idd = 0
+            vrank = list(vrank)
+            for nod in nodims:
+                if not gen.searchInlist(onc.dimensions, nod):
+                    print errormsg
+                    print '  ' +fname+ ": file '" + ncfile + "' does not have " +    \
+                      "dimension: '" + dn + "' !!"
+                    dimns = list(onc.dimensions)
+                    dimns.sort()
+                    print '    available ones:', dimns
+                    quit(-1)
+
+                if gen.searchInlist(vdimns, nod):
+                    stvdn.append(nod)
+                    stvid.append(vshape[idd])
+                    stvardims.remove(nod)
+                    vrank.pop(idd)
+                    masknodims.append(slice(0,vshape[idd]))
+                idd = idd + 1
+            vrank = np.array(vrank)
+
+            slcmin = np.ones(tuple(allNslices+stvid), dtype=ovar.dtype)*gen.fillValueF
+            slcmax = np.ones(tuple(allNslices+stvid), dtype=ovar.dtype)*gen.fillValueF
+            slcmean = np.ones(tuple(allNslices+stvid), dtype=ovar.dtype)*gen.fillValueF
+            slcstd= np.ones(tuple(allNslices+stvid), dtype=ovar.dtype)*gen.fillValueF
+            newvard = slicedimn + stvdn
+
+        newvarN = onewnc.createVariable(varn+'Nsliced', 'f', tuple(newvard),         \
+          fill_value=gen.fillValueF)
+        basicvardef(newvarN, varn+'Nsliced','number of values in slice of ' + varn + \
+          'sliced by '+ slcvarnS, vu)
+        add_varattrs(onc, onewnc, [varn], [varn+'Nsliced'])
+
+        newvarn = onewnc.createVariable(varn+'minsliced', 'f', tuple(newvard),       \
+          fill_value=gen.fillValueF)
+        basicvardef(newvarn, varn+'minsliced', 'minimum value of '+varn+'sliced by '+\
+          slcvarnS, vu)
+        add_varattrs(onc, onewnc, [varn], [varn+'minsliced'])
+
+        newvarx = onewnc.createVariable(varn+'maxsliced', 'f', tuple(newvard),       \
+          fill_value=gen.fillValueF)
+        basicvardef(newvarx, varn+'maxsliced', 'maximum value of '+varn+'sliced by '+\
+          slcvarnS, vu)
+        add_varattrs(onc, onewnc, [varn], [varn+'maxsliced'])
+
+        newvarm = onewnc.createVariable(varn+'meansliced', 'f', tuple(newvard),      \
+          fill_value=gen.fillValueF)
+        basicvardef(newvarm, varn+'meansliced', 'mean value of '+varn+'sliced by ' + \
+          slcvarnS, vu)
+        add_varattrs(onc, onewnc, [varn], [varn+'meansliced'])
+
+        newvars = onewnc.createVariable(varn+'stdsliced', 'f', tuple(newvard),       \
+          fill_value=gen.fillValueF)
+        basicvardef(newvars, varn+'stdsliced', 'standard deviation of '+varn+        \
+          'sliced by ' + slcvarnS, vu)
+        add_varattrs(onc, onewnc, [varn], [varn+'stdsliced'])
+
+        # Obviously we get memory problems... proceed slide by slide to fill the file
+        #newvarNmasked = np.ones(tuple(newvarshape), dtype=ovar.dtype)*gen.fillValueF
+        newvarNmasked = np.ones(tuple(vshape), dtype=ovar.dtype)*gen.fillValueF
+
+        # Masking!
+        allslices = gen.all_consecutive_combs(allNslices)
+        Nallslices = len(allslices)
+
+        for iallslc in range(Nallslices):
+            islcN = allslices[iallslc]
+            slcnewvar = []
+            iv = 0
+            for vslcn in slcvarns:
+                masksv = maskvarslcs[vslcn]
+                slcnewvar.append(islcN[iv])
+                if iv == 0: newmask = masksv[islcN[iv]]
+                else: newmask = newmask+masksv[islcN[iv]]
+                iv = iv + 1
+
+            for iid in vshape:
+               slcnewvar.append(slice(0,iid))
+
+            mavals = ma.array(varv, mask=newmask)
+            #newvarNmasked[tuple(slcnewvar)] = mavals.filled(gen.fillValueF)
+            newvarNmasked[:] = mavals.filled(gen.fillValueF)
+
+            # This can get quite huge !!!
+            #newvar[tuple(slcnewvar)] = newvarNmasked[:]
+
+            manewvar = ma.masked_equal(newvarNmasked, gen.fillValueF)
+ 
+            if len(masknodims) > 0: 
+                islcst = list(islcN) + list(masknodims)
+                if iallslc == 0:
+                    print '  variable slice statisitcs not using:', nodims
+                    print '    size statistics variables:', allNslices+stvid
+                    print '    slice statistics variables:', islcst
+                    print '    running shapes of the variable:', vrank
+            else: islcst = islcN
+
+            newvarN[tuple(islcst)]=np.sum(~manewvar.mask, axis=tuple(vrank))
+
+            slcmin[tuple(islcst)] = np.min(manewvar, axis=tuple(vrank))
+            slcmax[tuple(islcst)] = np.max(manewvar, axis=tuple(vrank))
+            slcmean[tuple(islcst)] = np.mean(manewvar, axis=tuple(vrank))
+            slcstd[tuple(islcst)] = np.std(manewvar, axis=tuple(vrank))
+
+        newvarn[:] = slcmin
+        newvarx[:] = slcmax
+        newvarm[:] = slcmean
+        newvars[:] = slcstd
+
+        onewnc.sync()
+        # Adding dimension-variables
+        for dn in vdimns:
+            if not dimvars.has_key(dn):
+                print errormsg
+                print '  ' + fname + ": no dimension-variable provided for " +       \
+                  "dimension '" + dn + "' !!"
+                print '    provided ones _______'
+                gen.printing_dictionary(dimvars)
+                quit(-1)
+            if not onewnc.variables.has_key(dimvars[dn]): 
+                add_vars(onc, onewnc, [dimvars[dn]])
+                onewnc.sync()
+
+    # Add global attributes
+    add_global_PyNCplot(onewnc, 'nc_var_tools', fname, '1.0')
+    add_globattrs(onc,onewnc,'all')
+    onewnc.sync()
+    onc.close()
+    onewnc.close()
+    print fname + ": successfull written of file '" + fname + ".nc' !!"
+
+    return
+
+values='lon,286.,323.6,4.;lat,-63.,19.,4.;orog,250.,7000.,500.@time|time:lon|lon:lat|lat@time'
+compute_slices_stats_areaweighted(values, '/media/lluis/ExtDiskC_ext3/DATA/estudios/Andes/DATA/concatenated/historical/tasmin/tasmin_Amon_ACCESS1-0_historical_r1i1p1_185001-200512_Andes_19600101000000-19900101000000.nc', 'tasmin')
 
 #quit()