Changeset 2463 in lmdz_wrf for trunk/tools

Timestamp:

Apr 23, 2019, 8:56:35 PM (7 years ago)

Author:

lfita

Message:

Adding:

• `fix_CFdates': Fixing CF-time values with wrong setting

File:

• trunk/tools/generic_tools.py (modified) (18 diffs)

trunk/tools/generic_tools.py

@@ -118 +118 @@
 # files_folder_HMT: Function to retrieve a list of files from a folder [fold] and files with [head]*[middle]*[tail]
 # fill_Narray: Function to fill a n-dimensional array with an arrary of lesser rank
+# fix_CFdates: Fixing CF-time values with wrong setting
+# gen_impose_gregorian: Function to impose gregorian calendar to a series of times with a 
+#   non-standard calendar independently from datetime (imposes yr > 1900)
 # gen_leapdays: Function to provide how many leap days are between two years avoiding datetime
 # generate_CFtimes: Function to generate CFtimes for a given period, frequency and units
@@ -519 +522 @@
     return combos
 
+def sign(number):
+    """ Function to provide the sign of a number
+      num: number to provide the sign
+    >>> sign(7109)
+    1
+    >>> sign(-3.43512)
+    -1.0
+    """
+    fname = 'sign'
+
+    sign = np.abs(number)/number
+
+    return sign
+
+##
+####
+######
+####### Temporal functions #######
+#####
+###
+#
 def days_month(year,month):
     """ Function to give the number of days of a month of a given year
@@ -592 +616 @@
 
     return Nleapdays
+
+def juliandate(yr,mo,dd,hr,mi,ss):
+    """ Function to provide the julian date (in days since -4712 January 1st 00 UTC) 
+          for any greogiran date
+      refdate: 1 January 4713 BC (= -4712 January 1), Greenwich mean noon (= 12h UT).
+    FROM: https://aa.usno.navy.mil/faq/docs/JD_Formula.php
+      after 1990 edition of the U.S. Naval Observatory's Almanac for Computers
+    >>> juliandate(1877,8,11,7,30,0)
+    2406843.3125
+    >>> juliandate(877,8,11,7,30,0)
+    2041601.3125
+    >>> juliandate(1970,1,1,0,0,0)
+    2440588.0
+    >>> juliandate(1978,1,1,0,0,0)
+    2443510.0
+    >>> juliandate(1978,7,21,15,0,0)
+    2443712.625
+    >>> juliandate(2000,1,1,0,0,0)
+    2451545.0
+    """
+    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)
+
+    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
+
+def juliandate_gregorian(juliand):
+    """ Function to convert from julian date (days since -4712 January 1st 00 UTC) 
+         to the equivalent greogiran date
+      juliand: julian date to convert
+    FROM: https://aa.usno.navy.mil/faq/docs/JD_Formula.php
+      after 1990 edition of the U.S. Naval Observatory's Almanac for Computers
+    >>> juliandate_gregorian(2406843.3125)
+    [1877, 8.0, 11, 7, 30, 0.0]
+    >>> juliandate_gregorian(2041601.3125)
+    [877, 8.0, 11, 7, 30, 0.0]
+    >>> juliandate_gregorian(2440588.0)
+    [1970, 1.0, 1, 0, 0, 0.0]
+    >>> juliandate_gregorian(2443510.0)
+    [1978, 1.0, 1, 0, 0, 0.0]
+    >>> juliandate_gregorian(2443712.625)
+    [1978, 7.0, 22, 15, 0, 0.0]
+    """
+    fname = 'juliandate_gregorian'
+
+    L = int(juliand+68569)
+    N = 4*L/146097
+    L = L-(146097*N+3)/4
+    I = 4000*(L+1)/1461001
+    L = L-1461*I/4+31
+    J = 80*L/2447
+    K = L-2447*J/80
+    L = J/11
+    J = J+2.-12.*L
+    I = 100*(N-49)+I+L
+
+    yr = I
+    mo = J
+    dd = K
+
+    # hr, mi, ss
+    ut = juliand - int(juliand)
+    utsecs = ut*24.*3600.
+    hh = int(utsecs/3600.)
+    mi = int((utsecs - hh*3600.)/60.)
+    ss = utsecs - hh*3600. - mi*60.
+
+    return [yr, mo, dd, hh, mi, ss]
 
 def datetimeStr_datetime(StringDT):
@@ -1957 +2074 @@
       Tunits: CF units of time [Tunits] since [DateRef]
       cal: non standard calendar
-      >>> impose_gregorian(np.range(800,10000,100), 'days since 1949-12-01 00:00:00',
+      >>> impose_gregorian(np.arange(800,10000,100), 'days since 1949-12-01 00:00:00',
         'noleap')
       [ 800. 901.]
@@ -2058 +2175 @@
       Tunits: CF units of time [Tunits] since [DateRef]
       cal: non standard calendar
-      >>> gen_impose_gregorian(np.range(800,10000,100), 'days since 1850-01-01 00:00:00',
+      >>> gen_impose_gregorian(np.arange(800,1000,100), 'days since 1850-01-01 00:00:00',
         'noleap')
-      [ 800. 901.]
+      [ 801. 901.]
       >>> gen_impose_gregorian(-7277.460938, 'days since 1949-12-01 00:00:00', 'noleap')
       [-7272.460938]
@@ -2080 +2197 @@
         Srefdate = lTunits[2] + ' ' + lTunits[3]
     refdate = [int(Srefdate[0:4]), int(Srefdate[5:7]), int(Srefdate[8:10]),          \
-      int(Srefdate[11:13]), int(Srefdate[14,16]), int(Srefdate[17:19]) ]
+      int(Srefdate[11:13]), int(Srefdate[14:16]), int(Srefdate[17:19]) ]
 
     if type(timev) == type(range(2)):
@@ -2101 +2218 @@
 
     if Tu == 'weeks':
-        timevsecs = timev * 24. * 7 * 3600.
+        timevdays = timev * 7
     elif Tu == 'days':
-        timevsecs = timev * 24. * 3600.
+        timevdays = timev * 1.
     elif Tu == 'hours':
-        timevsecs = timev * 3600.
+        timevdays = timev / 24.
     elif Tu == 'minutes':
-        timevsecs = timev * 60.
+        timevdays = timev / (24.*60.)
     elif Tu == 'seconds':
-        timevsecs = timev * 1.
+        timevdays = timev / (24.*3600.)
     else:
         print errormsg
@@ -2118 +2235 @@
     newtimes = np.zeros((dimt), dtype=np.float64)
 
-    refsecs = juliandate(refdate[0], refdate[1], refdate[2], refdate[3], refdate[4], refdate[5])
-    refsecs = refsecs*3600.
+    refdays = juliandate(refdate[0], refdate[1], refdate[2], refdate[3], refdate[4], refdate[5])
 
     yrref = refdate[0]
     if cal == 'noleap' or cal == '365d':
         for it in range(dimt):
-            # Some issues in implementations/machines ...
-            # FROM: https://stackoverflow.com/questions/33566939/timedelta-error-with-numpy-longdouble-dtype
-            tv = int(timevsecs[it])
-            deltaT = dt.timedelta(seconds=tv)
-            datetime = refdate + deltaT
-            yr = datetime.year
+            tdays = refdays + timevdays[it]
+            gdatev = juliandate_gregorian(tdays)
+
+            yr = gdatev[0]
+            mon = gdatev[1]
             if yr > yrref:
-                Nleapd = leapdays(yrref, yr)
+                Nleapd = gen_leapdays(yrref, yr)
             else:
-                Nleapd = leapdays(yr, yrref)
-            if datetime.month <= 2 and days_month(yr, 2) > 28: Nleapd = Nleapd - 1
-            newtimes[it] = timevsecs[it] + Nleapd*24*3600.
+                Nleapd = gen_leapdays(yr, yrref)
+            if mon <= 2 and leapyr(yr): Nleapd = Nleapd - 1
+
+            newtimes[it] = timevdays[it] + Nleapd
     else:
         print errormsg
@@ -2144 +2260 @@
     # Re-transforming to the original units
     if Tu == 'weeks':
-        newtimes = newtimes/(24.*7*3600.)
+        newtimes = newtimes/(7.)
     elif Tu == 'days':
-        newtimes = newtimes/(24.*3600.)
+        newtimes = newtimes*1.
     elif Tu == 'hours':
-        newtimes = newtimes/(3600.)
+        newtimes = newtimes*60.
     elif Tu == 'minutes':
-        newtimes = newtimes/(60.)
+        newtimes = newtimes*24.*3600.
 
     return newtimes
@@ -15954 +16070 @@
 #        print '  ' , '(', inds[i,j,0], ',', inds[i,j,1], ')'
 
-def sign(number):
-    """ Function to provide the sign of a number
-      num: number to provide the sign
-    >>> sign(7109)
-    1
-    >>> sign(-3.43512)
-    -1.0
-    """
-    fname = 'sign'
-
-    sign = np.abs(number)/number
-
-    return sign
-
-def juliandate(yr,mo,dd,hr,mi,ss):
-    """ Function to provide the julian date (in days since -4712 January 1st 00 UTC) 
-          for any greogiran date
-      refdate: 1 January 4713 BC (= -4712 January 1), Greenwich mean noon (= 12h UT).
-    FROM: https://aa.usno.navy.mil/faq/docs/JD_Formula.php
-      after 1990 edition of the U.S. Naval Observatory's Almanac for Computers
-    >>> juliandate(1877,8,11,7,30,0)
-    2406843.3125
-    >>> juliandate(877,8,11,7,30,0)
-    2041601.3125
-    >>> juliandate(1970,1,1,0,0,0)
-    2440588.0
-    >>> juliandate(1978,1,1,0,0,0)
-    2443510.0
-    >>> juliandate(1978,7,21,15,0,0)
-    2443712.625
-    """
-    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)
-
-    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
-
-def juliandate_gregorian(juliand):
-    """ Function to convert from julian date (days since -4712 January 1st 00 UTC) 
-         to the equivalent greogiran date
-      juliand: julian date to convert
-    FROM: https://aa.usno.navy.mil/faq/docs/JD_Formula.php
-      after 1990 edition of the U.S. Naval Observatory's Almanac for Computers
-    >>> juliandate_gregorian(2406843.3125)
-    [1877, 8.0, 11, 7, 30, 0.0]
-    >>> juliandate_gregorian(2041601.3125)
-    [877, 8.0, 11, 7, 30, 0.0]
-    >>> juliandate_gregorian(2440588.0)
-    [1970, 1.0, 1, 0, 0, 0.0]
-    >>> juliandate_gregorian(2443510.0)
-    [1978, 1.0, 1, 0, 0, 0.0]
-    >>> juliandate_gregorian(2443712.625)
-    [1978, 7.0, 22, 15, 0, 0.0]
-    """
-    fname = 'juliandate_gregorian'
-
-    L = int(juliand+68569)
-    N = 4*L/146097
-    L = L-(146097*N+3)/4
-    I = 4000*(L+1)/1461001
-    L = L-1461*I/4+31
-    J = 80*L/2447
-    K = L-2447*J/80
-    L = J/11
-    J = J+2.-12.*L
-    I = 100*(N-49)+I+L
-
-    yr = I
-    mo = J
-    dd = K
-
-    # hr, mi, ss
-    ut = juliand - int(juliand)
-    utsecs = ut*24.*3600.
-    hh = int(utsecs/3600.)
-    mi = int((utsecs - hh*3600.)/60.)
-    ss = utsecs - hh*3600. - mi*60.
-
-    return [yr, mo, dd, hh, mi, ss]
+#print juliandate_gregorian(2453408.5)
+#print juliandate_gregorian(2453683.5)
+
+#quit()
 
 def fix_CFdates(timevals, origcftimeu, origcal,                                      \
@@ -16070 +16086 @@
     >>> fix_CFdates(times, 'days since 1850-01-01 00:00:00', 'noleap', 
       'days since 1949-12-01 00:00:00')
-    [-36477.5 -36448.  -36418.5 -36388.  -36357.5 -36327.  -36296.5 -36265.5
-     -36235.  -36204.5 -36174. ]
+    [-36478.5 -36449.  -36419.5 -36389.  -36358.5 -36328.  -36297.5 -36266.5
+      -36236.  -36205.5 -36175. ]
+    >>> times = [56649.5, 56680, 56710.5, 56741, 56771.5, 56802.5, 56833, 56863.5, 56894, 56924.5]
+    >>> fix_CFdates(times, 'days since 1850-01-01 00:00:00', 'noleap', 
+      'days since 1949-12-01 00:00:00')
+    [ 20193.5  20224.   20254.5  20285.   20315.5  20346.5  20377.   20407.5
+      20438.   20468.5]
     """
     import datetime as dt
@@ -16085 +16106 @@
     origSrefdate = origtv[2]
 
-    yrref=origSrefdate[0:4]
-    monref=origSrefdate[5:7]
-    dayref=origSrefdate[8:10]
+    yrref = int(origSrefdate[0:4])
+    monref = int(origSrefdate[5:7])
+    dayref = int(origSrefdate[8:10])
 
     if int(yrref) < 1970: origsing = -1.
@@ -16110 +16131 @@
         secref = 0
 
-    seconds = juliandate(yrref,monref,dayref,horref,minref,secref)
-    seconds = seconds*24.*3600.
+    origdays = juliandate(yrref,monref,dayref,horref,minref,secref)
 
     # new
@@ -16118 +16138 @@
     newSrefdate = newtv[2]
 
-    yrnew=newSrefdate[0:4]
-    monnew=newSrefdate[5:7]
-    daynew=newSrefdate[8:10]  
+    yrnew = int(newSrefdate[0:4])
+    monnew = int(newSrefdate[5:7])
+    daynew = int(newSrefdate[8:10])  
 
     if int(yrnew) < 1970: newsing = -1.
@@ -16143 +16163 @@
         secnew = 0
 
-    seconds = juliandate(yrnew,monnew,daynew,hornew,minnew,secnew)
+    newdays = juliandate(yrnew,monnew,daynew,hornew,minnew,secnew)
 
     # Using shell date `coreutil' (a relic of the past...)
@@ -16150 +16170 @@
     #newsecondsdiff = np.float64(seconds.communicate()[0].replace('\n',''))
 
-    # If the calendar is not gregorian, it needs to be fixed from origSrefdate!
-    if origcal == 'noleap' or origcal == '365d':
-        print '  ' + fname + ": found non-standard calendar '" + origcal + "' !!"
-        # Imposing to gregorian calendar
-        newtimevals = gen_impose_gregorian(timevals, origcftimeu, origcal)
-
-    # Seconds difference respect new date
-    secsdiff = origsecondsdiff - newsecondsdiff
-
     if origtunits == 'years':
-        secsdtr = 365 * 24 * 3600.
+        daysdtr = 365.
     elif origtunits == 'weeks':
-        secsdtr = 7 * 24 * 3600.
+        daysdtr = 7.
     elif origtunits == 'days':
-        secsdtr = 24 * 3600.
+        daysdtr = 1.
     elif origtunits == 'hours':
-        secsdtr = 3600.
+        daysdtr = 1./24.
     elif origtunits == 'minutes':
-        secsdtr = 60.
+        daysdtr = 1./(24.*60.)
     elif origtunits == 'seconds':
-        secsdtr = 1.
+        daysdtr = 1./(24.*3600.)
     else:
         print errormsg
@@ -16177 +16188 @@
         quit(-1)
 
+    # If the calendar is not gregorian, it needs to be fixed from origSrefdate!
+    if origcal == 'noleap' or origcal == '365d':
+        print '  ' + fname + ": found non-standard calendar '" + origcal + "' !!"
+        # Imposing to gregorian calendar
+        timevals = gen_impose_gregorian(timevals, origcftimeu, origcal)     
+
+    # Seconds difference respect new date
+    daysdiff = origdays - newdays
+
     if type(timevals) != type(np.ones((2), dtype=int)):
         newtimevals = np.array(timevals, dtype=np.float64)
-        newtimevals = newtimevals*secsdtr
+        newtimevals = newtimevals*daysdtr
     else:
-        newtimevals = timevals*secsdtr
-
-    newtimevals = newtimevals + origsecondsdiff
-    newtimevals = newtimevals - newsecondsdiff
+        newtimevals = timevals*daysdtr
+
+    # Re-scaling to the absolute reference
+    newtimevals = newtimevals + origdays
+
+    # Re-scaling to the CF-time reference
+    newtimevals = newtimevals - newdays
 
     # Re-scaling to the new time units
     if newtunits == 'years':
-        secsdtr = 365 * 24 * 3600.
+        daysdtr = 1. / (365.)
     elif newtunits == 'weeks':
-        secsdtr = 7 * 24 * 3600.
+        daysdtr = 1. / 7.
     elif newtunits == 'days':
-        secsdtr = 24 * 3600.
+        daysdtr = 1.
     elif newtunits == 'hours':
-        secsdtr = 3600.
+        daysdtr = 24.
     elif newtunits == 'minutes':
-        secsdtr = 60.
+        daysdtr = 24. * 60.
     elif newtunits == 'seconds':
-        secsdtr = 1.
+        daysdtr = 24. * 3600.
     else:
         print errormsg
@@ -16204 +16227 @@
         print '    available ones: ', availtu
         quit(-1)
-    newtimevals = newtimevals / secsdtr
+    newtimevals = newtimevals * daysdtr
 
     return newtimevals
 
 #times = [15.5, 45, 74.5, 105, 135.5, 166, 196.5, 227.5, 258, 288.5, 319]
-#times = [56649.5, 56680, 56710.5, 56741, 56771.5, 56802.5, 56833, 56863.5, 56894, 56924.5]
-#print fix_CFdates(times, 'days since 1850-01-01 00:00:00', 'noleap', \
-#  'days since 1949-12-01 00:00:00')
+times = [56649.5, 56680, 56710.5, 56741, 56771.5, 56802.5, 56833, 56863.5, 56894, 56924.5]
+print fix_CFdates(times, 'days since 1850-01-01 00:00:00', 'noleap', \
+  'days since 1949-12-01 00:00:00')
 
 #quit()