Changeset 2462 in lmdz_wrf for trunk/tools/generic_tools.py

Timestamp:

Apr 23, 2019, 7:43:12 PM (6 years ago)

Author:

lfita

Message:

Adding:

• `gen_leapdays': Function to provide how many leap days are between two years avoiding datetime

File:

• trunk/tools/generic_tools.py (modified) (10 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
+# 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
 # get_configuration: Function to get the configuration from an ASCII external file
 # get_right_CFtimeunits: Function to get the right CFtime units from any given format
 # get_specdictionary_HMT: Function to get specific values from a dictionary by selcting that keys with H*M*T
 # getting_fixedline: Function to get the values from a line of text with fixed lenght of different values
-# generate_CFtimes: Function to generate CFtimes for a given period, frequency and units
 # grib_CFequiv: Function to provide the CF name of a GRIB variable code number
 # inf_operSlist: Function to provide information from a string as a list separated by 
@@ -142 +143 @@
 # juliandate: Function to provide the julian date (in days since -4712 January 1st 00 UTC) 
 #   for any greogiran date
+# juliandate_gregorian: Function to convert from julian date (days since -4712 
+#   January 1st 00 UTC) to the equivalent greogiran date
 # latex_fig_array: Function to add an array of figures to an existing tex file
 # latex_text: Function to transform a text to LaTeX following style rules
@@ -546 +549 @@
     for iyr in range(yeari,yeare+1):
         if days_month(iyr,2) == 29: Nleapdays=Nleapdays+1
+
+    return Nleapdays
+
+def leapyr(yr):
+    """ Function to determine if it is a leap year (independently of datetime)
+    FROM: https://en.wikipedia.org/wiki/Leap_year
+      yr: year to check
+    >>> leapyr(2019)
+    False
+    >>> leapyr(1976)
+    True
+    >>> leapyr(1800)
+    False
+    >>> leapyr(1850)
+    False
+    """
+    fname = 'leapyr'
+
+    isleap = False
+    if np.mod(yr,4) == 0: isleap = True
+    if np.mod(yr,100) == 0: isleap = False
+    if np.mod(yr,400) == 0: isleap = True
+
+    return isleap
+
+def gen_leapdays(yeari, yeare):
+    """ Function to provide how many leap days are between two years avoiding datetime
+      yeari: initial year from [yeari] Jan 01
+      yeare: end year until [yearf] Dec 12
+    >>> gen_leapdays(1970,1980)
+    3
+    >>> gen_leapdays(1976,1976)
+    1
+    >>> gen_leapdays(1850,1976)
+    31
+    """
+    fname = 'gen_leapdays'
+
+    Nleapdays = 0
+    for iyr in range(yeari,yeare+1):
+        if leapyr(iyr): Nleapdays=Nleapdays+1
 
     return Nleapdays
@@ -1976 +2020 @@
 
     yrref = refdate.year
+    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
+            if yr > yrref:
+                Nleapd = 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.
+    else:
+        print errormsg
+        print '  ' + fname + ": calendar '" + cal + "' not ready !!"
+        print '    available ones:', availcalendar
+        quit(-1)
+        
+    # Re-transforming to the original units
+    if Tu == 'weeks':
+        newtimes = newtimes/(24.*7*3600.)
+    elif Tu == 'days':
+        newtimes = newtimes/(24.*3600.)
+    elif Tu == 'hours':
+        newtimes = newtimes/(3600.)
+    elif Tu == 'minutes':
+        newtimes = newtimes/(60.)
+
+    return newtimes
+
+def gen_impose_gregorian(timev, Tunits, cal):
+    """ Function to impose gregorian calendar to a series of times with a 
+        non-standard calendar independently from datetime (imposes yr > 1900)
+      timev: list of values
+      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',
+        'noleap')
+      [ 800. 901.]
+      >>> gen_impose_gregorian(-7277.460938, 'days since 1949-12-01 00:00:00', 'noleap')
+      [-7272.460938]
+    """
+    fname = 'gen_impose_gregorian'
+
+    availTu = ['weeks', 'days', 'hours', 'minutes', 'seconds']
+    availcalendar = ['noleap', '365d']
+
+    lTunits = Tunits.split(' ')
+    NTunits = len(lTunits)
+    Tu = lTunits[0]
+
+    if NTunits < 4: 
+        print infmsg
+        print '  ' + fname + ": CF-time units without time, adding '00:00:00' !!"
+        Srefdate = lTunits[2] + ' 00:00:00'
+    else:
+        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]) ]
+
+    if type(timev) == type(range(2)):
+        dimt = len(timev)
+        timev = np.array(timev)
+    elif type(timev) == type(np.arange(2)):
+        dimt = timev.shape[0]
+    elif type(timev) == type(int(2)) or type(timev) == type(np.float(2.)) or         \
+      type(timev) == type(float(2.)) or type(timev) == type(np.float32(2.)) or       \
+      type(timev) == type(np.float64(2.)):
+        dimt = 1
+        timev = np.ones(dimt, dtype=type(timev))*timev
+    else:
+        print errormsg
+        print '  ' + fname + ": type of time values: '", type(timev), "' 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)
+
+    if Tu == 'weeks':
+        timevsecs = timev * 24. * 7 * 3600.
+    elif Tu == 'days':
+        timevsecs = timev * 24. * 3600.
+    elif Tu == 'hours':
+        timevsecs = timev * 3600.
+    elif Tu == 'minutes':
+        timevsecs = timev * 60.
+    elif Tu == 'seconds':
+        timevsecs = timev * 1.
+    else:
+        print errormsg
+        print '  ' + fname + ": CF time units '" + Tu + "' not ready !!"
+        print '    available ones:', availTu
+        quit(-1)
+
+    newtimes = np.zeros((dimt), dtype=np.float64)
+
+    refsecs = juliandate(refdate[0], refdate[1], refdate[2], refdate[3], refdate[4], refdate[5])
+    refsecs = refsecs*3600.
+
+    yrref = refdate[0]
     if cal == 'noleap' or cal == '365d':
         for it in range(dimt):
@@ -15808 +15954 @@
 #        print '  ' , '(', inds[i,j,0], ',', inds[i,j,1], ')'
 
-def leapyr(yr):
-    """ Function to determine if it is a leap year (independently of datetime)
-    FROM: https://en.wikipedia.org/wiki/Leap_year
-      yr: year to check
-    >>> leapyr(2019)
-    False
-    >>> leapyr(1976)
-    True
-    >>> leapyr(1800)
-    False
-    >>> leapyr(1850)
-    False
-    """
-    fname = 'leapyr'
-
-    isleap = False
-    if np.mod(yr,4) == 0: isleap = True
-    if np.mod(yr,100) == 0: isleap = False
-    if np.mod(yr,400) == 0: isleap = True
-
-    return isleap
-
 def sign(number):
     """ Function to provide the sign of a number
@@ -15891 +16015 @@
 
     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 fix_CFdates(timevals, origcftimeu, origcal,                                      \
@@ -15943 +16110 @@
         secref = 0
 
-    addsecs = horref*3600 + minref*60 + secref
-    addsS = str(addsecs)
-
-    # Using shell date `coreutil'
-    seconds = sub.Popen('/bin/date +%s -u -d"' + origSrefdate[0:10] + ' ' + addsS +  \
-      ' seconds"', shell=True, stdout=sub.PIPE)
-    origsecondsdiff = np.float64(seconds.communicate()[0].replace('\n',''))
+    seconds = juliandate(yrref,monref,dayref,horref,minref,secref)
+    seconds = seconds*24.*3600.
 
     # new
@@ -15956 +16118 @@
     newSrefdate = newtv[2]
 
-    yrref=newSrefdate[0:4]
-    monref=newSrefdate[5:7]
-    dayref=newSrefdate[8:10]  
-
-    if int(yrref) < 1970: newsing = -1.
+    yrnew=newSrefdate[0:4]
+    monnew=newSrefdate[5:7]
+    daynew=newSrefdate[8:10]  
+
+    if int(yrnew) < 1970: newsing = -1.
     else: newsing = 1.
 
 # Does reference date contain a time value [YYYY]-[MM]-[DD] [HH]:[MI]:[SS]
 ##
-    trefT = newcftimeu.find(':')
-    if not trefT == -1:
+    tnewT = newcftimeu.find(':')
+    if not tnewT == -1:
         if len(newtv) == 3:
-            horref=int(newSrefdate[11:13])
-            minref=int(newSrefdate[14:16])
-            secref=int(newSrefdate[17:19])
+            hornew=int(newSrefdate[11:13])
+            minnew=int(newSrefdate[14:16])
+            secnew=int(newSrefdate[17:19])
         else:
             newSreftime = newtv[3]
-            horref=int(newSreftime[0:2])
-            minref=int(newSreftime[3:5])
-            secref=int(newSreftime[6:8])
+            hornew=int(newSreftime[0:2])
+            minnew=int(newSreftime[3:5])
+            secnew=int(newSreftime[6:8])
     else:
-        horref = 0
-        minref = 0
-        secref = 0
-
-    addsecs = horref*3600 + minref*60 + secref
-    addsS = str(addsecs)
-
-    # Using shell date `coreutil'
-    seconds = sub.Popen('/bin/date +%s -u -d"' + newSrefdate[0:10] + ' ' + addsS +   \
-      ' seconds"', shell=True, stdout=sub.PIPE)
-    newsecondsdiff = np.float64(seconds.communicate()[0].replace('\n',''))
+        hornew = 0
+        minnew = 0
+        secnew = 0
+
+    seconds = juliandate(yrnew,monnew,daynew,hornew,minnew,secnew)
+
+    # Using shell date `coreutil' (a relic of the past...)
+    #seconds = sub.Popen('/bin/date +%s -u -d"' + newSrefdate[0:10] + ' ' + addsS +   \
+    #  ' seconds"', shell=True, stdout=sub.PIPE)
+    #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
@@ -16015 +16182 @@
     else:
         newtimevals = timevals*secsdtr
+
     newtimevals = newtimevals + origsecondsdiff
     newtimevals = newtimevals - newsecondsdiff
@@ -16038 +16206 @@
     newtimevals = newtimevals / secsdtr
 
-    if origcal == 'noleap3' or origcal == '365d':
-        print '  ' + fname + ": found non-standard calendar '" + origcal + "' !!"
-        # Imposing to gregorian calendar
-        newtimevals = impose_gregorian(newtimevals, newcftimeu, origcal)
-
     return newtimevals