UTide
UTide copied to clipboard
lat=0 issue:
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%
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
I can push a PR if you want me to