dankelley / oce

R package for oceanographic processing
http://dankelley.github.io/oce/
GNU General Public License v3.0
144 stars 42 forks source link

tidem() and t_tide do not produce identical results #1653

Open richardsc opened 4 years ago

richardsc commented 4 years ago

There's a lot of info in #1621 , which lead to the following exploration of tidem() and t_tide on a particular dataset (but there's not reason to think that the results are specific to that dataset):

compare_t_tide_with_tidem.pdf

Below is the result of a 10 constituent fit using tidem(): image

and below the same but for t_tide: image image

Note that the fit amplitudes are all quite similar, as are the phases -- except for the M4 phase (200.1 vs 177.15).

Why is this?

dankelley commented 4 years ago

I'll look into this. Please note that I created https://github.com/dankelley/oce-issues/tree/master/16xx/1653 and put the dataset there, along with some code I wrote to test, and your .pdf. If you could drop your Rmd into that directory, it might be helpful to me.

The fact that only one constituent is affected tells me that it might be some kind of satellite / inference thing relating to grouped constituents. This is explained in the Foreman documents, referred to in ?tidem ... which I'll need to reread, I think.

Since I'll need to study the Foreman documents, and the tricky code that does some of this satellite/inference work, my guess is that this issue will take a while. However, your very detailed report will make that work a lot easier, so thanks for that!!

Also, Clark, if you could put the matlab output file in the directory, that would be great, too.

richardsc commented 4 years ago

I'll post the Rmd and the matlab output -- likely tomorrow when I'm back in the office.

richardsc commented 4 years ago

Here is the Rmd:

compare_t_tide_with_tidem.zip

And a zip with all the files required (Matlab scripts and saved files):

compare_t_tide_with_tidem.zip

dankelley commented 4 years ago

Further to our f2f regarding possible problems with phase inference for small-amplitude constituents, the code

## Test whether phase drifts for small-amplitude constituents
library(oce)
data(tidedata)
dt <- 15 * 60                          # 15 minutes
t <- seq(as.POSIXct(Sys.time()), length.out=30*86400/dt, by=dt)
tt <- as.numeric(t)
## frequencies are in cycles/hour
M2 <- tidedata$const$freq[tidedata$const$name == "M2"] / 3600 * (2 * pi)
S2 <- tidedata$const$freq[tidedata$const$name == "S2"] / 3600 * (2 * pi)
M4 <- tidedata$const$freq[tidedata$const$name == "M4"] / 3600 * (2 * pi)
aM2 <- 1
aS2 <- 0.5
for (aM4 in 10^seq(1, -5)) {
    eta <- aM2*sin(M2*tt) + aM2*cos(M2*tt) + aS2*sin(S2*tt) + aS2*cos(S2*tt) + aM4*sin(M4*tt) + aM4*cos(M4*tt)
    sl <- as.sealevel(eta, t)
    capture.output(m <- tidem(sl))
    ## summary(m)
    phase <- m[["phase"]]
    name <- m[["name"]]
    cat(sprintf("aM4=%.5f yields phase=%.3f\n", aM4, phase[name=="S4"]))
}

produces

aM4=10.00000 yields phase=349.673
aM4=1.00000 yields phase=174.271
aM4=0.10000 yields phase=171.876
aM4=0.01000 yields phase=171.806
aM4=0.00100 yields phase=171.826
aM4=0.00010 yields phase=171.777
aM4=0.00001 yields phase=171.795

I'm not sure I understand this, but thought it worth posting the results, to have them on the record.

dankelley commented 4 years ago

I dumped the CR zipfile into https://github.com/dankelley/oce-issues/tree/master/16xx/1653/cr and excerpted the table comparing results into https://github.com/dankelley/oce-issues/blob/master/16xx/1653/dk/1653dk_01.R so I could look in more detail. That R script has enough comments to explain, so I won't replicate what I did in words here.

The output of https://github.com/dankelley/oce-issues/blob/master/16xx/1653/dk/1653dk_01.R is at Appendix 1 below. There's a lot I could look at, but my eye was drawn to the first big phase difference at MO3. I see big phase differences nearby (in freq) at MK3 and SK3.

Now, take a look at Appendix 2, which is a table (figure, really) from Foreman (1981). This deals with frequencies that are near each other. This time-series is 768 hours long, so the bottom row of the table is irrelevant (i.e. the Rayleigh criterion means we do not fit for SO3). But look at the others. Note that M3 is related (in this sense) to MO3, MK3 and SK3.

And (strokes 2-day old beard) this triplet, MO3, MK3 and SK3 have phase differences in the hundreds of degrees.

I need to do some more reading on this topic. (I realize that I am not explaining this well here, but this is something @richardsc and I were discussing at our last f2f, so I think there's no real need to inarticulately state something here that is probably very clear in Foreman's document.)

References

Appendix 1

    Name     Freq    msnr phaseMismatch timeShiftHours
1     MM 0.001512    1.70          0.03          0.055
2    MSF 0.002822    0.41          0.02          0.020
3  *ALP1 0.034397    2.80          0.27          0.022
4    2Q1 0.035706    1.30         -0.30         -0.023
5    *Q1 0.037219    5.00          0.10          0.007
6    *O1 0.038731  150.00         -0.03         -0.002
7    NO1 0.040269    1.50          0.48          0.033
8    *K1 0.041781  220.00          0.06          0.004
9     J1 0.043293    0.87         -0.66         -0.042
10   OO1 0.044831    0.68          0.48          0.030
11  UPS1 0.046343    0.46          1.00          0.060
12 *EPS2 0.076177    2.40          0.09          0.003
13  *MU2 0.077689    4.60         -0.07         -0.003
14   *N2 0.078999  110.00          0.00          0.000
15   *M2 0.080511 1300.00         -0.03         -0.001
16   *L2 0.082024   15.00         -0.02         -0.001
17   *S2 0.083333   38.00         -0.04         -0.001
18  ETA2 0.085074    1.10         -0.22         -0.007
19  *MO3 0.119242    5.40        271.58          6.327
20    M3 0.120767    0.85          0.66          0.015
21  *MK3 0.122292    3.50        122.52          2.783
22   SK3 0.125114    0.12        109.09          2.422
23  *MN4 0.159511    2.60         23.88          0.416
24   *M4 0.161023   14.00         22.90          0.395
25  *SN4 0.162333    2.60         12.82          0.219
26   MS4 0.163845    1.60         11.58          0.196
27    S4 0.166667    0.52          0.56          0.009
28  2MK5 0.202804    1.00        135.69          1.859
29 *2SK5 0.208447    2.90       -249.15         -3.320
30  2MN6 0.240022    1.10         35.18          0.407
31   *M6 0.241534    4.00         34.47          0.396
32 *2MS6 0.244356    2.40         23.15          0.263
33  2SM6 0.247178    0.49         12.43          0.140
34  3MK7 0.283315    0.26       -216.01         -2.118
35   *M8 0.322046    2.60         45.74          0.395

Appendix 2

foreman1981_table4