Changeset 2352 in lmdz_wrf

Timestamp:

Feb 19, 2019, 5:48:44 PM (6 years ago)

Author:

lfita

Message:

Adding:

• `usefile_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 using pre-calculated area-weights from an existing file

Location:

trunk/tools

Files:

• nc_var.py (modified) (4 diffs)
• nc_var_tools.py (modified) (2 diffs)

trunk/tools/nc_var.py

@@ -72 +72 @@
 ## e.g. # nc_var.py -o compute_slices_stats_areaweighted -S 'lat,-63.,19.,2.;orog,500.,7000.,500.;rangefaces,fixed,-2.5|-0.5,-0.5|0.5,0.5|2.5@time|time:lon|lon:lat|lat@time@lon|lon_bnds,lat|lat_bnds@lon|lon_bnds,lat|lat_bnds@lat,lon@time' -f /media/lluis/ExtDiskC_ext3/DATA/estudios/Andes/DATA/concatenated/historical/tasmin/tasmin_Amon_ACCESS1-0_historical_r1i1p1_185001-200512_Andes_19600101000000-19900101000000.nc -v tasmin
 ## e.g. # nc_var.py -o except_fillValue -S 'orog:range,500.,7000.:None' -f /home/lluis/estudios/Andes/cmip_data/fx/orog_fx_ACCESS1-0_historical_r0i0p0.nc -v 'all'
+## e.g. # nc_var.py -o usefile_compute_slices_stats_areaweighted -S 'compute_slices_stats_areaweighted.nc@XLAT_M-HGT_M-rangefaces@Time|WRFtime:west_east|XLONG_M:south_north|XLAT_M:land_cat|INTrange@Time,land_cat' -f 'geo_em.d01.nc' -v 'LANDUSEF'
 
 from optparse import OptionParser
@@ -202 +203 @@
 # TimeSplitmat: Function to transform a file with CFtimes to a matrix [Nyear,Nmonth,Nday,Nhour,Nminute,Nsecond] 
 # timemean: Function to retrieve a time mean series from a multidimensional variable of a file
+# usefile_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 using pre-calculated area-weights from an existing file
 # valmod: Function to modify the value of a variable
 # valmod_dim: Function to modify the value of a variable at a given dimension and value
@@ -251 +256 @@
   'splitfile_dim', 'statcompare_files',                                              \
   'subbasin', 'submns', 'subyrs', 'temporal_stats', 'TimeInf', 'time_reset',         \
-  'TimeSplitmat', 'timemean', 'valmod', 'valmod_dim','varaddattrk', 'varaddattr',    \
+  'TimeSplitmat', 'timemean', 'usefile_compute_slices_stats_areaweighted',           \
+  'valmod', 'valmod_dim','varaddattrk', 'varaddattr',                                \
   'varaddref',                                                                       \
   'var_creation', 'varout', 'varoutold', 'varrmattr', 'varrm', 'VarVal_FillValue',   \
@@ -554 +560 @@
 elif oper == 'timemean':
     ncvar.timemean(opts.values, opts.ncfile, opts.varname)
+elif oper == 'usefile_compute_slices_stats_areaweighted':
+    ncvar.usefile_compute_slices_stats_areaweighted(opts.values, opts.ncfile,        \
+      opts.varname)
 elif oper == 'valmod':
     ncvar.valmod(opts.values, opts.ncfile, opts.varname)

trunk/tools/nc_var_tools.py

@@ -187 +187 @@
 # TimeSplitmat: Function to transform a file with CFtimes to a matrix [Nyear,Nmonth,Nday,Nhour,Nminute,Nsecond]
 # uploading_timestep: Function to upload a given time-step of a variable inside a netCDF
+# usefile_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 using pre-calculated area-weights from an existing file
 # valmod: Function to modify the value of a variable
 # valmod_dim: Function to modify the value of a variable at a given dimension and value
@@ -29724 +29728 @@
 #except_fillValue('orog:range,500.,7000.:None:value,0', '/home/lluis/estudios/Andes/cmip_data/fx/orog_fx_ACCESS1-0_historical_r0i0p0.nc', 'all')
 
+
+def usefile_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 using 
+        pre-calculated area-weights from an existing file
+      values=[areaweightsfile]@[slicen]@[dimvars]@[slicestatsdim]
+        [areaweightsfile]: file with the pre-calculated area-weights
+        [slicen]: name of the slice
+        [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. 
+            'WRFtime': to compute CF-compilant times from WRF's variable 'Times'
+            'INTrange': variable as dimension-length series of consecutive integers 
+        [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 = 'usefile_compute_slices_stats_areaweighted'
+
+    # Statistics from values within final slice
+    statns = ['min', 'max', 'mean', 'mean2', 'stddev', 'median', 'count']
+    Lstatns = {'min': 'space minimum', 'max': 'space maximum', 'mean': 'space mean', \
+      'mean2': 'space quadratic mean', 'stddev': 'spcae standard deviation',         \
+      'median': 'median', 'count': 'amount of grids'}
+    Ustatns = {'min': 'vu', 'max': 'vu', 'mean': 'vu', 'mean2': 'vu', 'stddev': 'vu',\
+      'median': 'vu', 'count': 'counts'}
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print usefile_compute_slices_stats_areaweighted.__doc__
+        quit()
+
+    expectargs = '[areaweightsfile]@[slicen]@[dimvars]@[slicestatsdim]'
+    gen.check_arguments(fname, values, expectargs, '@')
+
+    areaweightsfile = values.split('@')[0]
+    slicen = values.split('@')[1]
+    dimvars0 = values.split('@')[2]
+    slicestatsdim = values.split('@')[3].split(',')
+
+    dimvars = gen.stringS_dictvar(dimvars0, Dc=':', DVs='|')
+
+    if not os.path.isfile(areaweightsfile):
+        print errormsg
+        print '  '+fname+ ": pre calculated area-weights file '" + areaweightsfile + \
+          "' does not exist !!"
+        quit(-1)
+
+    awonc = NetCDFFile(areaweightsfile, 'r')
+
+    if variable == 'all':
+        varns = awonc.variables.keys()
+    else:
+        varns = variable.split(',')
+
+    vsurnames = ['Ngrid', 'gridin', 'areagridpercen', 'area']
+
+    invarn = slicen + vsurnames[0]
+
+    if not awonc.variables.has_key(invarn):
+        print errormsg
+        print '  '+fname+ ": pre calculated area-weights file '" + areaweightsfile + \
+          "' does not have slice '" + invarn + "' !!"
+        Varns = list(awonc.variables)
+        Varns.sort()
+        vvns = []
+        for vn in Varns:
+            Lvn = len(vn)
+            if vn[Lvn-5:Lvn] == 'Ngrid': vvns.append(vn)
+        print '  available ones:', vn
+        quit(-1)
+
+    osN = awonc.variables[slicen+vsurnames[0]]
+    osin = awonc.variables[slicen+vsurnames[1]]
+    osgpa = awonc.variables[slicen+vsurnames[2]]
+    osga = awonc.variables[slicen+vsurnames[3]]
+
+    sN = osN[:]
+    sin = osin[:]
+    sgpa = osgpa[:]
+
+    # removing monotones. Fortran only accepts rank <= 7 arrays !!
+    sN, Ndims = gen.remove_monotones(sN, osN.dimensions)
+    sin, indims = gen.remove_monotones(sin, osin.dimensions)
+    sgpa, gpadims = gen.remove_monotones(sgpa, osga.dimensions)
+
+    ssh = sN.shape
+    dimsslice = list(Ndims)
+    shapeslice = list(sN.shape)
+
+    sNt = sN.transpose()
+    #sint = sin.transpose()
+    sint0 = sin.transpose()
+    ## We are going to use them in Fortran thus (dx,dy) and 1st:1 !!!!
+    sint = np.ones((sint0.shape), dtype=np.float)*gen.fillValueF
+    sint[...,0] = sint0[...,1]
+    sint[...,1] = sint0[...,0]
+    sint = sint + 1
+    spt = sgpa.transpose()
+
+    print '  ' + infmsg
+    print '    final slice:', Ndims, 'shape:', sN.shape
+    onc = NetCDFFile(ncfile, 'r')
+
+    onewnc = NetCDFFile(fname + '.nc', 'w')
+
+    # Adding slicing dimensions and variables to the output file
+    slcvn = slicen.split('-')
+    for slcvn in slcvn:
+        vn = 'slice_' + slcvn
+        if not onewnc.variables.has_key(vn): add_vars(awonc, onewnc, [vn])
+        onewnc.sync()
+    if awonc.variables.has_key('slice_fake'): add_vars(awonc, onewnc, ['slice_fake'])
+            
+    awonc.close()
+
+    for varn in varns:
+        if not onc.variables.has_key(varn):
+            print errormsg
+            print '  '+fname+ ": file '" + ncfile + "' does not have variable '" + varn + "' !!"
+            Varns = list(onc.variables)
+            Varns.sort()
+            print '  available ones:', Varns
+            quit(-1)
+        print '    ' + varn + ' ...'
+
+        ovar = onc.variables[varn]
+        vdimns = list(ovar.dimensions)
+        vshape = list(ovar.shape)
+        vu = ovar.units
+
+        # Looking for possible new dimensions and dim-variables
+        vardims = []
+        varshape = []
+        idd = 0
+        for dnn in vdimns:
+            vdimn = dimvars[dnn]
+            if vdimn == 'WRFtime':
+                vdimn = 'time'
+                if not onewnc.variables.has_key('time'): 
+                    odimvar = onc.variables['Times']
+                    timewrfv = odimvar[:]
+                    tvals, urefvals = compute_WRFtime(timewrfv,                      \
+                      refdate='19491201000000', tunitsval='minutes')
+                    tvals = np.array(tvals)
+                    dimtwrf = len(tvals)
+                    if not gen.searchInlist(onewnc.dimensions,vdimn):
+                        newdim = onewnc.createDimension(vdimn, None)
+                    if not gen.searchInlist(onewnc.variables,vdimn):
+                        newvar = onewnc.createVariable(vdimn, 'f', ('time'))
+                        newvar[:] = tvals
+                        basicvardef(newvar, 'time', 'Time', urefvals)
+                        newvar.setncattr('calendar','gregorian')
+                        onewnc.sync()
+                if gen.searchInlist(slicestatsdim, dnn):
+                    vardims.append(vdimn)
+                    varshape.append(dimtwrf)
+            elif vdimn == 'INTrange':
+                vdimn = dnn + '_intrange'
+                if not onewnc.variables.has_key(vdimn): 
+                    Lvdim = len(onc.dimensions[dnn])
+                    if not gen.searchInlist(onewnc.dimensions,dnn):
+                        newdim = onewnc.createDimension(dnn, Lvdim)
+                    if not gen.searchInlist(onewnc.variables,vdimn):
+                        print "    creation of '" + vdimn + "' "
+                        newvar = onewnc.createVariable(vdimn, 'i', (dnn))
+                        newvar[:] = np.arange(Lvdim)
+                        basicvardef(newvar, vdimn, vdimn, '1')
+                        onewnc.sync()
+                if gen.searchInlist(slicestatsdim, dnn):
+                    vardims.append(dnn)
+                    varshape.append(Lvdim)
+            else:
+                if gen.searchInlist(slicestatsdim, dnn):
+                    vardims.append(vdimn)
+                    varshape.append(len(onc.dimensions[dnn]))
+                    if not gen.searchInlist(onewnc.dimensions,dnn): 
+                        add_dims(onc, onewnc, [dnn])
+                    if not gen.searchInlist(onewnc.variables,vdimn): add_vars(onc,   \
+                      onewnc, [vdimn])
+            idd = idd + 1
+
+        newvardims = vardims + dimsslice
+        newvarshape = varshape + shapeslice
+
+        print 'Lluis newvardims:', newvardims, 'newvarshape:', newvarshape
+        print 'Lluis file dims:', onewnc.dimensions.keys()
+        # Generic newvar slice
+        newvarslice = []
+        for idd in range(len(newvarshape)):
+            newvarslice.append(slice(0,newvarshape[idd]))
+
+        onewvars = {}
+        for statn in statns:
+            Lstn = Lstatns[statn]
+            if Ustatns[statn] == 'vu': vvu = vu
+            else: 
+                if statn == 'mean2': vvu = vu + '2'
+                else: vvu = Ustatns[statn]
+            newvar= onewnc.createVariable(varn+'sliced'+statn,'f', tuple(newvardims),\
+              fill_value = gen.fillValueF)
+            basicvardef(newvar, varn+'sliced'+statn, Lstn + ' ' + varn +             \
+              ' within slice by '+ slicen, vvu)
+            onewvars[statn] = newvar
+            onewnc.sync()
+
+        vout = None
+        voutt = None
+        varvt = None
+
+        # Checking for size...
+        Nvdims = len(vshape)
+        newvarslc = {}
+        varslc = {}
+        if Nvdims == 4:
+            vSIZE = np.prod(ovar.shape[:])*8.
+            if vSIZE > gen.oneGB*3.:
+                print '  ' + infmsg
+                print '    variable size: ', vSIZE / gen.oneGB,' larger than 3GB !!'
+                print '    splitting time evolution to compute'
+                vHSIZE = np.prod(ovar.shape[1:3])*8.
+                print '    var. horizontal size:', vHSIZE, 'dimt:', ovar.shape[0]
+                tsteps = gen.oneGB*3. / vHSIZE 
+                tslices = range(0,ovar.shape[0],tsteps)
+                print '    time-slices:', tslices
+                for itt in tslcices[1,len(tslices)]:
+                    iit = tslices[itt-1]
+                    eet = tslices[itt]
+                    varslc[itt-1] = [slice(iit,eet), slice(0,vshape[1]),             \
+                      slice(0,vbshape[2]), slice(0,vbshape[3])]
+                    newvarslc[itt-1] = [slice(iit,eet)] + newvarslice[1:4]
+            else:
+                tslices = range(0,ovar.shape[0])
+                varslc[0] = [slice(0,vshape[0]), slice(0,vshape[1]),                 \
+                  slice(0,vshape[2]), slice(0,vshape[3])]
+                newvarslc[0] =  newvarslice
+        elif Nvdims == 3:
+            vSIZE = np.prod(ovar.shape[:])*8.
+            if vSIZE > gen.oneGB*3.:
+                print '  ' + infmsg
+                print '    variable size: ', vSIZE / gen.oneGB,' larger than 3GB !!'
+                print '    splitting time evolution to compute'
+                vHSIZE = np.prod(ovar.shape[1:3])*8.
+                print '    var. horizontal size:', vHSIZE, 'dimt:', ovar.shape[0]
+                tsteps = gen.oneGB*3. / vHSIZE 
+                tslices = range(0,ovar.shape[0],tsteps)
+                print '    time-slices:', tslices
+                for itt in tslcices[1,len(tslices)]:
+                    iit = tslices[itt-1]
+                    eet = tslices[itt]
+                    varslc[itt-1] = [slice(iit,eet), slice(0,vshape[1]),             \
+                      slice(0,vbshape[2])]
+                    newvarslc[itt-1] = [slice(iit,eet)] + newvarslice[1:3]
+            else:
+                tslices = range(0,ovar.shape[0])
+                varslc[0] = [slice(0,vshape[0]), slice(0,vshape[1]),                 \
+                  slice(0,vshape[2])]
+                newvarslc[0] =  newvarslice
+        else:
+            slicevv = []
+            for iid in range(Nvdims):
+                slicevv.append(slice(0,vshape[0]))
+            varslc[0] = slicevv
+            newvarslc[0] = newvarslice
+
+        if Nvdims == 4:
+            d0 = ovar.shape[0]
+            d1 = ovar.shape[1]
+            d2 = ovar.shape[2]
+            d3 = ovar.shape[3]
+        elif Nvdims == 3:
+            d0 = ovar.shape[0]
+            d1 = ovar.shape[1]
+            d2 = ovar.shape[2]
+        elif Nvdims == 2:
+            d0 = ovar.shape[0]
+            d1 = ovar.shape[1]
+        elif Nvdims == 1:
+            d0 = ovar.shape[0]
+            if slicespacedim.count(vdimns[0]) != 1:
+                print errormsg
+                print '  ' + fname + ": there is no way to be able to compute " +    \
+                  "sliced statistics for variable '" + varn + "' because it " +      \
+                  "does not have any horizontal spatial dimension:", slicespacedim,  \
+                  "!!"
+                quit(-1)
+            ind = slicespacedim.index(vdimns[0])
+        else:
+            print errormsg
+            print '  ' + fname + ": rank ", Nvdims, " of variable '" + varn +        \
+              "' not ready !!"
+            quit(-1)
+
+        print '    computing statistics ...'
+
+        for itt in varslc.keys():
+           nslcv = newvarslc[itt]
+           slcv = varslc[itt]
+           varv = ovar[tuple(slcv)]
+
+           if Nvdims == 4:
+               d0 = varv.shape[0]
+               d1 = varv.shape[1]
+               d2 = varv.shape[2]
+               d3 = varv.shape[3]
+           elif Nvdims == 3:
+               d0 = varv.shape[0]
+               d1 = varv.shape[1]
+               d2 = varv.shape[2]
+           elif Nvdims == 2:
+               d0 = varv.shape[0]
+               d1 = varv.shape[1]
+           elif Nvdims == 1:
+               d0 = varv.shape[0]
+               if slicespacedim.count(vdimns[0]) != 1:
+                   print errormsg
+                   print '  ' + fname + ": there is no way to be able to compute " + \
+                     "sliced statistics for variable '" + varn + "' because it " +   \
+                     "does not have any horizontal spatial dimension:",              \
+                     slicespacedim, "!!"
+                   quit(-1)
+               ind = slicespacedim.index(vdimns[0])
+           else:
+               print errormsg
+               print '  ' + fname + ": rank ", Nvdims, " of variable '" + varn +     \
+                 "' not ready !!"
+               quit(-1)
+
+           varvt = varv.transpose()
+           if Nvdims == 3 and len(sN.shape) == 4:
+               voutt=fsci.module_scientific.multi_spaceweightstats_in3drk3_slc3v4(   \
+                 varin=varvt, ngridsin=sNt, gridsin=sint, percentages=spt,           \
+                 di1=d2, di2=d1, di3=d0, ds1=ssh[3], ds2=ssh[2], ds3=ssh[1],         \
+                 ds4=ssh[0], maxngridsin=sin.shape[1])
+           elif Nvdims == 4 and len(sN.shape) == 3:
+               voutt=fsci.module_scientific.multi_spaceweightstats_in4drk3_4_slc3v3( \
+                 varin=varvt, ngridsin=sNt, gridsin=sint, percentages=spt,           \
+                 di1=d3, di2=d2, di3=d1, di4=d0, ds1=ssh[2], ds2=ssh[1], ds3=ssh[0], \
+                 maxngridsin=sin.shape[1])
+           elif Nvdims == 3 and len(sN.shape) == 3:
+               voutt=fsci.module_scientific.multi_spaceweightstats_in3drk3_slc3v3(   \
+                 varin=varvt, ngridsin=sNt, gridsin=sint, percentages=spt,           \
+                 di1=d2, di2=d1, di3=d0, ds1=ssh[2], ds2=ssh[1], ds3=ssh[0],         \
+                 maxngridsin=sin.shape[1])
+           elif Nvdims == 2 and len(sN.shape) == 3:
+               voutt=fsci.module_scientific.multi_spaceweightstats_in2drkno_slc3v3(  \
+                 varin=varvt, ngridsin=sNt, gridsin=sint, percentages=spt,           \
+                 di1=d1, di2=d0, ds1=ssh[2], ds2=ssh[1], ds3=ssh[0],                 \
+                 maxngridsin=sin.shape[1])
+           elif Nvdims == 1 and len(sN.shape) == 3:
+               voutt=fsci.module_scientific.multi_spaceweightstats_in1drkno_slc3v3(  \
+                 varin=varvt, idv=ind, ngridsin=sNt, gridsin=sint, percentages=spt,  \
+                 di1=d0, ds1=ssh[2], ds2=ssh[1], ds3=ssh[0], maxngridsin=sin.shape[1])
+           else:
+               print errormsg
+               print '  ' + fname + ": variable rank:", Nvdims, 'combined with',     \
+                 ' slice rank:', len(sN.shape), 'not ready!!'
+               print '    available ones _______'
+               print ' var_rank: 3   slice_rank: 4'
+               print ' var_rank: 4   slice_rank: 3'
+               print ' var_rank: 3   slice_rank: 3'
+               print ' var_rank: 2   slice_rank: 3'
+               print ' var_rank: 1   slice_rank: 3'
+               quit(-1)
+
+           vout = voutt.transpose()
+           ist = 0
+           for statn in statns:
+               newvar = onewvars[statn]
+               rightvals = np.where(vout[6] == 0, gen.fillValueF, vout[ist])
+               newvar[tuple(nslcv)] = rightvals
+               onewnc.sync()
+               ist = ist + 1
+           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='compute_slices_stats_areaweighted.nc:XLAT_M-HGT_M-rangefaces'
+#usefile_compute_slices_stats_areaweighted(values, '/media/lluis/ExtDiskC_ext4/bkup/llamp/estudios/dominios/SA5k/geo_em.d01.nc', 'LANDUSEF')
+
 #quit()