Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fixes for 'lat=0' and 0 time series issues #62

Open
wants to merge 4 commits into
base: master
Choose a base branch
from

Conversation

apatlpo
Copy link

@apatlpo apatlpo commented Jun 21, 2018

No description provided.

Duu = Puu[c] * Duu / np.trace(Duu)
Duu = Puu[c] * Duu
if np.trace(Duu) != 0. :
Duu = Duu / np.trace(Duu)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't compute np.trace twice and we can do the Duu division in place:

trace = np.trace(Duu)
if trace != 0. :
    Duu /= trace

if abs(lat) < 5:
if lat == 0:
lat = 5
elif abs(lat) < 5:
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM

@ocefpaf
Copy link
Collaborator

ocefpaf commented Jun 21, 2018

xref.: #60 and #61

@@ -252,7 +252,9 @@ def _confidence(coef, cnstit, opt, t, e, tin, elor, xraw, xmod, W, m, B,
varcov_mCw[c, :2, :2] = Duu

if not opt.white:
Duu = Puu[c] * Duu / np.trace(Duu)
Duu = Puu[c] * Duu
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is a trailing space after Duu.

@@ -252,7 +252,9 @@ def _confidence(coef, cnstit, opt, t, e, tin, elor, xraw, xmod, W, m, B,
varcov_mCw[c, :2, :2] = Duu

if not opt.white:
Duu = Puu[c] * Duu / np.trace(Duu)
Duu = Puu[c] * Duu
if np.trace(Duu) != 0. :
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is a whitespace before :.

Duu = Puu[c] * Duu
trace = np.trace(Duu)
if trace != 0.:
Duu = Duu / trace
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Any reason to not use the compact in place Duu /= trace?

@ocefpaf
Copy link
Collaborator

ocefpaf commented Jun 21, 2018

I'll leave this here so @efiring and @wesleybowman can take a second look before merging.

Thanks @apatlpo for the PR!

@efiring
Copy link
Collaborator

efiring commented Jun 21, 2018

Two real problems have been identified, but based on a quick look, I'm not sure the solutions are good ones.

  • For the zero-data situation, the solution might be OK; what is its net result? Does it leave zero-confidence limits, as it should? Are all outputs as expected in this pathological zero-data case?
  • For the zero-latitude situation, the original low-lat code doesn't make sense to me, so I suspect a more fundamental fix may be in order. How can it be correct to pin low lats to +-5? Especially when a blow-up at lat==0 suggests sensitivity to low latitudes?

An oddity related to the latter point is that the code is using conditional statements involving only values in constants, which are of course constant, so no conditionals should be needed.

@efiring
Copy link
Collaborator

efiring commented Jun 22, 2018

I think one of my potential objections is invalid, but the main one--the odd behavior within 5 degrees of the equator--still needs attention. I see that the behavior is inherited from t_tide, which includes some explanatory comments but still doesn't make sense to me. It implies a discontinuity at the equator.

@apatlpo
Copy link
Author

apatlpo commented Jun 23, 2018

Here are the explanatory comments:

  if nargin==4, % If we have a latitude, get nodal corrections.

    % Apparently the second-order terms in the tidal potential go to zero
    % at the equator, but the third-order terms do not. Hence when trying
    % to infer the third-order terms from the second-order terms, the 
    % nodal correction factors blow up. In order to prevent this, it is 
    % assumed that the equatorial forcing is due to second-order forcing 
    % OFF the equator, from about the 5 degree location. Latitudes are 
    % hence (somewhat arbitrarily) forced to be no closer than 5 deg to 
    % the equator, as per note in Foreman.

    if isfinite(lat) & (abs(lat)<5); lat=sign(lat).*5; end

    slat=sin(pi.*lat./180);
    % Satellite amplitude ratio adjustment for latitude. 

Maybe there is more in Foreman's code

It may be more simple to ask also experts.
I am attending a meeting next week and Richard Ray should be there, i'll ask him about this.

@apatlpo
Copy link
Author

apatlpo commented Jul 2, 2018

After a brief discussion and a lengthier email, Richard Ray's opinion is that ''bringing latitude dependence in the calculation of nodal modulations, as well as in the inference of minor tides is not a good idea''.

I asked him for a good reference in order to know what we are supposed to do instead, here is his reply:

Section 4.2.2 of the Pugh-Woodworth 2014 sea level book is ok.
(plus 4.2.3 on shallow water tides).

I've always liked Doodson & Warburg (1941) - Admiralty Manual of Tides,
which a few years ago was still in print.

@efiring
Copy link
Collaborator

efiring commented Jul 2, 2018

@apatlpo Thank you, I hope he can provide a reference we can work with.

Copy link
Collaborator

@efiring efiring left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

  1. Is it clear that the problem of zero data occurs only in the non-white-noise case, so this is the ideal place to intercept it? Regardless, I think it might be best to just raise an exception when obviously useless data is provided. A sanity check could be done early on, as an alternative to checking the trace of Duu.
  2. There must be a better way to handle the singularity at lat=0. We just have to find it. It's possible a horrible ad-hoc fix something like the one here will be necessary temporarily. Rather than having a jump at the equator, I suspect it would make more sense to have a linear transition in rr values from 5N to 5S. But I don't understand the math and physics behind rr right now, so I don't know how to handle it.

@@ -252,7 +252,10 @@ def _confidence(coef, cnstit, opt, t, e, tin, elor, xraw, xmod, W, m, B,
varcov_mCw[c, :2, :2] = Duu

if not opt.white:
Duu = Puu[c] * Duu / np.trace(Duu)
Duu = Puu[c] * Duu
trace = np.trace(Duu)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is not the same as the original even in the normal case. The order of the two lines above would need to be swapped.

@efiring
Copy link
Collaborator

efiring commented Sep 9, 2020

In a week or two I expect to be able to have a look at the Pugh-Woodworth book recommended by Richard Ray; I hope it will have a clear explanation of a replacement for the problematic calculation we have now. @DanCodigaMWRA, do you still work with tides? If so, do you have some ideas about this? Or know someone with the expertise and willingness to help out?

@apatlpo
Copy link
Author

apatlpo commented Sep 9, 2020

I don't have the time unfortunately to dive into this, sorry about this.

I would like to point out the pangeo-pytide library though.
A solution to the problem at hand may have been found there.

@DanCodiga
Copy link

In a week or two I expect to be able to have a look at the Pugh-Woodworth book recommended by Richard Ray; I hope it will have a clear explanation of a replacement for the problematic calculation we have now. @DanCodigaMWRA, do you still work with tides? If so, do you have some ideas about this? Or know someone with the expertise and willingness to help out?

@efiring thanks for asking. I do still work on tides a bit. Mostly just helping folks who are using the Matlab UTide functions.

As for the specific topic here, latitudes b/w -5 and +5, I don't have much insight and would have to dig around to try to come up to speed. The quoted documentation from t-tide shown above seems to attribute it to Foreman. He was helpful to me back when I was putting UTide together, would be worth getting his input if possible.

p.s. When working on tides I am @DanCodiga (not @DanCodigaMWRA).

p.s.s. I do aspire to revamp UTide as it has been almost 10 years now. I will probably start with a survey of all users to learn what would be the most helpful next step, i.e. update the Matlab version or the python version (or build an R version?); help you folks improve the python version here; add new features; add better tutorial/documentation; etc. (If it's python-focused I would see it as a chance for me to finally get up to speed with using github too...)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants