dankelley / oce

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

as.tidem() cannot handle tidal constituent 'M1' #1592

Closed dankelley closed 5 years ago

dankelley commented 5 years ago

I am trying to use as.tidem() on a dataset for Chuuk, Moen Island. I think the data file, communicated to me personally, might be from xtide. At any rate, in its comments at the start, it states the source as http://tidesandcurrents.noaa.gov/ and the file creation date 2009.

This file has a constituent named "M1", but but as.tidem() does not handle M1, so we get a not-very-informative error message. I'll make that more informative in a little while.

But the bigger issue is why oce does not handle M1. As a reminder, oce gets constituent names from Rich Pawlowicz's t-tide matlab package, which I think gets them from Mike Foreman's fortran code. So, did I miss something? Or did Rich and Mike miss something? Or is there another name for M1? I will look into the last of these, first, because I've encountered multiple names before.

At https://en.wikipedia.org/wiki/Theory_of_tides I see M1 named as "Smaller lunar elliptic diurnal", with stated speed of 14.4920521 (and thus freq=speed/360 of 0.0402557). So, to check whether oce has this speed under another name, I did

library(oce)
data(tidedata)
D<-data.frame(name=tidedata$const$name,
    freq=tidedata$const$freq,
    speed=tidedata$const$freq*360)
D$dev <- D$speed - 14.4920521
D[which(abs(D$dev)<0.5),]

which tells me

   name       freq    speed          dev
14 TAU1 0.03895881 14.02517 -0.466879204
15 BET1 0.04004044 14.41456 -0.077495392
16  NO1 0.04026859 14.49669  0.004641884
17 CHI1 0.04047097 14.56955  0.077495444
18  PI1 0.04143851 14.91786  0.425812580
19   P1 0.04155259 14.95893  0.466879256

so nothing actually matches. The closest is NO1, but shouldn't we know these to a lot of digits?

Next steps: look in the t-tide and Foreman docs to see if M1 is mentioned.

dankelley commented 5 years ago

I just checked Foreman's doc (ref 1) and he does not list M1 in the right 'speed' range, as shown by the following screenshot (left: oce constituents near this speed; right: Foreman's document, appendix 1 on page 36, showing name, freq, and comparitor-name)

Screen Shot 2019-08-17 at 8 12 26 am

References

  1. Foreman, M. G. G. “Manual for Tidal Currents Analysis and Prediction.” Pacific Marine Science Report. British Columbia, Canada: Institute of Ocean Sciences, Patricia Bay, 1978.
dankelley commented 5 years ago

Oh, I see in the data file that it is public domain, so I will reproduce it below. I still do not know its provenance.

# Harmonic constants from web snapshot taken 2009-12-22
# Datum from benchmark sheet, publication date 2004-10-21
# BEGIN HOT COMMENTS
# country: F.S.M.
# source: http://tidesandcurrents.noaa.gov/
# restriction: Public domain
# station_id_context: NOS
# station_id: 1840000
# date_imported: 20091223
# datum: Mean Lower Low Water
# confidence: 10
# !units: feet
# !longitude: 151.8467
# !latitude: 7.4467
Chuuk, Moen Island, E. Caroline Islands, F.S.M.
+00:00 :Pacific/Truk
0.8596 feet
J1              0.0350   81.20
K1              0.6170   64.80
K2              0.0990  143.20
L2              0.0060   73.60
M1              0.0280   54.40
M2              0.2370  129.20
M3              0.0140   30.60
x 0 0
x 0 0
x 0 0
N2              0.0540  185.00
2N2             0.0070  240.70
O1              0.3890   44.00
OO1             0.0160   85.60
P1              0.1990   64.20
Q1              0.0750   34.20
2Q1             0.0100   23.20
R2              0.0030  147.40
S1              0.0130  253.00
S2              0.3440  147.50
x 0 0
x 0 0
T2              0.0200  147.50
LDA2            0.0020  137.70
MU2             0.0190  133.40
NU2             0.0100  177.40
RHO1            0.0140   35.10
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
MF              0.0420   16.70
x 0 0
x 0 0
SA              0.1080   62.60
SSA             0.1040   25.70
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
x 0 0
dankelley commented 5 years ago

I have decided not to make the oce functions handle the M1 constituent. I think it's a bad idea to add things that are not in the Foreman (1978) documentation and not in the T_TIDE code.

I am not even convinced that M1 is distinct from the constituents that we have already, evidence for which is the following quote, from the docs for ?as.tidem:

Warnings are issued for any constituent name that is not in this list; as of the late summer of 2019, the only example seen in practice is M1, which according to Wikipedia (2019) has frequency 0.0402557, which is very close to that of NO1, i.e. 0.04026859, perhaps explaining why Foreman (1978) did not handle this constituent. A warning is issued if this or any other unhandled constituent is provided in the name argument to as.tidem(). Who knows, maybe NO1 is another name for M1. We already have several synonyms being handled by oce.

dankelley commented 5 years ago

As an example of how this works, the "develop" branch, commit 93163248d23af59019c4eb09af6cc773c46d2bb4, now ignores M1, making warnings. And the test code has multiple ways to ensure that unknown constituents are being identified and that warnings/errors are issued. So I think it's OK to close this issue, and will do so soon.

dankelley commented 5 years ago

Just to finish off this discussion, the code below makes the graph as shown. If I knew that this file format was "official" (documented, and widely used in open-source software), I might consider coding this into a function to read constituents. One thing that makes me disinclined to do that, however, is that the file format does not seem to state explicitly when the data for the harmonic analysis were measured. Without that, I don't see how to set a proper tRef ... and in the code below, I just fake it, to get a rough idea of what the tides look like.

library(oce)
lines <- readLines("Chuuk-Harmonics-2009.txt")
location <- lines[grep("^#", lines, invert=TRUE)[1]]
lat <- as.numeric(scan(text=lines[grep("latitude:", lines)], what="character", quiet=TRUE)[3])
unit <- scan(text=lines[grep("units:", lines)], what="character", quiet=TRUE)[3]
factor <- if (unit == "feet") .3 else 1
con <- lines[grep("^[A-Z]{1-3}[1-9]{1}", lines)]
T <- read.table(text=con, header=FALSE,
                col.names=c("name", "amp", "phase"),
                stringsAsFactors=FALSE)
if (unit == "feet") {
    T$amp <- 0.3048 * T$amp
    warning("converting amplitude from metres to feet\n")
}
tRef <- as.POSIXct("2019-07-17 00:00:00", tz="UTC")
days <- 28
time <- tRef + seq(0, days*86400, 3600)
tide <- as.tidem(tRef=tRef, latitude=lat, name=T$name, amplitude=T$amp, phase=T$phase)
summary(tide)
eta <- predict(tide, newdata=time)
if (!interactive()) png("chuuk.png")
oce.plot.ts(time, eta, type="l", drawTimeRange=FALSE)
grid()
mtext(location)
if (!interactive()) dev.off()

chuuk

dankelley commented 5 years ago

Here's a ref on M1, which touches on the history: https://www.ocean-sci.net/15/431/2019/os-15-431-2019.pdf