Changeset 1857 in lmdz_wrf for trunk

Timestamp:

Mar 25, 2018, 7:51:05 PM (7 years ago)

Author:

lfita

Message:

Working version with 'time_bnds' generation (only tested for WRFtimes)

Location:

trunk/tools

Files:

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

trunk/tools/nc_var.py

@@ -19 +19 @@
 ## e.g. # nc_var.py -o reconstruct_matrix_from_vector -f cruncep_halfdeg_1958.nc -S 'nav_lon:nav_lat:-90.:-7.:-67.:15.:latlon:0.5:0.5:0.05' -v all 
 ## e.g. # nc_var.py -o nc_var.py -o CFmorzization -S 'X|west_east|XLONG,Y|south_north|XLAT,T|Time|WRFtime,Z|bottom_top|ZNU:GlobalAttr1995.inf:proj1995.inf' -f ~/PY/wrfout_d01_1995-01-01_00\:00\:00 -v QFX@instantaneous@None
+## e.g.# nc_var.py -o CFmorzization -S 'X|west_east|XLONG,Y|south_north|XLAT,T|Time|WRFtime,Z|bottom_top|ZNU:GlobalAttr1995.inf:proj1995.inf' -f ~/PY/wrfout_d01_1995-01-01_00\:00\:00 -v 'RAINC@time:!accumulation@WRFtime_bnds|h|3'
 
 ## e.g. ccrc468-17 # ./nc_var.py -v time -f 123/CCRC_NARCliM_Sydney_All_1990-1999_pr10max.nc -o out -S 1:-1

trunk/tools/nc_var_tools.py

@@ -48 +48 @@
 # compute_tevolboxtraj: Function to compute tevolboxtraj: temporal evolution at a given point along a box following a trajectory
 # compute_tevolboxtraj_radialsec: Function to compute tevolboxtraj_radialsec: temporal evolution at a given point along a number of radii at a given frequency following a trajectory
+# compute_time_bnds: Function to copmute CFtime_bnd from time variable and a given period
 # compute_WRFtime: function to copmute CFtimes from WRFtime variable
 # compute_WRFtime_bnds: Function to copmute CFtime_bnd from WRFtime variable and a given period
@@ -22496 +22497 @@
   refdate='19491201000000', tunitsval='minutes'):
     """ Function to copmute CFtime_bnd from WRFtime variable and a given period
-      refdate= [YYYYMMDDMIHHSS] format of reference date
-      tunitsval= CF time units
       timewrfv= matrix string values of WRF 'Times' variable
       period= [tunit],[quantity] period to compute the bounds as:
@@ -22506 +22505 @@
        'begperiod': WRF times represent the beginning of the period
        'centperiod': WRF times represent the center of the period
-       'endtperiod': WRF times represent the end of the period
+       'endperiod': WRF times represent the end of the period
+      refdate= [YYYYMMDDMIHHSS] format of reference date
+      tunitsval= CF time units
     """
     import datetime as dt
@@ -22652 +22653 @@
               tunitsval)
             WRFtime[itb] = gen.realdatetime1_CFcompilant(newTmat, refdate, tunitsval)
+    else:
+        print errormsg
+        print '  ' + fname + ": kind of WRF time variable: '" + kindWRFt +           \
+          "' not ready!!"
+        print '     available ones:', ['begperiod']
             
     tunits = tunitsval + ' since ' + refdateS
@@ -22671 +22677 @@
 
 #onc.close()
+
+def compute_time_bnds(timev, tunits, period, kindt='begperiod',                      \
+  refdate='19491201000000', tunitsval='minutes'):
+    """ Function to copmute CFtime_bnd from time variable and a given period
+      timev= values of time variable
+      tunits= CF-time units of the input time-values
+      period= [tunit],[quantity] period to compute the bounds as:
+        [tunit]: unit of time: 'c' century, 'y' year, 'm' month, 'w' week, 'd' day, 
+          'h' hour, 'i' minute, 's' second
+        [quantity] = amount (integer) of [tunit] to cover a period of time_bnds
+      kindt: type of time values
+       'begperiod': times represent the beginning of the period
+       'centperiod': times represent the center of the period
+       'endperiod': times represent the end of the period
+      refdate= [YYYYMMDDMIHHSS] format of reference date
+      tunitsval= CF time units
+    """
+    import datetime as dt
+    fname = 'compute_time_bnds'
+
+    tunit = period.split(',')[0]
+    quantity = int(period.split(',')[1])
+
+    yrref=refdate[0:4]
+    monref=refdate[4:6]
+    dayref=refdate[6:8]
+    horref=refdate[8:10]
+    minref=refdate[10:12]
+    secref=refdate[12:14]
+
+    refdateS = yrref + '-' + monref + '-' + dayref + ' ' + horref + ':' + minref +   \
+      ':' + secref
+
+    dimt = timev.shape[0]
+
+    # Times period
+    begTmat= gen.datetimeStr_conversion(str(timev[0]),'cfTime,'+tunits, 'matYmdHMS')
+    endTmat= gen.datetimeStr_conversion(str(timev[dimt-1]),'cfTime,'+tunits,'matYmdHMS')
+    CFbeg = gen.realdatetime1_CFcompilant(begTmat, refdate, tunitsval)
+    CFend = gen.realdatetime1_CFcompilant(endTmat, refdate, tunitsval)
+    begT = dt.datetime(begTmat[0], begTmat[1], begTmat[2], begTmat[3], begTmat[4],   \
+      begTmat[5])
+    endT = dt.datetime(endTmat[0], endTmat[1], endTmat[2], endTmat[3], endTmat[4],   \
+      endTmat[5])
+    dFCt = CFend - CFbeg
+
+    # Units for the variable 'time' to seconds
+    availtunitvals = ['centuries', 'years', 'days', 'hours', 'minutes', 'seconds']
+    if tunitsval == 'centuries':
+        tsecsuv = 100 * 365 * 24 * 3600.
+        dFCt = dFCt*100*365*24*3600
+    elif tunitsval == 'years':
+        tsecsuv = 365 * 24 * 3600.
+        dFCt = dFCt*365*24*3600
+    elif tunitsval == 'weeks':
+        tsecsuv = 7 * 24 * 3600.
+        dFCt = dFCt*7*24*3600
+    elif tunitsval == 'days':
+        tsecsuv = 24 * 3600.
+        dFCt = dFCt*24*3600
+    elif tunitsval == 'hours':
+        tsecsuv = 3600.
+        dFCt = dFCt*3600
+    elif tunitsval == 'minutes':
+        tsecsuv = 60.
+        dFCt = dFCt*60
+    elif tunitsval == 'seconds':
+        tsecsuv = 1.
+    else:
+        print errormsg
+        print '  ' + fname + ":' tunitsvals= '" + tunitsval + "' not ready!!"
+        print '    available ones: ', availtunitvals
+        quit(-1)
+
+    # Quantity of time_bnds
+    availtunits = ['c', 'y', 'm', 'd', 'h', 'i', 's']
+    if tunit == 'c':
+        dtbnds = quantity * 100 * 365 * 24 * 3600.
+        begTBmat = [int(begTmat[0]/100)*100, 1, 1, 0, 0, 0]
+        endTBmat = [int(endTmat[0]/100)*100+100, 1, 1, 0, 0, 0]
+        dTBmat = [quantity*100, 0, 0, 0, 0, 0]
+    elif tunit == 'y':
+        dtbnds = quantity * 365 * 24 * 3600.
+        begTBmat = [begTmat[0], 1, 1, 0, 0, 0]
+        endTBmat = [endTmat[0]+1, 1, 1, 0, 0, 0]
+        dTBmat = [quantity, 0, 0, 0, 0, 0]
+    elif tunit == 'm':
+        iyr = begTmat[0]
+        imo = begTmat[1]
+        eyr = endTmat[0]
+        emo = endTmat[1]
+        begTBmat = [begTmat[0], begTmat[1], 1, 0, 0, 0]
+        endTBmat = [endTmat[0], endTmat[1]+1, 1, 0, 0, 0]
+        dTBmat = [0, quantity, 0, 0, 0, 0]
+    elif tunit == 'w':
+        dtbnds = quantity * 7 * 24 * 3600.
+        begTBmat = [begTmat[0], begTmat[1], begTmat[2], 0, 0, 0]
+        endTBmat = [endTmat[0], endTmat[1], endTmat[2]+7, 0, 0, 0]
+        dTBmat = [0, 0, quantity*7, 0, 0, 0]
+    elif tunit == 'd':
+        dtbnds = quantity * 24 * 3600.
+        begTBmat = [begTmat[0], begTmat[1], begTmat[2], 0, 0, 0]
+        endTBmat = [endTmat[0], endTmat[1], endTmat[2]+1, 0, 0, 0]
+        dTBmat = [0, 0, quantity, 0, 0, 0]
+    elif tunit == 'h':
+        dtbnds = quantity * 3600.
+        begTBmat = [begTmat[0], begTmat[1], begTmat[2], begTmat[3], 0, 0]
+        endTBmat = [endTmat[0], endTmat[1], endTmat[2], endTmat[3]+1, 0, 0]
+        dTBmat = [0, 0, 0, quantity, 0, 0]
+    elif tunit == 'i':
+        dtbnds = quantity * 60.
+        begTBmat = [begTmat[0], begTmat[1], begTmat[2], begTmat[3], begTmat[4], 0]
+        endTBmat = [endTmat[0], endTmat[1], endTmat[2], endTmat[3], endTmat[4]+1, 0]
+        dTBmat = [0, 0, 0, 0, quantity, 0]
+    elif tunit == 's':
+        dtbnds = quantity * 1.
+        begTBmat = [begTmat[0], begTmat[1], begTmat[2], begTmat[3], begTmat[4],      \
+          begTmat[5]]
+        endTBmat = [endTmat[0], endTmat[1], endTmat[2], endTmat[3], endTmat[4],      \
+          endTmat[5]+1]
+        dTBmat = [0, 0, 0, 0, 0, quantity]
+    else:
+        print errormsg
+        print '  ' + fname + ":' tunit= '" + tunit + "' not ready!!"
+        print '    available ones: ', availtunits
+        quit(-1)
+
+    Nbnds = int((dFCt + dtbnds)/ dtbnds)
+    if Nbnds < 2:
+        print errormsg
+        print '  ' + fname + ": from: ", begTmat, 'to', endTmat, 'found', Nbnds,     \
+          'time_bnds of', quantity,'*',tunit, '!!'
+        print '    revise parameters !!'
+        quit(-1)
+    print '  ' + fname + ' from: ', begTmat, 'to', endTmat, '(', dFCt, ' seconds )', \
+      'found', Nbnds, 'time_bnds of', quantity,'*',tunit, '(', dtbnds, 's)'
+
+    begTB = dt.datetime(begTBmat[0], begTBmat[1], begTBmat[2], begTBmat[3],          \
+      begTBmat[4], begTBmat[5])
+    endTB = dt.datetime(endTBmat[0], endTBmat[1], endTBmat[2], endTBmat[3],          \
+      endTBmat[4], endTBmat[5])
+
+    time_bnds = np.zeros((Nbnds,2), dtype=np.float)
+    time = np.zeros((Nbnds), dtype=np.float)
+    # For the time_bnds
+    # In the way that the beginning/ending of the boundaries is computed there is not
+    #   distinction among 'beg/cent/end' time-kind values
+    for itb in range(Nbnds):
+        newTBi = begTB + dt.timedelta(seconds=dtbnds*itb)
+        newTBe = begTB + dt.timedelta(seconds=dtbnds*(itb+1))
+        newT = begTB + dt.timedelta(seconds=dtbnds*itb+dtbnds/2.)
+        newTBmati = [newTBi.year, newTBi.month, newTBi.day, newTBi.hour,             \
+          newTBi.minute, newTBi.second]
+        newTBmate = [newTBe.year, newTBe.month, newTBe.day, newTBe.hour,             \
+          newTBe.minute, newTBe.second]
+        newTmat = [newT.year, newT.month, newT.day, newT.hour, newT.minute,          \
+          newT.second]
+        time_bnds[itb,0] = gen.realdatetime1_CFcompilant(newTBmati, refdate,         \
+          tunitsval)
+        time_bnds[itb,1] = gen.realdatetime1_CFcompilant(newTBmate, refdate,         \
+          tunitsval)
+        time[itb] = gen.realdatetime1_CFcompilant(newTmat, refdate, tunitsval)
+            
+    tunits = tunitsval + ' since ' + refdateS
+
+    return time_bnds, time, tunits
+
+#timevals = np.array([23712570, 23712750, 23712930, 23713110, 23713290, 23713470,     \
+#    23713650,    \
+#    23713830, 23714010, 23714190, 23714370, 23714550, 23714730, 23714910,            \
+#    23715090, 23715270, 23715450, 23715630, 23715810, 23715990, 23716170,            \
+#    23716350, 23716530, 23716710, 23716890, 23717070, 23717250, 23717430,            \
+#    23717610, 23717790, 23717970, 23718150, 23718330, 23718510, 23718690,            \
+#    23718870, 23719050, 23719230, 23719410, 23719590, 23719770, 23719950,            \
+#    23720130, 23720310, 23720490, 23720670, 23720850, 23721030, 23721210,            \
+#    23721390, 23721570, 23721750, 23721930, 23722110, 23722290, 23722470])
+#tbnds, tc, tu = compute_time_bnds(timevals, 'minutes since 1949-12-01 00:00:00',     \
+#  'd,1', kindt='begperiod', refdate='19491201000000', tunitsval='minutes')
+
+#itb = gen.netCDFdatetime_realdatetime(tu, 'standard', tbnds[:,0])
+#etb = gen.netCDFdatetime_realdatetime(tu, 'standard', tbnds[:,1])
+#ctb = gen.netCDFdatetime_realdatetime(tu, 'standard', tc[:])
+
+#for it in range(len(tc)):
+#    print itb[it], ctb[it], etb[it]
+#quit()
 
 def CFmorzization(values, ncfile, variable):
@@ -23204 +23396 @@
                     tbvals, tvals, urefvals = compute_WRFtime_bnds(timewrfv, period, \
                       kindWRFt='begperiod', refdate=refT, tunitsval='minutes')
+                    # Re-setting time values to the center of the bounds
+                    ovtime[:] = tvals[:]
+                elif vtb[0:10] == 'Ctime_bnds':
+                    vtbvals = vtb.split('|')
+                    if len(vtbvals) != 4:
+                        print errormsg
+                        print '  ' + fname + ": to compute 'time_bnds' from a " +    \
+                          "given variable '[tvarn]', one needs to provide the name "+\
+                          "of the variable and the period of the bounds as: " +      \
+                          "bounds as: 'Ctime_bnds'|[tvarn]|[tunit]|[quantity]"
+                        print "    values passed: '" + vtb
+                        onc.close()
+                        quit(-1)
+                    if not gen.searchInlist(onc.variables.keys(), vtbvals[1]):
+                        print errormsg
+                        print '  ' + fname + ": file does not have variable '" +     \
+                          vtbvals[1] + "' to be used to compute 'time_bnds'"
+                        print '    available ones:', onc.variables.keys()
+
+                    period = vtbvals[2] + ',' + vtbvals[3]
+                    print infmsg
+                    print '      ' + fname + ": creation of variable 'time_bnds' " + \
+                      "from variable '" + tvtbvals[1] + "' using period '" + period +\
+                      "' !!"
+                    odimvar = onc.variables[tvtbvals[1]]
+                    timev = odimvar[:]
+                    ttu = odimvar.getncattr('units')
+
+                    refT = gen.datetimeStr_conversion(gattr['basetime'],             \
+                      'Y-m-dTH:M:SZ', 'YmdHMS')
+                    tbvals, tvals, urefvals = compute_time_bnds(timev, ttu, period,  \
+                      kindWRFt='begperiod', refdate=refT, tunitsval='minutes')
+                    # Re-setting time values to the center of the bounds
+                    ovtime[:] = tvals[:]
                 else:
                     if not onc.variables.has_key(vardimn):
@@ -23213 +23439 @@
                         quit(-1)
                     tbvals = onc.variables[vardimn]
-                # Re-setting time values to the center of the bounds
-                ovtime[:] = tvals[:]
-                onewnc.sync()
 
                 # creation of time_bnds variable
@@ -23222 +23445 @@
                 newvar = onewnc.createVariable('time_bnds', 'f8', ('time', 'nv'))
                 newvar[:] = tbvals[:]
+                onewnc.sync()
 
         # Adding projection variable