Issue with tidem function

[Pra713]

Short summary of problem

Running tidem function with constituents gives me an error message when I include "UPS1" or exclude "Z0" constituents.

Details (optional)

Despite my using "UPS1" as a constituent, it gives me an error saying that "UPSI" is not a valid constituent. Despite my trying to exclude "Z0" as a constituent, it gives me an error saying that "Z0" is not a valid constituent. They are both valid constituents accepted by the function. In the case of "UPS1", the error makes it seem like I'm inputting "UPSI" which I am not, and actually suggest "UPS1" as a constituent among others.

What you did

I downloaded tide data from NOAA website, https://tidesandcurrents.noaa.gov/. After some preprocessing, I tried running tidem for selected constituents. I have attached a link to my Google drive folder with the code as well as the preprocessed data here:
https://drive.google.com/drive/folders/1YAMLGC3_J1wATueLmlQz4Wlk9Ja4NR9S?usp=sharing

The code (comments removed due to formatting issue, but the code in link is fully commented):

library(oce)
df_data <- readRDS('df_data.RDS')
constituents_29 <- c("MSF", "2Q1", "Q1",  "O1",  "NO1", "K1",  "J1",  "OO1", "UPS1", "N2",  "M2",  "S2",  "ETA2", "MO3", "M3",  "MK3", "SK3", "MN4", "M4",  "MS4", "S4",  "2MK5", "2SK5", "2MN6", "M6",  "2MS6", "2SM6", "3MK7", "M8")
t_result <- tidem(df_data$datetime, df_data$verified_m, constituents = constituents_29)
constituents_exlude <- c("standard", "-Z0", "-SA", "-SSA", "-MSM", "-MM", "-MF", "-ALP1", "-SIG1", "-RHO1", "-TAU1", "-BET1", "-CHI1", "-PI1", "-P1", "-S1", "-PSI1", "-PHI1", "-THE1", "-SO1", "-OQ2", "-EPS2", "-2N2", "-MU2", "-NU2", "-GAM2", "-H1", "-H2", "-MKS2", "-LDA2", "-L2", "-T2", "-R2", "-K2", "-MSN2", "-SO3", "-SN4", "-MK4", "-SK4", "-2MK6", "-MSK6")
t_result <- tidem(df_data$datetime, df_data$verified_m, constituents = constituents_exlude)

What you expected to happen

The code to run and give me a fitted tide model.

What happened

I instead got the following error when I tried to include "UPS1" as a constituent:
"'UPSI' is not a known tidal constituent"

I got the following error when I tried to exclude "Z0" as a contituent:
"Z0is not in the list of constituents currently under study"

How urgent is this?

Yes, this error can wait a few days.

Output from sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] oce_1.1-1      gsw_1.0-5      testthat_2.2.1

loaded via a namespace (and not attached):
[1] compiler_3.6.1 magrittr_1.5   R6_2.4.0       tools_3.6.1    yaml_2.2.0     Rcpp_1.0.2    
[7] rlang_0.4.0   

[richardsc]

Hi @Pra713

Thanks for opening the issue, and nice job producing a clear and reproducible issue report!

@dankelley is better in the weeds of the tidal fit stuff than I am, so I'll let him comment too, but a couple thoughts:

• Typically the approach for doing tidal harmonic analysis is to let the algorithm figure out what components are appropriate to try fitting. For example, depending on how long your time series is there may be components that can't be fit (period is too long). Also, the fitting procedure can figure out which components have a statistically significant fit, and identifies those.

• Not sure why UPS1 exists in t_tide and not in tidem(), as they are supposedly both based on Foreman 1978. But anyway the source of that error is easy to fix -- don't try and include a component that doesn't exist! 😄 (@dankelley will weigh in on this I'm sure).

• As for excluding Z0, I suspect that you can't exclude it because it isn't actually a harmonic component. It's essentially just the offset from zero for the elevation time series.

[Pra713]

Hi Richard,

Thanks for the detailed response.

• While that is what I would do for research purposes, I am currently taking a course for grad school where we use MATLAB's t_tide for the analysis. But I prefer to code in R. And for the sake of consistency, I wanted to make sure the constituents matched the ones used by t_tide.

• tidem() does have UPS1. But for some reason, when I input UPS1, it seems to wrongly read that I am inputting UPSI and gives the error. Or I might be missing something. Here is the full error message I get when I include UPS1 as a constituent (relevant parts bolded):

Error in tidem(df_data$datetime, df_data$verified_m, constituents = constituents_29) :
'UPSI' is not a known tidal constituent; try one of: Z0 SA SSA MSM MM MSF MF ALP1 2Q1 SIG1 Q1 RHO1 O1 TAU1 BET1 NO1 CHI1 PI1 P1 S1 K1 PSI1 PHI1 THE1 J1 2PO1 SO1 OO1 UPS1 ST36 2NS2 ST37 ST1 OQ2 EPS2 ST2 ST3 O2 2N2 MU2 SNK2 N2 NU2 ST4 OP2 GAM2 H1 M2 H2 MKS2 ST5 ST6 LDA2 L2 2SK2 T2 S2 R2 K2 MSN2 ETA2 ST7 2SM2 ST38 SKM2 2SN2 NO3 MO3 M3 NK3 SO3 MK3 SP3 SK3 ST8 N4 3MS4 ST39 MN4 ST9 ST40 M4 ST10 SN4 KN4 MS4 MK4 SL4 S4 SK4 MNO5 2MO5 3MP5 MNK5 2MP5 2MK5 MSK5 3KM5 2SK5 ST11 2NM6 ST12 2MN6 ST13 ST41 M6 MSN6 MKN6 ST42 2MS6 2MK6 NSK6 2SM6 MSK6 S6 ST14 ST15 M7 ST16 3MK7 ST17 ST18 3MN8 ST19 M8 ST20 ST21 3MS8 3MK8 ST22 ST23 ST24 ST25 ST26 4MK9 ST27 ST28 M10 ST29 ST30 ST31 ST32 ST33 M12 ST34 ST35
In addition: Warning message:
In constituentNameFix(constituents) :
constituent name switched from T-TIDE 'UPSI' to Foreman (1977) 'UPS1'

• That is good to know and makes perfect sense. One suggestion I would, however, give, is that it would have been helpful if this was provided in the help documentation for the function (?tidem).

[dankelley]

Hi. Sorry to be coming to this late in the day (um, early in the day, in my timezone). This is documented fairly clearly in the output from ?tidem, e.g. the screenshots:

Screenshot 1

Screenshot 2

[dankelley]

Re

I got the following error when I tried to exclude "Z0" as a contituent:
"Z0is not in the list of constituents currently under study"

I'll fix that up sometime today or tomorrow, by checking specifically whether Z0 is in the exclusion list, and issuing a more helpful warning if so. I've never seen negative numbers reported by a tide gauge, so I think removing it could have weird consequences for real-life data, but I won't make it an error, because someone might use tidem() on fake data, and those fake data might not have a constant offset from zero.

[Pra713]

I've never seen negative numbers reported by a tide gauge, so I think removing it could have weird consequences for real-life data, but I won't make it an error, because someone might use tidem() on fake data, and those fake data might not have a constant offset from zero.

I got this data from the NOAA tides and currents website. I believe that the water level data can be positive or negative since it simply depends on the datum chosen. So while the data reported by a tide gauge may never be negative, it can be negative when adjusted for some datums (eg. MHW, MSL, etc).

Screenshot (water level plot from NOAA Tides website with MSL as datum):

[dankelley]

Interesting, re the datum in the noaa dataset. Thanks.

[Pra713]

Screenshot 2

I have tried both UPS1 and UPSI as an input to the constituents to get the exact same error. UPS1 is a valid constituent in tidem() since running the code with standard constituents gives the fitting parameters for UPS1.

Screenshot 1 (using UPSI):

Screenshot 2 (using UPS1):

Screenshot 3 (using standard constituents; note that UPS1 is among the fitted constituents):

[dankelley]

I'll look at this over the weekend.

[Pra713]

Thanks!

[Pra713]

I'll fix that up sometime today or tomorrow, by checking specifically whether Z0 is in the exclusion list, and issuing a more helpful warning if so. I've never seen negative numbers reported by a tide gauge, so I think removing it could have weird consequences for real-life data, but I won't make it an error, because someone might use tidem() on fake data, and those fake data might not have a constant offset from zero.

And yes, the datum,Z0, is a fixed known value (in my study at least). So, allowing this parameter to be excluded from the constituents list in tidem() would be ideal.

[dankelley]

This may work now, in the "develop" branch. See the commit note, pasted below. However, I do not know whether the results are okay for the test contributed by @Pra713 ... all I know is that the fairly extensive tests in https://github.com/dankelley/oce/blob/develop/tests/testthat/test_tidem.R still pass, with this change.

@Pra713 -- if you have T-TIDE installed, perhaps you can do a test to see whether tidem() recovers the amplitude and phase. If not, I'd appreciate a comment that gives me 10 digits on both of these, for the T-TIDE case, so I can use that as a target for further tests.

Thanks for your patience on this. Last week was busy with teaching work.

commit 57b9ec1
Author: dankelley [email protected]
Date: Fri Nov 8 15:00:51 2019 -0400

tidem() improvements for issue 1621

See https://github.com/dankelley/oce/issues/1621 and note:

1. The renaming from T-TIDE was wrong (reversed) for "UPS1".

2. Now, we have a new user-visible function named
tidemConstituentNameFix() that does renaming. This will make it easier
for users to debug problems, and also submit patches (if more naming
differences are uncovered).

3. The internal functioning of Z0 is changed. The tests in
tests/testthat/test_tidem.R still work properly, which means that we are
still reproducing the T_TIDE results (to the extent that we check them).
Also, those tests require that amplitude and phase match values from
previous versions of oce, so we know that this is okay also.
Finally, there are tests that predict() works sensibly in that file.
Given all of this, I am fairly confident that things are OK, at least
for tidem() and predict().  (I have not checked into as.tidem() yet,
though.) The only thing that might affect a user (well, at least a usr
who is doing unrecommended investigation of the internals of a tidem
object) is that the cos() terms are stored before the sin() terms, in
building the columns of the regression matrix. However, the form of that
matrix is not documented, and so any user that is looking that deep into
the data object is working in unrecommended territory. Sticking to the
interface, that is predict() and summary(), ought to be okay, according
to the tests.

[Pra713]

This may work now, in the "develop" branch. See the commit note, pasted below. However, I do not know whether the results are okay for the test contributed by @Pra713 ... all I know is that the fairly extensive tests in https://github.com/dankelley/oce/blob/develop/tests/testthat/test_tidem.R still pass, with this change.

@Pra713 -- if you have T-TIDE installed, perhaps you can do a test to see whether tidem() recovers the amplitude and phase. If not, I'd appreciate a comment that gives me 10 digits on both of these, for the T-TIDE case, so I can use that as a target for further tests.

I am fairly busy today with my coursework and TA work, but I'll be sure to check it in a day or two. I can compare the results from T-TIDE and tidem() to see if there are any discrepancies. I will post the results here as soon as I have them.

[Pra713]

I'm sorry but I do not have much experience with Github, and installing packages from Github. I seem to be unable to install oce from github using the code provided. I tried uninstalling the previous version of oce, but still the install_github functions throws an error.

library(devtools)
install_github('dankelley/oce', ref='develop')

The error I get is:
Error: Failed to install 'oce' from GitHub:
(converted from warning) cannot remove prior installation of package ‘digest’

I do not have 'digest' package installed.

Using the 'oce' package from CRAN still throws the same error messages.

[Pra713]

Also, I managed to go into the source code form tidem and fix the 'UPSI' issue. Running the code gives me different amplitude and phase values compared to 't_tide'. They are very close, but slightly different. As a result, the output predicted water-level values are slightly different too. I'm assuming this is because 'Z0' amplitude is one of the parameters in tidem while it is 0 in t_tide.

[dankelley]

1. I am not sure what that 'digest' error is about. Sorry, but it's hard to advise on building, apart from saying that I've not heard of that error before.

2. The oce on CRAN will definitely still produce the old error messages. CRAN does not look kindly on frequent updates, and we would definitely not update it soon after adding a new feature like this.

3. I don't follow the statement about Z0 being 0 in T_TIDE. It might help if you could post more details.

Some other notes:

1. We have tests in the build test suite that still pass, which suggests that the tidem matches T_TIDE in some published test cases (see

   oce/tests/testthat/test_tidem.R

   Line 160 in bc31bfc

   test_that("Foreman (1977 App 7.3) and T-TIDE (Pawlowciz 2002 Table 1) test",

   ) and also code that checks whether the new version of tidem matches the old version. However, this does not argue strongly that tidem is "right" (see notes in the test suite about whether T_TIDE is "right" compared with Foreman...), though.

2. @richardsc might have some advice on building, and perhaps also on T_TIDE.

[richardsc]

I also don't know what the digest issue is about. A quick google turns up some similar issues:

https://stackoverflow.com/questions/36533880/cannot-install-flexdashboard-package/36533946#36533946

Maybe try some of what is in there and see if it helps? Particularly making sure that you try the install from a "fresh" R session (i.e. one that doesn't have other packages loaded).

Failing that, Dan can provide you with a pre-compiled binary (for Windows or OSX) that you can install directly.

It would help greatly if you could post both the data that you are doing the fit to, as well as the code that you are running, especially for the "very close, but slightly different" test that you did.

[Pra713]

Yes, I definitely want to add more details and post my codes and results, but I've just been so busy with Grad school this weekend, and probably until this Thursday.

Regarding that 'digest' issue, I had have a few issues before that, which I fixed after Googling. I mainly had to uninstall a couple of previously installed packages. And I did start from a "fresh" R session each time. This 'digest' issue was the one where I got stuck.

[Pra713]

Here are my results for tidem vs Matlab T-Tide. I am using August 15 to September 15 data from NOAA Tides and Currents website.
Link to data:
Raw Data in csv format
Cleaned Data in R's RDS format

################################################

Constituents:

MATLAB (Note that x0 here corresponds to the z0 in tidem)

R

################################################

Plots

Observed VS Predicted (MATLAB)

Observed VS Predicted (R)

Predicted (MATLAB) VS Predicted (R)

[dankelley]

These do not look "very close" (to quote an early characterization). For the R case, you are not reporting the output of summary(tidem(...)), and that is what we need to see, for a comparison.

To repeat the request of @richardsc, can you please state the code you used for the matlab and R cases, and also whether you are using the latest commit of the "develop" branch of oce? Thanks.

[dankelley]

Further on the data,

> load("~/Dropbox/df_data.RDS")
Error in load("~/Dropbox/df_data.RDS") : 
  bad restore file magic number (file may be corrupted) -- no data loaded
In addition: Warning message:
file ‘df_data.RDS’ has magic number 'X'
  Use of save versions prior to 2 is deprecated 

so I cannot read your RDS (perhaps because I have a newer R installed). I can try your .csv file, but you'll need to tell me (by posting your code!) which of the two sea-level columns you are using.

[Pra713]

I'll have to do some cleanup to post my full code and everything else required. This was just a quick and dirty run with no formatting. I can upload everything this Friday as I do not have time right now.

Regarding the data,
The following code will work:
df_data <- readRDS('~/Dropbox/df_data.RDS')
The column I used is "verified_m" if you are using the RDS file or "Verified (m)" if you are using the csv file.

[Pra713]

And apologies for the "very close" quote, now that I look at it more closely, they are not that close. Although the top constituents (constituents with the highest amplitudes) do have comparable amplitudes. But beyond that, they are quite different, which is concerning.

[dankelley]

I wrote R code to read it (using
commit 5d4077f of the "develop" branch of oce), and am showing both the code and its results below.

CODE

library(oce)
e <- read.csv("CO-OPS_8720218_wl.csv")
t <- as.POSIXct(paste(e$Date, e$Time..GMT.), tz="UTC")
e <- e$Verified..m.
sl <- as.sealevel(e, t)
summary(tidem(sl))

OUTPUT

Note: the tidal record is too short to fit for constituents:  SA SSA MSM MF SIG1 RHO1 TAU1 BET1 CHI1 PI1 P1 S1 PSI1 PHI1 THE1 SO1 OQ2 2N2 NU2 GAM2 H1 H2 MKS2 LDA2 T2 R2 K2 MSN2 SO3 MK4 SK4 2MK6 MSK6 
tidem summary
-------------

Call:
tidem(t = sl)
RMS misfit to data:  0.1381428 

Fitted Model:
         Freq Amplitude Phase       p    
Z0   0.000000  0.230167   0.0 < 2e-16 ***
MM   0.001512  0.158694  28.8 < 2e-16 ***
MSF  0.002822  0.067187 358.3 2.8e-11 ***
ALP1 0.034397  0.009362  30.3 0.01158 *  
2Q1  0.035706  0.004546 132.1 0.17815    
Q1   0.037219  0.012508 159.9 0.00464 ** 
O1   0.038731  0.060395 202.5 3.5e-09 ***
NO1  0.040269  0.005807 344.1 0.39586    
K1   0.041781  0.079132 229.0 < 2e-16 ***
J1   0.043293  0.004263 219.1 0.22700    
OO1  0.044831  0.003579 266.1 0.27172    
UPS1 0.046343  0.002460 107.4 0.50203    
EPS2 0.076177  0.024778 355.7 0.00012 ***
MU2  0.077689  0.035582  17.1 < 2e-16 ***
N2   0.078999  0.184835  15.5 < 2e-16 ***
M2   0.080511  0.666149  29.6 < 2e-16 ***
L2   0.082024  0.070815  35.4 < 2e-16 ***
S2   0.083333  0.099766  49.4 < 2e-16 ***
ETA2 0.085074  0.017726 307.5 0.06386 .  
MO3  0.119242  0.011012 311.1 0.00088 ***
M3   0.120767  0.004126 201.6 0.20877    
MK3  0.122292  0.008852 169.7 0.24359    
SK3  0.125114  0.001332 144.4 0.68104    
MN4  0.159511  0.015209 204.1 0.00311 ** 
M4   0.161023  0.036210 199.3 6.3e-08 ***
SN4  0.162333  0.015240 358.9 0.44842    
MS4  0.163845  0.011808 211.4 0.00339 ** 
S4   0.166667  0.005324 290.4 0.21890    
2MK5 0.202804  0.001819 243.1 0.59182    
2SK5 0.208447  0.003879  29.9 0.26156    
2MN6 0.240022  0.004262 196.2 0.33584    
M6   0.241534  0.008063 233.1 0.01863 *  
2MS6 0.244356  0.006855 247.2 0.12272    
2SM6 0.247178  0.002769 273.3 0.58167    
3MK7 0.283315  0.000578  83.6 0.88729    
M8   0.322046  0.002084  52.1 0.51530    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
* Processing Log
    - 2019-11-13 22:16:36.644 UTC: `create 'tidem' object`
    - 2019-11-13 22:16:36.644 UTC: `tidem(t = sl)`

[stale]

This issue has been automatically marked as stale because no comments have been made in a month. Unless comments are made in the next week (or the 'pinned' label is applied), this issue will be closed automatically in a week.

[richardsc]

Pinning because I was in the middle of starting to look into this with Matlab (but which is only licensed to use on a specific network) but haven't finished yet.

[dankelley]

@richardsc still working on this? (I am wondering, because we are on a countdown to a CRAN release, so I'm looking to see if there are any issues that we need to fix beforehand.)

[richardsc]

I'm working on this! Agree would be good to sort out before CRAN submission

[richardsc]

It took me longer to do this evaluation than I planned, and I've learned a lot about tidem() and t_tide along the way.

The short answer, is that I think that there are no fundamental issues with tidem() when compared with t_tide. Both do an essentially equivalent job with prediction when they are called in a similar manner. Some important points:

1. Both tidem() and t_tide use the Rayleigh criterion for determining which constituents to try and fit, and therefore both fit to the same set if not provided with a specified list (the constituents= arg in tidem() and the 'rayleigh' arg in t_tide.

2. However, the defaults for prediction are that t_tide will only use fit constituents with an SNR of 2 or greater. So, the t_tide prediction will differ (slightly) from that obtained by predict(tidem(x)) because of the constituents with low p values. This can be changed by setting 'synthesis', 0 in the t_tide arguments.

3. There are still some numerical differences between the approaches -- for amplitude it is mostly on the order of ~5e-5, but for phase it can be quite large (10's of degrees in some cases). The reason for this is not clear. I plan to open another issue so that @dankelley and I can fiddle with it when we have time to see what is going on.

The notes that were used for this are contained in the pdf below:
compare_t_tide_with_tidem.pdf

The following are a few summary plots showing the results of the comparison (which use the dataset that started this issue):

Difference between tidem() amplitudes and t_tide amplitudes:

t_tide phase (with errorbars) plotted against tidem phase:

tidem prediction minus t_tide prediction:

Note that the rms error between the two predictions and the original sealevel time series are:

• tidem --> 0.1381338
• t_tide --> 0.1381341

So, in terms of how "good" the prediction is, they are essentially the same.

I'm going to close this now, and will open another issue to discuss why there are differences in the fit coefficients. Thanks again @Pra713 for opening this, and I hope this is helpful for you in using tidem() and R instead of t_tide in the future (the truth is t_tide is fantastic -- sometimes I do wish that we had a translation of it in R instead of an alternative approach).