R/tides.R
In dankelley/oce: Analysis of Oceanographic Data
Defines functions
webtide
predict.tidem
tidem
tidemConstituentNameFix
tidemAstron
tidemVuf
as.tidem
Documented in
as.tidem
predict.tidem
tidem
tidemAstron
tidemConstituentNameFix
tidemVuf
webtide
# vim:textwidth=128:expandtab:shiftwidth=4:softtabstop=4
#' Class to Store Tidal Models
#'
#' This class stores tidal-constituent models.
#'
#' @templateVar class tidem
#'
#' @templateVar dataExample {}
#'
#' @templateVar metadataExample {}
#'
#' @template slot_summary
#'
#' @template slot_put
#'
#' @template slot_get
#'
#' @author Dan Kelley
#' @family functions that plot oce data
#' @family things related to tides
setClass("tidem", contains = "oce")
setMethod(
f = "initialize",
signature = "tidem",
definition = function(.Object, ...) {
.Object <- callNextMethod(.Object, ...)
.Object@metadata$version <- ""
.Object@processingLog$time <- presentTime()
.Object@processingLog$value <- "create 'tidem' object"
return(.Object)
}
)
#' Tidal Constituent Information
#'
#' The `tidedata` dataset contains Tide-constituent information that is
#' use by [tidem()] to fit tidal models. `tidedata` is a list
#' containing
#' \describe{
#' \item{`const`}{
#' a list containing vectors
#' `name` (a string with constituent name),
#' `freq` (the frequency, in cycles per hour),
#' `kmpr` (a string naming the comparison constituent, blank if there is none),
#' `ikmpr` (index of comparison constituent, or `0` if there is none),
#' `df` (frequency difference between constituent and its
#' comparison, used in the Rayleigh criterion),
#' `d1` through `d6` (the first through sixth Doodson numbers),
#' `semi`,
#' `nsat` (number of satellite constituents),
#' `ishallow`,
#' `nshallow`,
#' `doodsonamp`,
#' and
#' `doodsonspecies`.
#' }
#' \item{`sat`}{
#' a list containing vectors
#' `deldood`,
#' `phcorr`,
#' `amprat`,
#' `ilatfac`,
#' and
#' `iconst`.
#' }
#' \item{`shallow`}{
#' a list containing vectors
#' `iconst`,
#' `coef`,
#' and
#' `iname`.
#' }
#' }
#' Apart from the use of `d1` through `d6`, the naming and content
#' follows `T_TIDE` (see Pawlowicz et al. 2002), which in turn builds upon
#' the analysis of Foreman (1978).
#'
#' @name tidedata
#'
#' @docType data
#'
#' @author Dan Kelley
#'
#' @references
#'
#' Foreman, M. G. G., 1978. Manual for Tidal Currents Analysis and Prediction.
#' Pacific Marine Science Report. British Columbia, Canada: Institute of Ocean
#' Sciences, Patricia Bay.
#'
#' Pawlowicz, Rich, Bob Beardsley, and Steve Lentz, 2002. Classical tidal
#' harmonic analysis including error estimates in MATLAB using `T_TIDE`.
#' Computers and Geosciences, 28, 929-937.
#'
#' @source The data come from the `tide3.dat` file of the `T_TIDE`
#' package (Pawlowicz et al., 2002), and derive from Appendices provided by
#' Foreman (1978). The data are scanned using \file{tests/tide.R} in this
#' package, which also performs some tests using `T_TIDE` values as a
#' reference.
#'
#' @family things related to tides
NULL
#' Tidal Current Dataset
#'
#' The `tidalCurrent` dataset contains tidal velocities reported in
#' Foreman's (1978) report (reference 1) on his Fortran code for the analysis of
#' tidal currents and provided in an associated webpage (reference 2).
#' Here, `tidalCurrent` is data frame containing
#' * `time` a POSIXct time.
#' * `u` the eastward component of velocity in m/s.
#' * `v` the northward component of velocity in m/s.
#'
#' @name tidalCurrent
#'
#' @docType data
#'
#' @examples
#' library(oce)
#' data(tidalCurrent)
#' par(mfrow = c(2, 1))
#' oce.plot.ts(tidalCurrent$time, tidalCurrent$u, ylab = "u [m/s]")
#' abline(h = 0, col = 2)
#' oce.plot.ts(tidalCurrent$time, tidalCurrent$v, ylab = "v [m/s]")
#' abline(h = 0, col = 2)
#'
#' @author Dan Kelley (reformatting data provided by Michael Foreman)
#'
#' @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.
#' 2. \code{https://www.dfo-mpo.gc.ca/science/documents/data-donnees/tidal-marees/tidpack.zip}
#'
#' @source The data come from the `tide8.dat` and `tide9.dat` files provided
#' at reference 2.
#'
#' @family things related to tides
NULL
#' Summarize a tidem Object
#'
#' By default, all fitted constituents are plotted, but it is quite useful to
#' set e.g. p=0.05 To see just those constituents that are significant at the 5
#' percent level.
#' Note that the p values are estimated as the average of the p values for the
#' sine and cosine components at a given frequency.
#'
#' @param object an object of class [tidem], as created by
#' [as.tidem()] or [tidem()].
#'
#' @param p optional value of the maximum p value for the display of an
#' individual coefficient. If not given, all coefficients are shown.
#'
#' @param constituent optional character vector holding the names
#' of constituents on which to focus.
#' @template tideconst
#'
#' @param \dots further arguments passed to or from other methods.
#'
#' @return `NULL`
#'
#' @author Dan Kelley
#'
#' @section Sample of Usage:
#' \preformatted{
#' library(oce)
#' data(sealevel)
#' tide <- tidem(sealevel)
#' summary(tide)
#' }
#'
#' @family things related to tides
setMethod(
f = "summary",
signature = "tidem",
definition = function(object, p = 1.0, constituent, ...) {
debug <- if ("debug" %in% names(list(...))) list(...)$debug else 0
version <- object@metadata$version
ok <- object@data$p <= p | version == 3L
haveP <- any(!is.na(object@data$p))
if (missing(constituent)) {
fit <- data.frame(
Const = object@data$const[ok],
Name = object@data$name[ok],
Freq = object@data$freq[ok],
Amplitude = object@data$amplitude[ok],
Phase = object@data$phase[ok],
p = object@data$p[ok]
)
if (debug) {
cat("For missing(constituent) case, fit is:\n")
print(fit)
}
} else {
i <- NULL
bad <- NULL
for (iconst in seq_along(constituent)) {
w <- which(object@data$name == constituent[iconst])
if (length(w) == 1L) {
i <- c(i, w)
} else {
bad <- c(bad, iconst)
}
}
if (length(bad)) {
warning("the following constituents are not handled: \"",
paste(constituent[bad], collapse = "\", \""), "\"\n",
sep = ""
)
}
if (length(i) == 0) {
stop("In summary,tidem-method() : no known constituents were provided", call. = FALSE)
}
i <- unique(i)
fit <- data.frame(
Const = object@data$const[i],
Name = object@data$name[i],
Freq = object@data$freq[i],
Amplitude = object@data$amplitude[i],
Phase = object@data$phase[i],
p = object@data$p[i]
)
if (debug) {
cat("For !missing(constituent) case, fit is:\n")
print(fit)
}
}
cat("tidem summary\n-------------\n")
if (version != "3") {
cat("\nCall:\n")
cat(paste(deparse(object[["call"]]), sep = "\n", collapse = "\n"), "\n", sep = "")
cat("RMS misfit to data: ", sqrt(var(object[["model"]]$residuals)), "\n")
cat("\nFitted Model:\n")
f <- fit[3:6]
if (debug > 0) {
cat("fit:\n")
print(fit)
cat("f:\n")
print(f)
}
rownames(f) <- as.character(fit[, 2])
if (haveP) {
printCoefmat(f,
digits = 3,
signif.stars = getOption("show.signif.stars"),
signif.legend = TRUE,
P.values = TRUE, has.Pvalue = TRUE, ...
)
} else {
printCoefmat(f[, -4], digits = 3)
}
} else {
cat("\nSupplied Model:\n")
f <- fit[3:5]
rownames(f) <- as.character(fit[, 2])
printCoefmat(f, digits = 3)
}
processingLogShow(object)
invisible(NULL)
}
)
#' Extract Something From a tidem Object
#'
#' @param x a [tidem-class] object.
#'
#' @section Details of the Specialized Method:
#'
#' * If `i` is `"?"`, then the return value is a list
#' containing four items, each of which is a character vector
#' holding the names of things that can be accessed with `[[`.
#' The `data` and `metadata` items hold the names of
#' entries in the object's data and metadata
#' slots, respectively. Note that `metadataDerived`
#' holds only `""`, because no derived metadata values
#' are defined for `tidem` objects.
#'
#' * If `i` is `"frequency"` or `"freq"`, then a vector of
#' constituent frequencies is returned.
#'
#' * If `i` is `"amplitude"` then a vector of constituent amplitudes
#' is returned.
#'
#' * If `i` is `"phase"` then a vector of constituent phases
#' is returned.
#'
#' * If `i` is `"constituents"` then a data frame holding constituent
#' name, frequency, amplitude and phase is returned.
#'
#' * If `i` is a vector of constituent names, then the return
#' value is as for `"constituents"`, except that only the named
#' those constituents are returned.
#'
#' @template sub_subTemplate
#'
#' @family things related to tides
setMethod(
f = "[[",
signature(x = "tidem", i = "ANY", j = "ANY"),
definition = function(x, i, j, ...) {
# UNUSED const <- x@data$name
if (i[1] == "?") {
return(list(
metadata = sort(names(x@metadata)),
metadataDerived = NULL,
data = c(sort(names(x@data)), x@data$name),
dataDerived = "frequency"
))
}
if (i[1] == "frequency") {
return(x@data$freq)
} else if (i[1] == "constituents" || sum(i %in% x@data$name) > 0L) {
all <- with(
x@data,
data.frame(
name = name, freq = freq,
amplitude = amplitude, phase = phase
)
)
if (i[[1]] == "constituents") {
return(all)
} else {
return(all[all$name %in% i, ])
}
}
callNextMethod() # [[
}
)
#' Replace Parts of a tidem Object
#'
#' @param x a [tidem-class] object.
#'
#' @template sub_subsetTemplate
#'
#' @family things related to tides
setMethod(
f = "[[<-",
signature(x = "tidem", i = "ANY", j = "ANY"),
definition = function(x, i, j, ..., value) {
callNextMethod(x = x, i = i, j = j, ... = ..., value = value) # [[<-
}
)
#' Plot a tidem Object
#'
#' Plot a summary diagram for a tidal fit.
#'
#' @param x a [tidem-class] object.
#'
#' @param which integer flag indicating plot type, 1 for stair-case spectral, 2
#' for spike spectral.
#'
#' @param constituents character vector holding the names of constituents that are
#' to be drawn and labelled. If `NULL`, then no constituents will be shown.
#' @template tideconst
#'
#' @param sides an integer vector of length equal to that of `constituents`,
#' designating the side on which the constituent labels are to be drawn. As in
#' all R graphics, the value `1` indicates the bottom of the plot, and
#' `3` indicates the top. If `sides=NULL`, the default, then all labels
#' are drawn at the top. Any value of `sides` that is not either 1 or 3
#' is converted to 3.
#'
#' @param col a character vector naming colors to be used for `constituents`.
#' Ignored if `sides=3`. Repeated to be of the same length as
#' `constituents`, otherwise.
#'
#' @param log if set to "`x`", the frequency axis will be logarithmic.
#'
#' @param mgp 3-element numerical vector to use for `par(mgp)`, and also
#' for `par(mar)`, computed from this. The default is tighter than the R
#' default, in order to use more space for the data and less for the axes.
#'
#' @param mar value to be used with `[`par`]("mar")`.
#'
#' @param \dots optional arguments passed to plotting functions, not all
#' of which are obeyed. For example, if \dots contains `type`, that value will be
#' ignored because it is set internally, according to the value of `which`.
#'
#' @section Sample of Usage:
#' \preformatted{
#' library(oce)
#' data(sealevel)
#' tide <- tidem(sealevel)
#' plot(tide)
#' }
#'
#' @section Historical note:
#' An argument named `labelIf` was removed in July 2016,
#' because it was discovered never to have worked as documented, and
#' because the more useful argument `constituents` had been added.
#'
#' @author Dan Kelley
#'
#' @family functions that plot oce data
#'
#' @aliases plot.tidem
#'
#' @family things related to tides
setMethod(
f = "plot",
signature = signature("tidem"),
definition = function(x,
which = 1,
constituents = c("SA", "O1", "K1", "M2", "S2", "M4"),
sides = NULL,
col = "blue",
log = "",
mgp = getOption("oceMgp"),
mar = c(mgp[1] + 1, mgp[1] + 1, mgp[2] + 0.25, mgp[2] + 1),
...) {
dots <- list(...)
dotsNames <- names(dots)
data("tidedata", package = "oce", envir = environment())
tidedata <- get("tidedata") # , pos=globalenv())
drawConstituent <- function(name = "M2", side = 3, col = "blue", adj = NULL, cex = 0.8) {
w <- which(tidedata$const$name == name)
if (!length(w)) {
warning("constituent \"", name, "\" is unknown")
return()
}
frequency <- tidedata$const$freq[w]
abline(v = frequency, col = col, lty = "dotted")
if (par("usr")[1] < frequency && frequency <= par("usr")[2]) {
if (is.null(adj)) {
mtext(name, side = side, at = frequency, col = col, cex = cex)
} else {
mtext(name, side = side, at = frequency, col = col, cex = cex, adj = adj)
}
}
}
opar <- par(no.readonly = TRUE)
lw <- length(which)
if (lw > 1) {
on.exit(par(opar))
}
par(mgp = mgp, mar = mar)
frequency <- x@data$freq[-1] # trim z0
amplitude <- x@data$amplitude[-1]
# name <- x@data$name[-1]
# Place in frequency order (required for cumulative plot), because user
# might have given constituents in any order.
order <- order(frequency)
frequency <- frequency[order]
# name <- name[order]
amplitude <- amplitude[order]
# print(data.frame(name=name, period=1/frequency, amplitude=amplitude))
if (!is.null(constituents)) {
sides <- if (is.null(sides)) {
rep(3, length(constituents))
} else {
rep(sides, length.out = length(constituents))
}
col <- rep(col, length.out = length(constituents))
}
sides[sides != 1 & sides != 3] <- 3 # default to top
cex <- if ("cex" %in% names(dots)) dots$cex else 0.8
# We specify these things directly but they are not parameters of this
# function, so we must remove them from '...' before passing that to plot().
if ("cex" %in% dotsNames) dots$cex <- NULL
if ("type" %in% dotsNames) dots$type <- NULL
for (w in 1:lw) {
# message("w=", w, "; which[w]=", which[w])
if (which[w] == 2) {
plot(frequency, amplitude,
type = "n",
xlab = "Frequency [ cph ]", ylab = "Amplitude [ m ]", log = log, ...
)
segments(frequency, 0, frequency, amplitude)
for (i in seq_along(constituents)) {
drawConstituent(constituents[i], side = sides[i], col = col[i], cex = cex)
}
} else if (which[w] == 1) {
plot(frequency, cumsum(amplitude),
xlab = "Frequency [ cph ]", ylab = "Amplitude [ m ]", log = log, type = "s", ...
)
for (i in seq_along(constituents)) {
drawConstituent(constituents[i], side = sides[i], col = col[i], cex = cex)
}
} else {
stop("unknown value of which ", which, "; should be 1 or 2")
}
}
if (!all(is.na(pmatch(names(list(...)), "main")))) title(...)
}
)
#' Create tidem Object From Fitted Harmonic Data
#'
#' This function takes a set of tidal constituent amplitudes
#' and phases, and constructs a return value of similar form
#' to that returned by [tidem()]. Its purpose is to enable
#' predictions based on published constituent amplitudes
#' and phases. Since [as.tidem()] does not account for a
#' reference height, it is the user's responsible to account
#' for this after a prediction is made using [predict.tidem()].
#'
#' All the constituent names used by [tidem()] are permitted here,
#' *except* for `"Z0"` (see \dQuote{Description} regarding reference
#' height).
#' To get a list of constituent names, please consult Foreman (1978),
#' or type the following in an R console:
#' \preformatted{
#' data(tidedata)
#' data.frame(name=tidedata$const$name, freq=tidedata$const$freq)
#' }
#'
#' In addition to the above, [as.tidem()] can handle NOAA names
#' of constituents. For the most part, these match oce names, but
#' there are 4 exceptions: NOAA names
#' `"LAM2"`, `"M1"`, `"RHO"`, and `"2MK3"` are converted to oce names
#' `"LDA2"`, `"NO1"`, `"RHO1"`, and `"MO3"`. The name mapping was
#' inferred by matching frequencies; for these constituents, the
#' fractional mismatch in frequencies was under 4e-8;
#' see Reference 5 for more details.
#' A message is printed if these name conversions are required
#' in the particular use of [as.tidem()].
#'
#' Apart from the standard oce names and this set of NOAA synonyms,
#' any other constituent name is reported in a warning message.
#'
#' @param tRef a POSIXt value indicating the mean time of the
#' observations used to develop the harmonic model. This is rounded
#' to the nearest hour in [as.tidem()], to match the behaviour
#' of [tidem()].
#'
#' @param latitude numerical value indicating the latitude of the
#' observations that were used to create the harmonic model. This
#' is needed for nodal-correction procedures carried out
#' by [tidemVuf()].
#'
#' @param name character vector holding names of constituents.
#'
#' @param amplitude,phase numeric vectors of constituent amplitudes
#' and phases. These must be of the same length as `name`.
#'
#' @param frequency,speed optional numeric vectors giving the frequencies
#' of the constituents (in cycles per hour) or the analogous speeds
#' (in degrees per hour). Only one of these may be given, and a conversion
#' is done from the latter to the former, if required. If the frequencies
#' are thus specified, then these are used instead of the frequencies that
#' oce normally used, as defined in `data(tideconst)`. A warning will
#' be issued if the absolute value of the relative frequency mismatch for any
#' given component exceeds 1e-6, and this will occur for any NOAA tables
#' containing the SA component, for which this relative mismatch
#' is approximately 4e-5 (see reference 5).
#'
#' @template debugTemplate
#'
#' @return An object of [tidem-class], with only minimal
#' contents.
#'
#' @section Known issues:
#' There are two known differences between [tidem()] and the Matlab
#' `T_TIDE` package, as listed in references 3 and 4.
#'
#' @examples
#'
#' # Example 1: show agreement with tidem()
#' data(sealevelTuktoyaktuk)
#' # 'm0' is model fitted by tidem()
#' m0 <- tidem(sealevelTuktoyaktuk)
#' p0 <- predict(m0, sealevelTuktoyaktuk[["time"]])
#' m1 <- as.tidem(
#' mean(sealevelTuktoyaktuk[["time"]]), sealevelTuktoyaktuk[["latitude"]],
#' m0[["name"]], m0[["amplitude"]], m0[["phase"]]
#' )
#' # Test agreement with tidem() result, by comparing predicted sealevels.
#' p1 <- predict(m1, sealevelTuktoyaktuk[["time"]])
#' stopifnot(max(abs(p1 - p0), na.rm = TRUE) < 1e-10)
#'
#' # Example 2: See the effect of dropping weak constituents
#' m0[["name"]][which(m0[["amplitude"]] > 0.05)]
#' h <- "
#' name amplitude phase
#' Z0 1.98061875 0.000000
#' MM 0.21213065 263.344739
#' MSF 0.15605629 133.795004
#' O1 0.07641438 74.233130
#' K1 0.13473817 81.093134
#' OO1 0.05309911 235.749693
#' N2 0.08377108 44.521462
#' M2 0.49041340 77.703594
#' S2 0.22023705 137.475767"
#' coef <- read.table(text = h, header = TRUE)
#' m2 <- as.tidem(
#' mean(sealevelTuktoyaktuk[["time"]]),
#' sealevelTuktoyaktuk[["latitude"]],
#' coef$name, coef$amplitude, coef$phase
#' )
#' p2 <- predict(m2, sealevelTuktoyaktuk[["time"]])
#' par(mfrow = c(3, 1))
#' oce.plot.ts(sealevelTuktoyaktuk[["time"]], p0)
#' ylim <- par("usr")[3:4] # to match scales in other panels
#' oce.plot.ts(sealevelTuktoyaktuk[["time"]], p1, ylim = ylim)
#' oce.plot.ts(sealevelTuktoyaktuk[["time"]], p2, ylim = ylim)
#'
#' @references
#'
#' 1. Foreman, M. G. G., 1978. Manual for Tidal Currents Analysis and Prediction.
#' Pacific Marine Science Report. British Columbia, Canada: Institute of Ocean
#' Sciences, Patricia Bay.
#'
#' 2. Wikipedia, "Theory of Tides." https://en.wikipedia.org/wiki/Theory_of_tides
#' Downloaded Aug 17, 2019.
#'
#' 3. Github issue 1653 "tidem() and t_tide do not produce identical results"
#' (https://github.com/dankelley/oce/issues/1653)
#'
#' 4. Github issue 1654 "predict(tidem()) uses all constituents, unlike T_TIDE"
#' (https://github.com/dankelley/oce/issues/1654)
#'
#' 5. Github issue 2143 "mismatch in oce/NOAA frequency of SA tidal constituent"
#' (https://github.com/dankelley/oce/issues/2143)
#'
#' @family things related to tides
as.tidem <- function(tRef, latitude, name, amplitude, phase, frequency, speed, debug = getOption("oceDebug")) {
oceDebug(debug, "as.tidem() {\n", sep = "", unindent = 1)
if (missing(tRef)) {
stop("tRef must be given")
}
if (missing(latitude)) {
stop("latitude must be given")
}
if (missing(name)) {
stop("name must be given")
}
if (missing(amplitude)) {
stop("amplitude must be given")
}
if (missing(phase)) {
stop("phase must be given")
}
gaveFrequency <- !missing(frequency)
gaveSpeed <- !missing(speed)
if (gaveSpeed) {
if (gaveFrequency) {
stop("cannot give both speed and frequency")
}
frequency <- speed / 360.0
gaveFrequency <- TRUE
}
nname <- length(name)
if (nname != length(amplitude)) {
stop("lengths of name and amplitude must be equal but they are ", nname, " and ", length(amplitude))
}
if (nname != length(phase)) {
stop("lengths of name and phase must be equal but they are ", nname, " and ", length(phase))
}
# convert from NOAA name to oce name
# NOAA: LAM2 M1 RHO 2MK3
# oce: LDA2 NO1 RHO1 MO3
tmpNOAA <- c("LAM2", "M1", "RHO", "2MK3")
tmpOce <- c("LDA2", "NO1", "RHO1", "MO3")
for (tmpi in seq_along(tmpNOAA)) {
if (tmpNOAA[tmpi] %in% name) {
name[name == tmpNOAA[tmpi]] <- tmpOce[tmpi]
message("converted NOAA name \"", tmpNOAA[tmpi], "\" to oce name \"", tmpOce[tmpi], "\"")
}
}
data("tidedata", package = "oce", envir = environment())
tidedata <- get("tidedata")
tRef <- numberAsPOSIXct(3600 * round(as.numeric(tRef, tz = "UTC") / 3600), tz = "UTC")
oceDebug(debug, "input head(name): ", paste(head(name), collapse = " "), "\n")
oceDebug(debug, "input head(phase): ", paste(head(phase), collapse = " "), "\n")
oceDebug(debug, "input head(amplitude): ", paste(head(amplitude), collapse = " "), "\n")
freq <- rep(NA, nname)
indices <- rep(NA, nname)
ibad <- NULL
if (gaveSpeed) {
oceDebug(debug, "computing frequency from velocity\n")
frequency <- speed / 360.0
gaveFrequency <- TRUE
}
for (i in seq_along(name)) {
oceDebug(debug, "adjusting amplitude and phase for constituent \"", name[i], "\"\n", sep = "")
j <- which(tidedata$const$name == name[i])
oceDebug(debug, " inferred j=", j, " from constituent name\n", sep = "")
if (length(j)) {
vuf <- tidemVuf(tRef, j = j, latitude = latitude)
oceDebug(debug, " inferred vuf=", deparse(vuf), "\n")
indices[i] <- j
amplitude[i] <- amplitude[i] * vuf$f
phase[i] <- phase[i] - (vuf$v + vuf$u) * 360
if (gaveFrequency) {
freq[i] <- frequency[i]
afd <- abs(frequency[i] - tidedata$const$freq[j]) / tidedata$const$freq[j]
if (afd > 1e-6) {
msg <- sprintf("%s: absolute fractional mismatch between specified and built-in frequency is %.3e\n", name[i], afd)
warning(msg)
}
oceDebug(debug, msg)
} else {
freq[i] <- tidedata$const$freq[j]
}
} else {
ibad <- c(ibad, i)
}
}
if (length(ibad)) {
warning("the following constituents are not handled: \"",
paste(name[ibad], collapse = "\", \""), "\"\n",
sep = ""
)
indices <- indices[-ibad]
name <- name[-ibad]
amplitude <- amplitude[-ibad]
phase <- phase[-ibad]
freq <- freq[-ibad]
}
oceDebug(debug, "after vuf correction, head(name): ", paste(head(name), collapse = " "), "\n")
oceDebug(debug, "after vuf correction, head(phase): ", paste(head(phase), collapse = " "), "\n")
oceDebug(debug, "after vuf correction, head(amplitude): ", paste(head(amplitude), collapse = " "), "\n")
oceDebug(debug, "} # as.tidem()\n", sep = "", unindent = 1)
phase <- phase %% 360
res <- new("tidem")
res@data <- list(
tRef = tRef,
const = indices,
name = name,
freq = freq,
amplitude = amplitude,
phase = phase,
p = rep(NA, length(name))
)
res@metadata$version <- 3
res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep = "", collapse = ""))
oceDebug(debug, "} # as.tidem()\n", sep = "", unindent = 1)
res
}
#' Nodal Modulation Calculations for Tidal Analyses
#'
#' Carry out nodal modulation calculations for [tidem()]. This function is based directly
#' on `t_vuf` in the `T_TIDE` Matlab package (Pawlowicz et al., 2002),
#' which inherits from the Fortran code described by Foreman (1978).
#'
#' @param t a single time in [POSIXct()] format, with timezone `"UTC"`.
#'
#' @param j integer vector, giving indices of tidal constituents to use.
#'
#' @param latitude optional numerical value containing the latitude in degrees North.
#' If not provided, `u` in the return value will be a vector consisting of
#' repeated 0 value, and `f` will be a vector of repeated 1 value.
#'
#' @return A `list` containing
#' items named `v`, `u` and `f` as described in the `T_TIDE` documentation,
#' as well in Pawlowicz et al. (2002) and Foreman (1978).
#'
#' @author Dan Kelley translated this from the `t_vuf` function
#' of the `T_TIDE` Matlab package (see Pawlowicz et al. 2002).
#'
#' @examples
#' # Look up values for the M2 constituent in Halifax Harbour, Canada.
#' library(oce)
#' data("tidedata")
#' j <- with(tidedata$const, which(name == "M2"))
#' tidemVuf(t = as.POSIXct("2008-01-22 18:50:24"), j = j, lat = 44.63)
#'
#' @references
#' * Foreman, M. G. G., 1978. Manual for Tidal Currents Analysis and Prediction.
#' Pacific Marine Science Report. British Columbia, Canada: Institute of Ocean
#' Sciences, Patricia Bay.
#'
#' * Pawlowicz, Rich, Bob Beardsley, and Steve Lentz, 2002. Classical tidal
#' harmonic analysis including error estimates in MATLAB using `T_TIDE`.
#' Computers and Geosciences, 28, 929-937.
#'
#' @family things related to tides
tidemVuf <- function(t, j, latitude = NULL) {
debug <- 0
if (length(t) > 1) {
stop("t must be a single POSIXct item")
}
data("tidedata", package = "oce", envir = environment())
tidedata <- get("tidedata") # , pos=globalenv())
a <- tidemAstron(t)
if (debug > 0) print(a)
doodson <- cbind(
tidedata$const$d1,
tidedata$const$d2,
tidedata$const$d3,
tidedata$const$d4,
tidedata$const$d5,
tidedata$const$d6
)
oceDebug(
debug,
"doodson[1,]=", doodson[1, ], "\n",
"doodson[2,]=", doodson[2, ], "\n",
"doodson[3,]=", doodson[3, ], "\n"
)
v <- doodson %*% a$astro + tidedata$const$semi
oceDebug(debug, "tidedata$const$semi[", j, "]=", tidedata$const$semi[j], "\n")
v <- v - trunc(v)
oceDebug(debug, "v[1:3]=", v[1:3], "\n")
if (!is.null(latitude) && !is.na(latitude)) {
if (abs(latitude) < 5) {
latitude <- sign(latitude) * 5
}
slat <- sin(pi * latitude / 180)
k <- which(tidedata$sat$ilatfac == 1)
rr <- tidedata$sat$amprat
rr[k] <- rr[k] * 0.36309 * (1.0 - 5.0 * slat * slat) / slat
k <- which(tidedata$sat$ilatfac == 2)
rr[k] <- rr[k] * 2.59808 * slat
uu <- tidedata$sat$deldood %*% a$astro[4:6] + tidedata$sat$phcorr
uu <- uu - trunc(uu)
oceDebug(debug, "uu[1:3]=", uu[1:3], "\n")
nsat <- length(tidedata$sat$iconst)
# nfreq <- length(tidedata$const$numsat)
# loop, rather than make a big matrix
oceDebug(
debug,
"tidedata$sat$iconst=", tidedata$sat$iconst, "\n",
"length(sat$iconst)=", length(tidedata$sat$iconst), "\n"
)
fsum.vec <- vector("numeric", nsat)
u.vec <- vector("numeric", nsat)
for (isat in 1:nsat) {
oceDebug(debug, "isat=", isat, "\n")
use <- tidedata$sat$iconst == isat
fsum.vec[isat] <- 1 + sum(rr[use] * exp(1i * 2 * pi * uu[use]))
u.vec[isat] <- Arg(fsum.vec[isat]) / 2 / pi
if (isat == 8 && debug > 0) {
cat("TEST at isat=8:\n")
cat("fsum.vec[", isat, "]=", fsum.vec[isat], " (EXPECT 1.18531604917590 - 0.08028013402313i)\n")
cat("u.vec[ ", isat, "]=", u.vec[isat], " (EXPECT -0.01076294959868)\n")
}
}
oceDebug(
debug,
"uvec[", j, "]=", u.vec[j], "\n",
"fsum.vec[", j, "]=", fsum.vec[j], "\n"
)
f <- abs(fsum.vec)
u <- Arg(fsum.vec) / 2 / pi
oceDebug(debug, "f=", f, "\n")
oceDebug(debug, "u=", u, "\n")
for (k in which(!is.na(tidedata$const$ishallow))) {
ik <- tidedata$const$ishallow[k] + 0:(tidedata$const$nshallow[k] - 1)
f[k] <- prod(f[tidedata$shallow$iname[ik]]^abs(tidedata$shallow$coef[ik]))
u[k] <- sum(u[tidedata$shallow$iname[ik]] * tidedata$shallow$coef[ik])
v[k] <- sum(v[tidedata$shallow$iname[ik]] * tidedata$shallow$coef[ik])
if (debug > 0 && k < 28) {
cat("k=", k, "f[k]=", f[k], " u[k]=", u[k], "v[k]=", v[k], "\n")
}
}
u <- u[j]
v <- v[j]
f <- f[j]
} else {
v <- v[j]
u <- rep(0, length(j))
f <- rep(1, length(j))
}
if (length(v) < length(u)) {
warning("trimming u and f to get same length as v -- this is a bug")
u <- u[seq_along(v)]
f <- f[seq_along(v)]
}
list(v = v, u = u, f = f)
}
#' Astronomical Calculations for tidem
#'
#' Do some astronomical calculations for [tidem()]. This function is based directly
#' on `t_astron` in the `T_TIDE` Matlab package (see Pawlowicz et al. 2002),
#' which inherits from the Fortran code described by Foreman (1978).
#'
#' @param t Either a time in `POSIXct` format (with `"UTC"` timezone,
#' or else odd behaviours may result),
#' or an integer. In the second case, it is converted to a time with
#' [numberAsPOSIXct()], using `tz="UTC"`.
#'
#' @return A `list` containing items named
#' `astro` and `ader` (see the `T_TIDE` documentation).
#'
#' @author Dan Kelley translated this from the `t_astron` function
#' of the `T_TIDE` Matlab package (see Pawlowicz et al. 2002).
#'
#' @examples
#' tidemAstron(as.POSIXct("2008-01-22 18:50:24"))
#'
#' @references
#' * Foreman, M. G. G., 1978. Manual for Tidal Currents Analysis and Prediction.
#' Pacific Marine Science Report. British Columbia, Canada: Institute of Ocean
#' Sciences, Patricia Bay.
#'
#' * Pawlowicz, Rich, Bob Beardsley, and Steve Lentz, 2002. Classical tidal
#' harmonic analysis including error estimates in MATLAB using `T_TIDE`.
#' Computers and Geosciences, 28, 929-937.
#'
#' @family things related to tides
tidemAstron <- function(t) {
if (length(t) > 1) {
stop("t must be a single POSIXct item")
}
debug <- FALSE
if (is.numeric(t)) {
t <- numberAsPOSIXct(t, tz = "UTC")
}
d <- as.numeric(difftime(t, ISOdatetime(1899, 12, 31, 12, 0, 0, tz = "UTC"), units = "days"))
D <- d / 10000
a <- matrix(c(1.0, d, D^2, D^3), ncol = 1)
oceDebug(debug, "d=", formatC(d, digits = 10), "D=", D, "a=", a, "\n")
scHcPcNpPp <-
matrix(
c(
270.434164, 13.1763965268, -0.0000850, 0.000000039,
279.696678, 0.9856473354, 0.00002267, 0.000000000,
334.329556, 0.1114040803, -0.0007739, -0.00000026,
-259.183275, 0.0529539222, -0.0001557, -0.000000050,
281.220844, 0.0000470684, 0.0000339, 0.000000070
),
nrow = 5, ncol = 4, byrow = TRUE
)
astro <- ((scHcPcNpPp %*% a) / 360) %% 1
oceDebug(debug, "astro=", astro, "\n")
rem <- as.numeric(difftime(t, trunc.POSIXt(t, units = "days"), tz = "UTC", units = "days"))
oceDebug(debug, "rem2=", rem, "\n")
tau <- rem + astro[2, 1] - astro[1, 1]
astro <- c(tau, astro)
da <- matrix(c(0.0, 1.0, 2e-4 * D, 3e-4 * D^2), nrow = 4, ncol = 1)
ader <- (scHcPcNpPp %*% da) / 360
dtau <- 1 + ader[2, 1] - ader[1, 1]
ader <- c(dtau, ader)
list(astro = astro, ader = ader)
}
#' Change Tidal Constituent Name from T-TIDE to Foreman Convention
#'
#' This is used by [tidem()] to permit users to specify constituent names in either
#' the T-TIDE convention (see Pawlowicz et al. 2002) or Foreman convention
#' (see Foreman (1978). There are only two such instances:
#' `"MS"`, which gets translated to `"M8"`, and `"UPSI"`,
#' which gets translated to `"UPS1"`.
#'
#' @param names a vector of character values, holding constituent names
#'
#' @param debug an integer controlling warnings. If this is zero, then no warnings
#' are issued during processing; otherwise, as is the default, warnings are
#' issued for each conversion that is required.
#'
#' @return A vector of character values of tidal constituent names, in the Foreman naming convention.
#'
#' @references
#' Foreman, M. G. G., 1978. Manual for Tidal Currents Analysis and Prediction.
#' Pacific Marine Science Report. British Columbia, Canada: Institute of Ocean
#' Sciences, Patricia Bay.
#'
#' Pawlowicz, Rich, Bob Beardsley, and Steve Lentz, 2002. Classical tidal
#' harmonic analysis including error estimates in MATLAB using `T_TIDE`.
#' Computers and Geosciences, 28, 929-937.
tidemConstituentNameFix <- function(names, debug = 1) {
if ("MS" %in% names) {
if (debug) {
warning("constituent name switched from T-TIDE 'MS' to Foreman (1978) 'M8'")
}
names[names == "MS"] <- "M8"
}
if ("-MS" %in% names) {
if (debug) {
warning("removed-constituent name switched from T-TIDE 'MS' to Foreman (1978) 'M8'")
}
names[names == "-MS"] <- "-M8"
}
if ("UPSI" %in% names) {
if (debug) {
warning("constituent name switched from T-TIDE 'UPSI' to Foreman (1978) 'UPS1'")
}
names[names == "UPSI"] <- "UPS1"
}
if ("-UPSI" %in% names) {
if (debug) {
warning("removed-constituent name switched from T-TIDE 'UPSI' to Foreman (1978) 'UPS1'")
}
names[names == "-UPSI"] <- "-UPS1"
}
names
}
#' Fit a Tidal Model to a Timeseries
#'
#' The fit is done in terms of sine and cosine components at the indicated
#' tidal frequencies (after possibly pruning to satisfy the Rayleigh criterion,
#' as explained in phase 2 of the procedure outlined in
#' \dQuote{Details}), with the amplitude and phase being
#' calculated from the resultant coefficients on the sine and cosine terms.
#' The scheme was devised for hourly data; for other sampling schemes,
#' see \dQuote{Application to non-hourly data}.
#'
#' A summary of constituents used by [tidem()] may be found with:
#' \preformatted{
#' data(tidedata)
#' print(tidedata$const)
#' }
#'
#' A multi-stage procedure is used to establish the list of tidal
#' constituents to be used in the fit.
#'
#' *Phase 1: initial selection.*
#'
#' The initial list tidal constituents (prior to the pruning phase described
#' below) to be used in the analysis are specified as follows;
#' see \dQuote{Constituent Naming Convention}.
#'
#' 1. If `constituents` is not provided, then the constituent
#' list will be made up of the 69 constituents designated by Foreman as "standard".
#' These include astronomical frequencies and some shallow-water frequencies,
#' and are as follows: `c("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", "SO1", "OO1", "UPS1",
#' "OQ2", "EPS2", "2N2", "MU2", "N2", "NU2", "GAM2", "H1", "M2", "H2", "MKS2",
#' "LDA2", "L2", "T2", "S2", "R2", "K2", "MSN2", "ETA2", "MO3", "M3", "SO3",
#' "MK3", "SK3", "MN4", "M4", "SN4", "MS4", "MK4", "S4", "SK4", "2MK5", "2SK5",
#' "2MN6", "M6", "2MS6", "2MK6", "2SM6", "MSK6", "3MK7", "M8")`.
#'
#' 2. If the first item in `constituents` is the string
#' `"standard"`, then a provisional list is set up as in Case 1, and then
#' the (optional) rest of the elements of `constituents` are examined, in
#' order. Each of these constituents is based on the name of a tidal
#' constituent in the Foreman (1978) notation. (To get the list, execute
#' `data(tidedata)` and then execute `cat(tideData$name)`.) Each
#' named constituent is added to the existing list, if it is not already there.
#' But, if the constituent is preceded by a minus sign, then it is removed
#' from the list (if it is already there). Thus, for example, a user
#' might specify `constituents=c("standard", "-M2", "ST32")` to remove the M2
#' constituent and add the ST32 constituent.
#'
#' 3. If the first item is not `"standard"`, then the list of
#' constituents is processed as in Case 2, but without starting with the
#' standard list. As an example, `constituents=c("K1", "M2")` would fit
#' for just the K1 and M2 components. (It would be strange to use a minus sign
#' to remove items from the list, but the function allows that.)
#'
#' In each of the above cases, the list is reordered in frequency prior to the
#' analysis, so that the results of [summary,tidem-method()] will be in a
#' familiar form.
#'
#' *Phase 2: pruning based on Rayleigh's criterion.*
#'
#' Once the initial constituent list is determined during Phase 1,
#' `tidem` applies the Rayleigh criterion, which holds that
#' constituents of frequencies \eqn{f_1}{f1} and \eqn{f_2}{f2} cannot be
#' resolved unless the time series spans a time interval of at least
#' \eqn{rc/(f_1-f_2)}{rc/(f1-f2)}. The details of the decision of which
#' constituent to prune is fairly complicated, involving decisions
#' based on a comparison table. The details of this process are provided
#' by Foreman (1978).
#'
#' *Phase 3: optionally overruling Rayleigh's criterion.*
#'
#' The pruning list from phase 2 can be overruled by the user. This
#' is done by retaining any constituents that the user has named
#' in the `constituents` parameter. This action
#' was added on 2017-12-27, to make `tidem` behave the same
#' way as the Foreman (1978) code, as illustrated in his
#' Appendices 7.2 and 7.3. (As an aside, his Appendix 7.3 has some errors,
#' e.g. the frequency for the 2SK5 constituent is listed there (p58) as
#' 0.20844743, but it is listed as 0.2084474129 in his Appendix 7.1 (p41).
#' For this reason, the frequency comparison is relaxed to a `tol`
#' value of `1e-7` in a portion of the oce test suite
#' (see `tests/testthat/test_tidem.R` in the source).
#'
#' *Comments on phases 2 and 3*
#'
#' A specific example may be of help in understanding the removal of unresolvable
#' constituents. For example, the `data(sealevel)` dataset is of length
#' 6718 hours, and this is too short to resolve the full list of constituents,
#' with the conventional (and, really, necessary) limit of `rc=1`.
#' From Table 1 of Foreman (1978), this timeseries is too short to resolve the
#' `SA` constituent, so that `SA` will not be in the resultant.
#' Similarly, Table 2 of Foreman (1978) dictates the removal of
#' `PI1`, `S1` and `PSI1` from the list. And, finally,
#' Table 3 of Foreman (1978) dictates the removal of
#' `H1`, `H2`, `T2` and `R2`, and since that document
#' suggests that `GAM2` be subsumed into `H1`,
#' then if `H1` is already being deleted, then
#' `GAM2` will also be deleted.
#'
#' @param t either a vector of times or a [sealevel-class] object
#' (as created with [read.sealevel()] or [as.sealevel()]).
#' In the first case, `x` must be provided. In the second
#' case, though, any `x` that is provided will be ignored,
#' because sealevel objects contain both `time` and water
#' `elevation`, and the latter is used for `x`.
#'
#' @param x an optional numerical vector holding something that varies with
#' time. This is ignored if `t` is a [sealevel-class] object,
#' because it is inferred automatically, using `t[["elevation"]]`.
#'
#' @param constituents an optional character vector holding the names
#' of tidal constituents to which the fit is done; see \dQuote{Details}
#' and \dQuote{Constituent Naming Convention}.
#'
#' @template tideconst
#'
#' @param infer a list of constituents to be inferred from
#' fitted constituents according to the method outlined
#' in Section 2.3.4 of Foreman (1978).
#' If `infer` is `NULL`, the default, then
#' no such inferences are made. Otherwise, some constituents
#' are computed based on other constituents, instead of being
#' determined by regression at the proper frequency.
#' If provided, `infer` must be a list containing
#' four elements:
#' `name`, a vector of strings naming the constituents to be
#' inferred; `from`, a vector of strings naming the fitted
#' constituents used as the sources for those inferences (these
#' source constituents are added to the regression list, if they
#' are not already there);
#' `amp`, a numerical vector of factors to be applied to the
#' source amplitudes; and `phase`, a numerical vector of angles,
#' in degrees, to be subtracted from the source phases. For example,
#' Following Foreman (1998), if any of the `name` items
#' have already been computed, then the suggested inference is ignored,
#' and the already-computed values are used.
#' \preformatted{
#' infer=list(name=c("P1","K2"),
#' from=c("K1", "S2"),
#' amp=c(0.33093, 0.27215),
#' phase=c(-7.07, -22.4)
#' }
#' means that the amplitude of `P1` will be set as 0.33093 times the calculated amplitude
#' of `K1`, and that the `P1` phase will be set to the `K1` phase,
#' minus an offset of `-7.07` degrees.
#' (This example is used in the Foreman (1978) discussion of a
#' Fortran analysis code and also in Pawlowicz et al. (2002) discussion
#' of the T_TIDE Matlab code.
#' Rounded to the 0.1mm resolution of values reported in Foreman (1978)
#' and Pawlowicz et al. (2002),
#' the `tidem` results have root-mean-square amplitude difference
#' to Foreman's (1978) Appendix 7.3 of 0.06mm; by comparison,
#' the results in Table 1 of Pawlowicz et al. (2002) agree with Foreman's
#' results to RMS difference 0.04mm.)
#'
#' @param latitude if provided, the latitude of the observations. If not
#' provided, `tidem` will try to infer this from the first argument,
#' if it is a [sealevel-class] object.
#'
#' @param rc the value of the coefficient in the Rayleigh criterion.
#'
#' @param regress function to be used for regression, by default
#' [lm()], but could be for example `rlm` from the
#' `MASS` package.
#'
#' @template debugTemplate
#'
#' @return An object of [tidem-class], consisting of
#' \item{const}{constituent number, e.g. 1 for `Z0`, 1 for `SA`,
#' etc.} \item{model}{the regression model} \item{name}{a vector of constituent
#' names, in non-subscript format, e.g. "`M2`".} \item{frequency}{a vector
#' of constituent frequencies, in inverse hours.} \item{amplitude}{a vector of
#' fitted constituent amplitudes, in metres.} \item{phase}{a vector of fitted
#' constituent phase. NOTE: The definition of phase is likely to change as
#' this function evolves. For now, it is phase with respect to the first data
#' sample.} \item{p}{a vector containing a sort of p value for each
#' constituent. This is calculated as the average of the p values for the
#' sine() and cosine() portions used in fitting; whether it makes any sense is
#' an open question.}
#'
#' @section Application to non-hourly data:
#'
#' The framework on which [tidem()] rests on the assumption of data
#' that have been sampled on a 1-hour interval (see e.g. Foreman, 1977).
#' Since regression (as opposed to spectral analysis) is used to infer
#' the amplitude and phase of tidal constituents, data gaps do not pose
#' a serious problem. Sampling intervals under an hour are also not a
#' problem. However, trying to use [tidem()] on time series that are
#' sampled at uniform intervals that exceed 1 hour can lead to results
#' that are difficult to interpret. For example, some drifter data are
#' sampled at a 6-hour interval. This makes it impossible to fit for the
#' S4 component (which has exactly 6 cycles per day), because the method
#' works by constructing sine and cosine series at tidal frequencies and
#' using these as the basis for regression. Each of these series will have
#' a constant value through the constructed time, and regression cannot handle
#' that (in addition to a constant-value constructed series that is used to fit
#' for the Z0 constituent). [tidem()] tries to handle such problems by examining
#' the range of the constructed sine and cosine time-series, omitting any
#' constituents that yield near-constant values in either of these. Messages are
#' issued if this problem is encountered. This prevents failure of the regression,
#' and the predictions of the regression seem to represent the data reasonably well,
#' but the inferred constituent amplitudes are not physically reasonable. Cautious
#' use of [tidem()] to infer individual constituents might be warranted, but
#' users must be aware that the results will be difficult to interpret. The tool
#' is simply not designed for this use.
#'
#' @section Bugs:
#'
#' \enumerate{
#' \item This function is not fully developed yet, and both the
#' form of the call and the results of the calculation may change.
#'
#' \item The reported `p` value may make no sense at all, and it might be
#' removed in a future version of this function. Perhaps a significance level
#' should be presented, as in the software developed by both Foreman and
#' Pawlowicz.
#'
#' }
#'
#' @section Constituent Naming Convention:
#'
#' `tidem` uses constituent names that follow the convention
#' set by Foreman (1978). This convention is slightly different
#' from that used in the T-TIDE package of Pawlowicz et al.
#' (2002), with Foreman's `UPS1` and `M8` becoming
#' `UPSI` and `MS` in T-TIDE. To permit the use of either notation,
#' [tidem()] uses [tidemConstituentNameFix()] to
#' convert from T-TIDE names to the
#' Foreman names, issuing warnings when doing so.
#'
#' @section Agreement with `T_TIDE` results:
#'
#' The `tidem` amplitude and phase results, obtained with
#' \preformatted{
#' tidem(sealevelTuktoyaktuk, constituents=c("standard", "M10"),
#' infer=list(name=c("P1", "K2"),
#' from=c("K1", "S2"),
#' amp=c(0.33093, 0.27215),
#' phase=c(-7.07, -22.40)))
#' }
#' match the `T_TIDE` values listed in
#' Table 1 of Pawlowicz et al. (2002),
#' after rounding amplitude and phase to 4 and 2 digits past
#' the decimal place, respectively, to match the format of
#' that table.
#'
#' @author Dan Kelley
#'
#' @references
#'
#' Foreman, M G., 1977 (revised 1996).
#' Manual for Tidal Heights Analysis and Prediction.
#' Pacific Marine Science Report 77-10.
#' British Columbia, Canada: Institute of Ocean Sciences, Patricia Bay.
#'
#' Foreman, M. G. G., 1978.
#' Manual for Tidal Currents Analysis and Prediction.
#' Pacific Marine Science Report 78-6.
#' British Columbia, Canada: Institute of Ocean Sciences, Patricia Bay,
#'
#' Foreman, M. G. G., Neufeld, E. T., 1991.
#' Harmonic tidal analyses of long time series.
#' International Hydrographic Review, 68 (1), 95-108.
#'
#' Leffler, K. E. and D. A. Jay, 2009. Enhancing tidal harmonic analysis:
#' Robust (hybrid) solutions. Continental Shelf Research, 29(1):78-88.
#'
#' Pawlowicz, Rich, Bob Beardsley, and Steve Lentz, 2002. Classical tidal
#' harmonic analysis including error estimates in MATLAB using `T_TIDE`.
#' Computers and Geosciences, 28, 929-937.
#'
#' @examples
#' library(oce)
#' # The demonstration time series from Foreman (1978),
#' # also used in T_TIDE (Pawlowicz, 2002).
#' data(sealevelTuktoyaktuk)
#' tide <- tidem(sealevelTuktoyaktuk)
#' summary(tide)
#'
#' # AIC analysis
#' extractAIC(tide[["model"]])
#'
#' # Fake data at M2
#' library(oce)
#' data("tidedata")
#' M2 <- with(tidedata$const, freq[name == "M2"])
#' t <- seq(0, 10 * 86400, 3600)
#' eta <- sin(M2 * t * 2 * pi / 3600)
#' sl <- as.sealevel(eta)
#' m <- tidem(sl)
#' summary(m)
#'
#' @family things related to tides
tidem <- function(
t, x, constituents, infer = NULL, latitude = NULL,
rc = 1, regress = lm, debug = getOption("oceDebug")) {
oceDebug(debug, "tidem(t, x,\n", sep = "", unindent = 1)
oceDebug(debug, " constituents=",
if (missing(constituents)) {
"(missing)"
} else {
paste("c(\"", paste(constituents, collapse = "\", \""), "\")", sep = "")
},
",\n",
sep = "", unindent = 1
)
oceDebug(debug, " latitude=", if (is.null(latitude)) "NULL" else latitude, ",\n", sep = "", unindent = 1)
oceDebug(debug, " rc=", rc, ",\n", sep = "", unindent = 1)
oceDebug(debug, " debug=", debug, ") {\n", sep = "", unindent = 1)
cl <- match.call()
if (missing(t)) {
stop("must supply 't', either a vector of times or a sealevel object")
}
if (inherits(t, "sealevel")) {
sl <- t
t <- sl[["time"]]
x <- sl[["elevation"]]
if (is.null(latitude)) {
latitude <- sl[["latitude"]]
}
} else {
if (missing(x)) {
stop("must supply 'x', since the first argument is not a sealevel object")
}
if (inherits(x, "POSIXt")) {
warning("tidem() switching first 2 args to permit old-style usage")
tmp <- x
x <- t
t <- tmp
}
if (length(x) != length(t)) {
stop("lengths of 'x' and 't' must match, but they are ", length(x), " and ", length(t), " respectively")
}
if (inherits(t, "POSIXt")) {
t <- as.POSIXct(t)
} else {
stop("t must be a vector of POSIXt times")
}
sl <- as.sealevel(x, t)
}
# Check infer extensively, to prevent weird errors for e.g. an improperly-named
# constituent.
data("tidedata", package = "oce", envir = environment())
tidedata <- get("tidedata") # , pos=globalenv())
tc <- tidedata$const
ntc <- length(tc$name)
# 'infer' sanity-check
if (!is.null(infer)) {
if (!is.list(infer)) {
stop("infer must be a list")
}
if (length(infer) != 4) {
stop("infer must hold 4 elements")
}
if (!all.equal(sort(names(infer)), c("amp", "from", "name", "phase"))) {
stop("infer must contain 'name', 'from', 'amp', and 'phase', and nothing else")
}
if (!is.character(infer$name)) {
stop("infer$name must be a vector of character strings")
}
infer$name <- tidemConstituentNameFix(infer$name)
if (!is.character(infer$from)) {
stop("infer$from must be a vector of character strings")
}
infer$from <- tidemConstituentNameFix(infer$from)
if (length(infer$name) != length(infer$from)) {
stop("lengths of infer$name and infer$from must be equal")
}
if (length(infer$name) != length(infer$amp)) {
stop("lengths of infer$name and infer$amp must be equal")
}
if (length(infer$name) != length(infer$phase)) {
stop("lengths of infer$name and infer$phase must be equal")
}
for (n in infer$name) {
if (!(n %in% tc$name)) {
stop("infer$name value \"", n, "\" is not a known tidal constituent; see const$name in ?tidedata")
}
}
for (n in infer$from) {
if (!(n %in% tc$name)) {
stop("infer$from value \"", n, "\" is not a known tidal constituent; see const$name in ?tidedata")
}
}
} # 'infer' sanity-check
# The arguments seem to be OK, so start the actual analysis now.
startTime <- t[1]
endTime <- tail(t, 1)
years <- as.numeric(difftime(endTime, startTime, units = "secs")) / 86400 / 365.25
if (years > 18.6) {
warning("Time series spans 18.6 years, but tidem() is ignoring this important fact")
}
standard <- tc$ikmpr > 0
addedConstituents <- NULL
if (missing(constituents)) {
# Default 'name', 'freq', 'kmpr' and 'indices'; drop Z0 because we infer it directly.
name <- tc$name[standard]
freq <- tc$freq[standard]
kmpr <- tc$kmpr[standard]
indices <- seq(1:ntc)[standard] # NB. Z0 need not be dropped; we work with indices later anyway
oceDebug(debug, "starting with ", length(name), " default constituents: ", paste(name, collapse = " "), sep = "", "\n")
} else {
# Build up 'name'; later, infer 'indices' and thereby 'freq' and 'kmpr'.
name <- NULL
for (i in seq_along(constituents)) {
if (constituents[i] == "standard") {
# must be first!
if (i != 1) {
stop("\"standard\" must occur first in constituents list")
}
name <- tc$name[standard]
} else {
constituents <- tidemConstituentNameFix(constituents)
if (substr(constituents[i], 1, 1) == "-") {
# Case 1: removal. Require a valid name, and warn if not in the list already.
nameRemove <- substr(constituents[i], 2, nchar(constituents[i]))
if (1 != sum(tc$name == nameRemove)) {
stop(
"\"", nameRemove, "\" is not a known tidal constituent; try one of: ",
paste(tc$name, collapse = "\" \"")
)
}
remove <- which(name == nameRemove)
oceDebug(debug > 1, "removed \"", nameRemove, "\"\n", sep = "")
if (0 == length(remove)) {
warning("\"", nameRemove, "\" is not in the list of constituents currently under study", sep = "")
} else {
name <- name[-remove]
}
} else {
# Case 2: addition. Require a valid name, and ignore repeat requests.
add <- which(tc$name == constituents[i])
if (1 != length(add)) {
stop("\"", constituents[i], "\" is not a known tidal constituent (line 1093)")
}
if (!(constituents[i] %in% name)) {
name <- c(name, tc$name[add])
addedConstituents <- c(addedConstituents, add)
}
}
}
}
}
oceDebug(debug, "will fit for ", length(name), " constituents: ", paste(name, collapse = " "), "\n", sep = "")
# We let users add "Z0" as a constituent, but we remove it now since the
# regression will have an intercept and that becomes Z0.
fitForZ0 <- "Z0" %in% name
oceDebug(debug, "fitForZ0=", fitForZ0, "\n")
# All of the names should be valid from the above, but to protect against code changes,
# we check to be sure.
if (any(!(name %in% tc$name))) {
bad <- NULL
for (n in name) {
if (!(n %in% tc$name)) {
bad <- c(bad, n)
}
}
stop("unknown constituent names: ", paste(bad, collapse = " "), " (please report this error to developer)")
}
# Infer indices from the names, sort them as in the tidal-constituent
# table, tc, and then look up freq and kmpr from that table.
indices <- sort(unlist(lapply(name, function(name) which(tc$name == name))))
name <- tc$name[indices]
freq <- tc$freq[indices]
kmpr <- tc$kmpr[indices]
nc <- length(name)
# Check Rayleigh criterion
interval <- diff(range(as.numeric(sl@data$time), na.rm = TRUE)) / 3600 # in hours
oceDebug(debug, "interval=", interval, " hours\n")
dropTerm <- NULL
for (i in 1:nc) {
cc <- which(tc$name == kmpr[i])
if (length(cc)) {
cannotFit <- (interval * abs(freq[i] - tc$freq[cc])) < rc
oceDebug(debug, "i=", i, ", name=", name[i], ", kmpr[", i, "]=", kmpr[i], ", cannotFit=", cannotFit, "\n", sep = "")
if (cannotFit) {
dropTerm <- c(dropTerm, i)
}
}
}
oceDebug(debug, "before trimming constituents for Rayleigh condition, name[1:", length(name), "]=", paste(name, collapse = " "), sep = "", "\n")
if (length(dropTerm) > 0) {
# Bookmark 1A (see also 1B: link up variables)
message("Note: the tidal record is too short to fit for constituents: ", paste(name[dropTerm], collapse = ", "))
indices <- indices[-dropTerm]
name <- name[-dropTerm]
freq <- freq[-dropTerm]
kmpr <- kmpr[-dropTerm]
}
oceDebug(debug, "after trimming constituents for Rayleigh condition, name[1:", length(name), "]=", paste(name, collapse = " "), sep = "", "\n")
# Ensure that any added constituents are in the list, i.e. prevent
# the Rayleigh criterion from trimming them. (Before work on
# issue 1350, they would simply be dropped if they failed the Rayleigh
# criterion. Although that was a sensible choice, it was decided
# on 2017-12-27, whilst working on issue 1350, to make tidem() do the
# the same thing as the Foreman 1978 code as exemplified in his appendices
# 7.2 and 7.3.)
if (length(addedConstituents)) {
oceDebug(debug, "addedConstituents=", paste(addedConstituents, collapse = " "), "\n")
for (a in addedConstituents) {
# Avoid adding constituents that we already have
if (!(tc$name[a] %in% name)) {
if (debug > 0) {
message("ADDING:")
message(" tc$name[a=", a, "]=\"", tc$name[a], "\"", sep = "")
message(" tc$freq[a=", a, "]=\"", tc$freq[a], "\"", sep = "")
message(" tc$kmpr[a=", a, "]=\"", tc$kmpr[a], "\"", sep = "")
}
indices <- c(indices, which(tc$name == name[a]))
name <- c(name, tc$name[a])
freq <- c(freq, tc$freq[a])
kmpr <- c(kmpr, tc$kmpr[a])
}
}
}
oceDebug(debug, "after adding new constituents, ", vectorShow(name))
oceDebug(debug, "after adding new constituents, ", vectorShow(freq))
# Ensure that we fit for any infer$from constituents, *regardless* of whether
# those consitituents are permitted by the Rayleigh criterion.
if (!is.null(infer)) {
for (n in c(infer$from)) {
if (!(n %in% name)) {
a <- which(tc$name == n)
indices <- c(indices, a)
name <- c(name, tc$name[a])
freq <- c(freq, tc$freq[a])
kmpr <- c(kmpr, tc$kmpr[a])
message("fitting for infer$from=", n, ", even though the Rayleigh Criterion would exclude it")
}
}
}
oceDebug(debug, "after handling 'infer', ", vectorShow(name))
oceDebug(debug, "after handling 'infer', ", vectorShow(freq))
# sort constituents by index (which, among other things, ensures that Z0 is at the start, if it exists)
oindices <- order(indices)
indices <- indices[oindices]
name <- name[oindices]
freq <- freq[oindices]
kmpr <- kmpr[oindices]
nc <- length(name)
oceDebug(debug, "after reordering indices, ", vectorShow(name))
oceDebug(debug, "after reordering indices, ", vectorShow(freq))
rm(oindices) # clean up namespace
if (0L == nc) {
stop("cannot fit for any constituents")
}
elevation <- sl[["elevation"]]
time <- sl[["time"]]
nt <- length(elevation)
x <- array(dim = c(nt, 2 * nc))
oceDebug(debug, vectorShow(nc))
oceDebug(debug, vectorShow(dim(x)))
x[, 1] <- rep(1, nt)
pi <- 4 * atan2(1, 1)
rpd <- atan2(1, 1) / 45 # radians per degree
# tRef <- ISOdate(1899, 12, 31, 12, 0, 0, tz="UTC") # was this ever used?
# tRef <- centralTime # used previous to "dk" branch early 2018
tRef <- numberAsPOSIXct(3600 * round(mean(as.numeric(time, tz = "UTC")) / 3600), tz = "UTC")
hour2pi <- 2 * pi * (as.numeric(time) - as.numeric(tRef)) / 3600
oceDebug(debug, "tRef=", tRef, ", nc=", nc, ", length(name)=", length(name), "\n")
# The sameCriterion was added in January 2023, after this portion of tidem()
# had been stable for a decade or more. Also, I am trying to make tidem()
# do something that (I think) neither Foreman's code nor t-tide does.
# Therefore, I am putting a long comment here! Please also see
# https://github.com/dankelley/oce/issues/2034 to learn where this idea of
# using sameCriterion got started.
#
# In the code below, we keep track of any constituents for which the
# constructed C and S vectors will be constant (or nearly so). This came up
# because I was trying to analyse data sampled on a 6-hour interval. That
# works out to match the S4 period (approximately, but sampling times are
# not given to an infinite number of digits and neither is the S4 period in
# this code). Imagine a perfect match. Then, S and C would both be constant
# across the times used in the loop. And that's bad because the regression
# already has a constant column (for Z0). The lm() computation then gives
# NA values for coefficients, p value, etc. for the repeat. (Luckily, it
# doesn't just say the matrix cannot be inverted and die!) The approximate
# criterion is based on a test case in the context of the local machine
# epsilon. I use a criterion of 0.01 because for "good" cases (and my
# 6-hour test file for drifter data) have C and S span of very nearly 2, but for "bad" cases, it is
# of order e-7 or so. Any dividing line would do, I think, but maybe for short records
# the span might not get to be -1 to +1 and so I am choosing 1e-2 as a criterion. This value
# might need to be revisited.
iBad <- NULL
icBad <- NULL
sdCriterion <- 1e-2
oceDebug(debug, vectorShow(sdCriterion))
# 20230122 danS <- danC <- NULL
for (i in 1:nc) {
# 20230122 oceDebug(debug+1, "setting ", i, "-th coefficient (name=", name[i], " period=", 1/freq[i], " h)", "\n", sep="")
ft <- freq[i] * hour2pi
C <- cos(ft)
S <- sin(ft)
sdS <- sd(S)
sdC <- sd(C)
# 20230122 danS <- c(danS, sdS) # FIXME: remove
# 20230122 danC <- c(danC, sdC) # FIXME: remove
# 20230122 oceDebug(debug+1, sprintf(" sdC %.4g; sdS %.4g\n", sdC, sdS))
# Find whether anything is uncomputable. We ignore Z0 because that is handled later.
if (name[i] != "Z0" && (sdS < sdCriterion || sdC < sdCriterion)) {
oceDebug(debug, " ** uncomputable at ", name[i], " (period ", 1 / freq[i], "h) **\n")
icBad <- c(icBad, 1 + 2 * (i - 1))
icBad <- c(icBad, 2 + 2 * (i - 1))
iBad <- c(iBad, i)
}
x[, 1 + 2 * (i - 1)] <- C
x[, 2 + 2 * (i - 1)] <- S
}
name2 <- matrix(rbind(paste(name, "_C", sep = ""), paste(name, "_S", sep = "")), nrow = length(name), ncol = 2)
dim(name2) <- c(2 * length(name), 1)
colnames(x) <- name2
#<< cat("The following should be omitted: ", vectorShow(name2[uncomputable], n=20))
#<< cat("Their indices are", vectorShow(uncomputable, n=20))
oceDebug(debug, "cleaning up 'x' matrix prior to doing regression\n")
# Remove problematic constutuents; see https://github.com/dankelley/oce/issues/2034
if (length(iBad)) {
message("Note: the sampling interval is too coarse to fit for constituents: ", paste(name[iBad], collapse = ", "))
# 20230122 browser()
# 20230122 A<-colnames(x)
# 20230122 B<-name2
# 20230122 C<-name
indices <- indices[-iBad]
name <- name[-iBad]
freq <- freq[-iBad]
kmpr <- kmpr[-iBad]
nc <- nc - length(iBad)
oceDebug(debug, "Removing: ", paste(name2[icBad], collapse = ", "))
x <- x[, -icBad]
name2 <- name2[-icBad]
# 20230122 # A check with 6-h case (sandbox/dk/tidem/04_restrict_coefficients.R)
# 20230122 print(A[!A%in%colnames(x)]) # should be S2_C S2_S S4_C S4_S
# 20230122 print(B[!B%in%name2]) # should be S2_C S2_S S4_C S4_S
# 20230122 print(C[!C%in%name]) # should be S2 S4
}
# Remove the sine() part of the Z0 constituent, which makes no sense for a constant.
if ("Z0_S" %in% colnames(x)) {
x <- x[, -which("Z0_S" == colnames(x))]
oceDebug(debug, "model has Z0, so trimming the sin(freq*time) column\n")
}
oceDebug(debug, "about to do regression\n")
oceDebug(debug, vectorShow(colnames(x)))
model <- regress(elevation ~ x - 1, na.action = na.exclude)
#>> browser()
if (debug > 0) {
cat("regression worked OK; the results are as follows:\n")
print(summary(model))
}
coef <- model$coefficients
p.all <- if (4L == ncol(summary(model)$coefficients)) {
summary(model)$coefficients[, 4L]
} else {
rep(NA, length = 1L + nc)
}
amplitude <- phase <- p <- vector("numeric", length = nc)
oceDebug(debug, vectorShow(nc))
oceDebug(debug, vectorShow(phase))
oceDebug(debug, vectorShow(name))
ic <- 1
for (i in seq_len(nc)) {
# Z0 has Z0_C but no Z0_S (since zero is a degenerate regression variable)
if (name[i] == "Z0") {
if (i != 1) {
stop("Z0 should be at the start of the regression coefficients. Please report this to developer.")
}
j <- which(tidedata$const$name == name[i])
vuf <- tidemVuf(tRef, j = j, latitude = latitude)
amplitude[i] <- coef[ic]
phase[i] <- 0
p[i] <- p.all[ic]
oceDebug(debug, "processed coefs at i=", i, ", ic=", ic,
", name=", name[i], ", f=", vuf$f, ", angle adj=", (vuf$u + vuf$v) * 360,
", amplitude=", amplitude[i], ", phase=", phase[i], ", p=", p[i], "\n",
sep = ""
)
ic <- ic + 1 # only skip forward 1 since Z0 takes 1 column (contrast below)
} else {
C <- coef[ic] # coefficient on cos(t)
S <- coef[ic + 1] # coefficient on sin(t)
amplitude[i] <- sqrt(S^2 + C^2)
# Calculate phase from the coefficients on sin() and cos(). Generally,
# cos(t - phase) == cos(phase)*cos(t) + sin(phase)*sin(t)
# By the definition of the regression model, we have
# cos(t - phase) == c * cos(t) + s * sin(t)
# and thus phase is defined by
# tan(phase) == s/c
phase[i] <- atan2(S, C)
# Adjust amplitude phase, as in ~/src/foreman/tide12_r2.f:405
j <- which(tidedata$const$name == name[i])
vuf <- tidemVuf(tRef, j = j, latitude = latitude)
amplitude[i] <- amplitude[i] / vuf$f
p[i] <- 0.5 * (p.all[ic + 1] + p.all[ic])
#<<cat("i=", i, ", ic=", ic, ", name=\"", name[i], "\", p.all[ic]=",
# p.all[ic], ", p.all[ic+1]=", p.all[ic+1], ", we set p[i]=", p[i],
# " (the mean)\n", sep="")
oceDebug(debug, "processed coefs at i=", i, ", ic=", ic, ", name=", name[i],
", S=", S, ", C=", C, ", f=", vuf$f, ", angle adj=", (vuf$u + vuf$v) * 360, ",
amplitude=", amplitude[i], ", phase=", phase[i], ", p=", p[i], "\n",
sep = ""
)
ic <- ic + 2 # skip forward 2 since non-Z0 takes 2 columns (contrast above)
}
}
oceDebug(debug, vectorShow(phase))
phase <- phase * 180 / pi
phase <- ifelse(phase < -360, 720 + phase, phase)
phase <- ifelse(phase < 0, 360 + phase, phase)
# Do Greenwich phase correction, if `infer` is TRUE
C <- unlist(lapply(name, function(n) which(n == tidedata$const$name)))
vuf <- tidemVuf(tRef, j = C, latitude = latitude)
oceDebug(debug, vectorShow(freq))
oceDebug(debug, vectorShow(phase))
oceDebug(debug, vectorShow(vuf$u))
oceDebug(debug, vectorShow(vuf$v))
phase <- phase + (vuf$v + vuf$u) * 360
phase <- ifelse(phase < 0, phase + 360, phase)
phase <- ifelse(phase > 360, phase - 360, phase)
# Handle (optional) inferred constituents. We know that
# this list is well-formed because of extensive tests near
# the start of this function.
if (!is.null(infer)) {
if (debug > 0) {
cat("BEFORE inference:\n")
print(data.frame(name = name, freq = round(freq, 6), amplitude = round(amplitude, 4)))
}
for (n in seq_along(infer$name)) {
oceDebug(debug, "n=", n, "; handling inferred constituent ", infer$name[n], "\n")
iname <- which(tc$name == infer$name[n])[1]
oceDebug(debug, "infer$name[\"", n, "\"]=\"", infer$name[n], "\" yields iname=", iname, "\n", sep = "")
oceDebug(debug, "iname=", iname, "\n")
# NOTE: we know that infer$from has been fitted for, because we forced it to be,
# irrespective of the Rayleight Criterion. Still, we test here, in case the code
# gets changed later.
if (infer$from[n] %in% name) {
ifrom <- which(name == infer$from[n])[1]
if (infer$name[n] %in% name) {
# message("name already in list")
iname <- which(name == infer$name[n])[1]
# Update, skipping 'indices', 'name' and 'freq', since they are already OK.
amplitude[iname] <- infer$amp[n] * amplitude[ifrom]
phase[iname] <- phase[ifrom] - infer$phase[n]
p[iname] <- p[ifrom]
oceDebug(debug, "replace existing ", name[iname], " based on ", name[ifrom], " (", freq[ifrom], " cph)\n", sep = "")
warning(
"inferring \"", infer$name[n],
"\" which is already included in the regression. Foreman says to skip it; unsure on what T_TIDE does\n"
)
} else {
# We append new values at the end, knowing that they will get
# shifted back to their proper positions when we reorder the
# whole solution, after handling these inferences.
#
# The first step is to adjust the amp and phase of infer$from; this
# is done based on formulae in Foreman (1978) sec 2.3.4. It looks
# as though t_tide.m on and after about line 472 is doing a
# similar thing, although the numbers do not agree exactly,
# as shown in issue 1351, code 1351c.R.
#
# Foreman (1978) sec 2.3.4.
# Notation: suffices "1" and "2" refer to "from" and "name" here.
i1 <- which(tc$name == infer$from[n])[1]
i2 <- which(tc$name == infer$name[n])[1]
oceDebug(
debug, "tRef=", format(tRef, "%Y-%m-%d %H:%M:%S"),
", i1=", i1, ", i2=", i2, ", lat=", latitude, "\n"
)
vuf12 <- tidemVuf(tRef, c(i1, i2), latitude = latitude)
# vuf2 <- tidemVuf(tRef, i2, latitude=latitude)
f1 <- vuf12$f[1]
f2 <- vuf12$f[2]
oceDebug(debug, "f1=", f1, ", f2=", f2, "\n")
# FIXME: what is unit of u and v? t_tide.m:482 suggests it is degrees
# Foreman's tide12_r2.f:399 suggests U and V are in cycles,
# and this is consistent with Pawlowicz's t_tide.m:451
# We convert vu1 and vu2 to be in degrees, as t_tide.m does
vu1 <- (vuf12$v[1] + vuf12$u[1]) * 360
vu2 <- (vuf12$v[2] + vuf12$u[2]) * 360
oceDebug(debug, "vu1=", vu1, ", vu2=", vu2, "\n")
sigma1 <- tc$freq[i1]
sigma2 <- tc$freq[i2]
oceDebug(debug, "sigma1=", sigma1, ", sigma2=", sigma2, "\n")
# tmp is pi*N*(sigma2-sigma1) in Foreman
tmp <- pi * interval * (sigma2 - sigma1)
r12 <- infer$amp[n]
# FIXME: sign for Foreman?
zeta <- infer$phase[n]
S <- r12 * (f2 / f1) * sin(tmp) * sin(rpd * (vu2 - vu1 + zeta)) / tmp
C <- 1 + r12 * (f2 / f1) * sin(tmp) * cos(rpd * (vu2 - vu1 + zeta)) / tmp
oceDebug(debug, "tmp=", tmp, ", S=", S, ", C=", C, ", sqrt(S^2+C^2)=", sqrt(S^2 + C^2), "\n")
oceDebug(debug, infer$from[n], "amplitude, old=", amplitude[ifrom], ", new=", amplitude[ifrom] / sqrt(S^2 + C^2), "\n")
amplitude[ifrom] <- amplitude[ifrom] / sqrt(S^2 + C^2)
oceDebug(debug, infer$from[n], "phase, old=", phase[ifrom], ", new=", phase[ifrom] + atan2(S, C) / rpd, "\n")
phase[ifrom] <- phase[ifrom] + atan2(S, C) / rpd
# End of Foreman 1978 inference calculation. Now we can define 'name' i.t.o. 'from'
iname <- which(tc$name == infer$name[n])[1]
oceDebug(debug, "Below is inference for ", infer$name[n], " (index=", iname, ")\n")
indices <- c(indices, iname)
name <- c(name, infer$name[n])
freq <- c(freq, tc$freq[iname])
amplitudeInferred <- infer$amp[n] * amplitude[ifrom]
phaseInferred <- phase[ifrom] - infer$phase[n]
oceDebug(debug, " ", infer$name[n], "inferred amplitude=", amplitudeInferred, "\n")
oceDebug(debug, " ", infer$name[n], "inferred phase=", phaseInferred, "\n")
amplitude <- c(amplitude, amplitudeInferred)
phase <- c(phase, phaseInferred)
p <- c(p, p[ifrom])
oceDebug(debug, " create ", infer$name[n], " (index=", iname, ", ",
tc$freq[iname], " cph) based on ", name[ifrom],
" (index ", ifrom, ", ", freq[ifrom], " cph)\n",
sep = ""
)
}
} else {
stop("Internal error (please report): cannot infer ", infer$name[n], " from ", infer$from[n], " because the latter was not computed")
}
}
# reorder by original position in tc
o <- order(indices)
indices <- indices[o]
stopifnot(length(o) == length(name))
stopifnot(length(o) == length(freq))
stopifnot(length(o) == length(amplitude))
stopifnot(length(o) == length(phase))
stopifnot(length(o) == length(p))
name <- name[o]
freq <- freq[o]
amplitude <- amplitude[o]
phase <- phase[o]
p <- p[o]
rm(o)
if (debug > 0) {
cat("AFTER inference\n")
print(data.frame(name = name, freq = round(freq, 5), amplitude = round(amplitude, 4)))
}
}
phase <- phase %% 360
res <- new("tidem")
o <- order(freq)
res@data <- list(
model = model,
call = cl,
tRef = tRef,
const = indices[o],
name = name[o],
freq = freq[o],
amplitude = amplitude[o],
phase = phase[o],
p = p[o]
)
res@metadata$rc <- rc
res@metadata$version <- 2
res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep = "", collapse = ""))
oceDebug(debug, "} # tidem()\n", sep = "", unindent = 1)
res
}
#' Predict a Tidal Signal
#'
#' This creates a time-series of predicted tides, based on a
#' tidal model object that was created by [as.tidem()] or [tidem()].
#'
#' All the tidal constituents that are stored in `object` are used,
#' not just those that are statistically significant or that have amplitude
#' exceeding any particular value. In this respect, [predict.tidem()]
#' follows a pattern established by e.g. [predict.lm()]. Note that
#' the constituents in `object` are straightforward if it
#' was constructed with [as.tidem()], but considerably more complicated
#' for [tidem()], and so the documentation for the latter ought to
#' be studied closely, especially with regard to the Rayleigh criterion.
#'
#' @param object a [tidem-class] object.
#'
#' @param newdata vector of POSIXt times at which to make the
#' prediction. For models created with [tidem()],
#' the `newdata` argument is optional, and if it is not provided, then
#' the predictions are at the observation times given to
#' [tidem()]. However, `newdata` is required if [as.tidem()]
#' had been used to create `object`.
#'
#' @param \dots optional arguments passed on to children.
#'
#' @return A vector of predictions.
#'
#' @examples
#'
#' # Show non-tidal sealevel signal in Halifax Harbour during
#' # the year 2002. The spike resulted from Hurricane Juan.
#' library(oce)
#' data(sealevel)
#' time <- sealevel[["time"]]
#' elevation <- sealevel[["elevation"]]
#' prediction <- tidem(sealevel) |> predict()
#' oce.plot.ts(time, elevation - prediction)
#'
#' @section Sample of Usage:
#' \preformatted{
#' # prediction at specified times
#' data(sealevel)
#' m <- tidem(sealevel)
#' # Check fit over 2 days (interpolating to finer timescale)
#' look <- 1:48
#' time <- sealevel[["time"]]
#' elevation <- sealevel[["elevation"]]
#' oce.plot.ts(time[look], elevation[look])
#' # 360s = 10 minute timescale
#' t <- seq(from=time[1], to=time[max(look)], by=360)
#' lines(t, predict(m, newdata=t), col="red")
#' legend("topright", col=c("black","red"),
#' legend=c("data","model"),lwd=1)
#' }
#'
#' @author Dan Kelley
#'
#' @family things related to tides
predict.tidem <- function(object, newdata, ...) {
dots <- list(...)
debug <- if ("debug" %in% names(dots)) dots$debug else 0
oceDebug(debug, "predict.tidem() {\n", sep = "", unindent = 1)
if (!missing(newdata) && !inherits(newdata, "POSIXt")) {
stop("newdata must be of class POSIXt")
}
version <- object@metadata$version
if (!is.null(version) && version == 3) {
oceDebug(debug, "object@metadata$version is 3, so assuming the object was created by as.tidem()\n")
if (missing(newdata)) {
stop("must supply newdata because object was created with as.tidem()")
}
hour2pi <- 2 * pi * (as.numeric(newdata) - as.numeric(object[["tRef"]])) / 3600
oceDebug(debug, vectorShow(hour2pi))
# message("head(hour2pi): ", paste(head(hour2pi), collapse=" "))
nc <- length(object@data$name)
res <- rep(0, length(hour2pi))
for (i in seq_len(nc)) {
oceDebug(debug, "accounting for constituent[", i, "] = ", object@data$name[i], "\n", sep = "")
omega.t <- object@data$freq[i] * hour2pi
a <- object@data$amplitude[i] * sin(2 * pi * object@data$phase[i] / 360)
b <- object@data$amplitude[i] * cos(2 * pi * object@data$phase[i] / 360)
res <- res + a * sin(omega.t) + b * cos(omega.t)
}
} else if (!is.null(version) && version == 2) {
oceDebug(debug, "object@metadata$version is 2, so assuming the object was created by tidem()\n")
if (!missing(newdata) && !is.null(newdata)) {
oceDebug(debug, "newdata was provided\n")
freq <- object@data$freq
name <- object@data$name
nc <- length(freq)
tt <- as.numeric(as.POSIXct(newdata, tz = "UTC"))
nt <- length(tt)
x <- array(dim = c(nt, 2 * nc))
x[, 1] <- rep(1, nt)
hour2pi <- 2 * pi * (as.numeric(tt) - as.numeric(object[["tRef"]])) / 3600
for (i in 1:nc) {
omega.t <- freq[i] * hour2pi
x[, 2 * i - 1] <- cos(omega.t)
x[, 2 * i] <- sin(omega.t)
}
colnames(x) <- matrix(rbind(paste(name, "_C", sep = ""), paste(name, "_S", sep = "")), nrow = length(name), ncol = 2)
if ("Z0_S" %in% colnames(x)) {
x <- x[, -which("Z0_S" == colnames(x))]
oceDebug(debug, "model has Z0, so trimming the sin(freq*time) column\n")
}
res <- as.numeric(predict(object@data$model, newdata = list(x = x), ...))
} else {
oceDebug(debug, "newdata was not provided\n")
res <- as.numeric(predict(object@data$model, ...))
}
} else {
if (!("version" %in% names(object@metadata))) {
warning("prediction is being made based on an old object; it may be wrong\n")
}
res <- as.numeric(predict(object@data$model, ...))
}
oceDebug(debug, "} # predict.tidem()\n", sep = "", unindent = 1)
res
}
#' Get a Tidal Prediction From a WebTide Database
#'
#' Get a tidal prediction from a WebTide database. This only
#' works if the standalone WebTide application is installed,
#' and if it is installed in a standard location. The details
#' of installation are not within the oce purview.
#'
#' There are two methods of using this function.
#' *Case 1:* `action="map"`. In this case, if
#' `plot` is `FALSE`, a list is returned, containing
#' all the `node`s in the selected database, along with all
#' the `latitude`s and `longitude`s. This value is
#' also returned (silently) if `plot` is true, but in that case,
#' a plot is drawn to indicate the node locations. If `latitude` and
#' `longitude` are given, then the node nearest that spot is indicated on
#' the map; otherwise, if `node` is given, then the location of that
#' node is indicated. There is also a special case: if `node` is negative
#' and `interactive()` is `TRUE`,
#' then [locator()] is called, and the node nearest the spot
#' where the user clicks the mouse is indicated in the plot and in the
#' return value.
#'
#' *Case 2:* `action="predict"`. If `plot` is `FALSE`,
#' then a list is returned, indicating `time`, predicted
#' `elevation`, velocity components `u` and `v`,
#' `node` number, the name of the `basedir`, and
#' the `region`. If `plot` is `TRUE`, this list is returned
#' silently, and time-series plots are drawn for elevation, u, and v.
#'
#' Naturally, `webtide` will not work unless WebTide has been installed on
#' the computer.
#'
#' @param action An indication of the action, either `action="map"` to
#' draw a map or `action="predict"` to get a prediction; see
#' \dQuote{Details}.
#'
#' @param longitude,latitude optional location at which prediction is required (ignored if
#' `node` is given).
#'
#' @param node optional integer relating to a node in the database. If `node`
#' is given, then neither `latitude` nor `longitude` may be given.
#' If `node` is positive, then specifies indicates the node. If it is negative,
#' [locator()] is called so that the user can click (once) on the map, after
#' which the node is displayed on the map.
#'
#' @param time a vector of times (in the UTC timezone)
#' at which prediction is to be made.
#' If not supplied, this will be the week starting at the present time,
#' computed with [presentTime()], with a 15 minute increment.
#'
#' @param basedir directory containing the `WebTide` application.
#'
#' @param region database region, given as a directory name in the WebTide
#' directory. For example, `h3o` is for Halifax Harbour, `nwatl` is
#' for the northwest Atlantic, and `sshelf` is for the Scotian Shelf and
#' Gulf of Maine.
#'
#' @param plot boolean indicating whether to plot.
#'
#' @param tformat optional argument passed to [oce.plot.ts()], for
#' plot types that call that function. (See [strptime()] for the
#' format used.)
#'
#' @template debugTemplate
#'
#' @param \dots optional arguments passed to plotting functions. A common
#' example is to set `xlim` and `ylim`, to focus a map region.
#'
#' @return The value depends on `action`:
#'
#' * If `action="map"` the return value is a
#' list containing the index of the nearest node, along with the
#' `latitude` and `longitude` of that node. If
#' `plot` is `FALSE`, this value is returned invisibly.
#'
#' * If `action="predict"`, the return value is a list containing a vector
#' of times (`time`), as well as vectors of the predicted `elevation`
#' in metres and the predicted horizontal components of velocity, `u` and
#' `v`, along with the `node` number, and the `basedir` and
#' `region` as supplied to this function. If `plot` is `FALSE`,
#' this value is returned invisibly.
#'
#' @source The WebTide software may be downloaded for free at the
#' Department of Fisheries and Oceans (Canada) website at
#' `http://www.bio.gc.ca/science/research-recherche/ocean/webtide/index-en.php`
#' (checked February 2016 and May 2017).
#'
#' @section Caution:
#' WebTide is not an open-source application, so the present function was
#' designed based on little more than guesses about the WebTide file structure.
#' Users should be on the lookout for odd results.
#'
#' @section Sample of Usage:
#' \preformatted{
#' # needs WebTide at the system level
#' library(oce)
#' # 1. prediction at Halifax NS
#' longitude <- -63.57
#' latitude <- 44.65
#' prediction <- webtide("predict", longitude=longitude, latitude=latitude)
#' mtext(paste0("prediction at ", latitude, "N and ", longitude, "E"), line=0.75, side=3)
#' # 2. map
#' webtide(lon=-63.57,lat=44.65,xlim=c(-64,-63),ylim=c(43.0,46))
#' }
#'
#' @author Dan Kelley
#'
#' @family things related to tides
webtide <- function(
action = c("map", "predict"),
longitude, latitude, node, time,
basedir = getOption("webtide"),
region = "nwatl",
plot = TRUE, tformat, debug = getOption("oceDebug"), ...) {
debug <- max(0, min(floor(debug), 2))
oceDebug(debug, "webtide(action=\"", action, "\", ...)\n",
sep = "", unindent = 1, style = "bold"
)
rpd <- atan2(1, 1) / 45 # radians per degree
action <- match.arg(action)
nodeGiven <- !missing(node)
longitudeGiven <- !missing(longitude)
latitudeGiven <- !missing(latitude)
path <- paste(basedir, "/data/", region, sep = "")
# 2016-02-03: it seems that there are several possibilities for this filename.
suffices <- c(".nod", "ll.nod", "_ll.nod")
nodFiles <- paste(region, suffices, sep = "")
triangles <- NULL
for (nodFile in nodFiles) {
if (1 == length(list.files(path = path, pattern = nodFile))) {
triangles <- read.table(paste(path, nodFile, sep = "/"), col.names = c("triangle", "longitude", "latitude"), encoding = "latin1")
oceDebug(debug, "found webtide information in \"", nodFile, "\"\n", sep = "")
break
} else {
oceDebug(debug, "looked for webtide information in \"", nodFile, "\" but this file does not exist\n", sep = "")
}
}
if (is.null(triangles)) {
stop("cannot find WebTide data file; rerun with debug=1 to see the searched list")
}
if (action == "map") {
if (plot) {
asp <- 1 / cos(rpd * mean(range(triangles$latitude, na.rm = TRUE)))
par(mfrow = c(1, 1), mar = c(3, 3, 2, 1), mgp = c(2, 0.7, 0))
plot(triangles$longitude, triangles$latitude,
pch = 2, cex = 1 / 4, lwd = 1 / 8,
asp = asp, xlab = "", ylab = "", ...
)
# Try for a coastline of well-suite resolution, if we have ocedata installed.
usr <- par("usr")
bestcoastline <- coastlineBest(lonRange = usr[1:2], latRange = usr[3:4], debug = debug - 1)
oceDebug(debug, "coastlineBest() suggests using", bestcoastline, "as the coastline\n")
if (bestcoastline == "coastlineWorld") {
data(list = bestcoastline, package = "oce", envir = environment())
coastlineWorld <- get("coastlineWorld")
} else {
if (requireNamespace("ocedata", quietly = TRUE)) {
data(list = bestcoastline, package = "ocedata", envir = environment())
oceDebug(debug, "Using", bestcoastline, "from the ocedata package.\n")
coastlineWorld <- get(bestcoastline)
} else {
data(list = "coastlineWorld", package = "oce", envir = environment())
oceDebug(debug, "The ocedata package is not available, so using", bestcoastline, "from oce\n")
coastlineWorld <- get("coastlineWorld")
}
}
polygon(coastlineWorld[["longitude"]], coastlineWorld[["latitude"]], col = "tan")
# use lon and lat, if node not given
if (!nodeGiven && longitudeGiven && latitudeGiven) {
closest <- which.min(geodDist(triangles$longitude, triangles$latitude, longitude, latitude))
node <- triangles$triangle[closest]
}
if (nodeGiven && node < 0 && interactive()) {
point <- locator(1)
node <- which.min(geodDist(triangles$longitude, triangles$latitude, point$x, point$y))
}
if (missing(node)) {
node <- triangles$number
longitude <- triangles$longitude
latitude <- triangles$latitude
} else {
if (is.finite(node)) {
node <- triangles$triangle[node]
longitude <- triangles$longitude[node]
latitude <- triangles$latitude[node]
points(longitude, latitude, pch = 20, cex = 2, col = "blue")
legend("topleft",
pch = 20, pt.cex = 2, cex = 3 / 4, col = "blue", bg = "white",
legend = sprintf("node %.0f %.3fN %.3fE", node, latitude, longitude)
)
}
}
oceDebug(debug, "} # webtide()\n", sep = "", unindent = 1, style = "bold")
return(invisible(list(node = node, latitude = latitude, longitude = longitude)))
} else {
node <- triangles$triangle
longitude <- triangles$longitude
latitude <- triangles$latitude
oceDebug(debug, "} # webtide()\n", sep = "", unindent = 1, style = "bold")
return(list(node = node, latitude = latitude, longitude = longitude))
}
} else if (action == "predict") {
if (missing(time)) {
time <- seq.POSIXt(from = presentTime(), by = "15 min", length.out = 7 * 4 * 24)
} # Q: what about timezone?
if (missing(node)) {
if (missing(longitude) || missing(latitude)) {
stop("'longitude' and 'latitude' must be given unless 'node' is given")
}
node <- which.min(geodDist(triangles$longitude, triangles$latitude, longitude, latitude))
} else {
latitude <- triangles$latitude[node]
longitude <- triangles$longitude[node]
}
oceDebug(debug, latitude, "N ", longitude, "E -- use node ", node,
" at ", triangles$latitude[node], "N ", triangles$longitude[node], "E\n",
sep = ""
)
constituentse <- dir(path = path, pattern = "*.s2c")
constituentsuv <- dir(path = path, pattern = "*.v2c")
nconstituents <- length(constituentse)
period <- ampe <- phasee <- ampu <- phaseu <- ampv <- phasev <- vector("numeric", length(nconstituents))
data("tidedata", package = "oce", envir = environment())
tidedata <- get("tidedata") # , pos=globalenv())
for (i in 1:nconstituents) {
twoLetter <- substr(constituentse[i], 1, 2)
C <- which(twoLetter == tidedata$const$name)
period[i] <- 1 / tidedata$const$freq[C]
# Elevation file contains one entry per node, starting with e.g.:
# tri
# period 23.934470 (hours) first harmonic
# 260.000000 (days)
# 1 0.191244 223.820954
# 2 0.188446 223.141200
coneFile <- paste(path, constituentse[i], sep = "/")
cone <- read.table(coneFile, skip = 3, encoding = "latin1")[node, ]
ampe[i] <- cone[[2]]
phasee[i] <- cone[[3]]
conuvFile <- paste(path, constituentsuv[i], sep = "/")
conuv <- read.table(conuvFile, skip = 3, encoding = "latin1")[node, ]
ampu[i] <- conuv[[2]]
phaseu[i] <- conuv[[3]]
ampv[i] <- conuv[[4]]
phasev[i] <- conuv[[5]]
oceDebug(debug, coneFile, sprintf("%s ", twoLetter),
sprintf("%4.2fh", period[i]),
sprintf(" (node %d) ", node),
sprintf("%4.4fm ", ampe[i]), sprintf("%3.3fdeg", phasee[i]), "\n",
sep = ""
)
}
elevation <- u <- v <- rep(0, length(time))
# NOTE: tref is the *central time* for tidem()
tRef <- ISOdate(1899, 12, 31, 12, 0, 0, tz = "UTC")
h <- (as.numeric(time) - as.numeric(tRef)) / 3600
tRef <- 3600 * round(mean(as.numeric(time)) / 3600)
for (i in 1:nconstituents) {
twoLetter <- substr(constituentse[i], 1, 2)
C <- which(twoLetter == tidedata$const$name)
vuf <- tidemVuf(tRef, j = C, latitude = latitude)
phaseOffset <- (vuf$u + vuf$v) * 360
# NOTE: phase is *subtracted* here, but *added* in tidem()
elevation <- elevation + ampe[i] * cos((360 * h / period[i] - phasee[i] + phaseOffset) * rpd)
#> lines(time, elevation, col=i,lwd=3) # Debug
u <- u + ampu[i] * cos((360 * h / period[i] - phaseu[i] + phaseOffset) * rpd)
v <- v + ampv[i] * cos((360 * h / period[i] - phasev[i] + phaseOffset) * rpd)
oceDebug(debug, sprintf("%s ", twoLetter),
sprintf("%4.2fh ", period[i]),
sprintf("%4.4fm ", ampe[i]), sprintf("%3.3fdeg", phasee[i]), "\n",
sep = ""
)
}
if (plot) {
par(mfrow = c(3, 1))
oce.plot.ts(time, elevation,
type = "l", xlab = "", ylab = resizableLabel("elevation"),
main = sprintf("node %.0f %.3fN %.3fE", node, latitude, longitude),
tformat = tformat
)
abline(h = 0, lty = "dotted", col = "gray")
oce.plot.ts(time, u,
type = "l", xlab = "", ylab = resizableLabel("u"),
drawTimeRange = FALSE, tformat = tformat
)
abline(h = 0, lty = "dotted", col = "gray")
oce.plot.ts(time, v,
type = "l", xlab = "", ylab = resizableLabel("v"),
drawTimeRange = FALSE, tformat = tformat
)
abline(h = 0, lty = "dotted", col = "gray")
oceDebug(debug, "} # webtide()\n", sep = "", unindent = 1, style = "bold")
return(invisible(list(
time = time, elevation = elevation, u = u, v = v,
node = node, basedir = basedir, region = region
)))
} else {
oceDebug(debug, "} # webtide()\n", sep = "", unindent = 1, style = "bold")
return(list(
time = time, elevation = elevation, u = u, v = v,
node = node, basedir = basedir, region = region
))
}
}
}
dankelley/oce documentation built on April 28, 2024, 1:47 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions webtide predict.tidem tidem tidemConstituentNameFix tidemAstron tidemVuf as.tidem
Documented in as.tidem predict.tidem tidem tidemAstron tidemConstituentNameFix tidemVuf webtide
# vim:textwidth=128:expandtab:shiftwidth=4:softtabstop=4
#' Class to Store Tidal Models
#'
#' This class stores tidal-constituent models.
#'
#' @templateVar class tidem
#'
#' @templateVar dataExample {}
#'
#' @templateVar metadataExample {}
#'
#' @template slot_summary
#'
#' @template slot_put
#'
#' @template slot_get
#'
#' @author Dan Kelley
#' @family functions that plot oce data
#' @family things related to tides
setClass("tidem", contains = "oce")
setMethod(
f = "initialize",
signature = "tidem",
definition = function(.Object, ...) {
.Object <- callNextMethod(.Object, ...)
.Object@metadata$version <- ""
.Object@processingLog$time <- presentTime()
.Object@processingLog$value <- "create 'tidem' object"
return(.Object)
}
)
#' Tidal Constituent Information
#'
#' The `tidedata` dataset contains Tide-constituent information that is
#' use by [tidem()] to fit tidal models. `tidedata` is a list
#' containing
#' \describe{
#' \item{`const`}{
#' a list containing vectors
#' `name` (a string with constituent name),
#' `freq` (the frequency, in cycles per hour),
#' `kmpr` (a string naming the comparison constituent, blank if there is none),
#' `ikmpr` (index of comparison constituent, or `0` if there is none),
#' `df` (frequency difference between constituent and its
#' comparison, used in the Rayleigh criterion),
#' `d1` through `d6` (the first through sixth Doodson numbers),
#' `semi`,
#' `nsat` (number of satellite constituents),
#' `ishallow`,
#' `nshallow`,
#' `doodsonamp`,
#' and
#' `doodsonspecies`.
#' }
#' \item{`sat`}{
#' a list containing vectors
#' `deldood`,
#' `phcorr`,
#' `amprat`,
#' `ilatfac`,
#' and
#' `iconst`.
#' }
#' \item{`shallow`}{
#' a list containing vectors
#' `iconst`,
#' `coef`,
#' and
#' `iname`.
#' }
#' }
#' Apart from the use of `d1` through `d6`, the naming and content
#' follows `T_TIDE` (see Pawlowicz et al. 2002), which in turn builds upon
#' the analysis of Foreman (1978).
#'
#' @name tidedata
#'
#' @docType data
#'
#' @author Dan Kelley
#'
#' @references
#'
#' Foreman, M. G. G., 1978. Manual for Tidal Currents Analysis and Prediction.
#' Pacific Marine Science Report. British Columbia, Canada: Institute of Ocean
#' Sciences, Patricia Bay.
#'
#' Pawlowicz, Rich, Bob Beardsley, and Steve Lentz, 2002. Classical tidal
#' harmonic analysis including error estimates in MATLAB using `T_TIDE`.
#' Computers and Geosciences, 28, 929-937.
#'
#' @source The data come from the `tide3.dat` file of the `T_TIDE`
#' package (Pawlowicz et al., 2002), and derive from Appendices provided by
#' Foreman (1978). The data are scanned using \file{tests/tide.R} in this
#' package, which also performs some tests using `T_TIDE` values as a
#' reference.
#'
#' @family things related to tides
NULL
#' Tidal Current Dataset
#'
#' The `tidalCurrent` dataset contains tidal velocities reported in
#' Foreman's (1978) report (reference 1) on his Fortran code for the analysis of
#' tidal currents and provided in an associated webpage (reference 2).
#' Here, `tidalCurrent` is data frame containing
#' * `time` a POSIXct time.
#' * `u` the eastward component of velocity in m/s.
#' * `v` the northward component of velocity in m/s.
#'
#' @name tidalCurrent
#'
#' @docType data
#'
#' @examples
#' library(oce)
#' data(tidalCurrent)
#' par(mfrow = c(2, 1))
#' oce.plot.ts(tidalCurrent$time, tidalCurrent$u, ylab = "u [m/s]")
#' abline(h = 0, col = 2)
#' oce.plot.ts(tidalCurrent$time, tidalCurrent$v, ylab = "v [m/s]")
#' abline(h = 0, col = 2)
#'
#' @author Dan Kelley (reformatting data provided by Michael Foreman)
#'
#' @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.
#' 2. \code{https://www.dfo-mpo.gc.ca/science/documents/data-donnees/tidal-marees/tidpack.zip}
#'
#' @source The data come from the `tide8.dat` and `tide9.dat` files provided
#' at reference 2.
#'
#' @family things related to tides
NULL
#' Summarize a tidem Object
#'
#' By default, all fitted constituents are plotted, but it is quite useful to
#' set e.g. p=0.05 To see just those constituents that are significant at the 5
#' percent level.
#' Note that the p values are estimated as the average of the p values for the
#' sine and cosine components at a given frequency.
#'
#' @param object an object of class [tidem], as created by
#' [as.tidem()] or [tidem()].
#'
#' @param p optional value of the maximum p value for the display of an
#' individual coefficient. If not given, all coefficients are shown.
#'
#' @param constituent optional character vector holding the names
#' of constituents on which to focus.
#' @template tideconst
#'
#' @param \dots further arguments passed to or from other methods.
#'
#' @return `NULL`
#'
#' @author Dan Kelley
#'
#' @section Sample of Usage:
#' \preformatted{
#' library(oce)
#' data(sealevel)
#' tide <- tidem(sealevel)
#' summary(tide)
#' }
#'
#' @family things related to tides
setMethod(
f = "summary",
signature = "tidem",
definition = function(object, p = 1.0, constituent, ...) {
debug <- if ("debug" %in% names(list(...))) list(...)$debug else 0
version <- object@metadata$version
ok <- object@data$p <= p | version == 3L
haveP <- any(!is.na(object@data$p))
if (missing(constituent)) {
fit <- data.frame(
Const = object@data$const[ok],
Name = object@data$name[ok],
Freq = object@data$freq[ok],
Amplitude = object@data$amplitude[ok],
Phase = object@data$phase[ok],
p = object@data$p[ok]
)
if (debug) {
cat("For missing(constituent) case, fit is:\n")
print(fit)
}
} else {
i <- NULL
bad <- NULL
for (iconst in seq_along(constituent)) {
w <- which(object@data$name == constituent[iconst])
if (length(w) == 1L) {
i <- c(i, w)
} else {
bad <- c(bad, iconst)
}
}
if (length(bad)) {
warning("the following constituents are not handled: \"",
paste(constituent[bad], collapse = "\", \""), "\"\n",
sep = ""
)
}
if (length(i) == 0) {
stop("In summary,tidem-method() : no known constituents were provided", call. = FALSE)
}
i <- unique(i)
fit <- data.frame(
Const = object@data$const[i],
Name = object@data$name[i],
Freq = object@data$freq[i],
Amplitude = object@data$amplitude[i],
Phase = object@data$phase[i],
p = object@data$p[i]
)
if (debug) {
cat("For !missing(constituent) case, fit is:\n")
print(fit)
}
}
cat("tidem summary\n-------------\n")
if (version != "3") {
cat("\nCall:\n")
cat(paste(deparse(object[["call"]]), sep = "\n", collapse = "\n"), "\n", sep = "")
cat("RMS misfit to data: ", sqrt(var(object[["model"]]$residuals)), "\n")
cat("\nFitted Model:\n")
f <- fit[3:6]
if (debug > 0) {
cat("fit:\n")
print(fit)
cat("f:\n")
print(f)
}
rownames(f) <- as.character(fit[, 2])
if (haveP) {
printCoefmat(f,
digits = 3,
signif.stars = getOption("show.signif.stars"),
signif.legend = TRUE,
P.values = TRUE, has.Pvalue = TRUE, ...
)
} else {
printCoefmat(f[, -4], digits = 3)
}
} else {
cat("\nSupplied Model:\n")
f <- fit[3:5]
rownames(f) <- as.character(fit[, 2])
printCoefmat(f, digits = 3)
}
processingLogShow(object)
invisible(NULL)
}
)
#' Extract Something From a tidem Object
#'
#' @param x a [tidem-class] object.
#'
#' @section Details of the Specialized Method:
#'
#' * If `i` is `"?"`, then the return value is a list
#' containing four items, each of which is a character vector
#' holding the names of things that can be accessed with `[[`.
#' The `data` and `metadata` items hold the names of
#' entries in the object's data and metadata
#' slots, respectively. Note that `metadataDerived`
#' holds only `""`, because no derived metadata values
#' are defined for `tidem` objects.
#'
#' * If `i` is `"frequency"` or `"freq"`, then a vector of
#' constituent frequencies is returned.
#'
#' * If `i` is `"amplitude"` then a vector of constituent amplitudes
#' is returned.
#'
#' * If `i` is `"phase"` then a vector of constituent phases
#' is returned.
#'
#' * If `i` is `"constituents"` then a data frame holding constituent
#' name, frequency, amplitude and phase is returned.
#'
#' * If `i` is a vector of constituent names, then the return
#' value is as for `"constituents"`, except that only the named
#' those constituents are returned.
#'
#' @template sub_subTemplate
#'
#' @family things related to tides
setMethod(
f = "[[",
signature(x = "tidem", i = "ANY", j = "ANY"),
definition = function(x, i, j, ...) {
# UNUSED const <- x@data$name
if (i[1] == "?") {
return(list(
metadata = sort(names(x@metadata)),
metadataDerived = NULL,
data = c(sort(names(x@data)), x@data$name),
dataDerived = "frequency"
))
}
if (i[1] == "frequency") {
return(x@data$freq)
} else if (i[1] == "constituents" || sum(i %in% x@data$name) > 0L) {
all <- with(
x@data,
data.frame(
name = name, freq = freq,
amplitude = amplitude, phase = phase
)
)
if (i[[1]] == "constituents") {
return(all)
} else {
return(all[all$name %in% i, ])
}
}
callNextMethod() # [[
}
)
#' Replace Parts of a tidem Object
#'
#' @param x a [tidem-class] object.
#'
#' @template sub_subsetTemplate
#'
#' @family things related to tides
setMethod(
f = "[[<-",
signature(x = "tidem", i = "ANY", j = "ANY"),
definition = function(x, i, j, ..., value) {
callNextMethod(x = x, i = i, j = j, ... = ..., value = value) # [[<-
}
)
#' Plot a tidem Object
#'
#' Plot a summary diagram for a tidal fit.
#'
#' @param x a [tidem-class] object.
#'
#' @param which integer flag indicating plot type, 1 for stair-case spectral, 2
#' for spike spectral.
#'
#' @param constituents character vector holding the names of constituents that are
#' to be drawn and labelled. If `NULL`, then no constituents will be shown.
#' @template tideconst
#'
#' @param sides an integer vector of length equal to that of `constituents`,
#' designating the side on which the constituent labels are to be drawn. As in
#' all R graphics, the value `1` indicates the bottom of the plot, and
#' `3` indicates the top. If `sides=NULL`, the default, then all labels
#' are drawn at the top. Any value of `sides` that is not either 1 or 3
#' is converted to 3.
#'
#' @param col a character vector naming colors to be used for `constituents`.
#' Ignored if `sides=3`. Repeated to be of the same length as
#' `constituents`, otherwise.
#'
#' @param log if set to "`x`", the frequency axis will be logarithmic.
#'
#' @param mgp 3-element numerical vector to use for `par(mgp)`, and also
#' for `par(mar)`, computed from this. The default is tighter than the R
#' default, in order to use more space for the data and less for the axes.
#'
#' @param mar value to be used with `[`par`]("mar")`.
#'
#' @param \dots optional arguments passed to plotting functions, not all
#' of which are obeyed. For example, if \dots contains `type`, that value will be
#' ignored because it is set internally, according to the value of `which`.
#'
#' @section Sample of Usage:
#' \preformatted{
#' library(oce)
#' data(sealevel)
#' tide <- tidem(sealevel)
#' plot(tide)
#' }
#'
#' @section Historical note:
#' An argument named `labelIf` was removed in July 2016,
#' because it was discovered never to have worked as documented, and
#' because the more useful argument `constituents` had been added.
#'
#' @author Dan Kelley
#'
#' @family functions that plot oce data
#'
#' @aliases plot.tidem
#'
#' @family things related to tides
setMethod(
f = "plot",
signature = signature("tidem"),
definition = function(x,
which = 1,
constituents = c("SA", "O1", "K1", "M2", "S2", "M4"),
sides = NULL,
col = "blue",
log = "",
mgp = getOption("oceMgp"),
mar = c(mgp[1] + 1, mgp[1] + 1, mgp[2] + 0.25, mgp[2] + 1),
...) {
dots <- list(...)
dotsNames <- names(dots)
data("tidedata", package = "oce", envir = environment())
tidedata <- get("tidedata") # , pos=globalenv())
drawConstituent <- function(name = "M2", side = 3, col = "blue", adj = NULL, cex = 0.8) {
w <- which(tidedata$const$name == name)
if (!length(w)) {
warning("constituent \"", name, "\" is unknown")
return()
}
frequency <- tidedata$const$freq[w]
abline(v = frequency, col = col, lty = "dotted")
if (par("usr")[1] < frequency && frequency <= par("usr")[2]) {
if (is.null(adj)) {
mtext(name, side = side, at = frequency, col = col, cex = cex)
} else {
mtext(name, side = side, at = frequency, col = col, cex = cex, adj = adj)
}
}
}
opar <- par(no.readonly = TRUE)
lw <- length(which)
if (lw > 1) {
on.exit(par(opar))
}
par(mgp = mgp, mar = mar)
frequency <- x@data$freq[-1] # trim z0
amplitude <- x@data$amplitude[-1]
# name <- x@data$name[-1]
# Place in frequency order (required for cumulative plot), because user
# might have given constituents in any order.
order <- order(frequency)
frequency <- frequency[order]
# name <- name[order]
amplitude <- amplitude[order]
# print(data.frame(name=name, period=1/frequency, amplitude=amplitude))
if (!is.null(constituents)) {
sides <- if (is.null(sides)) {
rep(3, length(constituents))
} else {
rep(sides, length.out = length(constituents))
}
col <- rep(col, length.out = length(constituents))
}
sides[sides != 1 & sides != 3] <- 3 # default to top
cex <- if ("cex" %in% names(dots)) dots$cex else 0.8
# We specify these things directly but they are not parameters of this
# function, so we must remove them from '...' before passing that to plot().
if ("cex" %in% dotsNames) dots$cex <- NULL
if ("type" %in% dotsNames) dots$type <- NULL
for (w in 1:lw) {
# message("w=", w, "; which[w]=", which[w])
if (which[w] == 2) {
plot(frequency, amplitude,
type = "n",
xlab = "Frequency [ cph ]", ylab = "Amplitude [ m ]", log = log, ...
)
segments(frequency, 0, frequency, amplitude)
for (i in seq_along(constituents)) {
drawConstituent(constituents[i], side = sides[i], col = col[i], cex = cex)
}
} else if (which[w] == 1) {
plot(frequency, cumsum(amplitude),
xlab = "Frequency [ cph ]", ylab = "Amplitude [ m ]", log = log, type = "s", ...
)
for (i in seq_along(constituents)) {
drawConstituent(constituents[i], side = sides[i], col = col[i], cex = cex)
}
} else {
stop("unknown value of which ", which, "; should be 1 or 2")
}
}
if (!all(is.na(pmatch(names(list(...)), "main")))) title(...)
}
)
#' Create tidem Object From Fitted Harmonic Data
#'
#' This function takes a set of tidal constituent amplitudes
#' and phases, and constructs a return value of similar form
#' to that returned by [tidem()]. Its purpose is to enable
#' predictions based on published constituent amplitudes
#' and phases. Since [as.tidem()] does not account for a
#' reference height, it is the user's responsible to account
#' for this after a prediction is made using [predict.tidem()].
#'
#' All the constituent names used by [tidem()] are permitted here,
#' *except* for `"Z0"` (see \dQuote{Description} regarding reference
#' height).
#' To get a list of constituent names, please consult Foreman (1978),
#' or type the following in an R console:
#' \preformatted{
#' data(tidedata)
#' data.frame(name=tidedata$const$name, freq=tidedata$const$freq)
#' }
#'
#' In addition to the above, [as.tidem()] can handle NOAA names
#' of constituents. For the most part, these match oce names, but
#' there are 4 exceptions: NOAA names
#' `"LAM2"`, `"M1"`, `"RHO"`, and `"2MK3"` are converted to oce names
#' `"LDA2"`, `"NO1"`, `"RHO1"`, and `"MO3"`. The name mapping was
#' inferred by matching frequencies; for these constituents, the
#' fractional mismatch in frequencies was under 4e-8;
#' see Reference 5 for more details.
#' A message is printed if these name conversions are required
#' in the particular use of [as.tidem()].
#'
#' Apart from the standard oce names and this set of NOAA synonyms,
#' any other constituent name is reported in a warning message.
#'
#' @param tRef a POSIXt value indicating the mean time of the
#' observations used to develop the harmonic model. This is rounded
#' to the nearest hour in [as.tidem()], to match the behaviour
#' of [tidem()].
#'
#' @param latitude numerical value indicating the latitude of the
#' observations that were used to create the harmonic model. This
#' is needed for nodal-correction procedures carried out
#' by [tidemVuf()].
#'
#' @param name character vector holding names of constituents.
#'
#' @param amplitude,phase numeric vectors of constituent amplitudes
#' and phases. These must be of the same length as `name`.
#'
#' @param frequency,speed optional numeric vectors giving the frequencies
#' of the constituents (in cycles per hour) or the analogous speeds
#' (in degrees per hour). Only one of these may be given, and a conversion
#' is done from the latter to the former, if required. If the frequencies
#' are thus specified, then these are used instead of the frequencies that
#' oce normally used, as defined in `data(tideconst)`. A warning will
#' be issued if the absolute value of the relative frequency mismatch for any
#' given component exceeds 1e-6, and this will occur for any NOAA tables
#' containing the SA component, for which this relative mismatch
#' is approximately 4e-5 (see reference 5).
#'
#' @template debugTemplate
#'
#' @return An object of [tidem-class], with only minimal
#' contents.
#'
#' @section Known issues:
#' There are two known differences between [tidem()] and the Matlab
#' `T_TIDE` package, as listed in references 3 and 4.
#'
#' @examples
#'
#' # Example 1: show agreement with tidem()
#' data(sealevelTuktoyaktuk)
#' # 'm0' is model fitted by tidem()
#' m0 <- tidem(sealevelTuktoyaktuk)
#' p0 <- predict(m0, sealevelTuktoyaktuk[["time"]])
#' m1 <- as.tidem(
#' mean(sealevelTuktoyaktuk[["time"]]), sealevelTuktoyaktuk[["latitude"]],
#' m0[["name"]], m0[["amplitude"]], m0[["phase"]]
#' )
#' # Test agreement with tidem() result, by comparing predicted sealevels.
#' p1 <- predict(m1, sealevelTuktoyaktuk[["time"]])
#' stopifnot(max(abs(p1 - p0), na.rm = TRUE) < 1e-10)
#'
#' # Example 2: See the effect of dropping weak constituents
#' m0[["name"]][which(m0[["amplitude"]] > 0.05)]
#' h <- "
#' name amplitude phase
#' Z0 1.98061875 0.000000
#' MM 0.21213065 263.344739
#' MSF 0.15605629 133.795004
#' O1 0.07641438 74.233130
#' K1 0.13473817 81.093134
#' OO1 0.05309911 235.749693
#' N2 0.08377108 44.521462
#' M2 0.49041340 77.703594
#' S2 0.22023705 137.475767"
#' coef <- read.table(text = h, header = TRUE)
#' m2 <- as.tidem(
#' mean(sealevelTuktoyaktuk[["time"]]),
#' sealevelTuktoyaktuk[["latitude"]],
#' coef$name, coef$amplitude, coef$phase
#' )
#' p2 <- predict(m2, sealevelTuktoyaktuk[["time"]])
#' par(mfrow = c(3, 1))
#' oce.plot.ts(sealevelTuktoyaktuk[["time"]], p0)
#' ylim <- par("usr")[3:4] # to match scales in other panels
#' oce.plot.ts(sealevelTuktoyaktuk[["time"]], p1, ylim = ylim)
#' oce.plot.ts(sealevelTuktoyaktuk[["time"]], p2, ylim = ylim)
#'
#' @references
#'
#' 1. Foreman, M. G. G., 1978. Manual for Tidal Currents Analysis and Prediction.
#' Pacific Marine Science Report. British Columbia, Canada: Institute of Ocean
#' Sciences, Patricia Bay.
#'
#' 2. Wikipedia, "Theory of Tides." https://en.wikipedia.org/wiki/Theory_of_tides
#' Downloaded Aug 17, 2019.
#'
#' 3. Github issue 1653 "tidem() and t_tide do not produce identical results"
#' (https://github.com/dankelley/oce/issues/1653)
#'
#' 4. Github issue 1654 "predict(tidem()) uses all constituents, unlike T_TIDE"
#' (https://github.com/dankelley/oce/issues/1654)
#'
#' 5. Github issue 2143 "mismatch in oce/NOAA frequency of SA tidal constituent"
#' (https://github.com/dankelley/oce/issues/2143)
#'
#' @family things related to tides
as.tidem <- function(tRef, latitude, name, amplitude, phase, frequency, speed, debug = getOption("oceDebug")) {
oceDebug(debug, "as.tidem() {\n", sep = "", unindent = 1)
if (missing(tRef)) {
stop("tRef must be given")
}
if (missing(latitude)) {
stop("latitude must be given")
}
if (missing(name)) {
stop("name must be given")
}
if (missing(amplitude)) {
stop("amplitude must be given")
}
if (missing(phase)) {
stop("phase must be given")
}
gaveFrequency <- !missing(frequency)
gaveSpeed <- !missing(speed)
if (gaveSpeed) {
if (gaveFrequency) {
stop("cannot give both speed and frequency")
}
frequency <- speed / 360.0
gaveFrequency <- TRUE
}
nname <- length(name)
if (nname != length(amplitude)) {
stop("lengths of name and amplitude must be equal but they are ", nname, " and ", length(amplitude))
}
if (nname != length(phase)) {
stop("lengths of name and phase must be equal but they are ", nname, " and ", length(phase))
}
# convert from NOAA name to oce name
# NOAA: LAM2 M1 RHO 2MK3
# oce: LDA2 NO1 RHO1 MO3
tmpNOAA <- c("LAM2", "M1", "RHO", "2MK3")
tmpOce <- c("LDA2", "NO1", "RHO1", "MO3")
for (tmpi in seq_along(tmpNOAA)) {
if (tmpNOAA[tmpi] %in% name) {
name[name == tmpNOAA[tmpi]] <- tmpOce[tmpi]
message("converted NOAA name \"", tmpNOAA[tmpi], "\" to oce name \"", tmpOce[tmpi], "\"")
}
}
data("tidedata", package = "oce", envir = environment())
tidedata <- get("tidedata")
tRef <- numberAsPOSIXct(3600 * round(as.numeric(tRef, tz = "UTC") / 3600), tz = "UTC")
oceDebug(debug, "input head(name): ", paste(head(name), collapse = " "), "\n")
oceDebug(debug, "input head(phase): ", paste(head(phase), collapse = " "), "\n")
oceDebug(debug, "input head(amplitude): ", paste(head(amplitude), collapse = " "), "\n")
freq <- rep(NA, nname)
indices <- rep(NA, nname)
ibad <- NULL
if (gaveSpeed) {
oceDebug(debug, "computing frequency from velocity\n")
frequency <- speed / 360.0
gaveFrequency <- TRUE
}
for (i in seq_along(name)) {
oceDebug(debug, "adjusting amplitude and phase for constituent \"", name[i], "\"\n", sep = "")
j <- which(tidedata$const$name == name[i])
oceDebug(debug, " inferred j=", j, " from constituent name\n", sep = "")
if (length(j)) {
vuf <- tidemVuf(tRef, j = j, latitude = latitude)
oceDebug(debug, " inferred vuf=", deparse(vuf), "\n")
indices[i] <- j
amplitude[i] <- amplitude[i] * vuf$f
phase[i] <- phase[i] - (vuf$v + vuf$u) * 360
if (gaveFrequency) {
freq[i] <- frequency[i]
afd <- abs(frequency[i] - tidedata$const$freq[j]) / tidedata$const$freq[j]
if (afd > 1e-6) {
msg <- sprintf("%s: absolute fractional mismatch between specified and built-in frequency is %.3e\n", name[i], afd)
warning(msg)
}
oceDebug(debug, msg)
} else {
freq[i] <- tidedata$const$freq[j]
}
} else {
ibad <- c(ibad, i)
}
}
if (length(ibad)) {
warning("the following constituents are not handled: \"",
paste(name[ibad], collapse = "\", \""), "\"\n",
sep = ""
)
indices <- indices[-ibad]
name <- name[-ibad]
amplitude <- amplitude[-ibad]
phase <- phase[-ibad]
freq <- freq[-ibad]
}
oceDebug(debug, "after vuf correction, head(name): ", paste(head(name), collapse = " "), "\n")
oceDebug(debug, "after vuf correction, head(phase): ", paste(head(phase), collapse = " "), "\n")
oceDebug(debug, "after vuf correction, head(amplitude): ", paste(head(amplitude), collapse = " "), "\n")
oceDebug(debug, "} # as.tidem()\n", sep = "", unindent = 1)
phase <- phase %% 360
res <- new("tidem")
res@data <- list(
tRef = tRef,
const = indices,
name = name,
freq = freq,
amplitude = amplitude,
phase = phase,
p = rep(NA, length(name))
)
res@metadata$version <- 3
res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep = "", collapse = ""))
oceDebug(debug, "} # as.tidem()\n", sep = "", unindent = 1)
res
}
#' Nodal Modulation Calculations for Tidal Analyses
#'
#' Carry out nodal modulation calculations for [tidem()]. This function is based directly
#' on `t_vuf` in the `T_TIDE` Matlab package (Pawlowicz et al., 2002),
#' which inherits from the Fortran code described by Foreman (1978).
#'
#' @param t a single time in [POSIXct()] format, with timezone `"UTC"`.
#'
#' @param j integer vector, giving indices of tidal constituents to use.
#'
#' @param latitude optional numerical value containing the latitude in degrees North.
#' If not provided, `u` in the return value will be a vector consisting of
#' repeated 0 value, and `f` will be a vector of repeated 1 value.
#'
#' @return A `list` containing
#' items named `v`, `u` and `f` as described in the `T_TIDE` documentation,
#' as well in Pawlowicz et al. (2002) and Foreman (1978).
#'
#' @author Dan Kelley translated this from the `t_vuf` function
#' of the `T_TIDE` Matlab package (see Pawlowicz et al. 2002).
#'
#' @examples
#' # Look up values for the M2 constituent in Halifax Harbour, Canada.
#' library(oce)
#' data("tidedata")
#' j <- with(tidedata$const, which(name == "M2"))
#' tidemVuf(t = as.POSIXct("2008-01-22 18:50:24"), j = j, lat = 44.63)
#'
#' @references
#' * Foreman, M. G. G., 1978. Manual for Tidal Currents Analysis and Prediction.
#' Pacific Marine Science Report. British Columbia, Canada: Institute of Ocean
#' Sciences, Patricia Bay.
#'
#' * Pawlowicz, Rich, Bob Beardsley, and Steve Lentz, 2002. Classical tidal
#' harmonic analysis including error estimates in MATLAB using `T_TIDE`.
#' Computers and Geosciences, 28, 929-937.
#'
#' @family things related to tides
tidemVuf <- function(t, j, latitude = NULL) {
debug <- 0
if (length(t) > 1) {
stop("t must be a single POSIXct item")
}
data("tidedata", package = "oce", envir = environment())
tidedata <- get("tidedata") # , pos=globalenv())
a <- tidemAstron(t)
if (debug > 0) print(a)
doodson <- cbind(
tidedata$const$d1,
tidedata$const$d2,
tidedata$const$d3,
tidedata$const$d4,
tidedata$const$d5,
tidedata$const$d6
)
oceDebug(
debug,
"doodson[1,]=", doodson[1, ], "\n",
"doodson[2,]=", doodson[2, ], "\n",
"doodson[3,]=", doodson[3, ], "\n"
)
v <- doodson %*% a$astro + tidedata$const$semi
oceDebug(debug, "tidedata$const$semi[", j, "]=", tidedata$const$semi[j], "\n")
v <- v - trunc(v)
oceDebug(debug, "v[1:3]=", v[1:3], "\n")
if (!is.null(latitude) && !is.na(latitude)) {
if (abs(latitude) < 5) {
latitude <- sign(latitude) * 5
}
slat <- sin(pi * latitude / 180)
k <- which(tidedata$sat$ilatfac == 1)
rr <- tidedata$sat$amprat
rr[k] <- rr[k] * 0.36309 * (1.0 - 5.0 * slat * slat) / slat
k <- which(tidedata$sat$ilatfac == 2)
rr[k] <- rr[k] * 2.59808 * slat
uu <- tidedata$sat$deldood %*% a$astro[4:6] + tidedata$sat$phcorr
uu <- uu - trunc(uu)
oceDebug(debug, "uu[1:3]=", uu[1:3], "\n")
nsat <- length(tidedata$sat$iconst)
# nfreq <- length(tidedata$const$numsat)
# loop, rather than make a big matrix
oceDebug(
debug,
"tidedata$sat$iconst=", tidedata$sat$iconst, "\n",
"length(sat$iconst)=", length(tidedata$sat$iconst), "\n"
)
fsum.vec <- vector("numeric", nsat)
u.vec <- vector("numeric", nsat)
for (isat in 1:nsat) {
oceDebug(debug, "isat=", isat, "\n")
use <- tidedata$sat$iconst == isat
fsum.vec[isat] <- 1 + sum(rr[use] * exp(1i * 2 * pi * uu[use]))
u.vec[isat] <- Arg(fsum.vec[isat]) / 2 / pi
if (isat == 8 && debug > 0) {
cat("TEST at isat=8:\n")
cat("fsum.vec[", isat, "]=", fsum.vec[isat], " (EXPECT 1.18531604917590 - 0.08028013402313i)\n")
cat("u.vec[ ", isat, "]=", u.vec[isat], " (EXPECT -0.01076294959868)\n")
}
}
oceDebug(
debug,
"uvec[", j, "]=", u.vec[j], "\n",
"fsum.vec[", j, "]=", fsum.vec[j], "\n"
)
f <- abs(fsum.vec)
u <- Arg(fsum.vec) / 2 / pi
oceDebug(debug, "f=", f, "\n")
oceDebug(debug, "u=", u, "\n")
for (k in which(!is.na(tidedata$const$ishallow))) {
ik <- tidedata$const$ishallow[k] + 0:(tidedata$const$nshallow[k] - 1)
f[k] <- prod(f[tidedata$shallow$iname[ik]]^abs(tidedata$shallow$coef[ik]))
u[k] <- sum(u[tidedata$shallow$iname[ik]] * tidedata$shallow$coef[ik])
v[k] <- sum(v[tidedata$shallow$iname[ik]] * tidedata$shallow$coef[ik])
if (debug > 0 && k < 28) {
cat("k=", k, "f[k]=", f[k], " u[k]=", u[k], "v[k]=", v[k], "\n")
}
}
u <- u[j]
v <- v[j]
f <- f[j]
} else {
v <- v[j]
u <- rep(0, length(j))
f <- rep(1, length(j))
}
if (length(v) < length(u)) {
warning("trimming u and f to get same length as v -- this is a bug")
u <- u[seq_along(v)]
f <- f[seq_along(v)]
}
list(v = v, u = u, f = f)
}
#' Astronomical Calculations for tidem
#'
#' Do some astronomical calculations for [tidem()]. This function is based directly
#' on `t_astron` in the `T_TIDE` Matlab package (see Pawlowicz et al. 2002),
#' which inherits from the Fortran code described by Foreman (1978).
#'
#' @param t Either a time in `POSIXct` format (with `"UTC"` timezone,
#' or else odd behaviours may result),
#' or an integer. In the second case, it is converted to a time with
#' [numberAsPOSIXct()], using `tz="UTC"`.
#'
#' @return A `list` containing items named
#' `astro` and `ader` (see the `T_TIDE` documentation).
#'
#' @author Dan Kelley translated this from the `t_astron` function
#' of the `T_TIDE` Matlab package (see Pawlowicz et al. 2002).
#'
#' @examples
#' tidemAstron(as.POSIXct("2008-01-22 18:50:24"))
#'
#' @references
#' * Foreman, M. G. G., 1978. Manual for Tidal Currents Analysis and Prediction.
#' Pacific Marine Science Report. British Columbia, Canada: Institute of Ocean
#' Sciences, Patricia Bay.
#'
#' * Pawlowicz, Rich, Bob Beardsley, and Steve Lentz, 2002. Classical tidal
#' harmonic analysis including error estimates in MATLAB using `T_TIDE`.
#' Computers and Geosciences, 28, 929-937.
#'
#' @family things related to tides
tidemAstron <- function(t) {
if (length(t) > 1) {
stop("t must be a single POSIXct item")
}
debug <- FALSE
if (is.numeric(t)) {
t <- numberAsPOSIXct(t, tz = "UTC")
}
d <- as.numeric(difftime(t, ISOdatetime(1899, 12, 31, 12, 0, 0, tz = "UTC"), units = "days"))
D <- d / 10000
a <- matrix(c(1.0, d, D^2, D^3), ncol = 1)
oceDebug(debug, "d=", formatC(d, digits = 10), "D=", D, "a=", a, "\n")
scHcPcNpPp <-
matrix(
c(
270.434164, 13.1763965268, -0.0000850, 0.000000039,
279.696678, 0.9856473354, 0.00002267, 0.000000000,
334.329556, 0.1114040803, -0.0007739, -0.00000026,
-259.183275, 0.0529539222, -0.0001557, -0.000000050,
281.220844, 0.0000470684, 0.0000339, 0.000000070
),
nrow = 5, ncol = 4, byrow = TRUE
)
astro <- ((scHcPcNpPp %*% a) / 360) %% 1
oceDebug(debug, "astro=", astro, "\n")
rem <- as.numeric(difftime(t, trunc.POSIXt(t, units = "days"), tz = "UTC", units = "days"))
oceDebug(debug, "rem2=", rem, "\n")
tau <- rem + astro[2, 1] - astro[1, 1]
astro <- c(tau, astro)
da <- matrix(c(0.0, 1.0, 2e-4 * D, 3e-4 * D^2), nrow = 4, ncol = 1)
ader <- (scHcPcNpPp %*% da) / 360
dtau <- 1 + ader[2, 1] - ader[1, 1]
ader <- c(dtau, ader)
list(astro = astro, ader = ader)
}
#' Change Tidal Constituent Name from T-TIDE to Foreman Convention
#'
#' This is used by [tidem()] to permit users to specify constituent names in either
#' the T-TIDE convention (see Pawlowicz et al. 2002) or Foreman convention
#' (see Foreman (1978). There are only two such instances:
#' `"MS"`, which gets translated to `"M8"`, and `"UPSI"`,
#' which gets translated to `"UPS1"`.
#'
#' @param names a vector of character values, holding constituent names
#'
#' @param debug an integer controlling warnings. If this is zero, then no warnings
#' are issued during processing; otherwise, as is the default, warnings are
#' issued for each conversion that is required.
#'
#' @return A vector of character values of tidal constituent names, in the Foreman naming convention.
#'
#' @references
#' Foreman, M. G. G., 1978. Manual for Tidal Currents Analysis and Prediction.
#' Pacific Marine Science Report. British Columbia, Canada: Institute of Ocean
#' Sciences, Patricia Bay.
#'
#' Pawlowicz, Rich, Bob Beardsley, and Steve Lentz, 2002. Classical tidal
#' harmonic analysis including error estimates in MATLAB using `T_TIDE`.
#' Computers and Geosciences, 28, 929-937.
tidemConstituentNameFix <- function(names, debug = 1) {
if ("MS" %in% names) {
if (debug) {
warning("constituent name switched from T-TIDE 'MS' to Foreman (1978) 'M8'")
}
names[names == "MS"] <- "M8"
}
if ("-MS" %in% names) {
if (debug) {
warning("removed-constituent name switched from T-TIDE 'MS' to Foreman (1978) 'M8'")
}
names[names == "-MS"] <- "-M8"
}
if ("UPSI" %in% names) {
if (debug) {
warning("constituent name switched from T-TIDE 'UPSI' to Foreman (1978) 'UPS1'")
}
names[names == "UPSI"] <- "UPS1"
}
if ("-UPSI" %in% names) {
if (debug) {
warning("removed-constituent name switched from T-TIDE 'UPSI' to Foreman (1978) 'UPS1'")
}
names[names == "-UPSI"] <- "-UPS1"
}
names
}
#' Fit a Tidal Model to a Timeseries
#'
#' The fit is done in terms of sine and cosine components at the indicated
#' tidal frequencies (after possibly pruning to satisfy the Rayleigh criterion,
#' as explained in phase 2 of the procedure outlined in
#' \dQuote{Details}), with the amplitude and phase being
#' calculated from the resultant coefficients on the sine and cosine terms.
#' The scheme was devised for hourly data; for other sampling schemes,
#' see \dQuote{Application to non-hourly data}.
#'
#' A summary of constituents used by [tidem()] may be found with:
#' \preformatted{
#' data(tidedata)
#' print(tidedata$const)
#' }
#'
#' A multi-stage procedure is used to establish the list of tidal
#' constituents to be used in the fit.
#'
#' *Phase 1: initial selection.*
#'
#' The initial list tidal constituents (prior to the pruning phase described
#' below) to be used in the analysis are specified as follows;
#' see \dQuote{Constituent Naming Convention}.
#'
#' 1. If `constituents` is not provided, then the constituent
#' list will be made up of the 69 constituents designated by Foreman as "standard".
#' These include astronomical frequencies and some shallow-water frequencies,
#' and are as follows: `c("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", "SO1", "OO1", "UPS1",
#' "OQ2", "EPS2", "2N2", "MU2", "N2", "NU2", "GAM2", "H1", "M2", "H2", "MKS2",
#' "LDA2", "L2", "T2", "S2", "R2", "K2", "MSN2", "ETA2", "MO3", "M3", "SO3",
#' "MK3", "SK3", "MN4", "M4", "SN4", "MS4", "MK4", "S4", "SK4", "2MK5", "2SK5",
#' "2MN6", "M6", "2MS6", "2MK6", "2SM6", "MSK6", "3MK7", "M8")`.
#'
#' 2. If the first item in `constituents` is the string
#' `"standard"`, then a provisional list is set up as in Case 1, and then
#' the (optional) rest of the elements of `constituents` are examined, in
#' order. Each of these constituents is based on the name of a tidal
#' constituent in the Foreman (1978) notation. (To get the list, execute
#' `data(tidedata)` and then execute `cat(tideData$name)`.) Each
#' named constituent is added to the existing list, if it is not already there.
#' But, if the constituent is preceded by a minus sign, then it is removed
#' from the list (if it is already there). Thus, for example, a user
#' might specify `constituents=c("standard", "-M2", "ST32")` to remove the M2
#' constituent and add the ST32 constituent.
#'
#' 3. If the first item is not `"standard"`, then the list of
#' constituents is processed as in Case 2, but without starting with the
#' standard list. As an example, `constituents=c("K1", "M2")` would fit
#' for just the K1 and M2 components. (It would be strange to use a minus sign
#' to remove items from the list, but the function allows that.)
#'
#' In each of the above cases, the list is reordered in frequency prior to the
#' analysis, so that the results of [summary,tidem-method()] will be in a
#' familiar form.
#'
#' *Phase 2: pruning based on Rayleigh's criterion.*
#'
#' Once the initial constituent list is determined during Phase 1,
#' `tidem` applies the Rayleigh criterion, which holds that
#' constituents of frequencies \eqn{f_1}{f1} and \eqn{f_2}{f2} cannot be
#' resolved unless the time series spans a time interval of at least
#' \eqn{rc/(f_1-f_2)}{rc/(f1-f2)}. The details of the decision of which
#' constituent to prune is fairly complicated, involving decisions
#' based on a comparison table. The details of this process are provided
#' by Foreman (1978).
#'
#' *Phase 3: optionally overruling Rayleigh's criterion.*
#'
#' The pruning list from phase 2 can be overruled by the user. This
#' is done by retaining any constituents that the user has named
#' in the `constituents` parameter. This action
#' was added on 2017-12-27, to make `tidem` behave the same
#' way as the Foreman (1978) code, as illustrated in his
#' Appendices 7.2 and 7.3. (As an aside, his Appendix 7.3 has some errors,
#' e.g. the frequency for the 2SK5 constituent is listed there (p58) as
#' 0.20844743, but it is listed as 0.2084474129 in his Appendix 7.1 (p41).
#' For this reason, the frequency comparison is relaxed to a `tol`
#' value of `1e-7` in a portion of the oce test suite
#' (see `tests/testthat/test_tidem.R` in the source).
#'
#' *Comments on phases 2 and 3*
#'
#' A specific example may be of help in understanding the removal of unresolvable
#' constituents. For example, the `data(sealevel)` dataset is of length
#' 6718 hours, and this is too short to resolve the full list of constituents,
#' with the conventional (and, really, necessary) limit of `rc=1`.
#' From Table 1 of Foreman (1978), this timeseries is too short to resolve the
#' `SA` constituent, so that `SA` will not be in the resultant.
#' Similarly, Table 2 of Foreman (1978) dictates the removal of
#' `PI1`, `S1` and `PSI1` from the list. And, finally,
#' Table 3 of Foreman (1978) dictates the removal of
#' `H1`, `H2`, `T2` and `R2`, and since that document
#' suggests that `GAM2` be subsumed into `H1`,
#' then if `H1` is already being deleted, then
#' `GAM2` will also be deleted.
#'
#' @param t either a vector of times or a [sealevel-class] object
#' (as created with [read.sealevel()] or [as.sealevel()]).
#' In the first case, `x` must be provided. In the second
#' case, though, any `x` that is provided will be ignored,
#' because sealevel objects contain both `time` and water
#' `elevation`, and the latter is used for `x`.
#'
#' @param x an optional numerical vector holding something that varies with
#' time. This is ignored if `t` is a [sealevel-class] object,
#' because it is inferred automatically, using `t[["elevation"]]`.
#'
#' @param constituents an optional character vector holding the names
#' of tidal constituents to which the fit is done; see \dQuote{Details}
#' and \dQuote{Constituent Naming Convention}.
#'
#' @template tideconst
#'
#' @param infer a list of constituents to be inferred from
#' fitted constituents according to the method outlined
#' in Section 2.3.4 of Foreman (1978).
#' If `infer` is `NULL`, the default, then
#' no such inferences are made. Otherwise, some constituents
#' are computed based on other constituents, instead of being
#' determined by regression at the proper frequency.
#' If provided, `infer` must be a list containing
#' four elements:
#' `name`, a vector of strings naming the constituents to be
#' inferred; `from`, a vector of strings naming the fitted
#' constituents used as the sources for those inferences (these
#' source constituents are added to the regression list, if they
#' are not already there);
#' `amp`, a numerical vector of factors to be applied to the
#' source amplitudes; and `phase`, a numerical vector of angles,
#' in degrees, to be subtracted from the source phases. For example,
#' Following Foreman (1998), if any of the `name` items
#' have already been computed, then the suggested inference is ignored,
#' and the already-computed values are used.
#' \preformatted{
#' infer=list(name=c("P1","K2"),
#' from=c("K1", "S2"),
#' amp=c(0.33093, 0.27215),
#' phase=c(-7.07, -22.4)
#' }
#' means that the amplitude of `P1` will be set as 0.33093 times the calculated amplitude
#' of `K1`, and that the `P1` phase will be set to the `K1` phase,
#' minus an offset of `-7.07` degrees.
#' (This example is used in the Foreman (1978) discussion of a
#' Fortran analysis code and also in Pawlowicz et al. (2002) discussion
#' of the T_TIDE Matlab code.
#' Rounded to the 0.1mm resolution of values reported in Foreman (1978)
#' and Pawlowicz et al. (2002),
#' the `tidem` results have root-mean-square amplitude difference
#' to Foreman's (1978) Appendix 7.3 of 0.06mm; by comparison,
#' the results in Table 1 of Pawlowicz et al. (2002) agree with Foreman's
#' results to RMS difference 0.04mm.)
#'
#' @param latitude if provided, the latitude of the observations. If not
#' provided, `tidem` will try to infer this from the first argument,
#' if it is a [sealevel-class] object.
#'
#' @param rc the value of the coefficient in the Rayleigh criterion.
#'
#' @param regress function to be used for regression, by default
#' [lm()], but could be for example `rlm` from the
#' `MASS` package.
#'
#' @template debugTemplate
#'
#' @return An object of [tidem-class], consisting of
#' \item{const}{constituent number, e.g. 1 for `Z0`, 1 for `SA`,
#' etc.} \item{model}{the regression model} \item{name}{a vector of constituent
#' names, in non-subscript format, e.g. "`M2`".} \item{frequency}{a vector
#' of constituent frequencies, in inverse hours.} \item{amplitude}{a vector of
#' fitted constituent amplitudes, in metres.} \item{phase}{a vector of fitted
#' constituent phase. NOTE: The definition of phase is likely to change as
#' this function evolves. For now, it is phase with respect to the first data
#' sample.} \item{p}{a vector containing a sort of p value for each
#' constituent. This is calculated as the average of the p values for the
#' sine() and cosine() portions used in fitting; whether it makes any sense is
#' an open question.}
#'
#' @section Application to non-hourly data:
#'
#' The framework on which [tidem()] rests on the assumption of data
#' that have been sampled on a 1-hour interval (see e.g. Foreman, 1977).
#' Since regression (as opposed to spectral analysis) is used to infer
#' the amplitude and phase of tidal constituents, data gaps do not pose
#' a serious problem. Sampling intervals under an hour are also not a
#' problem. However, trying to use [tidem()] on time series that are
#' sampled at uniform intervals that exceed 1 hour can lead to results
#' that are difficult to interpret. For example, some drifter data are
#' sampled at a 6-hour interval. This makes it impossible to fit for the
#' S4 component (which has exactly 6 cycles per day), because the method
#' works by constructing sine and cosine series at tidal frequencies and
#' using these as the basis for regression. Each of these series will have
#' a constant value through the constructed time, and regression cannot handle
#' that (in addition to a constant-value constructed series that is used to fit
#' for the Z0 constituent). [tidem()] tries to handle such problems by examining
#' the range of the constructed sine and cosine time-series, omitting any
#' constituents that yield near-constant values in either of these. Messages are
#' issued if this problem is encountered. This prevents failure of the regression,
#' and the predictions of the regression seem to represent the data reasonably well,
#' but the inferred constituent amplitudes are not physically reasonable. Cautious
#' use of [tidem()] to infer individual constituents might be warranted, but
#' users must be aware that the results will be difficult to interpret. The tool
#' is simply not designed for this use.
#'
#' @section Bugs:
#'
#' \enumerate{
#' \item This function is not fully developed yet, and both the
#' form of the call and the results of the calculation may change.
#'
#' \item The reported `p` value may make no sense at all, and it might be
#' removed in a future version of this function. Perhaps a significance level
#' should be presented, as in the software developed by both Foreman and
#' Pawlowicz.
#'
#' }
#'
#' @section Constituent Naming Convention:
#'
#' `tidem` uses constituent names that follow the convention
#' set by Foreman (1978). This convention is slightly different
#' from that used in the T-TIDE package of Pawlowicz et al.
#' (2002), with Foreman's `UPS1` and `M8` becoming
#' `UPSI` and `MS` in T-TIDE. To permit the use of either notation,
#' [tidem()] uses [tidemConstituentNameFix()] to
#' convert from T-TIDE names to the
#' Foreman names, issuing warnings when doing so.
#'
#' @section Agreement with `T_TIDE` results:
#'
#' The `tidem` amplitude and phase results, obtained with
#' \preformatted{
#' tidem(sealevelTuktoyaktuk, constituents=c("standard", "M10"),
#' infer=list(name=c("P1", "K2"),
#' from=c("K1", "S2"),
#' amp=c(0.33093, 0.27215),
#' phase=c(-7.07, -22.40)))
#' }
#' match the `T_TIDE` values listed in
#' Table 1 of Pawlowicz et al. (2002),
#' after rounding amplitude and phase to 4 and 2 digits past
#' the decimal place, respectively, to match the format of
#' that table.
#'
#' @author Dan Kelley
#'
#' @references
#'
#' Foreman, M G., 1977 (revised 1996).
#' Manual for Tidal Heights Analysis and Prediction.
#' Pacific Marine Science Report 77-10.
#' British Columbia, Canada: Institute of Ocean Sciences, Patricia Bay.
#'
#' Foreman, M. G. G., 1978.
#' Manual for Tidal Currents Analysis and Prediction.
#' Pacific Marine Science Report 78-6.
#' British Columbia, Canada: Institute of Ocean Sciences, Patricia Bay,
#'
#' Foreman, M. G. G., Neufeld, E. T., 1991.
#' Harmonic tidal analyses of long time series.
#' International Hydrographic Review, 68 (1), 95-108.
#'
#' Leffler, K. E. and D. A. Jay, 2009. Enhancing tidal harmonic analysis:
#' Robust (hybrid) solutions. Continental Shelf Research, 29(1):78-88.
#'
#' Pawlowicz, Rich, Bob Beardsley, and Steve Lentz, 2002. Classical tidal
#' harmonic analysis including error estimates in MATLAB using `T_TIDE`.
#' Computers and Geosciences, 28, 929-937.
#'
#' @examples
#' library(oce)
#' # The demonstration time series from Foreman (1978),
#' # also used in T_TIDE (Pawlowicz, 2002).
#' data(sealevelTuktoyaktuk)
#' tide <- tidem(sealevelTuktoyaktuk)
#' summary(tide)
#'
#' # AIC analysis
#' extractAIC(tide[["model"]])
#'
#' # Fake data at M2
#' library(oce)
#' data("tidedata")
#' M2 <- with(tidedata$const, freq[name == "M2"])
#' t <- seq(0, 10 * 86400, 3600)
#' eta <- sin(M2 * t * 2 * pi / 3600)
#' sl <- as.sealevel(eta)
#' m <- tidem(sl)
#' summary(m)
#'
#' @family things related to tides
tidem <- function(
t, x, constituents, infer = NULL, latitude = NULL,
rc = 1, regress = lm, debug = getOption("oceDebug")) {
oceDebug(debug, "tidem(t, x,\n", sep = "", unindent = 1)
oceDebug(debug, " constituents=",
if (missing(constituents)) {
"(missing)"
} else {
paste("c(\"", paste(constituents, collapse = "\", \""), "\")", sep = "")
},
",\n",
sep = "", unindent = 1
)
oceDebug(debug, " latitude=", if (is.null(latitude)) "NULL" else latitude, ",\n", sep = "", unindent = 1)
oceDebug(debug, " rc=", rc, ",\n", sep = "", unindent = 1)
oceDebug(debug, " debug=", debug, ") {\n", sep = "", unindent = 1)
cl <- match.call()
if (missing(t)) {
stop("must supply 't', either a vector of times or a sealevel object")
}
if (inherits(t, "sealevel")) {
sl <- t
t <- sl[["time"]]
x <- sl[["elevation"]]
if (is.null(latitude)) {
latitude <- sl[["latitude"]]
}
} else {
if (missing(x)) {
stop("must supply 'x', since the first argument is not a sealevel object")
}
if (inherits(x, "POSIXt")) {
warning("tidem() switching first 2 args to permit old-style usage")
tmp <- x
x <- t
t <- tmp
}
if (length(x) != length(t)) {
stop("lengths of 'x' and 't' must match, but they are ", length(x), " and ", length(t), " respectively")
}
if (inherits(t, "POSIXt")) {
t <- as.POSIXct(t)
} else {
stop("t must be a vector of POSIXt times")
}
sl <- as.sealevel(x, t)
}
# Check infer extensively, to prevent weird errors for e.g. an improperly-named
# constituent.
data("tidedata", package = "oce", envir = environment())
tidedata <- get("tidedata") # , pos=globalenv())
tc <- tidedata$const
ntc <- length(tc$name)
# 'infer' sanity-check
if (!is.null(infer)) {
if (!is.list(infer)) {
stop("infer must be a list")
}
if (length(infer) != 4) {
stop("infer must hold 4 elements")
}
if (!all.equal(sort(names(infer)), c("amp", "from", "name", "phase"))) {
stop("infer must contain 'name', 'from', 'amp', and 'phase', and nothing else")
}
if (!is.character(infer$name)) {
stop("infer$name must be a vector of character strings")
}
infer$name <- tidemConstituentNameFix(infer$name)
if (!is.character(infer$from)) {
stop("infer$from must be a vector of character strings")
}
infer$from <- tidemConstituentNameFix(infer$from)
if (length(infer$name) != length(infer$from)) {
stop("lengths of infer$name and infer$from must be equal")
}
if (length(infer$name) != length(infer$amp)) {
stop("lengths of infer$name and infer$amp must be equal")
}
if (length(infer$name) != length(infer$phase)) {
stop("lengths of infer$name and infer$phase must be equal")
}
for (n in infer$name) {
if (!(n %in% tc$name)) {
stop("infer$name value \"", n, "\" is not a known tidal constituent; see const$name in ?tidedata")
}
}
for (n in infer$from) {
if (!(n %in% tc$name)) {
stop("infer$from value \"", n, "\" is not a known tidal constituent; see const$name in ?tidedata")
}
}
} # 'infer' sanity-check
# The arguments seem to be OK, so start the actual analysis now.
startTime <- t[1]
endTime <- tail(t, 1)
years <- as.numeric(difftime(endTime, startTime, units = "secs")) / 86400 / 365.25
if (years > 18.6) {
warning("Time series spans 18.6 years, but tidem() is ignoring this important fact")
}
standard <- tc$ikmpr > 0
addedConstituents <- NULL
if (missing(constituents)) {
# Default 'name', 'freq', 'kmpr' and 'indices'; drop Z0 because we infer it directly.
name <- tc$name[standard]
freq <- tc$freq[standard]
kmpr <- tc$kmpr[standard]
indices <- seq(1:ntc)[standard] # NB. Z0 need not be dropped; we work with indices later anyway
oceDebug(debug, "starting with ", length(name), " default constituents: ", paste(name, collapse = " "), sep = "", "\n")
} else {
# Build up 'name'; later, infer 'indices' and thereby 'freq' and 'kmpr'.
name <- NULL
for (i in seq_along(constituents)) {
if (constituents[i] == "standard") {
# must be first!
if (i != 1) {
stop("\"standard\" must occur first in constituents list")
}
name <- tc$name[standard]
} else {
constituents <- tidemConstituentNameFix(constituents)
if (substr(constituents[i], 1, 1) == "-") {
# Case 1: removal. Require a valid name, and warn if not in the list already.
nameRemove <- substr(constituents[i], 2, nchar(constituents[i]))
if (1 != sum(tc$name == nameRemove)) {
stop(
"\"", nameRemove, "\" is not a known tidal constituent; try one of: ",
paste(tc$name, collapse = "\" \"")
)
}
remove <- which(name == nameRemove)
oceDebug(debug > 1, "removed \"", nameRemove, "\"\n", sep = "")
if (0 == length(remove)) {
warning("\"", nameRemove, "\" is not in the list of constituents currently under study", sep = "")
} else {
name <- name[-remove]
}
} else {
# Case 2: addition. Require a valid name, and ignore repeat requests.
add <- which(tc$name == constituents[i])
if (1 != length(add)) {
stop("\"", constituents[i], "\" is not a known tidal constituent (line 1093)")
}
if (!(constituents[i] %in% name)) {
name <- c(name, tc$name[add])
addedConstituents <- c(addedConstituents, add)
}
}
}
}
}
oceDebug(debug, "will fit for ", length(name), " constituents: ", paste(name, collapse = " "), "\n", sep = "")
# We let users add "Z0" as a constituent, but we remove it now since the
# regression will have an intercept and that becomes Z0.
fitForZ0 <- "Z0" %in% name
oceDebug(debug, "fitForZ0=", fitForZ0, "\n")
# All of the names should be valid from the above, but to protect against code changes,
# we check to be sure.
if (any(!(name %in% tc$name))) {
bad <- NULL
for (n in name) {
if (!(n %in% tc$name)) {
bad <- c(bad, n)
}
}
stop("unknown constituent names: ", paste(bad, collapse = " "), " (please report this error to developer)")
}
# Infer indices from the names, sort them as in the tidal-constituent
# table, tc, and then look up freq and kmpr from that table.
indices <- sort(unlist(lapply(name, function(name) which(tc$name == name))))
name <- tc$name[indices]
freq <- tc$freq[indices]
kmpr <- tc$kmpr[indices]
nc <- length(name)
# Check Rayleigh criterion
interval <- diff(range(as.numeric(sl@data$time), na.rm = TRUE)) / 3600 # in hours
oceDebug(debug, "interval=", interval, " hours\n")
dropTerm <- NULL
for (i in 1:nc) {
cc <- which(tc$name == kmpr[i])
if (length(cc)) {
cannotFit <- (interval * abs(freq[i] - tc$freq[cc])) < rc
oceDebug(debug, "i=", i, ", name=", name[i], ", kmpr[", i, "]=", kmpr[i], ", cannotFit=", cannotFit, "\n", sep = "")
if (cannotFit) {
dropTerm <- c(dropTerm, i)
}
}
}
oceDebug(debug, "before trimming constituents for Rayleigh condition, name[1:", length(name), "]=", paste(name, collapse = " "), sep = "", "\n")
if (length(dropTerm) > 0) {
# Bookmark 1A (see also 1B: link up variables)
message("Note: the tidal record is too short to fit for constituents: ", paste(name[dropTerm], collapse = ", "))
indices <- indices[-dropTerm]
name <- name[-dropTerm]
freq <- freq[-dropTerm]
kmpr <- kmpr[-dropTerm]
}
oceDebug(debug, "after trimming constituents for Rayleigh condition, name[1:", length(name), "]=", paste(name, collapse = " "), sep = "", "\n")
# Ensure that any added constituents are in the list, i.e. prevent
# the Rayleigh criterion from trimming them. (Before work on
# issue 1350, they would simply be dropped if they failed the Rayleigh
# criterion. Although that was a sensible choice, it was decided
# on 2017-12-27, whilst working on issue 1350, to make tidem() do the
# the same thing as the Foreman 1978 code as exemplified in his appendices
# 7.2 and 7.3.)
if (length(addedConstituents)) {
oceDebug(debug, "addedConstituents=", paste(addedConstituents, collapse = " "), "\n")
for (a in addedConstituents) {
# Avoid adding constituents that we already have
if (!(tc$name[a] %in% name)) {
if (debug > 0) {
message("ADDING:")
message(" tc$name[a=", a, "]=\"", tc$name[a], "\"", sep = "")
message(" tc$freq[a=", a, "]=\"", tc$freq[a], "\"", sep = "")
message(" tc$kmpr[a=", a, "]=\"", tc$kmpr[a], "\"", sep = "")
}
indices <- c(indices, which(tc$name == name[a]))
name <- c(name, tc$name[a])
freq <- c(freq, tc$freq[a])
kmpr <- c(kmpr, tc$kmpr[a])
}
}
}
oceDebug(debug, "after adding new constituents, ", vectorShow(name))
oceDebug(debug, "after adding new constituents, ", vectorShow(freq))
# Ensure that we fit for any infer$from constituents, *regardless* of whether
# those consitituents are permitted by the Rayleigh criterion.
if (!is.null(infer)) {
for (n in c(infer$from)) {
if (!(n %in% name)) {
a <- which(tc$name == n)
indices <- c(indices, a)
name <- c(name, tc$name[a])
freq <- c(freq, tc$freq[a])
kmpr <- c(kmpr, tc$kmpr[a])
message("fitting for infer$from=", n, ", even though the Rayleigh Criterion would exclude it")
}
}
}
oceDebug(debug, "after handling 'infer', ", vectorShow(name))
oceDebug(debug, "after handling 'infer', ", vectorShow(freq))
# sort constituents by index (which, among other things, ensures that Z0 is at the start, if it exists)
oindices <- order(indices)
indices <- indices[oindices]
name <- name[oindices]
freq <- freq[oindices]
kmpr <- kmpr[oindices]
nc <- length(name)
oceDebug(debug, "after reordering indices, ", vectorShow(name))
oceDebug(debug, "after reordering indices, ", vectorShow(freq))
rm(oindices) # clean up namespace
if (0L == nc) {
stop("cannot fit for any constituents")
}
elevation <- sl[["elevation"]]
time <- sl[["time"]]
nt <- length(elevation)
x <- array(dim = c(nt, 2 * nc))
oceDebug(debug, vectorShow(nc))
oceDebug(debug, vectorShow(dim(x)))
x[, 1] <- rep(1, nt)
pi <- 4 * atan2(1, 1)
rpd <- atan2(1, 1) / 45 # radians per degree
# tRef <- ISOdate(1899, 12, 31, 12, 0, 0, tz="UTC") # was this ever used?
# tRef <- centralTime # used previous to "dk" branch early 2018
tRef <- numberAsPOSIXct(3600 * round(mean(as.numeric(time, tz = "UTC")) / 3600), tz = "UTC")
hour2pi <- 2 * pi * (as.numeric(time) - as.numeric(tRef)) / 3600
oceDebug(debug, "tRef=", tRef, ", nc=", nc, ", length(name)=", length(name), "\n")
# The sameCriterion was added in January 2023, after this portion of tidem()
# had been stable for a decade or more. Also, I am trying to make tidem()
# do something that (I think) neither Foreman's code nor t-tide does.
# Therefore, I am putting a long comment here! Please also see
# https://github.com/dankelley/oce/issues/2034 to learn where this idea of
# using sameCriterion got started.
#
# In the code below, we keep track of any constituents for which the
# constructed C and S vectors will be constant (or nearly so). This came up
# because I was trying to analyse data sampled on a 6-hour interval. That
# works out to match the S4 period (approximately, but sampling times are
# not given to an infinite number of digits and neither is the S4 period in
# this code). Imagine a perfect match. Then, S and C would both be constant
# across the times used in the loop. And that's bad because the regression
# already has a constant column (for Z0). The lm() computation then gives
# NA values for coefficients, p value, etc. for the repeat. (Luckily, it
# doesn't just say the matrix cannot be inverted and die!) The approximate
# criterion is based on a test case in the context of the local machine
# epsilon. I use a criterion of 0.01 because for "good" cases (and my
# 6-hour test file for drifter data) have C and S span of very nearly 2, but for "bad" cases, it is
# of order e-7 or so. Any dividing line would do, I think, but maybe for short records
# the span might not get to be -1 to +1 and so I am choosing 1e-2 as a criterion. This value
# might need to be revisited.
iBad <- NULL
icBad <- NULL
sdCriterion <- 1e-2
oceDebug(debug, vectorShow(sdCriterion))
# 20230122 danS <- danC <- NULL
for (i in 1:nc) {
# 20230122 oceDebug(debug+1, "setting ", i, "-th coefficient (name=", name[i], " period=", 1/freq[i], " h)", "\n", sep="")
ft <- freq[i] * hour2pi
C <- cos(ft)
S <- sin(ft)
sdS <- sd(S)
sdC <- sd(C)
# 20230122 danS <- c(danS, sdS) # FIXME: remove
# 20230122 danC <- c(danC, sdC) # FIXME: remove
# 20230122 oceDebug(debug+1, sprintf(" sdC %.4g; sdS %.4g\n", sdC, sdS))
# Find whether anything is uncomputable. We ignore Z0 because that is handled later.
if (name[i] != "Z0" && (sdS < sdCriterion || sdC < sdCriterion)) {
oceDebug(debug, " ** uncomputable at ", name[i], " (period ", 1 / freq[i], "h) **\n")
icBad <- c(icBad, 1 + 2 * (i - 1))
icBad <- c(icBad, 2 + 2 * (i - 1))
iBad <- c(iBad, i)
}
x[, 1 + 2 * (i - 1)] <- C
x[, 2 + 2 * (i - 1)] <- S
}
name2 <- matrix(rbind(paste(name, "_C", sep = ""), paste(name, "_S", sep = "")), nrow = length(name), ncol = 2)
dim(name2) <- c(2 * length(name), 1)
colnames(x) <- name2
#<< cat("The following should be omitted: ", vectorShow(name2[uncomputable], n=20))
#<< cat("Their indices are", vectorShow(uncomputable, n=20))
oceDebug(debug, "cleaning up 'x' matrix prior to doing regression\n")
# Remove problematic constutuents; see https://github.com/dankelley/oce/issues/2034
if (length(iBad)) {
message("Note: the sampling interval is too coarse to fit for constituents: ", paste(name[iBad], collapse = ", "))
# 20230122 browser()
# 20230122 A<-colnames(x)
# 20230122 B<-name2
# 20230122 C<-name
indices <- indices[-iBad]
name <- name[-iBad]
freq <- freq[-iBad]
kmpr <- kmpr[-iBad]
nc <- nc - length(iBad)
oceDebug(debug, "Removing: ", paste(name2[icBad], collapse = ", "))
x <- x[, -icBad]
name2 <- name2[-icBad]
# 20230122 # A check with 6-h case (sandbox/dk/tidem/04_restrict_coefficients.R)
# 20230122 print(A[!A%in%colnames(x)]) # should be S2_C S2_S S4_C S4_S
# 20230122 print(B[!B%in%name2]) # should be S2_C S2_S S4_C S4_S
# 20230122 print(C[!C%in%name]) # should be S2 S4
}
# Remove the sine() part of the Z0 constituent, which makes no sense for a constant.
if ("Z0_S" %in% colnames(x)) {
x <- x[, -which("Z0_S" == colnames(x))]
oceDebug(debug, "model has Z0, so trimming the sin(freq*time) column\n")
}
oceDebug(debug, "about to do regression\n")
oceDebug(debug, vectorShow(colnames(x)))
model <- regress(elevation ~ x - 1, na.action = na.exclude)
#>> browser()
if (debug > 0) {
cat("regression worked OK; the results are as follows:\n")
print(summary(model))
}
coef <- model$coefficients
p.all <- if (4L == ncol(summary(model)$coefficients)) {
summary(model)$coefficients[, 4L]
} else {
rep(NA, length = 1L + nc)
}
amplitude <- phase <- p <- vector("numeric", length = nc)
oceDebug(debug, vectorShow(nc))
oceDebug(debug, vectorShow(phase))
oceDebug(debug, vectorShow(name))
ic <- 1
for (i in seq_len(nc)) {
# Z0 has Z0_C but no Z0_S (since zero is a degenerate regression variable)
if (name[i] == "Z0") {
if (i != 1) {
stop("Z0 should be at the start of the regression coefficients. Please report this to developer.")
}
j <- which(tidedata$const$name == name[i])
vuf <- tidemVuf(tRef, j = j, latitude = latitude)
amplitude[i] <- coef[ic]
phase[i] <- 0
p[i] <- p.all[ic]
oceDebug(debug, "processed coefs at i=", i, ", ic=", ic,
", name=", name[i], ", f=", vuf$f, ", angle adj=", (vuf$u + vuf$v) * 360,
", amplitude=", amplitude[i], ", phase=", phase[i], ", p=", p[i], "\n",
sep = ""
)
ic <- ic + 1 # only skip forward 1 since Z0 takes 1 column (contrast below)
} else {
C <- coef[ic] # coefficient on cos(t)
S <- coef[ic + 1] # coefficient on sin(t)
amplitude[i] <- sqrt(S^2 + C^2)
# Calculate phase from the coefficients on sin() and cos(). Generally,
# cos(t - phase) == cos(phase)*cos(t) + sin(phase)*sin(t)
# By the definition of the regression model, we have
# cos(t - phase) == c * cos(t) + s * sin(t)
# and thus phase is defined by
# tan(phase) == s/c
phase[i] <- atan2(S, C)
# Adjust amplitude phase, as in ~/src/foreman/tide12_r2.f:405
j <- which(tidedata$const$name == name[i])
vuf <- tidemVuf(tRef, j = j, latitude = latitude)
amplitude[i] <- amplitude[i] / vuf$f
p[i] <- 0.5 * (p.all[ic + 1] + p.all[ic])
#<<cat("i=", i, ", ic=", ic, ", name=\"", name[i], "\", p.all[ic]=",
# p.all[ic], ", p.all[ic+1]=", p.all[ic+1], ", we set p[i]=", p[i],
# " (the mean)\n", sep="")
oceDebug(debug, "processed coefs at i=", i, ", ic=", ic, ", name=", name[i],
", S=", S, ", C=", C, ", f=", vuf$f, ", angle adj=", (vuf$u + vuf$v) * 360, ",
amplitude=", amplitude[i], ", phase=", phase[i], ", p=", p[i], "\n",
sep = ""
)
ic <- ic + 2 # skip forward 2 since non-Z0 takes 2 columns (contrast above)
}
}
oceDebug(debug, vectorShow(phase))
phase <- phase * 180 / pi
phase <- ifelse(phase < -360, 720 + phase, phase)
phase <- ifelse(phase < 0, 360 + phase, phase)
# Do Greenwich phase correction, if `infer` is TRUE
C <- unlist(lapply(name, function(n) which(n == tidedata$const$name)))
vuf <- tidemVuf(tRef, j = C, latitude = latitude)
oceDebug(debug, vectorShow(freq))
oceDebug(debug, vectorShow(phase))
oceDebug(debug, vectorShow(vuf$u))
oceDebug(debug, vectorShow(vuf$v))
phase <- phase + (vuf$v + vuf$u) * 360
phase <- ifelse(phase < 0, phase + 360, phase)
phase <- ifelse(phase > 360, phase - 360, phase)
# Handle (optional) inferred constituents. We know that
# this list is well-formed because of extensive tests near
# the start of this function.
if (!is.null(infer)) {
if (debug > 0) {
cat("BEFORE inference:\n")
print(data.frame(name = name, freq = round(freq, 6), amplitude = round(amplitude, 4)))
}
for (n in seq_along(infer$name)) {
oceDebug(debug, "n=", n, "; handling inferred constituent ", infer$name[n], "\n")
iname <- which(tc$name == infer$name[n])[1]
oceDebug(debug, "infer$name[\"", n, "\"]=\"", infer$name[n], "\" yields iname=", iname, "\n", sep = "")
oceDebug(debug, "iname=", iname, "\n")
# NOTE: we know that infer$from has been fitted for, because we forced it to be,
# irrespective of the Rayleight Criterion. Still, we test here, in case the code
# gets changed later.
if (infer$from[n] %in% name) {
ifrom <- which(name == infer$from[n])[1]
if (infer$name[n] %in% name) {
# message("name already in list")
iname <- which(name == infer$name[n])[1]
# Update, skipping 'indices', 'name' and 'freq', since they are already OK.
amplitude[iname] <- infer$amp[n] * amplitude[ifrom]
phase[iname] <- phase[ifrom] - infer$phase[n]
p[iname] <- p[ifrom]
oceDebug(debug, "replace existing ", name[iname], " based on ", name[ifrom], " (", freq[ifrom], " cph)\n", sep = "")
warning(
"inferring \"", infer$name[n],
"\" which is already included in the regression. Foreman says to skip it; unsure on what T_TIDE does\n"
)
} else {
# We append new values at the end, knowing that they will get
# shifted back to their proper positions when we reorder the
# whole solution, after handling these inferences.
#
# The first step is to adjust the amp and phase of infer$from; this
# is done based on formulae in Foreman (1978) sec 2.3.4. It looks
# as though t_tide.m on and after about line 472 is doing a
# similar thing, although the numbers do not agree exactly,
# as shown in issue 1351, code 1351c.R.
#
# Foreman (1978) sec 2.3.4.
# Notation: suffices "1" and "2" refer to "from" and "name" here.
i1 <- which(tc$name == infer$from[n])[1]
i2 <- which(tc$name == infer$name[n])[1]
oceDebug(
debug, "tRef=", format(tRef, "%Y-%m-%d %H:%M:%S"),
", i1=", i1, ", i2=", i2, ", lat=", latitude, "\n"
)
vuf12 <- tidemVuf(tRef, c(i1, i2), latitude = latitude)
# vuf2 <- tidemVuf(tRef, i2, latitude=latitude)
f1 <- vuf12$f[1]
f2 <- vuf12$f[2]
oceDebug(debug, "f1=", f1, ", f2=", f2, "\n")
# FIXME: what is unit of u and v? t_tide.m:482 suggests it is degrees
# Foreman's tide12_r2.f:399 suggests U and V are in cycles,
# and this is consistent with Pawlowicz's t_tide.m:451
# We convert vu1 and vu2 to be in degrees, as t_tide.m does
vu1 <- (vuf12$v[1] + vuf12$u[1]) * 360
vu2 <- (vuf12$v[2] + vuf12$u[2]) * 360
oceDebug(debug, "vu1=", vu1, ", vu2=", vu2, "\n")
sigma1 <- tc$freq[i1]
sigma2 <- tc$freq[i2]
oceDebug(debug, "sigma1=", sigma1, ", sigma2=", sigma2, "\n")
# tmp is pi*N*(sigma2-sigma1) in Foreman
tmp <- pi * interval * (sigma2 - sigma1)
r12 <- infer$amp[n]
# FIXME: sign for Foreman?
zeta <- infer$phase[n]
S <- r12 * (f2 / f1) * sin(tmp) * sin(rpd * (vu2 - vu1 + zeta)) / tmp
C <- 1 + r12 * (f2 / f1) * sin(tmp) * cos(rpd * (vu2 - vu1 + zeta)) / tmp
oceDebug(debug, "tmp=", tmp, ", S=", S, ", C=", C, ", sqrt(S^2+C^2)=", sqrt(S^2 + C^2), "\n")
oceDebug(debug, infer$from[n], "amplitude, old=", amplitude[ifrom], ", new=", amplitude[ifrom] / sqrt(S^2 + C^2), "\n")
amplitude[ifrom] <- amplitude[ifrom] / sqrt(S^2 + C^2)
oceDebug(debug, infer$from[n], "phase, old=", phase[ifrom], ", new=", phase[ifrom] + atan2(S, C) / rpd, "\n")
phase[ifrom] <- phase[ifrom] + atan2(S, C) / rpd
# End of Foreman 1978 inference calculation. Now we can define 'name' i.t.o. 'from'
iname <- which(tc$name == infer$name[n])[1]
oceDebug(debug, "Below is inference for ", infer$name[n], " (index=", iname, ")\n")
indices <- c(indices, iname)
name <- c(name, infer$name[n])
freq <- c(freq, tc$freq[iname])
amplitudeInferred <- infer$amp[n] * amplitude[ifrom]
phaseInferred <- phase[ifrom] - infer$phase[n]
oceDebug(debug, " ", infer$name[n], "inferred amplitude=", amplitudeInferred, "\n")
oceDebug(debug, " ", infer$name[n], "inferred phase=", phaseInferred, "\n")
amplitude <- c(amplitude, amplitudeInferred)
phase <- c(phase, phaseInferred)
p <- c(p, p[ifrom])
oceDebug(debug, " create ", infer$name[n], " (index=", iname, ", ",
tc$freq[iname], " cph) based on ", name[ifrom],
" (index ", ifrom, ", ", freq[ifrom], " cph)\n",
sep = ""
)
}
} else {
stop("Internal error (please report): cannot infer ", infer$name[n], " from ", infer$from[n], " because the latter was not computed")
}
}
# reorder by original position in tc
o <- order(indices)
indices <- indices[o]
stopifnot(length(o) == length(name))
stopifnot(length(o) == length(freq))
stopifnot(length(o) == length(amplitude))
stopifnot(length(o) == length(phase))
stopifnot(length(o) == length(p))
name <- name[o]
freq <- freq[o]
amplitude <- amplitude[o]
phase <- phase[o]
p <- p[o]
rm(o)
if (debug > 0) {
cat("AFTER inference\n")
print(data.frame(name = name, freq = round(freq, 5), amplitude = round(amplitude, 4)))
}
}
phase <- phase %% 360
res <- new("tidem")
o <- order(freq)
res@data <- list(
model = model,
call = cl,
tRef = tRef,
const = indices[o],
name = name[o],
freq = freq[o],
amplitude = amplitude[o],
phase = phase[o],
p = p[o]
)
res@metadata$rc <- rc
res@metadata$version <- 2
res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep = "", collapse = ""))
oceDebug(debug, "} # tidem()\n", sep = "", unindent = 1)
res
}
#' Predict a Tidal Signal
#'
#' This creates a time-series of predicted tides, based on a
#' tidal model object that was created by [as.tidem()] or [tidem()].
#'
#' All the tidal constituents that are stored in `object` are used,
#' not just those that are statistically significant or that have amplitude
#' exceeding any particular value. In this respect, [predict.tidem()]
#' follows a pattern established by e.g. [predict.lm()]. Note that
#' the constituents in `object` are straightforward if it
#' was constructed with [as.tidem()], but considerably more complicated
#' for [tidem()], and so the documentation for the latter ought to
#' be studied closely, especially with regard to the Rayleigh criterion.
#'
#' @param object a [tidem-class] object.
#'
#' @param newdata vector of POSIXt times at which to make the
#' prediction. For models created with [tidem()],
#' the `newdata` argument is optional, and if it is not provided, then
#' the predictions are at the observation times given to
#' [tidem()]. However, `newdata` is required if [as.tidem()]
#' had been used to create `object`.
#'
#' @param \dots optional arguments passed on to children.
#'
#' @return A vector of predictions.
#'
#' @examples
#'
#' # Show non-tidal sealevel signal in Halifax Harbour during
#' # the year 2002. The spike resulted from Hurricane Juan.
#' library(oce)
#' data(sealevel)
#' time <- sealevel[["time"]]
#' elevation <- sealevel[["elevation"]]
#' prediction <- tidem(sealevel) |> predict()
#' oce.plot.ts(time, elevation - prediction)
#'
#' @section Sample of Usage:
#' \preformatted{
#' # prediction at specified times
#' data(sealevel)
#' m <- tidem(sealevel)
#' # Check fit over 2 days (interpolating to finer timescale)
#' look <- 1:48
#' time <- sealevel[["time"]]
#' elevation <- sealevel[["elevation"]]
#' oce.plot.ts(time[look], elevation[look])
#' # 360s = 10 minute timescale
#' t <- seq(from=time[1], to=time[max(look)], by=360)
#' lines(t, predict(m, newdata=t), col="red")
#' legend("topright", col=c("black","red"),
#' legend=c("data","model"),lwd=1)
#' }
#'
#' @author Dan Kelley
#'
#' @family things related to tides
predict.tidem <- function(object, newdata, ...) {
dots <- list(...)
debug <- if ("debug" %in% names(dots)) dots$debug else 0
oceDebug(debug, "predict.tidem() {\n", sep = "", unindent = 1)
if (!missing(newdata) && !inherits(newdata, "POSIXt")) {
stop("newdata must be of class POSIXt")
}
version <- object@metadata$version
if (!is.null(version) && version == 3) {
oceDebug(debug, "object@metadata$version is 3, so assuming the object was created by as.tidem()\n")
if (missing(newdata)) {
stop("must supply newdata because object was created with as.tidem()")
}
hour2pi <- 2 * pi * (as.numeric(newdata) - as.numeric(object[["tRef"]])) / 3600
oceDebug(debug, vectorShow(hour2pi))
# message("head(hour2pi): ", paste(head(hour2pi), collapse=" "))
nc <- length(object@data$name)
res <- rep(0, length(hour2pi))
for (i in seq_len(nc)) {
oceDebug(debug, "accounting for constituent[", i, "] = ", object@data$name[i], "\n", sep = "")
omega.t <- object@data$freq[i] * hour2pi
a <- object@data$amplitude[i] * sin(2 * pi * object@data$phase[i] / 360)
b <- object@data$amplitude[i] * cos(2 * pi * object@data$phase[i] / 360)
res <- res + a * sin(omega.t) + b * cos(omega.t)
}
} else if (!is.null(version) && version == 2) {
oceDebug(debug, "object@metadata$version is 2, so assuming the object was created by tidem()\n")
if (!missing(newdata) && !is.null(newdata)) {
oceDebug(debug, "newdata was provided\n")
freq <- object@data$freq
name <- object@data$name
nc <- length(freq)
tt <- as.numeric(as.POSIXct(newdata, tz = "UTC"))
nt <- length(tt)
x <- array(dim = c(nt, 2 * nc))
x[, 1] <- rep(1, nt)
hour2pi <- 2 * pi * (as.numeric(tt) - as.numeric(object[["tRef"]])) / 3600
for (i in 1:nc) {
omega.t <- freq[i] * hour2pi
x[, 2 * i - 1] <- cos(omega.t)
x[, 2 * i] <- sin(omega.t)
}
colnames(x) <- matrix(rbind(paste(name, "_C", sep = ""), paste(name, "_S", sep = "")), nrow = length(name), ncol = 2)
if ("Z0_S" %in% colnames(x)) {
x <- x[, -which("Z0_S" == colnames(x))]
oceDebug(debug, "model has Z0, so trimming the sin(freq*time) column\n")
}
res <- as.numeric(predict(object@data$model, newdata = list(x = x), ...))
} else {
oceDebug(debug, "newdata was not provided\n")
res <- as.numeric(predict(object@data$model, ...))
}
} else {
if (!("version" %in% names(object@metadata))) {
warning("prediction is being made based on an old object; it may be wrong\n")
}
res <- as.numeric(predict(object@data$model, ...))
}
oceDebug(debug, "} # predict.tidem()\n", sep = "", unindent = 1)
res
}
#' Get a Tidal Prediction From a WebTide Database
#'
#' Get a tidal prediction from a WebTide database. This only
#' works if the standalone WebTide application is installed,
#' and if it is installed in a standard location. The details
#' of installation are not within the oce purview.
#'
#' There are two methods of using this function.
#' *Case 1:* `action="map"`. In this case, if
#' `plot` is `FALSE`, a list is returned, containing
#' all the `node`s in the selected database, along with all
#' the `latitude`s and `longitude`s. This value is
#' also returned (silently) if `plot` is true, but in that case,
#' a plot is drawn to indicate the node locations. If `latitude` and
#' `longitude` are given, then the node nearest that spot is indicated on
#' the map; otherwise, if `node` is given, then the location of that
#' node is indicated. There is also a special case: if `node` is negative
#' and `interactive()` is `TRUE`,
#' then [locator()] is called, and the node nearest the spot
#' where the user clicks the mouse is indicated in the plot and in the
#' return value.
#'
#' *Case 2:* `action="predict"`. If `plot` is `FALSE`,
#' then a list is returned, indicating `time`, predicted
#' `elevation`, velocity components `u` and `v`,
#' `node` number, the name of the `basedir`, and
#' the `region`. If `plot` is `TRUE`, this list is returned
#' silently, and time-series plots are drawn for elevation, u, and v.
#'
#' Naturally, `webtide` will not work unless WebTide has been installed on
#' the computer.
#'
#' @param action An indication of the action, either `action="map"` to
#' draw a map or `action="predict"` to get a prediction; see
#' \dQuote{Details}.
#'
#' @param longitude,latitude optional location at which prediction is required (ignored if
#' `node` is given).
#'
#' @param node optional integer relating to a node in the database. If `node`
#' is given, then neither `latitude` nor `longitude` may be given.
#' If `node` is positive, then specifies indicates the node. If it is negative,
#' [locator()] is called so that the user can click (once) on the map, after
#' which the node is displayed on the map.
#'
#' @param time a vector of times (in the UTC timezone)
#' at which prediction is to be made.
#' If not supplied, this will be the week starting at the present time,
#' computed with [presentTime()], with a 15 minute increment.
#'
#' @param basedir directory containing the `WebTide` application.
#'
#' @param region database region, given as a directory name in the WebTide
#' directory. For example, `h3o` is for Halifax Harbour, `nwatl` is
#' for the northwest Atlantic, and `sshelf` is for the Scotian Shelf and
#' Gulf of Maine.
#'
#' @param plot boolean indicating whether to plot.
#'
#' @param tformat optional argument passed to [oce.plot.ts()], for
#' plot types that call that function. (See [strptime()] for the
#' format used.)
#'
#' @template debugTemplate
#'
#' @param \dots optional arguments passed to plotting functions. A common
#' example is to set `xlim` and `ylim`, to focus a map region.
#'
#' @return The value depends on `action`:
#'
#' * If `action="map"` the return value is a
#' list containing the index of the nearest node, along with the
#' `latitude` and `longitude` of that node. If
#' `plot` is `FALSE`, this value is returned invisibly.
#'
#' * If `action="predict"`, the return value is a list containing a vector
#' of times (`time`), as well as vectors of the predicted `elevation`
#' in metres and the predicted horizontal components of velocity, `u` and
#' `v`, along with the `node` number, and the `basedir` and
#' `region` as supplied to this function. If `plot` is `FALSE`,
#' this value is returned invisibly.
#'
#' @source The WebTide software may be downloaded for free at the
#' Department of Fisheries and Oceans (Canada) website at
#' `http://www.bio.gc.ca/science/research-recherche/ocean/webtide/index-en.php`
#' (checked February 2016 and May 2017).
#'
#' @section Caution:
#' WebTide is not an open-source application, so the present function was
#' designed based on little more than guesses about the WebTide file structure.
#' Users should be on the lookout for odd results.
#'
#' @section Sample of Usage:
#' \preformatted{
#' # needs WebTide at the system level
#' library(oce)
#' # 1. prediction at Halifax NS
#' longitude <- -63.57
#' latitude <- 44.65
#' prediction <- webtide("predict", longitude=longitude, latitude=latitude)
#' mtext(paste0("prediction at ", latitude, "N and ", longitude, "E"), line=0.75, side=3)
#' # 2. map
#' webtide(lon=-63.57,lat=44.65,xlim=c(-64,-63),ylim=c(43.0,46))
#' }
#'
#' @author Dan Kelley
#'
#' @family things related to tides
webtide <- function(
action = c("map", "predict"),
longitude, latitude, node, time,
basedir = getOption("webtide"),
region = "nwatl",
plot = TRUE, tformat, debug = getOption("oceDebug"), ...) {
debug <- max(0, min(floor(debug), 2))
oceDebug(debug, "webtide(action=\"", action, "\", ...)\n",
sep = "", unindent = 1, style = "bold"
)
rpd <- atan2(1, 1) / 45 # radians per degree
action <- match.arg(action)
nodeGiven <- !missing(node)
longitudeGiven <- !missing(longitude)
latitudeGiven <- !missing(latitude)
path <- paste(basedir, "/data/", region, sep = "")
# 2016-02-03: it seems that there are several possibilities for this filename.
suffices <- c(".nod", "ll.nod", "_ll.nod")
nodFiles <- paste(region, suffices, sep = "")
triangles <- NULL
for (nodFile in nodFiles) {
if (1 == length(list.files(path = path, pattern = nodFile))) {
triangles <- read.table(paste(path, nodFile, sep = "/"), col.names = c("triangle", "longitude", "latitude"), encoding = "latin1")
oceDebug(debug, "found webtide information in \"", nodFile, "\"\n", sep = "")
break
} else {
oceDebug(debug, "looked for webtide information in \"", nodFile, "\" but this file does not exist\n", sep = "")
}
}
if (is.null(triangles)) {
stop("cannot find WebTide data file; rerun with debug=1 to see the searched list")
}
if (action == "map") {
if (plot) {
asp <- 1 / cos(rpd * mean(range(triangles$latitude, na.rm = TRUE)))
par(mfrow = c(1, 1), mar = c(3, 3, 2, 1), mgp = c(2, 0.7, 0))
plot(triangles$longitude, triangles$latitude,
pch = 2, cex = 1 / 4, lwd = 1 / 8,
asp = asp, xlab = "", ylab = "", ...
)
# Try for a coastline of well-suite resolution, if we have ocedata installed.
usr <- par("usr")
bestcoastline <- coastlineBest(lonRange = usr[1:2], latRange = usr[3:4], debug = debug - 1)
oceDebug(debug, "coastlineBest() suggests using", bestcoastline, "as the coastline\n")
if (bestcoastline == "coastlineWorld") {
data(list = bestcoastline, package = "oce", envir = environment())
coastlineWorld <- get("coastlineWorld")
} else {
if (requireNamespace("ocedata", quietly = TRUE)) {
data(list = bestcoastline, package = "ocedata", envir = environment())
oceDebug(debug, "Using", bestcoastline, "from the ocedata package.\n")
coastlineWorld <- get(bestcoastline)
} else {
data(list = "coastlineWorld", package = "oce", envir = environment())
oceDebug(debug, "The ocedata package is not available, so using", bestcoastline, "from oce\n")
coastlineWorld <- get("coastlineWorld")
}
}
polygon(coastlineWorld[["longitude"]], coastlineWorld[["latitude"]], col = "tan")
# use lon and lat, if node not given
if (!nodeGiven && longitudeGiven && latitudeGiven) {
closest <- which.min(geodDist(triangles$longitude, triangles$latitude, longitude, latitude))
node <- triangles$triangle[closest]
}
if (nodeGiven && node < 0 && interactive()) {
point <- locator(1)
node <- which.min(geodDist(triangles$longitude, triangles$latitude, point$x, point$y))
}
if (missing(node)) {
node <- triangles$number
longitude <- triangles$longitude
latitude <- triangles$latitude
} else {
if (is.finite(node)) {
node <- triangles$triangle[node]
longitude <- triangles$longitude[node]
latitude <- triangles$latitude[node]
points(longitude, latitude, pch = 20, cex = 2, col = "blue")
legend("topleft",
pch = 20, pt.cex = 2, cex = 3 / 4, col = "blue", bg = "white",
legend = sprintf("node %.0f %.3fN %.3fE", node, latitude, longitude)
)
}
}
oceDebug(debug, "} # webtide()\n", sep = "", unindent = 1, style = "bold")
return(invisible(list(node = node, latitude = latitude, longitude = longitude)))
} else {
node <- triangles$triangle
longitude <- triangles$longitude
latitude <- triangles$latitude
oceDebug(debug, "} # webtide()\n", sep = "", unindent = 1, style = "bold")
return(list(node = node, latitude = latitude, longitude = longitude))
}
} else if (action == "predict") {
if (missing(time)) {
time <- seq.POSIXt(from = presentTime(), by = "15 min", length.out = 7 * 4 * 24)
} # Q: what about timezone?
if (missing(node)) {
if (missing(longitude) || missing(latitude)) {
stop("'longitude' and 'latitude' must be given unless 'node' is given")
}
node <- which.min(geodDist(triangles$longitude, triangles$latitude, longitude, latitude))
} else {
latitude <- triangles$latitude[node]
longitude <- triangles$longitude[node]
}
oceDebug(debug, latitude, "N ", longitude, "E -- use node ", node,
" at ", triangles$latitude[node], "N ", triangles$longitude[node], "E\n",
sep = ""
)
constituentse <- dir(path = path, pattern = "*.s2c")
constituentsuv <- dir(path = path, pattern = "*.v2c")
nconstituents <- length(constituentse)
period <- ampe <- phasee <- ampu <- phaseu <- ampv <- phasev <- vector("numeric", length(nconstituents))
data("tidedata", package = "oce", envir = environment())
tidedata <- get("tidedata") # , pos=globalenv())
for (i in 1:nconstituents) {
twoLetter <- substr(constituentse[i], 1, 2)
C <- which(twoLetter == tidedata$const$name)
period[i] <- 1 / tidedata$const$freq[C]
# Elevation file contains one entry per node, starting with e.g.:
# tri
# period 23.934470 (hours) first harmonic
# 260.000000 (days)
# 1 0.191244 223.820954
# 2 0.188446 223.141200
coneFile <- paste(path, constituentse[i], sep = "/")
cone <- read.table(coneFile, skip = 3, encoding = "latin1")[node, ]
ampe[i] <- cone[[2]]
phasee[i] <- cone[[3]]
conuvFile <- paste(path, constituentsuv[i], sep = "/")
conuv <- read.table(conuvFile, skip = 3, encoding = "latin1")[node, ]
ampu[i] <- conuv[[2]]
phaseu[i] <- conuv[[3]]
ampv[i] <- conuv[[4]]
phasev[i] <- conuv[[5]]
oceDebug(debug, coneFile, sprintf("%s ", twoLetter),
sprintf("%4.2fh", period[i]),
sprintf(" (node %d) ", node),
sprintf("%4.4fm ", ampe[i]), sprintf("%3.3fdeg", phasee[i]), "\n",
sep = ""
)
}
elevation <- u <- v <- rep(0, length(time))
# NOTE: tref is the *central time* for tidem()
tRef <- ISOdate(1899, 12, 31, 12, 0, 0, tz = "UTC")
h <- (as.numeric(time) - as.numeric(tRef)) / 3600
tRef <- 3600 * round(mean(as.numeric(time)) / 3600)
for (i in 1:nconstituents) {
twoLetter <- substr(constituentse[i], 1, 2)
C <- which(twoLetter == tidedata$const$name)
vuf <- tidemVuf(tRef, j = C, latitude = latitude)
phaseOffset <- (vuf$u + vuf$v) * 360
# NOTE: phase is *subtracted* here, but *added* in tidem()
elevation <- elevation + ampe[i] * cos((360 * h / period[i] - phasee[i] + phaseOffset) * rpd)
#> lines(time, elevation, col=i,lwd=3) # Debug
u <- u + ampu[i] * cos((360 * h / period[i] - phaseu[i] + phaseOffset) * rpd)
v <- v + ampv[i] * cos((360 * h / period[i] - phasev[i] + phaseOffset) * rpd)
oceDebug(debug, sprintf("%s ", twoLetter),
sprintf("%4.2fh ", period[i]),
sprintf("%4.4fm ", ampe[i]), sprintf("%3.3fdeg", phasee[i]), "\n",
sep = ""
)
}
if (plot) {
par(mfrow = c(3, 1))
oce.plot.ts(time, elevation,
type = "l", xlab = "", ylab = resizableLabel("elevation"),
main = sprintf("node %.0f %.3fN %.3fE", node, latitude, longitude),
tformat = tformat
)
abline(h = 0, lty = "dotted", col = "gray")
oce.plot.ts(time, u,
type = "l", xlab = "", ylab = resizableLabel("u"),
drawTimeRange = FALSE, tformat = tformat
)
abline(h = 0, lty = "dotted", col = "gray")
oce.plot.ts(time, v,
type = "l", xlab = "", ylab = resizableLabel("v"),
drawTimeRange = FALSE, tformat = tformat
)
abline(h = 0, lty = "dotted", col = "gray")
oceDebug(debug, "} # webtide()\n", sep = "", unindent = 1, style = "bold")
return(invisible(list(
time = time, elevation = elevation, u = u, v = v,
node = node, basedir = basedir, region = region
)))
} else {
oceDebug(debug, "} # webtide()\n", sep = "", unindent = 1, style = "bold")
return(list(
time = time, elevation = elevation, u = u, v = v,
node = node, basedir = basedir, region = region
))
}
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.