Changeset 2476 in lmdz_wrf

Timestamp:

Apr 26, 2019, 9:21:12 PM (6 years ago)

Author:

lfita

Message:

Adding:

• `from360d_reg': Function to transform from 360d = 12 * 30d calendar to a regular one

Location:

trunk/tools

Files:

• generic_tools.py (modified) (8 diffs)
• module_generic.f90 (modified) (1 diff)

trunk/tools/generic_tools.py

@@ -16 +16 @@
 from cStringIO import StringIO
 import sys
+import module_ForGen as gen
 
 main = 'nc_var_tools.py'
@@ -123 +124 @@
 # fill_Narray: Function to fill a n-dimensional array with an arrary of lesser rank
 # fix_CFdates: Fixing CF-time values with wrong setting
+# from360d_reg: Function to transform from 360d = 12 * 30d calendar to a regular one
 # fromanydate_CFYmd: Function to transform from any string format of date to Y-m-d used in 
 #   CF-conventions
@@ -638 +640 @@
     2443510.0
     >>> juliandate(1978,7,21,15,0,0)
-    2443712.625
+    2443711.625
     >>> juliandate(2000,1,1,0,0,0)
     2451545.0
@@ -644 +646 @@
     fname = 'juliandate'
 
-    sec1 = dd - 32075.
-    sec2 = 1461.*(yr+4800.+(mo-14)/12.)/4.
-    sec3 = 367.*(mo-2-(mo-14)/12.*12)/12.
-    sec4 = 3.*((yr+4900.+(mo-14.)/12.)/100.)/4.
- 
-    juliand = int(sec1 + sec2 + sec3 - sec4)
+    I= int(yr)
+    J= int(mo)
+    K= int(dd)
+
+    # With python it is not working!!
+    #sec1 = K - 32075.
+    #sec2 = 1461.*(I+4800.+(J-14.)/12.)/4.
+    #sec3 = 367.*(J-2.-(J-14.)/12.*12.)/12.
+    #sec4 = 3.*((I+4900.+(J-14.)/12.)/100.)/4.
+    #juliand = int(sec1 + sec2 + sec3 - sec4)
+    juliand = gen.module_generic.jd(I,J,K)
 
     juliand = juliand*1. + (hr + mi/60. + ss/3600.)/24.
-
-    # Alternative (without full agreement, just kept in case)
-    #K = yr + 0
-    #M = mo + 0
-    #I = dd + 0
-    #UT = hr*1. + mi/60. + ss/3600.
-
-    #sec1 = 367.*K
-    #sec2 = int((7.*(K+int((M+9)/12.)))/4.)
-    #sec3 = int((275.*M)/9.)
-    #sec4 = I*1.
-    #sec5 = 1721013.5
-    #sec6 = UT/24.
-    #sec7 = 0.5*sign(100.*K+M*1.-190002.5)
-    #sec8 = 0.5
-
-    #juliand = sec1 - sec2 + sec3 + sec4 + sec5 + sec6 - sec7 + sec8
-    #print '  juliand:', juliand
 
     return juliand
@@ -687 +676 @@
     >>> juliandate_gregorian(2443510.0)
     [1978, 1.0, 1, 0, 0, 0.0]
-    >>> juliandate_gregorian(2443712.625)
-    [1978, 7.0, 22, 15, 0, 0.0]
+    >>> juliandate_gregorian(2443711.625)
+    [1978, 7.0, 21, 15, 0, 0.0]
     """
     fname = 'juliandate_gregorian'
@@ -703 +692 @@
     I = 100*(N-49)+I+L
 
-    yr = I
-    mo = J
-    dd = K
+    yr = int(I)
+    mo = int(J)
+    dd = int(K)
 
     # hr, mi, ss
@@ -16473 +16462 @@
     """ Classs to provide information of CFtime units
       cftimeu: CF time units to provide information for [Timeu] since [RefDate]
+      * attributes:
+        refdate: reference date in Y-m-d H:M:S format
+        refdateYMHdms = reference date in YmdHMS format
+        refdatemat = reference date in matricial format
+        refdatejuliand = reference date as julian date (days since -4712 January 1st 
+          00 UTC)
+        refdateDT = reference date as datetime object
+        Tunits = origtunits
+        utoweeks = Factor to convert values from original temporal units to weeks
+        utodays = Factor to convert values from original temporal units to days
+        utohours = Factor to convert values from original temporal units to hours
+        utominutes = Factor to convert values from original temporal units to minutes
+        utoseconds = Factor to convert values from original temporal units to seconds
     """
 
     def __init__(self, cftimeu):
+        import datetime as dt
         fname = 'CFtimeU_inf'
+
+        availtu = ['weeks', 'days', 'hours', 'minutes', 'seconds']
+
         origtv = cftimeu.split(' ')
         origtunits = origtv[0]
@@ -16509 +16515 @@
         secrefS = str(secref).zfill(2)
 
+        # Transforming factors
+        if origtunits == 'weeks':
+            self.utoweeks = 1. 
+            self.utodays = 7.
+            self.utohours = 7. * 24.
+            self.utominutes = 7. * 24. * 60.
+            self.utoseconds = 7. * 24. * 3600.
+        elif origtunits == 'days':
+            self.utoweeks = 1. / 7.
+            self.utodays = 1.
+            self.utohours = 24.
+            self.utominutes = 24. * 60.
+            self.utoseconds = 24. * 3600.
+        elif origtunits == 'hours':
+            self.utoweeks = 1. / ( 7. * 24.)
+            self.utodays = 1. / 24.
+            self.utohours = 1.
+            self.utominutes = 60.
+            self.utoseconds = 3600.
+        elif origtunits == 'minutes':
+            self.utoweeks = 1. / ( 7. * 24. * 60.)
+            self.utodays = 1. / ( 24. * 60.)
+            self.utohours = 1. / 60.
+            self.utominutes = 1.
+            self.utoseconds = 60.
+        elif origtunits == 'seconds':
+            self.utoweeks = 1. / ( 7. * 24. * 3600.)
+            self.utodays = 1. / ( 24. * 3600.)
+            self.utohours = 1. / 3600.
+            self.utominutes = 1. / 60.
+            self.utoseconds = 1.
+        else:
+            print errormsg
+            print '  ' + fname + ": new time-units '" + origtunits + "' not ready !!"
+            print '    available ones: ', availtu
+            quit(-1)
+
         self.refdate = origSrefdate + ' ' + horrefS + ':' + minrefS + ':' + secrefS
         self.refdateYMHdms = origSrefdate[0:4] + origSrefdate[5:7] +                 \
           origSrefdate[8:10] + horrefS + ':' + minrefS + ':' + secrefS
         self.refdatemat = [yrref, monref, dayref, horref, minref, secref]
+        self.refdatejuliand = juliandate(yrref,monref,dayref,horref,minref,secref)
+        if yrref >= 1900: self.refdateDT = dt.datetime(yrref,monref,dayref,horref,   \
+          minref,secref)
+        else: self.refdateDT = None
+
         self.Tunits = origtunits
 
-def from360d_reg(tvals, tunits):
+def from360d_reg(tvals, tunits, kind='right'):
     """ Function to transform from 360d = 12 * 30d calendar to a regular one
       tvals: temporal values
       tunits: units of the values in CF format [Tunits] since [RefDate]
-    """
+      kind: kind of transformation to apply
+        'right': only transforming the right days of each month
+        'monpercen': transforming using the day as the percentage of month
+        'spread': spreading the dates to a regular year leaving the last 5 days empty
+    >>> from360d_reg(np.arange(15.5, 120+15.5,30), 'days since 1970-01-01 00:00:00', kind='right')
+    [  15.5   46.5   74.5  105.5]
+    >>> from360d_reg(np.arange(15.5, 120+15.5,30), 'days since 1970-01-01 00:00:00', kind='monpercen')
+    [  15.5   45.5   74.5  105.5]
+    """
+    import numpy.ma as ma
+    import datetime as dt
+
     fname = 'from360d_reg'
 
+    availkinds = ['right', 'monpercen']
+
     inftunits = CFtimeU_inf(tunits)
 
     refdate = inftunits.refdatemat
+
+    if type(tvals) == type(range(2)):
+        dimt = len(tvals)
+        tvals = np.array(tvals)
+        dT = tvals[1] - tvals[0]
+    elif type(tvals) == type(np.arange(2)):
+        dimt = tvals.shape[0]
+        dT = tvals[1] - tvals[0]
+    elif type(tvals) == type(int(2)) or type(tvals) == type(np.float(2.)) or         \
+      type(tvals) == type(float(2.)) or type(tvals) == type(np.float32(2.)) or       \
+      type(tvals) == type(np.float64(2.)):
+        dimt = 1
+        tvals = np.ones(dimt, dtype=type(tvals))*tvals
+        dT = 30.
+    else:
+        print errormsg
+        print '  ' + fname + ": type of time values: '", type(tvals), "' not ready !!"
+        print '    available ones:', type(range(2)), type(np.arange(2)),type(int(2)),\
+          type(np.float(2.)), type(float(2.)), type(np.float32(2.)),                 \
+          type(np.float64(2.))
+        quit(-1)
+
+    # Transforming to days
+    dtvals = tvals*inftunits.utodays
+    dyear = 12.*30
     
-
+    newtvals = tvals*0.
+
+    if kind == 'right':
+        # Here only the 'correct' dates will be transformed (right day of the month)
+        for it in range(dimt):
+
+            # Getting year, month and day
+            Nyears = int(tvals[it] / dyear)
+            yd = refdate[0] + Nyears
+            jd = tvals[it] - Nyears*dyear
+            md = int(jd/30) + 1
+            dd = int(jd - (md-1)*30)
+
+            # Getting hour, minute, second
+            tt = jd - int(jd)
+            tsecs = tt*24.*3600.
+            hh = int(tsecs/(3600.))
+            mi = int((tsecs - hh*3600.)/60.)
+            ss = tsecs - hh*3600. - mi*60.
+
+            # Number days of the month in the regular calendar
+            Nd_regmon = days_month(yd,md)
+            if dd > Nd_regmon: mewtvals[it] = fillvalueF
+            else:
+                juliand = juliandate(yd,md,dd,hh,mi,ss)
+                newtvals[it] = juliand - inftunits.refdatejuliand + 1
+                # Transforming to original units
+                newtvals[it] = newtvals[it]/inftunits.utodays
+    elif kind == 'monpercen':
+        # Here only the 'correct' dates will be transformed (day as percentage of the month)
+        for it in range(dimt):
+
+            # Getting year, month and day
+            Nyears = int(tvals[it] / dyear)
+            yd = int(refdate[0] + Nyears)
+            jd = tvals[it] - Nyears*dyear
+            md = int(jd/30) + 1
+            dd = int(jd - (md-1)*30)
+
+            # Getting hour, minute, second
+            tt = jd - int(jd)
+            tsecs = tt*24.*3600.
+            hh = int(tsecs/(3600.))
+            mi = int((tsecs - hh*3600.)/60.)
+            ss = tsecs - hh*3600. - mi*60.
+
+            # Number days of the month in the regular calendar
+            Nd_regmon = days_month(yd,md)
+            dd = int(dd*Nd_regmon/30.)
+            ddt = dt.datetime(yd,md,dd,hh,mi,int(ss)) - inftunits.refdateDT
+            juliand = juliandate(yd,md,dd,hh,mi,ss)
+            newtvals[it] = juliand - inftunits.refdatejuliand + 1
+            # Transforming to original units
+            newtvals[it] = newtvals[it]/inftunits.utodays
+    else:
+        print errormsg
+        print '  ' + fname + ": kind of transformtion '" + kind + "' not ready !!"
+        print '    available ones:', availkinds
+        quit(-1)
+
+    if dimt == 1:
+        newtt = newtvals[0] + 0.
+        newtvals = None
+        newtvals = newtt + 0.
+    else:
+        newtvals = ma.masked_equal(newtvals, fillValueF)
+    
     return newtvals
 

trunk/tools/module_generic.f90

@@ -1160 +1160 @@
   K= day
 
-  PRINT *,1461*(I+4800+(J-14)/12)/4, 1461*(I+4800+(J-14)/12)/4.
-
   JD= K-32075+1461*(I+4800+(J-14)/12)/4+367*(J-2-(J-14)/12*12)/12-3*((I+4900+(J-14)/12)/100)/4