Changeset 2281 in lmdz_wrf for trunk

Timestamp:

Jan 23, 2019, 9:02:40 PM (6 years ago)

Author:

lfita

Message:

Adding:

• `ovar_reducedims': Function to reduce a variable removing some of its dimensions

File:

• trunk/tools/nc_var_tools.py (modified) (15 diffs)

trunk/tools/nc_var_tools.py

@@ -133 +133 @@
 # operdim: Function to operate along a series of dimensions
 # ovar_onc: Function to copy an object variable to a nother netcdf object
+# ovar_reducedims: Function to reduce a variable removing some of its dimensions
 # Partialmap_Entiremap: Function to transform from a partial global map (e.g.: only land points) to an entire one Coincidence of points is done throughout a first guess from fractions of the total domain of search
 # Partialmap_EntiremapFor: Function to transform from a partial global map (e.g.: only land points) to an entire one usinf Fortran code
@@ -27690 +27691 @@
 #area_weighted('yes:min,max,mean,stddev',fvals,'ct_values')
 
+def ovar_reducedims(ovar, dimns):
+    """ Function to reduce a variable removing some of its dimensions
+      ovar: Object variable to reduce
+      dimns: list of dimensions to remove from the variable
+    """
+    fname = 'ovar_reducedims'
+ 
+    dimv = ovar.dimensions
+    
+    slicevar = []
+    dimvar = []
+    shapevar = []
+    iid = 0
+    for idn in dimv:
+        if gen.searchInlist(dimns,idn): slicevar.append(0)
+        else:
+            ldim = ovar.shape[iid]
+            slicevar.append(slice(0,ldim))
+            dimvar.append(idn)
+            shapevar.append(ldim)
+        iid = iid + 1
+
+    newvar = ovar[tuple(slicevar)]
+
+    return newvar, dimvar, shapevar
+
 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]
+      values=[slicevarsinf]@[dimvars]@[sliceremovedim]@[slicebndsdim]@[slicestatsdim]
         [slicevarsinf] ';' separated list of [varn1],[minvar1],[maxvar1],
           [slcevar1]; ...;[varnN],[minvarN],[maxvarN],[slcevarN]
@@ -27709 +27736 @@
           be able to compute the area-weights or it will be assumed that dimension
           has no extension
+        [sliceremovedim]: ',' list of dimensions to remove from the slicing variables
+          ('None' for any)
+        [slicebndsdim]: ',' [dimn]|[dimbnd] list of dimensions and their associated 
+          bounds variables from the slicing variables ('None' for any)
         [slicestatsdim]: ',' list of name of dimensions to do not use for computing 
           statistics of variables at each slice ('any' for not leaving any 
@@ -27722 +27753 @@
         quit()
 
-    expectargs = '[slicevarsinf]@[dimvars]@[slicestatsdim]'
+    expectargs = '[slicevarsinf]@[dimvars]@[sliceremovedim]@[slicebndsdim]@' +       \
+      '[slicestatsdim]'
     gen.check_arguments(fname, values, expectargs, '@')
 
     varvalues0 = values.split('@')[0]
     dimvars0 = values.split('@')[1]
-    slicestatsdim = values.split('@')[2]
+    sliceremovedim = values.split('@')[2]
+    slicebndsdim0 = values.split('@')[2]
+    slicestatsdim = values.split('@')[4].split(',')
 
     varvalues = varvalues0.split(';')
+
+    if sliceremovedim != 'None':
+        sliceremovedim = gen.str_list(sliceremovedim, ',')
+    else:
+        sliceremovdim = None
+    if slicebndsdim0 != 'None':
+        slicebndsdim0 = gen.str_list(slicebndsdim0, ',')
+        slicebndsdim = {}
+        for slcn in slicebndsdim0:
+            dn = slcn.split('|')[0]
+            bndn = slcn.split('|')[1]
+            slicebndsdim[dn] = bndn
+    else:
+        sliceremovdim = None
 
     varvals = {}
@@ -27772 +27820 @@
     # Dictionary with the result of slicing a dimension-variable with bounds
     vardimbndslice = {}
+    # Dimension variables of 2D
+    dimvar2D = []
 
     # Preparing slices and output file
@@ -27806 +27856 @@
         slcvalsc = np.zeros((Nslices), dtype=np.float)
         slcvals = np.zeros((Nslices,2), dtype=np.float)
-        varv = ovar[:]
+
+        # Removing undesired dimensions from slicing variable
+        if sliceremovedim is not None:
+            varv, rmdims, rmshape = ovar_reducedims(ovar, sliceremovedim)
+        else:
+            varv = ovar[:]
+            rmdims = dimnv
+            rmshape = list(ovar.shape)
+
         for islc in range(Nslices):
             slcvalsc[islc] = slices[islc]
@@ -27814 +27872 @@
             #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
+            print islc, ':', slcvals[islc,:], '=', np.sum(~slcvar[islc,])
 
         slicesinf[varn] = [Nslices, dimnv, slcvar, slcvalsc, slcvals]
@@ -27842 +27901 @@
           'variable ' + varn, vu)
 
-        slcdims = ['slice_'+varn] + dimnv
-        slcshape = [Nslices] + list(ovar.shape)
+        slcdims = ['slice_'+varn] + rmdims
+        slcshape = [Nslices] + rmshape
         slcv = np.zeros(tuple(slcshape), dtype=int)
         for islc in range(Nslices):
@@ -27856 +27915 @@
         # Looking for boundaries
         varbnds = []
+#        dn = dimvars[varn]
         dn = varn
         vdarattrs = ovar.ncattrs()
@@ -27863 +27923 @@
             boundsv = ovar.getncattr('bounds')
             print '    ' + infmsg
-            print '    ' +fname+ ": slicing variable '" +dn+ "' with bounds !!"
+            print '    ' +fname+ ": slicing variable '"+varn+ "' with dimension '" + \
+              dn + "' bounds !!"
             print '      bounds found:', boundsv
             print '      getting them to retrieve the slices'
@@ -27917 +27978 @@
                     dyget = 1
 
+                elif len(dv.shape) == 2:
+                    dimvar2D.append(dn)
+                    continue
+                # Issue with WRF's XLONG(time, sout_north, west_east) ....
                 else:
                     print errormsg
                     print '  '+fname+": rank ", dv.shape, " of slicing varibale '" + \
                       dn + "' not ready !!"
-                    print '    function only works with 1D slicing variables'
+                    print '    function only works with 1D and 2D slicing variables'
                     quit(-1)
 
@@ -27951 +28016 @@
             print '    directly using grid point values'
 
+            #if len(varv.shape) > 1: 
+            #    print errormsg
+            #    print '  ' + fname + ": wrong shape for slicing vriable '" + varn + "' !!"
+            #    print '    function only ready to slice 1D varables and it has:',        \
+            #      varv.shape
+            #    quit(-1)
+
             Ngridsin = np.zeros((1,Nslices), dtype=int)           
             for islc in range(Nslices):
@@ -27960 +28032 @@
             for islc in range(Nslices):
                 unmasked = slcvar[islc,]
-                indices = gen.multi_index_mat(unmasked,False)
+                unma = np.where(unmasked, 0., 1.)
+                unmat = unma.transpose()
+                if (len(unma.shape) == 1):
+                  indices = gen.multi_index_mat(unmat,1.)
+                  Nt = len(indices)
+                elif (len(unma.shape) == 2):
+                    drng = unma.shape[0]*unma.shape[1]
+                    Nt,indices=fsci.module_scientific.multi_index_mat2drk(d12=drng,  \
+                      mat=unmat, value=1., d1=unma.shape[1], d2=unma.shape[0])
+                elif (len(unma.shape) == 3):
+                    drng = unma.shape[0]*unma.shape[1]*unma.shape[2]
+                    Nt,indices=fsci.module_scientific.multi_index_mat3drk(d123=drng, \
+                      mat=unmat, value=1., d1=unma.shape[2], d2=unma.shape[1],       \
+                      d3=unma.shape[0])
+                elif (len(unma.shape) == 4):
+                    drng = unma.shape[0]*unma.shape[1]*unma.shape[2]*unma.shape[3]
+                    Nt,indices=fsci.module_scientific.multi_index_mat4drk(d1234=drng,\
+                      mat=unmat, value=1., d1=unma.shape[3], d2=unma.shape[2],       \
+                      d3=unma.shape[1], d4=unma.shape[0])
+                else:
+                  rS = str(len(unma.shape))
+                  print errormsg
+                  print '  ' + fname + ": 'multi_index_mat" + rS + "drk' not ready !!"
+                  quit(-1)
+
                 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]
+                    gridsin[0,0:Ngridsin[0,islc],0,islc] = ainds[0,0:Nt]
+                    gridsin[1,0:Ngridsin[0,islc],0,islc] = ainds[1,0:Nt]
                     percens[0:Ngridsin[0,islc],0,islc] = 1.
 
@@ -28002 +28098 @@
                 pnewvar[0:Ngridsin[j,i],i]= percens[0:Ngridsin[j,i],j,i]
         onewnc.sync()
+
+#    # Slicing 2D dimension-variables
+#    if len(dimvar2D) != 0:
+#        ovar = onc.variables[varn]
+#        vu = ovar.units
+#        dimnv = list(ovar.dimensions)
+
+#        varbndns = []
+#        dimvals = {}
+#        dimslcs = {}
+#        for dn in dimnv:
+#            # 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.replace(' ','').split(',')
+
+#                for vn in bndvarns:
+#                    if not gen.searchInlist(varbnds, vn): varbnds.append(vn)
+                
+#        # Slicing following boundaries
+#        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]
+#            slcvalsc = varslcv[3]
+#            slcvals = varslcv[4]
+#            dimslcs[dn] = slcvals        
+
+#         if len(varbndns) != 0:
+#             for boundsv in varbnds:
+#                ovarbnds = onc.variables[boundsv]
+#                varbnds = ovarbnds[:]
+#                dimvals[
+
+
 
     for varn in varns:
@@ -28306 +28449 @@
     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')
+#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')
+values='XLONG,-74.,-36.4,4.;XLAT,-63.,19.,4.;HGT,250.,7000.,500.@Time|WRFtime:' +    \
+  'west_east|XLONG:south_north|XLAT@south_north|lat_bnds,west_east|lon_bnds@Time@Time'
+compute_slices_stats_areaweighted(values, '/home/lluis/PY/wrfout_d01_1995-01-01_00:00:00', 'T2')
 
 #quit()