Changeset 2165 in lmdz_wrf for trunk

Timestamp:

Oct 4, 2018, 10:44:41 PM (7 years ago)

Author:

lfita

Message:

Adding:

• `temporal_stats': Function to compute temporal statistics. Rank of the variables are preserved along other non-temporal dimensions

Location:

trunk/tools

Files:

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

trunk/tools/nc_var.py

@@ -60 +60 @@
 ## e.g. # nc_var.py -o curve_section -f /home/lluis/PY/test.nc -S 'gridline,x,y,8.,8.,16.,16.,32' -v all
 ## e.g. # nc_var.py -o merge_files -S 'plev|plev,time|time:merged.nc' -f 'ncobs/AliceSprings/snd_94326_198201-198201.nc,ncobs/AliceSprings/snd_94326_198202-198202.nc' -v 'plev,time'
+## e.g. # nc_var.py -o temporal_stats -S 'Time:WRFtime:day@1@min,LTday@-3@1@min:bottom_top@ZNU,south_north@XLAT,west_east@XLONG' -f wrfout_d01_1995-01-01_00:00:00 -v T2,Q2
 
 from optparse import OptionParser
@@ -171 +172 @@
 # submns: Function to retrieve a series of months from a file
 # subyrs: Function to retrieve a series of years from a file
+# temporal_stats: Function to compute temporal statistics. Rank of the variables are 
+#   preserved along other non-temporal dimensions
 # TimeInf: Function to print all the information from the variable time
 # time_reset: Function to re-set the time axis of a file
@@ -220 +223 @@
   'setvar_nc', 'sorttimesmat', 'spacemean', 'SpatialWeightedMean',                   \
   'splitfile_dim', 'statcompare_files',                                              \
-  'subbasin', 'submns', 'subyrs', 'TimeInf', 'time_reset',                           \
+  'subbasin', 'submns', 'subyrs', 'temporal_stats', 'TimeInf', 'time_reset',         \
   'TimeSplitmat', 'timemean', 'valmod', 'valmod_dim','varaddattrk', 'varaddattr',    \
   'varaddref',                                                                       \
@@ -497 +500 @@
 elif oper == 'subyrs':
     ncvar.subyrs(opts.values, opts.ncfile, opts.varname)
+elif oper == 'temporal_stats':
+    ncvar.temporal_stats(opts.values, opts.ncfile, opts.varname)
 elif oper == 'TimeInf':
     ncvar.TimeInf(opts.ncfile, opts.varname)

trunk/tools/nc_var_tools.py

@@ -158 +158 @@
 # subyrs: Function to retrieve a series of years from a file
 # testvarinfile: Check existence of a variable inside a netCDF file
+# temporal_stats: Function to compute temporal statistics. Rank of the variables are 
+#   preserved along other non-temporal dimensions
 # TimeInf: Function to print all the information from the variable time
 # time_information: Function to provide information about variable time
@@ -17663 +17665 @@
     dimvarns = {}
     for vardn in dimvars:
+
         dimn = vardn.split('@')[0]
         dimvn = vardn.split('@')[1]
@@ -25269 +25272 @@
     return [minval, maxval, meanval, stdval], [isw, jsw, ine, jne]
 
+def temporal_stats(values, ncfile, variable):
+    """ Function to compute temporal statistics. Rank of the variables are preserved 
+          along other non-temporal dimensions
+      values=[timedimn]:[timevarn]:[timestats]:[dimndimvarn]
+        [timedimn]: name of the dimension time
+        [timevarn]: name of the variable with CFtimes ('WRFtime' for times in WRF 
+          model format)
+        [timestats]:  ',' separated list of [period]@[amount]@[stat] periods and 
+            statistics
+          [period]: period to compute the statistics
+            'year': statistics within full years [01/01 00:00:00 - 12/31 23:59:59]
+            'month': statistics within full months [01 00:00:00 - 
+              [31,30,28/29] 23:59:59]
+            'day': statistics within full days [00:00:00 - 23:59:59]
+            'dayLT'@[UTCshift]: statistics within full local-time ([UTshift] local 
+              shift in hours from UTC) days [00:00:00+UTCshift - 23:59:59+UTCshift]
+            'hour': statistics within full hours [00:00 - 59:59]
+            'minute': statistics within full minutes [00 - 59]
+          [amount]: quantity of periods to group
+          [statn]: statistics to compute
+            'min': minimum value within [amunt]*[period]
+            'max': maximum value within [amunt]*[period]
+            'mean': mean value within [amunt]*[period]
+            'mean2': quadratic mean value within [amunt]*[period]
+            'std': stanard deviation value within [amunt]*[period]
+            'pecentiles'@[Npercen]: [Npercen] percentiles within [amunt]*[period]      
+          [dimndimvarn]: ',' list of [dimn]@[dimvarn] dimension and its associated 
+            variable dimension
+      ncfile= file from which provide the information
+      variable= ',' list of values to compute statistics ('all' for all variables 
+        only that ones with [timedimn])
+    """
+    fname = 'temporal_stats'
+
+    Lperiodn = {'year': 'yearly', 'month': 'monthly', 'day': 'daily',                \
+      'LTday': 'local-time daily', 'hour': 'hourly', 'minute': 'minutely'}
+    Lstatn = {'min': 'minimum', 'max': 'maximum', 'mean': 'mean',                    \
+      'mean2': 'quadratic mean', 'std': 'standard deviation',                        \
+      'percentiles': 'percentiles'}
+    availcftunits = ['days', 'hours', 'minutes', 'seconds']
+
+    if values == 'h':
+        print fname + '_____________________________________________________________'
+        print tepmoral_stats.__doc__
+        quit()
+
+    expectargs = '[timedimn]:[timevarn]:[timestats]:[dimndimvarn]'
+    gen.check_arguments(fname, values, expectargs, ':')
+
+    timedimn = values.split(':')[0]
+    timevarn = values.split(':')[1]
+    timestatsS = gen.str_list(values.split(':')[2], ',')
+    dimndimvarn = gen.str_list(values.split(':')[3], ',')
+
+    # Getting dim/dimvar correspondences
+    dimvarns = {}
+    for dvn in dimndimvarn:
+        dn = dvn.split('@')[0]
+        dvn = dvn.split('@')[1]
+        dimvarns[dn] = dvn
+        
+    onc = NetCDFFile(ncfile, 'r')
+
+    # Getting time values
+    if not gen.searchInlist(onc.dimensions, timedimn):
+        print errormsg
+        print '  ' +fname+ ": file '" + ncfile + "' does not have dimension time '"+ \
+          timedimn + "' !!"
+        dimns = onc.dimensions
+        dimns.sort()
+        print '    avaialble ones:', dimns
+        quit(-1)
+
+    if timevarn != 'WRFtime' and not onc.variables.has_key(timevarn):
+        print errormsg
+        print '  ' +fname+ ": file '" + ncfile + "' does not have variable time '" + \
+          timevarn + "' !!"
+        varns = onc.variables.keys()
+        varns.sort()
+        print '    avaialble ones:', varns
+        quit(-1)
+
+    if timevarn != 'WRFtime':
+        otimev = onc.variables[timevarn]
+        timev = otimev[:]
+        timeu = otimev.units
+        if gen.searchInlist(otimev.ncattrs(), 'calendar'): timec = otimev.calendar
+        else:
+            print warnsmg
+            print '    ' + fname + ": variable time without calendar attribute !"
+            print "    assuming 'standard/greogorian'"
+            timec = 'standard'
+    else:
+        otimev = onc.variables['Times']
+        timewrfv = otimev[:]
+        timev, timeu = compute_WRFtime(timewrfv, refdate='19491201000000',           \
+          tunitsval='minutes')
+        timec = 'standard'
+
+    tunits = timeu.split(' ')[0]
+    trefdate = timeu.split(' ')[2] + '_' + timeu.split(' ')[3]
+
+    # Getting periods
+    timestats = {}
+    iTst = 0
+    periods = []
+    amounts = []
+    stats = []
+    for timest in timestatsS:
+        timestatsv = []
+        period = timest.split('@')[0]
+        timestatsv.append(period)
+        if period == 'LTday':
+            UTCshift = int(timest.split('@')[1])
+            periods.append(period + '(' + str(UTCshift) + ')')
+            if tunits == 'days': UTCshift = UTCshift / 24.
+            elif tunits == 'hours': UTCshift = UTCshift / 1.
+            elif tunits == 'minutes': UTCshift = UTCshift * 60.
+            elif tunits == 'seconds': UTCshift = UTCshift * 3600.
+            else:
+                print errormsg
+                print '  ' + fname + ": CF temporal units '" + tunits +              \
+                  "' not ready to transform 'UTCshift' !!"
+                print '    available ones:', availcftunits
+                quit(-1)
+            timestatsv.append(UTCshift)
+            iamount = 2
+        else:
+            periods.append(period)
+            iamount = 1
+
+        amount = int(timest.split('@')[iamount])
+        amounts.append(amount)
+
+        timestatsv.append(amount)
+        statn = timest.split('@')[iamount+1]
+        timestatsv.append(statn)
+        if statn == 'percentiles': 
+            Npercen = int(timest.split('@')[iamount+2])
+            timestatsv.append(Npercen)
+            stats.append(statn + '(' + str(Npercen) + ')')
+        else:
+            stats.append(statn)
+
+        timestats[iTst] = timestatsv
+
+        iTst = iTst + 1
+
+    NTsts = iTst
+    print infmsg
+    print '    ' + fname + ': ', NTsts, 'Temporal statistics to compute _______'
+    for iTst in range(NTsts):
+        print timestats[iTst]
+
+    # Getting tempral time-slices for each period
+    statslices = {}
+    for iTst in range(NTsts):
+        timestatv = timestats[iTst]
+        period = timestatv[0]
+        if period == 'LTday':
+            UTCshift = timestatv[1]
+            timevc = timev + UTCshift
+            period = 'day'
+            iamount = 2
+        else:
+            timevc = timev + 0
+            iamount = 1
+        amount = timestatv[iamount]
+
+        timeslice, Ntimeslice = gen.time_slices(timevc, timeu, timec, period, amount)
+        statslices[iTst] = [Ntimeslice, timeslice]
+ 
+    # Creation of output file
+    ofile = fname + '.nc'
+    onewnc = NetCDFFile(ofile, 'w')
+
+    # Creation of dimensions
+    newdim = onewnc.createDimension('Tstat', NTsts)
+    newdim = onewnc.createDimension('Lstring', 256)
+    newdim = onewnc.createDimension('bounds', 2)
+
+    # Creation of variables
+    newvar = onewnc.createVariable('period', 'c', ('Tstat', 'Lstring'))
+    newvals = writing_str_nc(newvar, periods, 256)
+    basicvardef(newvar, 'period', 'Name of the time-period', '-')
+
+    newvar = onewnc.createVariable('statn', 'c', ('Tstat', 'Lstring'))
+    newvals = writing_str_nc(newvar, stats, 256)
+    basicvardef(newvar, 'statn', 'Name of the time-statistics', '-')
+
+    newvar = onewnc.createVariable('amount', 'i', ('Tstat'))
+    newvar[:] = amounts
+    basicvardef(newvar, 'amount', 'amount of the time-period', '1')
+
+    # Computing statistics
+    if variable == 'all':
+        varnames = onc.variables.keys()
+    else:
+        varnames = variable.split(',')
+
+    # Removing non-temporal variables
+    vartdim = {}
+    for varn in varnames:
+        ovar = onc.variables[varn]
+        if not gen.searchInlist(ovar.dimensions, timedimn): 
+            varnames.remove(varn)
+        else:
+            vartdim[varn] = gen.index_vec(ovar.dimensions, timedimn)
+       
+    for iTst in range(NTsts):
+        statinf = timestats[iTst]
+        sliceinf = statslices[iTst]
+
+        Ntslc = sliceinf[0]
+        tslcs = sliceinf[1]
+
+        periodn = statinf[0]
+        if periodn == 'LTday':
+            UTCshift = statinf[1]
+            iamount = 2
+        else:
+            iamount = 1
+        amount = statinf[iamount]
+        statn = statinf[iamount+1]
+
+        if statn == 'percentiles': Npercen = statinf[iamount+2]
+
+        print '    ' + fname + ": period '" + periodn + "' amount:", amount,         \
+          'statistics:', statn, " ..."
+
+        # Adding period dimension and variables
+        newdim = onewnc.createDimension(periodn+'_time', Ntslc)
+        newvart = onewnc.createVariable(periodn+'_time', 'f8',                       \
+          (periodn+'_time'))
+        basicvardef(newvart, periodn+'_time', 'time for ' + periodn, timeu)
+        newvart.setncattr('calendar', timec)
+        newvart.setncattr('bounds', periodn+'_time_bounds')
+        newvarb = onewnc.createVariable(periodn+'_time_bounds', 'f8',                \
+          (periodn+'_time', 'bounds'))
+        basicvardef(newvarb, periodn+'_time_bounds', 'time boundaries for '+ periodn,\
+          timeu)
+        newvarb.setncattr('calendar', timec)
+        for it in range(Ntslc):
+            timeslcev = tslcs[it]
+            islc = timevc[timeslcev[0]]
+            eslc = timevc[timeslcev[1]]
+
+            newvart[it] = (eslc + islc) / 2.
+            newvarb[it,:] = [islc, eslc]
+        onewnc.sync()
+
+        for varn in varnames:
+            print '      variable:', varn
+            ovar = onc.variables[varn]
+            vardims = ovar.dimensions
+
+            # Shaping statistics variable
+            shapetstat = []
+            dimntstat = []
+            iid = 0
+            for dn in vardims:
+                if dn == timedimn:
+                    dimntstat.append(periodn+'_time')
+                    shapetstat.append(Ntslc)
+                    itdim = iid
+                else:
+                    dimntstat.append(dn)
+                    Ldim = len(onc.dimensions[dn])
+                    shapetstat.append(Ldim)
+                    # inserting dimension and its variable
+                    if not gen.searchInlist(onewnc.dimensions, dn):
+                        od = onc.dimensions[dn]
+                        Ldim = len(od)
+                        if od.isunlimited(): newdim = onewnc.createDimension(dn,None)
+                        else: newdim = onewnc.createDimension(dn,Ldim)
+
+                    dimvn = dimvarns[dn]
+                    if not onewnc.variables.has_key(dimvn): 
+                        add_vars(onc,onewnc,[dimvn])
+                            
+                iid = iid + 1
+
+            if amount == 1:
+                newvarn = varn+periodn+statn
+            else:
+                newvarn = varn+str(amount)+periodn+statn
+
+            newvar = onewnc.createVariable(newvarn, 'f', tuple(dimntstat),          \
+              fill_value=gen.fillValueR)
+            varunits = get_varunits(ovar)
+            if statn == 'std': varunits = '(' + varunits + ')**2)'
+            CFvarattr = gen.variables_values(varn)
+            Lname = Lstatn[statn] + ' ' + Lperiodn[periodn] + ' ' +                  \
+              CFvarattr[4].replace('|',' ')
+            add_varattrs(onc,onewnc,[varn],[newvarn])
+            basicvardef(newvar, newvarn, Lname, varunits)
+            newvar.setncattr('cell_method', 'time: ' + Lstatn[statn])
+            newvar.setncattr('bounds', periodn+'_time_bounds')
+            if periodn == 'LTday': set_attributek(newvar, 'UTCshift', UTCshift, 'I')
+
+            tstat = np.ones(tuple(shapetstat), dtype=np.float)*gen.fillValueR
+
+            # Looping into the period
+            for islc in range(Ntslc):
+                timeslcv = tslcs[islc]
+                
+                timeslc = []
+                shapetimeslc = []
+                tstatslc = []
+                for dn in vardims:
+                    if dn == timedimn:
+                        timeslc.append(slice(timeslcv[0],timeslcv[1],timeslcv[2]))
+                        shapetimeslc.append(timeslcv[1]-timeslcv[0]+1)
+                        tstatslc.append(islc)
+                    else:
+                        Ldim = len(onc.dimensions[dn])
+                        timeslc.append(slice(0,Ldim))
+                        shapetimeslc.append(Ldim)
+                        tstatslc.append(slice(0,Ldim))
+
+                tvals = ovar[tuple(timeslc)]
+                if statn == 'min':
+                    tstat[tuple(tstatslc)] = np.min(tvals, axis=itdim)
+                elif statn == 'max':
+                    tstat[tuple(tstatslc)] = np.max(tvals, axis=itdim)
+                elif statn == 'mean':
+                    tstat[tuple(tstatslc)] = np.mean(tvals, axis=itdim)
+                elif statn == 'mean2':
+                    tstat[tuple(tstatslc)] = np.mean(tvals**2, axis=itdim)
+                elif statn == 'std':
+                    tstat[tuple(tstatslc)] = np.std(tvals, axis=itdim)
+                else:
+                    print errormsg
+                    print '  ' + fname + ": statisitcs '" + statn + "' not ready !!"
+                    statnor = Lstatn.keys()
+                    statnor.sort()
+                    print '    available ones:', statnor
+                    quit(-1)
+
+            newvar[:] = tstat
+            onewnc.sync()
+
+    # Adding PyNCplot
+    add_global_PyNCplot(onewnc, 'nc_var_tools', fname, '1.0')
+    onewnc.sync()
+
+    add_globattrs(onc, onewnc, 'all')
+    onewnc.sync()
+
+    onewnc.close()
+    print fname + ": successfull creation of file '" + ofile + "' !!"
+
+    return
+
+#vals = 'Time:WRFtime:day@1@min:bottom_top@ZNU,south_north@XLAT,west_east@XLONG'
+#temporal_stats(vals, '/home/lluis/PY/wrfout_d01_1995-01-01_00:00:00', 'T2')
+
 #quit()