UTide icon indicating copy to clipboard operation
UTide copied to clipboard

lat=0 issue:

Open apatlpo opened this issue 6 years ago • 2 comments

Proving lat=0 leads to a crash with the following example:

    import utide
    import numpy as np
    import matplotlib.dates as mdates
    import datetime, dateutil

    print(utide.__version__)

    time = np.arange(0.,30.*86400,60.)
    t0 = datetime.datetime(2011,9,10)
    dtime = [t0+dateutil.relativedelta.relativedelta(seconds=t) for t in time]
    mtime = mdates.date2num(dtime)

    tdat = utide.ut_constants['const']
    iM2 = tdat['name'].tolist().index('M2')
    omega = tdat['freq'][iM2]
    
    eta = 1.*np.cos(2.*np.pi*omega*time/3600.)

    coef = utide.solve(mtime, eta, lat=0., method='ols', conf_int='MC')

With the following error message:

(equinox) r2i2n0 datahome/aponte% python tests.py
0.2.2
/home1/datahome/aponte/.miniconda3/envs/equinox/lib/python3.6/site-packages/utide/harmonics.py:132: RuntimeWarning: divide by zero encountered in double_scalars
  rr[j] *= 0.36309 * (1.0 - 5.0 * slat**2)/slat
/home1/datahome/aponte/.miniconda3/envs/equinox/lib/python3.6/site-packages/numpy/core/_methods.py:32: RuntimeWarning: invalid value encountered in reduce
  return umr_sum(a, axis, dtype, out, keepdims)
solve: matrix prep ... solution ... Traceback (most recent call last):
  File "tests.py", line 57, in <module>
    test_utide()
  File "tests.py", line 46, in test_utide
    coef = utide.solve(mtime, eta, lat=0., method='ols', conf_int='MC')
  File "/home1/datahome/aponte/.miniconda3/envs/equinox/lib/python3.6/site-packages/utide/_solve.py", line 200, in solve
    coef = _solv1(t, u, v, lat, **compat_opts)
  File "/home1/datahome/aponte/.miniconda3/envs/equinox/lib/python3.6/site-packages/utide/_solve.py", line 291, in _solv1
    m = np.linalg.lstsq(B, xraw, rcond=None)[0]
  File "/home1/datahome/aponte/.miniconda3/envs/equinox/lib/python3.6/site-packages/numpy/linalg/linalg.py", line 2031, in lstsq
    0, work, lwork, rwork, iwork, 0)
ValueError: On entry to DLASCL parameter number 4 had an illegal value
(equinox) r2i2n0 datahome/aponte%

apatlpo avatar Jun 21 '18 13:06 apatlpo

Substituying https://github.com/wesleybowman/UTide/blob/82fd9b79e7968f1ad30cf6fce7fdb7fdaa11f53f/utide/harmonics.py#L125 by:

if abs(lat) == 0:
      lat = 5
elif abs(lat) < 5:
      lat = np.sign(lat) * 5

apatlpo avatar Jun 21 '18 13:06 apatlpo

I can push a PR if you want me to

apatlpo avatar Jun 21 '18 13:06 apatlpo