# vim:textwidth=80:expandtab:shiftwidth=4:softtabstop=4
#' Class to Store Sealevel Data
#'
#' This class stores sealevel data, e.g. from a tide gauge.
#'
#' @templateVar class sealevel
#'
#' @templateVar dataExample The key items stored in this slot are `time` and `elevation`.
#'
# nolint start (long lines)
#' @templateVar metadataExample An example of the former might be the location at which a `sealevel` measurement was made, stored in `longitude` and `latitude`, and of the latter might be `filename`, the name of the data source.
# nolint end (long lines)
#'
#' @template slot_summary
#'
#' @template slot_put
#'
#' @template slot_get
#'
#' @author Dan Kelley
#'
#' @family classes provided by oce
#' @family things related to sealevel data
setClass("sealevel", contains = "oce")
#' Sample sealevel Data (Halifax Harbour)
#'
#' This sample sea-level dataset is the 2003 record from Halifax Harbour in
#' Nova Scotia, Canada. For reasons that are not mentioned on the data archive
#' website, the record ends on the 8th of October.
#'
#' See [predict.tidem()] for an example that reveals the storm surge
#' that resulted from Hurricane Juan, in this year.
#'
#' @name sealevel
#'
#' @docType data
#'
#' @author Dan Kelley
#'
#' @source The data were created as \preformatted{ sealevel <-
#' read.oce("490-01-JAN-2003_slev.csv") sealevel <- oce.edit(sealevel,
#' "longitude", -sealevel[["longitude"]], reason="Fix longitude hemisphere") }
#' where the csv file was downloaded from reference 1. Note the correction of longitude
#' sign, which is required because the data file has no indication that this is
#' the western hemisphere.
#'
#' @references
#' 1. Fisheries and Oceans Canada \code{http://www.meds-sdmm.dfo-mpo.gc.ca/isdm-gdsi/index-eng.html}
#'
#' @family datasets provided with oce
#' @family things related to sealevel data
NULL
#' Sample sealevel Data (Tuktoyaktuk)
#'
#' This sea-level dataset is provided with in Appendix 7.2 of Foreman (1977)
#' and also with the `T_TIDE` package (Pawlowicz et al., 2002). It results
#' from measurements made in 1975 at Tuktoyaktuk, Northwest Territories,
#' Canada.
#'
#' The data set contains 1584 points, some of which have NA for sea-level
#' height.
#'
#' Although Foreman's Appendix 7.2 states that times are in Mountain standard
#' time, the timezone is set to `UTC` in the present case, so that the
#' results will be similar to those he provides in his Appendix 7.3.
#'
#' @name sealevelTuktoyaktuk
#'
#' @docType data
#'
#' @references Foreman, M. G. G., 1977. Manual for tidal heights analysis and
#' prediction. Pacific Marine Science Report 77-10, Institute of Ocean
#' Sciences, Patricia Bay, Sidney, BC, 58pp.
#'
#' 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 were based on the `T_TIDE` dataset, which in turn
#' seems to be based on Appendix 7.2 of Foreman (1977). Minor editing was on
#' file format, and then the `sealevelTuktoyaktuk` object was created
#' using [as.sealevel()].
#'
## @examples
## \donttest{
## library(oce)
## data(sealevelTuktoyaktuk)
## time <- sealevelTuktoyaktuk[["time"]]
## elevation <- sealevelTuktoyaktuk[["elevation"]]
## oce.plot.ts(time, elevation, type="l", ylab="Height [m]", ylim=c(-2, 6))
## legend("topleft", legend=c("Tuktoyaktuk (1975)","Detided"),
## col=c("black","red"),lwd=1)
## tide <- tidem(sealevelTuktoyaktuk)
## detided <- elevation - predict(tide)
## lines(time, detided, col="red")
## }
#'
#' @section Historical note:
#' Until Jan 6, 2018, the time in this dataset had been increased
#' by 7 hours. However, this alteration was removed on this date,
#' to make for simpler comparison of amplitude and phase output with
#' the results obtained by Foreman (1977) and Pawlowicz et al. (2002).
#'
#' @family datasets provided with oce
#' @family things related to sealevel data
NULL
setMethod(
f = "initialize",
signature = "sealevel",
definition = function(.Object, elevation, time, ...) {
.Object <- callNextMethod(.Object, ...)
if (!missing(elevation)) {
.Object@data$elevation <- elevation
}
if (!missing(time)) {
.Object@data$time <- time
}
.Object@processingLog$time <- presentTime()
.Object@processingLog$value <- "create 'sealevel' object"
return(.Object)
}
)
#' Summarize a sealevel Object
#'
#' Summarizes some of the data in a sealevel object.
#'
#' @param object A [sealevel-class] object.
#'
#' @param \dots further arguments passed to or from other methods.
#'
#' @return A matrix containing statistics of the elements of the `data`
#' slot.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' data(sealevel)
#' summary(sealevel)
#'
#' @family things related to sealevel data
setMethod(
f = "summary",
signature = "sealevel",
definition = function(object, ...) {
cat("Sealevel Summary\n----------------\n\n")
showMetadataItem(object, "stationNumber", "number: ")
showMetadataItem(object, "version", "version: ")
showMetadataItem(object, "stationName", "name: ")
showMetadataItem(object, "region", "region: ")
showMetadataItem(object, "deltat", "sampling delta-t: ", postlabel = " hour")
cat("* Location: ", latlonFormat(object@metadata$latitude,
object@metadata$longitude,
digits = 5
), "\n")
showMetadataItem(object, "year", "year: ")
ndata <- length(object@data$elevation)
cat("* number of observations: ", ndata, "\n")
cat("* \" non-missing: ", sum(!is.na(object@data$elevation)), "\n")
invisible(callNextMethod()) # summary
}
)
#' @title Subset a sealevel Object
#'
#' @description
#' This function is somewhat analogous to [subset.data.frame()], but
#' subsetting is only permitted by time.
#'
#' @param x a [sealevel-class] object.
#'
#' @param subset a condition to be applied to the `data` portion of
#' `x`.
#'
#' @param \dots ignored.
#'
#' @return A new `sealevel` object.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' data(sealevel)
#' plot(sealevel)
#' plot(subset(sealevel, time < mean(range(sealevel[["time"]]))))
#'
#' @family things related to sealevel data
#' @family functions that subset oce objects
setMethod(
f = "subset",
signature = "sealevel",
definition = function(x, subset, ...) {
res <- new("sealevel")
res@metadata <- x@metadata
res@processingLog <- x@processingLog
for (i in seq_along(x@data)) {
r <- eval(expr = substitute(expr = subset, env = environment()), envir = x@data, enclos = parent.frame(2))
r <- r & !is.na(r)
res@data[[i]] <- x@data[[i]][r]
}
names(res@data) <- names(x@data)
subsetString <- paste(deparse(substitute(expr = subset, env = environment())), collapse = " ")
res@processingLog <- processingLogAppend(res@processingLog, paste("subset.sealevel(x, subset=", subsetString, ")", sep = ""))
res
}
)
#' Extract Something From a sealevel Object
#'
#' @param x a [sealevel-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. The `dataDerived`
#' and `metadataDerived` items are each NULL, because
#' no derived values are defined by [sealevel-class] objects.
#'
#' * In many cases, the focus will be on variations of
#' sealevel elevation over time, so it is common to use e.g.
#' `x[["time"]]` and `x[["elevation"]]` to retrieve vectors
#' of these quantities. Another common task is to retrieve
#' the location of the observations, using e.g. `x[["longitude"]]`
#' and `x[["latitude"]]`.
#'
#' @template sub_subTemplate
#'
#' @author Dan Kelley
#'
#' @family things related to sealevel data
setMethod(
f = "[[",
signature(x = "sealevel", i = "ANY", j = "ANY"),
definition = function(x, i, j, ...) {
if (i == "?") {
return(list(
metadata = sort(names(x@metadata)),
metadataDerived = NULL,
data = sort(names(x@data)),
dataDerived = NULL
))
}
callNextMethod() # [[
}
)
#' @title Replace Parts of a sealevel Object
#'
#' @param x a [sealevel-class] object.
#'
#' @template sub_subsetTemplate
#'
#' @family things related to sealevel data
setMethod(
f = "[[<-",
signature(x = "sealevel", i = "ANY", j = "ANY"),
definition = function(x, i, j, ..., value) {
callNextMethod(x = x, i = i, j = j, ... = ..., value = value) # [[<-
}
)
setValidity(
"sealevel",
function(object) {
ndata <- length(object@data)
lengths <- vector("numeric", ndata)
for (i in 1:ndata) {
lengths[i] <- length(object@data[[i]])
}
if (var(lengths) != 0) {
cat("lengths of data elements are unequal\n")
return(FALSE)
} else {
return(TRUE)
}
}
)
#' Coerce Data Into a sealevel Object
#'
#' Coerces a dataset (minimally, a sequence of times and heights) into a
#' sealevel dataset.
#' The arguments are based on the standard data format, as were described in a
#' file formerly available at reference 1.
#'
#' @param elevation a list of sea-level heights in metres, in an hourly
#' sequence.
#'
#' @param time optional list of times, in POSIXct format. If missing, the list
#' will be constructed assuming hourly samples, starting at 0000-01-01
#' 00:00:00.
#'
#' @param header a character string as read from first line of a standard data
#' file.
#'
#' @param stationNumber three-character string giving station number.
#'
#' @param stationVersion single character for version of station.
#'
#' @param stationName the name of station (at most 18 characters).
#'
#' @param region the name of the region or country of station (at most 19
#' characters).
#'
#' @param year the year of observation.
#'
#' @param longitude the longitude in decimal degrees, positive east of
#' Greenwich.
#'
#' @param latitude the latitude in decimal degrees, positive north of the
#' equator.
#'
#' @param GMTOffset offset from GMT, in hours.
#'
#' @param decimationMethod a coded value, with 1 meaning filtered, 2 meaning a
#' simple average of all samples, 3 meaning spot readings, and 4 meaning some
#' other method.
#'
#' @param referenceOffset ?
#'
#' @param referenceCode ?
#'
#' @param deltat optional interval between samples, in hours (as for the
#' [ts()] timeseries function). If this is not provided, and `t`
#' can be understood as a time, then the difference between the first two times
#' is used. If this is not provided, and `t` cannot be understood as a
#' time, then 1 hour is assumed.
#'
#' @return A [sealevel-class] object (for details, see [read.sealevel()]).
#'
#' @author Dan Kelley
#'
#' @seealso The documentation for the [sealevel-class] class explains the
#' structure of sealevel objects, and also outlines the other functions dealing
#' with them.
#'
#' @references `http://ilikai.soest.hawaii.edu/rqds/hourly.fmt` (this link
#' worked for years but failed at least temporarily on December 4, 2016).
#'
#' @examples
#' library(oce)
#'
#' # Construct a year of M2 tide, starting at the default time
#' # 0000-01-01T00:00:00.
#' h <- seq(0, 24 * 365)
#' elevation <- 2.0 * sin(2 * pi * h / 12.4172)
#' sl <- as.sealevel(elevation)
#' summary(sl)
#'
#' # As above, but start at the Y2K time.
#' time <- as.POSIXct("2000-01-01") + h * 3600
#' sl <- as.sealevel(elevation, time)
#' summary(sl)
#' @family things related to sealevel data
as.sealevel <- function(
elevation, time, header = NULL,
stationNumber = NA, stationVersion = NA, stationName = NULL,
region = NULL, year = NA, longitude = NA, latitude = NA, GMTOffset = NA,
decimationMethod = NA, referenceOffset = NA, referenceCode = NA, deltat) {
if (missing(elevation)) {
stop("must supply sealevel height, elevation, in metres")
}
if (inherits(elevation, "POSIXt")) {
stop("elevation must be a numeric vector, not a time vector")
}
res <- new("sealevel")
n <- length(elevation)
if (missing(time)) {
# construct hourly from time "zero"
start <- as.POSIXct("0000-01-01 00:00:00", tz = "UTC")
time <- as.POSIXct(start + seq(0, n - 1, 1) * 3600, tz = "UTC")
if (is.na(GMTOffset)) {
GMTOffset <- 0
} # FIXME: do I want to do this?
} else {
time <- as.POSIXct(time, tz = "UTC")
}
if (missing(deltat)) {
deltat <- as.numeric(difftime(time[2], time[1], units = "hours"))
}
if (is.na(deltat) || deltat <= 0) {
deltat <- 1
}
res@metadata$filename <- ""
res@metadata$header <- header
res@metadata$year <- year
res@metadata$stationNumber <- stationNumber
res@metadata$stationVersion <- stationVersion
res@metadata$stationName <- stationName
res@metadata$region <- region
res@metadata$latitude <- latitude
res@metadata$longitude <- longitude
res@metadata$GMTOffset <- GMTOffset
res@metadata$decimationMethod <- decimationMethod
res@metadata$referenceOffset <- referenceOffset
res@metadata$referenceCode <- referenceCode
res@metadata$units <- list(elevation = list(unit = expression(m), scale = ""))
res@metadata$n <- length(t)
res@metadata$deltat <- deltat
res@data$elevation <- elevation
res@data$time <- time
res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep = "", collapse = ""))
res
}
#' @title Plot a sealevel Object
#'
#' @description
#' Creates a plot for a sea-level dataset, in one of two varieties. Depending
#' on the length of `which`, either a single-panel or multi-panel plot is
#' drawn. If there is just one panel, then the value of `par` used in
#' `plot,sealevel-method` is retained upon exit, making it convenient to add to
#' the plot. For multi-panel plots, `par` is returned to the value it had
#' before the call.
#'
#' @param x a [sealevel-class] object.
#'
#' @param which a numerical or string vector indicating desired plot types,
#' with possibilities 1 or `"all"` for a time-series of all the elevations, 2 or
#' `"month"` for a time-series of just the first month, 3 or
#' `"spectrum"` for a power spectrum (truncated to frequencies below 0.1
#' cycles per hour, or 4 or `"cumulativespectrum"` for a cumulative
#' integral of the power spectrum.
#'
#' @param drawTimeRange boolean that applies to panels with time as the
#' horizontal axis, indicating whether to draw the time range in the top-left
#' margin of the plot.
#'
#' @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 marginsAsImage logical value indicating whether to put a
#' wide margin to the right of time-series plots, matching the space
#' used up by a palette in an [imagep()] plot.
#'
#' @param grid logical value, indicating whether to draw a grid with
#' [grid()].
#'
#' @param xlim,ylim optional limits for axes. If not supplied,
#' reasonable choices will be made
#'
#' @param xaxs,yaxs axis-limit parameters, as for standard graphics.
#' The default is to make the time axis extend to the edges of
#' the box, but to make the y axis have some space above and below
#' the range of the data.
#'
#' @param debug a flag that turns on debugging, if it exceeds 0.
#'
#' @param \dots optional arguments passed to plotting functions.
#'
#' @return None.
#'
#' @author Dan Kelley
#'
#' @seealso The documentation for the [sealevel-class] class explains the
#' structure of sealevel objects, and also outlines the other functions dealing
#' with them.
#'
#' @section Historical Note:
#' Until 2020-02-06, sea-level plots had the mean value removed, and indicated
#' with a tick mark and margin note on the right-hand side of the plot.
#' This behaviour was confusing. The change did not go through the usual
#' deprecation process, because the margin-note behaviour had not
#' been documented.
#'
#' @references The example refers to Hurricane Juan, which caused a great deal
#' of damage to Halifax in 2003. Since this was in the era of the digital
#' photo, a casual web search will uncover some spectacular images of damage,
#' from both wind and storm surge.
#' Landfall, within 30km of this sealevel
#' gauge, was between 00:10 and 00:20 Halifax local time on Monday, Sept 29,
#' 2003.
#'
#' @examples
#' library(oce)
#' data(sealevel)
#' # local Halifax time is UTC + 4h
#' juan <- as.POSIXct("2003-09-29 00:15:00", tz = "UTC") + 4 * 3600
#' plot(sealevel, which = 1, xlim = juan + 86400 * c(-7, 7))
#' abline(v = juan, col = "red")
#'
#' @family functions that plot oce data
#' @family things related to sealevel data
#'
#' @aliases plot.sealevel
setMethod(
f = "plot",
signature = signature("sealevel"),
definition = function(x, which = 1:3,
drawTimeRange = getOption("oceDrawTimeRange"),
mgp = getOption("oceMgp"),
mar = c(mgp[1] + 0.5, mgp[1] + 1.5, mgp[2] + 1, mgp[2] + 3 / 4),
marginsAsImage = FALSE,
grid = TRUE,
xlim, ylim, xaxs = "i", yaxs = "r",
debug = getOption("oceDebug"),
...) {
oceDebug(debug, "plot.sealevel(..., mar=c(", paste(mar, collapse = ", "), "), ...) START\n", sep = "", unindent = 1)
xlimGiven <- !missing(xlim)
ylimGiven <- !missing(ylim)
titlePlot <- function(x) {
title <- ""
if (!is.null(x@metadata$stationNumber) || !is.null(x@metadata$stationName) || !is.null(x@metadata$region)) {
# "Station"
# title <- paste(title, gettext("Station ", domain = "R-oce"),
title <- paste(
if (!is.na(x@metadata$stationNumber)) x@metadata$stationNumber else "",
if (!is.null(x@metadata$stationName)) x@metadata$stationName else "",
if (!is.null(x@metadata$region)) x@metadata$region else ""
)
oceDebug(debug, paste0("stage 1. title=\"", title, "\"\n"))
}
if (!is.na(x@metadata$latitude) && !is.na(x@metadata$longitude)) {
title <- paste(title, latlonFormat(x@metadata$latitude, x@metadata$longitude), sep = "")
}
oceDebug(debug, paste0("stage 2. title=\"", title, "\"\n"))
title <- trimws(title)
oceDebug(debug, paste0("stage 3. title=\"", title, "\"\n"))
if (nchar(title) > 0) {
title <- paste0(gettext("Station ", domain = "R-oce"), title)
oceDebug(debug, paste0("stage 4. title=\"", title, "\"\n"))
mtext(side = 3, title, adj = 1, cex = 2 / 3)
}
}
drawConstituent <- function(frequency = 0.0805114007, label = "M2", col = "darkred", side = 1) {
abline(v = frequency, col = col)
mtext(label, side = side, at = frequency, col = col, cex = 3 / 4 * par("cex"))
}
drawConstituents <- function() {
drawConstituent(0.0387306544, "O1", side = 1)
# draw.constituent(0.0416666721, "S1", side=3)
drawConstituent(0.0417807462, "K1", side = 3)
drawConstituent(0.0789992488, "N2", side = 1)
drawConstituent(0.0805114007, "M2", side = 3)
drawConstituent(0.0833333333, "S2", side = 1)
}
if (!inherits(x, "sealevel")) {
stop("method is only for objects of class '", "sealevel", "'")
}
opar <- par(no.readonly = TRUE)
par(mgp = mgp, mar = mar)
lw <- length(which)
if (marginsAsImage) {
scale <- 0.7
w <- (1.5 + par("mgp")[2]) * par("csi") * scale * 2.54 + 0.5
if (lw > 1) {
layout(matrix(1:(2 * lw), nrow = lw, byrow = TRUE), widths = rep(c(1, lcm(w)), lw))
}
} else {
if (lw > 1) {
layout(cbind(1:lw))
}
}
if (lw > 1) on.exit(par(opar))
# tidal constituents (in cpd):
# http://www.soest.hawaii.edu/oceanography/dluther/HOME/Tables/Kaw.htm
num.NA <- sum(is.na(x@data$elevation))
par(mgp = mgp)
# par(mar=c(mgp[1],mgp[1]+2.5,mgp[2]+0.5,mgp[2]+1))
par(mar = mar)
# ?used? n <- length(x@data$elevation) # do not trust value in metadata
oceDebug(debug, "which:", which, "\n")
which2 <- oce.pmatch(which, list(all = 1, month = 2, spectrum = 3, cumulativespectrum = 4))
oceDebug(debug, "which2:", which2, "\n")
for (w in seq_along(which2)) {
oceDebug(debug, "plotting for code which2[", w, "] = ", which2[w], "\n", sep = "")
if (which2[w] == 1) {
xlim <- if (xlimGiven) xlim else (range(x@data$time, na.rm = TRUE))
ylim <- if (ylimGiven) ylim else (range(x@data$elevation, na.rm = TRUE))
plot(x@data$time, x@data$elevation,
xlab = "",
ylab = resizableLabel("elevation"),
type = "l",
lwd = 0.5, axes = FALSE,
xlim = xlim, ylim = ylim, xaxs = xaxs, yaxs = yaxs,
...
)
tics <- oce.axis.POSIXct(1, x@data$time, drawTimeRange = drawTimeRange, cex.axis = 1, debug = debug - 1)
box()
titlePlot(x)
yax <- axis(2)
if (grid) {
abline(h = yax, col = "darkgray", lty = "dotted")
abline(v = tics, col = "darkgray", lty = "dotted")
}
} else if (which2[w] == 2) {
# sample month
from <- trunc(x@data$time[1], "day")
to <- from + 28 * 86400 # 28 days
look <- from <= x@data$time & x@data$time <= to
xx <- x
for (i in seq_along(x@data)) {
xx@data[[i]] <- x@data[[i]][look]
}
if (any(is.finite(xx@data$elevation))) {
xlim <- if (xlimGiven) xlim else (range(x@data$time, na.rm = TRUE))
ylim <- if (ylimGiven) ylim else (range(x@data$elevation, na.rm = TRUE))
atWeek <- seq(from = from, to = to, by = "week")
atDay <- seq(from = from, to = to, by = "day")
plot(xx@data$time, xx@data$elevation,
xlab = "",
ylab = resizableLabel("elevation"),
type = "l", xaxs = "i",
xlim = xlim, ylim = ylim, xaxs = xaxs, yaxs = yaxs,
axes = FALSE
)
oce.axis.POSIXct(1, xx@data$time, drawTimeRange = drawTimeRange, cex.axis = 1, debug = debug - 1)
yax <- axis(2)
box()
if (grid) {
abline(h = yax, col = "lightgray", lty = "dotted")
abline(v = atWeek, col = "darkgray", lty = "dotted")
abline(v = atDay, col = "lightgray", lty = "dotted")
}
} else {
plot(0:1, 0:1, type = "n", xlab = "", ylab = "", axes = FALSE)
box()
text(0.5, 0.5, "Cannot show first month, since all data are NA then")
}
} else if (which2[w] == 3) { # "spectrum"
if (num.NA == 0) {
Elevation <- ts(x@data$elevation, start = 1, deltat = x@metadata$deltat)
s <- spectrum(Elevation - mean(Elevation), plot = FALSE, log = "y", demean = TRUE, detrend = TRUE)
par(mar = c(mgp[1] + 1.25, mgp[1] + 1.5, mgp[2] + 0.25, mgp[2] + 3 / 4))
xlim <- if (xlimGiven) xlim else c(0, 0.1)
ylim <- if (ylimGiven) {
ylim
} else {
range(subset(
s$spec,
xlim[1] <= s$freq & s$freq <= xlim[2]
), na.rm = TRUE)
}
plot(s$freq, s$spec,
xlab = resizableLabel("frequency cph"),
ylab = resizableLabel("spectral density m2/cph"),
# [m^2/cph]",
xlim = xlim, ylim = ylim, xaxs = xaxs, yaxs = yaxs,
type = "l", log = "y",
...
)
if (grid) {
grid(col = "darkgray", lty = "dotted")
}
drawConstituents()
} else {
plot(0:1, 0:1, type = "n", xlab = "", ylab = "", axes = FALSE)
box()
text(0.5, 0.5, "Some elevations are NA, so cannot calculate the spectrum")
}
} else if (which2[w] == 4) { # "cumulativespectrum"
if (num.NA == 0) {
# ?used? n <- length(x@data$elevation)
Elevation <- ts(x@data$elevation, start = 1, deltat = x@metadata$deltat)
s <- spectrum(Elevation - mean(Elevation), plot = FALSE, log = "y", demean = TRUE, detrend = TRUE)
nCumSpec <- length(s$spec)
cumSpec <- sqrt(cumsum(s$spec) / nCumSpec)
par(mar = c(mgp[1] + 1.25, mgp[1] + 2.5, mgp[2] + 0.25, mgp[2] + 0.25))
xlim <- if (xlimGiven) xlim else c(0, 0.1)
ylim <- if (ylimGiven) ylim else range(cumSpec, na.rm = TRUE)
plot(s$freq, cumSpec,
type = "l",
xlab = resizableLabel("frequency cph"),
ylab = expression(paste(integral(Gamma, 0, f), " df [m]")),
xlim = xlim, ylim = ylim, xaxs = xaxs, yaxs = yaxs
)
if (grid) {
grid(col = "darkgray", lty = "dotted")
}
drawConstituents()
} else {
warning("cannot draw sealevel spectum, because the series contains missing values")
}
} else {
stop("unrecognized value of which: ", which[w])
}
if (marginsAsImage) {
# blank plot, to get axis length same as for images
omar <- par("mar")
par(mar = c(mar[1], 1 / 4, mgp[2] + 1 / 2, mgp[2] + 1))
plot(1:2, 1:2, type = "n", axes = FALSE, xlab = "", ylab = "")
par(mar = omar)
}
}
oceDebug(debug, "END plot.sealevel()\n", unindent = 1)
invisible(NULL)
}
)
#' Read a sealevel File
#'
#' Read a data file holding sea level data. BUG: the time vector assumes GMT,
#' regardless of the GMT.offset value.
#'
#' This function starts by scanning the first line of the file, from which it
#' determines whether the file is in one of two known formats: type 1, the
#' format used at the Hawaii archive centre, and type 2, the
#' comma-separated-value format used by the Marine Environmental Data Service.
#' The file type is inferred by examination of its first line. If that contains
#' the string `Station_Name` the file is of type 2. If
#' the file is in neither of these formats, the user might wish to scan it
#' directly, and then to use [as.sealevel()] to create a
#' `sealevel` object.
#' The Hawaii archive site at
#' `http://ilikai.soest.hawaii.edu/uhslc/datai.html` at one time provided a graphical
#' interface for downloading sealevel data in Type 1, with format that was once
#' described at `http://ilikai.soest.hawaii.edu/rqds/hourly.fmt` (although that link
#' was observed to no longer work, on December 4, 2016).
#' Examination of data retrieved from what seems to be a replacement Hawaii server
#' (https://uhslc.soest.hawaii.edu/data/?rq) in September 2019 indicated that the
#' format had been changed to what is called Type 3 by `read.sealevel`.
#' Web searches did not uncover documentation on this format, so the
#' decoding scheme was developed solely through examination of
#' data files, which means that it might be not be correct.
#' The MEDS repository (\code{http://www.isdm-gdsi.gc.ca/isdm-gdsi/index-eng.html})
#' provides Type 2 data.
#'
#' @param file a connection or a character string giving the name of the file
#' to load. See Details for the types of files that are recognized.
#'
#' @param tz time zone. The default value, `oceTz`, is set to `UTC`
#' at setup. (If a time zone is present in the file header, this will
#' supercede the value given here.)
#'
#' @template encodingTemplate
#'
#' @param processingLog if provided, the action item to be stored in the log.
#' (Typically only provided for internal calls; the default that it provides is
#' better for normal calls by a user.)
#' @template debugTemplate
#'
#' @return A [sealevel-class] object.
#'
#' @author Dan Kelley
#'
#' @family things related to sealevel data
read.sealevel <- function(
file, tz = getOption("oceTz"), encoding = "latin1",
processingLog, debug = getOption("oceDebug")) {
if (missing(file)) {
stop("must supply 'file'")
}
if (is.character(file)) {
if (!file.exists(file)) {
stop("cannot find file \"", file, "\"")
}
if (0L == file.info(file)$size) {
stop("empty file \"", file, "\"")
}
}
oceDebug(debug, "read.sealevel(file=\"", file, "\", ...) START\n", sep = "", unindent = 1)
filename <- "?"
if (is.character(file)) {
filename <- fullFilename(file)
file <- file(file, "r", encoding = encoding)
on.exit(close(file))
}
if (!inherits(file, "connection")) {
stop("argument `file' must be a character string or connection")
}
if (!isOpen(file)) {
filename <- "(connection)"
open(file, "r", encoding = encoding)
on.exit(close(file))
}
firstLine <- readLines(file, n = 1)
header <- firstLine
oceDebug(debug, "header (first line in file): '", header, "'\n", sep = "")
pushBack(firstLine, file)
stationNumber <- NA
stationVersion <- NA
stationName <- NULL
region <- NULL
year <- NA
latitude <- NA
longitude <- NA
GMTOffset <- NA
decimationMethod <- NA
referenceOffset <- NA
referenceCode <- NA
res <- new("sealevel")
if (substr(firstLine, 1, 12) == "Station_Name") {
oceDebug(debug, "File is of format 1 (e.g. as in MEDS archives)\n")
# Station_Name,HALIFAX
# Station_Number,490
# Latitude_Decimal_Degrees,44.666667
# Longitude_Decimal_Degrees,63.583333
# Datum,CD
# Time_Zone,AST
# SLEV=Observed Water Level
# Obs_date,SLEV
# 01/01/2001 12:00 AM,1.82,
headerLength <- 8
header <- readLines(file, n = headerLength)
if (debug > 0) {
print(header)
}
stationName <- strsplit(header[1], ",")[[1]][2]
stationNumber <- as.numeric(strsplit(header[2], ",")[[1]][2])
latitude <- as.numeric(strsplit(header[3], ",")[[1]][2])
longitude <- as.numeric(strsplit(header[4], ",")[[1]][2])
tz <- strsplit(header[6], ",")[[1]][2] # needed for get GMT offset
GMTOffset <- GMTOffsetFromTz(tz)
oceDebug(debug, "about to read data\n")
x <- read.csv(file, header = FALSE, stringsAsFactors = FALSE, encoding = encoding)
oceDebug(debug, "... finished reading data\n")
if (length(grep("[0-9]{4}/", x$V1[1])) > 0) {
oceDebug(debug, "Date format is year/month/day hour:min with hour in range 1:24\n")
time <- strptime(as.character(x$V1), "%Y/%m/%d %H:%M", "UTC") + 3600 * GMTOffset
} else {
oceDebug(debug, "Date format is day/month/year hour:min AMPM with hour in range 1:12 and AMPM indicating whether day or night\n")
time <- strptime(as.character(x$V1), "%d/%m/%Y %I:%M %p", "UTC") + 3600 * GMTOffset
}
elevation <- as.numeric(x$V2)
oceDebug(
debug, "tz=", tz, "so GMTOffset=", GMTOffset, "\n",
"first pass has time string:", as.character(x$V1)[1], "\n",
"first pass has time start:", format(time[1]), " ", attr(time[1], "tzone"), "\n"
)
year <- as.POSIXlt(time[1])$year + 1900
} else {
oceDebug(debug, "File is of type 2 or 3\n")
d <- readLines(file)
n <- length(d)
header <- d[1]
if (grepl("LAT=", header) && grepl("LONG=", header) && grepl("TIMEZONE", header)) {
# URL
# http://uhslc.soest.hawaii.edu/woce/h275.dat
# is a sample file, which starts as below (with quote marks added):
# '275HALIFAX 1895 LAT=44 40.0N LONG=063 35.0W TIMEZONE=GMT '
# '275HALIFAX 189501011 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999'
# '275HALIFAX 189501012 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999'
oceDebug(debug, "type 3 (format inferred/guessed from e.g. http://uhslc.soest.hawaii.edu/woce/h275.dat)\n")
stationNumber <- strtrim(header, 3)
oceDebug(debug, " stationNumber=\"", stationNumber, "\"\n", sep = "")
longitudeString <- gsub("^.*LONG=([ 0-9.]*[EWew]).*$", "\\1", header)
latitudeString <- gsub("^.*LAT=([ 0-9.]*[NSns]).*$", "\\1", header)
oceDebug(debug, " longitudeString=\"", longitudeString, "\"\n", sep = "")
oceDebug(debug, " latitudeString=\"", latitudeString, "\"\n", sep = "")
longitudeSplit <- strsplit(longitudeString, split = "[ \t]+")[[1]]
longitudeDegree <- longitudeSplit[1]
longitudeMinute <- longitudeSplit[2]
oceDebug(debug, " longitudeDegree=\"", longitudeDegree, "\"\n", sep = "")
oceDebug(debug, " longitudeMinute=\"", longitudeMinute, "\"\n", sep = "")
longitudeSign <- if (grepl("[wW]", longitudeMinute)) -1 else 1
oceDebug(debug, " longitudeSign=", longitudeSign, "\n")
longitudeMinute <- gsub("[EWew]", "", longitudeMinute)
oceDebug(debug, " longitudeMinute=\"", longitudeMinute, "\" after removing EW suffix\n", sep = "")
longitude <- longitudeSign * (as.numeric(longitudeDegree) + as.numeric(longitudeMinute) / 60)
oceDebug(debug, " longitude=", longitude, "\n")
latitudeSplit <- strsplit(latitudeString, split = "[ \t]+")[[1]]
latitudeDegree <- latitudeSplit[1]
latitudeMinute <- latitudeSplit[2]
latitudeSign <- if (grepl("[sS]", latitudeMinute)) -1 else 1
oceDebug(debug, " latitudeSign=", latitudeSign, "\n")
latitudeMinute <- gsub("[SNsn]", "", latitudeMinute)
oceDebug(debug, " latitudeMinute=\"", latitudeMinute, "\" after removing NS suffix\n", sep = "")
latitude <- latitudeSign * (as.numeric(latitudeDegree) + as.numeric(latitudeMinute) / 60)
oceDebug(debug, " latitude=", latitude, "\n")
# Remove interspersed year boundaries (which look like the first line).
d2 <- d[!grepl("LAT=.*LONG=", d)]
start <- 1 + which(strsplit(header, "")[[1]] == " ")[1]
d3 <- substr(d2, start, 1000L)
# Fix problem where the month is sometimes e.g. " 1" instead of "01"
d4 <- gsub("^([1-9][0-9]{3}) ", "\\10", d3)
# Fix problem where the day is sometimes e.g. " 1" instead of "01"
d5 <- gsub("^([1-9][0-9]{3}[0-9]{2}) ", "\\10", d4)
n <- length(d5)
# Now we have as below. But the second block sometimes has " " for "0", so we
# need to fix that.
# 275HALIFAX 189501011 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
twelve <- seq(1, 12, 1)
elevation <- rep(NA, 12 * n)
time <- rep(NA, 12 * n)
lastDayPortion <- NULL # value defined at i=1, checked at i>2, so the initial value is immaterial
for (i in 1:n) {
sp <- strsplit(d5[i], "[ ]+")[[1]]
target.index <- 12 * (i - 1) + twelve
if (length(sp) != 13) {
stop("cannot parse tokens on line \"", d2[i], "\"\n", sep = "")
}
elevation[target.index] <- as.numeric(sp[2:13])
dayPortion <- as.numeric(substr(sp[1], 9, 9))
if (i == 1) {
startDay <- as.POSIXct(strptime(paste(substr(sp[1], 1, 8), "00:00:00"), "%Y%m%d"), tz = tz)
oceDebug(debug, " startDay=", startDay, "\n")
} else {
if (dayPortion == 1) {
if (i > 2 && lastDayPortion != 2) {
stop("non-alternating day portions on data line ", i)
}
} else if (dayPortion == 2) {
if (i > 2 && lastDayPortion != 1) {
stop("non-alternating day portions on data line ", i)
}
} else {
stop("day portion is ", dayPortion, " but must be 1 or 2, on data line", i)
}
}
lastDayPortion <- dayPortion
time[target.index] <- as.POSIXct(sp[1], format = "%Y%m%d", tz = "UTC") + 3600 * (seq(0, 11) + 12 * (dayPortion - 1))
}
elevation[elevation == 9999] <- NA
elevation <- elevation / 1000 # convert mm to m
time <- numberAsPOSIXct(time, tz = "UTC") # guess on timezone
} else {
oceDebug(debug, "type 2 (an old Hawaii format, inferred from documentation)\n")
stationNumber <- substr(header, 1, 3)
stationVersion <- substr(header, 4, 4)
stationName <- substr(header, 6, 23)
stationName <- sub("[ ]*$", "", stationName)
region <- substr(header, 25, 43)
region <- sub("[ ]*$", "", region)
year <- substr(header, 45, 48)
latitudeStr <- substr(header, 50, 55) # degrees,minutes,tenths,hemisphere
latitude <- as.numeric(substr(latitudeStr, 1, 2)) + (as.numeric(substr(latitudeStr, 3, 5))) / 600
if (tolower(substr(latitudeStr, 6, 6)) == "s") latitude <- -latitude
longitudeStr <- substr(header, 57, 63) # degrees,minutes,tenths,hemisphere
longitude <- as.numeric(substr(longitudeStr, 1, 3)) + (as.numeric(substr(longitudeStr, 4, 6))) / 600
if (tolower(substr(longitudeStr, 7, 7)) == "w") longitude <- -longitude
GMTOffset <- substr(header, 65, 68) # hours,tenths (East is +ve)
oceDebug(debug, "GMTOffset=", GMTOffset, "\n")
decimationMethod <- substr(header, 70, 70) # 1=filtered 2=average 3=spot readings 4=other
referenceOffset <- substr(header, 72, 76) # add to values
referenceCode <- substr(header, 77, 77) # add to values
units <- substr(header, 79, 80)
oceDebug(debug, "units=", units, "\n")
if (nchar(units) == 0) {
warning("no units can be inferred from the file, so assuming \"mm\"")
} else {
if (units != "mm" && units != "MM") {
stop("require units to be \"mm\" or \"MM\", not \"", units, "\"")
}
}
elevation <- array(NA_real_, 12 * (n - 1))
twelve <- seq(1, 12, 1)
lastDayPortion <- -1 # ignored; prevents undefined warning in code analysis
for (i in 2:n) {
sp <- strsplit(d[i], "[ ]+")[[1]]
target.index <- 12 * (i - 2) + twelve
elevation[target.index] <- as.numeric(sp[4:15])
dayPortion <- as.numeric(substr(sp[3], 9, 9))
if (i == 2) {
startDay <- as.POSIXct(strptime(paste(substr(sp[3], 1, 8), "00:00:00"), "%Y%m%d"), tz = tz)
} else {
if (dayPortion == 1) {
if (i > 2 && lastDayPortion != 2) {
stop("non-alternating day portions on data line ", i)
}
} else if (dayPortion == 2) {
if (i > 2 && lastDayPortion != 1) {
stop("non-alternating day portions on data line ", i)
}
} else {
stop("day portion is ", dayPortion, " but must be 1 or 2, on data line", i)
}
}
lastDayPortion <- dayPortion
}
time <- as.POSIXct(startDay + 3600 * (seq(0, 12 * (n - 1) - 1)), tz = tz)
elevation[elevation == 9999] <- NA
if (tolower(units) == "mm") {
elevation <- elevation / 1000
} else {
stop("require units to be MM")
}
}
}
res@metadata$filename <- filename
res@metadata$header <- header
res@metadata$year <- year
res@metadata$stationNumber <- stationNumber
res@metadata$stationVersion <- stationVersion
res@metadata$stationName <- stationName
res@metadata$region <- region
res@metadata$latitude <- latitude
res@metadata$longitude <- longitude
res@metadata$GMTOffset <- GMTOffset
res@metadata$decimationMethod <- decimationMethod
res@metadata$referenceOffset <- referenceOffset
res@metadata$referenceCode <- referenceCode
res@metadata$units <- list(elevation = list(unit = expression(m), scale = ""))
res@metadata$n <- length(time)
# deltat is in hours
res@metadata$deltat <- if (res@metadata$n > 1) (as.numeric(time[2]) - as.numeric(time[1])) / 3600 else 0
res@data$elevation <- elevation
res@data$time <- time
res@processingLog <- processingLogAppend(
res@processingLog,
paste("read.sealevel(file=\"", filename, "\", tz=\"", tz, "\")", sep = "", collapse = "")
)
oceDebug(debug, "END read.sealevel()\n", unindent = 1)
res
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.