Changeset 1823 in lmdz_wrf for trunk/tools/nc_var_tools.py

Timestamp:

Mar 20, 2018, 1:07:45 AM (7 years ago)

Author:

lfita

Message:

Adding:

• `compute_WRFtime_bnds': Function to copmute CFtime_bnd from WRFtime variable and a given period

File:

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

trunk/tools/nc_var_tools.py

@@ -48 +48 @@
 # 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_WRFtime: function to copmute CFtimes from WRFtime variable
+# compute_WRFtime_bnds: Function to copmute CFtime_bnd from WRFtime variable and a given period
 # curve_section: Function to provide a section of a file following a given curve
 # DataSetSection: Function to get a section (values along a dimension) of a given data-set
@@ -21191 +21192 @@
     return
 
+
+def compute_WRFtime_bnds(timewrfv, period, kindWRFt='begperiod',                    \
+  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:
+        [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
+     kindWRFt: type of time from WRF
+       '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
+    """
+    import datetime as dt
+    fname = 'compute_WRFtime_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 = timewrfv.shape[0]
+
+    # Times period
+    begTmat = gen.datetimeStr_conversion(timewrfv[0,:], 'WRFdatetime', 'matYmdHMS')
+    endTmat= gen.datetimeStr_conversion(timewrfv[dimt-1,:], 'WRFdatetime','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)
+    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])
+
+    WRFtime_bnds = np.zeros((Nbnds,2), dtype=np.float)
+    WRFtime = np.zeros((Nbnds), dtype=np.float)
+    # For the time_bnds
+    if kindWRFt == 'begperiod':
+        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]
+            WRFtime_bnds[itb,0] = gen.realdatetime1_CFcompilant(newTBmati, refdate,  \
+              tunitsval)
+            WRFtime_bnds[itb,1] = gen.realdatetime1_CFcompilant(newTBmate, refdate,  \
+              tunitsval)
+            WRFtime[itb] = gen.realdatetime1_CFcompilant(newTmat, refdate, tunitsval)
+            
+    tunits = tunitsval + ' since ' + refdateS
+
+    return WRFtime_bnds, WRFtime, tunits
+
+#onc = NetCDFFile('/home/lluis/PY/wrfout_d01_1995-01-01_00:00:00', 'r')
+#otime = onc.variables['Times']
+#wtimev = otime[:]
+#tbnds, tc, tu = compute_WRFtime_bnds(wtimev, 'i,1', kindWRFt='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]
+
+#onc.close()
 #quit()