Merge pull request #103 from Cambridge-ICCS/python-fcDatetick

Python translation of fcDatetick (plus some extra testing for datenum and a fix)

oneflux_steps/ustar_cp/python-pisaac/cpdEvaluateUStarTh4Season20100901.py

@@ -10,7 +10,7 @@
 # local modules
 from fcDatevec import mydatevec
 from fcDoy import mydoy
-from fcDatenum import mydatenum
+from fcDatenum import datenum
 from fcBin import cpdBin
 from cpdFindChangePoint20100901 import cpdFindChangePoint20100901
 from myprctile import myprctile
@@ -71,7 +71,7 @@ def cpdEvaluateUStarTh4Season20100901(t, NEE, uStar, T, fNight, fPlot, cSiteYr):
     Y, M, D, h, m, s = mydatevec(t)
     # mydatevec returns float because it uses numpy.nan
     iYr = int(numpy.median(Y))
-    dn = mydatenum(iYr, 12, 31, 12, 0, 0)
+    dn = datenum(iYr, 12, 31, 12, 0, 0)
     EndDOY = mydoy(dn)
     nPerDay = int(round(1 / numpy.median(numpy.diff(t))))
     nSeasons = 4

oneflux_steps/ustar_cp/python-pisaac/fcDoy.py

@@ -1,6 +1,6 @@
 import numpy
 from fcDatevec import mydatevec
-from fcDatenum import mydatenum
+from fcDatenum import datenum
 
 def mydoy(t):
     """
@@ -21,6 +21,6 @@ def mydoy(t):
     yr0 = 2000
     Y, M, D, h, m, s = mydatevec(t)
     Y = Y + yr0
-    tt = mydatenum(int(Y), int(M), int(D), 0, 0, 0)
-    d = numpy.floor(tt - mydatenum(int(Y) - 1, 12, 31, 0, 0, 0))
+    tt = datenum(int(Y), int(M), int(D), 0, 0, 0)
+    d = numpy.floor(tt - datenum(int(Y) - 1, 12, 31, 0, 0, 0))
     return d

oneflux_steps/ustar_cp_python/__init__.py

@@ -1,8 +1,11 @@
 # Import single-function-per-file modules
 from .fcx2colvec import fcx2colvec
 from .fcx2rowvec import fcx2rowvec
+from .fcDatetick import fcDatetick
+from .fcDoy import fcDoy
 from .fcr2Calc import fcr2Calc
 from .cpdFmax2pCp3 import interpolate_FmaxCritical, cpdFmax2pCp3
+from .fcDatenum import datenum
 from .fcBin import fcBin
 
 from os.path import dirname, basename, isfile, join

oneflux_steps/ustar_cp_python/fcDatenum.py

@@ -1,8 +1,10 @@
 from datetime import datetime as dt
 from datetime import timedelta as td
-from math import ceil
+from math import ceil, floor
+from numpy import vectorize
+import numpy as np
 
-def mydatenum(Y, M, D):
+def datenum(Y : int | np.ndarray, M : int | np.ndarray, D : int| np.ndarray) -> int | np.ndarray:
     """
     Convert date to serial date number.
 
@@ -18,20 +20,30 @@ def mydatenum(Y, M, D):
     Returns:
     int: Serial date number corresponding to the input date.
 
+    Or all of these could be NumPy arrays of integers
+
     Notes:
     - Year 0 is treated as a leap year.
     - If the day (D) is 0, it is adjusted to the last day of the previous month.
     - If the month (M) is 0, it is adjusted to January.
     - Negative years are treated as if they are in the year 1 AD, with adjustments for leap years.
 
     Examples:
-    >>> mydatenum(2023, 10, 5)
+    >>> datenum(2023, 10, 5)
     739164
-    >>> mydatenum(0, 3, 1)
+    >>> datenum(0, 3, 1)
     61
-    >>> mydatenum(-1, 12, 31)
+    >>> datenum(-1, 12, 31)
     0
+    >>> datenum(0,24,31)
+    731
+    >>> datenum(1,24,31)
+    1096
     """
+    # mimic MATLAB's ability to handle scalar or vector inputs
+    if (hasattr(Y, "__len__") and len(Y) > 0 and hasattr(Y[0], "__len__")):
+        # Input is 2-Dimensional, so vectorise ourselves
+        return vectorize(datenum)(Y,M,D)
 
     adjustment = td()
     # A zero day means we need to subtract one day
@@ -41,12 +53,13 @@ def mydatenum(Y, M, D):
     # A zero month is interpreted as 1
     if M == 0:
       M = 1
-
-    # Treat year 0 (Y = 0) as a leap year, so we must
-    # add the leap day from Y = 0 if we are after the
-    # first leap day (Feb 29th of Year 0)
-    if (Y >= 1) | (Y == 0) & (M > 2):
-      adjustment = adjustment + td(1)
+    elif M > 12:
+      # If the month is greater than 12, we need to add the
+      # number of years to the year and adjust the month
+      Y = Y + (floor((M - 1) / 12))
+      M = ((M - 1) % 12) + 1
+      if M == 0:
+        M = 1
 
     if Y < 0:
         # If the year is negative, treat it as if we are in
@@ -58,8 +71,11 @@ def mydatenum(Y, M, D):
         # plus leap year corrections
         dn = dn + ceil(Y / 4)
     else:
-        d = dt(Y + 1, M, D, 0, 0, 0)
-        dn = d.toordinal()
+        # Adjust the year forwards by 1 to AD
+        # then by 3 to get the correct leap year calculate
+        # later compenating back by subtracting 3 years of non-leap years
+        d = dt(Y + 4, M, D, 0, 0, 0)
+        dn = d.toordinal() - (365*3)
 
     # turn adjustment (timedelta) in a number of days
     return dn + adjustment.days

oneflux_steps/ustar_cp_python/fcDatetick.py

@@ -0,0 +1,96 @@
+from oneflux_steps.ustar_cp_python.utils import floor, ceil, dot, arange, xlim, unique
+from oneflux_steps.ustar_cp_python.fcDatevec import fcDatevec
+from oneflux_steps.ustar_cp_python.fcDatenum import datenum
+import numpy as np
+import matplotlib.pyplot as plt
+import datetime
+
+def fcDatetick(t : float | np.ndarray, sFrequency : str, iDateStr : int, fLimits : float):
+    """
+    Generate date ticks for a plot based on the given time vector and frequency.
+
+    This function generates date ticks for a plot based on the provided time vector `t`,
+    the frequency `sFrequency`, the date string format `iDateStr`, and the plot limits `fLimits`.
+    It replicates some minimal behavior from the previous codebase for plotting purposes.
+
+    Parameters:
+    t (array-like): Time vector.
+    sFrequency (str): Frequency of the ticks. Possible values are "Dy", "14Dy", "Mo".
+    iDateStr (int): Date string format.
+    fLimits (float): Limits of the plot.
+
+    Returns:
+    None
+
+    Notes:
+    - This function is *not* used in the rest of the Python code.
+    - It is *not* thoroughly tested and only replicates some minimal behavior from the previous codebase for plotting.
+
+    Examples:
+    >>> t = list(range(0, 49))
+    >>> sFrequency = "Mo"
+    >>> iDateStr = 4
+    >>> fLimits = 1
+    >>> fcDatetick(t, sFrequency, iDateStr, fLimits)
+    """
+
+    y, m, d, h, mn, s = fcDatevec(t)
+    iYrs = unique(y)
+    iSerMos = dot((y - 1), 12) + m
+    iSerMo1 = min(iSerMos)
+    iSerMo2 = max(iSerMos)
+    nSerMos = iSerMo2 - iSerMo1 + 1
+    xDates = np.array([])
+
+    match (sFrequency):
+      case "Dy":
+        xDates = t[::48]
+
+      case "14Dy":
+        iYr1 = floor(iSerMo1 / 12) + 1
+        iMo1 = iSerMo1 % 12
+        if iMo1 == 0:
+            iMo1 = 12
+            iYr1 = iYr1 - 1
+        for iDy in arange(1, 15, 14).reshape(-1):
+            xDates = [datenum(iYr1, int(month), iDy) for month in arange(iMo1, (iMo1 + nSerMos))]
+
+      case "Mo":
+        iYr1 = floor(iSerMo1 / 12) + 1
+        iMo1 = iSerMo1 % 12
+        if iMo1 == 0:
+            iMo1 = 12
+            iYr1 = iYr1 - 1
+        # define xDates as the array given my mapping over
+        # every month in arange(iMo1, (iMo1 + nSerMos)) 
+        # and applying `datenum(iYr1, month, 1)` to it
+        xDates = [datenum(iYr1, int(month), 1) for month in arange(iMo1, (iMo1 + nSerMos))]
+        # oneflux_steps/ustar_cp_refactor_wip/fcDatetick.m:36
+        # # # 			datestr(xDates)
+        # # # 			datestr([min(t) max(t)])
+        # # # 			pause;
+
+    xDates = unique(xDates)
+
+    # Set current `x` access to have values given by xDates
+    plt.gca().set_xticks(xDates)
+    # set label to empty
+    plt.gca().set_xticklabels([])
+    if iDateStr > 0:
+        # compute a datestring for each xDates based on iDateStr
+        # TODO: this is different to the original MATLAB code which was
+        # cDates = datestr(xDates, iDateStr)
+
+        # convert from a date floating-point
+        # ordinal to a datetime object
+        cDates = [datetime.datetime.fromordinal(floor(t) + 1).strftime("%Y-%m-%d") for t in xDates]
+
+        plt.gca().set_xticklabels(cDates)
+
+    if fLimits == 1:
+        xlim([floor(min(xDates)), ceil(max(xDates))])
+        # Turn on the grid and box using matplotlib
+        plt.grid("on")
+        plt.box("on")
+
+    return None

oneflux_steps/ustar_cp_python/fcDatevec.py

@@ -24,6 +24,7 @@ def fcDatevec(t : numpy.ndarray) -> tuple:
     if (hasattr(t, "__len__") and len(t) > 0 and hasattr(t[0], "__len__")):
         # Input is 2-Dimensional, so vectorise ourselves
         return numpy.vectorize(fcDatevec)(t)
+
     t = numpy.asarray(t)
 
     # Quantise the input to the granularity of 0.0001 seconds in a day

oneflux_steps/ustar_cp_python/fcDoy.py

@@ -1,5 +1,5 @@
 from oneflux_steps.ustar_cp_python.fcDatevec import fcDatevec
-from oneflux_steps.ustar_cp_python.fcDatenum import mydatenum
+from oneflux_steps.ustar_cp_python.fcDatenum import datenum
 import numpy as np
 
 def fcDoy(t=None):
@@ -43,9 +43,9 @@ def convert(y, m, d):
         return np.vectorize(convert)(y, m, d)
       else:
         # Convert the date to a datetime object
-        tt = mydatenum(int(y), int(m), int(d)) # datenum(y, m, d)
+        tt = datenum(int(y), int(m), int(d)) # datenum(y, m, d)
         # Subtract the last day of the previous year
-        d = np.floor(tt - mydatenum(int(y - 1), 12, 31))
+        d = np.floor(tt - datenum(int(y - 1), 12, 31))
         return d
 
     d = np.vectorize(convert)(y, m, d)

oneflux_steps/ustar_cp_python/utils.py

@@ -232,17 +232,23 @@ def any(a):
     return np.any(a)
 
 
+def arange_column(start, stop, step=1, **kwargs):
+    """
+    >>> a=arange(1,10) # 1:10
+    >>> size(a)
+    matlabarray([[ 1, 10]])
+    """
+    expand_value = 1 if step > 0 else -1
+    return np.arange(start, stop + expand_value, step, **kwargs).reshape(1, -1),
+
 def arange(start, stop, step=1, **kwargs):
     """
     >>> a=arange(1,10) # 1:10
     >>> size(a)
     matlabarray([[ 1, 10]])
     """
     expand_value = 1 if step > 0 else -1
-    return matlabarray(
-        np.arange(start, stop + expand_value, step, **kwargs).reshape(1, -1),
-        **kwargs,
-    )
+    return np.arange(start, stop + expand_value, step, **kwargs)
 
 
 def concat(args, axis=1):
@@ -816,7 +822,7 @@ def unique(a):
     """
     Return the unique elements of an array.
     """
-    return matlabarray(np.unique(np.asarray(a)))
+    return np.unique(np.asarray(a))
 
 
 def interp1(x, v, xq, method):