# vim:textwidth=80:expandtab:shiftwidth=4:softtabstop=4
#' A Package for Oceanographic Analysis
#'
#' @description
#' The oce package provides functions for working with
#' Oceanographic data, for calculations that are specific
#' to Oceanography, and for producing graphics that
#' match the conventions of the field.
#'
#' @details
#'
#' Over a dozen specialized data types are handled by oce,
#' with generic plots and summaries for each, along with
#' the specialized functions needed for typical Oceanographic
#' analysis.
#'
#' See [oce-class] for a summary of the class structure
#' and links to documentation for the many subclasses of
#' oce objects, each aligned with a class of instrument or
#' or type of dataset. For a more task-oriented approach,
#' see the several vignettes that are provided with oce,
#' and a book
#' (Kelley, Dan E. Oceanographic Analysis with R. New York: Springer-Verlag, 2018.
#' https://link.springer.com/book/10.1007/978-1-4939-8844-0) written
#' by one of the oce co-authors.
#'
#' @section Specialized Functions:
#' A key function is [read.oce()], which will attempt
#' to read Oceanographic data in raw format. This uses
#' [oceMagic()] to try to detect the file type,
#' based on the file name and contents. If this detection
#' is not possible, users will need to go beyond [read.oce()],
#' using a more specialized function, e.g. [read.ctd()] for CTD files,
#' [read.ctd.sbe()] for Teledyne-Seabird files, etc.
#'
#' @section Generic Methods:
#' A list of the generic methods in oce is provided by
#' [methods]`(class="oce")`; a few that are used frequently
#' are as follows.
#' * `[[` Finds the value of an item in the object's
#' `metadata` or `data` slot. If the item does
#' not exist, but can be calculated from the other items,
#' then the calculated value is returned. As an example of the
#' latter, consider the built-in `ctd` dataset, which does
#' not contain potential temperature, "`theta`". Using
#' `ctd[["theta"]]` therefore causes [swTheta()]
#' to be called, to calculate `theta`.
#' See \link{[[,oce-method} or type `?"[[,oce-method"`
#' to learn more about general functioning, or a specialized
#' method like \link{[[,ctd-method} for CTD data, etc.
#' * `[[<-` Alters the named item in the object's `metadata` or
#' `data` slot. If the item does not exist, it is created.
#' See \link{[[<-,oce-method} or type `?"[[<-,oce-method"`
#' to learn more about the general methodology, or a specialized
#' method like \link{[[<-,ctd-method} for CTD data, etc.
#' * `summary()` Displays some information about the object named as an
#' argument, including a few elements from its `metadata` slot
#' and some statistics of the contents of its `data` slot.
#' See \link{summary,oce-method} or type `?"summary,oce-method"`
#' to learn more about general functioning, or a specialized
#' method like \link{summary,ctd-method} for CTD data, etc.
#' * `subset()` Takes a subset of an oce object.
#' See \link{subset,oce-method} or type `?"subset,oce-method"`
#' to learn more about general functioning, or a specialized
#' method like \link{subset,ctd-method} for CTD data, etc.
#'
#' @name oce
#' @docType package
#' @keywords internal
"_PACKAGE"
NULL
#################################################################
# DEVELOPER NOTE: Update this whenever function status changes. #
#################################################################
#' Deprecated and Defunct Elements of the oce Package
#'
#' Certain functions and function arguments are still provided for
#' compatibility with older versions of \CRANpkg{oce}, but will be removed soon.
#' The \CRANpkg{oce} scheme for removing functions is similar to that used by
#' `Bioconductor`: items are marked as "deprecated" in one release, marked as
#' "defunct" in the next, and removed in the next after that. This goal is to provide a
#' gentle migration path for users who keep their packages reasonably
#' up-to-date.
#'
#' The following are marked "deprecated" in the present CRAN release of \CRANpkg{oce}.
#' Please use the replacement functions as listed below. The upcoming CRAN
#' release of \CRANpkg{oce} will mark these as "defunct", which is the
#' last step before outright removal.
#'
#' \tabular{lllll}{
#' **Deprecated** \tab **Replacement** \tab **Deprecated** \tab **Defunct** \tab **Removed** \cr
# `byteToBinary(x,"endian")` \tab [rawToBits()] \tab 1.1-1 \tab 1.1-3 \tab 1.1-4 \cr
#' }
#'
#' The following are marked "defunct", so calling them in the
#' the present version produces an error message that hints at a replacement
#' function. Once a function is marked "defunct" on one CRAN release, it will
#' be slated for outright deletion in some subsequent release.
#'
#' \tabular{lll}{
#' **Defunct** \tab **Replacement** \tab **Version**\cr
#' }
#'
#' The following functions were removed after having been marked as "deprecated"
#' in at least one CRAN release, and possibly as "defunct" in at least
#' one CRAN release. (The version number in the table is the first
#' version to lack the named function.)
#'
#' \tabular{lll}{
#' **Function** \tab **Replacement** \tab **Version**\cr
#' `addColumn()` \tab [oceSetData()] \tab 1.1-2 \cr
#' `ctdAddColumn()` \tab [oceSetData()] \tab 1.1-2 \cr
#' `ctdUpdateHeader()` \tab [oceSetMetadata()] \tab 1.1-2 \cr
#' `findInOrdered()` \tab [findInterval()] \tab 1.1-2 \cr
#' `makeSection()` \tab [as.section()] \tab 0.9.24 \cr
#' `mapMeridians()` \tab [mapGrid()] \tab 1.1-2 \cr
#' `mapZones()` \tab [mapGrid()] \tab 1.1-2 \cr
#' `oce.as.POSIXlt()` \tab [lubridate::parse_date_time()] \tab 1.1-2 \cr
#' `renameData()` \tab [oceRenameData()] \tab 1.7-9 \cr
#' `trimString()` \tab [trimws()] \tab 1.8-2 \cr
#' }
#'
#' Several \dQuote{oce} function arguments are considered "defunct", which
#' means they will be removed in the next CRAN release. They are as follows.
#'
#' * The `fill` argument of [mapPlot()] was confusing
#' to users, so it was designated as deprecated in June 2016.
#' (The confusion stemmed from subtle differences between
#' [plot()] and [polygon()], and the problem is that
#' [mapPlot()] can use either of these functions, according
#' to whether coastlines are to be filled.)
#' The functionality is preserved, in the `col` argument.
#'
#' @aliases oce-defunct
#'
#' @name oce-deprecated
#'
#' @seealso The \dQuote{Bioconductor} scheme for removing functions is
#' described at
#' `https://www.bioconductor.org/developers/how-to/deprecation/` and it is
#' extended here to function arguments.
NULL
#' Coerce Something Into an oce Object
#'
#' @details
#' This function is limited and not intended for common use.
#' In most circumstances, users should employ a function such
#' as [as.ctd()] to construct specialized oce sub-classes.
#'
#' `as.oce` creates an oce object from data contained within its
#' first argument, which may be a list, a data frame, or an object
#' of [oce-class]. (In the last case, `x` is
#' simply returned, without modification.)
#'
#' If `x` is a list containing items named `longitude` and
#' `latitude`, then [as.coastline()] is called (with
#' the specified \dots value) to create a coastline object.
#'
#' If `x` is a list created by `read_odf()` from the (as
#' yet unreleased) ODF package developed by the Bedford Institute of
#' Oceanography, then [ODF2oce()] is called (with
#' no arguments other than the first) to calculate a return value.
#' If the sub-class inference made by [ODF2oce()] is
#' incorrect, users should call that function directly, specifying
#' a value for its `coerce` argument.
#'
#' If `x` has not been created by `read_odf()`, then the names
#' of the items it contains are examined, and used to try to infer
#' the proper return value. There
#' are only a few cases (although more may be added if there is
#' sufficient user demand). The cases are as follows.
#'
#' * If `x` contains items named `temperature`,
#' `pressure` and either `salinity` or `conductivity`,
#' then an object of type [ctd-class] will be returned.
#'
#' * If `x` contains columns named `longitude` and `latitude`,
#' but no other columns, then an object of class [coastline-class]
#' is returned.
#'
#' @param x an item containing data. This may be data frame, list, or an oce object.
#'
#' @param \dots optional extra arguments, passed to conversion functions
#' [as.coastline()] or [ODF2oce()], if these are used.
#'
#' @return An [oce-class] object.
as.oce <- function(x, ...) {
if (inherits(x, "oce")) {
names <- names(x)
return(if ("EVENT_HEADER" %in% names) ODF2oce(x) else x)
}
if (!is.list(x) && !is.data.frame(x)) {
stop("x must be a list, a data frame, or an oce object")
}
names <- names(x)
if ("temperature" %in% names && "pressure" %in% names && "salinity" %in% names) {
# Assume it's a CTD
res <- as.ctd(salinity = x$salinity, temperature = x$temperature, pressure = x$pressure)
# Add other columns
for (name in names) {
if (name != "temperature" && name != "pressure" && name != "salinity") {
res <- oceSetData(res, name = name, value = x[[name]])
}
}
return(res)
}
if ("longitude" %in% names && "latitude" %in% names && length(names) == 2) {
# Assume it's a coastline
return(as.coastline(longitude = x$longitude, latitude = x$latitude))
}
# Not sure what it is, so make a generic oce object. We
# have no way to guess the unit.
res <- new("oce")
for (name in names) {
res <- oceSetData(res, name = name, value = x[[name]])
}
res
}
#' Replace the Heading for One Instrument With That of Another
#'
#' Replace the heading angles in one oce object with that from another,
#' possibly with a constant adjustment.
#'
#' @param b object holding data from an instrument whose heading is bad, but
#' whose other data are good.
#'
#' @param g object holding data from an instrument whose heading is good, and
#' should be interpolated to the time base of `b`.
#'
#' @param add an angle, in degrees, to be added to the heading.
#'
#' @return A copy of `b`, but with `b$data$heading` replaced with
#' heading angles that result from linear interpolation of the headings in
#' `g`, and then adding the angle `add`.
#'
#' @author Dan Kelley
useHeading <- function(b, g, add = 0) {
if (!"heading" %in% names(b@data)) {
stop("'b' does not have any heading data (in b@data$heading)")
}
if (!"time" %in% names(b@data)) {
stop("'b' does not have any time data (in b@data$time)")
}
if (!"heading" %in% names(g@data)) {
stop("'g' does not have any heading data (in g@data$heading)")
}
if (!"time" %in% names(g@data)) {
stop("'g' does not have any time data (in g@data$time)")
}
res <- b
t0 <- as.numeric(g@data$time[1])
if (is.na(t0)) {
stop("need first element of g@data$time to be non-NA")
}
b.t <- as.numeric(b@data$time) - t0 # FIXME: what if heading in tsSlow?
g.t <- as.numeric(g@data$time) - t0 # FIXME: what if heading in tsSlow?
res@data$heading <- approx(x = g.t, y = g@data$heading, xout = b.t)$y + add
res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep = "", collapse = ""))
res
}
#' Window an oce Object by Time or Distance
#'
#' Windows `x` on either time or distance, depending on the value of
#' `which`. In each case, values of `start` and `end` may be
#' integers, to indicate a portion of the time or distance range. If
#' `which` is `"time"`, then the `start` and `end` values
#' may also be provided as POSIX times, or character strings indicating times
#' (in time zone given by the value of `getOption("oceTz")`).
#' Note that [subset()] may be more useful than this function.
#'
#' @param x an [oce-class] object.
#'
#' @param start the start time (or distance) of the time (or space) region of
#' interest. This may be a single value or a vector.
#'
#' @param end the end time (or distance) of the time (or space) region of
#' interest. This may be a single value or a vector.
#'
#' @param frequency not permitted yet.
#'
#' @param deltat not permitted yet
#'
#' @param extend not permitted yet
#'
#' @param which string containing the name of the quantity on which sampling is
#' done. Possibilities are `"time"`, which applies the windowing on the
#' `time` entry of the `data` slot, and `"distance"`, for the
#' `distance` entry (for those objects, such as `adp`, that have this
#' entry).
#'
#' @param indexReturn boolean flag indicating whether to return a list of the
#' "kept" indices for the `time` entry of the `data` slot, as well as
#' the `timeSlow` entry, if there is one.. Either of these lists will be
#' `NULL`, if the object lacks the relevant items.
#'
#' @param debug a flag that turns on debugging.
#'
#' @param \dots ignored
#'
#' @return Normally, this is new `oce` object. However, if
#' `indexReturn=TRUE`, the return value is two-element list containing
#' items named `index` and `indexSlow`, which are the indices for the
#' `time` entry of the `data` slot (and the `timeSlow`, if it
#' exists).
#'
#' @author Dan Kelley
#'
#' @seealso [subset()] provides more flexible selection of subsets.
#'
#' @examples
#' library(oce)
#' data(adp)
#' plot(adp)
#' early <- window(adp, start = "2008-06-26 00:00:00", end = "2008-06-26 12:00:00")
#' plot(early)
#' bottom <- window(adp, start = 0, end = 20, which = "distance")
#' plot(bottom)
window.oce <- function(
x, start = NULL, end = NULL, frequency = NULL, deltat = NULL, extend = FALSE,
which = c("time", "distance"), indexReturn = FALSE,
debug = getOption("oceDebug"), ...) {
oceDebug(debug, "window.oce(...,start=",
paste(format(start), collapse = ","), ",end=",
paste(format(end), collapse = ","),
",indexReturn=", indexReturn, ",...) {\n",
unindent = 1
)
if (extend) {
stop("cannot handle extend=TRUE yet")
}
if (!is.null(frequency)) {
stop("cannot handle frequency yet")
}
if (!is.null(deltat)) {
stop("cannot handle deltat yet")
}
if (is.null(start)) {
stop("must provide start")
}
if (is.null(end)) {
stop("must provide end")
}
oceDebug(debug, "class of (x) is: ", paste(class(x), collapse = ","), "\n")
res <- x
which <- match.arg(which)
nstart <- length(start)
if (nstart != length(end)) {
stop("lengths of 'start' and 'end' must match")
}
if (which == "time") {
oceDebug(debug, "windowing by time\n")
if (!("time" %in% names(x@data))) {
warning("oce object has no @data$time vector, so window is returning it unaltered")
return(x)
}
if (is.character(start)) {
start <- as.POSIXct(start, tz = getOption("oceTz"))
}
if (is.character(end)) {
end <- as.POSIXct(end, tz = getOption("oceTz"))
}
oceDebug(debug, "tz of start:", attr(start, "tzone"), "\n")
oceDebug(debug, "tz of end:", attr(end, "tzone"), "\n")
oceDebug(debug, "tz of data$time:", attr(res@data$time, "tzone"), "\n")
nstart <- length(start)
ntime <- length(x@data$time)
keep <- rep(FALSE, ntime)
haveSlow <- "timeSlow" %in% names(x@data)
keepSlow <- if (haveSlow) rep(FALSE, length(x@data$timeSlow)) else NULL
for (w in 1:nstart) {
keep <- keep | (start[w] <= res@data$time & res@data$time <= end[w])
if (haveSlow) {
keepSlow <- keepSlow | (start[w] <= res@data$timeSlow & res@data$timeSlow <= end[w])
}
oceDebug(debug, "data window (start=", format(start[w]), ", end=", format(end[w]), ") retains", sum(keep) / ntime * 100, "percent\n")
}
if (indexReturn) {
res <- list(index = keep, indexSlow = keepSlow)
return(res)
} else {
for (name in names(res@data)) {
if ("distance" == name) {
next
}
if (length(grep("^time", name)) || is.vector(res@data[[name]])) {
if (1 == length(agrep("Slow$", name))) {
oceDebug(debug, "windowing 'slow' variable data$", name, "\n", sep = "")
res@data[[name]] <- x@data[[name]][keepSlow]
} else {
oceDebug(debug, "windowing data@", name, "\n", sep = "")
res@data[[name]] <- x@data[[name]][keep]
}
} else if (is.matrix(res@data[[name]])) {
oceDebug(debug, "windowing data@", name, ", which is a matrix\n", sep = "")
res@data[[name]] <- x@data[[name]][keep, ]
} else if (is.array(res@data[[name]])) {
oceDebug(debug, "windowing data@", name, ", which is an array\n", sep = "")
res@data[[name]] <- x@data[[name]][keep, , , drop = FALSE]
}
}
}
} else if (which == "distance") {
if (!inherits(x, "adp")) {
warning("window(..., which=\"distance\") only works for objects of class adp")
return(x)
}
if (!("distance" %in% names(x@data))) {
warning("oce object has no @data$s$distance vector, so window is returning it unaltered")
return(x)
}
oceDebug(debug, "windowing an ADP object by distance\n")
# FIXME: make it work on sections, on CTD, etc.
keep <- start <= res@data$distance & res@data$distance < end
res@data$distance <- x@data$distanc[keep]
for (name in names(res@data)) {
if (is.array(res@data[[name]]) && 3 == length(dim(x@data[[name]]))) {
res@data[[name]] <- res@data[[name]][, keep, ]
}
}
} else {
stop("unknown value of which \"", which, "\"") # cannot get here
}
# sync up some metadata that might have been altered
if (inherits(x, "adp")) {
res@metadata$numberOfSamples <- dim(res@data$v)[1]
res@metadata$numberOfCells <- dim(res@data$v)[2]
} else if (inherits(x, "adv")) {
res@metadata$numberOfSamples <- dim(res@data$v)[1]
}
oceDebug(debug, "} # window.oce()\n", unindent = 1)
res
}
#' Extract The Start of an Oce Object
#'
#' @templateVar headOrTail head
#' @template head_or_tail
#'
#' @param x an [oce-class] object.
#'
#' @param n Number of elements to extract, as for [head()].
#'
#' @param ... ignored
#'
#' @seealso [tail.oce()], which yields the end of an `oce` object.
#'
#' @author Dan Kelley
head.oce <- function(x, n = 6L, ...) {
headOrTail(x = x, n = n, headTail = head, ...)
}
headOrTail <- function(x, n = 6L, headTail = head, ...) {
# Using headTail (a) lets us handle both head.oce and
# tail.oce and also handles both positive and negative n.
res <- x
if (inherits(x, "adp") || inherits(x, "adv")) {
look <- headTail(seq_len(dim(x@data$v)[1]), n)
for (name in names(x@data)) {
if ("distance" == name) {
next
} else {
if (name == "time") {
res@data[[name]] <- headTail(x@data[[name]], n)
} else if (is.vector(x@data[[name]])) {
res@data[[name]] <- x@data[[name]][look]
} else if (is.matrix(x@data[[name]])) {
res@data[[name]] <- x@data[[name]][look, ]
} else if (is.array(x@data[[name]])) {
res@data[[name]] <- x@data[[name]][look, , ]
} else {
res@data[[name]] <- x@data[[name]][look] # for reasons unknown, 'time' is not a vector
}
}
}
} else if (inherits(x, "amsr")) {
looklon <- headTail(seq_along(x@metadata$longitude), n)
looklat <- headTail(seq_along(x@metadata$latitude), n)
res@metadata$longitude <- x@metadata$longitude[looklon]
res@metadata$latitude <- x@metadata$latitude[looklat]
for (n in names(res@data)) {
if (is.matrix(res@data[[n]])) {
res@data[[n]] <- x@data[[n]][looklon, looklat]
}
}
} else if (inherits(x, "argo")) {
for (name in names(x@metadata)) {
if (name %in% c("direction", "juldQc", "positionQc")) {
# select characters in a string
look <- headTail(seq_len(nchar(x@metadata[[name]])), n)
res@metadata[[name]] <- substr(x@metadata[[name]], look[1], tail(look, 1))
} else if (name == "flags") {
look <- headTail(seq_len(dim(x@metadata$flags[[1]])[2]), n)
for (fname in names(x@metadata$flags)) {
res@metadata$flags[[fname]] <- x@metadata$flags[[fname]][, look]
}
} else if (is.vector(x@metadata[[name]])) {
res@metadata[[name]] <- headTail(x@metadata[[name]], n)
} else if (is.matrix(x@metadata[[name]])) {
look <- headTail(seq_len(dim(x@metadata[[name]])[2]), n)
res@metadata[[name]] <- x@metadata[[name]][, look]
}
}
for (name in names(x@data)) {
if (is.vector(x@data[[name]])) {
res@data[[name]] <- headTail(x@data[[name]], n)
} else if (is.matrix(x@data[[name]])) {
look <- headTail(seq_len(dim(x@data[[name]])[2]), n)
res@data[[name]] <- x@data[[name]][, look]
} else if (name == "time") {
# for reasons unknown, time is not a vector
res@data[[name]] <- headTail(x@data[[name]], n)
} else {
warning("ignoring data item: '", name, "'")
}
}
} else if (inherits(x, "ctd")) {
for (name in names(x@data)) {
res@data[[name]] <- headTail(x@data[[name]], n)
}
for (fname in names(x@metadata$flags)) {
res@metadata$flags[[fname]] <- headTail(x@metadata$flags[[fname]], n)
}
} else if (inherits(x, "coastline")) {
for (name in c("longitude", "latitude")) {
res@data[[name]] <- headTail(x@data[[name]], n)
}
} else if (inherits(x, "echosounder")) {
look <- headTail(seq_along(x@data$latitude), n = n)
for (name in c("longitude", "latitude", "time")) {
res@data[[name]] <- x@data[[name]][look]
}
res@data$a <- x@data$a[look, ]
} else if (inherits(x, "g1sst")) {
# not in the test suite, because the files are too big, but
# looks fine manually on the following test file
# f <- dc.g1sst(2015, 10, 14, -66, -60, 41, 46)
# d <- read.g1sst(f)
looklon <- headTail(seq_along(x@metadata$longitude), n)
looklat <- headTail(seq_along(x@metadata$latitude), n)
res@metadata$longitude <- x@metadata$longitude[looklon]
res@metadata$latitude <- x@metadata$latitude[looklat]
res@data$SST <- x@data$SST[looklon, looklat]
} else if (inherits(x, "lisst")) {
look <- headTail(seq_along(x@data[[1]]), n = n)
for (name in names(x@data)) {
res@data[[name]] <- x@data[[name]][look]
}
} else if (inherits(x, "met")) {
look <- headTail(seq_along(x@data$temperature), n = n)
for (name in names(x@data)) {
res@data[[name]] <- x@data[[name]][look]
}
} else if (inherits(x, "section")) {
look <- headTail(seq_along(x@metadata$latitude), n = n)
for (name in c("stationId", "longitude", "latitude", "time")) {
res@metadata[[name]] <- x@metadata[[name]][look]
}
res@data$station <- x@data$station[look]
} else if (inherits(x, "topo")) {
looklon <- headTail(seq_along(x@data$longitude), n)
looklat <- headTail(seq_along(x@data$latitude), n)
res@data$longitude <- x@data$longitude[looklon]
res@data$latitude <- x@data$latitude[looklat]
res@data$z <- x@data$z[looklon, looklat]
} else if (inherits(x, "landsat")) {
# Actually, handling this would not be too hard, but there are some
# interlocked metadata elements that require some thought.
warning("head.oce() and tail.oce() cannot handle landsat, so returning it unaltered\n")
} else {
# The following general code works on the following objects (and is tested
# for them, in tests/testthat/test_oce.R)
#
# cm
# ladp
# lobo
# rsk
# sealevel
for (name in names(x@data)) {
# For reasons I don't understand, the 'time' items are not vectors.
if ((is.vector(x@data[[name]]) && !is.list(x@data[[name]])) || (name == "time" || name == "tstamp")) {
res@data[[name]] <- headTail(x@data[[name]], n)
} else if (is.data.frame(x@data[[name]])) {
res@data[[name]] <- headTail(x@data[[name]], n)
} else if (is.null(x@data[[name]])) {
} else {
warning("ignoring '", name, "' because it is not a vector\n")
}
}
for (fname in names(x@metadata$flags)) {
if (is.list(x@metadata$flags[[fname]])) {
#> message("x@metadata$flags$", fname, " is a list")
for (fname2 in names(x@metadata$flags[[fname]])) {
res@metadata$flags[[fname]][[fname2]] <- headTail(x@metadata$flags[[fname]][[fname2]], n)
}
} else {
res@metadata$flags[[fname]] <- headTail(x@metadata$flags[[fname]], n)
}
}
}
res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep = "", collapse = ""))
res
}
#' Extract the End of an Oce Object
#'
#' @templateVar headOrTail tail
#' @template head_or_tail
#'
#' @param x an [oce-class] object.
#'
#' @param n Number of elements to extract, as for [tail()].
#'
#' @param ... ignored
#'
#' @seealso [head.oce()], which yields the start of an `oce` object.
#'
#' @author Dan Kelley
tail.oce <- function(x, n = 6L, ...) {
headOrTail(x = x, n = n, headTail = tail, ...)
}
#' Draw a Polar Plot
#'
#' Creates a crude polar plot.
#'
#' @param r radii of points to plot.
#'
#' @param theta angles of points to plot, in degrees.
#'
#' @param debug a flag that turns on debugging. Set to 1 to get a moderate
#' amount of debugging information, or to 2 to get more.
#'
#' @param \dots optional arguments passed to the lower-level plotting
#' functions.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' r <- rnorm(50, mean = 2, sd = 0.1)
#' theta <- runif(50, 0, 360)
#' plotPolar(r, theta)
plotPolar <- function(r, theta, debug = getOption("oceDebug"), ...) {
oceDebug(debug, "plotPolar(...)\n", unindent = 1)
if (missing(r)) {
stop("must supply 'r'")
}
if (missing(theta)) {
stop("must supply 'theta'")
}
thetaRad <- theta * atan2(1, 1) / 45
x <- r * cos(thetaRad)
y <- r * sin(thetaRad)
R <- 1.2 * max(r, na.rm = TRUE)
# Rpretty <- pretty(c(0, R))
plot.new()
plot.window(c(-R, R), c(-R, R), asp = 1)
points(x, y, ...)
xa <- axis(1, pos = 0)
abline(v = 0)
th <- seq(0, 2 * atan2(1, 1) * 4, length.out = 100)
for (radius in xa[xa > 0]) {
lines(radius * cos(th), radius * sin(th))
}
abline(h = 0)
abline(v = 0)
oceDebug(debug, "} # plotPolar()\n", unindent = 1)
}
#' Interpolate 1D Data with UNESCO or Reiniger-Ross Algorithm
#'
#' Interpolate one-dimensional data using schemes that permit curvature but
#' tends minimize extrema that are not well-indicated by the data.
#'
#' Setting `method="rr"` yields the weighted-parabola algorithm of
#' Reiniger and Ross (1968). For procedure is as follows. First, the
#' interpolant for any `xout` value that is outside the range of `x`
#' is set to NA. Next, linear interpolation is used for any `xout` value
#' that has only one smaller neighboring `x` value, or one larger
#' neighboring value. For all other values of `xout`, the 4 neighboring
#' points `x` are sought, two smaller and two larger. Then two parabolas
#' are determined, one from the two smaller points plus the nearest larger
#' point, and the other from the nearest smaller point and the two larger
#' points. A weighted sum of these two parabolas provides the interpolated
#' value. Note that, in the notation of Reiniger and Ross (1968), this
#' algorithm uses `m`=2 and `n`=1. (A future version of this routine
#' might provide the ability to modify these values.)
#'
#' Setting `method="unesco"` yields the method that is used by the U.S.
#' National Oceanographic Data Center. It is described in pages 48-50 of
#' reference 2; reference 3 presumably contains the same information but it is
#' not as easily accessible. The method works as follows.
#'
#' * If there are data above 5m depth, then the surface value is taken to
#' equal to the shallowest recorded value.
#'
#' * Distance bounds are put on the four neighboring points, and the
#' Reiniger-Ross method is used for interpolated points with sufficiently four
#' close neighbors. The bounds are described in table 15 of reference 2 only
#' for so-called standard depths; in the present instance they are transformed
#' to the following rules. Inner neighbors must be within 5m for data above
#' 10m, 50m above 250m 100m above 900m, 200m above 2000m, or within 1000m
#' otherwise. Outer neighbors must be within 200m above 500m, 400m above
#' 1300m, or 1000m otherwise. If two or more points meet these criteria,
#' Lagrangian interpolation is used. If not, `NA` is used as the
#' interpolant.
#'
#' After these rules are applied, the interpolated value is compared with the
#' values immediately above and below it, and if it is outside the range,
#' simple linear interpolation is used.
#'
#' @param x the independent variable (z or p, usually).
#'
#' @param y the dependent variable.
#'
#' @param xout the values of the independent variable at which interpolation is
#' to be done.
#'
#' @param method method to use. See \dQuote{Details}.
#'
#' @return A vector of interpolated values, corresponding to the `xout`
#' values and equal in number.
#'
#' @author Dan Kelley
#'
#' @references
#' 1. R.F. Reiniger and C.K. Ross, 1968. A method of interpolation with
#' application to oceanographic data. *Deep Sea Research*, **15**,
#' 185-193.
#'
#' 2. Daphne R. Johnson, Tim P. Boyer, Hernan E. Garcia, Ricardo A.
#' Locarnini, Olga K. Baranova, and Melissa M. Zweng, 2011. World Ocean
#' Database 2009 Documentation. NODC Internal report 20. Ocean Climate
#' Laboratory, National Oceanographic Data Center. Silver Spring, Maryland.
#'
#' 3. UNESCO, 1991. Processing of oceanographic station data, 138 pp.,
#' Imprimerie des Presses Universitaires de France, United Nations Educational,
#' Scientific and Cultural Organization, France.
#'
#' @aliases oce.approx
#'
#' @examples
#' library(oce)
#' if (require(ocedata)) {
#' data(RRprofile)
#' zz <- seq(0, 2000, 2)
#' plot(RRprofile$temperature, RRprofile$depth, ylim = c(500, 0), xlim = c(2, 11))
#' # Contrast two methods
#' a1 <- oce.approx(RRprofile$depth, RRprofile$temperature, zz, "rr")
#' a2 <- oce.approx(RRprofile$depth, RRprofile$temperature, zz, "unesco")
#' lines(a1, zz)
#' lines(a2, zz, col = "red")
#' legend("bottomright", lwd = 1, col = 1:2, legend = c("rr", "unesco"), cex = 3 / 4)
#' }
oceApprox <- function(x, y, xout, method = c("rr", "unesco")) {
method <- match.arg(method)
if (missing(x)) {
stop("must supply x")
}
if (missing(y)) {
stop("must supply y")
}
lx <- length(x)
ly <- length(y)
if (lx != ly) {
stop("length of x (", lx, ") and y (", ly, ") must agree")
}
ok <- !is.na(x) & !is.na(y)
x <- x[ok]
y <- y[ok]
o <- order(x)
x <- x[o]
y <- y[o]
keep <- c(TRUE, 0 != diff(x))
x <- x[keep]
y <- y[keep]
if (missing(xout)) {
xout <- seq(min(x, na.rm = TRUE), max(x, na.rm = TRUE), length.out = lx)
} else {
if (any(is.na(xout))) {
stop("must not have any NA values in xout")
}
}
do_oceApprox(x = x, y = y, xout = xout, method = pmatch(method, c("unesco", "rr")))
}
oce.approx <- oceApprox
#' Draw a Stick Plot
#'
#' The arrows are drawn with directions on the graph that match the directions
#' indicated by the `u` and `v` components. The arrow size is set
#' relative to the units of the `y` axis, according to the value of
#' `yscale`, which has the unit of `v` divided by the unit of
#' `y`.
#' The interpretation of diagrams produced by `plotSticks` can be
#' difficult, owing to overlap in the arrows. For this reason, it is often
#' a good idea to smooth `u` and `v` before using this function.
#'
#' @param x x coordinates of stick origins.
#'
#' @param y y coordinates of stick origins. If not supplied, 0 will be used;
#' if length is less than that of x, the first number is repeated and the rest
#' are ignored.
#'
#' @param u x component of stick length.
#'
#' @param v y component of stick length.
#'
#' @param yscale scale from u and v to y (see \dQuote{Description}).
#'
#' @param add boolean, set `TRUE` to add to an existing plot.
#'
#' @param length value to be provided to [arrows()]; here, we set a
#' default that is smaller than normally used, because these plots tend to be
#' crowded in oceanographic applications.
#'
#' @param mgp 3-element numerical vector to use for [`par`]`("mgp")`. Note
#' that the default `mar` is computed from the `mgp` value.
#' 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 xlab,ylab labels for the plot axes. The default is not to label them.
#'
#' @param col color of sticks, in either numerical or character format. This is
#' made to have length matching that of `x` by a call to [rep()],
#' which can be handy in e.g. colorizing a velocity field by day.
#'
#' @param \dots graphical parameters passed down to [arrows()]. It
#' is common, for example, to use smaller arrow heads than [arrows()]
#' uses; see \dQuote{Examples}.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#'
#' # Flow from a point source
#' n <- 16
#' x <- rep(0, n)
#' y <- rep(0, n)
#' theta <- seq(0, 2 * pi, length.out = n)
#' u <- sin(theta)
#' v <- cos(theta)
#' plotSticks(x, y, u, v, xlim = c(-2, 2), ylim = c(-2, 2))
#' rm(n, x, y, theta, u, v)
#'
#' # Oceanographic example
#' data(met)
#' t <- met[["time"]]
#' u <- met[["u"]]
#' v <- met[["v"]]
#' p <- met[["pressure"]]
#' oce.plot.ts(t, p)
#' plotSticks(t, 99, u, v, yscale = 25, add = TRUE)
plotSticks <- function(x, y, u, v, yscale = 1, add = FALSE, length = 1 / 20,
mgp = getOption("oceMgp"),
mar = c(mgp[1] + 1, mgp[1] + 1, 1, 1 + par("cex")),
xlab = "", ylab = "", col = 1, ...) {
pin <- par("pin")
page.ratio <- pin[2] / pin[1]
if (missing(x)) {
stop("must supply x")
}
nx <- length(x)
if (missing(y)) {
y <- rep(0, nx)
}
if (length(y) < nx) {
y <- rep(y[1], nx)
}
if (missing(u)) {
stop("must supply u")
}
if (missing(v)) {
stop("must supply v")
}
n <- length(x)
if (length(y) != n) {
stop("lengths of x and y must match, but they are ", n, " and ", length(y))
}
if (length(u) != n) {
stop("lengths of x and u must match, but they are ", n, " and ", length(u))
}
if (length(v) != n) {
stop("lenghts of x and v must match, but they are ", n, " and ", length(v))
}
col <- rep(col, length.out = n)
par(mar = mar, mgp = mgp)
if (!add) {
plot(range(x), range(y), type = "n", xlab = xlab, ylab = ylab, ...)
}
usr <- par("usr")
yrxr <- (usr[4] - usr[3]) / (usr[2] - usr[1])
warn <- options("warn")$warn # FIXME: fails to quieten arrows()
options(warn = 0)
ok <- !is.na(u) & !is.na(v) & (u^2 + v^2) > 0
arrows(as.numeric(x[ok]),
y[ok],
(as.numeric(x[ok]) + u[ok] / yscale / yrxr * page.ratio),
(y[ok] + v[ok] / yscale),
length = length, col = col[ok],
...
)
options(warn = warn)
}
#' Add a Grid to an Existing Oce Plot
#'
#' @details
#' For plots not created by oce functions, or for missing `xat` and `yat`,
#' this is the same as a call to [grid()] with missing `nx` and
#' `ny`. However, if `xat` is the return value from certain oce functions,
#' a more sophisticated grid is constructed. The problem with [grid()] is
#' that it cannot handle axes with non-uniform grids, e.g. those with time axes
#' that span months of differing lengths.
#'
#' As of early February 2015, `oce.grid` handles `xat` produced as the
#' return value from the following functions: [imagep()] and
#' [oce.plot.ts()], [plot,adp-method()],
#' [plot,echosounder-method()], and [plotTS()].
#' It makes no sense to try to use `oce.grid` for multipanel oce plots,
#' e.g. the default plot from [plot,adp-method()].
#'
#' @examples
#' library(oce)
#' i <- imagep(volcano)
#' oce.grid(i, lwd = 2)
#'
#' data(sealevel)
#' i <- oce.plot.ts(sealevel[["time"]], sealevel[["elevation"]])
#' oce.grid(i, col = "red")
#'
#' data(ctd)
#' i <- plotTS(ctd)
#' oce.grid(i, col = "red")
#'
#' data(adp)
#' i <- plot(adp, which = 1)
#' oce.grid(i, col = "gray", lty = 1)
#'
#' data(echosounder)
#' i <- plot(echosounder)
#' oce.grid(i, col = "pink", lty = 1)
#'
#' @param xat either a list of x values at which to draw the grid, or the return value from an oce plotting function
#'
#' @param yat a list of y values at which to plot the grid (ignored if `gx` was a return value from an oce plotting function)
#'
#' @param col color of grid lines (see [par()])
#'
#' @param lty type for grid lines (see [par()])
#'
#' @param lwd width for grid lines (see [par()])
oce.grid <- function(xat, yat, col = "lightgray", lty = "dotted", lwd = par("lwd")) {
if (missing(xat) && missing(yat)) {
grid(col = col, lty = lty, lwd = lwd)
} else {
if (is.list(xat)) {
# following over-rides the args
yat <- xat$yat
xat <- xat$xat
}
if (!missing(xat)) {
abline(v = xat, col = col, lty = lty, lwd = lwd)
}
if (!missing(yat)) {
abline(h = yat, col = col, lty = lty, lwd = lwd)
}
}
}
#' Oce Variant of plot.ts
#'
#' Plot a time-series, obeying the timezone and possibly drawing the range in
#' the top-left margin.
#'
#' Depending on the version of R, the standard [plot()] and
#' [plot.ts()] routines will not obey the time zone of the data.
#' This routine gets around that problem. It can also plot the time range in
#' the top-left margin, if desired; this string includes the timezone, to
#' remove any possible confusion.
#' The time axis is drawn with [oce.axis.POSIXct()].
#'
#' @param x the times of observations. If this is not a [POSIXt] object, then an
#' attempt is made to convert it to one using [as.POSIXct()].
#'
#' @param y the observations.
#'
#' @param type plot type, `"l"` for lines, `"p"` for points.
#'
#' @param xlim optional limit for x axis. This has an additional effect,
#' beyond that for conventional R functions: it effectively windows the data,
#' so that autoscaling will yield limits for y that make sense within the
#' window.
#'
#' @param ylim optional limit for y axis.
#'
#' @param log a character value that must be either empty (the default) for linear
#' `y` axis, or `"y"` for logarithmic `y` axis. (Unlike
#' [plot.default()] etc., `oce.plot.ts` does not permit
#' logarithmic time, or `x` axis.)
#'
#' @param logStyle a character value that indicates how to draw the y axis, if
#' `log="y"`. If it is `"r"` (the default) then the conventional R style is used,
#' in which a logarithmic transform connects y values to position on the "page"
#' of the plot device, so that tics will be nonlinearly spaced, but not
#' organized by integral powers of 10. However, if it is `"decade"`, then
#' the style will be that used in the scientific literature, in which large
#' tick marks are used for integral powers of 10, with smaller tick marks
#' at integral multiples of those powers, and with labels that use exponential
#' format for values above 100 or below 0.01. The value of `logStyle` is passed
#' to [oceAxis()], which draws the axis.
#'
#' @param drawTimeRange an optional indication of whether/how to draw a time range,
#' in the top-left margin of the plot; see [oce.axis.POSIXct()] for details.
#'
#' @param simplify an integer value that indicates
#' whether to speed up `type="l"` plots by replacing the data
#' with minimum and maximum values within a subsampled time mesh.
#' This can speed up plots of large datasets (e.g. by factor 20 for 10^7 points),
#' sometimes with minor changes in appearance.
#' This procedure is skipped if `simplify` is `NA` or
#' if the number of visible data points is less than 5 times `simplify`.
#' Otherwise, `oce.plot.ts` creates `simplify` intervals ranging across
#' the visible time range. Intervals with under 2 finite
#' `y` data are ignored. In the rest, `y` values
#' are replaced with their range, and `x` values are replaced
#' with the repeated midpoint time. Thus, each retained sub-interval
#' has exactly 2 data points.
#' A warning is printed if this replacement is done.
#' The default value of `simplify` means that cases with
#' under 2560 visible points are plotted conventionally.
#'
#' @param fill boolean, set `TRUE` to fill the curve to zero (which it
#' does incorrectly if there are missing values in `y`).
#' @param col The colours for points (if `type=="p"`) or lines (if `type=="l"`).
#' For the `type="p"` case,
#' if there are fewer `col` values than there are `x` values, then the `col` values
#' are recycled in the standard fashion.
#' For the `type="l"` case, the line is plotted in the first colour specified.
#'
#' @param pch character code, used if `type=="p"`.
#' If there are fewer `pch` values than there are `x` values, then the `pch` values
#' are recycled in the standard fashion.
#' See [points()] for the possible values for `pch`.
#'
#' @param cex numeric character expansion factor for points on plots, ignored unless
#' `type` is `"p"`. This may be a single number, applied to all points, or
#' a vector of numbers to be applied to the points in sequence. If there are
#' fewer `pch` values than there are `x` values, then the `pch` values are recycled
#' in the standard fashion. See [par()] for more on `cex`.
#'
#' @param cex.axis,cex.lab,cex.main numeric character expansion factors for axis numbers,
#' axis names and plot titles; see [par()].
#'
#' @param flipy Logical, with `TRUE` indicating that the graph
#' should have the y axis reversed, i.e. with smaller values at
#' the bottom of the page.
#'
#' @param xlab name for x axis; defaults to `""`.
#'
#' @param ylab name for y axis; defaults to the plotted item.
#'
#' @param xaxs control x axis ending; see [`par`]`("xaxs")`.
#'
#' @param yaxs control y axis ending; see [`par`]`("yaxs")`.
#'
#' @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")` to set margins.
#' The default value uses significantly tighter margins than is the norm in R,
#' which gives more space for the data. However, in doing this, the existing
#' `par("mar")` value is ignored, which contradicts values that may have
#' been set by a previous call to [drawPalette()]. To get plot with
#' a palette, first call [drawPalette()], then call
#' `oce.plot.ts` with `mar=par("mar")`.
#'
#' @param main title of plot.
#'
#' @param despike boolean flag that can turn on despiking with
#' [despike()].
#'
#' @param axes boolean, set to `TRUE` to get axes plotted
#'
#' @param tformat optional format for labels on the time axis
#'
#' @param marginsAsImage boolean indicating whether to set the right-hand
#' margin to the width normally taken by an image drawn with
#' [imagep()].
#'
#' @param grid if `TRUE`, a grid will be drawn for each panel. (This
#' argument is needed, because calling [grid()] after doing a
#' sequence of plots will not result in useful results for the individual
#' panels.
#'
#' @param grid.col color of grid
#'
#' @param grid.lty line type of grid
#'
#' @param grid.lwd line width of grid
#'
#' @param debug a flag that turns on debugging. Set to 1 to get a moderate
#' amount of debugging information, or to 2 to get more.
#'
#' @param \dots graphical parameters passed down to [plot()].
#'
#' @return A list is silently returned, containing `xat` and `yat`,
#' values that can be used by [oce.grid()] to add a grid to the plot.
#'
#' @author Dan Kelley and Clark Richards
#'
#' @examples
#' library(oce)
#' t0 <- as.POSIXct("2008-01-01", tz = "UTC")
#' t <- seq(t0, length.out = 48, by = "30 min")
#' y <- sin(as.numeric(t - t0) * 2 * pi / (12 * 3600))
#' oce.plot.ts(t, y, type = "l", xaxs = "i")
#' # Show how col, pch and cex get recycled
#' oce.plot.ts(t, y,
#' type = "p", xaxs = "i",
#' col = 1:3, pch = c(rep(1, 6), rep(20, 6)), cex = sqrt(1:6)
#' )
#' # Trimming x; note the narrowing of the y view
#' oce.plot.ts(t, y, type = "p", xlim = c(t[6], t[12]))
#' # Flip the y axis
#' oce.plot.ts(t, y, flipy = TRUE)
oce.plot.ts <- function(
x, y, type = "l", xlim, ylim, log = "", logStyle = "r", flipy = FALSE, xlab, ylab,
drawTimeRange, simplify = 2560, fill = FALSE, col = par("col"), pch = par("pch"),
cex = par("cex"), cex.axis = par("cex.axis"), cex.lab = par("cex.lab"), cex.main = par("cex.main"),
xaxs = par("xaxs"), yaxs = par("yaxs"),
mgp = getOption("oceMgp"),
mar = c(mgp[1] + if (nchar(xlab) > 0) 1.5 else 1, mgp[1] + 1.5, mgp[2] + 1, mgp[2] + 3 / 4),
main = "",
despike = FALSE,
axes = TRUE, tformat,
marginsAsImage = FALSE,
grid = FALSE, grid.col = "lightgray", grid.lty = "dotted", grid.lwd = par("lwd"),
debug = getOption("oceDebug"),
...) {
if (is.function(x)) {
stop("x cannot be a function")
}
if (!inherits(x, "POSIXt")) {
x <- as.POSIXct(x)
}
if (missing(xlab)) {
xlab <- ""
}
if (missing(ylab)) {
ylab <- deparse(substitute(expr = y, env = environment()))
}
if (missing(drawTimeRange)) {
drawTimeRange <- getOption("oceDrawTimeRange", TRUE)
}
# ocex <- par("cex")
# par(cex=cex)
debug <- min(debug, 4)
oceDebug(debug, "oce.plot.ts(...,debug=", debug, ",",
argShow(type),
argShow(flipy),
argShow(log), "\n",
sep = "", unindent = 1, style = "bold"
)
oceDebug(debug,
" ",
argShow(mar),
argShow(mgp), "\n",
sep = "", style = "bold"
)
oceDebug(debug,
" ",
argShow(cex),
argShow(cex.axis),
argShow(cex.lab),
"...) {\n",
sep = "", style = "bold"
)
if (!is.logical(flipy)) {
stop("flipy must be TRUE or FALSE")
}
if (!log %in% c("", "y")) {
stop("log must be either an empty string or \"y\", not \"", log, "\" as given")
}
if (!is.na(simplify) && (!is.numeric(simplify) || simplify <= 0)) {
stop("simplify must be NA or a positive number, but it is ", simplify)
}
# oceDebug(debug, "length(x)", length(x), "; length(y)", length(y), "\n")
# oceDebug(debug, "marginsAsImage=", marginsAsImage, ")\n")
oceDebug(debug, "x has timezone", attr(x[1], "tzone"), "\n")
# Repeat col, pch and cex to the right length, for possible trimming later.
drawingPoints <- type == "p" || type == "o" || type == "b"
if (drawingPoints) {
# BOOKMARK 1a: alter 1b if more added or subtracted here
nx <- length(x)
col <- rep(col, length.out = nx)
pch <- rep(pch, length.out = nx)
cex <- rep(cex, length.out = nx)
oceDebug(debug, "made col, pch and cex of length ", nx, " to match length(x)\n")
}
pc <- paletteCalculations(maidiff = rep(0, 4))
oceDebug(debug, as.character(dput(pc)), "\n", style = "red")
par(mgp = mgp, mar = mar)
args <- list(...)
xlimGiven <- !missing(xlim)
# Trim data to visible window (to improve speed and to facilitate 'simplify' calculation)
if (xlimGiven) {
if (2 != length(xlim)) {
stop("'xlim' must be of length 2, but it is of length ", length(xlim))
}
if (xlim[2] <= xlim[1]) {
stop(
"the elements of xlim must be in order, but they are ",
format(xlim[1]), " and ", format(xlim[2]), ", respectively"
)
}
# Comment-out next line for issue 1508, since trim_ts
# fails if times are NA.
# ends <- trim_ts(x, xlim, 0.04)
dx <- diff(as.numeric(xlim))
keep <- (xlim[1] - dx * 0.04) < x & x < (xlim[2] + dx * 0.04)
x <- x[keep]
y <- y[keep]
if (drawingPoints) {
# BOOKMARK 1b: alter 1a if more added or subtracted here
col <- col[keep]
pch <- pch[keep]
cex <- cex[keep]
}
}
if (length(y) == 1) {
y <- rep(y, length(x))
}
if (despike) {
y <- despike(y)
}
if (marginsAsImage) {
# FIXME: obey their mar?
the.mai <- pc$mai0
the.mai <- clipmin(the.mai, 0) # just in case
oceDebug(debug, "the.mai=", vectorShow(the.mai))
par(mai = the.mai)
drawPalette(mai = rep(0, 4))
}
# Find data ranges. Note that x is a time, so it's either going to be NA or
# sensible; thus the na.rm argument to range() is suitable for trimming bad
# values. However, for y, we emulate plot(), by trimming (and warning).
if ("y" %in% log) {
yBAD <- (!is.finite(y)) | y <= 0.0
nyBAD <- sum(yBAD)
if (nyBAD > 0L) {
warning(nyBAD, " y value <= 0 omitted from logarithmic oce.plot.ts\n")
x <- x[!yBAD]
y <- y[!yBAD]
}
}
# Handle 'simplify' argument
nx <- length(x)
if (type == "l" && is.numeric(simplify) && nx > (5L * simplify)) {
warning("simplifying a large dataset; set simplify=NA to see raw data\n")
xgrid <- seq(min(x, na.rm = TRUE), max(x, na.rm = TRUE), length.out = simplify)
df <- data.frame(x, y)
# Tests N=1e8 suggest split(X,findInterval()) is 2X faster than split(X,cut())
#>>> dfSplit <- split(df, cut(df$x, xgrid))
dfSplit <- split(df, as.factor(findInterval(df$x, xgrid)))
# Compute within the sub-intervals, inserting NAs when no data there
tz <- attr(x, "tzone") # cause gridded x to inherit timezone from original x
x <- rep(unname(sapply(dfSplit, function(DF) if (length(DF$x) > 2) mean(DF$x, na.rm = TRUE) else NA)), each = 2)
x <- numberAsPOSIXct(x, tz = tz)
ymin <- unname(sapply(dfSplit, function(DF) if (length(DF$y) > 2) min(DF$y, na.rm = TRUE) else NA))
ymax <- unname(sapply(dfSplit, function(DF) if (length(DF$y) > 2) max(DF$y, na.rm = TRUE) else NA))
y <- as.vector(rbind(ymin, ymax))
# Remove any segments for which min and max could not be computed
bad <- !is.finite(y)
x <- x[!bad]
y <- y[!bad]
}
xrange <- range(x, na.rm = TRUE)
yrange <- range(y, finite = TRUE)
maybeflip <- function(y) if (flipy) rev(sort(y)) else y
if (!is.finite(yrange[1])) {
plot(xrange, c(0, 1),
axes = FALSE, xaxs = xaxs, yaxs = yaxs,
xlim = if (xlimGiven) xlim else xrange,
ylim = if (missing(ylim)) maybeflip(range(y, na.rm = TRUE)) else maybeflip(ylim),
xlab = xlab, ylab = ylab, type = "n", log = log, col = col, pch = pch, cex = cex
)
oce.axis.POSIXct(1, drawTimeRange = FALSE)
box()
mtext("bad data", side = 3, line = -1, cex = cex)
warning("no valid data for '", ylab, "'", sep = "")
oceDebug(debug, "} # oce.plot.ts()\n", unindent = 1, style = "bold")
return()
} else {
if (fill) {
xx <- c(x[1], x, x[length(x)])
yy <- c(0, y, 0)
plot(x, y,
axes = FALSE, xaxs = xaxs, yaxs = yaxs,
xlim = if (xlimGiven) xlim else range(x, na.rm = TRUE),
ylim = if (missing(ylim)) maybeflip(range(y, na.rm = TRUE)) else maybeflip(ylim),
xlab = "", ylab = "",
type = type, col = col, cex = cex, pch = pch, log = log, ...
)
mtext(xlab, side = 1, cex = cex.lab * par("cex"), line = mgp[1])
mtext(ylab, side = 2, cex = cex.lab * par("cex"), line = mgp[1])
fillcol <- if ("col" %in% names(args)) args$col else "lightgray" # FIXME: should this be a formal argument?
do.call(polygon, list(x = xx, y = yy, col = fillcol))
} else {
plot(x, y,
axes = FALSE, xaxs = xaxs, yaxs = yaxs,
xlim = if (missing(xlim)) NULL else xlim,
ylim = if (missing(ylim)) maybeflip(range(y, na.rm = TRUE)) else maybeflip(ylim),
xlab = "", ylab = "",
type = type, col = col, cex = cex, cex.axis = cex.axis, cex.lab = cex.lab, pch = pch, log = log, ...
)
# mtext(paste("TEST: xlab at mgp[1]", xlab), side=1, cex=cex.lab*par("cex"), line=mgp[1])
mtext(ylab, side = 2, cex = cex.lab * par("cex"), line = mgp[1])
}
xat <- NULL
yat <- NULL
if (axes) {
xaxt <- list(...)["xaxt"]
drawxaxis <- !is.null(xaxt) && xaxt != "n"
yaxt <- list(...)["yaxt"]
drawyaxis <- !is.null(yaxt) && yaxt != "n"
if (drawxaxis) {
xlabs <- oce.axis.POSIXct(1,
x = x, drawTimeRange = drawTimeRange, main = main,
mgp = mgp,
xlim = if (missing(xlim)) range(x) else xlim,
cex.axis = cex.axis, cex.main = cex.main, cex.lab = cex.lab,
tformat = tformat,
debug = debug - 1
)
xat <- xlabs
oceDebug(debug, "drawing x axis; set xat=c(", paste(xat, collapse = ","), ")", "\n", sep = "")
}
if (FALSE && grid) { # FIXME: don't w do this below? Why here also?
lwd <- par("lwd")
if (drawxaxis) {
abline(v = xlabs, col = "lightgray", lty = "dotted", lwd = lwd)
}
yaxp <- par("yaxp")
if (log == "y") {
abline(h = axTicks(2), col = "lightgray", lty = "dotted", lwd = lwd)
} else {
abline(
h = seq(yaxp[1], yaxp[2], length.out = 1 + yaxp[3]),
col = "lightgray", lty = "dotted", lwd = lwd
)
}
}
box()
if (drawyaxis) {
if (log != "y" || logStyle == "r") {
axis(2, cex.axis = cex.axis, cex = cex.axis)
yat <- axis(4, labels = FALSE)
} else if (log == "y" && logStyle == "decade") {
yat <- oceAxis(2, logStyle = logStyle, cex.axis = cex.axis, cex = cex.axis)
oceAxis(4, logStyle = logStyle, labels = FALSE)
} else if (log == "y") {
stop("if log=\"y\", then logStyle must be \"r\" or \"decade\", not \"", logStyle, "\"")
}
}
}
if (grid) {
if (log == "y") {
at <- axTicks(2)
lat <- log10(at)
powerOfTen <- abs(lat - round(lat)) < 1e-4
abline(
h = if (any(powerOfTen)) at[powerOfTen] else at,
col = grid.col, lty = grid.lty, lwd = grid.lwd
)
} else {
abline(h = axTicks(2), col = grid.col, lty = grid.lty, lwd = grid.lwd)
}
abline(v = axTicks(1), col = grid.col, lty = grid.lty, lwd = grid.lwd)
}
oceDebug(debug, "} # oce.plot.ts()\n", unindent = 1, style = "bold")
invisible(list(xat = xat, yat = yat))
}
} # oce.plot.ts()
#' Trim an oce File
#'
#' Create an oce file by copying the first `n` data chunks of another such file.
#' This can be useful in supplying small sample files for bug reports. Only
#' a few file types (as inferred with [oceMagic()]) are permitted.
#'
#' @param infile name of an AD2CP source file.
#'
#' @param n integer indicating the number of data chunks to keep. The default is
#' to keep 100 chunks, a common good choice for sample files.
#'
#' @param outfile optional name of the new file to be created. If this is not
#' supplied, a default is used, by adding `_trimmed` to the base filename, e.g.
#' for an AD2CP file named `"a.ad2cp"`, the constructed value of `outfile` will
#' be `a_trimmed.ad2cp`.
#'
#' @param debug an integer value indicating the level of debugging. If
#' this is 1L, then a brief indication is given of the processing steps. If it
#' is > 1L, then information is given about each data chunk, which can yield
#' very extensive output.
#'
#' @return `oceFileTrim()` returns the name of the output file, either provided
#' in the `outfile` parameter or constructed by this function.
#'
#' @family functions that trim data files
#'
#' @section Sample of Usage:
#' \preformatted{
#' # Can only be run by the developer, since it uses a private file.
#' f <- "~/Dropbox/oce_secret_data/ad2cp/byg_trimmed.ad2cp"
#' if (file.exists(f)) {
#' oceFileTrim(f, 10L) # this file holds 100 data segments
#' }
#' }
#' @author Dan Kelley
oceFileTrim <- function(infile, n = 100L, outfile, debug = getOption("oceDebug")) {
magic <- oceMagic(infile)
allowed <- c("adp/nortek/ad2cp", "adp/rdi")
if (!magic %in% allowed) {
stop("this file, of type ", magic, ", is not understood; only files of the following types are allowed: ", paste(allowed, collapse = ", "))
}
switch(magic,
"adp/nortek/ad2cp" = adpAd2cpFileTrim(infile, n, outfile, debug),
"adp/rdi" = adpRdiFileTrim(infile, n, outfile, debug)
)
}
#' Edit an Oce Object
#'
#' Edit an element of an oce object, inserting a note in the processing
#' log of the returned object.
#'
#' There are several ways to use this function.
#'
#' 1. If both an `item` and `value` are supplied, then
#' either the object's metadata or data slot may be altered. There are
#' two ways in which this can be done.
#'
#' * Case 1A. If the `item` string does not contain an
#' `@` character, then the `metadata` slot is examined
#' for an entry named `item`, and that is modified if so.
#' Alternatively, if `item` is found in `metadata`, then
#' that value is modified. However, if `item` is not found in
#' either `metadata` or `data`, then an error is reported
#' (see 1B for how to add something that does not yet exist).
#'
#' * Case 1B. If the `item` string contains
#' the `@` character, then the text to the left of that character
#' must be either `"metadata"` or `"data"`, and it names the slot
#' in which the change is done. In contrast with case 1A, this will
#' *create* a new item, if it is not already in existence.
#'
#'
#' 2. If `item` and `value` are not supplied, then `action` must
#' be supplied. This is a character string specifying some action to be
#' performed on the object, e.g. a manipulation of a column. The action must
#' refer to the object as `x`, as in Example 2.
#'
#' In any case, a log entry is stored in the object, to document the change.
#' Indeed, this is the main benefit to using this function, instead of altering
#' the object directly. The log entry will be most useful if it contains a
#' brief note on the `reason` for the change, and the name of the
#' `person` doing the work.
#'
#' @aliases oce.edit
#'
#' @param x an [oce-class] object. The exact action of [oceEdit()] depends
#' on the sub-class of `x`.
#'
#' @param item if supplied, a character string naming an item in the object's
#' `metadata` or `data` slot, the former being checked first.
#' An exception is if `item` starts with `"data@"` or
#' `"metadata@"`, in which case the named slot is updated with a changed
#' value of the contents of `item` after the `@` character.
#'
#' @param value new value for `item`, if both supplied.
#'
#' @param action optional character string containing R code to carry out some
#' action on the object.
#'
#' @param reason character string giving the reason for the change.
#'
#' @param person character string giving the name of person making the change.
#'
#' @param debug an integer that specifies a level of debugging, with 0 or less
#' indicating no debugging, and 1 or more indicating debugging.
#'
#' @return A [oce-class] object, altered appropriately, and with a log item
#' indicating the nature of the alteration.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' data(ctd)
#' # Example 1: change latitude
#' ctd2 <- oceEdit(ctd,
#' item = "latitude", value = 47.8879,
#' reason = "illustration", person = "Dan Kelley"
#' )
#' # Example 2: add 0.1 dbar to pressure
#' ctd3 <- oceEdit(ctd, action = "x@@data$pressure<-x@@data$pressure+0.1")
oceEdit <- function(x, item, value, action, reason = "", person = "", debug = getOption("oceDebug")) {
oceDebug(debug, "oceEdit() {\n", unindent = 1)
if (!inherits(x, "oce")) {
stop("method is only for oce objects")
}
if (missing(item) && missing(value) && missing(action)) {
x@processingLog <- processingLogAppend(x@processingLog, paste(deparse(match.call()), sep = "", collapse = ""))
oceDebug(debug, "} # oceEdit()\n", unindent = 1)
return(x)
}
slot <- NULL
if (!missing(item)) {
if (missing(value)) {
stop("must supply a value")
}
# Split out the slotname, if any.
if (length(grep("@", item))) {
slot <- gsub("@.*$", "", item)
if (slot != "metadata" && slot != "data") {
stop("slot must be 'metadata' or 'data'")
}
item <- gsub("^.*@", "", item)
}
if (missing(value)) {
stop("must supply a 'value' for this 'item'")
}
if (inherits(x, "adv")) {
oceDebug(debug, "object is an ADV\n")
hpr <- 0 < length(grep("heading|pitch|roll", item))
if (hpr) {
# FIXME: I think this violates the 1A rule on creating new data,
# FIXME: but I am retaining this since it's years old.
# FIXME: why are adp and adv handled differently, anyway? Is
# FIXME: this a fast/slow variable issue?
x@data[[item]] <- value
} else {
if (!is.null(slot)) {
slot(x, slot)[[item]] <- value
} else if (item %in% names(x@metadata)) {
x@metadata[[item]] <- value
} else if (item %in% names(x@data)) {
x@data[[item]] <- value
} else {
stop("nothing named '", item, "' in object's metadata or data")
}
}
} else if (inherits(x, "adp")) {
oceDebug(debug, "object is an ADP\n")
hpr <- 0 < length(grep("heading|pitch|roll", item))
if (hpr) {
# FIXME: I think this violates the 1A rule on creating new data,
# FIXME: but I am retaining this since it's years old.
# FIXME: why are adp and adv handled differently, anyway? Is
# FIXME: this a fast/slow variable issue?
x@data[[item]] <- value
} else {
if (!is.null(slot)) {
slot(x, slot)[[item]] <- value
} else if (item %in% names(x@metadata)) {
x@metadata[[item]] <- value
} else if (item %in% names(x@data)) {
x@data[[item]] <- value
} else {
stop("nothing named '", item, "' in object's metadata or data")
}
}
} else if ("instrumentType" %in% names(x@metadata) && x@metadata$instrumentType == "aquadopp-hr") {
# FIXME: what if S4?
oceDebug(debug, "About to try editing AQUADOPP ...\n")
hpr <- 0 < length(grep("heading|pitch|roll", item)) # FIXME: possibly aquadopp should have tsSlow
x@data[[item]] <- value
if (hpr) {
x@data[[item]] <- value
} else {
if (!is.null(slot)) {
slot(x, slot)[[item]] <- value
} else if (item %in% names(x@metadata)) {
x@metadata[[item]] <- value
} else if (item %in% names(x@data)) {
x@data[[item]] <- value
} else {
stop("nothing named '", item, "' in object's metadata or data")
}
}
oceDebug(debug, "...AQUADOPP edited\n")
} else {
oceDebug(debug, "general object; item=\"", item, "\"; slot=\"", slot, "\"\n", sep = "")
if (!is.null(slot)) {
slot(x, slot)[[item]] <- value
} else if (item %in% names(x@metadata)) {
x@metadata[[item]] <- value
} else if (item %in% names(x@data)) {
x@data[[item]] <- value
} else {
stop("nothing named '", item, "' in object's metadata or data")
}
}
} else if (!missing(action)) {
t <- try(eval(parse(text = action)), silent = TRUE)
if (inherits(t, "try-error")) {
stop("cannot parse and evaluate the supplied action")
}
} else {
stop("must supply either an 'item' plus a 'value', or an 'action'")
}
x@processingLog <- processingLogAppend(x@processingLog, paste(deparse(match.call()), sep = "", collapse = ""))
oceDebug(debug, "} # oceEdit()\n", unindent = 1)
x
}
oce.edit <- oceEdit
#' Write the Data Portion of Object to a File
#'
#' The output has a line containing the names of the columns in `x$data`,
#' each enclosed in double quotes. After that line are lines for the data
#' themselves. The default is to separate data items by a single space
#' character, but this can be altered by using a `sep` argument in the
#' `...` list; see [utils::write.table()].
#'
#' This function is little more than a thin wrapper around
#' [utils::write.table()], the only difference being that row names
#' are omitted here, making for a file format that is more conventional in
#' Oceanography.
#'
#' @param x an [oce-class] object.
#'
#' @param file file name, as passed to [utils::write.table()]. Use
#' `""` to get a listing in the terminal window.
#'
#' @param ... optional arguments passed to [utils::write.table()].
#'
#' @return The value returned by [utils::write.table()].
#'
#' @author Dan Kelley
#'
#' @seealso `[utils::write.table()], which does the actual work.
oce.write.table <- function(x, file = "", ...) {
if (!inherits(x, "oce")) {
stop("method is only for oce objects")
}
if (!("row.names" %in% names(list(...)))) {
write.table(x@data, file, row.names = FALSE, ...)
} else {
write.table(x@data, file, ...)
}
}
#' Standard Oceanographic Depths
#'
#' This returns a vector of numbers that build upon the shorter lists
#' provided in Chapter 10 of reference 1 and the more modern World
#' Ocean Atlases (e.g. reference 2).
#' With the default call,
#' i.e. with `n=0`, the result is
#' \code{c(0, 10, 20, 30, 40, 50, 75, 100, 125, 150, 200, 250,
#' seq(300, 1500, by=100), 1750, seq(2000, 10000, by=500))}.
#' For higher values of `n`, progressively more and more values
#' are added between each pair in this sequence.
#' See the documentation for
#' [sectionGrid()] for how `standardDepths` can be used
#' in gridding data for section plots.
#'
#' @param n Integer specifying the number of subdivisions to insert between
#' each of the stated levels. For exmple, setting `n=1` puts a 5m level
#' between the 0 and 10m levels, and `n=2` puts 3.33 and 6.66 between
#' 0 and 10m.
#'
#' @return A vector of depths that are more closely spaced for small values,
#' i.e. a finer grid near the ocean surface.
#'
#' @examples
#' depth <- standardDepths()
#' depth1 <- standardDepths(1)
#' plot(depth, depth)
#' points(depth1, depth1, col = 2, pch = 20, cex = 1 / 2)
#'
#' @author Dan Kelley
#'
#' @references
#' 1. Sverdrup, H U, Martin W Johnson, and Richard H Fleming. The Oceans,
#' Their Physics, Chemistry, and General Biology. New York: Prentice-Hall, 1942.
# `https://ark.cdlib.org/ark:/13030/kt167nb66r`
# next worked well most of the time, but I got a failure on 2021-08-07 and
# I just don't see the point in retaining a link that will *ever* fail, given
# the possible consequence in terms of CRAN.
#' \code{https://publishing.cdlib.org/ucpressebooks/view?docId=kt167nb66r}
#'
#' 2.Locarnini, R. A., A. V. Mishonov, J. I. Antonov, T. P. Boyer,
#' H. E. Garcia, O. K. Baranova, M. M. Zweng, D. R. Johnson, and
#' S. Levitus. \dQuote{World Ocean Atlas 2009 Temperature.}
#' US Government printing Office, 2010.
standardDepths <- function(n = 0) {
d <- c(
0, 10, 20, 30, 40, 50, 75, 100, 125, 150, 200, 250,
seq(300, 1500, by = 100), 1750, seq(2000, 10000, by = 500)
)
n <- as.integer(n)
if (n < 0) {
stop("cannot have negative n")
}
if (n == 0) {
return(d)
}
ul <- cbind(head(d, -1), tail(d, -1))
res <- NULL
for (i in seq_len(nrow(ul))) {
res <- c(res, approx(c(0, 1), ul[i, ], seq(0, 1, by = 1 / (n + 1)))$y)
}
unique(res) # remove duplicates from one "l" being the next "u"
}
#' Find the Type of an Oceanographic Data File
#'
#' `oceMagic` tries to infer the file type, based on the data
#' within the file, the file name, or a combination of the two.
#'
#' `oceMagic` was previously called `oce.magic`, but that
#' alias was removed in version 0.9.24; see [oce-defunct].
#'
#' @param file a connection or a character string giving the name of the file
#' to be checked.
#'
#' @template encodingTemplate
#'
#' @param debug an integer, set non-zero to turn on debugging. Higher values
#' indicate more debugging.
#'
#' @return A character string indicating the file type, or `"unknown"`, if
#' the type cannot be determined. If the result contains `"/"` characters,
#' these separate a list describing the file type, with the first element being
#' the general type, the second element being the manufacturer, and the third
#' element being the manufacturer's name for the instrument. For example,
#' `"adp/nortek/aquadopp"` indicates a acoustic-doppler profiler made by
#' NorTek, of the model type called Aquadopp.
#'
#' @author Dan Kelley
#'
#' @seealso This is used mainly by [read.oce()].
oceMagic <- function(file, encoding = "latin1", debug = getOption("oceDebug")) {
filename <- file
oceDebug(debug, paste("oceMagic(file=\"", filename, "\") {\n", sep = ""), unindent = 1, style = "bold", sep = "")
isdir <- file.info(file)$isdir
if (is.finite(isdir) && isdir) {
tst <- file.info(paste(file, "/", file, "_MTL.txt", sep = ""))$isdir
if (!is.na(tst) && !tst) {
oceDebug(debug, "} # oceMagic returning landsat\n", unindent = 1)
return("landsat")
} else {
stop("please supply a file name, not a directory name")
}
}
if (is.character(file)) {
oceDebug(debug, "'file' is a character value\n")
if (grepl(".asc$", filename)) {
someLines <- readLines(file, n = 1L)
if (42 == length(strsplit(someLines[1], " ")[[1]])) {
oceDebug(debug, "} # oceMagic returning lisst\n", unindent = 1, style = "bold")
return("lisst")
}
}
if (grepl(".adr$", filename)) {
oceDebug(debug, "file names ends in .adr, so this is an adv/sontek/adr file.\n")
oceDebug(debug, "} # oceMagic returning adv/sontek/adr\n", unindent = 1)
return("adv/sontek/adr")
}
if (grepl(".rsk$", filename)) {
oceDebug(debug, "file names ends with \".rsk\", so this is an RBR/rsk file.\n")
oceDebug(debug, "} # oceMagic returning RBR/rsk\n", unindent = 1)
return("RBR/rsk")
}
if (grepl(".s4a.", filename)) {
oceDebug(debug, "file names contains \".s4a.\", so this is an interocean S4 file.\n")
oceDebug(debug, "} # oceMagic returning interocean/s4\n", unindent = 1)
return("interocean/s4")
}
if (grepl(".ODF$", filename, ignore.case = TRUE)) {
# in BIO files, the data type seems to be on line 14. Read more, for safety.
lines <- readLines(file, encoding = "latin1")
dt <- try(
{
grep("DATA_TYPE[ \t]*=", lines, perl = TRUE)
},
silent = TRUE
)
if (inherits(dt, "try-error") || length(t) < 1) {
stop("cannot infer type of ODF file")
}
subtype <- gsub("[',]", "", tolower(strsplit(lines[dt[1]], "=")[[1]][2]))
subtype <- gsub("^\\s*", "", subtype)
subtype <- gsub("\\s*$", "", subtype)
res <- paste(subtype, "odf", sep = "/")
oceDebug(debug, paste0("} # oceMagic() returning ", res, "\n"), unindent = 1, style = "bold")
return(res)
}
if (grepl(".WCT$", filename, ignore.case = TRUE)) {
# old-style WOCE
oceDebug(debug, "} # oceMagic returning ctd/woce/other\n", unindent = 1)
return("ctd/woce/other") # e.g. http://cchdo.ucsd.edu/data/onetime/atlantic/a01/a01e/a01ect.zip
}
if (grepl(".nc$", filename, ignore.case = TRUE)) {
# argo or netcdf?
if (requireNamespace("ncdf4", quietly = TRUE)) {
if (substr(filename, 1, 5) == "http:") {
stop(
"cannot open netcdf files over the web; try doing as follows\n download.file(\"",
filename, "\", \"", gsub(".*/", "", filename), "\")"
)
}
# NOTE: need to name ncdf4 package because otherwise R checks give warnings.
f <- ncdf4::nc_open(filename)
if ("DATA_TYPE" %in% names(f$var)) {
if (grepl("argo", ncdf4::ncvar_get(f, "DATA_TYPE"), ignore.case = TRUE)) {
oceDebug(debug, "} # oceMagic returning argo (upper-case style)\n", unindent = 1)
ncdf4::nc_close(f)
return("argo")
} else {
oceDebug(debug, "} # oceMagic returning netcdf (upper-case style)\n", unindent = 1)
ncdf4::nc_close(f)
return("netcdf")
}
} else if ("data_type" %in% names(f$var)) {
if (grepl("argo", ncdf4::ncvar_get(f, "data_type"), ignore.case = TRUE)) {
oceDebug(debug, "} # oceMagic returning argo (lower-case style)\n", unindent = 1)
ncdf4::nc_close(f)
return("argo")
} else {
oceDebug(debug, "} # oceMagic returning netcdf (lower-case style)\n", unindent = 1)
ncdf4::nc_close(f)
return("netcdf")
}
}
ncdf4::nc_close(f) # it's netcdf, but with no data_type
return("unknown")
} else {
stop("must install.packages(\"ncdf4\") to read a NetCDF file")
}
}
if (grepl(".xml$", filename, ignore.case = TRUE)) {
firstLine <- readLines(filename, 1L)
if (grepl(".weather.gc.ca", firstLine)) {
oceDebug(debug, "} # oceMagic returning met/xml2\n", unindent = 1)
return("met/xml2")
}
}
if (grepl(".osm.xml$", filename, ignore.case = TRUE)) {
oceDebug(debug, "} # oceMagic returning openstreetmap (xml style)\n", unindent = 1)
return("openstreetmap")
}
if (grepl(".osm$", filename, ignore.case = TRUE)) {
oceDebug(debug, "} # oceMagic returning openstreetmap (non xml style)\n", unindent = 1, style = "bold")
return("openstreetmap")
}
if (grepl(".gpx$", filename, ignore.case = TRUE)) {
oceDebug(debug, "} # oceMagic returning gpx (e.g. Garmin GPS data)\n", unindent = 1, style = "bold")
return("gpx")
}
if (grepl(".csv$", filename, ignore.case = TRUE)) {
con <- file(filename, "r", encoding = encoding)
someLines <- readLines(con, 30L) # , encoding="UTF-8-BOM")
close(con)
# print(someLines[1])
if (grepl("^SSDA Sea & Sun Technology", someLines[1])) {
return("ctd/ssda")
} else if (1L == length(grep("^.*\"WMO Identifier\",", someLines))) {
oceDebug(debug, "} # oceMagic returning met/csv1\n", unindent = 1, style = "bold")
return("met/csv1") # FIXME: may be other things too ...
} else if (grepl("^.*Longitude.*Latitude.*Station Name.*Climate ID.*Dew Point", someLines[1])) {
oceDebug(debug, "} # oceMagic returning met/csv2 or met/csv3\n", unindent = 1, style = "bold")
if (grepl("Time \\(LST\\)", someLines[1])) {
oceDebug(debug, "} # oceMagic returning met/csv2\n", unindent = 1, style = "bold")
return("met/csv3")
} else {
oceDebug(debug, "} # oceMagic returning met/csv3\n", unindent = 1, style = "bold")
return("met/csv2")
}
} else if (length(grep("^Station_Name,", someLines))) {
oceDebug(debug, "} # oceMagic returning sealevel\n", unindent = 1, style = "bold")
return("sealevel")
} else if (1L == length(grep("^CTD,", someLines))) {
oceDebug(debug, "} # oceMagic returning ctd/woce/exchange\n", unindent = 1, style = "bold")
return("ctd/woce/exchange")
} else if (1L == length(grep("^BOTTLE,", someLines))) {
oceDebug(debug, "} # oceMagic returning section\n", unindent = 1, style = "bold")
return("section")
} else {
return("unknown")
}
}
if (grepl(".edf$", filename, ignore.case = TRUE)) {
oceDebug(debug, "} # oceMagic returning xbt/edf\n", unindent = 1, style = "bold")
return("xbt/edf")
}
file <- file(file, "r", encoding = encoding)
}
if (!inherits(file, "connection")) {
stop("argument 'file' must be a character string or connection")
}
oceDebug(debug, "'file' is a connection\n")
if (!isOpen(file)) {
open(file, "r", encoding = encoding)
}
# Grab text at start of file.
lines <- readLines(file, n = 2L, skipNul = TRUE)
line <- lines[1]
line2 <- lines[2]
oceDebug(debug, "first line of file: ", line, "\n", sep = "")
oceDebug(debug, "second line of file: ", line2, "\n", sep = "")
close(file)
file <- file(filename, "rb")
bytes <- readBin(file, what = "raw", n = 4)
oceDebug(debug, paste("first two bytes in file: 0x", bytes[1], " and 0x", bytes[2], "\n", sep = ""))
on.exit(close(file))
if (bytes[1] == 0x00 && bytes[2] == 0x00 && bytes[3] == 0x27 && bytes[4] == 0x0a) {
oceDebug(debug, "this is a shapefile; see e.g. http://en.wikipedia.org/wiki/Shapefile\n }\n")
oceDebug(debug, "} # oceMagic returning shapefile\n", unindent = 1, style = "bold")
return("shapefile")
}
if (bytes[3] == 0xff && bytes[4] == 0xff) {
oceDebug(debug, "} # oceMagic returning echosounder\n", unindent = 1, style = "bold")
return("echosounder")
}
if (bytes[1] == 0x10 && bytes[2] == 0x02) {
# 'ADPManual v710.pdf' p83
if (96 == readBin(bytes[3:4], "integer", n = 1, size = 2, endian = "little")) {
oceDebug(debug, "this is adp/sontek (4 byte match)\n }\n")
} else {
oceDebug(debug, "this is adp/sontek (2 byte match, but bytes 3 and 4 should become integer 96)\n }\n")
}
oceDebug(debug, "} # oceMagic returning adp/sontek\n", unindent = 1, style = "bold")
return("adp/sontek")
}
if (bytes[1] == 0x7f && bytes[2] == 0x7f) {
oceDebug(debug, "} # oceMagic returning adp/rdi\n", unindent = 1, style = "bold")
return("adp/rdi")
}
if (bytes[1] == 0xa5 && bytes[2] == 0x05) {
# NorTek files require deeper inspection. Here, SIG stands for "System Integrator Guide",
# Dated Jue 2008 (Nortek Doc No PS100-0101-0608)
seek(file, 0)
oceDebug(debug, "This is probably a nortek file of some sort. Reading further to see for sure ...\n")
hardware.configuration <- readBin(file, what = "raw", n = 48) # FIXME: this hard-wiring is repeated elsewhere
if (hardware.configuration[1] != 0xa5 || hardware.configuration[2] != 0x05) {
return("unknown")
}
oceDebug(debug, "hardware.configuration[1:2]", hardware.configuration[1:2], "(expect 0xa5 0x05)\n")
head.configuration <- readBin(file, what = "raw", n = 224)
oceDebug(debug, "head.configuration[1:2]", head.configuration[1:2], "(expect 0xa5 0x04)\n")
if (head.configuration[1] != 0xa5 || head.configuration[2] != 0x04) {
return("unknown")
}
user.configuration <- readBin(file, what = "raw", n = 512)
oceDebug(debug, "user.configuration[1:2]", user.configuration[1:2], "(expect 0xa5 0x00)\n")
if (user.configuration[1] != 0xa5 || user.configuration[2] != 0x00) {
return("unknown")
}
nextTwoBytes <- readBin(file, what = "raw", n = 2)
oceDebug(
debug, "nextTwoBytes:", paste("0x", nextTwoBytes[1], sep = ""),
paste("0x", nextTwoBytes[2], sep = ""), "(e.g. 0x5 0x12 is adv/nortek/vector)\n"
)
if (nextTwoBytes[1] == 0xa5 && nextTwoBytes[2] == 0x12) {
oceDebug(debug, "} # oceMagic returning adv/nortek/vector\n", unindent = 1, style = "bold")
return("adv/nortek/vector")
}
if (nextTwoBytes[1] == 0xa5 && nextTwoBytes[2] == 0x01) {
oceDebug(debug, "these two bytes imply this is adp/nortek/aqudopp (see system-integrator-manual_jan2011.pdf Table 5.2)\n")
oceDebug(debug, "} # oceMagic returning adp/nortek/aquadopp\n", unindent = 1, style = "bold")
return("adp/nortek/aquadopp")
}
if (nextTwoBytes[1] == 0xa5 && nextTwoBytes[2] == 0x81) {
oceDebug(debug, "these two bytes imply this is adp/nortek/aqudopp (see N3015-023-Integrators-Guide-Classic_1220.pdf page 30)\n")
oceDebug(debug, "} # oceMagic returning adp/nortek/aquadoppPlusMagnetometer\n", unindent = 1, style = "bold")
return("adp/nortek/aquadoppPlusMagnetometer")
}
if (nextTwoBytes[1] == 0xa5 && nextTwoBytes[2] == 0x21) {
oceDebug(debug, "} # oceMagic returning adp/nortek/aquadoppProfiler\n", unindent = 1, style = "bold")
return("adp/nortek/aquadoppProfiler") # p37 SIG
}
if (nextTwoBytes[1] == 0xa5 && nextTwoBytes[2] == 0x2a) {
oceDebug(debug, "} # oceMagic returning adp/nortek/aquadoppHR\n", unindent = 1, style = "bold")
return("adp/nortek/aquadoppHR") # p38 SIG
}
stop("some sort of nortek ... two bytes are 0x", nextTwoBytes[1], " and 0x", nextTwoBytes[2], " but cannot figure out what the type is")
}
if (bytes[1] == 0xa5 && bytes[4] == 0x10) {
oceDebug(debug, "} # oceMagic returning adp/nortek/ad2cp\n", unindent = 1)
return("adp/nortek/ad2cp")
}
if (bytes[1] == 0x9b && bytes[2] == 0x00) {
warning(paste(
"Possibly this is an RDI CTD file. Oce cannot read such files yet, because\n",
" the author has not located file-format documents. If you get such documents\n",
" from RDI, please send them to dan.kelley@dal.ca so the format can be added."
))
return("possibly RDI CTD")
}
if (1 == length(grep("^CTD", line))) {
oceDebug(debug, "} # oceMagic returning ctd/woce/exchange\n", unindent = 1, style = "bold")
return("ctd/woce/exchange")
}
if (1 == length(grep("^EXPOCODE", line))) {
oceDebug(debug, "} # oceMagic returning ctd/woce/other\n", unindent = 1, style = "bold")
return("ctd/woce/other")
}
if (1 == length(grep("^\\s*ODF_HEADER", line))) {
oceDebug(debug, "} # oceMagic returning odf\n", unindent = 1, style = "bold")
return("odf")
}
if (grepl("^\\* Sea-Bird SBE", line) || grepl("^\\* Viking Buoy CTD file", line)) {
oceDebug(debug, "} # oceMagic returning ctd/sbe\n", unindent = 1, style = "bold")
return("ctd/sbe")
}
if (1 == length(grep("^%ITP", line))) {
oceDebug(debug, "} # oceMagic returning ctd/itp\n", unindent = 1, style = "bold")
return("ctd/itp")
}
if (1 == length(grep("^# -b", line))) {
oceDebug(debug, "} # oceMagic returning coastline\n", unindent = 1, style = "bold")
return("coastline")
}
if (1 == length(grep("^# Station_Name,", line))) {
oceDebug(debug, "} # oceMagic returning sealevel\n", unindent = 1, style = "bold")
return("sealevel")
}
if (1 == length(grep("^Station_Name,", line))) {
oceDebug(debug, "} # oceMagic returning sealevel\n", unindent = 1, style = "bold")
return("sealevel")
}
if (1 == length(grep("^[0-9][0-9][0-9][A-Z] ", line))) {
oceDebug(debug, "} # oceMagic returning sealevel\n", unindent = 1, style = "bold")
return("sealevel")
}
if (1 == length(grep("^NCOLS[ ]*[0-9]*[ ]*$", line))) {
oceDebug(debug, "} # oceMagic returning topo\n", unindent = 1, style = "bold")
return("topo")
}
if (1 == length(grep("^RBR TDR", line))) {
# FIXME: obsolete; to be removed Fall 2015
oceDebug(debug, "} # oceMagic returning RBR/dat\n", unindent = 1, style = "bold")
return("RBR/dat")
}
if (1 == length(grep("^Model=", line))) {
oceDebug(debug, "} # oceMagic returning RBR/txt\n", unindent = 1, style = "bold")
return("RBR/txt")
}
if (1 == length(grep("^BOTTLE", line))) {
oceDebug(debug, "} # oceMagic returning section\n", unindent = 1, style = "bold")
return("section")
}
oceDebug(debug, "this is unknown\n")
return("unknown")
}
#' Read an Oceanographic Data File
#'
#' Read an oceanographic data file, auto-discovering the file type from the
#' first line of the file.
#' This function tries to infer the file type from the first line, using
#' [oceMagic()]. If it can be discovered, then an
#' instrument-specific file reading function is called, with the `file`
#' and with any additional arguments being supplied.
#'
#' @param file a connection or a character string giving the name of the file
#' to load.
#'
#' @param ... arguments to be handed to whichever instrument-specific reading
#' function is selected, based on the header.
#'
#' @param encoding a character string giving the file encoding. This defaults
#' to `"latin1`", which seems to work for files available to the authors, but
#' be aware that a different setting may be required for files that contain
#' unusual accents or characters. (Try `"UTF-8"` if the default produces
#' errors.) Note that `encoding` is ignored in binary files, and also
#' in some text-based files, as well.
#'
#' @return An [oce-class] object of that is
#' specialized to the data type, e.g. [ctd-class],
#' if the data file contains `ctd` data.
#'
#' @author Dan Kelley
#'
#' @seealso The file type is determined by [oceMagic()]. If the file
#' type can be determined, then one of the following is called:
#' [read.ctd()], [read.coastline()]
#' [read.lobo()], [read.rsk()],
#' [read.sealevel()], etc.
#'
#' @examples
#' library(oce)
#' x <- read.oce(system.file("extdata", "ctd.cnv.gz", package = "oce"))
#' plot(x) # summary with TS and profiles
#' plotTS(x) # just the TS
read.oce <- function(file, ..., encoding = "latin1") {
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, "\"")
}
}
dots <- list(...)
debug <- if ("debug" %in% names(dots)) as.integer(dots$debug) else 0L
debug <- max(0L, debug)
if (missing(file)) {
stop("must supply 'file'")
}
if (is.character(file) && "http://" != substr(file, 1, 7) && !file.exists(file)) {
stop("In read.oce() : cannot open '", file, "' because there is no such file or directory", call. = FALSE)
}
if (is.character(file) && "http://" != substr(file, 1, 7) && 0 == file.info(file)$size) {
stop("empty file \"", file, "\"")
}
type <- oceMagic(file, debug = debug - 1)
oceDebug(debug, "read.oce(\"", as.character(file), "\") {\n", sep = "", unindent = 1, style = "bold")
if (is.character(file) && "http://" != substr(file, 1, 7) && 0 == file.info(file)$size) {
stop("empty file \"", file, "\"")
}
#> OLD: deparse is unhelpful if "file" is a variable in the calling code
#> OLD: processingLog <- paste(deparse(match.call()), sep="", collapse="")
processingLog <- paste("read.oce(\"", file, "\"", ifelse(length(list(...)), ", ...)", ")"), sep = "")
if (type == "shapefile") {
res <- read.coastline.shapefile(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "openstreetmap") {
res <- read.coastline.openstreetmap(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "echosounder") {
res <- read.echosounder(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "adp/rdi") {
res <- read.adp.rdi(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "adp/sontek") {
res <- read.adp.sontek(file, encoding = encoding, processingLog = processingLog, ...) # FIXME is pcadcp different?
} else if (type == "adp/nortek/aquadopp") {
res <- read.aquadopp(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "adp/nortek/aquadoppPlusMagnetometer") {
res <- read.aquadopp(file, type = "aquadoppPlusMagnetometer", encoding = encoding, processingLog = processingLog, ...)
} else if (type == "adp/nortek/aquadoppProfiler") {
res <- read.aquadoppProfiler(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "adp/nortek/aquadoppHR") {
res <- read.aquadoppHR(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "adp/nortek/ad2cp") {
res <- read.adp.ad2cp(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "adv/nortek/vector") {
res <- read.adv.nortek(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "adv/sontek/adr") {
res <- read.adv.sontek.adr(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "aml/txt") {
res <- read.ctd.aml(file, encoding = encoding, processingLog = processingLog, ...)
# FIXME need adv/sontek (non adr)
} else if (type == "interocean/s4") {
res <- read.cm.s4(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "ctd/sbe") {
res <- read.ctd.sbe(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "ctd/woce/exchange") {
res <- read.ctd.woce(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "ctd/woce/other") {
res <- read.ctd.woce.other(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "ctd/odf" || type == "mctd/odf" || type == "mvctd/odf") {
res <- read.ctd.odf(file, encoding = encoding, processingLog = processingLog, ...)
} else if (length(grep(".odf$", type))) {
res <- read.odf(file, encoding = encoding, ...)
} else if (type == "mtg/odf") {
# FIXME: document this data type
# Moored tide gauge: returns a data frame.
fromHeader <- function(key) {
i <- grep(key, lines)
if (length(i) < 1) {
""
} else {
gsub("\\s*$", "", gsub("^\\s*", "", gsub("'", "", gsub(",", "", strsplit(lines[i[1]], "=")[[1]][2]))))
}
}
lines <- readLines(file)
nlines <- length(lines)
headerEnd <- grep("-- DATA --", lines)
if (1 != length(headerEnd)) {
stop("found zero or multiple '-- DATA --' (end of header) lines in a mtg/odf file")
}
# header <- lines[1:headerEnd]
data <- lines[seq.int(headerEnd + 1, nlines)]
res <- read.table(text = data, header = FALSE, col.names = c("time", "temperature", "ptotal", "psea", "depth"), encoding = encoding)
res$time <- strptime(res$time, "%d-%B-%Y %H:%M:%S", tz = "UTC") # guess on timezone
missing_value <- -99.0 # FIXME: it's different for each column
res[res == missing_value] <- NA
attr(res, "scientist") <- fromHeader("CHIEF_SCIENTIST")
attr(res, "latitude") <- as.numeric(fromHeader("INITIAL_LATITUDE"))
attr(res, "longitude") <- as.numeric(fromHeader("INITIAL_LONGITUDE"))
attr(res, "cruise_name") <- fromHeader("CRUISE_NAME")
attr(res, "cruise_description") <- fromHeader("CRUISE_DESCRIPTION")
attr(res, "inst_type") <- fromHeader("INST_TYPE")
attr(res, "model") <- fromHeader("MODEL")
attr(res, "serial_number") <- fromHeader("SERIAL_NUMBER")
attr(res, "missing_value") <- missing_value
warning("Missing-value code for mtg/odf is hard-wired to -99, which will likely be wrong in other files")
warning("The format of mtg/odf objects is likely to change throughout April, 2015")
# if (type == "ctd/odv")
# return(read.ctd.odv(file, processingLog=processingLog, ...))
} else if (type == "ctd/itp") {
res <- read.ctd.itp(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "ctd/ssda") {
res <- read.ctd.ssda(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "gpx") {
res <- read.gps(file, type = "gpx", encoding = encoding, processingLog = processingLog, ...)
} else if (type == "coastline") {
res <- read.coastline(file, type = "mapgen", encoding = encoding, processingLog = processingLog, ...)
} else if (type == "argo") {
res <- read.argo(file, encoding = encoding, ...)
} else if (type == "lisst") {
res <- read.lisst(file, encoding = encoding)
} else if (type == "sealevel") {
res <- read.sealevel(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "topo") {
res <- read.topo(file)
} else if (type == "RBR/dat") { # FIXME: obsolete; to be removed by Fall 2015
res <- read.rsk(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "RBR/rsk") {
res <- read.rsk(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "RBR/txt") {
res <- read.rsk(file, encoding = encoding, processingLog = processingLog, type = "txt", ...)
} else if (type == "section") {
res <- read.section(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "ctd/woce/other") {
res <- read.ctd.woce.other(file, encoding = encoding, processingLog = processingLog, ...)
} else if (type == "landsat") {
res <- read.landsat(file, encoding = encoding, ...)
} else if (type == "netcdf") {
res <- read.netcdf(file, encoding = encoding, ...)
} else if (type == "met/csv1") {
res <- read.met(file, type = "csv1", encoding = encoding, ...)
} else if (type == "met/csv2") {
res <- read.met(file, type = "csv2", encoding = encoding, ...)
} else if (type == "met/csv3") {
res <- read.met(file, type = "csv2", encoding = encoding, ...) # FIXME: ok?
} else if (type == "met/xml2") {
res <- read.met(file, type = "xml2", encoding = encoding, ...)
} else if (type == "odf") {
res <- read.odf(file, encoding = encoding, ...)
} else if (type == "xbt/edf") {
res <- read.xbt.edf(file, encoding = encoding, ...)
} else {
stop("unknown file type \"", type, "\"")
}
oceDebug(debug, "} # read.oce()\n", unindent = 1, sep = "", style = "bold")
res
}
#' Draw an Axis, Possibly with Decade-style Logarithmic Scaling
#'
#' @param logStyle a character value that indicates how to draw the y axis, if
#' `log="y"`. If it is `"r"` (the default) then the conventional R style is used,
#' in which a logarithmic transform connects y values to position on the "page"
#' of the plot device, so that tics will be nonlinearly spaced, but not
#' organized by integral powers of 10. However, if it is `"decade"`, then
#' the style will be that used in the scientific literature, in which large
#' tick marks are used for integral powers of 10, with smaller tick marks
#' at integral multiples of those powers, and with labels that use exponential
#' format for values above 100 or below 0.01.
#' @param side an integer specifying which axis to draw, with 1 for bottom axis, 2 for left axis,
#' 3 for top axis, and 4 for right axis (as with [axis()]).
#' @param labels either a vector of character values used for labels or a logical value indicating
#' whether to draw such labels. The first form only works if the coordinate is not logarithmic,
#' and if `logStyle` is `"r"`.
#' @param \dots other graphical parameters, passed to [axis()].
#'
#' @return Numerical values at which tick marks were drawn (or would have been drawn, if `labels`
#' specified to draw them).
#'
#' @examples
#' library(oce)
#' Ra <- 10^seq(4, 10, 0.1)
#' Nu <- 0.085 * Ra^(1 / 3)
#' plot(Ra, Nu, log = "xy", axes = FALSE)
#' box()
#' oceAxis(1, logStyle = "decade")
#' oceAxis(2, logStyle = "decade")
#'
#' @author Dan Kelley
oceAxis <- function(side, labels = TRUE, logStyle = "r", ...) {
if (missing(side)) {
stop("in oceAxis() :\n argument \"side\" is missing, with no default", call. = FALSE)
}
if (length(side) != 1) {
stop("in oceAxis() :\n argument \"side\" must be a single number", call. = FALSE)
}
if (!(side %in% 1:4)) {
stop("in oceAxis() :\n argument \"side\" must be 1, 2, 3 or 4", call. = FALSE)
}
if (!(logStyle %in% c("r", "decade"))) {
stop("logStyle must be \"r\" or \"decade\", not \"", logStyle, "\"")
}
if (logStyle == "r") {
return(invisible(axis(side = side, labels = labels, ...)))
} else {
# use decade axis if previous plot() call made this coordinate be logarithmic
if (((side %in% c(1, 3)) && par("xlog")) || ((side %in% c(2, 4)) && par("ylog"))) {
usr <- if (side %in% c(1, 3)) par("usr")[1:2] else par("usr")[3:4]
lowerDecade <- floor(usr[1])
upperDecade <- floor(1 + usr[2])
smallTickAt <- NULL
bigTickAt <- NULL
bigTickLabel <- NULL
for (bigTick in lowerDecade:upperDecade) {
bigTickAt <- c(bigTickAt, bigTick)
bigTickLabel <- c(
bigTickLabel,
if (bigTick < -1L || bigTick > 1L) {
substitute(10^A, list(A = bigTick))
} else {
10^bigTick
}
)
smallTickAt <- c(smallTickAt, -1 + bigTick + log10(2:9))
smallTickAt <- c(smallTickAt, bigTick + log10(2:9))
}
bigTickInWindow <- usr[1] <= bigTickAt & bigTickAt <= usr[2]
bigTickAt <- bigTickAt[bigTickInWindow]
bigTickLabel <- as.expression(bigTickLabel[bigTickInWindow])
smallTickInWindow <- usr[1] <= smallTickAt & smallTickAt <= usr[2]
smallTickAt <- smallTickAt[smallTickInWindow]
rval <- axis(side = side, at = 10^bigTickAt, labels = if (labels) bigTickLabel else FALSE)
rug(side = side, x = 10^smallTickAt, tcl = 0.5 * par("tcl"), lwd = par("lwd"))
return(invisible(rval))
} else {
return(invisible(axis(side = side, labels = labels, ...)))
}
}
}
#' Create Two-Color Palette
#'
#' Create colors ranging between two specified limits, with white
#' in the middle.
#'
#' @aliases oceColorsTwo oce.colorsTwo
#'
#' @param n number of colors to generate.
#'
#' @param low,high numerical values (in range 0 to 1) specifying the hue
#' for the low and high ends of the color scale.
#'
#' @param smax numerical value (in range 0 to 1) for the color saturation.
#'
#' @param alpha numerical value (in ragne 0 to 1) for the alpha (transparency)
#' of the colors.
#'
#' @examples
#' library(oce)
#' imagep(volcano - mean(range(volcano)),
#' col = oceColorsTwo(128),
#' zlim = "symmetric", zlab = "oceColorsTwo"
#' )
#' @family things related to colors
oceColorsTwo <- function(n, low = 2 / 3, high = 0, smax = 1, alpha = 1) {
# code borrows heavily from cm.color()
if ((n <- as.integer(n[1])) > 0) {
even.n <- n %% 2 == 0
k <- n %/% 2
l1 <- k + 1 - even.n
l2 <- n - k + even.n
c(
if (l1 > 0) {
hsv(
h = low,
s = seq.int(smax, ifelse(even.n, 0.5 / k, 0), length.out = l1),
v = 1,
alpha = alpha
)
},
if (l2 > 1) {
hsv(
h = high,
s = seq.int(0, smax, length.out = l2)[-1],
v = 1,
alpha = alpha
)
}
)
} else {
character(0)
}
}
oce.colorsTwo <- oceColorsTwo
#' Create Colors in a GEBCO-like Scheme
#'
#' The colours were determined by examination of paper
#' charts printed during the GEBCO Fifth Edition era.
#' The hues range from dark blue to light blue, then
#' from light brown to dark brown. If used to show
#' topography in scheme centred on z=0, this means that
#' near-coastal regions are light in tone, with darker
#' colours representing both mountains and the deep sea.
#'
#' @aliases oceColorsGebco oce.colorsGebco
#'
#' @param n Number of colors to return
#'
#' @param region String indicating application region,
#' one of `"water"`, `"land"`, or `"both"`.
#'
#' @param type String indicating the purpose, one of `"fill"` or `"line"`.
#'
#' @param debug a flag that turns on debugging.
#'
#' @examples
#' library(oce)
#' imagep(volcano, col = oceColorsGebco(128, region = "both"))
#'
#' @family things related to colors
oceColorsGebco <- function(n = 9, region = c("water", "land", "both"), type = c("fill", "line"), debug = getOption("oceDebug")) {
oceDebug(debug, "oceColorsGebco(n=", n, ", region=\"", region, "\", type=\"", type, "\", debug=", debug, ")\n", sep = "", unindent = 1)
region <- match.arg(region)
type <- match.arg(type)
land <- c(
"#FEF1E0", "#FDE3C1", "#FBC784", "#F1C37A", "#E6B670", "#DCA865",
"#D19A5C", "#C79652", "#BD9248", "#B38E3E"
)
water <- c(
"#0F7CAB", "#2292B5", "#38A7BF", "#4FBBC9", "#68CDD4", "#83DEDE",
"#A0E8E4", "#BFF2EC", "#E1FCF7", "#F0FDFB"
)
if (type == "fill") {
# generate land colors by e.g. rgb(t(col2rgb(land[5])-1*c(10, 4, 10))/255)
# until 2020-12-14 land <- c("#FBC784", "#F1C37A", "#E6B670", "#DCA865", "#D19A5C",
# until 2020-12-14 "#C79652", "#BD9248", "#B38E3E", "#A98A34")
# until 2020-12-14 water <- rev(c("#E1FCF7", "#BFF2EC", "#A0E8E4", "#83DEDE", "#68CDD4",
# until 2020-12-14 "#4FBBC9", "#38A7BF", "#2292B5", "#0F7CAB"))
# land <- c("#FEF1E0", "#FDE3C1", "#FBC784", "#F1C37A", "#E6B670", "#DCA865",
# "#D19A5C", "#C79652", "#BD9248", "#B38E3E")
# water <- c("#0F7CAB", "#2292B5", "#38A7BF", "#4FBBC9", "#68CDD4", "#83DEDE",
# "#A0E8E4", "#BFF2EC", "#E1FCF7", "#F0FDFB")
land <- c(
"#FFF0DF", "#FFE9D0", "#FFE2C1", "#FDD6A6", "#FBC98A", "#F7C580", "#F2C37B", "#EDBE76", "#E8B872", "#E3B26D",
"#DEAB67", "#D9A563", "#D49E5E", "#CF995A", "#CA9755", "#C59550", "#C1934C", "#BC9147", "#B78F42", "#B38E3E"
)
water <- c(
"#0F7CAB", "#1886AF", "#2090B4", "#2B9AB9", "#35A4BD", "#40AEC2", "#4BB7C7", "#56C0CC", "#62C9D1", "#6FD1D6",
"#7BD9DB", "#89E0DF", "#96E4E2", "#A4E9E5", "#B3EEE9", "#C2F3ED", "#D2F7F2", "#E2FCF7", "#EBFDF9", "#F5FEFC"
)
} else {
oceDebug(debug, "type=\"line\"\n")
land <- c(
"#FBC784", "#F1C37A", "#E6B670", "#DCA865", "#D19A5C",
"#C79652", "#BD9248", "#B38E3E", "#A98A34"
)
water <- rev(c(
"#A4FCE3", "#72EFE9", "#4FE3ED", "#47DCF2", "#46D7F6",
"#3FC0DF", "#3FC0DF", "#3BB7D3", "#36A5C3"
))
}
if (region == "water") {
rgb.list <- col2rgb(water) / 255
l <- length(water)
r <- approx(1:l, rgb.list[1, 1:l], xout = seq(1, l, length.out = n))$y
g <- approx(1:l, rgb.list[2, 1:l], xout = seq(1, l, length.out = n))$y
b <- approx(1:l, rgb.list[3, 1:l], xout = seq(1, l, length.out = n))$y
res <- rgb(r, g, b)
} else if (region == "land") {
rgb.list <- col2rgb(land) / 255
l <- length(land)
r <- approx(1:l, rgb.list[1, 1:l], xout = seq(1, l, length.out = n))$y
g <- approx(1:l, rgb.list[2, 1:l], xout = seq(1, l, length.out = n))$y
b <- approx(1:l, rgb.list[3, 1:l], xout = seq(1, l, length.out = n))$y
res <- rgb(r, g, b)
} else {
# both
# See https://github.com/dankelley/oce/discussions/1756#discussioncomment-204754 for
# a discussion of adding some white 'ink' between the water and the land.
# ? rgb.list <- col2rgb(c(water, "#FFFFFF", "#FFFFFF", land)) / 255
# ? rgb.list <- col2rgb(c(water, "#FFFFFF", land)) / 255
# 20201214> rgb.list <- col2rgb(c(water, "#FFFFFF", land)) / 255
# 20201214> l <- ncol(rgb.list)
# 20201214> r <- approx(1:l, rgb.list[1, 1:l], xout=seq(1, l, length.out=n))$y
# 20201214> g <- approx(1:l, rgb.list[2, 1:l], xout=seq(1, l, length.out=n))$y
# 20201214> b <- approx(1:l, rgb.list[3, 1:l], xout=seq(1, l, length.out=n))$y
# I find it very difficult to see a difference between 'rgb' and 'Lab' spaces, and between
# 'linear' and 'spline' interpolations.
cr <- colorRamp(c(water, "#FFFFFF", land), bias = 1, space = "rgb", interpolate = "spline")(seq(0, 1, length.out = n))
res <- rgb(cr, maxColorValue = 255)
}
oceDebug(debug, "} # oceColorsGebco()", sep = "", unindent = 1)
res
}
oce.colorsGebco <- oceColorsGebco
#' Create Color Functions
#'
#' This function generates other functions that are used to specify colors.
#' It is used within oce to create [oceColorsTemperature()]
#' and its many cousins. Users may also find it helpful, for creating
#' custom color schemes (see \dQuote{Examples}).
#'
#' @param spec Specification of the color scheme. This may be a
#' character string, in which case it must be the name of an item stored
#' in `data(ocecolors)`, or either a 3-column data frame or
#' matrix, in which case the columns specify red, green and blue values
#' (in range from 0 to 1).
#'
#' @section Sample of Usage:
#' \preformatted{
#' # Update oxygen color scheme to latest matplotlib value.
#' library(oce)
#' oxy <- "https://raw.githubusercontent.com/matplotlib/cmocean/master/cmocean/rgb/oxy-rgb.txt"
#' oxyrgb <- read.table(oxy, header=FALSE)
#' oceColorsOxygenUpdated <- oceColorsClosure(oxyrgb)
#' par(mfrow=c(1, 2))
#' m <- matrix(1:256)
#' imagep(m, col=oceColorsOxygen, zlab="oxygen")
#' imagep(m, col=oceColorsOxygenUpdated, zlab="oxygenUpdated")
#' }
#' @family things related to colors
oceColorsClosure <- function(spec) {
function(n) {
if (is.character(spec)) {
data("ocecolors", package = "oce", envir = environment())
col <- get("ocecolors")[[spec]]
} else if (is.data.frame(spec) || is.matrix(spec)) {
col <- rgb(spec[, 1], spec[, 2], spec[, 3])
} else {
stop("oceColorsClosure(): first arg must be character, data frame, or 3-column matrix", call. = FALSE)
}
if (missing(n) || n <= 0) colorRampPalette(col) else colorRampPalette(col)(n)
}
}
#' Create Colors Similar to the Google Turbo Scheme
#'
#' This uses the coefficients published (with Apache license) by google,
#' as described by Mikhailo (2019).
#'
#' @aliases oce.colorsTurbo oceColorsTurbo
#'
#' @param n number of colors to create.
#'
#' @references
#' Mikhailo, Anton.
#' \dQuote{Turbo, An Improved Rainbow Colormap for Visualization.}
#' Google AI (blog), August 20, 2019.
#' `https://ai.googleblog.com/2019/08/turbo-improved-rainbow-colormap-for.html`
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' imagep(volcano,
#' col = oceColorsTurbo(128),
#' zlab = "oceColorsTurbo"
#' )
#'
#' @template colourBlindnessTemplate
#' @family things related to colors
oceColorsTurbo <- oce.colorsTurbo <- oceColorsClosure("turbo")
#' Create Colors Similar to the Matlab Viridis Scheme
#'
#' This is patterned on a \proglang{matlab}/\proglang{python} scheme that blends
#' from yellow to blue in a way that is designed to reproduce well
#' in black-and-white, and to be interpretable by those with
#' certain forms of color blindness. See the references for
#' notes about issues of colour blindness in computer graphics.
#' An alternative
#' to [oceColorsViridis] is provided in the \CRANpkg{viridis} package, as illustrated
#' in Example 2.
#'
#' @aliases oce.colorsViridis oceColorsViridis
#'
#' @param n number of colors to create.
#'
#' @template colourBlindnessTemplate
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' # Example 1: oceColorsViridis
#' imagep(volcano,
#' col = oceColorsViridis(128),
#' zlab = "oceColorsViridis"
#' )
#'
#' @family things related to colors
#'
#' @template colourBlindnessTemplate
oceColorsViridis <- oce.colorsViridis <- oceColorsClosure("viridis")
#' @templateVar colorItem CDOM
#' @templateVar colorItemUC CDOM
#' @template cmcolorTemplate
#'
#' @template colourBlindnessTemplate
oceColorsCDOM <- oce.colorsCDOM <- oceColorsClosure("cdom")
#' @templateVar colorItem chlorophyll
#' @templateVar colorItemUC Chlorophyll
#' @template cmcolorTemplate
#'
#' @template colourBlindnessTemplate
oceColorsChlorophyll <- oce.colorsChlorophyll <- oceColorsClosure("chlorophyll")
#' @templateVar colorItem density
#' @templateVar colorItemUC Density
#' @template cmcolorTemplate
#'
#' @template colourBlindnessTemplate
oceColorsDensity <- oce.colorsDensity <- oceColorsClosure("density")
#' @templateVar colorItem freesurface
#' @templateVar colorItemUC Freesurface
#' @template cmcolorTemplate
#'
#' @template colourBlindnessTemplate
oceColorsFreesurface <- oce.colorsFreesurface <- oceColorsClosure("freesurface")
#' @templateVar colorItem oxygen
#' @templateVar colorItemUC Oxygen
#' @template cmcolorTemplate
#'
#' @template colourBlindnessTemplate
oceColorsOxygen <- oce.colorsOxygen <- oceColorsClosure("oxygen")
#' @templateVar colorItem PAR
#' @templateVar colorItemUC PAR
#' @template cmcolorTemplate
#'
#' @template colourBlindnessTemplate
oceColorsPAR <- oce.colorsPAR <- oceColorsClosure("par")
#' @templateVar colorItem phase
#' @templateVar colorItemUC Phase
#' @template cmcolorTemplate
#'
#' @template colourBlindnessTemplate
oceColorsPhase <- oce.colorsPhase <- oceColorsClosure("phase")
#' @templateVar colorItem salinity
#' @templateVar colorItemUC Salinity
#' @template cmcolorTemplate
#'
#' @template colourBlindnessTemplate
oceColorsSalinity <- oce.colorsSalinity <- oceColorsClosure("salinity")
#' @templateVar colorItem temperature
#' @templateVar colorItemUC Temperature
#' @template cmcolorTemplate
#'
#' @template colourBlindnessTemplate
oceColorsTemperature <- oce.colorsTemperature <- oceColorsClosure("temperature")
#' @templateVar colorItem turbidity
#' @templateVar colorItemUC Turbidity
#' @template cmcolorTemplate
#'
#' @template colourBlindnessTemplate
oceColorsTurbidity <- oce.colorsTurbidity <- oceColorsClosure("turbidity")
#' @templateVar colorItem velocity
#' @templateVar colorItemUC Velocity
#' @template cmcolorTemplate
#'
#' @template colourBlindnessTemplate
oceColorsVelocity <- oce.colorsVelocity <- oceColorsClosure("velocity")
#' @templateVar colorItem vorticity
#' @templateVar colorItemUC Vorticity
#' @template cmcolorTemplate
#'
#' @template colourBlindnessTemplate
oceColorsVorticity <- oce.colorsVorticity <- oceColorsClosure("vorticity")
#' Create Colors Similar to the Matlab Jet Scheme
#' @aliases oceColorsJet oce.colorsJet oceColors9A oce.colors9A
#' @param n number of colors
#' @examples
#' library(oce)
#' imagep(volcano, col = oceColorsJet, zlab = "oceColorsJet")
#'
#' @template colourBlindnessTemplate
#' @family things related to colors
oceColorsJet <- function(n) {
if (missing(n) || n <= 0) {
colorRampPalette(c(
"#00007F", "blue", "#007FFF", "cyan",
"#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"
))
} else {
colorRampPalette(c(
"#00007F", "blue", "#007FFF", "cyan",
"#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"
))(n)
}
}
oce.colors9A <- oceColors9A <- oce.colorsJet <- oceColorsJet
#' Create Colors in a Red-Yellow-Blue Color Scheme
#'
#' The results are similar to those of [oceColorsJet()], but
#' with white hues in the centre, rather than green ones. The scheme
#' may be useful in displaying signed quantities, and thus is somewhat
#' analogous to [oceColorsTwo()], except that some viewers
#' may be able to distinguish more colors with `oceColors9B`.
#'
#' @aliases oceColors9B oce.colors9B
#'
#' @param n number of colors
#'
#' @examples
#' library(oce)
#' imagep(volcano,
#' col = oceColors9B(128),
#' zlab = "oceColors9B"
#' )
#'
#' @template colourBlindnessTemplate
#' @family things related to colors
oceColors9B <- function(n) {
if (missing(n) || n <= 0) {
colorRampPalette(c(
"#00007F", "blue", "#007FFF", "#22e4e7",
"white", "#ffe45e", "#FF7F00", "red", "#7F0000"
))
} else {
colorRampPalette(c(
"#00007F", "blue", "#007FFF", "#22e4e7",
"white", "#ffe45e", "#FF7F00", "red", "#7F0000"
))(n)
}
}
oce.colors9B <- oceColors9B
#' Create a Vector of Colors
#'
#' The available schemes are:
#' * `which=1` for a red-white-blue scheme.
#' * `which=2` for a red-yellow-blue scheme.
#' * `which=9.01`, `which="9A"` or `which="jet"` for [`oceColorsJet`]`(n)`.
#' * `which=9.02` or `which="9B"` for [`oceColors9B`]`(n)`.
#'
#' @param n number of colors to create
#'
#' @param which integer or character string indicating the palette
#' to use; see \dQuote{Details}.
#'
#' @aliases oce.colorsPalette oceColorsPalette
#'
#' @template colourBlindnessTemplate
#'
#' @family things related to colors
oceColorsPalette <- function(n, which = 1) {
if ((n <- as.integer(n[1])) > 0) {
if (which == 1) {
# Started with http://www.personal.psu.edu/cab38/ColorBrewer/ColorBrewer.html
# RdBu 11 divisions
# and then smoothed the components with smooth.spline(...,df=6)
rgb <- matrix(c(
103, 000, 026,
178, 024, 046,
214, 096, 072,
244, 165, 136,
253, 219, 195,
247, 247, 247,
209, 229, 238,
146, 197, 226,
067, 147, 184,
033, 102, 179,
005, 048, 97
), ncol = 3, byrow = TRUE) / 255
m <- dim(rgb)[1]
i <- 1:m
xout <- seq(1, m, length.out = n)
rev(rgb(
red = approx(i, rgb[, 1], xout, rule = 1)$y,
green = approx(i, rgb[, 2], xout, rule = 1)$y,
blue = approx(i, rgb[, 3], xout, rule = 1)$y,
alpha = 1
))
} else if (which == 2) {
# http://www.personal.psu.edu/cab38/ColorBrewer/ColorBrewer.html
m <- 11 # number of classes
r <- c(165, 215, 244, 253, 254, 255, 224, 171, 116, 69, 49) / 255
g <- c(0, 48, 109, 174, 224, 255, 243, 217, 173, 117, 54) / 255
b <- c(38, 39, 67, 97, 144, 191, 248, 233, 209, 180, 149) / 255
i <- 1:m
xout <- seq(1, m, length.out = n)
rev(rgb(
approx(i, r, xout, rule = 1)$y,
approx(i, g, xout, rule = 1)$y,
approx(i, b, xout, rule = 1)$y
))
} else if (which == 9.01 || which == "9A" || which == "jet") {
# jet, also known as 9A or 9.01
oceColorsJet(n)
} else if (which == 9.02 || which == "9B") {
oceColors9B(n)
} else {
stop("unknown which")
}
} else {
character(0)
}
}
oce.colorsPalette <- oceColorsPalette
#' Oce Version of axis.POSIXct
#'
#' A specialized variant of [axis.POSIXct()] that produces
#' results with less ambiguity in axis labels.
#'
#' The tick marks are set automatically based on examination of the time range on
#' the axis. The scheme was devised by constructing test cases with a typical plot
#' size and font size, and over a wide range of time scales. In some categories,
#' both small tick marks are interspersed between large ones.
#'
#' The user may set the format of axis numbers with the `tformat` argument.
#' If this is not supplied, the format is set based on the time span of the axis:
#'
#' * If this time span is less than a minute, the time axis labels are in
#' seconds (fractional seconds, if the interval is less than 2 seconds), with
#' leading zeros on small integers. (Fractional seconds are enabled with a trick:
#' the usual R format `"\%S"` is supplemented with a new format e.g.
#' `"\%.2S"`, meaning to use two digits after the decimal.)
#'
#' * If the time span exceeds a minute but is less than 1.5 days, the label
#' format is `"\%H:\%M:\%S"`.
#'
#' * If the time span exceeds 1.5 days but is less than 1 year, the format is
#' `"\%b \%d"` (e.g. Jul 15) and, again, the tick marks are set up for several
#' subcategories.
#'
#' * If the time span exceeds a year, the format is `"\%Y"`, i.e. the year
#' is displayed with 4 digits.
#'
#' It should be noted that this scheme differs from the R approach in several
#' ways. First, R writes day names for some time ranges, in a convention that is
#' seldom seen in the literature. Second, R will write nn:mm for both HH:MM and
#' MM:SS, an ambiguity that might confuse readers. Third, the use of both large
#' and small tick marks is not something that R does.
#'
#' Bear in mind that `tformat` may be set to alter the number format, but
#' that the tick mark scheme cannot (presently) be controlled.
#'
#' @param side as for [axis.POSIXct()].
#'
#' @param x as for [axis.POSIXct()].
#'
#' @param at as for [axis.POSIXct()].
#'
#' @param tformat as `format` for [axis.POSIXct()] for now, but
#' may eventually have new features for multiline labels, e.g. day on one line
#' and month on another.
#'
#' @param labels as for [axis.POSIXct()].
#'
#' @param drawTimeRange Optional indication of whether/how to draw the time range
#' in the margin on the side of the the plot opposite the time axis. If this is
#' not supplied, it defaults to the value returned by
#' [`getOption`]`("oceDrawTimeRange")`, and if that option is not set,
#' it defaults to `TRUE`. No time range is drawn if `drawTimeRange` is `FALSE`.
#' If it is `TRUE`, the range will be shown. This range refers to
#' range of the x axis (not the data). The format of the elements of that range is set by
#' [`getOption`]`("oceTimeFormat")` (or with the default value
#' of an empty string, if this option has not been set). The timezone will
#' be indicated if the time range is under a week. For preliminary work, it makes
#' sense to use `drawTimeRange=TRUE`, but for published work it can be better
#' to drop this label and indicate something about the time in the figure caption.
#'
#' @param drawFrequency boolean, `TRUE` to show the frequency of sampling
#' in the data
#'
#' @param abbreviateTimeRange boolean, `TRUE` to abbreviate the second
#' number in the time range, e.g. dropping the year if it is the same in the
#' first number.
#'
#' @param cex.axis,cex.lab,cex.main character expansion factors for axis numbers, axis names and plot titles; see [par()].
#'
#' @param mar value for `par(mar)` for axis
#'
#' @param mgp value for `par(mgp)` for axis
#'
#' @param main title of plot
#'
#' @param debug a flag that turns on debugging. Set to 1 to get a moderate
#' amount of debugging information, or to 2 to get more.
#'
#' @param \dots as for [axis.POSIXct()].
#'
#' @return A vector of times corresponding to axis ticks is returned silently.
#'
#' @author Dan Kelley
#'
#' @seealso This is used mainly by [oce.plot.ts()].
oce.axis.POSIXct <- function(
side, x, at, tformat, labels = TRUE,
drawTimeRange,
abbreviateTimeRange = FALSE, drawFrequency = FALSE,
cex.axis = par("cex.axis"), cex.lab = par("cex.lab"), cex.main = par("cex.main"),
mar = par("mar"),
mgp = par("mgp"),
main = "",
debug = getOption("oceDebug"), ...) {
oceDebug(debug, "oce.axis.POSIXct(..., debug=", debug, ",...) {\n", sep = "", unindent = 1, style = "bold")
oceDebug(debug, argShow(mar), "\n")
oceDebug(debug, argShow(mgp), "\n")
oceDebug(debug, "cex.axis=", cex.axis, ", cex.lab=", cex.lab, ", cex.main=", cex.main, "\n")
oceDebug(debug, vectorShow(x, "x"))
tformatGiven <- !missing(tformat)
if (missing(drawTimeRange)) {
drawTimeRange <- getOption("oceDrawTimeRange")
}
# This was written because axis.POSIXt in R version 2.8.x did not obey the
# time zone in the data. (Version 2.9.0 obeys the time zone.)
if (missing(x)) {
x <- numberAsPOSIXct(par("usr")[1:2])
}
dots <- list(...)
if ("xlim" %in% names(dots)) {
ok <- dots$xlim[1] <= x & x <= dots$xlim[2]
x <- x[ok]
}
mat <- missing(at) || is.null(at)
x <- if (!mat) as.POSIXct(at) else as.POSIXct(x)
range <- par("usr")[if (side %% 2) 1:2 else 3:4]
d <- range[2] - range[1] # time span, in seconds
z <- c(range, x[is.finite(x)])
attr(z, "tzone") <- attr(x, "tzone") # need this because c() makes it local time zone (!)
rr <- range + as.POSIXct("2000-01-20") - as.numeric(as.POSIXct("2000-01-20"))
attr(rr, "tzone") <- attr(x, "tzone")
oceDebug(
debug, "range=",
format.POSIXct(rr[1], "%Y-%m-%d %H:%M:%S", tz = "UTC"),
"UTC to ",
format.POSIXct(rr[2], "%Y-%m-%d %H:%M:%S", tz = "UTC"),
"UTC\n"
)
z.sub <- NULL # unlabelled tics may be set in some time ranges, e.g. hours, for few-day plots
oceDebug(debug, "d=", d, " (time range)\n")
if (d <= 2) {
oceDebug(debug, "Time range is under 2 sec\n")
# The time rounding will fail for very small time intervals;
# a wider range can be added easily.
t.start <- rr[1]
t.end <- rr[2]
span <- as.numeric(t.end) - as.numeric(t.start)
if (span > 1) {
round <- 0.5
} else if (span > 0.1) {
round <- 0.05
} else if (span > 0.01) {
round <- 0.005
} else if (span > 0.001) {
round <- 0.0005
} else if (span > 0.0001) {
round <- 0.00005
} else {
round <- 0.00001
}
t0 <- trunc(t.start, "sec")
t.start <- t0 + round * floor((as.numeric(t.start) - as.numeric(t0)) / round)
t.end <- t0 + round * floor((as.numeric(t.end) - as.numeric(t0)) / round)
z <- seq(t.start, t.end, by = round)
oceDebug(debug, vectorShow(z, "TIME RANGE is under 2 seconds; z="))
# BOOKMARK 1A
if (missing(tformat)) {
tformat <- "%.1S" # NOTE: this .1 is interpreted at BOOKMARK 1B
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else if (d <= 20) {
oceDebug(debug, "Time range is between 2 sec and 20 sec\n")
t.start <- trunc(rr[1] - 1, "secs")
t.end <- trunc(rr[2] + 1, "secs")
z <- seq(t.start, t.end, by = "1 sec")
oceDebug(debug, vectorShow(z))
if (missing(tformat)) {
tformat <- "%S"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else if (d <= 60) {
oceDebug(debug, "Time range is between 20 sec and 1 min\n")
t.start <- trunc(rr[1] - 1, "secs")
t.end <- trunc(rr[2] + 1, "secs")
z <- seq(t.start, t.end, by = "2 sec")
oceDebug(debug, vectorShow(z))
if (missing(tformat)) {
tformat <- "%S"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else if (d <= 60 * 3) {
oceDebug(debug, "Time range is between 1 min and 3 min\n")
t.start <- trunc(rr[1] - 60, "mins")
t.end <- trunc(rr[2] + 60, "mins")
z <- seq(t.start, t.end, by = "10 sec")
oceDebug(debug, vectorShow(z))
if (missing(tformat)) {
tformat <- "%M:%S"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else if (d <= 60 * 30) {
oceDebug(debug, "Time range is between 3 min and 30 min\n")
t.start <- trunc(rr[1] - 30, "mins")
t.end <- trunc(rr[2] + 30, "mins")
z <- seq(t.start, t.end, by = "min")
oceDebug(debug, vectorShow(z))
if (missing(tformat)) {
tformat <- "%M:%S"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else if (d <= 60 * 60) {
oceDebug(debug, "Time range is between 30 min and 1 hour\n")
t.start <- trunc(rr[1] - 30, "mins")
t.end <- trunc(rr[2] + 30, "mins")
z <- seq(t.start, t.end, by = "10 min")
oceDebug(debug, vectorShow(z))
if (missing(tformat)) {
tformat <- "%M:%S"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else if (d <= 60 * 60 * 2) {
oceDebug(debug, "Time range is between 1 and 2 hours\n")
t.start <- trunc(rr[1] - 30, "mins")
t.end <- trunc(rr[2] + 30, "mins")
z <- seq(t.start, t.end, by = "10 min")
oceDebug(debug, vectorShow(z))
if (missing(tformat)) {
tformat <- "%H:%M:%S"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else if (d <= 60 * 60 * 6) {
oceDebug(debug, "Time range is between 2 and 6 hours\n")
t.start <- trunc(rr[1], "hour")
t.end <- trunc(rr[2] + 3600, "hour")
z <- seq(t.start, t.end, by = "30 min")
oceDebug(debug, vectorShow(z))
if (missing(tformat)) {
tformat <- "%H:%M:%S"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else if (d <= 60 * 60 * 24 * 1.5) {
oceDebug(debug, "Time range is between 6 hours and 1.5 days\n")
t.start <- trunc(rr[1], "hour")
t.end <- trunc(rr[2] + 86400, "hour")
z <- seq(t.start, t.end, by = "2 hour")
z.sub <- seq(t.start, t.end, by = "hour")
oceDebug(debug, vectorShow(z))
if (missing(tformat)) {
tformat <- "%H:%M:%S"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else if (d <= 60 * 60 * 24 * 5) {
oceDebug(debug, "Time range is between 1.5 and 5 days\n")
t.start <- trunc(rr[1], "day")
t.end <- trunc(rr[2] + 86400, "day")
z <- seq(t.start, t.end, by = "day")
z.sub <- seq(t.start, t.end, by = "6 hour")
oceDebug(debug, vectorShow(z))
oceDebug(debug, vectorShow(z.sub))
if (missing(tformat)) {
tformat <- "%b %d"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else if (d <= 60 * 60 * 24 * 14) {
oceDebug(debug, "Time range is between 4 days and 2 weeks\n")
t.start <- trunc(rr[1], "day")
t.end <- trunc(rr[2] + 86400, "day")
z <- seq(t.start, t.end, by = "day")
z.sub <- seq(t.start, t.end, by = "12 hour")
oceDebug(debug, vectorShow(z))
oceDebug(debug, vectorShow(z.sub))
if (missing(tformat)) {
tformat <- "%b %d"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else if (d <= 60 * 60 * 24 * 31) {
oceDebug(debug, "Time range is between 2 weeks and 1 month (defined as 31 days)\n")
t.start <- trunc(rr[1], "day")
t.end <- trunc(rr[2] + 86400, "day")
z <- seq(t.start, t.end, by = "week")
z.sub <- seq(t.start, t.end, by = "day")
oceDebug(debug, vectorShow(z))
oceDebug(debug, vectorShow(z.sub))
if (missing(tformat)) {
tformat <- "%b %d"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else if (d < 60 * 60 * 24 * 31 * 2) {
oceDebug(debug, "Time range is between 1 and 2 months (defined as 31 days)\n")
t.start <- trunc(rr[1], "days")
t.end <- trunc(rr[2] + 86400, "days")
z <- seq(t.start, t.end, by = "week") # big ticks
z.sub <- seq(t.start, t.end, by = "day") # small ticks
oceDebug(debug, vectorShow(z))
oceDebug(debug, vectorShow(z.sub))
if (missing(tformat)) {
tformat <- "%b %d"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else if (d < 60 * 60 * 24 * 31 * 4) {
oceDebug(debug, "Time range is between 2 and 4 months (defined as 31 days)\n")
t.start <- trunc(rr[1], "days")
t.end <- trunc(rr[2] + 86400, "days")
z <- seq(t.start, t.end, by = "week")
oceDebug(debug, vectorShow(z))
if (missing(tformat)) {
tformat <- "%b %d"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else if (d < 1.1 * 60 * 60 * 24 * 365) {
oceDebug(debug, "Time range is between 4 months and 1 year\n")
rrl <- as.POSIXlt(rr)
rrl[1]$mday <- 1
rrl[2] <- rrl[2] + 31 * 86400
rrl[2]$mday <- 1
t.start <- trunc(rrl[1], "day")
t.end <- trunc(rrl[2] + 86400, "day")
z <- seq(t.start, t.end, by = "month")
oceDebug(debug, vectorShow(z))
if (missing(tformat)) {
tformat <- "%b %d"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else if (d < 3.1 * 60 * 60 * 24 * 365) {
oceDebug(debug, "Time range is between 1 and 3 years\n")
rrl <- as.POSIXlt(rr)
rrl[1]$mday <- 1
rrl[2] <- rrl[2] + 31 * 86400
rrl[2]$mday <- 1
t.start <- trunc(rrl[1], "days")
t.end <- trunc(rrl[2], "days")
z <- seq(t.start, t.end, by = "3 month")
oceDebug(debug, vectorShow(z))
z.sub <- seq(t.start, t.end, by = "month") # small ticks
oceDebug(debug, vectorShow(z.sub))
if (missing(tformat)) {
tformat <- "%Y %b"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
} else {
oceDebug(debug, "Time range is longer than 3 years\n")
class(z) <- c("POSIXt", "POSIXct")
tz <- attr(x, "tzone")
attr(z, "tzone") <- tz
zz <- unclass(as.POSIXlt(z, tz = tz))
zz$mday <- zz$wday <- zz$yday <- 1
zz$isdst <- -1
zz$mon <- zz$hour <- zz$min <- zz$sec <- 0
zz$year <- pretty(zz$year)
M <- length(zz$year)
zz <- lapply(zz, function(x) rep(x, length.out = M))
class(zz) <- c("POSIXt", "POSIXlt")
z <- as.POSIXct(zz, tz = tz)
attr(z, "tzone") <- attr(x, "tzone")
oceDebug(debug, vectorShow(z))
if (missing(tformat)) {
tformat <- "%Y"
oceDebug(debug, "automatic tformat=\"", tformat, "\"\n")
}
oceDebug(debug, vectorShow(z, "z="))
}
if (!mat) {
z <- x[is.finite(x)]
}
# FIXME: I was twiddling the numbers, to get more labels, but xaxs="r" fixes that.
twiddle <- 0.04 * diff(as.numeric(range)) # FIXME: do I need this anymore?
oceDebug(
debug, "range=",
format.POSIXct(rr[1], "%Y-%m-%d %H:%M:%S", tz = "UTC"),
"UTC to ",
format.POSIXct(rr[2], "%Y-%m-%d %H:%M:%S", tz = "UTC"),
"UTC\n"
)
keep <- range[1] <= (z + twiddle) & (z - twiddle) <= range[2]
# oceDebug(debug, vectorShow(keep, "keep"))
oceDebug(debug > 1, vectorShow(z, "z before keep"))
z <- z[keep]
oceDebug(debug > 1, vectorShow(z, "z after keep"))
if (!is.logical(labels)) {
labels <- labels[keep]
} else if (labels[1]) {
if (length(grep("[0-9]+S.*", tformat))) {
# BOOKMARK 1B a special trick to get fractional seconds (cf BOOKMARK 1A)
old <- options("digits.secs")$digits.secs
n <- as.numeric(gsub("^%.*\\.([0-9]*)S.*", "\\1", tformat))
options(digits.secs = n)
labels <- format(z) # "2016-01-01 hh:mm:ss.digits"
labels <- gsub("[0-9]{4}-[0-9]{2}-[0-9]{2} [0-9]{2}:[0-9]{2}:", "", labels)
options(digits.secs = old)
} else {
labels <- format(z, format = tformat)
}
} else if (!labels[1]) {
labels <- rep("", length(z))
}
oceDebug(debug, vectorShow(labels, n = -1))
oceDebug(debug, vectorShow(format(z), n = -1))
oceDebug(debug, vectorShow(z, n = -1))
if (drawTimeRange) {
time.range <- par("usr")[1:2] # axis, not data
class(time.range) <- c("POSIXt", "POSIXct")
attr(time.range, "tzone") <- attr(x, "tzone")[1]
time.range <- as.POSIXlt(time.range)
# time.range.data <- range(x, na.rm=TRUE)
# what was this for?# time.range[1] <- max(time.range[1], time.range.data[1], na.rm=TRUE)
# what was this for?# time.range[2] <- min(time.range[2], time.range.data[2], na.rm=TRUE)
if (!is.null(getOption("oceTimeFormat"))) {
tr1 <- format(time.range[1], getOption("oceTimeFormat"))
tr2 <- format(time.range[2], getOption("oceTimeFormat"))
} else {
tr1 <- format(time.range[1])
tr2 <- format(time.range[2])
}
if (abbreviateTimeRange) {
oceDebug(debug, "abbreviating time format\n")
if (time.range[1]$year == time.range[2]$year) {
tr2 <- substr(tr2, 6, nchar(tr2)) # remove the "YYYY-"
if (time.range[1]$mon == time.range[2]$mon) {
tr2 <- substr(tr2, 4, nchar(tr2)) # remove the "MM-"
if (time.range[1]$mday == time.range[2]$mday) {
tr2 <- substr(tr2, 4, nchar(tr2)) # remove the "DD-"
}
}
}
time.range <- as.POSIXct(time.range)
}
deltat <- mean(diff(as.numeric(x)), na.rm = TRUE)
# only show timezone if hours are shown
oceDebug(debug, "time.range[1]:", format(time.range[1]), "\n")
oceDebug(debug, "round(time.range[1], 'days'):", format(round(time.range[1], "days")), "\n")
oceDebug(debug, "time.range[2]:", format(time.range[2]), "\n")
oceDebug(debug, "round(time.range[2], 'days'):", format(round(time.range[2], "days")), "\n")
tzone <- c(attr(time.range[1], "tzone"), attr(time.range[2], "tzone"))
# Only show tzone if it is UTC (https://github.com/dankelley/oce/issues/1811)
tzone <- if (all(tzone == "UTC")) c("", "UTC") else rep("", 2)
# The below is not fool-proof, depending on how xlim might have been supplied; see
# https://bugs.r-project.org/bugzilla3/show_bug.cgi?id=14449
if (diff(as.numeric(time.range)) > 7 * 86400) {
label <- paste(tr1, tr2, sep = " to ")
} else {
label <- paste(tr1, tzone[1], " to ", tr2, tzone[2], sep = "")
}
if (drawFrequency && is.finite(1 / deltat)) {
label <- paste(label, "@", sprintf("%.4g Hz", 1 / deltat), sep = " ")
}
oceDebug(debug, "label=\"", label, "\" at cex.lab=", cex.lab, "\n", sep = "")
mtext(label, side = if (side == 1) 3 else 1, cex = cex.lab * par("cex"), adj = 0)
oceDebug(debug, "cex.axis=", cex.axis, "; par(\"cex\")=", par("cex"), "\n")
}
if (nchar(main) > 0) {
mtext(main, side = if (side == 1) 3 else 1, cex = cex.lab * par("cex"), adj = 1)
}
oceDebug(debug, vectorShow(z, "z="))
if (length(z.sub) > 0) {
axis(side, at = z.sub, line = 0, labels = FALSE, tcl = -0.25)
oceDebug(debug, vectorShow(z.sub, "z.sub="))
}
oceDebug(debug, vectorShow(labels, "labels="))
# ocex <- par("cex")
ocex.axis <- par("cex.axis")
ocex.lab <- par("cex.lab")
ocex.main <- par("cex.main")
omgp <- par("mgp")
par(cex.axis = cex.axis, cex.lab = cex.lab, cex.main = cex.main, mgp = mgp, tcl = -0.5)
# axis(side, at=z, line=0, labels=labels, cex=cex, cex.axis=cex.axis, cex.main=cex.main, mar=mar, mgp=mgp)
# If the user did gave tformat, shorten the strings for aesthetic reasons.
if (!tformatGiven) {
oceDebug(debug, "axis labels before shortenTimeString(): '", paste(labels, "', '"), "'\n")
labels <- shortenTimeString(labels, debug = debug - 1)
oceDebug(debug, "axis labels after shortenTimeString(): '", paste(labels, "', '"), "'\n")
}
axis(side, at = z, line = 0, labels = labels, mgp = mgp, cex.axis = cex.axis, cex.lab = cex.lab, cex.main = cex.main, ...)
par(cex.axis = ocex.axis, cex.lab = ocex.lab, cex.main = ocex.main, mgp = omgp)
oceDebug(debug, "} # oce.axis.POSIXct()\n", unindent = 1, style = "bold")
zzz <- as.numeric(z)
if (1 < length(zzz)) {
xaxp <- c(min(zzz, na.rm = TRUE), max(zzz, na.rm = TRUE), -1 + length(zzz))
par(xaxp = xaxp)
}
invisible(z) # FIXME: or z.sub?
} # oce.axis.POSIXct()
#' Convert a Numeric Time to Hour, Minute, and Second
#'
#' @param t a vector of factors or character strings, in the format 1200 for
#' 12:00, 0900 for 09:00, etc.
#'
#' @param default value to be used for the returned hour, minute and second if
#' there is something wrong with the input value (e.g. its length exceeds 4
#' characters, or it contains non-numeric characters)
#'
#' @return A list containing `hour`, `minute`, and `second`, the
#' last of which is always zero.
#'
#' @author Dan Kelley
#'
#' @examples
#' t <- c("0900", "1234")
#' numberAsHMS(t)
#'
#' @family things related to time
numberAsHMS <- function(t, default = 0) {
if ("factor" == class(t)) {
t <- as.character(t)
}
if (!is.character(t)) {
stop("can only handle strings or factors")
}
nc <- nchar(t)
t[nc == 0] <- "0000"
nc <- nchar(t)
t[nc == 1] <- paste("000", t[nc == 1], sep = "")
t[nc == 2] <- paste("00", t[nc == 2], sep = "")
t[nc == 3] <- paste("0", t[nc == 3], sep = "")
nc <- nchar(t)
warn <- options("warn")$warn
options(warn = -1)
try(
{
hour <- as.numeric(substr(t, 1, 2))
minute <- as.numeric(substr(t, 3, 4))
second <- rep(0, length(hour))
},
silent = TRUE
)
options(warn = warn)
bad <- nc > 4 | grepl("[^[0-9]]*", t)
hour[bad] <- default
minute[bad] <- default
second[bad] <- default
list(hour = hour, minute = minute, second = second)
}
#' Convert a Numeric Time to a POSIXct Time
#'
#' This converts numerical values into POSIXct times. There are many
#' schemes for doing this, with the `type` parameter being used
#' to select between them. See \dQuote{Details} for a listing, broken
#' down by scheme.
#'
#' The possible choices for `type` are as listed below.
#'
#' * `"unix"` handles Unix times, measured in seconds since the start
#' of the year 1970.
#'
#' * `"matlab"` handles Matlab times, measured in days since what
#' MathWorks (reference 1) calls ``January 0, 0000'' (i.e. \code{ISOdatetime(0, 1, 1, 0,
#' 0, 0)} in R notation).
#'
#' * `"gps"` handles the Global Positioning System convention. The
#' scheme is complicated, owing to hardware limitations of
#' GPS satellites. As illustrated in Example 3, `t` may be
#' a matrix with either 2 or 3 columns. In the 2-column format, the
#' first column holds the number of weeks after 1999-08-22, modulo
#' 1024 (approximately 19.6 years), and the second column (here
#' and also in the 3-column format) holds the number of seconds
#' in the referenced week, with leap seconds being handled with
#' the `leap` parameter. The modulo calculation is required
#' because GPS satellites dedicate only 10 bits to the
#' week number. The resultant ambiguity (e.g. a rollover in
#' 2019-04-07) is addressed in the 3-column format, in which the
#' last column holds the number of 1024-week rollover events
#' since 1980-01-06. Users should set this column to 0 for times
#' prior to 1999-08-22, to 1 for later times prior
#' to 2019-04-07, to 2 for later times prior to 2038-11-21,
#' etc. However, there will be an exception to this rule,
#' when satellites start dedicating 12 bits to the week value.
#' For such data, the third column will need to be 0 for all times
#' prior to 2137-01-06.
#'
#' * `"argo"` handles Argo times, measured in days since the start of
#' the year 1900.
#'
#' * `"excel"` handles Excel times, measured in days since the start of
#' the year 1900. (Note that excel incorrectly regards 1900 as a leap year,
#' so 1 day is subtracted from `t` unless the time is less than or equal
#' to 1900 Feb 28. Note that NA is returned for the day 60, which is
#' what excel codes for "Feb 29, 1900", the non-existing day that excel
#' accepts.
#'
#' * `"ncep1"` handles NCEP times, measured in hours since the start
#' of the year 1800.
#'
#' * `"ncep2"` handles NCEP times, measured in days since the start of
#' the year 1. (Note that, for reasons that are unknown at this time, a simple
#' R expression of this definition is out by two days compared with the UDUNITS
#' library, which is used by NCEP. Therefore, a two-day offset is applied. See
#' references 2 and 3.)
#'
#' * `"sas"` handles SAS times, indicated by `type="sas"`, have
#' origin at the start of 1960.
#'
#' * `"spss"` handles SPSS times, in seconds after 1582-10-14.
#'
#' * `"yearday"` handles a convention in which `t` is a
#' two-column matrix, with the first column being the year, and the second the
#' yearday (starting at 1 for the first second of January 1, to match the
#' convention used by Sea-Bird CTD software).
#'
#' * `"epic"` handles a convention used in the EPIC software library,
#' from the Pacific Marine Environmental Laboratory, in which `t` is a
#' two-column matrix, with the first column being the julian Day (as defined in
#' [julianDay()], for example), and with the second column being the
#' millisecond within that day. See reference 4.
#'
#' * `"vms"` handles a convention used in the VMS operating system and
#' for Modified Julian Day, in which `t` is the number of seconds
#' past 1859-11-17T00:00:00 UTC. See reference 5.
#'
#' @param t an integer corresponding to a time, in a way that depends on
#' `type`.
#'
#' @param type character value indicating the time type. The permitted values
#' are `"argo"`, `"epic"`, `"excel"`, `"gps"`, `"matlab"`, `"ncep1"`, `"ncep2"`,
#' `"sas"`, `"spss"`, `"unix"`, and `"yearday"`, the first of these being
#' the default.
#'
#' @param tz a string indicating the time zone, by default `"UTC"`.
#'
#' @param leap a logical value, TRUE by default, that applies only
#' if `type` is `"gps"`. If `leap` is TRUE, then the built-in
#' dataset named [.leap.seconds] is consulted to find
#' of the number of leap seconds between 1980 (when the
#' GPS program started) and the time computed from the other
#' parameters, and the return value is decreased accordingly
#' (see Example 3).
#'
#' @return A [POSIXct()] time vector.
#'
#' @author Dan Kelley
#'
#' @family things related to time
#'
#' @references
#' 1. Matlab times:
#' `https://www.mathworks.com/help/matlab/ref/datenum.html`
#'
#' 2. NCEP times: `https://psl.noaa.gov/data/gridded/faq.html`
#'
#' 3. Problem with NCEP times:
#' `https://github.com/dankelley/oce/issues/738`
#'
#' 4. EPIC times: software and manuals at `https://www.pmel.noaa.gov/epic/download/index.html#epslib`;
#' see also Denbo, Donald W., and Nancy N. Soreide. \dQuote{EPIC.} Oceanography 9 (1996).
#' \doi{10.5670/oceanog.1996.10}
#'
#' 5. VMS times: https://en.wikipedia.org/wiki/Epoch_(computing)
#'
#' 6. GPS times: https://www.labsat.co.uk/index.php/en/gps-time-calculator
#'
#'
#' @examples
#' # Example 1. default (unix)
#' numberAsPOSIXct(0)
#'
#' # Example 2. Matlab
#' numberAsPOSIXct(1, type = "matlab")
#'
#' # Example 3. GPS with default week rollover or with no rollover (Canada Day, year 2010)
#' numberAsPOSIXct(cbind(566, 345615), type = "gps")
#' numberAsPOSIXct(cbind(566, 345615, 1), type = "gps")
#' numberAsPOSIXct(cbind(1024 + 566, 345615, 0), type = "gps")
#' # Show how to deal with leap seconds (15 of them, in this case)
#' sum(as.POSIXct("1980-01-01") < .leap.seconds & .leap.seconds <= as.POSIXct("2010-07-01"))
#' -15 + numberAsPOSIXct(cbind(1024 + 566, 345615, 0), type = "gps", leap = FALSE)
#'
#' # Example 4. yearday
#' numberAsPOSIXct(cbind(2013, 1), type = "yearday") # start of 2013
#'
#' # Example 5. Epic time, one hour into Canada Day of year 2018. In computing the
#' # Julian day, note that this starts at noon.
#' jd <- julianDay(as.POSIXct("2018-07-01 12:00:00", tz = "UTC"))
#' numberAsPOSIXct(cbind(jd, 1e3 * 1 * 3600), type = "epic", tz = "UTC")
#'
#' # Example 6. Julian day, note that this starts at noon.
#' jd <- julianDay(as.POSIXct("2018-07-01 12:00:00", tz = "UTC"))
#' numberAsPOSIXct(cbind(jd, 1e3 * 1 * 3600), type = "epic", tz = "UTC")
#'
#' @family things related to time
numberAsPOSIXct <- function(t, type = "unix", tz = "UTC", leap = TRUE) {
typeAllowed <- c(
"argo", "epic", "excel", "gps", "matlab",
"ncep1", "ncep2", "sas", "spss", "unix", "yearday"
)
if (!type %in% typeAllowed) {
stop("type must be one of the following: \"", paste(typeAllowed, collapse = "\", \""), "\".", sep = "")
}
if (type == "unix") {
# We add something with a timezone, and then subtract it, as a trick to inherit the timezone
tref <- if (!is.null(tz)) as.POSIXct("2000-01-01", tz = tz) else as.POSIXct("2000-01-01")
return(tref + as.numeric(t) - as.numeric(tref))
} else if (type == "matlab") {
# R won't take a day "0", so subtract one
return(as.POSIXct(ISOdatetime(0000, 01, 01, 0, 0, 0, tz = tz) + 86400 * (t - 1)))
} else if (type == "yearday") {
if (2 != ncol(t)) {
stop("'t' must have two columns, one for year, the other for yearday")
}
return(ISOdatetime(t[, 1], 1, 1, 0, 0, 0, tz = tz) + (t[, 2] - 1) * 24 * 3600)
} else if (type == "argo") {
return(t * 86400 + as.POSIXct("1900-01-01 00:00:00", tz = "UTC"))
} else if (type == "excel") {
# We need a one-day offset if time is after Feb 28, 1900,
# because Excel thinks 1900 was a leap year. We can check for
# this by day count. See https://github.com/dankelley/oce/issues/1591
# for a discussion. Note that we return NA for day 60, because
# that is what excel produces for "Feb 29, 1900", which is in
# fact a non-existent day so I think we ought to inform oce users
# of that fact, with a NA.
offset <- ifelse(t < 61, 0, 1) # excel thinks 1900 is a leap year
rval <- 86400 * (t - 1 - offset) + as.POSIXct("1900-01-01 00:00:00", tz = "UTC")
rval[t == 60] <- NA
return(rval)
} else if (type == "ncep1") {
# hours since the start of 1800
return(t * 3600 + as.POSIXct("1800-01-01 00:00:00", tz = "UTC"))
} else if (type == "ncep2") {
# days since 1-1-1 00:00:0.0 (supposedly, but offset to match a test case; see
resOriginal <- t * 86400 + as.POSIXct("0001-01-01 00:00:00", tz = "UTC")
return(resOriginal - 2 * 86400) # kludge for ht of https://github.com/dankelley/oce/issues/738
} else if (type == "gps") {
if (!is.matrix(t) || !(ncol(t) %in% 2:3)) {
stop("for GPS times, 't' must be a matrix with either 2 or 3 columns")
}
if (ncol(t) == 2L) {
t <- cbind(t, 1)
}
# Account for leap seconds since the GPS start time in 1980 (for the present week wraparound grouping).
# 20171014 See http://en.wikipedia.org/wiki/Leap_second and other sources for a list. Updates can happen
# 20171014 # on June 30 and December 31 of any given year. The information below was last updated
# 20171014 # in January, 2017.
# 20171014 # leapsOLD <- as.POSIXct(strptime(c("1981-07-01", "1982-07-01", "1983-07-01", "1985-07-01", "1987-12-31",
# 20171014 # "1989-12-31", "1990-12-31", "1992-07-01", "1993-07-01", "1994-07-01",
# 20171014 # "1995-12-31", "1997-07-01", "1998-12-31", "2005-12-31", "2008-12-31",
# 20171014 # "2012-07-01", "2015-07-01", "2016-12-31"),
# 20171014 # format="%Y-%m-%d", tz="UTC"))
# 20171014 message("leapsOLD ", paste(leapsOLD, collapse=" "))
# 20171014 message("leaps A ", paste(leaps, collapse=" "))
# 20171014 message("leaps B ", paste(leaps, collapse=" "))
# 20171014 message("leaps C ", paste(leaps, collapse=" "))
week <- t[, 1] + 1024 * (t[, 3] - 1)
second <- t[, 2]
t <- as.POSIXct("1999-08-22 00:00:00", tz = "UTC") + 86400 * 7 * week + second
if (leap) {
leaps <- as.POSIXlt(.leap.seconds, tz = "UTC")
leapOffset <- sum(as.POSIXlt("1980-01-01 00:00:00", tz = "UTC") < leaps & leaps < as.POSIXlt(t, tz = "UTC"))
t <- t - leapOffset # subtrace 1 second for each leap-second year between 1980 and the computed time
}
} else if (type == "spss") {
t <- as.POSIXct(t, origin = "1582-10-14", tz = tz)
} else if (type == "sas") {
t <- as.POSIXct(t, origin = "1960-01-01", tz = tz)
} else if (type == "epic") {
if (!is.matrix(t) || dim(t)[2] != 2) {
stop("for epic times, 't' must be a two-column matrix, with first column the julian day, and second the millisecond within that day")
}
r <- do_epic_time_to_ymdhms(t[, 1], t[, 2])
t <- ISOdatetime(r$year, r$month, r$day, r$hour, r$minute, r$second, tz = tz)
} else if (type == "vms") {
t <- as.POSIXct(t, origin = "1858-11-17", tz = tz)
} else {
stop("unknown type '", type, "'")
}
t
}
#' Plot an Inset Diagram
#'
#' Adds an inset diagram to an existing plot. Note that if the inset is a map
#' or coastline, it will be necessary to supply `inset=TRUE` to prevent
#' the inset diagram from occupying the whole device width. After
#' `plotInset()` has been called, any further plotting will take place
#' within the inset, so it is essential to finish a plot before drawing an
#' inset.
#'
#' @param xleft location of left-hand of the inset diagram, in the existing
#' plot units. (PROVISIONAL FEATURE: this may also be `"bottomleft"`, to
#' put the inset there. Eventually, other positions may be added.)
#'
#' @param ybottom location of bottom side of the inset diagram, in the existing
#' plot units.
#'
#' @param xright location of right-hand side of the inset diagram, in the
#' existing plot units.
#'
#' @param ytop location of top side of the inset diagram, in the existing plot
#' units.
#'
#' @param expr An expression that draws the inset plot. This may be a single
#' plot command, or a sequence of commands enclosed in curly braces.
#'
#' @param mar margins, in line heights, to be used at the four sides of the
#' inset diagram. (This is often helpful to save space.)
#'
#' @param debug a flag that turns on debugging. Set to 1 to get a moderate
#' amount of debugging information, or to 2 to get more.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' # power law in linear and log form
#' x <- 1:10
#' y <- x^2
#' plot(x, y, log = "xy", type = "l")
#' plotInset(3, 1, 10, 8,
#' expr = plot(x, y, type = "l", cex.axis = 3 / 4, mgp = c(3 / 2, 1 / 2, 0)),
#' mar = c(2.5, 2.5, 1, 1)
#' )
#'
#' # CTD data with location
#' data(ctd)
#' plot(ctd, which = "TS")
#' plotInset(29.9, 2.7, 31, 10,
#' expr = plot(ctd,
#' which = "map",
#' coastline = "coastlineWorld",
#' span = 5000, mar = NULL, cex.axis = 3 / 4
#' )
#' )
plotInset <- function(
xleft, ybottom, xright, ytop, expr,
mar = c(2, 2, 1, 1),
debug = getOption("oceDebug")) {
opar <- par(no.readonly = TRUE)
mai <- par("mai") # bottom left top right
# omfg <- par("mfg") # original mfg
xLog <- par("xlog")
yLog <- par("ylog")
usr <- par("usr") # xmin xmax ymin ymax
fin <- par("fin") # figure width height
x2in <- function(x) {
if (xLog) {
mai[2] + (log10(x) - usr[1]) * (fin[1] - mai[2] - mai[4]) / (usr[2] - usr[1])
} else {
mai[2] + (x - usr[1]) * (fin[1] - mai[2] - mai[4]) / (usr[2] - usr[1])
}
}
y2in <- function(y) {
if (yLog) {
mai[1] + (log10(y) - usr[3]) * (fin[2] - mai[1] - mai[3]) / (usr[4] - usr[3])
} else {
mai[1] + (y - usr[3]) * (fin[2] - mai[1] - mai[3]) / (usr[4] - usr[3])
}
}
if (is.character(xleft)) {
if (xleft != "bottomleft") {
stop("only named position is \"bottomleft\"")
}
f1 <- 0.02
f2 <- 1 / 3
if (xLog) {
stop("cannot handle xlog yet")
}
xleft <- usr[1] + f1 * (usr[2] - usr[1])
xright <- usr[1] + f2 * (usr[2] - usr[1])
if (yLog) {
stop("cannot handle ylog yet")
}
ybottom <- usr[3] + f1 * (usr[4] - usr[3])
ytop <- usr[3] + f2 * (usr[4] - usr[3])
} else {
oceDebug(debug, "plotInset(xleft=", xleft, ", ybottom=", ybottom,
", xright=", xright, ", ytop=", ytop,
", ...) {\n",
sep = "", unindent = 1
)
}
oceDebug(debug, "par(\"mfg\")=", par("mfg"), "\n")
oceDebug(debug, "par(\"mai\")=", par("mai"), "\n")
oceDebug(debug, "par(\"usr\")=", par("usr"), "\n")
oceDebug(debug, "par(\"fin\")=", fin, "(figure width and height)\n")
nmai <- c(y2in(ybottom), x2in(xleft), fin[2] - y2in(ytop), fin[1] - x2in(xright))
oceDebug(debug, "nmai:", nmai, "\n")
if (any(nmai < 0)) {
warning("part of the inset is off the page")
}
nmai[nmai < 0] <- 0
if (nmai[1] < 0) nmai[1] <- fin[1]
if (nmai[2] < 0) nmai[2] <- fin[1]
if (nmai[3] > fin[2] - 0.2) nmai[3] <- fin[2] - 0.2
if (nmai[4] > fin[1] - 0.2) nmai[4] <- fin[1] - 0.2
oceDebug(debug, "nmai:", nmai, "(after trimming negatives)\n")
# mfg2 <- par("mfg")
par(new = TRUE, mai = nmai)
thismar <- par("mar")
par(mar = thismar + mar)
if (debug > 1) {
cat("\n\nBEFORE expr, PAR IS:\n")
str(par())
}
mfg <- par("mfg")
oceDebug(debug, "BEFORE expr, mfg=", mfg, "\n")
# Draw the inset plot (or perform any action, actually)
expr
if (mfg[1] == mfg[3] && mfg[2] == mfg[4]) {
oceDebug(debug, "setting new=FALSE; mfg=", mfg, "\n")
par(new = FALSE)
} else {
oceDebug(debug, "setting new=TRUE; mfg=", mfg, "\n")
par(new = TRUE)
}
# Reset some things that could have been set in the inset, and
# then adjust new appropriately.
par(usr = opar$usr, mai = opar$mai, cex = opar$cex, lwd = opar$lwd, lty = opar$lty, bg = opar$bg)
oceDebug(debug, "} # plotInset()\n", unindent = 1)
invisible(NULL)
}
#' Oce Version of as.POSIXct
#'
#' Each format in `timeFormats` is used in turn as the `format`
#' argument to [as.POSIXct()], and the first that produces a
#' non-`NA` result is used. If `timeFormats` is missing, the
#' following formats are tried, in the stated order:
#'
#' * `"\%b \%d \%Y \%H:\%M:\%S"` (e.g. `"Jul 1 2013 01:02:03"`)
#'
#' * `"\%b \%d \%Y"` (e.g. `"Jul 1 2013"`)
#'
#' * `"\%B \%d \%Y \%H:\%M:\%S"` (e.g. `"July 1 2013 01:02:03"`)
#'
#' * `"\%B \%d \%Y"` (e.g. `"July 1 2013"`)
#'
#' * `"\%d \%b \%Y \%H:\%M:\%S"` (e.g. `"1 Jul 2013 01:02:03"`)
#'
#' * `"\%d \%b \%Y"` (e.g. `"1 Jul 2013"`)
#'
#' * `"\%d \%B \%Y \%H:\%M:\%S"` (e.g. `"1 July 2013 01:02:03"`)
#'
#' * `"\%d \%B \%Y"` (e.g. `"1 July 2013"`)
#'
#' * `"\%Y-\%m-\%d \%H:\%M:\%S"` (e.g. `"2013-07-01 01:02:03"`)
#'
#' * `"\%Y-\%m-\%d"` (e.g. `"2013-07-01"`)
#'
#' * `"\%Y-\%b-\%d \%H:\%M:\%S"` (e.g. `"2013-July-01 01:02:03"`)
#'
#' * `"\%Y-\%b-\%d"` (e.g. `"2013-Jul-01"`)
#'
#' * `"\%Y-\%B-\%d \%H:\%M:\%S"` (e.g. `"2013-July-01 01:02:03"`)
#'
#' * `"\%Y-\%B-\%d"` (e.g. `"2013-July-01"`)
#'
#' * `"\%d-\%b-\%Y \%H:\%M:\%S"` (e.g. `"01-Jul-2013 01:02:03"`)
#'
#' * `"\%d-\%b-\%Y"` (e.g. `"01-Jul-2013"`)
#'
#' * `"\%d-\%B-\%Y \%H:\%M:\%S"` (e.g. `"01-July-2013 01:02:03"`)
#'
#' * `"\%d-\%B-\%Y"` (e.g. `"01-July-2013"`)
#'
#' * `"\%Y/\%b/\%d \%H:\%M:\%S"` (e.g. `"2013/Jul/01 01:02:03"`)
#'
#' * `"\%Y/\%b/\%d"` (e.g. `"2013/Jul/01"`)
#'
#' * `"\%Y/\%B/\%d \%H:\%M:\%S"` (e.g. `"2013/July/01 01:02:03"`)
#'
#' * `"\%Y/\%B/\%d"` (e.g. `"2013/July/01"`)
#'
#' * `"\%Y/\%m/\%d \%H:\%M:\%S"` (e.g. `"2013/07/01 01:02:03"`)
#'
#' * `"\%Y/\%m/\%d"` (e.g. `"2013/07/01"`)
#'
#' @param time Character string with an indication of the time.
#'
#' @param timeFormats Optional vector of time formats to use, as for [as.POSIXct()].
#'
#' @param tz Time zone.
#'
#' @return A time as returned by [as.POSIXct()].
#'
#' @author Dan Kelley
#'
#' @examples
#' decodeTime("July 1 2013 01:02:03")
#' decodeTime("Jul 1 2013 01:02:03")
#' decodeTime("1 July 2013 01:02:03")
#' decodeTime("1 Jul 2013 01:02:03")
#' decodeTime("2013-07-01 01:02:03")
#' decodeTime("2013/07/01 01:02:03")
#' decodeTime("2013/07/01")
#'
#' @family functions relating to time
decodeTime <- function(time, timeFormats, tz = "UTC") {
if (missing(timeFormats)) {
timeFormats <- c(
"%b %d %Y %H:%M:%S", "%b %d %Y", # Jul 1 2013
"%B %d %Y %H:%M:%S", "%B %d %Y", # July 1 2013
"%d %b %Y %H:%M:%S", "%d %b %Y", # 1 Jul 2013
"%d %B %Y %H:%M:%S", "%d %B %Y", # 1 July 2013
"%Y-%m-%d %H:%M:%S", "%Y-%m-%d", # 2013-07-01
"%Y-%b-%d %H:%M:%S", "%Y-%b-%d", # 2013-Jul-01
"%Y-%B-%d %H:%M:%S", "%Y-%B-%d", # 2013-July-01
"%d-%b-%Y %H:%M:%S", "%d-%b-%Y", # 01-Jul-2013
"%d-%B-%Y %H:%M:%S", "%d-%B-%Y", # 01-July-2013
"%Y/%m/%d %H:%M:%S", "%Y/%m/%d", # 2013/07/01
"%Y/%b/%d %H:%M:%S", "%Y/%b/%d", # 2013/Jul/01
"%Y/%B/%d %H:%M:%S", "%Y/%B/%d", # 2013/July/01
"%Y/%m/%d %H:%M:%S", "%Y/%m/%d"
) # 2013/07/01
}
res <- NA
for (format in timeFormats) {
res <- as.POSIXct(time, format = format, tz = tz)
if (!is.na(res)) {
break
}
}
res
}
#' Draw a Direction Field
#'
#' The direction field is indicated variously, depending on the value of
#' `type`:
#'
#' * For `type=1`, each indicator is drawn with a symbol, according to the
#' value of `pch` (either supplied globally, or as an element of the
#' `...` list) and of size `cex`, and color `col`. Then, a
#' line segment is drawn for each, and for this `lwd` and `col` may
#' be set globally or in the `...` list.
#'
#' * For `type=2`, the points are not drawn, but arrows are drawn instead
#' of the line segments. Again, `lwd` and `col` control the type of
#' the line.
#'
#' @param x,y coordinates at which velocities are specified. The
#' length of `x` and `y` depends on the form of `u`
#' and `v` (vectors or matrices).
#'
#' @param u,v velocity components in the x and y directions. Can be
#' either vectors with the same length as `x, y`, or
#' matrices, of dimension `length(x)` by `length(y)`.
#'
#' @param scalex,scaley scale to be used for the velocity arrows.
#' Exactly one of these must be specified. Arrows that have
#' `u^2+v^2=1` will have length `scalex` along the x
#' axis, or `scaley` along the y axis, according to which
#' argument is given.
#'
#' @param skip either an integer, or a two-element vector indicating
#' the number of points to skip when plotting arrows (for the
#' matrix `u, v` case). If a single value, the same
#' `skip` is applied to both the `x` and `y`
#' directions. If a two-element vector, specifies different values
#' for the `x` and `y` directions.
#'
#' @param length indication of *width* of arrowheads. The
#' somewhat confusing name of this argument is a consequence of
#' the fact that it is passed to [arrows()] for drawing
#' arrows. Note that the present default is smaller than the
#' default used by [arrows()].
#'
#' @param add if `TRUE`, the arrows are added to an existing
#' plot; otherwise, a new plot is started by calling
#' [plot()] with `x`, `y` and `type="n"`.
#' In other words, the plot will be very basic. In most cases, the
#' user will probably want to draw a diagram first, and `add`
#' the direction field later.
#'
#' @param type indication of the style of arrow-like indication of the
#' direction.
#'
#' @param col color of line segments or arrows; see [par()] for meaning
#'
#' @param pch,cex plot character and expansion factor, used for
#' `type=1`; see [par()] for meanings
#'
#' @param lwd,lty line width and type, used for `type=2`; see [par()] for meaning
#'
#' @param xlab,ylab `x` and `y` axis labels
#'
#' @param debug debugging value; set to a positive integer to get
#' debugging information.
#'
#' @param ... other arguments to be passed to plotting functions
#' (e.g. axis labels, etc).
#'
#' @return None.
#'
#' @examples
#' library(oce)
#' plot(c(-1.5, 1.5), c(-1.5, 1.5), xlab = "", ylab = "", type = "n")
#' drawDirectionField(
#' x = rep(0, 2), y = rep(0, 2),
#' u = c(1, 1), v = c(1, -1), scalex = 0.5, add = TRUE
#' )
#' plot(c(-1.5, 1.5), c(-1.5, 1.5), xlab = "", ylab = "", type = "n")
#' drawDirectionField(
#' x = rep(0, 2), y = rep(0, 2),
#' u = c(1, 1), v = c(1, -1), scalex = 0.5, add = TRUE, type = 2
#' )
#'
#' # 2D example
#' x <- seq(-2, 2, 0.1)
#' y <- x
#' xx <- expand.grid(x, y)[, 1]
#' yy <- expand.grid(x, y)[, 2]
#' z <- matrix(xx * exp(-xx^2 - yy^2), nrow = length(x))
#' gz <- grad(z, x, y)
#' drawDirectionField(x, y, gz$gx, gz$gy, scalex = 0.5, type = 2, len = 0.02)
#' oceContour(x, y, z, add = TRUE)
#'
#' @author Dan Kelley and Clark Richards
drawDirectionField <- function(
x, y, u, v, scalex, scaley, skip, length = 0.05, add = FALSE,
type = 1, col = par("fg"), pch = 1, cex = par("cex"),
lwd = par("lwd"), lty = par("lty"),
xlab = "", ylab = "",
debug = getOption("oceDebug"),
...) {
oceDebug(debug, "drawDirectionField(...) {\n", unindent = 1)
if (missing(x) || missing(y) || missing(u) || missing(v)) {
stop("must supply x, y, u, and v")
}
if ((missing(scalex) && missing(scaley)) || (!missing(scalex) && !missing(scaley))) {
stop("either 'scalex' or 'scaley' must be specified (but not both)")
}
if (is.matrix(u) && !is.matrix(v)) {
stop("if 'u' is a matrix, then 'v' must also be a matrix")
}
if (is.matrix(v) && !is.matrix(u)) {
stop("if 'v' is a matrix, then 'u' must also be a matrix")
}
if (is.matrix(u) && is.matrix(v)) {
oceDebug(debug, "u, v are matrices")
if ((dim(u)[1] != dim(v)[1]) && (dim(u)[2] != dim(v)[2])) {
stop("dimensions of u and v must match")
}
dim <- dim(u)
if (length(x) != dim[1] || length(y) != dim[2]) {
stop("lengths of x, y must match dimensions of u, v")
}
if (missing(skip)) {
skip <- c(1, 1)
} else {
if (length(skip) == 1) {
skip <- rep(skip, 2)
}
}
nx <- length(x)
ny <- length(y)
ix <- seq(1, nx, skip[1])
iy <- seq(1, ny, skip[2])
xx <- expand.grid(x[ix], y[iy])[, 1]
yy <- expand.grid(x[ix], y[iy])[, 2]
uu <- as.vector(u[ix, iy])
vv <- as.vector(v[ix, iy])
} else {
if (length(x) != length(y)) {
stop("lengths of x and y must match")
}
if (length(x) != length(u)) {
stop("lengths of x and u must match")
}
if (length(x) != length(v)) {
stop("lengths of x and v must match")
}
xx <- x
yy <- y
uu <- u
vv <- v
}
if (!add) {
plot(xx, yy, type = "n", xlab = xlab, ylab = ylab, ...)
}
usr <- par("usr")
pin <- par("pin")
# mai <- par("mai")
xPerInch <- diff(usr[1:2]) / pin[1]
yPerInch <- diff(usr[3:4]) / pin[2]
oceDebug(debug, "pin=", pin, "usr=", usr, "xPerInch=", xPerInch, "yPerInch=", yPerInch, "\n")
if (missing(scaley)) {
oceDebug(debug, "scaling for x\n")
uPerX <- 1 / scalex
vPerY <- uPerX * xPerInch / yPerInch
} else {
vPerY <- 1 / scaley
uPerX <- vPerY * yPerInch / xPerInch
oceDebug(debug, "scaling for y\n")
}
oceDebug(debug, "uPerX=", uPerX, "\n")
oceDebug(debug, "vPerY=", vPerY, "\n")
len <- sqrt((u / uPerX)^2 + (v / vPerY)^2)
ok <- len > 0.0
xx <- xx[ok]
yy <- yy[ok]
uu <- uu[ok]
vv <- vv[ok]
if (type == 1) {
points(xx, yy, col = col, pch = pch, cex = cex)
segments(x0 = xx, y0 = yy, x1 = xx + uu / uPerX, y1 = yy + vv / vPerY, col = col, lwd = lwd, lty = lty)
} else if (type == 2) {
arrows(x0 = xx, y0 = yy, x1 = xx + uu / uPerX, y1 = yy + vv / vPerY, length = length, col = col, lwd = lwd, lty = lty)
} else {
stop("unknown value of type ", type)
}
oceDebug(debug, "} # drawDirectionField\n", unindent = 1)
}
#' Oce Variant of contour
#'
#' This provides something analogous to [contour()], but with the
#' ability to flip x and y.
#' Setting `revy=TRUE` can be helpful if the `y` data represent
#' pressure or depth below the surface.
#'
#' @aliases oce.contour oceContour
#'
#' @param x values for x grid.
#'
#' @param y values for y grid.
#'
#' @param z matrix for values to be contoured. The first dimension of `z`
#' must equal the number of items in `x`, etc.
#'
#' @param revx set to `TRUE` to reverse the order in which the labels on
#' the x axis are drawn
#'
#' @param revy set to `TRUE` to reverse the order in which the labels on
#' the y axis are drawn
#'
#' @param add logical value indicating whether the contours should be added to
#' a pre-existing plot.
#'
#' @param tformat time format; if not supplied, a reasonable choice will be
#' made by [oce.axis.POSIXct()], which draws time axes.
#'
#' @param drawTimeRange logical, only used if the `x` axis is a time. If
#' `TRUE`, then an indication of the time range of the data (not the axis)
#' is indicated at the top-left margin of the graph. This is useful because
#' the labels on time axes only indicate hours if the range is less than a day,
#' etc.
#'
#' @param debug a flag that turns on debugging; set to 1 to information about
#' the processing.
#'
#' @param \dots optional arguments passed to plotting functions.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' data(topoWorld)
#' # coastline now, and in last glacial maximum
#' lon <- topoWorld[["longitude"]]
#' lat <- topoWorld[["latitude"]]
#' z <- topoWorld[["z"]]
#' oce.contour(lon, lat, z, levels = 0, drawlabels = FALSE)
#' oce.contour(lon, lat, z, levels = -130, drawlabels = FALSE, col = "blue", add = TRUE)
oce.contour <- function(
x, y, z, revx = FALSE, revy = FALSE, add = FALSE,
tformat, drawTimeRange = getOption("oceDrawTimeRange"),
debug = getOption("oceDebug"), ...) {
# dots <- list(...)
# dotsNames <- names(dots)
mustReverseX <- any(0 > diff(order(x)))
mustReverseY <- any(0 > diff(order(y)))
oceDebug(debug, "mustReverseX:", mustReverseX, "\n")
oceDebug(debug, "mustReverseY:", mustReverseY, "\n")
# perhaps get (x,y,z) from x, etc., trying to emulate contour()
if (missing(z)) {
if (!missing(x)) {
if (is.list(x)) {
z <- x$z
y <- x$y
x <- x$x
} else {
z <- x
x <- seq.int(0, 1, length.out = nrow(z))
}
} else {
stop("no 'z' matrix specified")
}
} else if (is.list(x)) {
y <- x$y
x <- x$x
}
zdim <- dim(z)
# store local values for the tricky cases of reversing axes etc
xx <- x
yy <- y
zz <- z
if (mustReverseX) {
xx <- rev(xx)
zz <- zz[seq.int(zdim[1], 1), ]
}
if (mustReverseY) {
yy <- rev(yy)
zz <- zz[, seq.int(zdim[2], 1)]
}
if (add) {
contour(xx, yy, zz, add = TRUE, ...)
} else {
if (revx) {
xx <- rev(xx)
}
if (revy) {
yy <- rev(yy)
}
contour(xx, yy, zz, axes = FALSE, ...)
# see src/library/graphics/R/contour.R
xaxp <- par("xaxp")
xat <- seq(xaxp[1], xaxp[2], length.out = xaxp[3])
xlabels <- format(xat)
yaxp <- par("yaxp")
yat <- seq(yaxp[1], yaxp[2], length.out = yaxp[3])
ylabels <- format(yat)
xIsTime <- inherits(x, "POSIXt") || inherits(x, "POSIXct") || inherits(x, "POSIXlt")
if (revx) {
if (xIsTime) {
oce.axis.POSIXct(
side = 1, x = x,
drawTimeRange = drawTimeRange,
# mar=mar, mgp=mgp,
tformat = tformat, debug = debug - 1
)
} else {
Axis(xx, side = 1, at = xat, labels = rev(xlabels))
}
} else {
if (xIsTime) {
oce.axis.POSIXct(
side = 1, x = x,
drawTimeRange = drawTimeRange,
# mar=mar, mgp=mgp,
tformat = tformat, debug = debug - 1
)
} else {
Axis(xx, side = 1)
}
}
if (revy) {
Axis(yy, side = 2, at = yat, labels = rev(ylabels))
} else {
Axis(yy, side = 2)
}
box()
}
}
oceContour <- oce.contour
#' Partial Matching of Strings or Numbers
#'
#' An extended version of [pmatch()] that allows `x` to be
#' numeric or string-based. As with [pmatch()], partial string
#' matches are handled.
#' This is a wrapper that is useful mainly for `which` arguments to
#' plotting functions.
#'
#' @aliases oce.pmatch
#'
#' @param x a code, or vector of codes. This may be numeric, in which case it
#' is simply returned without further analysis of the other arguments, or it
#' may be string-based, in which case [pmatch()] is used to find
#' numeric matches.
#'
#' @param table a list that maps strings to numbers; [pmatch()] is
#' used on `names(table)`. If the name contains characters that are
#' normally not permitted in a variable name, use quotes, e.g.
#' `list(salinity=1, temperature=2, "salinity+temperature"=3)`.
#'
#' @param nomatch value to be returned for cases of no match (passed to
#' [pmatch()].
#'
#' @param duplicates.ok code for the handling of duplicates (passed to
#' [pmatch()]).
#'
#' @return A number, or vector of numbers, corresponding to the matches.
#' Non-matches are indicated with `NA` values, or whatever value is given
#' by the `NA` argument.
#'
#' @author Dan Kelley
#'
#' @seealso Since [pmatch()] is used for the actual matching, its
#' documentation should be consulted.
#'
#' @examples
#' library(oce)
#' oce.pmatch(c("s", "at", "te"), list(salinity = 1, temperature = 3.1))
ocePmatch <- function(x, table, nomatch = NA_integer_, duplicates.ok = FALSE) {
# FIXME: do element by element, and extend as follows, to allow string numbers
# if (1 == length(grep("^[0-9]*$", ww))) which2[w] <- as.numeric(ww)
if (is.numeric(x)) {
return(x)
} else if (is.character(x)) {
nx <- length(x)
res <- NULL
for (i in 1:nx) {
if (1 == length(grep("^[0-9]*$", x[i]))) {
res <- c(res, as.numeric(x[i]))
} else {
ix <- pmatch(x[i], names(table), nomatch = nomatch, duplicates.ok = duplicates.ok)
# use unlist() to expand e.g. list(x=1:10)
res <- c(res, if (is.na(ix)) NA else unlist(table[[ix]]))
}
}
# m <- pmatch(x, names(table), nomatch=nomatch, duplicates.ok=duplicates.ok)
# return(as.numeric(table)[m])
return(as.numeric(res))
} else {
stop("'x' must be numeric or character")
}
}
oce.pmatch <- ocePmatch
#' Normalize a Spectrum
#'
#' This is a wrapper around the R [spectrum()] function, which
#' returns spectral values that are adjusted so that the integral of those
#' values equals the variance of the input `x`.
#'
#' @aliases oce.spectrum oceSpectrum
#'
#' @param x a univariate or multivariate time series, as for [spectrum()].
#'
#' @param \dots extra arguments passed on to [spectrum()].
#'
#' @return A spectrum that has values that integrate to the variance.
#'
#' @author Dan Kelley
#'
#' @seealso [spectrum()].
#'
#' @examples
#' x <- rnorm(1e3)
#' s <- spectrum(x, plot = FALSE)
#' ss <- oce.spectrum(x, plot = FALSE)
#' cat("variance of x=", var(x), "\n")
#' cat("integral of spectrum=", sum(s$spec) * diff(s$freq[1:2]), "\n")
#' cat("integral of oce.spectrum=", sum(ss$spec) * diff(ss$freq[1:2]), "\n")
oceSpectrum <- function(x, ...) {
args <- list(...)
want.plot <- FALSE
if ("plot" %in% names(args)) {
want.plot <- args$plot
args$plot <- FALSE
args$x <- x
res <- do.call(spectrum, args)
}
dt <- diff(res$freq[1:2])
normalize <- var(x) / (sum(res$spec) * dt)
res$spec <- normalize * res$spec
if (want.plot) {
plot(res)
}
invisible(res)
}
oce.spectrum <- oceSpectrum
#' Filter a Time Series
#'
#' Filter a time-series, possibly recursively
#'
#' The filter is defined as e.g. \eqn{y[i]=b[1]*x[i] + }{y[i]=b[1]*x[i] +
#' b[2]*x[i-1] + b[3]*x[i-2] + ... - a[2]*y[i-1] - a[3]*y[i-2] - a[4]*y[i-3] -
#' ...}\eqn{ b[2]*x[i-1] + b[3]*x[i-2] + ... - a[2]*y[i-1] - a[3]*y[i-2] -
#' }{y[i]=b[1]*x[i] + b[2]*x[i-1] + b[3]*x[i-2] + ... - a[2]*y[i-1] -
#' a[3]*y[i-2] - a[4]*y[i-3] - ...}\eqn{ a[4]*y[i-3] - ...}{y[i]=b[1]*x[i] +
#' b[2]*x[i-1] + b[3]*x[i-2] + ... - a[2]*y[i-1] - a[3]*y[i-2] - a[4]*y[i-3] -
#' ...}, where some of the illustrated terms will be omitted if the lengths of
#' `a` and `b` are too small, and terms are dropped at the start of
#' the time series where the index on `x` would be less than 1.
#'
#' By contrast with the [filter()] function of R, `oce.filter`
#' lacks the option to do a circular filter. As a consequence,
#' `oceFilter` introduces a phase lag. One way to remove this lag is to
#' run the filter forwards and then backwards, as in the \dQuote{Examples}.
#' However, the result is still problematic, in the sense that applying it in
#' the reverse order would yield a different result. (Matlab's `filtfilt`
#' shares this problem.)
#'
#' @aliases oce.filter
#'
#' @param x a vector of numeric values, to be filtered as a time series.
#'
#' @param a a vector of numeric values, giving the \eqn{a}{a} coefficients (see
#' \dQuote{Details}).
#'
#' @param b a vector of numeric values, giving the \eqn{b}{b} coefficients (see
#' \dQuote{Details}).
#'
#' @param zero.phase boolean, set to `TRUE` to run the filter forwards,
#' and then backwards, thus removing any phase shifts associated with the
#' filter.
#'
#' @return A numeric vector of the filtered results, \eqn{y}{y}, as denoted in
#' \dQuote{Details}.
#'
#' @note The first value in the `a` vector is ignored, and if
#' `length(a)` equals 1, a non-recursive filter results.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' par(mar = c(4, 4, 1, 1))
#' b <- rep(1, 5) / 5
#' a <- 1
#' x <- seq(0, 10)
#' y <- ifelse(x == 5, 1, 0)
#' f1 <- oceFilter(y, a, b)
#' plot(x, y, ylim = c(-0, 1.5), pch = "o", type = "b")
#' points(x, f1, pch = "x", col = "red")
#'
#' # remove the phase lag
#' f2 <- oceFilter(y, a, b, TRUE)
#' points(x, f2, pch = "+", col = "blue")
#'
#' legend("topleft",
#' col = c("black", "red", "blue"), pch = c("o", "x", "+"),
#' legend = c("data", "normal filter", "zero-phase filter")
#' )
#' mtext("note that normal filter rolls off at end")
oceFilter <- function(x, a = 1, b, zero.phase = FALSE) {
if (missing(x)) {
stop("must supply x")
}
if (missing(b)) {
stop("must supply b")
}
if (!zero.phase) {
return(do_oce_filter(x, a, b))
}
res <- do_oce_filter(x, a, b)
res <- rev(res)
res <- do_oce_filter(res, a, b)
return(rev(res))
}
oce.filter <- oceFilter
#' Smooth an oce Object
#'
#' Each data element is smoothed as a timeseries. For ADP data, this is done
#' along time, not distance. Time vectors, if any, are not smoothed. A good
#' use of `oce.smooth` is for despiking noisy data.
#'
#' @aliases oce.smooth
#'
#' @param x an [oce-class] object.
#'
#' @param \dots parameters to be supplied to [smooth()], which does
#' the actual work.
#'
#' @return An [oce-class] object that has been smoothed appropriately.
#'
#' @author Dan Kelley
#'
#' @seealso The work is done with [smooth()], and the `...`
#' arguments are handed to it directly by `oce.smooth`.
#'
#' @examples
#' library(oce)
#' data(ctd)
#' d <- oce.smooth(ctd)
#' plot(d)
oceSmooth <- function(x, ...) {
if (!inherits(x, "oce")) {
stop("method is only for oce objects")
}
res <- x
if (inherits(x, "adp")) {
stop("cannot smooth ADP objects (feel free to request this from the author)")
} else if (inherits(x, "adv")) {
for (name in names(x@data)) {
if (length(grep("^time", name))) {
next
}
if (is.vector(x@data[[name]])) {
oceDebug(debug, "smoothing x@data$", name, ", which is a vector\n", sep = "")
res@data[[name]] <- smooth(x@data[[name]], ...)
} else if (is.matrix(x@data[[name]])) {
for (j in seq_len(dim(x@data[[name]])[2])) {
oceDebug(debug, "smoothing x@data[[", name, ",", j, "]], which is a matrix\n", sep = "")
res@data[[name, j]] <- smooth(x@data[[name, j]], ...)
}
} else if (is.array(x@data[[name]])) {
dim <- dim(x@data[[name]])
for (k in seq_len(dim[2])) {
for (j in seq_len(dim[3])) {
oceDebug(debug, "smoothing x@data[[", name, ",", j, "]], which is an arry \n", sep = "")
res@data[[name, j, k]] <- smooth(x@data[[name, j, k]], ...)
}
}
}
}
warning("oce.smooth() has recently been recoded for 'adv' objects -- do not trust it yet!")
} else if (inherits(x, "ctd")) {
res <- x
for (name in names(x@data)) {
res@data[[name]] <- smooth(x@data[[name]], ...)
}
} else {
stop("smoothing does not work (yet) for objects of class ", paste(class(x), collapse = " "))
}
res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep = "", collapse = ""))
res
}
oce.smooth <- oceSmooth
#' Print a Debugging Message
#'
#' Print an indented debugging message.
#' Many oce functions decrease the `debug` level by 1 when they call other
#' functions, so the effect is a nesting, with more space for deeper function
#' level.
#'
#' @aliases oce.debug
#'
#' @param debug an integer, less than or equal to zero for no message, and
#' greater than zero for increasing levels of debugging. Values greater than 4
#' are treated like 4.
#'
#' @param \dots items to be supplied to [cat()], which does the
#' printing. Note that no newline will be printed unless \dots
#' contains a string with a newline character (as in the example).
#'
#' @param style either a string or a function. If a string,
#' it must be `"plain"` (the default) for plain text,
#' `"bold"`, `"italic"`, `"red"`, `"green"` or `"blue"` (with
#' obvious meanings). Note that none of these has any effect
#' for non-interactive use, because doing so would make it difficult
#' to work with R-markdown and similar documents that are to be run
#' through latex.
#'
#' If `style` is a function, it must prepend and postpend the text
#' with control codes, as in the cyan-coloured example; note that
#' \CRANpkg{crayon} provides many functions that work well for `style`.
#'
#' @param unindent integer giving the number of levels to un-indent,
#' e.g. for start and end lines from a called function.
#'
#' @param sep character to insert between elements of `...`, by passing
#' it to [cat()].
#'
#' @author Dan Kelley
#'
#' @examples
#' oceDebug(debug = 1, "Example", 1, "Plain text")
#' oceDebug(debug = 1, "Example", 2, "Bold", style = "bold")
#' oceDebug(debug = 1, "Example", 3, "Italic", style = "italic")
#' oceDebug(debug = 1, "Example", 4, "Red", style = "red")
#' oceDebug(debug = 1, "Example", 5, "Green", style = "green")
#' oceDebug(debug = 1, "Example", 6, "Blue", style = "blue")
#' mycyan <- function(...) paste("\033[36m", paste(..., sep = " "), "\033[0m", sep = "")
#' oceDebug(debug = 1, "Example", 7, "User-set cyan", style = mycyan)
oceDebug <- function(debug = 0, ..., style = "plain", unindent = 0, sep = "") {
catSpecial <- function(...) if (interactive()) cat(...)
debug <- if (debug > 4) 4 else max(0, floor(debug + 0.5))
if (debug > 0) {
n <- 5 - debug - unindent
if (is.character(style) && style == "plain") {
if (n > 0) {
cat(paste(rep(" ", n), collapse = ""))
}
cat(..., sep = sep)
} else if (is.character(style) && style == "bold") {
catSpecial("\033[1m")
if (n > 0) {
cat(paste(rep(" ", n), collapse = ""))
}
cat(..., sep = sep)
catSpecial("\033[0m")
} else if (is.character(style) && style == "italic") {
catSpecial("\033[3m")
if (n > 0) {
cat(paste(rep(" ", n), collapse = ""))
}
cat(..., sep = sep)
catSpecial("\033[0m")
} else if (is.character(style) && style == "red") {
catSpecial("\033[31m")
if (n > 0) {
cat(paste(rep(" ", n), collapse = ""))
}
cat(..., sep = sep)
catSpecial("\033[0m")
} else if (is.character(style) && style == "green") {
catSpecial("\033[32m")
if (n > 0) {
cat(paste(rep(" ", n), collapse = ""))
}
cat(..., sep = sep)
catSpecial("\033[0m")
} else if (is.character(style) && style == "blue") {
catSpecial("\033[34m")
if (n > 0) {
cat(paste(rep(" ", n), collapse = ""))
}
cat(..., sep = sep)
catSpecial("\033[0m")
} else if (is.function(style)) {
if (n > 0) {
cat(style(paste(rep(" ", n), collapse = "")))
}
cat(style(...), sep = sep)
} else { # fallback
if (n > 0) {
cat(paste(rep(" ", n), collapse = ""))
}
cat(..., sep = sep)
}
flush.console()
}
invisible(NULL)
}
oce.debug <- oceDebug
#' Version of as.raw() That Clips Data
#'
#' A version of as.raw() that clips data to prevent warnings
#'
#' Negative values are clipped to 0, while values above 255 are clipped to 255;
#' the result is passed to [as.raw()] and returned.
#'
#' @param x values to be converted to raw
#'
#' @return Raw values corresponding to `x`.
#'
#' @author Dan Kelley
#'
#' @examples
#' x <- c(-0.1, 0, 1, 255, 255.1)
#' data.frame(x, oce.as.raw(x))
oce.as.raw <- function(x) {
# prevent warnings from out-of-range with as.raw()
na <- is.na(x)
x[na] <- 0 # FIXME: what to do here?
x <- ifelse(x < 0, 0, x)
x <- ifelse(x > 255, 255, x)
x <- as.raw(x)
x
}
#' Convolve Two Time Series
#'
#' Convolve two time series, using a backward-looking method.
#' This function provides a straightforward convolution, which may be useful to
#' those who prefer not to use [convolve()] and `filter` in the
#' `stats` package.
#'
#' @aliases oce.convolve
#'
#' @param x a numerical vector of observations.
#'
#' @param f a numerical vector of filter coefficients.
#'
#' @param end a flag that controls how to handle the points of the `x`
#' series that have indices less than the length of `f`. If `end=0`,
#' the values are set to 0. If `end=1`, the original x values are used
#' there. If `end=2`, that fraction of the `f` values that overlap
#' with `x` are used.
#'
#' @return A vector of the convolution output.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' t <- 0:1027
#' n <- length(t)
#' signal <- ifelse(sin(t * 2 * pi / 128) > 0, 1, 0)
#' tau <- 10
#' filter <- exp(-seq(5 * tau, 0) / tau)
#' filter <- filter / sum(filter)
#' observation <- oce.convolve(signal, filter)
#' plot(t, signal, type = "l")
#' lines(t, observation, lty = "dotted")
oceConvolve <- function(x, f, end = 2) {
do_oce_convolve(x, f, end)
}
oce.convolve <- oceConvolve
#' Lowpass Digital Filtering
#'
#' The filter coefficients are constructed using standard definitions,
#' and then [stats::filter()] is
#' used to filter the data. This leaves `NA`
#' values within half the filter length of the ends of the time series, but
#' these may be replaced with the original `x` values, if the argument
#' `replace` is set to `TRUE`.
#'
#' @section Caution: This function was added in June of 2017,
#' and it may be extended during the rest of 2017. New arguments
#' may appear between `n` and `replace`, so users are
#' advised to call this function with named arguments, not positional
#' arguments.
#'
#' @param x a vector to be smoothed
#'
#' @param filter name of filter; at present, `"hamming"`, `"hanning"`, and `"boxcar"` are permitted.
#'
#' @param n length of filter (must be an odd integer exceeding 1)
#'
#' @param replace a logical value indicating whether points near the
#' ends of `x` should be copied into the end regions, replacing
#' the `NA` values that would otherwise be placed there by
#' [stats::filter()].
#'
#' @param coefficients logical value indicating whether to return
#' the filter coefficients, instead of the filtered values. In accordance
#' with conventions in the literature, the returned values are not
#' normalized to sum to 1, although of course that normalization
#' is done in the actual filtering.
#'
#' @return By default, `lowpass` returns a filtered version
#' of `x`, but if `coefficients` is `TRUE` then it
#' returns the filter coefficients.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' par(mfrow = c(1, 2), mar = c(4, 4, 1, 1))
#' coef <- lowpass(n = 5, coefficients = TRUE)
#' plot(-2:2, coef, ylim = c(0, 1), xlab = "Lag", ylab = "Coefficient")
#' x <- seq(-5, 5) + rnorm(11)
#' plot(1:11, x, type = "o", xlab = "time", ylab = "x and X")
#' X <- lowpass(x, n = 5)
#' lines(1:11, X, col = 2)
#' points(1:11, X, col = 2)
lowpass <- function(x, filter = "hamming", n, replace = TRUE, coefficients = FALSE) {
# .Call("hammingFilter", x, n)
if (missing(x) && !coefficients) {
stop("must supply x")
}
if (missing(n)) {
stop("must supply n")
}
if (n < 1) {
stop("n must be be an integer exceeding 1")
}
n2 <- n %/% 2 # half width
if (2 * n2 == n) {
stop("n must be an odd integer")
}
twopi <- 8 * atan2(1, 1)
ii <- n2 + seq.int(-n2, n2, 1)
if (filter == "hamming") {
f <- 0.54 - 0.46 * cos(twopi * ii / (n - 1))
} else if (filter == "hanning") {
f <- 0.5 * (1 - cos(twopi * ii / (n - 1)))
} else if (filter == "boxcar") {
f <- rep(1 / n, n)
} else {
stop("filter must be \"hanning\", \"hamming\", or \"boxcar\"")
}
if (coefficients) {
rval <- f
} else {
f <- f / sum(f)
rval <- as.numeric(stats::filter(x = x, filter = f, method = "convolution"))
if (replace) {
start <- seq.int(1, n2)
rval[start] <- x[start]
nx <- length(x)
end <- seq.int(nx - n2 + 1, nx)
rval[end] <- x[end]
}
}
rval
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.