Support arbitrary epoch for time inputs

Prior to this, time was assumed to be in the form of Matlab
datenums--but this is not what any typical Python code would
work with.  Now the default is matplotlib date numbers, and
any epoch can be specified.

doc/strategy.rst

@@ -38,32 +38,35 @@ Public interface
 The package is called `utide`, and presently has a very
 simple interface, with two function: `solve` and
 `reconstruct`.  These are simply English spellings of their
-slightly shortened Matlab counterparts.  Everything else
+Matlab counterparts.  Everything else
 should be considered to be private, regardless of whether it
 has a leading underscore.
 
-There is an overwhelming number of options; we might be able
-to find ways of making this interface friendlier.
+There is an overwhelming number of options, and many of the
+Matlab names are rather cryptic.  We have made some changes
+in the way options are specified, but the process is not
+complete.
 
-Options are being held internally in a `dict`.  We will
-switch to using a `Bunch` so as to provide both dictionary
-and attribute access syntax.
+Options are being held internally in a `Bunch` so as to
+provide both dictionary and attribute access syntax.
 
 Time
 ^^^^
-Presently, time inputs are assumed to be Matlab `datenum`
-arrays.  We need to make this more flexible, at the very
-least including the ability to handle time in days since a
-specified epoch. An array of Python datetime objects could
-be supported, but this is not top priority. At some point
-we will presumably handle the numpy datetime64 dtype, but we
-can wait until it has been reworked and is no longer in a
-semi-broken experimental state.  We will also need to
-investigate handling whatever Pandas produces.
+Time inputs are arrays of time in days relative to the epoch
+given in the `epoch` keyword argument.  In the Matlab version
+of utide these would be Matlab datenums; in the python version,
+using Matlab datenums requires `epoch = 'matlab'`.  The default is
+`epoch = 'python'`, corresponding to the `matplotlib` date
+numbers.  Any other epoch can be specified using either a
+string, like `'2015-01-01'`, or a Python standard library
+`datetime.datetime` or `datetime.date` instance.  The numpy
+`datetime64` dtype is not yet supported, nor are any Pandas
+constructs.
 
 Missing values
 ^^^^^^^^^^^^^^^
-The `t`, `u`, `v` inputs to `solve` now support any combination
+The `t`, `u`, `v` inputs to `solve` and the `t` input to `
+reconstruct` now support any combination
 of nans and masked array inputs to indicate missing values.
 
 The degree to which masked arrays will be used internally is

tests/test_astron.py

@@ -17,7 +17,9 @@
 
 
 def test_astron():
-    dns = [694086.000000000, 736094.552430556]
+    # Switch from Matlab epoch to python epoch.
+    dns = [t - 366 for t in [694086.000000000, 736094.552430556]]
+
     a_expected = np.array([[-0.190238223090586, -0.181487296022524],
                            [0.308043143259513, 0.867328798490917],
                            [0.117804920168928, 0.133410946894728],

tests/test_harmonics.py

@@ -21,6 +21,10 @@
 
 def test_FUV():
     x = loadbunch(fname, masked=False)
+    # Switch epoch from Matlab to Python
+    x.t -= 366
+    x.t0 -= 366
+
     for i, flag in enumerate(x.flags):
         F, U, V = FUV(x.t, x.t0, x.lind-1, x.lat, flag)
         print('i:', i, "ngflgs:", flag)

tests/test_normalize.py

@@ -0,0 +1,20 @@
+from __future__ import (absolute_import, division, print_function)
+
+import datetime
+
+import numpy as np
+from utide._time_conversion import _normalize_time
+
+
+def test_formats():
+    forms = [(np.array([693595.1]), 'python'),
+             (np.array([693961.1]), 'matlab'),
+             (np.array([2.1]), datetime.date(1899, 12, 29)),
+             (np.array([3.1]), datetime.datetime(1899, 12, 28)),
+             (np.array([2.6]), datetime.datetime(1899, 12, 28, 12)),
+             (np.array([2.1]), '1899-12-29')]
+
+    expected = _normalize_time(*forms[0])
+
+    for form in forms[1:]:
+        np.testing.assert_almost_equal(_normalize_time(*form), expected)

tests/test_uv.py

@@ -32,6 +32,7 @@ def fake_tide(t, M2amp, M2phase):
 @pytest.fixture
 def make_data():
     N = 500
+    np.random.seed(1234)
     t = date_range(start='2016-03-29', periods=N, freq='H')
     # Signal + some noise.
     u = fake_tide(np.arange(N), M2amp=2, M2phase=0) + np.random.randn(N)

utide/_reconstruct.py

@@ -3,9 +3,10 @@
 import numpy as np
 from .harmonics import ut_E
 from .utilities import Bunch
+from ._time_conversion import _normalize_time
 
 
-def reconstruct(t, coef, **opts):
+def reconstruct(t, coef, epoch='python', verbose=True, **opts):
     """
     Reconstruct a tidal signal.
 
@@ -15,9 +16,12 @@ def reconstruct(t, coef, **opts):
         Time in days since `epoch`.
     coef : `Bunch`
         Data structure returned by `utide.solve`
-    epoch : {string, int, `datetime.datetime`}, optional
-        Not implemented yet. It will default to the epoch
-        used for `coef`.
+    epoch : {string, `datetime.date`, `datetime.datetime`}, optional
+        Valid strings are 'python' (default); 'matlab' if `t` is
+        an array of Matlab datenums; or an arbitrary date in the
+        form 'YYYY-MM-DD'.  The default corresponds to the Python
+        standard library `datetime` proleptic Gregorian calendar,
+        starting with 1 on January 1 of year 1.
     verbose : {True, False}, optional
         True to enable output message (default). False turns off all
         messages.
@@ -31,7 +35,7 @@ def reconstruct(t, coef, **opts):
     """
 
     out = Bunch()
-    u, v = _reconstr1(t, coef, **opts)
+    u, v = _reconstr1(t, coef, epoch=epoch, verbose=verbose, **opts)
     if coef['aux']['opt']['twodim']:
         out.u, out.v = u, v
     else:
@@ -42,7 +46,7 @@ def reconstruct(t, coef, **opts):
 def _reconstr1(tin, coef, **opts):
 
     # Parse inputs and options.
-    t, opt = _rcninit(tin, **opts)
+    t, goodmask, opt = _rcninit(tin, **opts)
 
     if opt['RunTimeDisp']:
         print('reconstruct:', end='')
@@ -108,7 +112,7 @@ def _reconstr1(tin, coef, **opts):
 
     # Mean (& trend).
     u = np.nan * np.ones(tin.shape)
-    whr = ~np.isnan(tin)
+    whr = goodmask
     if coef['aux']['opt']['twodim']:
         v = np.nan * np.ones(tin.shape)
         if coef['aux']['opt']['notrend']:
@@ -127,7 +131,7 @@ def _reconstr1(tin, coef, **opts):
             u[whr] = np.real(fit) + coef['mean']
             u[whr] += coef['slope'] * (t-coef['aux']['reftime'])
 
-        v = []
+        v = None
 
     if opt['RunTimeDisp']:
         print('done.')
@@ -139,14 +143,24 @@ def _rcninit(tin, **opts):
 
     t = tin[:]
 
-    t[np.isnan(t)] = []
-    # t(isnan(t)) = []
+    # Supporting only 1-D arrays for now; we can add "group"
+    # support later.
+    if tin.ndim != 1:
+        raise ValueError("t must be a 1-D array")
+
+    # Step 0: apply epoch to time.
+    t = _normalize_time(tin, opts['epoch'])
+
+    # Step 1: remove invalid times from tin
+    t = np.ma.masked_invalid(t)
+    goodmask = ~np.ma.getmaskarray(t)
+    t = t.compressed()
+
     opt = {}
 
     opt['cnstit'] = False
     opt['minsnr'] = 2
     opt['minpe'] = 0
-    opt['RunTimeDisp'] = True
 
     for key, item in opts.items():
         # Be backward compatible with the MATLAB package syntax.
@@ -158,25 +172,4 @@ def _rcninit(tin, **opts):
         except KeyError:
             print('reconstruct: unrecognized input: {0}'.format(key))
 
-    # args = list(args)
-    # args = [string.lower() for string in args]
-
-    # Need an example of the args
-
-#    while ~isempty(args)
-#        switch(lower(args{1}))
-#            case 'cnstit'
-#                opt.cnstit = args{2};
-#                args(1:2) = [];
-#            case 'minsnr'
-#                opt.minsnr = args{2};
-#                args(1:2) = [];
-#            case 'minpe'
-#                opt.minpe = args{2};
-#                args(1:2) = [];
-#            otherwise
-#                error(['reconstruct: unrecognized input: ' args{1}]);
-#        end
-#    end
-
-    return t, opt
+    return t, goodmask, opt

utide/_solve.py

@@ -14,6 +14,7 @@
 from .utilities import Bunch
 from . import constit_index_dict
 from .robustfit import robustfit
+from ._time_conversion import _normalize_time
 
 default_opts = dict(constit='auto',
                     conf_int='linear',
@@ -26,7 +27,8 @@
                     Rayleigh_min=1,
                     robust_kw=dict(weight_function='cauchy'),
                     white=False,
-                    verbose=True
+                    verbose=True,
+                    epoch='python',
                     )
 
 
@@ -71,6 +73,7 @@ def _translate_opts(opts):
     oldopts.white = opts.white
     oldopts.newopts = opts  # So we can access new opts via the single "opt."
     oldopts['RunTimeDisp'] = opts.verbose
+    oldopts.epoch = opts.epoch
     return oldopts
 
 
@@ -88,8 +91,12 @@ def solve(t, u, v=None, lat=None, **opts):
         If `u` is a velocity component, `v` is the orthogonal component.
     lat : float
         Latitude in degrees; required.
-    epoch : {string, int, float, `datetime.datetime`}
-        Not implemented yet.
+    epoch : {string, `datetime.date`, `datetime.datetime`}, optional
+        Valid strings are 'python' (default); 'matlab' if `t` is
+        an array of Matlab datenums; or an arbitrary date in the
+        form 'YYYY-MM-DD'.  The default corresponds to the Python
+        standard library `datetime` proleptic Gregorian calendar,
+        starting with 1 on January 1 of year 1.
     constit : {'auto', array_like}, optional
         List of strings with standard letter abbreviations of
         tidal constituents; or 'auto' to let the list be determined
@@ -131,7 +138,7 @@ def solve(t, u, v=None, lat=None, **opts):
         residuals in the confidence limit estimates; if True,
         assume a white background spectrum.
     verbose : {True, False}, optional
-        True (default) turns on verbose output. False omits no messages.
+        True (default) turns on verbose output. False emits no messages.
 
     Note
     ----
@@ -336,6 +343,9 @@ def _slvinit(tin, uin, vin, lat, **opts):
 
     opt = Bunch(twodim=(vin is not None))
 
+    # Step 0: apply epoch to time.
+    tin = _normalize_time(tin, opts['epoch'])
+
     # Step 1: remove invalid times from tin, uin, vin
     tin = np.ma.masked_invalid(tin)
     uin = np.ma.masked_invalid(uin)

utide/_time_conversion.py

@@ -0,0 +1,39 @@
+"""
+Utility for allowing flexible time input.
+"""
+
+from __future__ import (absolute_import, division, print_function)
+
+import warnings
+from datetime import date, datetime
+
+try:
+    from datetime import timezone
+    have_tz = True
+except ImportError:
+    have_tz = False
+
+
+def _normalize_time(t, epoch):
+    if epoch == 'python':
+        return t
+    if epoch == 'matlab':
+        return t - 366
+    try:
+        epoch = datetime.strptime(epoch, "%Y-%m-%d")
+    except (TypeError, ValueError):
+        pass
+    if isinstance(epoch, date):
+        if not hasattr(epoch, 'time'):
+            return t + epoch.toordinal()
+        # It must be a datetime, which is also an instance of date.
+        if epoch.tzinfo is not None:
+            if have_tz:
+                epoch = epoch.astimezone(timezone.utc)
+            else:
+                warnings.warn("Timezone info in epoch is being ignored;"
+                              " UTC is assumed.")
+        ofs = (epoch.toordinal() + epoch.hour / 24 +
+               epoch.minute / 1440 + epoch.second / 86400)
+        return t + ofs
+    raise ValueError("Cannot parse epoch as string or date or datetime")

utide/astronomy.py

@@ -21,7 +21,9 @@ def ut_astron(jd):
     Parameters
     ----------
     jd : float, scalar or sequence
-        Time (UTC) in days since the Matlab datenum epoch.
+        Time (UTC) in days starting with 1 on 1 Jan. of the year 1
+        in the proleptic Gregorian calendar as in
+        `datetime.date.toordinal`.
 
     Returns
     -------
@@ -56,8 +58,10 @@ def ut_astron(jd):
 
     jd = np.atleast_1d(jd).flatten()
 
-    # datenum(1899,12,31,12,0,0)
-    daten = 693961.500000000
+    # Shift epoch to 1899-12-31 at noon:
+    # daten = 693961.500000000  Matlab datenum version
+
+    daten = 693595.5  # Python epoch is 366 days later than Matlab's
 
     d = jd - daten
     D = d / 10000

utide/confidence.py

@@ -134,7 +134,8 @@ def nearestSPD(A):
         # Tweaking strategy differs from D'Errico version.  It
         # is still a very small adjustment, but much larger than
         # his.
-        maxeig = np.linalg.eigvals(Ahat).max()
+        # Eigvals are or can be complex dtype, so take abs().
+        maxeig = np.abs(np.linalg.eigvals(Ahat)).max()
         Ahat[np.diag_indices(n)] += np.spacing(maxeig)
         # Normally no more than one adjustment will be needed.
         if k > 100:

utide/constituent_selection.py

@@ -41,7 +41,7 @@ def ut_cnstitsel(tref, minres, incnstit, infer):
     astro, ader = ut_astron(tref)
 
     ii = np.isfinite(const.ishallow)
-    const.freq[~ii] = np.dot(const.doodson[~ii, :], ader) / 24
+    const.freq[~ii] = np.dot(const.doodson[~ii, :], ader[:, 0]) / 24
 
     for k in ii.nonzero()[0]:
         ik = const.ishallow[k]+np.arange(const.nshallow[k])