R/misc.R

Defines functions labelWithUnit argShow approx3d angleRemap abbreviateVector gappyIndex pluralize asNumeric

Documented in angleRemap approx3d argShow gappyIndex labelWithUnit

# vim:textwidth=80:expandtab:shiftwidth=4:softtabstop=4

# Internal function, mainly used in plotting raw images

asNumeric <- function(x)
{
    if (is.numeric(x))
        return(x)
    rval <- as.numeric(x)
    if (is.array(x))
        dim(rval) <- dim(x)
    rval
}

# Internal function; use as e.g. oce:::pluralize().  Supply n and singular,
# but only give plural if it is different from singular with a 's' tacked on.
pluralize <- function(n, singular, plural)
{
    if (missing(plural))
        plural <- paste0(singular, "s")
    paste(n, if (n > 1L) plural else singular)
}

#' Create a Possibly Gappy Indexing Vector
#'
#' This is used internally to construct indexing arrays, mainly for adv and adp
#' functions, in which [readBin()] is used to access isolated regions within a
#' [raw] vector. The work is done in C++, for speed. Since this function is
#' designed for use within oce, it does not offer many safeguards on the
#' parameters, beyond detecting an overlapping situation that would occur if
#' `length` exceeded the space between `starts` elements.  Also, users ought
#' to be aware that the behaviour of `gappyIndex()` might change at any time;
#' simply stated, it is not intended for direct use except by the package
#' developers.
#'
#' For example, suppose data elements in a buffer named `buf` start at bytes
#' 1000 and 2000, and that the goal is to skip the first 4 bytes of each of
#' these sequences, and then to read the next 2 bytes as an unsigned 16-bit
#' integer. This could be accomplished as follows.
#'
#'```
#' library(oce)
#' buf <- readBin("filename", "raw", n=5000, size=1)
#' i <- gappyIndex(c(1000, 2000, 3000), 4, 2)
#' # i is 1004,1005, 2004,2005, 3004,3005
#' values <- readBin(buf[i], "integer", size=2, n=3, endian="little")
#'```
#'
#' @param starts integer vector of one or more values.
#'
#' @param offset integer value indicating the value to be added
#' to each of the `starts` value, as the beginning of the sequence.
#'
#' @param length integer value indicating the number of
#' elements of that sequence.
#'
#' @author Dan Kelley
gappyIndex <- function(starts, offset=0L, length=4L)
{
    if (missing(starts))
        stop("must provide 'starts', an integer vector")
    if (any(starts < 1L)) {
        w <- which(starts < 1L)[1]
        stop("'starts' must consist of positive values, but e.g. starts[", w, "] is ", starts[w])
    }
    if (length(offset) != 1L)
        stop("'offset' must be a single number")
    if (length(length) != 1L)
        stop("'length' must be a single number")
    if (offset < 0L)
        stop("'offset' must be non-negative")
    if (length < 1L)
        stop("'length' must be positive")
    do_gappy_index(starts, offset, length)
}

#OLD #' Create a Possibly Gappy Indexing Vector
#OLD #'
#OLD #' This is used internally to construct indexing arrays, mainly for adv and adp
#OLD #' functions, in which [readBin()] is used to access isolated regions within a
#OLD #' [raw] vector. The work is done in C++, for speed.
#OLD #'
#OLD #' For example, suppose data elements in a buffer named `buf` start at bytes
#OLD #' 1000 and 2000, and that the goal is to read the first 4 bytes of each of
#OLD #' these sequences as an unsigned, 4-byte integer.  This could be accomplished
#OLD #' as follows.  Note that the last argument to `gappyIndex()` is 3 in this
#OLD #' example, because the `0:3` has 4 elements.
#OLD #'
#OLD #'```
#OLD #' library(oce)
#OLD #' buf <- readBin("filename", "raw", n=5000, size=1)
#OLD #' i <- gappyIndex(c(1000, 2000), 0, 3) # Note the length(0:3) is 4.
#OLD #' # i is 1000,1001,1002,1003,2000,2001,2002,2003
#OLD #' readBin(buf[i], "integer", size=4, n=1, endian="little")
#OLD #'```
#OLD #'
#OLD #' @param starts integer vector of one or more values.
#OLD #'
#OLD #' @param from,to integer values controlling the generated sequence, as if
#OLD #' [seq(from,to)] had been used. See \dQuote{Details}.
#OLD #'
#OLD #' @author Dan Kelley
#OLD gappyIndex <- function(starts, from, to)
#OLD {
#OLD     if (missing(starts)) stop("must provide 'starts', an integer vector")
#OLD     if (any(starts < 1L)) stop("'starts' must consist of positive values")
#OLD     if (missing(from)) stop("must provide 'from', an integer value")
#OLD     if (missing(to)) stop("must provide 'to', an integer value")
#OLD     if (to <= from) stop("'from' must exceed 'to'")
#OLD     if (from < 0) stop("'from' must be positive")
#OLD     do_gappy_index(starts, from, to)
#OLD }

abbreviateVector <- function(x)
{
    if (1 >= length(x))
        return(x)
    ud <- unique(diff(x))
    if (1L == length(ud) && 1L == ud) paste(x[1], ":", tail(x, 1), sep="") else x
}


#' Convert angles from 0:360 to -180:180
#'
#' This is mostly used for instrument heading angles, in cases where the
#' instrument is aligned nearly northward, so that small variations in heading
#' (e.g. due to mooring motion) can yield values that swing from small angles
#' to large angles, because of the modulo-360 cut point.
#' The method is to use the cosine and sine of the angle in order to find "x"
#' and "y" values on a unit circle, and then to use [atan2()] to
#' infer the angles.
#'
#' @param theta an angle (in degrees) that is in the range from 0 to 360
#' degrees
#' @return A vector of angles, in the range -180 to 180.
#' @author Dan Kelley
#' @examples
#'
#' library(oce)
#' # fake some heading data that lie near due-north (0 degrees)
#' n <- 20
#' heading <- 360 + rnorm(n, sd=10)
#' heading <- ifelse(heading > 360, heading - 360, heading)
#' x <- 1:n
#' plot(x, heading, ylim=c(-10, 360), type="l", col="lightgray", lwd=10)
#' lines(x, angleRemap(heading))
angleRemap <- function(theta)
{
    toRad <- atan2(1, 1) / 45
    atan2(sin(toRad * theta), cos(toRad * theta)) / toRad
}

#' Trilinear interpolation in a 3D array
#'
#' Interpolate within a 3D array, using the trilinear approximation.
#'
#' Trilinear interpolation is used to interpolate within the `f` array,
#' for those (`xout`, `yout` and `zout`) triplets that are
#' inside the region specified by `x`, `y` and `z`.  Triplets
#' that lie outside the range of `x`, `y` or `z` result in
#' `NA` values.
#'
#' @param x vector of x values for grid (must be equi-spaced)
#' @param y vector of y values for grid (must be equi-spaced)
#' @param z vector of z values for grid (must be equi-spaced)
#' @param f matrix of rank 3, with the gridded values mapping to the `x`
#' values (first index of `f`), etc.
#' @param xout vector of x values for output.
#' @param yout vector of y values for output (length must match that of
#' `xout`).
#' @param zout vector of z values for output (length must match that of
#' `xout`).
#' @return A vector of interpolated values (or `NA` values), with length
#' matching that of `xout`.
#' @author Dan Kelley and Clark Richards
#'
#' @examples
#' # set up a grid
#' library(oce)
#' n <- 5
#' x <- seq(0, 1, length.out=n)
#' y <- seq(0, 1, length.out=n)
#' z <- seq(0, 1, length.out=n)
#' f <- array(1:n^3, dim=c(length(x), length(y), length(z)))
#' # interpolate along a diagonal line
#' m <- 100
#' xout <- seq(0, 1, length.out=m)
#' yout <- seq(0, 1, length.out=m)
#' zout <- seq(0, 1, length.out=m)
#' approx <- approx3d(x, y, z, f, xout, yout, zout)
#' # graph the results
#' plot(xout, approx, type="l")
#' points(xout[1], f[1, 1, 1])
#' points(xout[m], f[n, n, n])
approx3d <- function(x, y, z, f, xout, yout, zout)
{
    # Were all arguments given?
    if (missing(x))
        stop("must provide x")
    if (missing(y))
        stop("must provide y")
    if (missing(z))
        stop("must provide z")
    if (missing(f))
        stop("must provide f")
    if (missing(xout))
        stop("must provide xout")
    if (missing(yout))
        stop("must provide yout")
    if (missing(zout))
        stop("must provide zout")
    # Are there enough data to interpolate?
    if (length(x) < 2)
        stop("must have more than one x value")
    if (length(y) < 2)
        stop("must have more than one y value")
    if (length(z) < 2)
        stop("must have more than one z value")
    # Are the array dimensions consistent with x, y, and z?
    if (3 != length(dim(f)))
        stop("f must be an array with 3 dimensions")
    if (length(x) != dim(f)[1])
        stop("length of x and dim(f)[1] must agree, but they are ", length(x), " and ", dim(f)[1])
    if (length(y) != dim(f)[2])
        stop("length of y and dim(f)[2] must agree, but they are ", length(y), " and ", dim(f)[2])
    if (length(z) != dim(f)[3])
        stop("length of z and dim(f)[3] must agree, but they are ", length(z), " and ", dim(f)[3])
    # Are x, y and z equi-spaced?
    if (length(x) > 2) {
        equispaced <- function(a) sd(diff(a)) / mean(diff(a)) < 1e-5
        if (!equispaced(x))
            stop("x values must be equi-spaced")
        if (!equispaced(y))
            stop("y values must be equi-spaced")
        if (!equispaced(z))
            stop("z values must be equi-spaced")
    }
    do_approx3d(x, y, z, f, xout, yout, zout)
}


#' Show an argument to a function, e.g. for debugging
#'
#' @param x the argument
#'
#' @param nshow number of values to show at first (if length(x)> 1)
#'
#' @param last indicates whether this is the final argument to the function
#'
#' @param sep the separator between name and value
argShow <- function(x, nshow=4, last=FALSE, sep="=")
{
    if (missing(x))
        return("")
    name <- paste(substitute(expr=x, env=environment()))
    res <- ""
    nx <- length(x)
    if (missing(x)) {
        res <- "(missing)"
    } else {
        if (is.null(x)) {
            res <- NULL
        } else {
            if (is.function(x)) {
                res <- "(provided)"
            } else if (nx==1) {
                res <- if (is.character(x)) paste("\"", x[1], "\"", sep="") else x[1]
            } else {
                look <- seq.int(1L, min(nshow, nx)-1)
                res <- paste(format(x[look], digits=4), collapse=",")
                res <- paste(res, if (nx > nshow) ", ..., " else ",", x[nx], sep="")
            }
        }
    }
    res <- if (nx == 1) paste(name, "=", res,  sep="") else paste(name, "=c(", res, ")", sep="")
    if (!last) {
        res <- paste(res, ",", sep="")
    }
    res
}

#' Create label with unit
#'
#' `labelWithUnit` creates a label with a unit, for graphical
#' display, e.g. by \code{\link{plot,section-method}}.
#' The unit is enclosed in square brackets, although setting
#' `options(oceUnitBracket="(")` will cause parentheses to be
#' used, instead.
#'
#' If `name` is in a standard list, then alterations are made as appropriate,
#' e.g. `"SA"` or `"Absolute Salinity"` yields an S with subscript A; `"CT"` or
#' `"Conservative Temperature"` yields an upper-case Theta, `sigmaTheta`
#' yields a sigma with subscript theta, `sigma0` yields
#' sigma with subscript 0 (with similar for 1 through 4), `"N2"` yields "N" with
#' superscript 2, and `"pressure"` yields "p".
#' These basic hydrographic quantities have default units that will
#' be used if `unit` is not supplied (see \dQuote{Examples}).
#'
#' In addition to the above, several chemical names are recognized,
#' but no unit is guessed for them, because the oceanographic
#' community lacks agreed-upon standards.
#'
#' If `name` is not recognized, then it is simply repeated in the
#' return value.
#'
#' @param name character value naming a quantity.
#'
#' @param unit a list containing items `unit` and (optionally) `scale`, only the
#' first of which, an [expression()], is used.  If `unit` is not provided,
#' then a default will be used (see \dQuote{Details}).
#'
#' @return `labelWithUnit` returns a language object, created with [bquote()],
#' that that may supplied as a text string to [legend()], [mtext()], [text()],
#' etc.
#'
#' @examples
#' library(oce)
#' # 1. temperature has a predefined unit, but this can be overruled
#' labelWithUnit("temperature")
#' labelWithUnit("temperature",
#'     list(unit=expression(m/s), scale="erroneous"))
#' # 2. phosphate lacks a predefined unit
#' labelWithUnit("phosphate")
#' data(section)
#' labelWithUnit("phosphate",
#'     section[["station",1]][["phosphateUnit"]])
#'
#' @family functions that create labels
#'
#' @author Dan Kelley
labelWithUnit <- function(name, unit=NULL)
{
    u <- if (!is.null(unit) && length(unit$unit) > 0L) unit$unit[[1]] else "unitless"
    L <- if (getOption("oceUnitBracket", "[") == "[") " [" else " ("
    R <- if (getOption("oceUnitBracket", "[") == "[")  "]" else  ")"
    # Note that the code is alphabetical in the first item of
    # equivalents. Please follow that convention if adding new
    # entries. Also, use paste() to combine words, to prevent
    # text-editor reflow operations from inserting line breaks.
    #
    # There's no need to handle unitless quantities that don't
    # need renaming, since they are handled by the final else.
    if (name %in% c(paste("Absolute", "Salinity"), "SA")) {
        rval <- if (is.null(unit)) bquote(S[A]*.(L)*g/kg*.(R)) else bquote(S[A]*.(L)*.(u)*.(R))
    } else if (name %in% c(paste("Conservative", "Temperature"), "CT")) {
        rval <- if (is.null(unit)) bquote(Theta*.(L)*degree*C*.(R)) else bquote(Theta*.(L)*.(u)*.(R))
    } else if (name %in% c("density")) {
        rval <- if (is.null(unit)) bquote(rho*.(L)*kg/m^3*.(R)) else bquote(rho*.(L)*.(u)*.(R))
    } else if (name %in% c("depth")) {
        rval <- if (is.null(unit)) bquote(depth*.(L)*m*.(R)) else bquote(depth*.(L)*.(u)*.(R))
    } else if (name %in% "N2") {
        rval <- if (is.null(unit)) bquote(N^2*.(L)*s^-2*.(R)) else bquote(N^2*.(L)*.(u)*.(R))
    } else if (name == "nitrate") {
        rval <- if (is.null(unit)) bquote(NO[3]) else bquote(NO[3]*.(L)*.(u)*.(R))
    } else if (name == "nitrite") {
        rval <- if (is.null(unit)) bquote(NO[2]) else bquote(NO[2]*.(L)*.(u)*.(R))
    } else if (name == "NO2+NO3") {
        rval <- if (is.null(unit)) bquote(NO[2]+NO[3]) else bquote(NO[2]+NO[3]*.(L)*.(u)*.(R))
    } else if (name == "oxygen") {
        rval <- if (is.null(unit)) bquote(O[2]) else bquote(O[2]*.(L)*.(u)*.(R))
    } else if (name == "pressure") {
        rval <- if (is.null(unit)) bquote(p*.(L)*dbar*.(R)) else bquote(p*.(L)*.(u)*.(R))
    } else if (name == "phosphate") {
        rval <- if (is.null(unit)) bquote(PO[4]) else bquote(PO[4]*.(L)*.(u)*.(R))
    } else if (name %in% c(paste("potential", "temperature"), "theta")) {
        rval <- if (is.null(unit)) bquote(theta*.(L)*degree*C*.(R)) else bquote(theta*.(L)*.(u)*.(R))
    } else if (name == "Rrho") {
        rval <- expression(R[rho])
    } else if (name %in% c("salinity", "SP")) {
        rval <- expression("S")
    } else if (name == "sigma0") {
        rval <- if (is.null(unit)) bquote(sigma[0]*.(L)*kg/m^3*.(R)) else bquote(sigma[0]*.(L)*.(u)*.(R))
    } else if (name == "sigma1") {
        rval <- if (is.null(unit)) bquote(sigma[1]*.(L)*kg/m^3*.(R)) else bquote(sigma[1]*.(L)*.(u)*.(R))
    } else if (name == "sigma2") {
        rval <- if (is.null(unit)) bquote(sigma[2]*.(L)*kg/m^3*.(R)) else bquote(sigma[2]*.(L)*.(u)*.(R))
    } else if (name == "sigma3") {
        rval <- if (is.null(unit)) bquote(sigma[3]*.(L)*kg/m^3*.(R)) else bquote(sigma[3]*.(L)*.(u)*.(R))
    } else if (name == "sigma4") {
        rval <- if (is.null(unit)) bquote(sigma[4]*.(L)*kg/m^3*.(R)) else bquote(sigma[4]*.(L)*.(u)*.(R))
    } else if (name == "sigmaTheta") {
        rval <- if (is.null(unit)) bquote(sigma[theta]*.(L)*kg/m^3*.(R)) else bquote(sigma[theta]*.(L)*.(u)*.(R))
    } else if (name == "silicate") {
        rval <- if (is.null(unit)) bquote(SiO[4]) else bquote(SiO[4]*.(L)*.(u)*.(R))
    } else if (name == "spice") {
        rval <- if (is.null(unit)) bquote(spice*.(L)*kg/m^3*.(R)) else bquote(spice*.(L)*.(u)*.(R))
    } else if (name == "SR") {
        rval <- if (is.null(unit)) bquote(S[R]*.(L)*kg/m^3*.(R)) else bquote(S[R]*.(L)*.(u)*.(R))
    } else if (name == "Sstar") {
        rval <- if (is.null(unit)) bquote(S["*"]*.(L)*kg/m^3*.(R)) else bquote(S["*"]*.(L)*.(u)*.(R))
    } else if (name == "temperature") {
        # nolint start T_and_F_symbol_linter
        rval <- if (is.null(unit)) bquote(T*.(L)*degree*C*.(R)) else bquote(T*.(L)*.(u)*.(R))
        # nolint end T_and_F_symbol_linter
    } else if (name == "z") {
        rval <- if (is.null(unit)) bquote(z*.(L)*m*.(R)) else bquote(z*.(L)*.(u)*.(R))
    } else {
        rval <- name
    }
    rval
}

#' Get first finite value in a vector or array, or NULL if none
#' @param v A numerical vector or array.
firstFinite <- function(v)
{
    if (!is.vector(v))
        v <- as.vector(v)
    first <- which(is.finite(v))
    if (length(first) > 0) v[first[1]] else NULL
}

#' Read a World Ocean Atlas NetCDF File
#'
#' @param file character string naming the file
#' @param name of variable to extract. If not provided, an
#' error message is issued that lists the names of data in the file.
#' @param positive logical value indicating whether `longitude` should
#' be converted to be in the range from 0 to 360, with `name`
#' being shuffled accordingly. This is set to `FALSE` by default,
#' because the usual oce convention is for longitude to range between -180
#' to +180.
#'
#' @template encodingIgnoredTemplate
#'
#' @return A list containing vectors `longitude`, `latitude`,
#' `depth`, and an array with the specified name. If `positive`
#' is true, then `longitude` will be converted to range from 0
#' to 360, and the array will be shuffled accordingly.
#'
# @examples
#\dontrun{
# # Mean SST at 5-degree spatial resolution
# tmn <- read.woa("/data/woa13/woa13_decav_t00_5dv2.nc", "t_mn")
# imagep(tmn$longitude, tmn$latitude, tmn$t_mn[, , 1], zlab="SST")
#}
read.woa <- function(file, name, positive=FALSE, encoding=NA)
{
    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, "'")
    }
    if (!is.character(file))
        stop("'file' must be a character string")
    con <- ncdf4::nc_open(file)
    if (missing(name)) {
        varnames <- names(con$var)
        stop("must supply a name from the list: ", paste(varnames, collapse=", "))
        return(NULL)
    }
    longitude <- as.vector(ncdf4::ncvar_get(con, "lon"))
    latitude <- as.vector(ncdf4::ncvar_get(con, "lat"))
    depth <- as.vector(ncdf4::ncvar_get(con, "depth"))
    field <- ncdf4::ncvar_get(con, name)
    if (positive) {
        lon2 <- ifelse(longitude < 0, longitude + 360, longitude)
        i  <- order(lon2)
        longitude <- longitude[i]
        # Crude method to reorder field on first index, whether it is 2D, 3D or 4D,
        # although I'm not sure that any 4D items occur in the World Ocean Atlas.
        if (is.array(field)) {
            ndim <- length(dim(field))
            if (ndim == 2) {
                field <- field[i, ]
            } else if (ndim == 3) {
                field <- field[i, , ]
            } else if (ndim == 4) {
                field <- field[i, , , ]
            }
        }
    }
    rval <- list(longitude=longitude, latitude=latitude, depth=depth, field=field)
    names(rval) <- c(head(names(rval), -1), name)
    rval
}

# alphabetized functions END


# unalphabetized functions START

shortenTimeString <- function(t, debug=getOption("oceDebug"))
{
    tc <- as.character(t)
    oceDebug(debug, "shortenTimeString() {\n", sep="", unindent=1)
    oceDebug(debug, "A: '", paste(t, collapse="' '"), "'\n")
    tc <- gsub(" [A-Z]{3}$", "", tc) # remove timezone
    if (all(grepl("^[0-9]{4}", tc))) {
        # leading years
        years <- substr(tc, 1, 4)
        if (1 == length(unique(years))) {
            tc <- gsub("^[0-9]{4}", "", tc)
            tc <- gsub("^-", "", tc) # works for ISO dates
            oceDebug(debug, "B: '", paste(tc, collapse="' '"), "'\n", sep="")
        }
    } else if (any(grepl("[a-zA-Z]", tc))) {
        # Change e.g. 'Jul 01' to 'Jul' if all labels end in 01
        if (all(grepl("01\\s*$", tc))) {
            tc <- gsub(" 01\\s*$", "", tc)
            oceDebug(debug, "B: '", paste(tc, collapse="' '"), "'\n", sep="")
        }
    }
    oceDebug(debug, "C: '", paste(tc, collapse="' '"), "'\n", sep="")
    tc <- gsub("^\\s*", "", tc)
    tc <- gsub("\\s*$", "", tc)
    oceDebug(debug, "D: '", paste(tc, collapse="' '"), "'\n", sep="")
    oceDebug(debug, "}\n", unindent=1)
    tc
}

#' Convert each of a vector of strings from SNAKE_CASE to camelCase
#'
#' `snakeToCamel` converts "snake-case" characters such as `"NOVA_SCOTIA"`
#' to "camel-case" values, such as `"NovaScotia"`.  It was written for
#' use by [read.argo()], but it also may prove helpful in other contexts.
#'
#' The basic procedure is to chop the string up into substrings separated by
#' the underline character, then to upper-case the first letter of
#' all substrings except the first, and then to paste the substrings
#' together.
#'
#' However, there are exceptions.  First, any upper-case string that contains no
#' underlines is converted to lower case, but any mixed-case string with no
#' underlines is returned as-is (see the second example). Second, if
#' the `specialCases` argument contains `"QC"`, then the
#' `QC` is passed through directly (since it is an acronym) and
#' if the first letter of remaining text is upper-cased (contrast
#' see the four examples).
#'
#' @param s A vector of character values.
#'
#' @param specialCases A vector of character values that tell which
#' special-cases to apply, or `NULL` (the default) to turn off special
#' cases.  The only permitted special case at the moment is `"QC"` (see
#' \dQuote{Details}) but the idea of this argument is that other cases
#' can be added later, if needed.
#'
#' @return A vector of character values
#'
#' @examples
#' library(oce)
#' snakeToCamel("PARAMETER_DATA_MODE")   # "parameterDataMode"
#' snakeToCamel("PARAMETER")             # "parameter"
#' snakeToCamel("HISTORY_QCTEST")        # "historyQctest"
#' snakeToCamel("HISTORY_QCTEST", "QC")  # "historyQCTest"
#' snakeToCamel("PROFILE_DOXY_QC")       # "profileDoxyQc"
#' snakeToCamel("PROFILE_DOXY_QC", "QC") # "profileDoxyQC"
#' @author Dan Kelley
snakeToCamel <- function(s, specialCases=NULL)
{
    ns <- length(s)
    if ("QC" %in% specialCases) {
        s <- gsub("QCTEST", "Q_C_TEST", s) # for e.g. HISTORY_QCTEST
        s <- gsub("QC$", "Q_C", s)         # for e.g. PROFILE_DOXY_QC
        s <- gsub("Qc$", "QC", s)          # for e.g. positionQc (converted previously)
    }
    if (ns < 1)
        stop("'s' must be a vector of character values")
    res <- vector("character", length=ns)
    for (is in seq(1L, length(s))) {
        if (!grepl("_", s[is])) {
            # Handle the single-word case. If all upper-case, convert to lower,
            # but otherwise, leave it as it is.
            res[is] <- if (s[is] == toupper(s[is])) tolower(s[is]) else s[is]
        } else {
            # Handle the multi-word case. Start by making it lower case.
            s[is] <- tolower(s[is])
            # Now, split and then work through the words
            w <- strsplit(s[is], "_")[[1]]
            nw <- length(w)
            res[is] <- w[1]
            for (iw in seq(2L, nw)) {
                wl <- tolower(w[iw])
                res[is] <- paste0(res[is], toupper(substring(wl, 1, 1)), substring(wl, 2))
            }
        }
    }
    res
}


#' Decode units, from strings
#'
#' This is mainly intended for internal use within the package, e.g. by
#' [read.odf()], and so the list of string-to-unit mappings is not
#' documented, since developers can learn it from simple examination
#' of the code.  The focus of `unitFromString()` is on strings that are
#' found in oceanographic files available to the author, *not* on all
#' possible units.
#'
#' @param unit a character value indicating the unit. These
#' are matched according to rules developed to work with actual
#' data files, and so the list is not by any means exhaustive.
#'
#' @param scale a character value indicating the scale.  The default value
#' of `NULL` dictates that the scale is to be inferred from the unit. If
#' a non-`NULL` value is supplied, it will be used, even if it makes no sense
#' in relation to value of `unit`.
#'
#' @return A [list()] of two items: `unit` which is an
#' [expression()], and `scale`, which is a string.
#'
#' @examples
#' unitFromString("dbar")   # dbar (no scale)
#' unitFromString("deg c")  # modern temperature (ITS-90 scale)
#' @family functions that interpret variable names and units from headers
unitFromString <- function(unit, scale=NULL)
{
    if (length(unit) > 1L)
        stop("cannot work with a vector of strings")
    u <- trimws(unit)                  # remove any leading/trailing whitespace
    U <- toupper(u)                    # simplify some match tests
    #> message("unit=\"", unit, "\", scale=\"", scale, "\"")
    if (U == "" || U == "NONE" || U == "(NONE)")
        return(list(unit=expression(), scale=if (is.null(scale)) "" else scale))
    if (U == "10**3CELLS/L")
        return(list(unit=expression(10^3*cells/l), scale=if (is.null(scale)) "" else scale))
    if (U == "CODE")
        return(list(unit=expression(), scale=if (is.null(scale)) "" else scale))
    if (U == "COUNTS")
        return(list(unit=expression(), scale=if (is.null(scale)) "" else scale))
    if (U == "DB")
        # NOTE: this really should be decibel, but ODF files use it for decibar, sometimes
        return(list(unit=expression(dbar), scale=if (is.null(scale)) "" else scale))
    if (U == "DBAR" || U == "DECIBAR" || U == "DECIBARS")
        return(list(unit=expression(dbar), scale=if (is.null(scale)) "" else scale))
    if (U == "DEG C" || U == "DEGREES C")
        return(list(unit=expression(degree*C), scale=if (is.null(scale)) "ITS-90" else scale))
    if (U == "FMOL/KG")
        return(list(unit=expression(fmol/kg), scale=if (is.null(scale)) "" else scale))
    if (U == "G")
        return(list(unit=expression(g), scale=if (is.null(scale)) "" else scale))
    if (U == "GMT")
        return(list(unit=expression(), scale=if (is.null(scale)) "" else scale))
    if (U == "HPA")
        return(list(unit=expression(hPa), scale=if (is.null(scale)) "" else scale))
    if (U == "HZ")
        return(list(unit=expression(Hz), scale=if (is.null(scale)) "" else scale))
    if (U == "ITS-90" || U == "ITS-90 DEGC")
        return(list(unit=expression(degree*C), scale=if (is.null(scale)) "ITS-90" else scale))
    if (U == "IPTS-68" || U == "IPTS-68 DEGC" || U == "IPTS-68, DEG C")
        return(list(unit=expression(degree*C), scale=if (is.null(scale)) "IPTS-68" else scale))
    if (U == "ITS-68" || U == "ITS-68 DEGC" || U == "ITS-68, DEG C") # not an accepted unit, but seen in ODF files
        return(list(unit=expression(degree*C), scale=if (is.null(scale)) "IPTS-68" else scale))
    if (U == "KG/M^3" || U == "KG/M**3")
        return(list(unit=expression(kg/m^3), scale=if (is.null(scale)) "" else scale))
    if (U == "M" || U == "METER" || U == "METRE" || U == "METERS" || U == "METRES")
        return(list(unit=expression(m), scale=if (is.null(scale)) "" else scale))
    if (U == "M**3/KG" || U == "M^3/KG")
        return(list(unit=expression(m^3/kg), scale=if (is.null(scale)) "" else scale))
    if (U == "MA")
        return(list(unit=expression(ma), scale=if (is.null(scale)) "" else scale))
    if (U == "M/S" || U == "METER/SEC" || U == "M/S" || U == "METRE/SEC")
        return(list(unit=expression(m/s), scale=if (is.null(scale)) "" else scale))
    if (U == "MG/M^3" || U == "MG/M**3")
        return(list(unit=expression(mg/m^3), scale=if (is.null(scale)) "" else scale))
    if (U == "MICRON" || U == "MICRONS")
        return(list(unit=expression(mu*m), scale=if (is.null(scale)) "" else scale))
    if (grepl("^\\s*MICRO[ ]?MOL[E]?S/M(\\*){0,2}2/S(EC)?\\s*$", U))
        return(list(unit=expression(mu*mol/m/s), scale=if (is.null(scale)) "" else scale))
    if (U == "ML/L")
        return(list(unit=expression(ml/l), scale=if (is.null(scale)) "" else scale))
    if (U == "M/S" || U == "M/SEC")
        return(list(unit=expression(m/s), scale=if (is.null(scale)) "" else scale))
    if (U == "M^-1/SR")
        return(list(unit=expression(1/m/sr), scale=if (is.null(scale)) "" else scale))
    if (U == "MHO/CM" || U == "MHOS/CM") # 1 mho (archaic) = 1 Siemen (modern)
        return(list(unit=expression(S/cm), scale=if (is.null(scale)) "" else scale))
    if (U == "MHO/M" || U == "MHOS/M") # 1 mho (archaic) = 1 Siemen (modern)
        return(list(unit=expression(S/m), scale=if (is.null(scale)) "" else scale))
    if (U == "MHO/CM" || U == "MHOS/CM") # 1 mho (archaic) = 1 Siemen (modern)
        return(list(unit=expression(S/cm), scale=if (is.null(scale)) "" else scale))
    if (U == "MMHO") # 1 mho (archaic) = 1 Siemen (modern)
        return(list(unit=expression(mS), scale=if (is.null(scale)) "" else scale))
    if (U == "NBS SCALE")
        return(list(unit=expression(), scale=if (is.null(scale)) "NBS" else scale))
    if (U == "NTU")
        return(list(unit=expression(NTU), scale=if (is.null(scale)) "" else scale))
    if (U == "PPM")
        return(list(unit=expression(ppm), scale=if (is.null(scale)) "" else scale))
    if (U == "PSS-78" || U == "PSU")
        return(list(unit=expression(), scale=if (is.null(scale)) "PSS-78" else scale))
    if (U == "PMOL/KG")
        return(list(unit=expression(pmol/kg), scale=if (is.null(scale)) "" else scale))
    if (U == "PSU")
        return(list(unit=expression(), scale=if (is.null(scale)) "PSS-78" else scale))
    if (U == "ML/L")
        return(list(unit=expression(ml/l), scale=if (is.null(scale)) "" else scale))
    if (U == "S" || U == "SEC" || U == "SECOND")
        return(list(unit=expression(s), scale=if (is.null(scale)) "" else scale))
    if (u == "s/m")
        return(list(unit=expression(s/m), scale=if (is.null(scale)) "" else scale))
    if (u == "S/m")
        return(list(unit=expression(S/m), scale=if (is.null(scale)) "" else scale))
    if (u == "Total scale")
        return(list(unit=expression(), scale=if (is.null(scale)) "" else scale))
    if (u == "True degrees")
        return(list(unit=expression(degree), scale=if (is.null(scale)) "" else scale))
    if (u == "uA")
        return(list(unit=expression(mu*a), scale=if (is.null(scale)) "" else scale))
    # > stringi::stri_escape_unicode() indicates that Greek mu is "\u00b5"
    if (u == "\u00b5einsteins/s/m^2" || U == "UEINSTEINS/S/M**2" || U == "UEINSTEINS/S/M^2")
        return(list(unit=expression(mu*mol/m^2/s), scale=if (is.null(scale)) "" else scale))
    # > stringi::stri_escape_unicode() indicates that Greek mu is "\u00b5"
    if (u == "\u00b5M")
        return(list(unit=expression(mu*M), scale=if (is.null(scale)) "" else scale))
    if (U == "UMOL/KG")
        return(list(unit=expression(mu*mol/kg), scale=if (is.null(scale)) "" else scale))
    if (u == "ug/l")
        return(list(unit=expression(mu*g/l), scale=if (is.null(scale)) "" else scale))
    if (grepl("^mmol/m\\*\\*3$", unit, ignore.case=TRUE))
        return(list(unit=expression(mmol/m^3), scale=if (is.null(scale)) "" else scale))
    if (grepl("^mmol/m\\^3$", unit, ignore.case=TRUE))
        return(list(unit=expression(mmol/m^3), scale=if (is.null(scale)) "" else scale))
    if (U == "UMOL/KG")
        return(list(unit=expression(mmol/kg), scale=if (is.null(scale)) "" else scale))
    if (U == "UMOL/M**3")
        return(list(unit=expression(mu*mol/m^3), scale=if (is.null(scale)) "" else scale))
    if (U == "UMOL/M**2/S")
        return(list(unit=expression(mu*mol/m^2/s), scale=if (is.null(scale)) "" else scale))
    if (U == "UMOL PHOTONS/M2/S")
        return(list(unit=expression(mu*mol/m^2/s), scale=if (is.null(scale)) "" else scale))
    if (U == "UTC")
        return(list(unit=expression(), scale=if (is.null(scale)) "" else scale))
    if (U == "V")
        return(list(unit=expression(V), scale=if (is.null(scale)) "" else scale))
    if (U == "1/CM")
        return(list(unit=expression(1/cm), scale=if (is.null(scale)) "" else scale))
    if (U == "1/M")
        return(list(unit=expression(1/m), scale=if (is.null(scale)) "" else scale))
    if (U == "VOLT" || U == "VOLTS")
        return(list(unit=expression(V), scale=if (is.null(scale)) "" else scale))
    if (U == "%")
        return(list(unit=expression(percent), scale=if (is.null(scale)) "" else scale))
    # If none of the above worked, just try our best.
    return(list(unit=as.expression(unit), scale=if (is.null(scale)) "" else scale))
}

#' Rename duplicated character strings
#'
#' Append numeric suffices to character strings, to avoid repeats.
#' This is used by various data
#' input functions, to handle the fact that several oceanographic data
#' formats permit the reuse of variable names within a given file.
#'
#' @param strings Vector of character strings.
#'
#' @param style An integer giving the style. If `style`
#' is `1`, then e.g. a triplicate of `"a"` yields
#' `"a"`, `"a1"`, and `"a2"`.
#' If `style` is `2`, then the same input yields
#' `"a_001"`, `"a_002"`, and `"a_003"`.
#'
#' @return Vector of strings with repeats distinguished by suffix.
#'
#' @seealso Used by [read.ctd.sbe()] with `style=1` to
#' rename repeated data elements (e.g. for multiple temperature sensors)
#' in CTD data, and by [read.odf()] with `style=2` on
#' key-value pairs within ODF metadata.
#'
#' @examples
#' unduplicateNames(c("a", "b", "a", "c", "b"))
#' unduplicateNames(c("a", "b", "a", "c", "b"), style=2)
unduplicateNames <- function(strings, style=1)
{
    # Handle duplicated names
    if (style == 1) {
        for (i in seq_along(strings)) {
            w <- which(strings == strings[i])
            lw <- length(w)
            if (lw > 1) {
                w <- w[-1]
                strings[w] <- paste(strings[i], 1+seq.int(1, length(w)), sep="")
            }
        }
    } else if (style == 2) {
        for (i in seq_along(strings)) {
            w <- which(strings == strings[i])
            lw <- length(w)
            if (lw > 1) {
                suffices <- seq_len(lw)
                strings[w] <- sprintf("%s_%03d", strings[i], suffices)
            }
        }
    } else {
        stop("unknown style=", style, "; it must be 1 or 2")
    }
    strings
}

#' Calculate a rounded bound, rounded up to mantissa 1, 2, or 5
#'
#' @param x a single positive number
#'
#' @return for positive x, a value exceeding x that has mantissa 1, 2, or 5; otherwise, x
bound125 <- function(x)
{
    x <- x[1] # ignore all but first element
    if (x <= 0) {
        res <- x
    } else {
        exp10 <- 10^floor(log10(x))
        xx <- x / exp10
        m <- if (xx <= 1) 1 else if (xx <=2) 2 else if (xx <= 5) 5 else 10
        res <- m * exp10
    }
    res
}

#' Put longitude in the range from -180 to 180
#'
#' @param longitude in degrees East, possibly exceeding 180
#'
#' @return longitude in signed degrees East
#'
#' @seealso
#' [matrixShiftLongitude()] and [shiftLongitude()] are more
#' powerful relatives to `standardizeLongitude`.
standardizeLongitude <- function(longitude) ifelse(longitude > 180, longitude-360, longitude)

#' Try to associate data names with units, for use by summary()
#'
#' Note that the whole object is not being given as an argument;
#' possibly this will reduce copying and thus storage impact.
#'
#' @param names the names of data within an object
#'
#' @param units the units from metadata
#'
#' @return a vector of strings, with blank entries for data with unknown units
#'
#' @examples
#' library(oce)
#' data(ctd)
#' dataLabel(names(ctd@@data), ctd@@metadata$units)
dataLabel <- function(names, units)
{
    res <- names
    # message("in dataLabel()")
    if (!missing(units)) {
        # message("  dataLabel(); next line is names")
        # print(names)
        # message("  dataLabel(); next line is units")
        # print(units)
        unitsNames <- names(units)
        # message("  dataLabel(); next line is unitsNames")
        # print(unitsNames)
        for (i in seq_along(names)) {
            # message("  i: ", i, ", name: ", names[i])
            w <- which(unitsNames == names[i])
            if (length(w)) {
                # message("  we match a unit at index w=",  paste(w, collapse=" "))
                u <- units[w]
                if (!is.null(u)) {
                    if (is.character(u)) {
                        res[i] <- paste(res[i], " [", u, "]", sep="")
                    } else if (is.list(u)) {
                        res[i] <- paste(res[i], " [", u$unit[[1]], u$scale, "]", sep="")
                    }
                }
            }
        }
    }
    res <- gsub(" *\\[\\]", "", res)
    res
}

#' Capitalize first letter of each of a vector of words
#'
#' This is used in making labels for data names in some ctd functions
#'
#' @param w vector of character strings
#'
#' @return vector of strings patterned on `w` but with first letter
#' in upper case and others in lower case
titleCase <- function(w)
{
    unlist(lapply(seq_along(w),
            function(i) paste(toupper(substr(w[i], 1, 1)),
                tolower(substr(w[i], 2, nchar(w[i]))), sep="")))
}


#' Curl of 2D vector field
#'
#' Calculate the z component of the curl of an x-y vector field.
#'
#' The computed component of the curl is defined by \eqn{\partial }{dv/dx -
#' du/dy}\eqn{ v/\partial x - \partial u/\partial y}{dv/dx - du/dy} and the
#' estimate is made using first-difference approximations to the derivatives.
#' Two methods are provided, selected by the value of `method`.
#'
#' * For `method=1`, a centred-difference, 5-point stencil is used in
#' the interior of the domain.  For example, \eqn{\partial v/\partial x}{dv/dx}
#' is given by the ratio of \eqn{v_{i+1,j}-v_{i-1,j}}{v[i+1,j]-v[i-1,j]} to the
#' x extent of the grid cell at index \eqn{j}{j}. (The cell extents depend on
#' the value of `geographical`.)  Then, the edges are filled in with
#' nearest-neighbour values. Finally, the corners are filled in with the
#' adjacent value along a diagonal.  If `geographical=TRUE`, then `x`
#' and `y` are taken to be longitude and latitude in degrees, and the
#' earth shape is approximated as a sphere with radius 6371km.  The resultant
#' `x` and `y` are identical to the provided values, and the
#' resultant `curl` is a matrix with dimension identical to that of
#' `u`.
#'
#' * For `method=2`, each interior cell in the grid is considered
#' individually, with derivatives calculated at the cell center. For example,
#' \eqn{\partial v/\partial x}{dv/dx} is given by the ratio of
#' \eqn{0.5*(v_{i+1,j}+v_{i+1,j+1}) -
#' 0.5*(v_{i,j}+v_{i,j+1})}{0.5*(v[i+1,j]+v[i+1,j+1]) - 0.5*(v[i,j]+v[i,j+1])}
#' to the average of the x extent of the grid cell at indices \eqn{j}{j} and
#' \eqn{j+1}{j+1}. (The cell extents depend on the value of
#' `geographical`.)  The returned `x` and `y` values are the
#' mid-points of the supplied values. Thus, the returned `x` and `y`
#' are shorter than the supplied values by 1 item, and the returned `curl`
#' matrix dimensions are similarly reduced compared with the dimensions of
#' `u` and `v`.
#'
#' @param u matrix containing the 'x' component of a vector field
#'
#' @param v matrix containing the 'y' component of a vector field
#'
#' @param x the x values for the matrices, a vector of length equal to the
#' number of rows in `u` and `v`.
#'
#' @param y the y values for the matrices, a vector of length equal to the
#' number of cols in `u` and `v`.
#'
#' @param geographical logical value indicating whether `x` and `y`
#' are longitude and latitude, in which case spherical trigonometry is used.
#'
#' @param method A number indicating the method to be used to calculate the
#' first-difference approximations to the derivatives.  See \dQuote{Details}.
#'
#' @return A list containing vectors `x` and `y`, along with matrix
#' `curl`.  See \dQuote{Details} for the lengths and dimensions, for
#' various values of `method`.
#'
#' @section Development status.: This function is under active development as
#' of December 2014 and is unlikely to be stabilized until February 2015.
#'
#' @author Dan Kelley and Chantelle Layton
#'
#' @examples
#' library(oce)
#' # 1. Shear flow with uniform curl.
#' x <- 1:4
#' y <- 1:10
#' u <- outer(x, y, function(x, y) y/2)
#' v <- outer(x, y, function(x, y) -x/2)
#' C <- curl(u, v, x, y, FALSE)
#'
#' # 2. Rankine vortex: constant curl inside circle, zero outside
#' rankine <- function(x, y)
#' {
#'     r <- sqrt(x^2 + y^2)
#'     theta <- atan2(y, x)
#'     speed <- ifelse(r < 1, 0.5*r, 0.5/r)
#'     list(u=-speed*sin(theta), v=speed*cos(theta))
#' }
#' x <- seq(-2, 2, length.out=100)
#' y <- seq(-2, 2, length.out=50)
#' u <- outer(x, y, function(x, y) rankine(x, y)$u)
#' v <- outer(x, y, function(x, y) rankine(x, y)$v)
#' C <- curl(u, v, x, y, FALSE)
#' # plot results
#' par(mfrow=c(2, 2))
#' imagep(x, y, u, zlab="u", asp=1)
#' imagep(x, y, v, zlab="v", asp=1)
#' imagep(x, y, C$curl, zlab="curl", asp=1)
#' hist(C$curl, breaks=100)
#' @family things relating to vector calculus
curl <- function(u, v, x, y, geographical=FALSE, method=1)
{
    if (missing(u))
        stop("must supply u")
    if (missing(v))
        stop("must supply v")
    if (missing(x))
        stop("must supply x")
    if (missing(y))
        stop("must supply y")
    if (length(x) <= 1)
        stop("length(x) must exceed 1 but it is ", length(x))
    if (length(y) <= 1)
        stop("length(y) must exceed 1 but it is ", length(y))
    if (length(x) != nrow(u))
        stop("length(x) must equal nrow(u)")
    if (length(y) != ncol(u))
        stop("length(x) must equal ncol(u)")
    if (nrow(u) != nrow(v))
        stop("nrow(u) and nrow(v) must match")
    if (ncol(u) != ncol(v))
        stop("ncol(u) and ncol(v) must match")
    if (!is.logical(geographical))
        stop("geographical must be a logical quantity")
    method <- as.integer(round(method))
    if (1 == method) {
        res <- do_curl1(u, v, x, y, geographical)
    } else if (2 == method) {
        res <- do_curl2(u, v, x, y, geographical)
    } else {
        stop("method must be 1 or 2")
    }
    res
}


#' Calculate Range, Extended a Little, as is Done for Axes
#'
#' This is analogous to what is done as part of the R axis range calculation,
#' in the case where `xaxs="r"`.
#'
#' @param x a numeric vector.
#'
#' @param extend fraction to extend on either end
#'
#' @return A two-element vector with the extended range of `x`.
#'
#' @author Dan Kelley
rangeExtended <- function(x, extend=0.04) # extend by 4% on each end, like axes
{
    if (length(x) == 1) {
        x * c(1 - extend, 1 + extend)
    } else {
        r <- range(x, na.rm=TRUE)
        d <- diff(r)
        c(r[1] - d * extend, r[2] + d * extend)
    }
}

#' Extract (x, y, z) from (x, y, grid)
#'
#' Extract the grid points from a grid, returning columns.
#' This is useful for e.g. gridding large datasets, in which the first step
#' might be to use [binMean2D()], followed by
#' [interpBarnes()].
#'
#' @param x a vector holding the x coordinates of grid.
#'
#' @param y a vector holding the y coordinates of grid.
#'
#' @param grid a matrix holding the grid.
#'
#' @return A list containing three vectors: `x`, the grid x values,
#' `y`, the grid y values, and `grid`, the grid values.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' data(wind)
#' u <- interpBarnes(wind$x, wind$y, wind$z)
#' contour(u$xg, u$yg, u$zg)
#' U <- ungrid(u$xg, u$yg, u$zg)
#' points(U$x, U$y, col=oce.colorsViridis(100)[rescale(U$grid, rlow=1, rhigh=100)], pch=20)
ungrid <- function(x, y, grid)
{
    nrow <- nrow(grid)
    ncol <- ncol(grid)
    grid <- as.vector(grid) # by columns
    x <- rep(x, times=ncol)
    y <- rep(y, each=nrow)
    ok <- !is.na(grid)
    list(x=x[ok], y=y[ok], grid=grid[ok])
}


#' Draw error bars on an existing xy diagram
#'
#' @param x,y coordinates of points on the existing plot.
#'
#' @param xe,ye errors on x and y coordinates of points on the existing plot,
#' each either a single number or a vector of length identical to that of
#' the corresponding coordinate.
#'
#' @param percent boolean flag indicating whether `xe` and `ye` are
#' in terms of percent of the corresponding `x` and `y` values.
#'
#' @param style indication of the style of error bar.  Using `style=0`
#' yields simple line segments (drawn with [segments()]) and
#' `style=1` yields line segments with short perpendicular endcaps.
#'
#' @param length length of endcaps, for `style=1` only; it is passed to
#' [arrows()], which is used to draw that style of error bars.
#'
#' @param \dots graphical parameters passed to the code that produces the error
#' bars, e.g. to [segments()] for `style=0`.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' data(ctd)
#' S <- ctd[["salinity"]]
#' T <- ctd[["temperature"]]
#' plot(S, T)
#' errorbars(S, T, 0.05, 0.5)
errorbars <- function(x, y, xe, ye, percent=FALSE, style=0, length=0.025, ...)
{
    if (missing(x))
        stop("must supply x")
    if (missing(y))
        stop("must supply y")
    n <- length(x)
    if (n != length(y))
        stop("x and y must be of same length\n")
    if (missing(xe) && missing(ye))
        stop("must give either xe or ye")
    if (1 == length(xe))
        xe <- rep(xe, n) # FIXME probably gives wrong length
    if (1 == length(ye))
        ye <- rep(ye, n)
    if (!missing(xe)) {
        if (n != length(xe))
            stop("x and xe must be of same length\n")
        if (percent)
            xe <- xe * x / 100
        look <- xe != 0
        if (style == 0) {
            segments(x[look], y[look], x[look]+xe[look], y[look], ...)
            segments(x[look], y[look], x[look]-xe[look], y[look], ...)
        } else if (style == 1) {
            arrows(x[look], y[look], x[look] + xe[look], y[look], angle=90, length=length, ...)
            arrows(x[look], y[look], x[look] - xe[look], y[look], angle=90, length=length, ...)
        } else {
            stop("unknown value ", style, " of style; must be 0 or 1\n")
        }
    }
    if (!missing(ye)) {
        if (n != length(ye))
            stop("y and ye must be of same length\n")
        if (percent)
            ye <- ye * y / 100
        look <- ye != 0
        if (style == 0) {
            segments(x[look], y[look], x[look], y[look]+ye[look], ...)
            segments(x[look], y[look], x[look], y[look]-ye[look], ...)
        } else if (style == 1) {
            arrows(x[look], y[look], x[look], y[look] + ye[look], angle=90, length=length, ...)
            arrows(x[look], y[look], x[look], y[look] - ye[look], angle=90, length=length, ...)
        } else {
            stop("unknown value ", style, " of style; must be 0 or 1\n")
        }
    }
}

filterSomething <- function(x, filter)
{
    if (is.raw(x)) {
        x <- as.numeric(x)
        replace <- mean(x, na.rm=TRUE)
        x[is.na(x)] <- replace
        res <- as.integer(filter(x, filter))
        res <- ifelse(res < 0, 0, res)
        res <- ifelse(res > 255, 255, res)
        res <- as.raw(res)
    } else {
        replace <- mean(x, na.rm=TRUE)
        x[is.na(x)] <- replace
        res <- filter(x, filter)
    }
    res
}


#' Plot a Model-data Comparison Diagram
#'
#' Creates a diagram as described by Taylor (2001).  The graph is in the form
#' of a semicircle, with radial lines and spokes connecting at a focus point on
#' the flat (lower) edge.  The radius of a point on the graph indicates the
#' standard deviation of the corresponding quantity, i.e. x and the columns in
#' y.  The angle connecting a point on the graph to the focus provides an
#' indication of correlation coefficient with respect to x.  The ``east'' side
#' of the graph indicates \eqn{R=1}{R=1}, while \eqn{R=0}{R=0} is at the
#' "north" edge and \eqn{R=-1}{R=-1} is at the "west" side.  The `x`
#' data are indicated with a bullet on the graph, appearing on the lower edge
#' to the right of the focus at a distance indicating the standard deviation of
#' `x`.  The other points on the graph represent the columns of `y`,
#' coded automatically or with the supplied values of `pch` and
#' `col`.
#' The example shows two tidal models of the Halifax sealevel data, computed
#' with [tidem()] with just the M2 component and the S2 component;
#' the graph indicates that the M2 model is much better than the S2 model.
#'
#' @param x a vector of reference values of some quantity, e.g. measured over
#' time or space.
#'
#' @param y a matrix whose columns hold values of values to be compared with
#' those in x.  (If `y` is a vector, it is converted first to a one-column
#' matrix).
#'
#' @param scale optional scale, interpreted as the maximum value of standard
#' deviation.
#'
#' @param pch list of characters to plot, one for each column of `y`.
#'
#' @param col list of colors for points on the plot, one for each column of
#' `y`.
#'
#' @param labels optional vector of strings to use for labelling the points.
#'
#' @param pos optional vector of positions for labelling strings.  If not
#' provided, labels will be to the left of the symbols.
#'
#' @param \dots optional arguments passed by `plotTaylor` to more child
#' functions.
#'
#' @author Dan Kelley
#'
#' @references Taylor, Karl E., 2001.  Summarizing multiple aspects of model
#' performance in a single diagram, *J. Geophys. Res.*, 106:D7, 7183--7192.
#'
#' @examples
#' library(oce)
#' data(sealevel)
#' x <- sealevel[["elevation"]]
#' M2 <- predict(tidem(sealevel, constituents="M2"))
#' S2 <- predict(tidem(sealevel, constituents=c("S2")))
#' plotTaylor(x, cbind(M2, S2))
plotTaylor <- function(x, y, scale, pch, col, labels, pos, ...)
{
    if (missing(x))
        stop("must supply 'x'")
    if (missing(y))
        stop("must supply 'y'")
    if (is.vector(y))
        y <- matrix(y)
    ncol <- ncol(y)
    if (missing(pch))
        pch <- 1:ncol
    if (missing(col))
        col <- rep("black", ncol)
    haveLabels <- !missing(labels)
    if (missing(pos))
        pos <- rep(2, ncol)
    if (length(pos) < ncol)
        pos <- rep(pos[1], ncol)
    xSD <- sd(x, na.rm=TRUE)
    ySD <- sd(as.vector(y), na.rm=TRUE)
    if (missing(y))
        stop("must supply 'y'")
    halfArc <- seq(0, pi, length.out=200)
    # FIXME: use figure geometry, to avoid axis cutoff
    if (missing(scale))
        scale <- max(pretty(c(xSD, ySD)))
    plot.new()
    plot.window(c(-1.2, 1.2) * scale, c(0, 1.2) * scale, asp=1)
    #plot.window(c(-1.1, 1.1), c(0.1, 1.2), asp=1)
    sdPretty <- pretty(c(0, scale))
    for (radius in sdPretty)
        lines(radius * cos(halfArc), radius * sin(halfArc), col="gray")
    # spokes
    for (rr in seq(-1, 1, 0.2)) {
        lines(c(0, max(sdPretty)*cos(pi/2 + rr * pi / 2)),
            c(0, max(sdPretty)*sin(pi/2 + rr * pi / 2)), col="gray")
    }
    axisLabels <- format(sdPretty)
    axisLabels[1] <- paste(0)
    axis(1, pos=0, at=sdPretty, labels=axisLabels)
    # temporarily permit labels outside the platting zone
    xpdOld <- par("xpd")
    par(xpd=NA)
    m <- max(sdPretty)
    text(m, 0, "R=1", pos=4)
    text(0, m, "R=0", pos=3)
    text(-m, 0, "R=-1", pos=2)
    par(xpd=xpdOld)
    points(xSD, 0, pch=20, cex=1.5)
    for (column in seq_len(ncol(y))) {
        ySD <- sd(y[, column], na.rm=TRUE)
        R <- cor(x, y[, column])^2
        #cat("column=", column, "ySD=", ySD, "R=", R, "col=", col[column], "pch=", pch[column], "\n")
        xx <- ySD * cos((1 - R) * pi / 2)
        yy <- ySD * sin((1 - R) * pi / 2)
        points(xx, yy, pch=pch[column], lwd=2, col=col[column], cex=2)
        if (haveLabels) {
            #cat(labels[column], "at", pos[column], "\n")
            text(xx, yy, labels[column], pos=pos[column], ...)
        }
    }
}


#' Pretty lat/lon in deg, min, sec
#'
#' Round a geographical positions in degrees, minutes, and seconds
#' Depending on the range of values in `x`, rounding is done to degrees,
#' half-degrees, minutes, etc.
#'
#' @param x a series of one or more values of a latitude or longitude, in
#' decimal degrees
#'
#' @param debug set to a positive value to get debugging information during
#' processing.
#'
#' @return A vector of numbers that will yield good axis labels if
#' [formatPosition()] is used.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' formatPosition(prettyPosition(10+1:10/60+2.8/3600))
prettyPosition <- function(x, debug=getOption("oceDebug"))
{
    oceDebug(debug, "prettyPosition(...) {\n", sep="", unindent=1)
    r <- diff(range(x, na.rm=TRUE))
    oceDebug(debug, "range(x)=", range(x), ", r=", r, "\n")
    if (r > 5) {
        # D only
        res <- pretty(x)
    } else if (r > 1) {
        # round to 30 minutes
        res <- (1 / 2) * pretty(2 * x)
        oceDebug(debug, "case 1: res=", res, "\n")
    } else if (r > 30/60) {
        # round to 1 minute, with extras
        res <- (1 / 60) * pretty(60 * x, n=6)
        oceDebug("case 2: res=", res, "\n")
    } else if (r > 5/60) {
        # round to 1 minute
        res <- (1 / 60) * pretty(60 * x, 4)
        oceDebug(debug, "case 3: res=", res, "\n")
    } else if (r > 10/3600) {
        # round to 10 sec
        res <- (1 / 360) * pretty(360 * x)
        oceDebug(debug, "case 4: res=", res, "\n")
    } else {
        # round to seconds
        res <- (1 / 3600) * pretty(3600 * x)
        if (debug) cat("case 5: res=", res, "\n")
    }
    oceDebug(debug, "} # prettyPosition\n", unindent=1)
    res
}

smoothSomething <- function(x, ...)
{
    if (is.raw(x)) {
        x <- as.numeric(x)
        replace <- mean(x, na.rm=TRUE)
        x[is.na(x)] <- replace
        res <- as.integer(smooth(x, ...))
        res <- ifelse(res < 0, 0, res)
        res <- ifelse(res > 255, 255, res)
        res <- as.raw(res)
    } else {
        replace <- mean(x, na.rm=TRUE)
        x[is.na(x)] <- replace
        res <- smooth(x, ...)
    }
    res
}


#' Rescale values to lie in a given range
#'
#' This is helpful in e.g. developing a color scale for an image plot.  It is
#' not necessary that `rlow` be less than `rhigh`, and in fact
#' reversing them is a good way to get a reversed color scale for a plot.
#'
#' @param x a numeric vector.
#'
#' @param xlow `x` value to correspond to `rlow`.  If not given, it
#' will be calculated as the minimum value of `x`
#'
#' @param xhigh `x` value to correspond to `rhigh`.  If not given, it
#' will be calculated as the maximum value of `x`
#'
#' @param rlow value of the result corresponding to `x` equal to
#' `xlow`.
#'
#' @param rhigh value of the result corresponding to `x` equal to
#' `xhigh`.
#'
#' @param clip logical, set to `TRUE` to clip the result to the range
#' spanned by `rlow` and `rhigh`.
#'
#' @return A new vector, which has minimum `lim[1]` and maximum `lim[2]`.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' # Fake tow-yow data
#' t <- seq(0, 600, 5)
#' x <- 0.5 * t
#' z <- 50 * (-1 + sin(2 * pi * t / 360))
#' T <- 5 + 10 * exp(z / 100)
#' palette <- oce.colorsViridis(100)
#' zlim <- range(T)
#' drawPalette(zlim=zlim, col=palette)
#' plot(x, z, type="p", pch=20, cex=3,
#'      col=palette[rescale(T, xlow=zlim[1], xhigh=zlim[2], rlow=1, rhigh=100)])
rescale <- function(x, xlow, xhigh, rlow=0, rhigh=1, clip=TRUE)
{
    x <- as.numeric(x)
    finite <- is.finite(x)
    #r <- range(x, na.rm=TRUE)
    if (missing(xlow))
        xlow <- min(x, na.rm=TRUE)
    if (missing(xhigh))
        xhigh <- max(x, na.rm=TRUE)
    res <- rlow + (rhigh - rlow) * (x - xlow) / (xhigh - xlow)
    if (clip) {
        res <- ifelse(res < min(rlow, rhigh), rlow, res)
        res <- ifelse(res > max(rlow, rhigh), rhigh, res)
    }
    res[!finite] <- NA
    res
}


#' Adjust the time within Oce object
#'
#' This function compensates for drifting instrument clocks, according to
#' \eqn{t'=t + a + b (t-t0)}{t'=t + a + b*(t-t0)}, where \eqn{t'}{t'} is the
#' true time and \eqn{t}{t} is the time stored in the object.  A single check
#' on time mismatch can be described by a simple time offset, with a non-zero
#' value of `a`, a zero value of `b`, and an arbitrary value of
#' `t0`.  Checking the mismatch before and after an experiment yields
#' sufficient information to specify a linear drift, via `a`, `b`,
#' and `t0`.  Note that `t0` is just a convenience parameter, which
#' avoids the user having to know the "zero time" used in R and clarifies the
#' values of the other two parameters.  It makes sense for `t0` to have
#' the same timezone as the time within `x`.
#'
#' The returned object is computed by linear interpolation, using
#' [approx()] with `rule=2`, to avoid `NA` values at the
#' start or end.  The new time will be as given by the formula above. Note that
#' if the drift is large enough, the sampling rate will be changed.  It is a
#' good idea to start with an object that has an extended time range, so that,
#' after this is called, [subset()] can be used to trim to a desired
#' time range.
#'
#' @param x an [oce-class] object.
#'
#' @param a intercept (in seconds) in linear model of time drift (see
#' \dQuote{Details}).
#'
#' @param b slope (unitless) in linear model of time drift (unitless) (see
#' \dQuote{Details}).
#'
#' @param t0 reference time (in [POSIXct()] format) used in linear model of time
#' drift (see \dQuote{Details}).
#'
#' @param debug a flag that, if nonzero, turns on debugging.
#'
#' @return A new object, with time and other data adjusted.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' data(adv)
#' adv2 <- retime(adv,0,1e-4,as.POSIXct("2008-07-01 00:00:00", tz="UTC"))
#' plot(adv[["time"]], adv2[["time"]]-adv[["time"]], type="l")
retime <- function(x, a, b, t0, debug=getOption("oceDebug"))
{
    if (missing(x))
        stop("must give argument 'x'")
    if (missing(a))
        stop("must give argument 'a'")
    if (missing(b))
        stop("must give argument 'b'")
    if (missing(t0))
        stop("must give argument 't0'")
    oceDebug(debug, paste("retime.adv(x, a=", a, ", b=", b, ", t0=\"", format(t0), "\")\n"), sep="", unindent=1)
    res <- x
    oceDebug(debug, "retiming x@data$time")
    res@data$time <- x@data$time + a + b * (as.numeric(x@data$time) - as.numeric(t0))
    if ("timeSlow" %in% names(x@data)) {
        oceDebug(debug, "retiming x@data$timeSlow\n")
        res@data$timeSlow <- x@data$timeSlow + a + b * (as.numeric(x@data$timeSlow) - as.numeric(t0))
    }
    res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep="", collapse=""))
    oceDebug(debug, "} # retime.adv()\n", unindent=1)
    res
}


#' Calculate min, mean, and max values
#'
#' This is a simpler cousin of the standard [fivenum()] function,
#' used in [summary()] functions for `oce` objects.
#'
#' @param x a vector or matrix of numerical values.
#'
#' @return A character vector of three values: the minimum, the mean, the
#' maximum.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' threenum(1:10)
#' @section Historical note:
#' On Aug 5, 2019, the dimension was dropped as the fourth column, and
#' this function returned to the original intention (revealed by its
#' name).  Another change is that the function now returns numerical
#' results, leaving the task of setting the number of digits to
#' [summary()].
threenum <- function(x)
{
    dim <- dim(x)
    if (is.character(x) || is.null(x)) {
        res <- rep(NA, 3)
    } else if (is.list(x)) {
        if (2 == sum(c("lsb", "msb") %in% names(x))) {
            # e.g. landsat data
            x <- as.numeric(x$lsb) + 256 * as.numeric(x$msb)
            dim(x) <- dim
            res <- c(min(x, na.rm=TRUE), mean(x, na.rm=TRUE), max(x, na.rm=TRUE))
        } else {
            res <- rep(NA, 3)
        }
    } else if (is.raw(x)) {
        x <- as.numeric(x)
        dim(x) <- dim
        res <- c(min(x, na.rm=TRUE), mean(x, na.rm=TRUE), max(x, na.rm=TRUE))
    } else if (is.factor(x)) {
        res <- rep(NA, 3)
    } else if (0 < sum(!is.na(x))) {
        res <- c(min(x, na.rm=TRUE), mean(x, na.rm=TRUE), max(x, na.rm=TRUE))
    } else if (inherits(x, "POSIXt")) {
        res <- c(min(x, na.rm=TRUE), mean(x, na.rm=TRUE), max(x, na.rm=TRUE))
    } else {
        res <- rep(NA, 3)
    }
    res
}

normalize <- function(x)
{
    var <- var(x, na.rm=TRUE)
    if (var == 0) {
        rep(0, length(x))
    } else {
        (x - mean(x, na.rm=TRUE)) / sqrt(var)
    }
}


#' Detrend a set of observations
#'
#' Detrends `y` by subtracting a linear trend in `x`, to create
#' a vector that is zero for its first and last finite value.
#' If the second parameter (`y`) is missing, then `x` is
#' taken to be `y`, and a new `x` is constructed with
#' [seq_along()].  Any `NA` values are left as-is.
#'
#' A common application is to bring the end points of a time series
#' down to zero, prior to applying a digital filter. (See examples.)
#'
#' @param x a vector of numerical values.  If `y` is not given, then
#' `x` is taken for `y`.
#'
#' @param y an optional vector
#'
#' @return A list containing `Y`, the detrended version of `y`, and
#' the intercept `a` and slope `b` of the linear function of `x`
#' that is subtracted from `y` to yield `Y`.
#'
#' @author Dan Kelley
#'
#' @examples
#' x <- seq(0, 0.9 * pi, length.out=50)
#' y <- sin(x)
#' y[1] <- NA
#' y[10] <- NA
#' plot(x, y, ylim=c(0, 1))
#' d <- detrend(x, y)
#' points(x, d$Y, pch=20)
#' abline(d$a, d$b, col="blue")
#' abline(h=0)
#' points(x, d$Y + d$a + d$b * x, col="blue", pch="+")
detrend <- function(x, y)
{
    if (missing(x))
        stop("must give x")
    n <- length(x)
    if (missing(y)) {
        y <- x
        x <- seq_along(y)
    } else {
        if (length(y) != n)
            stop("x and y must be of same length, but they are ", n, " and ", length(y))
    }
    first <- which(is.finite(y))[1]
    last <- 1 + length(y) - which(is.finite(rev(y)))[1]
    if (x[first] == x[last])
        stop("the first and last x values must be distinct")
    b <- (y[first] - y[[last]]) / (x[first] - x[[last]])
    a <- y[first] - b * x[first]
    list(Y=y - (a+b*x), a=a, b=b)
}



#' Remove spikes from a time series
#'
#' The method identifies spikes with respect to a "reference" time-series, and
#' replaces these spikes with the reference value, or with `NA` according
#' to the value of `action`; see \dQuote{Details}.
#'
#' Three modes of operation are permitted, depending on the value of
#' `reference`.
#'
#' 1. For `reference="median"`, the first step is to linearly interpolate
#' across any gaps (spots where `x==NA`), using [approx()] with
#' `rule=2`. The second step is to pass this through
#' [runmed()] to get a running median spanning `k`
#' elements. The result of these two steps is the "reference" time-series.
#' Then, the standard deviation of the difference between `x`
#' and the reference is calculated.  Any `x` values that differ from
#' the reference by more than `n` times this standard deviation are considered
#' to be spikes.  If `replace="reference"`, the spike values are
#' replaced with the reference, and the resultant time series is
#' returned.  If `replace="NA"`, the spikes are replaced with `NA`,
#' and that result is returned.
#'
#' 2. For `reference="smooth"`, the processing is the same as for
#' `"median"`, except that [smooth()] is used to calculate the
#' reference time series.
#'
#' 3. For `reference="trim"`, the reference time series is constructed by
#' linear interpolation across any regions in which `x<min` or
#' `x>max`.  (Again, this is done with [approx()] with
#' `rule=2`.) In this case, the value of `n` is ignored, and the
#' return value is the same as `x`, except that spikes are replaced
#' with the reference series (if `replace="reference"` or with
#' `NA`, if `replace="NA"`.
#'
#' @param x a vector of (time-series) values, a list of vectors, a data frame,
#' or an [oce-class] object.
#'
#' @param reference indication of the type of reference time series to be used
#' in the detection of spikes; see \dQuote{Details}.
#'
#' @param n an indication of the limit to differences between `x` and the
#' reference time series, used for `reference="median"` or
#' `reference="smooth"`; see \dQuote{Details.}
#'
#' @param k length of running median used with `reference="median"`, and
#' ignored for other values of `reference`.
#'
#' @param min minimum non-spike value of `x`, used with
#' `reference="trim"`.
#'
#' @param max maximum non-spike value of `x`, used with
#' `reference="trim"`.
#'
#' @param replace an indication of what to do with spike values, with
#' `"reference"` indicating to replace them with the reference time
#' series, and `"NA"` indicating to replace them with `NA`.
#'
#' @param skip optional vector naming columns to be skipped. This is ignored if
#' `x` is a simple vector. Any items named in `skip` will be passed
#' through to the return value without modification.  In some cases,
#' `despike` will set up reasonable defaults for `skip`, e.g. for a
#' `ctd` object, `skip` will be set to \code{c("time", "scan",
#' "pressure")} if it is not supplied as an argument.
#'
#' @return A new vector in which spikes are replaced as described above.
#'
#' @author Dan Kelley
#'
#' @examples
#' n <- 50
#' x <- 1:n
#' y <- rnorm(n=n)
#' y[n/2] <- 10                    # 10 standard deviations
#' plot(x, y, type="l")
#' lines(x, despike(y), col="red")
#' lines(x, despike(y, reference="smooth"), col="darkgreen")
#' lines(x, despike(y, reference="trim", min=-3, max=3), col="blue")
#' legend("topright", lwd=1, col=c("black", "red", "darkgreen", "blue"),
#'        legend=c("raw", "median", "smooth", "trim"))
#'
#' # add a spike to a CTD object
#' data(ctd)
#' plot(ctd)
#' T <- ctd[["temperature"]]
#' T[10] <- T[10] + 10
#' ctd[["temperature"]] <- T
#' CTD <- despike(ctd)
#' plot(CTD)
despike <- function(x, reference=c("median", "smooth", "trim"), n=4, k=7, min=NA, max=NA,
    replace=c("reference", "NA"), skip)
{
    if (is.vector(x)) {
        x <- despikeColumn(x, reference=reference, n=n, k=k, min=min, max=max, replace=replace)
    } else {
        if (missing(skip)) {
            if (inherits(x, "ctd")) {
                skip <- c("time", "scan", "pressure")
            } else {
                skip <- NULL
            }
        }
        if (inherits(x, "oce")) {
            columns <- names(x@data)
            for (column in columns) {
                if (!(column %in% skip)) {
                    # check for NA column
                    if (all(is.na(x[[column]]))) {
                        warning(paste("Column", column, "contains only NAs. Skipping"))
                    } else {
                        x[[column]] <- despikeColumn(x[[column]],
                            reference=reference, n=n, k=k, min=min, max=max, replace=replace)
                    }
                }
            }
            x@processingLog <- processingLogAppend(x@processingLog, paste(deparse(match.call()), sep="", collapse=""))
        } else {
            columns <- names(x)
            for (column in columns) {
                if (!(column %in% skip)) {
                    if (all(is.na(x[[column]]))) {
                        warning(paste("Column", column, "contains only NAs. Skipping"))
                    } else {
                        x[[column]] <- despikeColumn(x[[column]],
                            reference=reference, n=n, k=k, min=min, max=max, replace=replace)
                    }
                }
            }
        }
    }
    x
}

despikeColumn <- function(x, reference=c("median", "smooth", "trim"), n=4, k=7, min=NA, max=NA,
    replace=c("reference", "NA"))
{
    reference <- match.arg(reference)
    replace <- match.arg(replace)
    gave.min <- !is.na(min)
    gave.max <- !is.na(max)
    nx <- length(x)
    # degap
    na <- is.na(x)
    if (sum(na) > 0) {
        i <- 1:nx
        x.gapless <- approx(i[!na], x[!na], i, rule=2)$y
    } else {
        x.gapless <- x
    }
    if (reference == "median" || reference == "smooth") {
        if (reference == "median") {
            x.reference <- runmed(x.gapless, k=k)
        } else {
            x.reference <- as.numeric(smooth(x.gapless))
        }
        distance <- abs(x.reference - x.gapless)
        stddev <- sqrt(var(distance))
        bad <- distance > n * stddev
        nbad <- sum(bad)
        if (nbad > 0) {
            if (replace == "reference") {
                x[bad] <- x.reference[bad]
            } else {
                x[bad] <- rep(NA, nbad)
            }
        }
    } else if (reference == "trim") {
        if (!gave.min || !gave.max)
            stop("must give min and max")
        bad <- !(min <= x & x <= max)
        nbad <- length(bad)
        if (nbad > 0) {
            i <- 1:nx
            if (replace == "reference") {
                x[bad] <- approx(i[!bad], x.gapless[!bad], i[bad], rule=2)$y
            } else {
                x[bad] <- NA
            }
        }
    } else {
        stop("unknown reference ", reference)
    }
    x
}



#' Substitute NA for data outside a range
#'
#' Substitute NA for data outside a range, e.g. to remove wild spikes in data.
#'
#' @param x vector of values
#'
#' @param min minimum acceptable value.  If not supplied, and if `max` is
#' also not supplied, a `min` of the 0.5 percentile will be used.
#'
#' @param max maximum acceptable value.  If not supplied, and if `min` is
#' also not supplied, a `min` of the 0.995 percentile will be used.
#'
#' @author Dan Kelley
#'
#' @examples
#'
#' ten.to.twenty <- rangeLimit(1:100, 10, 20)
rangeLimit <- function(x, min, max)
{
    if (missing(min) && missing(max)) {
        minmax <- quantile(x, 0.005, 0.995)
        min <- minmax[1]
        max <- minmax[2]
    }
    ifelse(max < x | x < min, NA, x)
}


#' Determine year from various abbreviations
#'
#' Various data files may contain various abbreviations for years.  For
#' example, 99 refers to 1999, and 8 refers to 2008.  Sometimes, even 108
#' refers to 2008 (the idea being that the "zero" year was 1900).  This
#' function deals with the three cases mentioned.  It will fail if someone
#' supplies 60, meaning year 2060 as opposed to 1960.
#'
#' @param year a year, or vector of years, possibly abbreviated
#'
#' @author Dan Kelley
#'
#' @examples
#' fullYear <- unabbreviateYear(c(99, 8, 108))
#' @family things related to time
unabbreviateYear <- function(year)
{
    # handle e.g. 2008 as 2008 (full year), 8 (year-2000 offset), or 108 (year 1900 offset)
    ifelse(year > 1800, year, ifelse(year > 50, year + 1900, year + 2000))
}



#' Unwrap an angle that suffers modulo-360 problems
#'
#' This is mostly used for instrument heading angles, in cases where the
#' instrument is aligned nearly northward, so that small variations in heading
#' (e.g. due to mooring motion) can yield values that swing from small angles
#' to large angles, because of the modulo-360 cut point.
#' The method is to use the cosine and sine of the angle, to construct "x" and
#' "y" values on a unit circle, then to find means and medians of x and y
#' respectively, and finally to use [atan2()] to infer the angles.
#'
#' @param angle an angle (in degrees) that is thought be near 360 degrees, with
#' added noise
#'
#' @return A list with two estimates: `mean` is based on an arithmetic
#' mean, and `median` is based on the median. Both are mapped to the range
#' 0 to 360.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' true <- 355
#' a <- true + rnorm(100, sd=10)
#' a <- ifelse(a > 360, a - 360, a)
#' a2 <- unwrapAngle(a)
#' par(mar=c(3, 3, 5, 3))
#' hist(a, breaks=360)
#' abline(v=a2$mean, col="blue", lty="dashed")
#' abline(v=true, col="blue")
#' mtext("true (solid)\n estimate (dashed)", at=true, side=3, col="blue")
#' abline(v=mean(a), col="red")
#' mtext("mean", at=mean(a), side=3, col="red")
unwrapAngle <- function(angle)
{
    toRad <- atan2(1, 1) / 45
    angle <- angle * toRad
    S <- sin(angle)
    C <- cos(angle)
    Smean <- mean(S, na.rm=TRUE)
    Smedian <- median(S, na.rm=TRUE)
    Cmean <- mean(C, na.rm=TRUE)
    Cmedian <- median(C, na.rm=TRUE)
    resMean <- atan2(Smean, Cmean)/toRad
    resMedian <- atan2(Smedian, Cmedian)/toRad
    resMean <- if (resMean < 0) resMean + 360 else resMean
    resMedian <- if (resMedian < 0) resMedian + 360 else resMedian
    list(mean=resMean, median=resMedian)
}


#' 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


#' Wrapper to give normalized 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


#' Show some values from a list, vector or matrix
#'
#' This is similar to [str()], but it shows data at the first and last of the
#' vector, which can be quite helpful in debugging.
#'
#' @param v the item to be summarized.  If this is a list of a vector of named
#' items, then information is provided for each element. Otherwise, information
#' is provided for the first and last `n` values.
#'
#' @param msg optional character value indicating a message to show,
#' introducing the vector.  If not provided, then
#' a message is created from `v`. If `msg` is a non-empty string,
#' then that string is pasted together with a colon (unless `msg` already
#' contains a colon), before pasting a summary of data values.
#'
#' @param postscript optional character value indicating an optional message
#' to append at the end of the return value.
#'
#' @param digits for numerical values of `v`, this is the number of digits
#' to use, in formatting the numbers with [format()]; otherwise,
#' `digits` is ignored.
#'
#' @param n number of elements to show at start and end. If `n`
#' is negative, then all the elements are shown.
#'
#' @param showNA logical value indicating whether to show the number
#' of `NA` values. This is done only if the output contains ellipses,
#' meaning that some values are skipped, because if all values are shown,
#' it will be perfectly obvious whether there are any `NA` values.
#'
#' @param showNewline logical value indicating whether to put a newline
#' character at the end of the output string.  The default, TRUE, is
#' convenient for printing, but using FALSE makes more sense if
#' the result is to be used with, e.g. [mtext()].
#'
#' @return A string ending in a newline character, suitable for
#' display with [cat()] or [oceDebug()].
#'
#' @examples
#' # List
#' limits <- list(low=0, high=1)
#' vectorShow(limits)
#'
#' # Vector of named items
#' planktonCount <- c(phytoplankton=100, zooplankton=20)
#' vectorShow(planktonCount)
#'
#' # Vector
#' vectorShow(pi)
#'
#' # Matrix
#' vectorShow(volcano)
#'
#' # Other arguments
#' knot2mps <- 0.5144444
#' vectorShow(knot2mps, postscript="knots per m/s")
#' vectorShow("January", msg="The first month is")
#'
#' @author Dan Kelley
vectorShow <- function(v, msg="", postscript="", digits=5L, n=2L, showNA=FALSE, showNewline=TRUE)
{
    startEnd <- function(v, n)
    {
        if (length(v) < 2L*n) {
            paste(v, collapse=", ")
        } else {
            paste0(paste(head(v, n), collapse=","), ",...,",
                paste(paste(tail(v, n)), collapse=","))
        }
    }
    dimv <- dim(v)
    nv <- length(v)
    if (!nchar(msg)) {
        msg <- deparse(substitute(expr=v, env=environment()))
    }
    if (is.list(v) || (is.vector(v) && length(names(v)) > 0L)) {
        names <- names(v)
        values <- unname(v)
        msg <- paste0(msg, ": ")
        nv <- length(v)
        for (iv in seq_len(nv)) {
            msg <- paste0(msg, names[iv], "=", startEnd(values[[iv]], n))
            if (iv < nv) {
                msg <- paste0(msg, ", ")
            }
        }
        if (showNewline) {
            msg <- paste0(msg, "\n")
        }
        return(msg)
    } else {
        if (!is.null(dimv)) {
            msg <- paste0(msg,
                paste0("[",
                    paste(unlist(lapply(dimv, function(x) paste0("1:", x))), collapse=", "),
                    "]"))
        } else if (nv > 1) {
            msg <- paste0(msg, paste0("[1:", nv, "]"))
        }
        if (nchar(msg) && !grepl(":$", msg)) {
            msg <- paste0(msg, ": ")
        }
    }
    res <- msg
    if (nv == 0) {
        res <- paste(res, "(empty vector)")
    } else {
        if (n < 0 || nv <= 2*n) {
            showAll <- TRUE
        } else {
            n <- floor(min(n, nv/2))
            showAll <- FALSE
        }
        if (is.numeric(v)) {
            if (showAll) {
                res <- paste(msg, paste(format(v, digits=digits), collapse=", "), sep="")
            } else {
                res <- paste(msg, paste(format(v[1:n], digits=digits), collapse=", "),
                    ", ..., ", paste(format(v[nv-seq.int(n-1, 0)], digits=digits), collapse=", "),
                    sep="")
                if (showNA)
                    res <- paste0(res, " (", sum(is.na(v)), " NA)")
            }
        } else {
            if (showAll) {
                res <- paste(msg, paste(v, collapse=", "), sep="")
            } else {
                res <- paste(msg, paste(v[1:n], collapse=", "),
                             ", ..., ", paste(v[nv-seq.int(n-1, 0)], collapse=", "), sep="")
                if (showNA)
                    res <- paste0(res, " (", sum(is.na(v)), " NA)")
            }
        }
    }
    if (nchar(postscript) > 0) {
        res <- paste(res, postscript)
    }
    if (showNewline) {
        res <- paste(res, "\n", sep="")
    }
    res
}


#' Full name of file, including path
#'
#' Determines the full name of a file, including the path.  Used by many
#' `read.X` routines, where `X` is the name of a class of object.
#' This is a wrapper around [normalizePath()], with warnings turned
#' off so that messages are not printed for unfound files (e.g. URLs).
#'
#' @param filename name of file
#'
#' @return Full file name
#'
#' @author Dan Kelley
fullFilename <- function(filename)
{
    warn <- options("warn")$warn
    options(warn=-1)
    res <- normalizePath(filename)
    options(warn=warn)
    res
}


#' Provide axis names in adjustable sizes
#'
#' Provide axis names in adjustable sizes, e.g. using T instead of Temperature,
#' and including units as appropriate.
#' Used by e.g. [plot,ctd-method()].
#'
#' @param item code for the label. The following common values are recognized:
#' `"absolute salinity"`, `"along-spine distance km"`, `"along-track distance km"`,
#' `"C"`, `"conductivity mS/cm"`, `"conductivity S/m"`, `"conservative temperature"`,
#' `"CT"`, `"depth"`, `"direction"`, `"distance"`, `"distance km"`, `"eastward"`,
#' `"elevation"`, `"fluorescence"`, `"frequency cph"`, `"heading"`, `"latitude"`,
#' `"longitude"`, `"N2"`, `"nitrate"`, `"nitrite"`, `"northward"`, `"oxygen"`,
#' `"oxygen mL/L"`, `"oxygen saturation"`, `"oxygen umol/kg"`, `"oxygen umol/L"`,
#' `"p"`, `"phosphate"`, `"pitch"`, `"roll"`, `"S"`, `"SA"`,
#' `"sigma0"`, `"sigma1"`, `"sigma2"`, `"sigma3"`, `"sigma4"`,
#' `"sigmaTheta"`,
#' `"silicate"`, `"sound speed"`, `"spectral density m2/cph"`, `"speed"`,
#' `"spice"`, `"T"`, `"theta"`, `"tritium"`, `"u"`, `"v"`, `"w"`, or `"z"`.
#' Other values may also be recognized, and if an unrecognized item is
#' given, then it is returned, unaltered.
#'
#' @param axis a string indicating which axis to use; must be `x` or
#' `y`.
#'
#' @param sep optional character string inserted between the unit and the
#' parentheses or brackets that enclose it. If not provided, then
#' [`getOption`]`("oceUnitSep")` is checked. If that exists, then it is
#' used as the separator; if not, no separator is used.
#'
#' @param unit optional unit to use, if the default is not satisfactory. This
#' might be the case if for example temperature was not measured in Celcius.
#'
#' @param debug optional debugging flag. Setting to 0 turns off debugging,
#' while setting to 1 causes some debugging information to be printed.
#'
#' @return A character string or expression, in either a long or a shorter
#' format, for use in the indicated axis at the present plot size.  Whether the
#' unit is enclosed in parentheses or square brackets is determined by the
#' value of `getOption("oceUnitBracket")`, which may be `"["` or
#' `"("`.  Whether spaces are used between the unit and these deliminators
#' is set by `psep` or [`getOption`]`("oceUnitSep")`.
#'
#' @family functions that create labels
#'
#' @author Dan Kelley
resizableLabel <- function(item, axis="x", sep, unit=NULL, debug=getOption("oceDebug"))
{
    oceDebug(debug, "resizableLabel(item=\"", item,
        "\", axis=\"", axis,
        "\", sep=\"", if (missing(sep)) "(missing)" else sep, "\", ...) {\n",
        sep="", unindent=1, style="bold")
    if (missing(item))
        stop("must provide 'item'")
    if (axis != "x" && axis != "y")
        stop("axis must be \"x\" or \"y\"")
    itemAllowed <- c("salinity", "S", "SA", "C", "CT", paste("conductivity", "mS/cm"),
        paste("conductivity", "S/m"), "T", "temperature", "theta", "sigmaTheta",
        paste("Conservative", "Temperature"), paste("Absolute", "Salinity"),
        "N2", "nitrate", "nitrite", "oxygen", paste("oxygen", "saturation"),
        paste("oxygen", "mL/L"), paste("oxygen", "umol/L"), paste("oxygen",
            "umol/kg"), "phosphate", "silicate", "tritium", "spice",
        "fluorescence", "p", "z", "distance", "distance km",
        paste("along-spine", "distance", "km"), paste("along-track", "distance",
            "km"), "heading", "pitch", "roll", "u", "v", "w", "speed",
        "direction", "eastward", "northward", "depth", "elevation", "latitude",
        "longitude", paste("frequency", "cph"), paste("sound", "speed"),
        paste("spectral", "density", "m2/cph"), "sigma0", "sigma1", "sigma2",
        "sigma3", "sigma4", "Sstar", "SR")

    # FIXME: if anything is added, run the next, and paste results into roxygen.
    # > A<-paste0("'",paste(sort(itemAllowed), collapse="'`, `'"),"'");A
    # NOTE: some hand-tweaking must be done to fix linebreaks and (preferably) to
    # change the ' into a ".
    if (!missing(unit)) {
        if (is.list(unit)) {
            unit <- unit[[1]] # second item is a scale
        }
        if (0 == length(unit) || 0 == nchar(unit)) {
            unit <- NULL
        }
    }
    # Previous to 2016-06-11, an error was reported if there was no match.
    itemAllowedMatch <- pmatch(item, itemAllowed)
    if (!is.na(itemAllowedMatch)) {
        item <- itemAllowed[itemAllowedMatch[1]]
    }
    if (getOption("oceUnitBracket") == "[") {
        L <- " ["
        R <- "]"
    } else {
        L <- " ("
        R <- ")"
    }
    if (missing(sep)) {
        tmp <- getOption("oceUnitSep")
        sep <- if (!is.null(tmp)) tmp else ""
    }
    L <- paste(L, sep, sep="")
    R <- paste(sep, R, sep="")
    if (item == "T" || item == "temperature") {
        var <- gettext("Temperature", domain="R-oce")
        if (is.null(unit)) {
            #message("no unit given for temperature")
            full <- bquote(.(var)*.(L)*degree*"C"*.(R))
            abbreviated <- bquote("T"*.(L)*degree*"C"*.(R))
        } else {
            #message("unit given for temperature")
            full <- bquote(.(var)*.(L)*.(unit[[1]])*.(R))
            abbreviated <- bquote("T"*.(L)*.(unit[[1]])*.(R))
        }
    } else if (item == paste("conductivity", "mS/cm")) {
        var <- gettext("Conductivity", domain="R-oce")
        full <- bquote(.(var)*.(L)*mS/cm*.(R))
        abbreviated <- bquote("C"*.(L)*mS/cm*.(R))
    } else if (item == paste("conductivity", "S/m")) {
        var <- gettext("Conductivity", domain="R-oce")
        full <- bquote(.(var)*.(L)*S/m*.(R))
        abbreviated <- bquote("C"*.(L)*S/m*.(R))
    } else if (item == "C") {
        # unitless form
        var <- gettext("Conductivity Ratio", domain="R-oce")
        unit <- gettext("unitless", domain="R-oce")
        full <- bquote(.(var)*.(L)*.(unit[[1]])*.(R))
        abbreviated <- bquote("C")
    } else if (item %in% c("CT",
        paste("conservative", "temperature"),
        paste("Conservative", "Temperature"))) {
        var <- gettext("Conservative Temperature", domain="R-oce")
        full <- bquote(.(var)*.(L)*degree*"C"*.(R))
        abbreviated <- bquote(Theta*.(L)*degree*"C"*.(R))
    } else if (item == "sigmaTheta") {
        var <- gettext("Potential density anomaly", domain="R-oce")
        full <- bquote(.(var)*.(L)*kg/m^3*.(R))
        abbreviated <- bquote(sigma[theta]*.(L)*kg/m^3*.(R))
    } else if (item == "sigma0") {
        var <- gettext("Potential density anomaly wrt surface", domain="R-oce")
        full <- bquote(.(var)*.(L)*kg/m^3*.(R))
        abbreviated <- bquote(sigma[0]*.(L)*kg/m^3*.(R))
    } else if (item == "sigma1") {
        var <- gettext("Potential density anomaly wrt 1000 dbar", domain="R-oce")
        full <- bquote(.(var)*.(L)*kg/m^3*.(R))
        abbreviated <- bquote(sigma[1]*.(L)*kg/m^3*.(R))
    } else if (item == "sigma2") {
        var <- gettext("Potential density anomaly wrt 2000 dbar", domain="R-oce")
        full <- bquote(.(var)*.(L)*kg/m^3*.(R))
        abbreviated <- bquote(sigma[2]*.(L)*kg/m^3*.(R))
    } else if (item == "sigma3") {
        var <- gettext("Potential density anomaly wrt 3000 dbar", domain="R-oce")
        full <- bquote(.(var)*.(L)*kg/m^3*.(R))
        abbreviated <- bquote(sigma[3]*.(L)*kg/m^3*.(R))
    } else if (item == "sigma4") {
        var <- gettext("Potential density anomaly wrt 4000 dbar", domain="R-oce")
        full <- bquote(.(var)*.(L)*kg/m^3*.(R))
        abbreviated <- bquote(sigma[4]*.(L)*kg/m^3*.(R))
    } else if (item %in% c("salinity", "SP")) {
        var <- "Salinity"
        abbreviated <- full <- bquote(.(var))
    } else if (item == "SR") {
        var <- "SR"
        abbreviated <- full <- bquote(S[R]*.(L)*kg/m^3*.(R))
    } else if (item == "Sstar") {
        var <- "Sstar"
        abbreviated <- full <- bquote(S["*"]*.(L)*kg/m^3*.(R))
    } else if (item == "theta") {
        var <- gettext("Potential Temperature", domain="R-oce")
        full <- bquote(.(var)*.(L)*degree*"C"*.(R))
        abbreviated <- bquote(theta*.(L)*degree*"C"*.(R))
    } else if (item == "tritium") {
        var <- gettext("Tritium", domain="R-oce")
        if (is.null(unit)) {
            full <- bquote(.(var)*.(L)*Tu*.(R))
            abbreviated <- bquote(phantom()^3*H*.(L)*Tu*.(R))
        } else {
            full <- bquote(.(var)*.(L)*.(unit[[1]])*.(R))
            abbreviated <- bquote(phantom()^3*H*.(L)*.(unit[[1]])*.(R))
        }
    } else if (item == "N2") {
        # full <- bquote("Square of Buoyancy Frequency"*.(L)*s^-2*.(R))
        full  <- bquote(N^2*.(L)*s^-2*.(R))
        abbreviated <- bquote(N^2*.(L)*s^-2*.(R))
    } else if (item == "nitrate") {
        var <- gettext("Nitrate", domain="R-oce")
        if (is.null(unit)) {
            full <- bquote(.(var)*.(L)*mu*mol/kg*.(R))
            abbreviated <- bquote(N*O[3]*.(L)*mu*mol/kg*.(R))
        } else {
            full <- bquote(.(var)*.(L)*.(unit[[1]])*.(R))
            abbreviated <- bquote(N*O[3]*.(L)*.(unit[[1]])*.(R))
        }
    } else if (item == "nitrite") {
        var <- gettext("Nitrite", domain="R-oce")
        if (is.null(unit)) {
            full <- bquote(.(var)*.(L)*mu*mol/kg*.(R))
            abbreviated <- bquote(N*O[2]*.(L)*mu*mol/kg*.(R))
        } else {
            full <- bquote(.(var)*.(L)*.(unit[[1]])*.(R))
            abbreviated <- bquote(N*O[2]*.(L)*.(unit[[1]])*.(R))
        }
    } else if (item == "oxygen") {
        var <- gettext("Oxygen", domain="R-oce")
        if (is.null(unit)) {
            full <- bquote(.(var))
            abbreviated <- bquote(O[2])
        } else {
            full <- bquote(.(var)*.(L)*.(unit[[1]])*.(R))
            abbreviated <- bquote(O[2]*.(L)*.(unit[[1]])*.(R))
        }
    } else if (item == paste("oxygen", "saturation")) {
        var <- gettext("Oxygen Saturation", domain="R-oce")
        full <- bquote(.(var))
        abbreviated <- bquote(O[2]*.(L)*percent*saturation*.(R))
    } else if (item ==  paste("oxygen", "mL/L")) {
        var <- gettext("Oxygen", domain="R-oce")
        full <- bquote(.(var)*.(L)*mL/L*.(R))
        abbreviated <- bquote(O[2]*.(L)*mL/L*.(R))
    } else if (item == paste("oxygen", "umol/L")) {
        var <- gettext("Oxygen", domain="R-oce")
        full <- bquote(.(var)*.(L)*mu*mol/L*.(R))
        abbreviated <- bquote(O[2]*.(L)*mu*mol/L*.(R))
    } else if (item == paste("oxygen", "umol/kg")) {
        var <- gettext("Oxygen", domain="R-oce")
        full <- bquote(.(var)*.(L)*mu*mol/kg*.(R))
        abbreviated <- bquote(O[2]*.(L)*mu*mol/kg*.(R))
    } else if (item == "phosphate") {
        var <- gettext("Phosphate", domain="R-oce")
        if (is.null(unit)) {
            full <- bquote(.(var)*.(L)*mu*mol/kg*.(R))
            abbreviated <- bquote(P*O[4]*.(L)*mu*mol/kg*.(R))
        } else {
            full <- bquote(.(var)*.(L)*.(unit[[1]])*.(R))
            abbreviated <- bquote(P*O[4]*.(L)*.(unit[[1]])*.(R))
        }
    } else if (item == paste("potential", "temperature")) {
        var <- gettext("Potential Temperature", domain="R-oce")
        if (is.null(unit)) {
            full <- bquote(.(var)*.(L)*degree*C*.(R))
            abbreviated <- bquote(theta*.(L)*degree*C*.(R))
        } else {
            full <- bquote(.(var)*.(L)*.(unit[[1]])*.(R))
            abbreviated <- bquote(theta*.(L)*.(unit[[1]])*.(R))
        }
    } else if (item == "pressure") {
        var <- gettext("Pressure", domain="R-oce")
        if (is.null(unit)) {
            full <- bquote(.(var)*.(L)*dbar*.(R))
            abbreviated <- bquote(theta*.(L)*dbar*.(R))
        } else {
            full <- bquote(.(var)*.(L)*.(unit[[1]])*.(R))
            abbreviated <- bquote(theta*.(L)*.(unit[[1]])*.(R))
        }
    } else if (item == "silicate") {
        var <- gettext("Silicate", domain="R-oce")
        if (is.null(unit)) {
            full <- bquote(.(var)*.(L)*mu*mol/kg*.(R))
            abbreviated <- bquote(Si*O[4]*.(L)*mu*mol/kg*.(R))
        } else {
            full <- bquote(.(var)*.(L)*.(unit[[1]])*.(R))
            abbreviated <- bquote(Si*O[4]*.(L)*.(unit[[1]])*.(R))
        }
    } else if (item == "fluorescence") {
        var <- gettext("Fluorescence", domain="R-oce")
        if (is.null(unit)) {
            # I've no idea what a 'standard' unit might be
            full <- bquote(.(var))
            abbreviated <- full
        } else {
            full <- bquote(.(var)*.(L)*.(unit[[1]])*.(R))
            abbreviated <- bquote("Fluor."*.(L)*.(unit[[1]])*.(R))
        }
    } else if (item == "spice") {
        var <- gettext("Spice", domain="R-oce")
        if (is.null(unit)) {
            full <- bquote(.(var)*.(L)*kg/m^3*.(R))
            abbreviated <- full
        } else {
            full <- bquote(.(var)*.(L)*.(unit[[1]])*.(R))
            abbreviated <- full
        }
    } else if (item == "S") {
        full <- gettext("Salinity", domain="R-oce")
        abbreviated <- expression(S)
    } else if (item %in% c("SA",
        paste("absolute", "salinity"),
        paste("Absolute", "Salinity"))) {
        var <- gettext("Absolute Salinity", domain="R-oce")
        full <- bquote(.(var)*.(L)*g/kg*.(R))
        abbreviated <- bquote(S[A]*.(L)*g/kg*.(R))
    } else if (item == "p") {
        var <- gettext("Pressure", domain="R-oce")
        full <- bquote(.(var)*.(L)*dbar*.(R))
        abbreviated <- bquote("p"*.(L)*dbar*.(R))
    } else if (item == "z") {
        var <- "z"
        abbreviated <- full <- bquote("z"*.(L)*m*.(R))
    } else if (item == "distance") {
        var <- gettext("Distance", domain="R-oce")
        abbreviated <- full <- bquote(.(var)*.(L)*m*.(R))
    } else if (item == "distance km") {
        var <- gettext("Distance", domain="R-oce")
        abbreviated <- full <- bquote(.(var)*.(L)*km*.(R))
    } else if (item == paste("along-spine", "distance", "km")) {
        var <- gettext("Along-spine Distance", domain="R-oce")
        abbreviated <- full <- bquote(.(var)*.(L)*km*.(R))
    } else if (item == paste("along-track", "distance", "km")) {
        var <- gettext("Along-track Distance", domain="R-oce")
        abbreviated <- full <- bquote(.(var)*.(L)*km*.(R))
    } else if (item == "heading") {
        var <- gettext("Heading", domain="R-oce")
        abbreviated <- full <- bquote(.(var)*.(L)*degree*.(R))
    } else if (item == "pitch") {
        var <- gettext("Pitch", domain="R-oce")
        abbreviated <- full <- bquote(.(var)*.(L)*degree*.(R))
    } else if (item == "roll") {
        var <- gettext("Roll", domain="R-oce")
        abbreviated <- full <- bquote(.(var)*.(L)*degree*.(R))
    } else if (item == "u" || item == "v" || item == "w") {
        abbreviated <- full <- bquote(.(item)*.(L)*m/s*.(R))
    } else if (item == "eastward") {
        var <- gettext("Eastward", domain="R-oce")
        abbreviated <- full <- bquote(.(var)*.(L)*m/s*.(R))
    } else if (item == "northward") {
        var <- gettext("Northward", domain="R-oce")
        abbreviated <- full <- bquote(.(var)*.(L)*m/s*.(R))
    } else if (item == "depth") {
        var <- gettext("Depth", domain="R-oce")
        abbreviated <- full <- bquote(.(var)*.(L)*m*.(R))
    } else if (item == "elevation") {
        var <- gettext("Elevation", domain="R-oce")
        abbreviated <- full <- bquote(.(var)*.(L)*m*.(R))
    } else if (item ==  "speed") {
        var <- gettext("Speed", domain="R-oce")
        abbreviated <- full <- bquote(.(var)*.(L)*m/s*.(R))
    } else if (item == "latitude") {
        var <- gettext("Latitude", domain="R-oce")
        # maybe add deg "N" "S" etc here, but maybe not (aesthetics)
        abbreviated <- full <- var
    } else if (item == "longitude") {
        var <- gettext("Longitude", domain="R-oce")
        # maybe add deg "E" "W" etc here, but maybe not (aesthetics)
        abbreviated <- full <- var
    } else if (item == paste("frequency", "cph")) {
        var <- gettext("Frequency", domain="R-oce")
        unit <- gettext("cph", domain="R-oce")
        abbreviated <- full <- bquote(.(var)*.(L)*.(unit[[1]])*.(R))
    } else if (item == paste("sound", "speed")) {
        var <- gettext("Sound Speed", domain="R-oce")
        unit <- gettext("m/s", domain="R-oce")
        abbreviated <- full <- bquote(.(var)*.(L)*.(unit[[1]])*.(R))
    } else if (item == paste("spectral", "density", "m2/cph")) {
        var <- gettext("Spectral density", domain="R-oce")
        full <- bquote(.(var)*.(L)*m^2/cph*.(R))
        var <- gettext("Spec. dens.", domain="R-oce")
        abbreviated <- bquote(.(var)*.(L)*m^2/cph*.(R))
    } else {
        oceDebug(debug, "unknown item=\"", item, "\"\n", sep="")
        if (is.null(unit)) {
            oceDebug(debug, "no unit given\n")
            #message("no unit given")
            full <- item
            abbreviated <- full
        } else {
            oceDebug(debug, "unit \"", unit, "\" given\n")
            full <- bquote(.(item)*.(L)*.(unit[[1]])*.(R))
            abbreviated <- full
        }
    }
    spaceNeeded <- strwidth(paste(full, collapse=""), "inches")
    whichAxis <- if (axis == "x") 1 else 2
    spaceAvailable <- abs(par("fin")[whichAxis])
    fraction <- spaceNeeded / spaceAvailable
    oceDebug(debug, "} # resizableLabel\n", unindent=1, style="bold")
    if (fraction < 1) full else abbreviated
}


#' Rotate velocity components within an oce object
#'
#' Alter the horizontal components of velocities in `adp`,
#' `adv` or `cm` objects, by applying a rotation about
#' the vertical axis.
#'
#' @param x an [adp-class], [adv-class], or [cm-class] object.
#'
#' @param angle The rotation angle about the z axis, in degrees counterclockwise.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' par(mfcol=c(2, 3))
#' # adp (acoustic Doppler profiler)
#' data(adp)
#' plot(adp, which="uv")
#' mtext("adp", side=3, line=0, adj=1, cex=0.7)
#' adpRotated <- rotateAboutZ(adp, 30)
#' plot(adpRotated, which="uv")
#' mtext("adp rotated 30 deg", side=3, line=0, adj=1, cex=0.7)
#' # adv (acoustic Doppler velocimeter)
#' data(adv)
#' plot(adv, which="uv")
#' mtext("adv", side=3, line=0, adj=1, cex=0.7)
#' advRotated <- rotateAboutZ(adv, 125)
#' plot(advRotated, which="uv")
#' mtext("adv rotated 125 deg", side=3, line=0, adj=1, cex=0.7)
#' # cm (current meter)
#' data(cm)
#' plot(cm, which="uv")
#' mtext("cm", side=3, line=0, adj=1, cex=0.7)
#' cmRotated <- rotateAboutZ(cm, 30)
#' plot(cmRotated, which="uv")
#' mtext("cm rotated 30 deg", side=3, line=0, adj=1, cex=0.7)
#'
#' @family things related to adp data
#' @family things related to adv data
#' @family things related to cm data
rotateAboutZ <- function(x, angle)
{
    if (missing(angle))
        stop("must supply angle")
    S <- sin(angle * pi / 180)
    C <- cos(angle * pi / 180)
    rotation <- matrix(c(C, S, -S, C), nrow=2)
    res <- x
    allowedClasses <- c("adp", "adv", "cm")
    if (!(class(x) %in% allowedClasses))
        stop("cannot rotate for class \"", class(x), "\"; try one of: \"",
            paste(allowedClasses, collapse="\" \""), "\")")
    if (inherits(x, "adp")) {
        if (is.ad2cp(x))
            stop("this function does not work yet for AD2CP data")
        if (x[["oceCoordinate"]] != "enu")
            stop("cannot rotate adp unless coordinate system is 'enu'; see ?toEnu or ?xyzToEnu")
        V <- x[["v"]]
        # Work through the bins, transforming a 3D array operation to a
        # sequence of 2D matrix operations.
        for (j in seq_len(dim(V)[2])) {
            uvr <- rotation %*% t(V[, j, 1:2])
            V[, j, 1] <- uvr[1, ]
            V[, j, 2] <- uvr[2, ]
        }
        res@data$v <- V
    } else if (inherits(x, "adv")) {
        if (x[["oceCoordinate"]] != "enu")
            stop("cannot rotate adv unless coordinate system is 'enu'; see ?toEnu or ?xyzToEnu")
        V <- x[["v"]]
        uvr <- rotation %*% t(V[, 1:2])
        V[, 1] <- uvr[1, ]
        V[, 2] <- uvr[2, ]
        res@data$v <- V
    } else if (inherits(x, "cm")) {
        uvr <- rotation %*% rbind(x@data$u, x@data$v)
        res@data$u <- uvr[1, ]
        res@data$v <- uvr[2, ]
    } else {
        stop("cannot rotate for class \"", class(x), "\"; try one of: \"",
            paste(allowedClasses, collapse="\" \""), "\". (internal error: please report)")
    }
    # Update processing log
    res@processingLog <- processingLogAppend(res@processingLog,
        paste("rotateAboutZ(x, angle=", angle, ")", sep=""))
    res
}

#' Format a latitude-longitude pair
#'
#' Format a latitude-longitude pair, using "S" for negative latitudes, etc.
#'
#' @param lat latitude in \eqn{^\circ}{deg}N north of the equator.
#'
#' @param lon longitude in \eqn{^\circ}{deg}N east of Greenwich.
#'
#' @param digits the number of significant digits to use when printing.
#'
#' @return A character string.
#'
#' @author Dan Kelley
#'
#' @seealso [latFormat()] and [lonFormat()].
latlonFormat <- function(lat, lon, digits=max(6, getOption("digits") - 1))
{
    n <- length(lon)
    res <- vector("character", n)
    for (i in 1:n) {
        if (is.na(lat[i]) || is.na(lon[i])) {
            res[i] <- "Lat and lon unknown"
        } else {
            res[i] <- paste(format(abs(lat[i]), digits=digits),
                if (lat[i] > 0) gettext("N", domain="R-oce") else gettext("S", domain="R-oce"),
                " ",
                format(abs(lon[i]), digits=digits),
                if (lon[i] > 0) gettext("E", domain="R-oce") else gettext("W", domain="R-oce"),
                sep="")
        }
    }
    res
}


#' Format a latitude
#'
#' Format a latitude, using "S" for negative latitude.
#'
#' @param lat latitude in \eqn{^\circ}{deg}N north of the equator.
#'
#' @param digits the number of significant digits to use when printing.
#'
#' @return A character string.
#'
#' @author Dan Kelley
#'
#' @seealso [lonFormat()] and [latlonFormat()].
latFormat <- function(lat, digits=max(6, getOption("digits") - 1))
{
    n <- length(lat)
    if (n < 1)
        return("")
    res <- vector("character", n)
    for (i in 1:n) {
        if (is.na(lat[i])) {
            res[i] <-  ""
        } else {
            res[i] <- paste(format(abs(lat[i]), digits=digits),
                gettext(if (lat[i] > 0.0) "N" else "S", domain="R-oce"),
                sep="")
        }
    }
    res
}


#' Format a longitude
#'
#' Format a longitude, using "W" for west longitude.
#'
#' @param lon longitude in \eqn{^\circ}{deg}N east of Greenwich.
#'
#' @param digits the number of significant digits to use when printing.
#'
#' @return A character string.
#'
#' @author Dan Kelley
#'
#' @seealso [latFormat()] and [latlonFormat()].
lonFormat <- function(lon, digits=max(6, getOption("digits") - 1))
{
    n <- length(lon)
    if (n < 1)
        return("")
    res <- vector("character", n)
    for (i in 1:n)
    if (is.na(lon[i])) {
        res[i] <- ""
    } else {
        res[i] <- paste(format(abs(lon[i]), digits=digits),
            gettext(if (lon[i] > 0.0) "E" else "W", domain="R-oce"),
            sep="")
    }
    res
}

#' Alter longitudes from -180:180 to 0:360 convention
#'
#' For numerical input, including vectors, matrices and arrays,
#' [lon360()] simply calls [ifelse()] to add 360 to any negative values. For
#' [section-class] objects, it changes `longitude` in the `metadata` slot
#' and then calls itself to handle the [ctd-class] objects stored as
#' as the entries in `station` within the `data` slot.
#' For this [ctd-class] object, and indeed for all non-[section-class]
#' objects, [lon360()] changes `longitude` values in the
#' `metadata` slot (if present) and also in the `data` slot (again, if
#' present). This function is not useful for dealing with coastline
#' data; see [coastlineCut()] for such data.
#'
#' @param x either a numeric vector or array, or an [oce-class] object.
#'
#' @examples
#' lon360(c(179, -179))
lon360 <- function(x)
{
    if (missing(x))
        stop("must provide 'x'")
    shift <- function(X) ifelse(X < 0, 360 + X, X)
    res <- x
    if (inherits(x, "section")) {
        # Shift section-level metadata.
        res@metadata$longitude <- shift(x@metadata$longitude)
        # Shift each station, by calling lon360() on that (CTD) object.
        for (istn in seq_along(x[["station"]])) {
            res@data$station[[istn]] <- lon360(x@data$station[[istn]])
        }
        res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep="", collapse=""))
    } else if (inherits(x, "oce")) {
        if ("longitude" %in% names(x[["metadata"]]))
            res@metadata$longitude <- shift(res@metadata$longitude)
        if ("longitude" %in% names(x[["data"]]))
            res@data$longitude <- shift(res@data$longitude)
        res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep="", collapse=""))
    } else if (is.numeric(x)) {
        res <- shift(x)
    } else {
        stop("x must be a numeric value of an oce object")
    }
    res
}

#' Determine time offset from timezone
#'
#' The data are from
#' `https://www.timeanddate.com/time/zones/` and were
#' hand-edited to develop this code, so there may be errors.  Also, note that
#' some of these contradict; if you examine the code, you'll see some
#' commented-out portions that represent solving conflicting definitions by
#' choosing the more common timezone abbreviation over a the less common one.
#'
#' @param tz a timezone, e.g. `UTC`.
#'
#' @return Number of hours in offset, e.g. `AST` yields 4.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' cat("Atlantic Standard Time is ", GMTOffsetFromTz("AST"), "hours after UTC")
#' @family things relating to time
GMTOffsetFromTz <- function(tz)
{
    # Data are from
    #   https://www.timeanddate.com/library/abbreviations/timezones/
    # and hand-edited, so there may be errors.  Also, note that some of these
    # contradict ... I've commented out conflicting definitions that I think
    # will come up most rarely in use, but perhaps something better should
    # be devised.  (Maybe this is not a problem.  Maybe only MEDS uses these,
    # as opposed to GMT offsets, and maybe they only work in 5 zones, anyway...)
    if (tz == "A")      return(-1) # Alpha Time Zone Military                UTC + 1 hour
    if (tz == "ACDT")   return(-10.5) # Aus. Central Daylight Time  Aus.     UTC + 10:30 hours
    if (tz == "ACST")   return(-9.5) # Aus. Central Standard Time Aus.ralia  UTC + 9:30 hours
    if (tz == "ADT")    return(3) # Atlantic Daylight Time  North America    UTC - 3 hours
    if (tz == "AEDT")   return(-11) # Aus. E Day. or Aus. E Sum Time Aus.    UTC + 11 hours
    if (tz == "AEST")   return(-10) # Australian Eastern Standard Time  Aus. UTC + 10 hours
    if (tz == "AKDT")   return(8) # Alaska Daylight Time    North America    UTC - 8 hours
    if (tz == "AKST")   return(9) # Alaska Standard Time    North America    UTC - 9 hours
    if (tz == "AST")    return(4) # Atlantic Standard Time  North America    UTC - 4 hours
    if (tz == "AWDT")   return(-9) # Australian Western Daylight Time Aus.   UTC + 9 hours
    if (tz == "AWST")   return(-8) # Australian Western Standard Time Aus.   UTC + 8 hours
    if (tz == "B")      return(-2) # Bravo Time Zone Military                UTC + 2 hours
    if (tz == "BST")    return(-1) # British Summer Time     Europe          UTC + 1 hour
    if (tz == "C")      return(-3) # Charlie Time Zone       Military        UTC + 3 hours
    #if (tz == "CDT")  return(-10.5) # Central Daylight Time   Australia    UTC + 10:30 hours
    if (tz == "CDT")    return(5) # Central Daylight Time   North America    UTC - 5 hours
    if (tz == "CEDT")   return(-2) # Central European Daylight Time  Europe  UTC + 2 hours
    if (tz == "CEST")   return(-2) # Central European Summer Time    Europe  UTC + 2 hours
    if (tz == "CET")    return(-1) # Central European Time   Europe          UTC + 1 hour
    #if (tz == "CST")  return(-10.5) # Central Summer Time     Australia    UTC + 10:30 hours
    #if (tz == "CST")  return(-9.5) # Central Standard Time   Australia     UTC + 9:30 hours
    if (tz == "CST")    return(6) # Central Standard Time   North America    UTC - 6 hours
    if (tz == "CXT")    return(-7) # Christmas Island Time   Australia       UTC + 7 hours
    if (tz == "D")      return(-4) # Delta Time Zone Military                UTC + 4 hours
    if (tz == "E")      return(-5) # Echo Time Zone  Military                UTC + 5 hours
    #if (tz == "EDT")  return(-11) # Eastern Daylight Time   Australia      UTC + 11 hours
    if (tz == "EDT")    return(4) # Eastern Daylight Time   North America    UTC - 4 hours
    if (tz == "EEDT")   return(-3) # Eastern European Daylight Time  Europe  UTC + 3 hours
    if (tz == "EEST")   return(-3) # Eastern European Summer Time    Europe  UTC + 3 hours
    if (tz == "EET")    return(-2) # Eastern European Time   Europe          UTC + 2 hours
    #if (tz == "EST")  return(-11) # Eastern Summer Time     Australia      UTC + 11 hours
    #if (tz == "EST")  return(-10) # Eastern Standard Time   Australia      UTC + 10 hours
    if (tz == "EST")    return(5) # Eastern Standard Time   North America    UTC - 5 hours
    if (tz == "F")      return(-6) # Foxtrot Time Zone       Military        UTC + 6 hours
    if (tz == "G")      return(-7) # Golf Time Zone  Military                UTC + 7 hours
    if (tz == "GMT")    return(0) # Greenwich Mean Time     Europe           UTC
    if (tz == "H")      return(-8) # Hotel Time Zone Military                UTC + 8 hours
    if (tz == "HAA")    return(3) # Heure Avancee de l'Atlantique N. Amer.   UTC - 3 hours
    if (tz == "HAC")    return(5) # Heure Avancee du Centre North America    UTC - 5 hours
    if (tz == "HADT")   return(9) # Hawaii-Aleutian Daylight Time N. Amer.   UTC - 9 hours
    if (tz == "HAE")    return(4) # Heure Avancee de l'Est  North America    UTC - 4 hours
    if (tz == "HAP")    return(7) # Heure Avancee du Pacifique  N. America   UTC - 7 hours
    if (tz == "HAR")    return(6) # Heure Avancee des Rocheuses N. America   UTC - 6 hours
    if (tz == "HAST")   return(10) # Hawaii-Aleutian Standard Time N. Amer.  UTC - 10 hours
    if (tz == "HAT")    return(2.5) # Heure Avancee de Terre-Neuve N. Amer.  UTC - 2:30 hours
    if (tz == "HAY")    return(8) # Heure Avancee du Yukon  North America    UTC - 8 hours
    if (tz == "HNA")    return(4) # Heure Normaee de l'Atlantique N. Amer.   UTC - 4 hours
    if (tz == "HNC")    return(6) # Heure Normale du Centre North America    UTC - 6 hours
    if (tz == "HNE")    return(5) # Heure Normale de l'Est  North America    UTC - 5 hours
    if (tz == "HNP")    return(8) # Heure Normale du Pacifique  N. America   UTC - 8 hours
    if (tz == "HNR")    return(7) # Heure Normale des Rocheuses N. America   UTC - 7 hours
    if (tz == "HNT")    return(3.5) # Heure Normale de Terre-Neuve N.Amer.   UTC - 3:30 hours
    if (tz == "HNY")    return(9) # Heure Normale du Yukon  North America    UTC - 9 hours
    if (tz == "I")      return(-9) # India Time Zone Military                UTC + 9 hours
    if (tz == "IST")    return(-1) # Irish Summer Time       Europe          UTC + 1 hour
    if (tz == "K")      return(-10) # Kilo Time Zone  Military               UTC + 10 hours
    if (tz == "L")      return(-11) # Lima Time Zone  Military               UTC + 11 hours
    if (tz == "M")      return(-12) # Mike Time Zone  Military               UTC + 12 hours
    if (tz == "MDT")    return(6) # Mountain Daylight Time  North America    UTC - 6 hours
    if (tz == "MESZ")   return(-2) # Mitteleuroaische Sommerzeit Europe      UTC + 2 hours
    if (tz == "MEZ")    return(-1) # Mitteleuropaische Zeit  Europe          UTC + 1 hour
    if (tz == "MST")    return(7) # Mountain Standard Time  North America    UTC - 7 hours
    if (tz == "N")      return(1) # November Time Zone      Military         UTC - 1 hour
    if (tz == "NDT")    return(2.5) # NFLD Daylight Time North America       UTC - 2:30 hours
    if (tz == "NFT")    return(-11.5) # Norfolk (Island) Time   Australia    UTC + 11:30 hours
    if (tz == "NST")    return(3.5) # NFLD Std. Time North America           UTC - 3:30 hours
    if (tz == "O")      return(1) # Oscar Time Zone Military                 UTC - 2 hours
    if (tz == "P")      return(3) # Papa Time Zone  Military                 UTC - 3 hours
    if (tz == "PDT")    return(7) # Pacific Daylight Time   North America    UTC - 7 hours
    if (tz == "PST")    return(8) # Pacific Standard Time   North America    UTC - 8 hours
    if (tz == "Q")      return(4) # Quebec Time Zone        Military         UTC - 4 hours
    if (tz == "R")      return(4) # Romeo Time Zone Military                 UTC - 5 hours
    if (tz == "S")      return(6) # Sierra Time Zone        Military         UTC - 6 hours
    if (tz == "T")      return(7) # Tango Time Zone Military                 UTC - 7 hours
    if (tz == "U")      return(8) # Uniform Time Zone       Military         UTC - 8 hours
    if (tz == "UTC")    return(0) # Coordinated Universal Time      Europe   UTC
    if (tz == "V")      return(9) # Victor Time Zone        Military         UTC - 9 hours
    if (tz == "W")      return(10) # Whiskey Time Zone       Military        UTC - 10 hours
    if (tz == "WDT")    return(-9) # Western Daylight Time   Australia       UTC + 9 hours
    if (tz == "WEDT")   return(-1) # Western European Daylight Time  Europe  UTC + 1 hour
    if (tz == "WEST")   return(-1) # Western European Summer Time    Europe  UTC + 1 hour
    if (tz == "WET")    return(0) # Western European Time   Europe  UTC
    #if (tz == "WST")  return(-9) # Western Summer Time     Australia       UTC + 9 hours
    if (tz == "WST")    return(-8) # Western Standard Time   Australia       UTC + 8 hours
    if (tz == "X")      return(11) # X-ray Time Zone Military                UTC - 11 hours
    if (tz == "Y")      return(12) # Yankee Time Zone        Military        UTC - 12 hours
    if (tz == "Z")      return(0) # Zulu Time Zone  Military                 UTC
}


#' Acceleration due to earth gravity
#'
#' Compute \eqn{g}{g}, the acceleration due to gravity, as a function of
#' latitude.
#'
#' Value not verified yet, except roughly.
#'
#' @param latitude Latitude in \eqn{^\circ}{deg}N or radians north of the
#' equator.
#'
#' @param degrees Flag indicating whether degrees are used for latitude; if set
#' to `FALSE`, radians are used.
#'
#' @return Acceleration due to gravity, in \eqn{m^2/s}{m^2/s}.
#'
#' @author Dan Kelley
#'
#' @references Gill, A.E., 1982. *Atmosphere-ocean Dynamics*, Academic
#' Press, New York, 662 pp.
#'
#' **Caution:** Fofonoff and Millard (1983 UNESCO) use a different formula.
#'
#' @examples
#' g <- gravity(45) # 9.8
gravity <- function(latitude=45, degrees=TRUE)
{
    if (degrees)
        latitude <- latitude * 0.0174532925199433
    9.780318 * (1.0 + 5.3024e-3 * sin(latitude)^2 - 5.9e-6 * sin(2*latitude)^2)
}


#' Make a digital filter
#'
#' The filter is suitable for use by [filter()],
#' [convolve()] or (for the `asKernal=TRUE` case) with
#' [kernapply()].  Note that [convolve()] should be faster
#' than [filter()], but it cannot be used if the time series has
#' missing values.  For the Blackman-Harris filter, the half-power frequency is
#' at `1/m` cycles per time unit, as shown in the \dQuote{Examples}
#' section.  When using [filter()] or [kernapply()] with
#' these filters, use `circular=TRUE`.
#'
#' @param type a string indicating the type of filter to use.  (See Harris
#' (1978) for a comparison of these and similar filters.)
#'
#' * `"blackman-harris"` yields a modified raised-cosine filter designated
#' as "4-Term (-92 dB) Blackman-Harris" by Harris (1978; coefficients given in
#' the table on page 65).  This is also called "minimum 4-sample Blackman
#' Harris" by that author, in his Table 1, which lists figures of merit as
#' follows: highest side lobe level -92dB; side lobe fall off -6 db/octave;
#' coherent gain 0.36; equivalent noise bandwidth 2.00 bins; 3.0-dB bandwidth
#' 1.90 bins; scallop loss 0.83 dB; worst case process loss 3.85 dB; 6.0-db
#' bandwidth 2.72 bins; overlap correlation 46 percent for 75\% overlap and 3.8
#' for 50\% overlap.  Note that the equivalent noise bandwidth is the width of
#' a spectral peak, so that a value of 2 indicates a cutoff frequency of
#' `1/m`, where `m` is as given below.
#'
#' * `"rectangular"` for a flat filter.  (This is just for convenience.  Note that
#' [`kernel`]`("daniell",....)` gives the same result, in kernel form.)
#' `"hamming"` for a Hamming filter (a raised-cosine that does not taper
#' to zero at the ends)
#'
#' * `"hann"` (a raised cosine that tapers to zero at the ends).
#'
#' @param m length of filter.  This should be an odd number, for any
#' non-rectangular filter.
#'
#' @param asKernel boolean, set to `TRUE` to get a smoothing kernel for
#' the return value.
#'
#' @return If `asKernel` is `FALSE`, this returns a list of filter
#' coefficients, symmetric about the midpoint and summing to 1.  These may be
#' used with [filter()], which should be provided with argument
#' `circular=TRUE` to avoid phase offsets.  If `asKernel` is
#' `TRUE`, the return value is a smoothing kernel, which can be applied to
#' a timeseries with [kernapply()], whose bandwidth can be determined
#' with [bandwidth.kernel()], and which has both print and plot
#' methods.
#'
#' @author Dan Kelley
#'
#' @references F. J. Harris, 1978.  On the use of windows for harmonic analysis
#' with the discrete Fourier Transform.  *Proceedings of the IEEE*, 66(1),
#' 51-83 (`http://web.mit.edu/xiphmont/Public/windows.pdf`.)
#'
#' @examples
#' library(oce)
#'
#' # 1. Demonstrate step-function response
#' y <- c(rep(1, 10), rep(-1, 10))
#' x <- seq_along(y)
#' plot(x, y, type="o", ylim=c(-1.05, 1.05))
#' BH <- makeFilter("blackman-harris", 11, asKernel=FALSE)
#' H <- makeFilter("hamming", 11, asKernel=FALSE)
#' yBH <- stats::filter(y, BH)
#' points(x, yBH, col=2, type="o")
#' yH <- stats::filter(y, H)
#' points(yH, col=3, type="o")
#' legend("topright", col=1:3, cex=2/3, pch=1,
#'        legend=c("input", "Blackman Harris", "Hamming"))
#'
#' # 2. Show theoretical and practical filter gain, where
#' #    the latter is based on random white noise, and
#' #    includes a particular value for the spans
#' #    argument of spectrum(), etc.
#'\dontrun{
#' # need signal package for this example
#' r <- rnorm(2048)
#' rh <- stats::filter(r, H)
#' rh <- rh[is.finite(rh)] # kludge to remove NA at start/end
#' sR <- spectrum(r, plot=FALSE, spans=c(11, 5, 3))
#' sRH <- spectrum(rh, plot=FALSE, spans=c(11, 5, 3))
#' par(mfrow=c(2, 1), mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0))
#' plot(sR$freq, sRH$spec/sR$spec, xlab="Frequency", ylab="Power Transfer",
#'      type="l", lwd=5, col="gray")
#' theory <- freqz(H, n=seq(0,pi,length.out=100))
#' # Note we must square the modulus for the power spectrum
#' lines(theory$f/pi/2, Mod(theory$h)^2, lwd=1, col="red")
#' grid()
#' legend("topright", col=c("gray", "red"), lwd=c(5, 1), cex=2/3,
#'        legend=c("Practical", "Theory"), bg="white")
#' plot(log10(sR$freq), log10(sRH$spec/sR$spec),
#'      xlab="log10 Frequency", ylab="log10 Power Transfer",
#'      type="l", lwd=5, col="gray")
#' theory <- freqz(H, n=seq(0,pi,length.out=100))
#' # Note we must square the modulus for the power spectrum
#' lines(log10(theory$f/pi/2), log10(Mod(theory$h)^2), lwd=1, col="red")
#' grid()
#' legend("topright", col=c("gray", "red"), lwd=c(5, 1), cex=2/3,
#'        legend=c("Practical", "Theory"), bg="white")
#'}
makeFilter <- function(type=c("blackman-harris", "rectangular", "hamming", "hann"), m, asKernel=TRUE)
{
    type <- match.arg(type)
    if (missing(m))
        stop("must supply 'm'")
    i <- seq(0, m - 1)
    if (type == "blackman-harris") {
        # See Harris (1978) table on p65
        if (m == 2 * floor(m/2)) {
            m <- m + 1
            warning("increased filter length by 1, to make it odd")
        }
        a <- c(0.35875, 0.488829, 0.14128, 0.01168) # 4-term (-92dB) coefficients
        ff <- pi * i / (m - 1)
        coef <- a[1] - a[2]*cos(2*ff) + a[3]*cos(4*ff) - a[4]*cos(6*ff)
    } else if (type == "rectangular") {
        coef <- rep(1 / m, m)
    } else if (type == "hamming") {
        coef <- 0.54 - 0.46 * cos(2 * pi * i / (m-1))
    } else if (type == "hann") {
        coef <- 0.50 - 0.50 * cos(2 * pi * i / (m-1))
    }
    coef <- coef / sum(coef)           # ensure unit sum
    if (!asKernel)
        return(coef)
    if (m == 2 * floor(m/2))
        stop("m must be odd")
    middle <- ceiling(m / 2)
    coef <- coef[middle:m]
    # the r=0 is to prevent code-analysis warning; it only applies to Fejer, which we do not use
    return(kernel(coef=coef, name=paste(type, "(", m, ")", sep=""), r=0))
}


#' 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


#' Grid data using Barnes algorithm
#'
#' The algorithm follows that described by Koch et al. (1983), except
#' that `interpBarnes` adds (1) the ability to
#' blank out the grid where data are
#' sparse, using the `trim` argument, and (2) the ability to
#' pre-grid, with the `pregrid` argument.
#'
#' @param x,y a vector of x and y locations.
#'
#' @param z a vector of z values, one at each (x,y) location.
#'
#' @param w a optional vector of weights at the (x,y) location.  If not
#' supplied, then a weight of 1 is used for each point, which means equal
#' weighting.  Higher weights give data points more influence. If `pregrid`
#' is `TRUE`, then any supplied value of `w` is ignored, and instead
#' each of the pregriddd points is given equal weight.
#'
#' @param xg,yg optional vectors defining the x and y grids.  If not supplied,
#' these values are inferred from the data, using e.g. `pretty(x, n=50)`.
#'
#' @param xgl,ygl optional lengths of the x and y grids, to be constructed with
#' [seq()] spanning the data range.  These values `xgl` are only
#' examined if `xg` and `yg` are not supplied.
#'
#' @param xr,yr optional values defining the x and y radii of the weighting ellipse.
#' If not supplied, these are calculated as the span of x
#' and y over the square root of the number of data.
#'
#' @param gamma grid-focussing parameter.  At each successive iteration, `xr` and
#' `yr` are reduced by a factor of `sqrt(gamma)`.
#'
#' @param iterations number of iterations.  Set this to 1 to perform just
#' one iteration, using the radii as described at `xr,yr` above.
#'
#' @param trim a number between 0 and 1, indicating the quantile of data weight
#' to be used as a criterion for blanking out the gridded value (using
#' `NA`).  If 0, the whole `zg` grid is returned.  If >0, any spots
#' on the grid where the data weight is less than the `trim`-th
#' [quantile()] are set to `NA`.  See examples.
#'
#' @param pregrid an indication of whether to pre-grid the data. If
#' `FALSE`, this is not done, i.e. conventional Barnes interpolation is
#' performed.  Otherwise, then the data are first averaged within grid cells
#' using [binMean2D()].  If `pregrid` is `TRUE` or
#' `4`, then this averaging is done within a grid that is 4 times finer
#' than the grid that will be used for the Barnes interpolation. Otherwise,
#' `pregrid` may be a single integer indicating the grid refinement (4
#' being the result if `TRUE` had been supplied), or a vector of two
#' integers, for the grid refinement in x and y. The purpose of using
#' `pregrid` is to speed processing on large datasets, and to remove
#' spatial bias (e.g. with a single station that is repeated frequently in an
#' otherwise seldom-sampled region).  A form of pregridding is done in the
#' World Ocean Atlas, for example.
#'
#' @param debug a flag that turns on debugging.  Set to 0 for no debugging
#' information, to 1 for more, etc; the value is reduced by 1 for each
#' descendent function call.
#'
#' @return A list containing: `xg`, a vector holding the x-grid);
#' `yg`, a vector holding the y-grid; `zg`, a matrix holding the
#' gridded values; `wg`, a matrix holding the weights used in the
#' interpolation at its final iteration; and `zd`, a vector of the same
#' length as `x`, which holds the interpolated values at the data points.
#'
#' @author Dan Kelley
#'
#' @seealso See [wind()].
#'
#' @references S. E.  Koch and M.  DesJardins and P. J. Kocin, 1983.  ``An
#' interactive Barnes objective map analysis scheme for use with satellite and
#' conventional data,'' *J.  Climate Appl.  Met.*, vol 22, p. 1487-1503.
#'
#' @examples
#' library(oce)
#'
#' # 1. contouring example, with wind-speed data from Koch et al. (1983)
#' data(wind)
#' u <- interpBarnes(wind$x, wind$y, wind$z)
#' contour(u$xg, u$yg, u$zg, labcex=1)
#' text(wind$x, wind$y, wind$z, cex=0.7, col="blue")
#' title("Numbers are the data")
#'
#' # 2. As 1, but blank out spots where data are sparse
#' u <- interpBarnes(wind$x, wind$y, wind$z, trim=0.1)
#' contour(u$xg, u$yg, u$zg, level=seq(0, 30, 1))
#' points(wind$x, wind$y, cex=1.5, pch=20, col="blue")
#'
#' # 3. As 1, but interpolate back to points, and display the percent mismatch
#' u <- interpBarnes(wind$x, wind$y, wind$z)
#' contour(u$xg, u$yg, u$zg, labcex=1)
#' mismatch <- 100 * (wind$z - u$zd) / wind$z
#' text(wind$x, wind$y, round(mismatch), col="blue")
#' title("Numbers are percent mismatch between grid and data")
#'
#' # 4. As 3, but contour the mismatch
#' mismatchGrid <- interpBarnes(wind$x, wind$y, mismatch)
#' contour(mismatchGrid$xg, mismatchGrid$yg, mismatchGrid$zg, labcex=1)
#'
#' # 5. One-dimensional example, smoothing a salinity profile
#' data(ctd)
#' p <- ctd[["pressure"]]
#' y <- rep(1, length(p)) # fake y data, with arbitrary value
#' S <- ctd[["salinity"]]
#' pg <- pretty(p, n=100)
#' g <- interpBarnes(p, y, S, xg=pg, xr=1)
#' plot(S, p, cex=0.5, col="blue", ylim=rev(range(p)))
#' lines(g$zg, g$xg, col="red")
interpBarnes <- function(x, y, z, w,
    xg, yg, xgl, ygl,
    xr, yr, gamma=0.5, iterations=2, trim=0,
    pregrid=FALSE,
    debug=getOption("oceDebug"))
{
    debug <- max(0, debug)
    oceDebug(debug, "interpBarnes(",
        argShow(x),
        argShow(y),
        argShow(z),
        argShow(w),
        argShow(xg),
        argShow(xg),
        argShow(yg),
        argShow(xgl),
        argShow(ygl),
        argShow(xr),
        argShow(yr),
        argShow(gamma),
        argShow(iterations),
        argShow(trim),
        argShow(pregrid, last=TRUE),
        ") {\n", unindent=1, sep="")
    if (!is.vector(x))
        stop("x must be a vector")
    n <- length(x)
    if (length(y) != n)
        stop("lengths of x and y disagree; they are ", n, " and ", length(y))
    if (length(z) != n)
        stop("lengths of x and z disagree; they are ", n, " and ", length(z))
    xrGiven <- !missing(xr)
    yrGiven <- !missing(yr)
    if (missing(w))
        w <- rep(1.0, length(x))
    if (missing(xg)) {
        if (missing(xgl)) {
            if (0 == diff(range(x, na.rm=TRUE))) {
                xg <- x[1]
            } else {
                xg <- pretty(x, n=50)
            }
        } else {
            xg <- seq(min(x, na.rm=TRUE), max(x, na.rm=TRUE), length.out=xgl)
        }
        oceDebug(debug, "computed ", vectorShow(xg))
    }
    if (missing(yg)) {
        if (missing(ygl)) {
            if (0 == diff(range(y, na.rm=TRUE))) {
                yg <- y[1]
            } else {
                yg <- pretty(y, n=50)
            }
        } else {
            yg <- seq(min(y, na.rm=TRUE), max(y, na.rm=TRUE), length.out=ygl)
        }
        oceDebug(debug, "computed ", vectorShow(yg))
    }
    if (!xrGiven) {
        xr <- diff(range(x, na.rm=TRUE)) / sqrt(n)
        if (xr == 0) {
            xr <- 1
        }
        oceDebug(debug, "computed xr=", xr, " based on data density\n")
    }
    if (!yrGiven) {
        yr <- diff(range(y, na.rm=TRUE)) / sqrt(n)
        if (yr == 0) {
            yr <- 1
        }
        oceDebug(debug, "computed yr=", yr, " based on data density\n")
    }
    # Handle pre-gridding (code not DRY but short enough to be ok)
    if (is.logical(pregrid)) {
        if (pregrid) {
            pregrid <- c(4, 4)
            oceDebug(debug, "pregrid: ", paste(pregrid, collapse=" "))
            pg <- binMean2D(x, y, z,
                xbreaks=seq(xg[1], tail(xg, 1), (xg[2]-xg[1]) / pregrid[1]),
                ybreaks=seq(yg[1], tail(yg, 1), (yg[2]-yg[1]) / pregrid[2]),
                flatten=TRUE)
            x <- pg$x
            y <- pg$y
            z <- pg$f
            w <- rep(1, length(x))
        }
    } else {
        if (!is.numeric(pregrid))
            stop("pregrid must be logical or a numeric vector")
        if (length(pregrid) < 0 || length(pregrid) > 2)
            stop("length(pregrid) must be 1 or 2")
        if (length(pregrid) == 1)
            pregrid <- rep(pregrid, 2)
        oceDebug(debug, "pregrid: ", paste(pregrid, collapse=" "))
        pg <- binMean2D(x, y, z,
            xbreaks=seq(xg[1], tail(xg, 1), (xg[2]-xg[1])/pregrid[1]),
            ybreaks=seq(yg[1], tail(yg, 1), (yg[2]-yg[1])/pregrid[2]),
            flatten=TRUE)
        x <- pg$x
        y <- pg$y
        z <- pg$f
    }
    for (i in seq_len(iterations)) {
        oceDebug(debug, "  Iteration ", i, ": use xr=", xr*gamma^((i-1)/2), " and yr=", yr*gamma^((i-1)/2), "\n")
    }
    ok <- !is.na(x) & !is.na(y) & !is.na(z) & !is.na(w)
    if (sum(ok) > 0) {
        g <- do_interp_barnes(x[ok], y[ok], z[ok], w[ok], xg, yg, xr, yr, gamma, iterations)
        if (trim >= 0 && trim <= 1) {
            bad <- g$wg < quantile(g$wg, trim, na.rm=TRUE)
            g$zg[bad] <- NA
        }
        rval <- list(xg=xg, yg=yg, zg=g$zg, wg=g$wg, zd=g$zd)
    } else {
        rval <- list(xg=xg, yg=yg,
            zg=matrix(NA, nrow=length(xg), ncol=length(yg)),
            wg=matrix(NA, nrow=length(xg), ncol=length(yg)),
            zd=rep(NA, length(x)))
    }
    oceDebug(debug, sprintf("filled %.3f%% of z matrix\n", 100*sum(is.finite(rval$zg))/prod(dim(rval$zg))))
    oceDebug(debug, "} # interpBarnes(...)\n", unindent=1, sep="")
    rval
}

#' Coriolis parameter on rotating earth
#'
#' Compute \eqn{f}{f}, the Coriolis parameter as a function of latitude
#' (see reference 1),
#' assuming earth siderial angular rotation rate
#' \eqn{omega}{omega}=7292115e-11 rad/s. See reference 1 for general notes, and
#' see reference 2 for comments on temporal variations
#' of \eqn{omega}{omega}.
#'
#' @param latitude Vector of latitudes in \eqn{^\circ}{deg}N or radians north of the equator.
#'
#' @param degrees Flag indicating whether degrees are used for latitude; if set
#' to `FALSE`, radians are used.
#'
#' @return Coriolis parameter, in radian/s.
#'
#' @author Dan Kelley
#'
#' @references
#' 1. Gill, A.E., 1982. *Atmosphere-ocean Dynamics*, Academic
#' Press, New York, 662 pp.
#'
#' 2. Groten, E., 2004: Fundamental Parameters and Current, 2004. Best
#'  Estimates of the Parameters of Common Relevance to Astronomy, Geodesy,
#'  and Geodynamics. Journal of Geodesy, 77:724-797.
#'  (downloaded from
#' `http://www.iag-aig.org/attach/e354a3264d1e420ea0a9920fe762f2a0/51-groten.pdf`
#' March 11, 2017).
#'
#' @examples
#' C <- coriolis(45) # 1e-4
coriolis <- function(latitude, degrees=TRUE)
{
    # Siderial day 86164.1 s.
    if (degrees) latitude <- latitude * 0.0174532925199433
    # http://www.iag-aig.org/attach/e354a3264d1e420ea0a9920fe762f2a0/51-groten.pdf 7292115e-11
    2 * 7292115e-11 * sin(latitude)
}


#' Correct for drift in instrument clock
#'
#' It is assumed that the instrument clock matches the real time at the start
#' of the sampling, and that the clock drifts linearly (i.e. is uniformly fast
#' or slow) over the sampling interval.  Linear interpolation is used to infer
#' the values of all variables in the `data` slot.  The data length is
#' altered in this process, e.g. a slow instrument clock (positive
#' `slowEnd`) takes too few samples in a given time interval, so
#' `undriftTime` will increase the number of data.
#'
#' @param x an [oce-class] object.
#'
#' @param slowEnd number of seconds to add to final instrument time, to get the
#' correct time of the final sample.  This will be a positive number, for a
#' "slow" instrument clock.
#'
#' @param tname Character string indicating the name of the time column in the
#' `data` slot of `x`.
#'
#' @return An object of the same class as `x`, with the `data` slot
#' adjusted appropriately.
#'
#' @author Dan Kelley
#'
#' @examples
#'\dontrun{
#' library(oce)
#' rbr011855 <- read.oce(
#'  "/data/archive/sleiwex/2008/moorings/m08/pt/rbr_011855/raw/pt_rbr_011855.dat")
#' d <- subset(rbr011855, time < as.POSIXct("2008-06-25 10:05:00"))
#' x <- undriftTime(d, 1)   # clock lost 1 second over whole experiment
#' summary(d)
#' summary(x)
#'}
undriftTime <- function(x, slowEnd = 0, tname="time")
{
    if (!inherits(x, "oce"))
        stop("method is only for oce objects")
    names <- names(x@data)
    if (!(tname %in% names))
        stop("no column named '", tname, "'; only found: ", paste(names, collapse=" "))
    res <- x
    time <- res@data[[tname]]
    nt <- length(time)
    if (nt < 2) {
        warning("too few data to to undrift time; returning object unaltered")
    } else {
        sampleInterval <- as.numeric(difftime(time[2], time[1], units="s"))
        nt <- length(time)
        nt.out <- floor(0.5 + nt + slowEnd / sampleInterval)
        time.out <- seq.POSIXt(from=time[1], by=sampleInterval, length.out=nt.out)
        i <- seq(from=1, by=1, length.out=nt)
        i.out <- seq(from=1, to=nt, length.out = nt.out)
        out <- data.frame(array(dim=c(nt.out, length(x@data))))
        names(out) <- names
        out[[tname]] <- time.out
        for (name in names) {
            if (name != tname) {
                yy <- approx(x=i, y=x@data[[name]], xout=i.out)$y
                out[[name]] <- yy
            }
        }
        res@data <- out
    }
    res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep="", collapse=""))
    res
}


#' Fill a gap in an oce object
#'
#' Sequences of `NA` values, are filled by linear interpolation between
#' the non-`NA` values that bound the gap.
#'
#' @param x an [oce-class] object.
#'
#' @param method to use; see \dQuote{Details}.
#'
#' @param rule integer controlling behaviour at start and end of `x`.  If
#' `rule=1`, `NA` values at the ends are left in the return value.
#' If `rule=2`, they are replaced with the nearest non-NA point.
#'
#' @return A new `oce` object, with gaps removed.
#'
#' @section Bugs:
#' 1. Eventually, this will be expanded to work
#' with any `oce` object.  But, for now, it only works for vectors that
#' can be coerced to numeric.
#' 2. If the first or last point is `NA`, then `x` is returned unaltered.
#' 3. Only method `linear` is permitted now.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' # Integers
#' x <- c(1:2, NA, NA, 5:6)
#' y <- fillGap(x)
#' print(data.frame(x, y))
#' # Floats
#' x <- x + 0.1
#' y <- fillGap(x)
#' print(data.frame(x, y))
fillGap <- function(x, method=c("linear"), rule=1)
{
    if (!is.numeric(x))
        stop("only works for numeric 'x'")
    method <- match.arg(method)
    class <- class(x)
    if (is.vector(x)) {
        #res <- .Call("fillgap1d", as.numeric(x), rule)
        res <- do_fill_gap_1d(x, rule)
    } else if (is.matrix(x))  {
        res <- x
        for (col in seq_len(ncol(x))) {
            res[, col] <- do_fill_gap_1d(x[, col], rule)
        }
        for (row in seq_len(nrow(x))) {
            res[row, ] <- do_fill_gap_1d(x[row, ], rule)
        }
    } else {
        stop("only works if 'x' is a vector or a matrix")
    }
    class(res) <-  class
    res
}


#' Smooth and Decimate, or Subsample, an Oce Object
#'
#' Later on, other methods will be added, and [ctdDecimate()] will be
#' retired in favour of this, a more general, function.  The filtering is done
#' with the [filter()] function of the stats package.
#'
#' @param x an [oce-class] object.
#'
#' @param by an indication of the subsampling.  If this is a single number,
#' then it indicates the spacing between elements of `x` that are
#' selected.  If it is two numbers (a condition only applicable if `x` is
#' an `echosounder` object, at present), then the first number indicates
#' the time spacing and the second indicates the depth spacing.
#'
#' @param to Indices at which to subsample.  If given, this over-rides
#' `by`.
#'
#' @param filter optional list of numbers representing a digital filter to be
#' applied to each variable in the `data` slot of `x`, before
#' decimation is done. If not supplied, then the decimation is done strictly by
#' sub-sampling.
#'
#' @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.
#'
#' @return An [oce-class] object that has been subsampled appropriately.
#'
#' @section Bugs: Only a preliminary version of this function is provided in
#' the present package.  It only works for objects of class `echosounder`,
#' for which the decimation is done after applying a running median filter and
#' then a boxcar filter, each of length equal to the corresponding component of
#' `by`.
#'
#' @author Dan Kelley
#'
#' @seealso Filter coefficients may be calculated using
#' [makeFilter()].  (Note that [ctdDecimate()] will be
#' retired when the present function gains equivalent functionality.)
#'
#' @examples
#' library(oce)
#' data(adp)
#' plot(adp)
#' adpDec <- decimate(adp, by=2, filter=c(1/4, 1/2, 1/4))
#' plot(adpDec)
decimate <- function(x, by=10, to, filter, debug=getOption("oceDebug"))
{
    if (!inherits(x, "oce"))
        stop("method is only for oce objects")
    oceDebug(debug, "in decimate(x, by=", by, ", to=", if (missing(to)) "unspecified" else to, "...)\n")
    res <- x
    do.filter <- !missing(filter)
    if ("time" %in% names(x@data)) {
        if (missing(to))
            to <- length(x@data$time[[1]])
        if (length(by) == 1) {
            # FIXME: probably should not be here
            select <- seq(from=1, to=to, by=by)
            oceDebug(debug, vectorShow(select, "select:"))
        }
    }
    if (inherits(x, "adp")) {
        oceDebug(debug, "decimate() on an ADP object\n")
        warning("decimate(adp) not working yet ... just returning the adp unchanged")
        return(res) # FIXME
        #nbeam <- dim(x@data$v)[3]
        for (name in names(x@data)) {
            oceDebug(debug, "decimating item named '", name, "'\n")
            if ("distance" == name)
                next
            if ("time" == name) {
                res@data[[name]] <- x@data[[name]][select]
            } else if (is.vector(x@data[[name]])) {
                oceDebug(debug, "subsetting x@data$", name, ", which is a vector\n", sep="")
                if (do.filter)
                    res@data[[name]] <- filterSomething(x@data[[name]], filter)
                res@data[[name]] <- res@data[[name]][select]
            } else if (is.matrix(x@data[[name]])) {
                dim <- dim(x@data[[name]])
                for (j in 1: dim[2]) {
                    oceDebug(debug, "subsetting x@data[[", name, ",", j, "]], which is a matrix\n", sep="")
                    if (do.filter)
                        res@data[[name]][, j] <- filterSomething(x@data[[name]][, j], filter)
                    res@data[[name]][, j] <- res@data[[name]][, j][select]
                }
            } else if (is.array(x@data[[name]])) {
                dim <- dim(x@data[[name]])
                #print(dim)
                for (k in seq_len(dim[2])) {
                    for (j in seq_len(dim[3])) {
                        oceDebug(debug, "subsetting x@data[[", name, "]][", j, ",", k, "], which is an array\n", sep="")
                        if (do.filter)
                            res@data[[name]][, j, k] <- filterSomething(x@data[[name]][, j, k], filter)
                        res@data[[name]][, j, k] <- res@data[[name]][, j, k][select]
                    }
                }
            }
        }
    } else if (inherits(x, "adv")) {
        # FIXME: the (newer) adp code is probably better than this ADV code
        oceDebug(debug, "decimate() on an ADV object\n")
        warning("decimate(adv) not working yet ... just returning the adv unchanged")
        return(res) # FIXME
        for (name in names(x@data)) {
            if ("time" == name) {
                res@data[[name]] <- x@data[[name]][select]
            } else if (is.vector(x@data[[name]])) {
                oceDebug(debug, "decimating x@data$", name, ", which is a vector\n", sep="")
                if (do.filter)
                    res@data[[name]] <- filterSomething(x@data[[name]], filter)
                res@data[[name]] <- res@data[[name]][select]
            } else if (is.matrix(x@data[[name]])) {
                dim <- dim(x@data[[name]])
                for (j in 1: dim[2]) {
                    oceDebug(debug, "decimating x@data[[", name, ",", j, "]], which is a matrix\n", sep="")
                    if (do.filter)
                        res@data[[name]][, j] <- filterSomething(x@data[[name]][, j], filter)
                    res@data[[name]][, j] <- res@data[[name]][, j][select]
                }
            } 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, "decimating x@data[[", name, ",", j, ",", k, "]], which is an array\n", sep="")
                        if (do.filter)
                            res@data[[name]][, j, k] <- filterSomething(x@data[[name]][, j, k], filter)
                        res@data[[name]][, j, k] <- res@data[[name]][, j, k][select]
                    }
                }
            } else {
                stop("item data[[", name, "]] is not understood; it must be a vector, a matrix, or an array")
            }
        }
    } else if (inherits(x, "ctd")) {
        warning("decimate(ctd) not working yet ... just returning the ctd unchanged")
        return(res) # FIXME
        if (do.filter)
            stop("cannot (yet) filter ctd data during decimation") # FIXME
        select <- seq(1, dim(x@data)[1], by=by)
        res@data <- x@data[select, ]
    } else if (inherits(x, "pt")) {
        warning("decimate(pt) not working yet ... just returning the pt unchanged")
        return(res) # FIXME
        if (do.filter)
            stop("cannot (yet) filter pt data during decimation") # FIXME
        for (name in names(res@data))
            res@data[[name]] <- x@data[[name]][select]
    } else if (inherits(x, "echosounder")) {
        oceDebug(debug, "decimate() on an 'echosounder' object\n")
        # use 'by', ignoring 'to' and filter'
        if (length(by) != 2)
            stop("length(by) must equal 2.  First element is width of boxcar in pings, second is width in depths")
        by <- as.integer(by)
        byPing <- by[1]
        kPing <- as.integer(by[1])
        if (0 == kPing%%2)
            kPing <- kPing + 1
        byDepth <- by[2]
        kDepth <- as.integer(by[2])
        if (0 == kDepth%%2)
            kDepth <- kDepth + 1
        if (byDepth > 1) {
            depth <- x[["depth"]]
            a <- x[["a"]]
            ncol <- ncol(a)
            nrow <- nrow(a)
            ii <- 1:ncol
            depth2 <- binAverage(ii, depth, 1, ncol, byDepth)$y
            a2 <- matrix(nrow=nrow(a), ncol=length(depth2))
            for (r in 1:nrow) {
                a2[r, ] <- binAverage(ii, runmed(a[r, ], kDepth), 1, ncol, byDepth)$y
            }
            res <- x
            res[["depth"]] <- depth2
            res[["a"]] <- a2
            x <- res # need for next step
        }
        if (byPing > 1) {
            #time <- x[["time"]]
            a <- x[["a"]]
            ncol <- ncol(a)
            nrow <- nrow(a)
            jj <- 1:nrow
            time2 <- binAverage(jj, as.numeric(x[["time"]]), 1, nrow, byPing)$y + as.POSIXct("1970-01-01 00:00:00", tz="UTC")
            a2 <- matrix(nrow=length(time2), ncol=ncol(a))
            for (c in 1:ncol) {
                a2[, c] <- binAverage(jj, runmed(a[, c], kPing), 1, nrow, byPing)$y
            }
            res <- x
            res[["time"]] <- time2
            res[["latitude"]] <- binAverage(jj, x[["latitude"]], 1, nrow, byPing)$y
            res[["longitude"]] <- binAverage(jj, x[["longitude"]], 1, nrow, byPing)$y
            res[["a"]] <- a2
        }
        # do depth, rows of matrix, time, cols of matrix
    } else if (inherits(x, "topo")) {
        oceDebug(debug, "Decimating a topo object")
        lonlook <- seq(1, length(x[["longitude"]]), by=by)
        latlook <- seq(1, length(x[["latitude"]]), by=by)
        res[["longitude"]] <- x[["longitude"]][lonlook]
        res[["latitude"]] <- x[["latitude"]][latlook]
        res[["z"]] <- x[["z"]][lonlook, latlook]
    } else if (inherits(x, "landsat")) {
        oceDebug(debug, "Decimating a landsat object with by=", by, "\n")
        for (i in seq_along(x@data)) {
            b <- x@data[[i]]
            if (is.list(b)) {
                dim <- dim(b$msb)
                if (!is.null(dim)) {
                    res@data[[i]]$msb <- b$msb[seq(1, dim[1], by=by), seq(1, dim[2], by=by)]
                }
                dim <- dim(b$lsb)
                res@data[[i]]$lsb <- b$lsb[seq(1, dim[1], by=by), seq(1, dim[2], by=by)]
            } else {
                dim <- dim(x@data[[i]])
                res@data[[i]] <- b[seq(1, dim[1], by=by), seq(1, dim[2], by=by)]
            }
        }
    } else {
        stop("decimation does not work (yet) for objects of class ", paste(class(x), collapse=" "))
    }
    if ("deltat" %in% names(x@metadata)) # FIXME: should handle for individual cases, not here
        res@metadata$deltat <- by * x@metadata$deltat
    res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep="", collapse=""))
    res
}



#' 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



#' Decode BCD to integer
#'
#' @param x a raw value, or vector of raw values, coded in binary-coded
#' decimal.
#'
#' @param endian character string indicating the endian-ness ("big" or
#' "little").  The PC/intel convention is to use "little", and so most data
#' files are in that format.
#'
#' @return An integer, or list of integers.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' twenty.five <- bcdToInteger(as.raw(0x25))
#' thirty.seven <- as.integer(as.raw(0x25))
bcdToInteger <- function(x, endian=c("little", "big"))
{
    endian <- match.arg(endian)
    x <- as.integer(x)
    byte1 <- as.integer(floor(x / 16))
    byte2 <- x - 16 * byte1
    if (endian=="little") 10*byte1 + byte2 else byte1 + 10*byte2
}


#' Format bytes as binary [defunct]
#'
#' **WARNING:** The `endian` argument will soon be removed
#' from this function; see \link{oce-defunct}.
#' This is because the actions for `endian="little"` made
#' no sense in practical work. The default value for `endian`
#' was changed to `"big"` on 2017 May 6.
#'
#' @param x an integer to be interpreted as a byte.
#'
#' @param endian character string indicating the endian-ness ("big" or
#' "little"). **WARNING:** This argument will be removed soon.
#'
#' @return A character string representing the bit strings for the elements of
#' `x`, in order of significance for the `endian="big"` case.
#' (The nibbles, or 4-bit sequences, are interchanged in the now-deprecated
#' `"little"` case.)
#' See \dQuote{Examples} for how this relates to the output from
#' \link{rawToBits}.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' # Note comparison with rawToBits():
#' a <- as.raw(0x0a)
#' byteToBinary(a, "big")        # "00001010"
#' as.integer(rev(rawToBits(a))) # 0 0 0 0 1 0 1 0
byteToBinary <- function(x, endian="big")
{
    if (endian != "big") {
        .Defunct("rawToBits",
                 msg="byteToBinary(.,'little') is disallowed and will be removed soon. See ?'oce-defunct'.")
    }
    # onebyte2binary <- function(x)
    # {
    #     c("0000", "0001", "0010", "0011",
    #       "0100", "0101", "0110", "0111",
    #       "1000", "1001", "1010", "1011",
    #       "1100", "1101", "1110", "1111")[x+1]
    # }
    # res <- NULL
    # if (class(x) == "raw")
    #     x <- readBin(x, "int", n=length(x), size=1, signed=FALSE)
    # for (i in seq_along(x)) {
    #     if (x[i] < 0) {
    #         res <- c(res, "??")
    #     } else {
    #         # FIXME: these are not bytes here; they are nibbles.  I don't think endian="little"
    #         # makes ANY SENSE at all.  2012-11-22
    #         byte1 <- as.integer(floor(x[i] / 16))
    #         byte2 <- x[i] - 16 * byte1
    #         #cat("input=",x[i],"byte1=",byte1,"byte2=",byte2,"\n")
    #         if (endian == "little")
    #             res <- c(res, paste(onebyte2binary(byte2), onebyte2binary(byte1), sep=""))
    #         else
    #             res <- c(res, paste(onebyte2binary(byte1), onebyte2binary(byte2), sep=""))
    #         #cat(" res=",res,"\n")
    #     }
    # }
    # res
    x <- as.raw(x)
    paste(ifelse(rev(rawToBits(x)==as.raw(0x01)), "1", "0"), collapse="")
}


#' Confidence interval in parenthetic notation
#'
#' Format a confidence interval in parenthetic notation.
#'
#' If a `model` is given, then `ci` is ignored, and a confidence
#' interval is calculated using [confint()] with `level` set to
#' 0.6914619.  This `level` corresponds to a range of plus or minus one
#' standard deviation, for the t distribution and a large number of degrees of
#' freedom (since `qt(0.6914619, 100000)` is 0.5).
#'
#' If `model` is missing, `ci` must be provided.  If it contains 3
#' elements, then first and third elements are taken as the range of the
#' confidence interval (which by convention should use the `level` stated
#' in the previous paragraph), and the second element is taken as the central
#' value.  Alternatively, if `ci` has 2 elements, they are taken to be
#' bounds of the confidence interval and their mean is taken to be the central
#' value.
#'
#' In the `+/-` notation, e.g. \eqn{a \pm b}{a +/- b} means that the true
#' value lies between \eqn{a-b}{a-b} and \eqn{a+b}{a+b} with a high degree of
#' certainty.  Mills et al. (1993, section 4.1 on page 83) suggest that
#' \eqn{b}{b} should be set equal to 2 times the standard uncertainty or
#' standard deviation.  JCGM (2008, section 7.2.2 on pages 25 and 26), however,
#' suggest that \eqn{b}{b} should be set to the standard uncertainty, while
#' also recommending that the \eqn{\pm}{+/-} notation be avoided altogether.
#'
#' The `parentheses` notation is often called the compact notation.  In
#' it, the digits in parentheses indicate the uncertainty in the corresponding
#' digits to their left, e.g. 12.34(3) means that the last digit (4) has an
#' uncertainty of 3.  However, as with the \eqn{\pm}{+/-} notation, different
#' authorities offer different advice on defining this uncertainty; Mills et
#' al. (1993, section 4.1 on page 83) provide an example in which the
#' parenthetic notation has the same value as the \eqn{\pm}{+/-} notation,
#' while JCM (2008, section 7.2.2 on pages 25 and 26) suggest halving the
#' number put in parentheses.
#'
#' The `foramtci` function is based on the JCM (2008) notation, i.e.
#' `formatCI(ci=c(8,12), style="+/-")` yields `"10+-2"`, and
#' `formatCI(ci=c(8,12), style="parentheses")` yields `"10(2)"`.
#'
#' **Note:** if the confidence range exceeds the value, the
#' `parentheses` format reverts to `+/-` format.
#'
#' @param ci optional vector of length 2 or 3.
#'
#' @param style string indicating notation to be used.
#'
#' @param model optional regression model, e.g. returned by [lm()] or
#' [nls()].
#'
#' @param digits optional number of digits to use; if not supplied,
#' [`getOption`[`("digits")` is used.
#'
#' @return If `ci` is given, the result is a character string with the
#' estimate and its uncertainty, in plus/minus or parenthetic notation.  If
#' `model` is given, the result is a 1-column matrix holding character
#' strings, with row names corresponding to the parameters of the model.
#'
#' @author Dan Kelley
#'
#' @references JCGM, 2008. *Evaluation of measurement data - Guide to the
#' expression of uncertainty in measurement (JCGM 100:2008)*, published by the
#' Joint Committee for Guides in Metrology,
#' http://www.bipm.org/en/publications/guides/gum.html (see section
#' 7.2.2 for a summary of notation, which shows equal values to the right of a
#' `+-` sign and in parentheses.
#'
#' I. Mills, T. Cvitas, K. Homann, N. Kallay, and K. Kuchitsu, 1993.
#' *Quantities, Units and Symbols in Physical Chemistry*, published
#' Blackwell Science for the International Union of Pure and Applied Chemistry.
#' (See section 4.1, page 83, for a summary of notation, which shows that a
#' value to the right of a `+-` sign is to be halved if put in
#'
#' @examples
#' x <- seq(0, 1, length.out=300)
#' y <- rnorm(n=300, mean=10, sd=1) * x
#' m <- lm(y~x)
#' print(formatCI(model=m))
formatCI <- function(ci, style=c("+/-", "parentheses"), model, digits=NULL)
{
    formatCI.one <- function(ci, style, digits=NULL)
    {
        debug <- FALSE
        if (missing(ci))
            stop("must supply ci")
        ci <- as.numeric(ci)
        if (length(ci) == 3) {
            x <- ci[2]
            ci <- ci[c(1, 3)]
        } else if (length(ci) == 2) {
            x <- mean(ci)
        } else {
            stop("ci must contain 2 or 3 elements")
        }
        sign <- sign(x)
        x <- abs(x)
        if (style == "+/-") {
            pm <- abs(diff(ci)/2)
            if (is.null(digits)) {
                paste(format(sign * x, digits=getOption("digits")), "+/-", format(pm, digits=getOption("digits")), sep="")
            } else {
                paste(format(sign * x, digits=digits), "+/-", format(pm, digits=digits), sep="")
            }
        } else {
            pm <- abs(diff(ci)/2)
            scale <- 10^floor(log10(pm))
            pmr <- round(pm / scale)
            if (pmr == 10) {
                pmr <- 1
                scale <- scale * 10
            }
            #scale <- 10^floor(log10(x))
            #x0 <- x / scale
            #ci0 <- ci / scale
            if (pm > x)
                return(paste(sign*x, "+/-", pm, sep=""))
            digits <- floor(log10(scale) + 0.1)
            if (digits < 0) {
                fmt <- paste("%.", abs(digits), "f", sep="")
            } else {
                fmt <- "%.f"
            }
            oceDebug(debug, "pm=", pm, ";pmr=", pmr, "; scale=", scale, "pm/scale=",
                pm/scale, "round(pm/scale)=", round(pm/scale), "\n", " x=", x,
                "; x/scale=", x/scale, "digits=", digits, "fmt=", fmt, "\n")
            paste(sprintf(fmt, sign*x), "(", pmr, ")", sep="")
        }
    }
    style <- match.arg(style)
    if (!missing(model)) {
        cm <- class(model)
        # > qt(0.6914619, 100000)
        # [1] 0.5
        if (cm == "lm" || cm == "nls") {
            ci <- confint(model, level=0.6914619)
            names <- dimnames(ci)[[1]]
            n <- length(names)
            res <- matrix("character", nrow=n, ncol=1)
            rownames(res) <- names
            colnames(res) <- "value"
            for (row in seq_len(dim(ci)[1])) {
                res[row, 1] <- formatCI.one(ci=ci[row, ], style=style, digits=digits)
            }
        }
        res
    } else {
        if (missing(ci)) {
            stop("must give either ci or model")
        }
        formatCI.one(ci=ci, style=style, digits=digits)
    }
}


#' Decode integer to corresponding ASCII code
#'
#' @param i an integer, or integer vector.
#'
#' @return A character, or character vector.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' A <- integerToAscii(65)
#' cat("A=", A, "\n")
integerToAscii <- function(i)
{
    c("", "\001", "\002", "\003", "\004", "\005", "\006", "\a", "\b",
      "\t", "\n", "\v", "\f", "\r", "\016", "\017", "\020", "\021",
      "\022", "\023", "\024", "\025", "\026", "\027", "\030", "\031",
      "\032", "\033", "\034", "\035", "\036", "\037", " ", "!", "\"",
      "#", "$", "%", "&", "'", "(", ")", "*", "+", ",", "-", ".", "/",
      "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", ":", ";", "<",
      "=", ">", "?", "@", "A", "B", "C", "D", "E", "F", "G", "H", "I",
      "J", "K", "L", "M", "N", "O", "P", "Q", "R", "S", "T", "U", "V",
      "W", "X", "Y", "Z", "[", "\\", "]", "^", "_", "`", "a", "b",
      "c", "d", "e", "f", "g", "h", "i", "j", "k", "l", "m", "n", "o",
      "p", "q", "r", "s", "t", "u", "v", "w", "x", "y", "z", "{", "|",
      "}", "~", "\177", "\x80", "\x81", "\x82", "\x83", "\x84", "\x85",
      "\x86", "\x87", "\x88", "\x89", "\x8a", "\x8b", "\x8c", "\x8d",
      "\x8e", "\x8f", "\x90", "\x91", "\x92", "\x93", "\x94", "\x95",
      "\x96", "\x97", "\x98", "\x99", "\x9a", "\x9b", "\x9c", "\x9d",
      "\x9e", "\x9f", "\xa0", "\xa1", "\xa2", "\xa3", "\xa4", "\xa5",
      "\xa6", "\xa7", "\xa8", "\xa9", "\xaa", "\xab", "\xac", "\xad",
      "\xae", "\xaf", "\xb0", "\xb1", "\xb2", "\xb3", "\xb4", "\xb5",
      "\xb6", "\xb7", "\xb8", "\xb9", "\xba", "\xbb", "\xbc", "\xbd",
      "\xbe", "\xbf", "\xc0", "\xc1", "\xc2", "\xc3", "\xc4", "\xc5",
      "\xc6", "\xc7", "\xc8", "\xc9", "\xca", "\xcb", "\xcc", "\xcd",
      "\xce", "\xcf", "\xd0", "\xd1", "\xd2", "\xd3", "\xd4", "\xd5",
      "\xd6", "\xd7", "\xd8", "\xd9", "\xda", "\xdb", "\xdc", "\xdd",
      "\xde", "\xdf", "\xe0", "\xe1", "\xe2", "\xe3", "\xe4", "\xe5",
      "\xe6", "\xe7", "\xe8", "\xe9", "\xea", "\xeb", "\xec", "\xed",
      "\xee", "\xef", "\xf0", "\xf1", "\xf2", "\xf3", "\xf4", "\xf5",
      "\xf6", "\xf7", "\xf8", "\xf9", "\xfa", "\xfb", "\xfc", "\xfd",
      "\xfe", "\xff")[i+1]
}


#' Earth magnetic declination, inclination, and intensity
#'
#' Implements the 12th and 13th generations of the
#' International Geomagnetic Reference Field
#' (IGRF), based on a reworked version of a Fortran program downloaded from a
#' NOAA website (see reference 1).
#'
#' The code (subroutines `igrf12syn` and `igrf13syn`) seem to have
#' been written by Susan Macmillan of the British Geological Survey.  Comments
#' in the source code `igrf13syn` (the current default used here)
#' indicate that its coefficients were agreed to in
#' December 2019 by the IAGA Working Group V-MOD.  Other comments in that code
#' suggest that the proposed application time interval is from years 1900 to 2025, inclusive,
#' but that only dates from 1945 to 2015 are to be considered definitive.
#'
#' @param longitude longitude in degrees east (negative for degrees west), as a
#' number, a vector, or a matrix.
#'
#' @param latitude latitude in degrees north, as a number, vector, or matrix.
#' The shape (length or dimensions) must conform to the dimensions of `longitude`.
#'
#' @param time The time at which the field is desired. This may be a
#' single value or a vector or matrix that is structured to match
#' `longitude` and `latitude`. The value may a decimal year,
#' a POSIXt time, or a Date time.
#'
#' @param version an integer that must be either 12 or 13, to specify
#' the version number of the formulae. Note that 13 became the default
#' on 2020 March 3, so to old code will need to specify `version=12`
#' to work as it did before that date.
#'
#' @return A list containing `declination`, `inclination`, and
#' `intensity`.
#'
#' @author Dan Kelley wrote the R code and a fortran wrapper to the
#' `igrf12.f` subroutine, which was written by Susan Macmillan of the
#' British Geological Survey and distributed ``without limitation'' (email from
#' SM to DK dated June 5, 2015).
#'
#' @section Historical Notes:
#' For about a decade, `magneticField` used the version 12 formulae provided
#' by IAGA, but the code was updated on March 3, 2020, to version 13.  Example
#' 3 shows that the differences in declination are typically under 2 degrees
#' (with 95 percent of the data lying between -1.7 and 0.7 degrees).
#'
#' @references
#' 1. The underlying Fortran code for version 12 is from `igrf12.f`, downloaded the NOAA
#' website (`https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html`) on June 7,
#' 2015. That for version 13 is `igrf13.f`, downloaded from the NOAA website
#' (`https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html` on March 3, 2020.
#' 2. Witze, Alexandra. \dQuote{Earth's Magnetic Field Is Acting up and Geologists Don't Know Why.}
#' Nature 565 (January 9, 2019): 143.
#' \doi{10.1038/d41586-019-00007-1}
#'
#' @examples
#' library(oce)
#' # 1. Today's value at Halifax NS
#' magneticField(-(63+36/60), 44+39/60, Sys.Date())
#'
#' # 2. World map of declination in year 2000.
#'\donttest{
#' data(coastlineWorld)
#' par(mar=rep(0.5, 4)) # no axes on whole-world projection
#' mapPlot(coastlineWorld, projection="+proj=robin", col="lightgray")
#' # Construct matrix holding declination
#' lon <- seq(-180, 180)
#' lat <- seq(-90, 90)
#' dec2000 <- function(lon, lat) {
#'     magneticField(lon, lat, 2000)$declination
#' }
#' dec <- outer(lon, lat, dec2000) # hint: outer() is very handy!
#' # Contour, unlabelled for small increments, labeled for
#' # larger increments.
#' mapContour(lon, lat, dec, col="blue", levels=seq(-180, -5, 5),
#'            lty=3, drawlabels=FALSE)
#' mapContour(lon, lat, dec, col="blue", levels=seq(-180, -20, 20))
#' mapContour(lon, lat, dec, col="red", levels=seq(5, 180, 5),
#'            lty=3, drawlabels=FALSE)
#' mapContour(lon, lat, dec, col="red", levels=seq(20, 180, 20))
#' mapContour(lon, lat, dec, levels=180, col="black", lwd=2, drawlabels=FALSE)
#' mapContour(lon, lat, dec, levels=0, col="black", lwd=2)
#'}
#'
#' # 3. Declination differences between versions 12 and 13
#'\donttest{
#' lon <- seq(-180, 180)
#' lat <- seq(-90, 90)
#' decDiff <- function(lon, lat) {
#'     old <- magneticField(lon, lat, 2020, version=13)$declination
#'     new <- magneticField(lon, lat, 2020, version=12)$declination
#'     new - old
#' }
#' decDiff <- outer(lon, lat, decDiff)
#' decDiff <- ifelse(decDiff > 180, decDiff - 360, decDiff)
#' # Overall (mean) shift -0.1deg
#' t.test(decDiff)
#' # View histogram, narrowed to small differences
#' par(mar=c(3.5, 3.5, 2, 2), mgp=c(2, 0.7, 0))
#' hist(decDiff, breaks=seq(-180, 180, 0.05), xlim=c(-2, 2),
#'      xlab="Declination difference [deg] from version=12 to version=13",
#'      main="Predictions for year 2020")
#' print(quantile(decDiff, c(0.025, 0.975)))
#' # Note that the large differences are at high latitudes
#' imagep(lon,lat,decDiff, zlim=c(-1,1)*max(abs(decDiff)))
#' lines(coastlineWorld[["longitude"]], coastlineWorld[["latitude"]])
#'}
#' @family things related to magnetism
magneticField <- function(longitude, latitude, time, version=13)
{
    if (missing(longitude) || missing(latitude) || missing(time))
        stop("must provide longitude, latitude, and time")
    dim <- dim(latitude)
    if (!all(dim == dim(longitude)))
        stop("dimensions of longitude and latitude must agree")
    n <- length(latitude)
    if (inherits(time, "Date"))
        time <- as.POSIXct(time, tz="UTC")
    if (inherits(time, "POSIXt")) {
        d <- as.POSIXlt(time, tz="UTC")
        year <- d$year+1900
        yearday <- d$yday
        time <- year + yearday / 365.25 # ignore leap year issue (formulae not daily)
    }
    if (length(time) == 1) {
        time <- rep(time, n)
    } else {
        if (!all(dim == dim(time)))
            stop("dimensions of latitude and time must agree")
    }
    if (!is.null(dim)) {
        dim(longitude) <- n
        dim(latitude) <- n
        dim(time) <- n
    }
    #isv <- 0
    #itype <- 1                          # geodetic
    #alt <- 0.0                          # altitude in km
    elong <- ifelse(longitude < 0, 360 + longitude, longitude)
    colat <- 90 - latitude
    iversion <- as.integer(version)
    if (!(iversion %in% c(12L, 13L))) {
        stop("version must be 12 or 13, but it is ", iversion)
    }
    # message("time:", time, " (", as.numeric(time), ")")
    r <- .Fortran("md_driver",
        as.double(colat), as.double(elong), as.double(time),
        as.integer(n),
        declination=double(n),
        inclination=double(n),
        intensity=double(n),
        as.integer(iversion))
    declination <- r$declination
    inclination <- r$inclination
    intensity <- r$intensity
    if (!is.null(dim)) {
        dim(declination) <- dim
        dim(inclination) <- dim
        dim(intensity) <- dim
    }
    list(declination=declination, inclination=inclination, intensity=intensity)
}


#' Locate byte sequences in a raw vector
#'
#' Find spots in a raw vector that match a given byte sequence.
#'
#' @param input a vector of raw (byte) values.
#'
#' @param b1 a vector of bytes to match (must be of length 2 or 3 at present;
#' for 1-byte, use [which()]).
#'
#' @param \dots additional bytes to match for (up to 2 permitted)
#'
#' @return List of the indices of `input` that match the start of the
#' `bytes` sequence (see example).
#'
#' @author Dan Kelley
#'
#' @examples
#' buf <- as.raw(c(0xa5, 0x11, 0xaa, 0xa5, 0x11, 0x00))
#' match <- matchBytes(buf, 0xa5, 0x11)
#' print(buf)
#' print(match)
matchBytes <- function(input, b1, ...)
{
    if (missing(input))
        stop("must provide \"input\"")
    if (missing(b1))
        stop("must provide at least one byte to match")
    #n <- length(input)
    dots <- list(...)
    lb <- 1 + length(dots)
    if (lb == 2) {
        .Call("match2bytes", as.raw(input), as.raw(b1), as.raw(dots[[1]]), FALSE)
    } else if (lb == 3) {
        .Call("match3bytes", as.raw(input), as.raw(b1), as.raw(dots[[1]]), as.raw(dots[[2]]))
    } else {
        stop("must provide 2 or 3 bytes")
    }
}


#' Rearrange areal matrix so Greenwich is near the centre
#'
#' Sometimes datasets are provided in matrix form, with first
#' index corresponding to longitudes ranging from 0 to 360.
#' `matrixShiftLongitude` cuts such matrices at
#' longitude=180, and swaps the pieces so that the dateline
#' is at the left of the matrix, not in the middle.
#'
#' @param m The matrix to be modified.
#'
#' @param longitude A vector containing the longitude in the 0-360 convention. If missing,
#' this is constructed to range from 0 to 360, with as many elements as the first index of `m`.
#'
#' @return A list containing `m` and `longitude`, both rearranged as appropriate.
#'
#' @seealso [shiftLongitude()] and [standardizeLongitude()].
matrixShiftLongitude <- function(m, longitude)
{
    if (missing(m))
        stop("must supply m")
    n <- dim(m)[1]
    if (missing(longitude))
        longitude <- seq.int(0, 360, length.out=n)
    if (n != length(longitude))
        stop("dim(m) and length(longitude) are incompatible")
    if (max(longitude, na.rm=TRUE) > 180) {
        cut <- which.min(abs(longitude-180))
        longitude <- c(longitude[seq.int(cut+1L, n)]-360, longitude[seq.int(1L, cut)])
        m <- m[c(seq.int(cut+1L, n), seq.int(1L, cut)), ]
    }
    list(m=m, longitude=longitude)
}


#' Smooth a Matrix
#'
#' The values on the edge of the matrix are unaltered.  For interior points,
#' the result is defined in terms in terms of the original as follows.
#' \eqn{r_{i,j} = (2 m_{i,j} + m_{i-1,j} + m_{i+1,j} + m_{i,j-1} +
#' m_{i,j+1})/6}{r_[i,j] = (2 m_[i,j] + m_[i-1,j] + m_[i+1,j] + m_[i,j-1] +
#' m_[i,j+1])/6}.  Note that missing values propagate to neighbours.
#'
#' @param m a matrix to be smoothed.
#'
#' @param passes an integer specifying the number of times the smoothing is to
#' be applied.
#'
#' @return A smoothed matrix.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' opar <- par(no.readonly = TRUE)
#' m <- matrix(rep(seq(0, 1, length.out=5), 5), nrow=5, byrow=TRUE)
#' m[3, 3] <- 2
#' m1 <- matrixSmooth(m)
#' m2 <- matrixSmooth(m1)
#' m3 <- matrixSmooth(m2)
#' par(mfrow=c(2, 2))
#' image(m,  col=rainbow(100), zlim=c(0, 4), main="original image")
#' image(m1, col=rainbow(100), zlim=c(0, 4), main="smoothed 1 time")
#' image(m2, col=rainbow(100), zlim=c(0, 4), main="smoothed 2 times")
#' image(m3, col=rainbow(100), zlim=c(0, 4), main="smoothed 3 times")
#' par(opar)
matrixSmooth <- function(m, passes=1)
{
    if (missing(m))
        stop("must provide matrix 'm'")
    storage.mode(m) <- "double"
    if (passes > 0) {
        for (pass in seq.int(1, passes, 1)) {
            m <- do_matrix_smooth(m)
        }
    } else {
        warning("matrixSmooth given passes<=0, so returning matrix unmodified")
    }
    m
}


#' Time interval as colon-separated string
#'
#' Convert a time interval to a colon-separated string
#'
#' @param sec length of time interval in seconds.
#'
#' @return A string with a colon-separated time interval.
#'
#' @author Dan Kelley
#'
#' @seealso See [ctimeToSeconds()], the inverse of this.
#'
#' @examples
#' library(oce)
#' cat("   10 s = ", secondsToCtime(10), "\n", sep="")
#' cat("   61 s = ", secondsToCtime(61), "\n", sep="")
#' cat("86400 s = ", secondsToCtime(86400), "\n", sep="")
#' @family things related to time
secondsToCtime <- function(sec)
{
    if (sec < 60)
        return(sprintf("00:00:%02d", sec))
    if (sec < 3600) {
        min <- floor(sec / 60)
        sec <- sec - 60 * min
        return(sprintf("00:%02d:%02d", min, sec))
    }
    hour <- floor(sec / 3600)
    sec <- sec - 3600 * hour
    min <- floor(sec / 60)
    sec <- sec - 60 * min
    return(sprintf("%02d:%02d:%02d", hour, min, sec))
}


#' Interpret a character string as a time interval
#'
#' Interpret a character string as a time interval
#' Strings are of the form MM:SS or HH:MM:SS.
#'
#' @param ctime a character string (see \dQuote{Details}.
#'
#' @return A numeric value, the number of seconds represented by the string.
#'
#' @author Dan Kelley
#'
#' @seealso See [secondsToCtime()], the inverse of this.
#'
#' @examples
#' library(oce)
#' cat("10      = ", ctimeToSeconds("10"), "s\n", sep="")
#' cat("01:04   = ", ctimeToSeconds("01:04"), "s\n", sep="")
#' cat("1:00:00 = ", ctimeToSeconds("1:00:00"), "s\n", sep="")
#' @family things related to time
ctimeToSeconds <- function(ctime)
{
    if (length(grep(":", ctime)) > 0) {
        parts <- as.numeric(strsplit(ctime, ":")[[1]])
        l <- length(parts)
        if (l == 1) s <- as.numeric(ctime)
        else if (l == 2) s <- parts[1] * 60 + parts[2]
        else if (l == 3) s <- parts[1] * 3600 + parts[2] * 60 + parts[3]
        else stop("cannot interpret \"time\"=", ctime, "as a time interval because it has more than 2 colons")
    } else {
        s <- as.numeric(ctime)
    }
    s
}

#showThrees <- function(x, indent="    ")
#{
#    if (!("threes" %in% names(x)))
#        stop("'x' has no item named 'threes'")
#    rownames <- rownames(x$threes)
#    colnames <- colnames(x$threes)
#    data.width <- max(nchar(colnames)) + 10
#    name.width <- max(nchar(rownames(x$threes))) + 4 # space for left-hand column
#    ncol <- length(colnames)
#    nrow <- length(rownames)
#    res <- indent
#    res <- paste(res, format(" ", width=1+name.width), collapse="")
#    res <- paste(res, paste(format(colnames, width=data.width, justify="right"), collapse=" "))
#    res <- paste(res, "\n", sep="")
#    digits <- max(5, getOption("digits") - 1)
#    for (irow in 1L:nrow) {
#        res <- paste(res, indent, format(rownames[irow], width=name.width), "  ", sep="") # FIXME: should not need the "  "
#        for (icol in 1L:ncol) {
#            res <- paste(res, format(x$threes[irow,icol], digits=digits, width=data.width, justify="right"), sep=" ")
#        }
#        res <- paste(res, "\n", sep="")
#    }
#    res
#}



#' 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


#' Show metadata item
#'
#' This is a helper function for various `summary` functions.
#'
#' @param object an [oce-class] object.
#'
#' @param name name of item
#'
#' @param label label to print before item
#'
#' @param postlabel label to print after item
#'
#' @param isdate boolean indicating whether the item is a time
#'
#' @param quote boolean indicating whether to enclose the item in quotes
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' data(ctd)
#' showMetadataItem(ctd, "ship", "ship")
showMetadataItem <- function(object, name, label="", postlabel="", isdate=FALSE, quote=FALSE)
{
    if (name %in% names(object@metadata)) {
        item <- object@metadata[[name]]
        if (is.null(item))
            return()
        if (is.na(item))
            return()
        if (is.character(item) && nchar(item) == 0)
            return()
        if (is.na(item))
            return()
        if (isdate)
            item <- format(item)
        if (quote)
            item <- paste("\"", item, "\"", sep="")
        cat(paste("* ", label, item, postlabel, "\n", sep=""))
    }
}


#' Trapezoidal Integration
#'
#' Estimate the integral of one-dimensional function using the trapezoidal
#' rule.
#'
#' @param x,y vectors of x and y values. In the normal case, these
#' vectors are both supplied, and of equal length. There are also two
#' special cases. First, if `y` is missing, then
#' `x` is taken to be `y`, and a new `x` is constructed
#' as [`seq_along`]`(y)1. Second, if `length(x)` is 1
#' and `length(y)` exceeds 1, then `x` is replaced by
#' `x*`[`seq_along`]`(y)`.
#'
#' @param type Flag indicating the desired return value (see \dQuote{Value}).
#'
#' @param xmin,xmax Optional numbers indicating the range of the integration.
#' These values may be used to restrict the range of integration, or to
#' extend it; in either case, [approx()] with `rule=2`
#' is used to create new x and y vectors.
#'
#' @return If `type="A"` (the default), a single value is returned,
#' containing the estimate of the integral of `y=y(x)`.  If
#' `type="dA"`, a numeric vector of the same length as `x`, of which
#' the first element is zero, the second element is the integral between
#' `x[1]` and `x[2]`, etc.  If `type="cA"`, the result is the
#' cumulative sum (as in [cumsum()]) of the values that would be
#' returned for `type="dA"`.  See \dQuote{Examples}.
#'
#' @section Bugs: There is no handling of `NA` values.
#'
#' @author Dan Kelley
#'
#' @examples
#' x <- seq(0, 1, length.out=10) # try larger length.out to see if area approaches 2
#' y <- 2*x + 3*x^2
#' A <- integrateTrapezoid(x, y)
#' dA <- integrateTrapezoid(x, y, "dA")
#' cA <- integrateTrapezoid(x, y, "cA")
#' print(A)
#' print(sum(dA))
#' print(tail(cA, 1))
#' print(integrateTrapezoid(diff(x[1:2]), y))
#' print(integrateTrapezoid(y))
integrateTrapezoid <- function(x, y, type=c("A", "dA", "cA"), xmin, xmax)
{
    type <- match.arg(type)
    if (missing(x))
        stop("must supply 'x'")
    if (missing(y)) {
        y <- x
        x <- seq_along(y)
    }
    if (length(x) == 1 && length(y) > 1) {
       x <- x * seq_along(y)
    }
    if (length(x) != length(y)) {
        stop("'x' and 'y' must be of same length")
    }
    xout <- x
    yout <- y
    # message("x: ", paste(x, collapse=" "))
    # message("y: ", paste(y, collapse=" "))
    if (!missing(xmin) || !missing(xmax)) {
        if (missing(xmin) || missing(xmax))
            stop("'xmin' and 'xmax' must both be supplied, if either is")
        if (xmin >= xmax)
            stop("'xmin' must be less than 'xmax'")

        if (xmin > max(x, na.rm=TRUE))
            stop("xmin must be less than max(x)")
        if (xmin < min(x, na.rm=TRUE)) {
            xout <- c(xmin, xout)
        } else {
            xout <- xout[xout >= xmin]
            xout <- c(xmin, xout)
        }
        if (xmax < min(x, na.rm=TRUE)) {
            stop("xmax must be greater than min(x)")
        }
        if (xmax > max(x, na.rm=TRUE)) {
            xout <- c(xout, xmax)
        } else {
            xout <- xout[xout < xmax]
            xout <- c(xout, xmax)
        }
        yout <- approx(x, y, xout, rule=2)$y
    }
    # I think we should be able to use trap(), which gets defined into
    # R/RcppExports.R but that doesn't seem to be put into the loadspace.
    # My guess is that the problem is because we are not doing exports in
    # the recommended (automatic) way, but I don't want to do exports that
    # way since things are ok now, and have been for years.
    # I don't see much point trying to figure this out, because we already
    # have things set up for a .Call() from before the switch from C to Cpp.
    #
    # NOTE: must run Rcpp::compileAttributes() after creating trap in
    # src/trap.cpp
    res <- do_trap(xout, yout, as.integer(switch(match.arg(type), A=0, dA=1, cA=2)))
    res
}


#' Calculate Matrix Gradient
#'
#' In the interior of the matrix, centred second-order differences are used to
#' infer the components of the grad.  Along the edges, first-order differences
#' are used.
#'
#' @param h a matrix of values
#'
#' @param x vector of coordinates along matrix columns (defaults to integers)
#'
#' @param y vector of coordinates along matrix rows (defaults to integers)
#'
#' @return A list containing \eqn{|\nabla h|}{abs(grad(h))} as `g`,
#' \eqn{\partial h/\partial x}{dh/dx} as `gx`,
#' and \eqn{\partial h/\partial y}{dh/dy} as `gy`,
#' each of which is a matrix of the same dimension as `h`.
#'
#' @author Dan Kelley, based on advice of Clark Richards, and mimicking a matlab function.
#'
#' @examples
#' # 1. Built-in volcano dataset
#' g <- grad(volcano)
#' par(mfrow=c(2, 2), mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0))
#' imagep(volcano, zlab="h")
#' imagep(g$g, zlab="|grad(h)|")
#' zlim <- c(-1, 1) * max(g$g)
#' imagep(g$gx, zlab="dh/dx", zlim=zlim)
#' imagep(g$gy, zlab="dh/dy", zlim=zlim)
#'
#' # 2. Geostrophic flow around an eddy
#' library(oce)
#' dx <- 5e3
#' dy <- 10e3
#' x <- seq(-200e3, 200e3, dx)
#' y <- seq(-200e3, 200e3, dy)
#' R <- 100e3
#' h <- outer(x, y, function(x, y) 500*exp(-(x^2+y^2)/R^2))
#' grad <- grad(h, x, y)
#' par(mfrow=c(2, 2), mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0))
#' contour(x,y,h,asp=1, main=expression(h))
#' f <- 1e-4
#' gprime <- 9.8 * 1 / 1024
#' u <- -(gprime / f) * grad$gy
#' v <-  (gprime / f) * grad$gx
#' contour(x, y, u, asp=1, main=expression(u))
#' contour(x, y, v, asp=1, main=expression(v))
#' contour(x, y, sqrt(u^2+v^2), asp=1, main=expression(speed))
#'
#' @family things relating to vector calculus
grad <- function(h, x=seq(0, 1, length.out=nrow(h)), y=seq(0, 1, length.out=ncol(h)))
{
    if (missing(h))
        stop("must give h")
    if (length(x) != nrow(h))
        stop("length of x (", length(x), ") must equal number of rows in h (", nrow(h), ")")
    if (length(y) != ncol(h))
        stop("length of y (", length(y), ") must equal number of cols in h (", ncol(h), ")")
    # ensure that all three args are double, so the C code won't misinterpret
    dim <- dim(h)
    h <- as.double(h)
    dim(h) <- dim
    x <- as.double(x)
    y <- as.double(y)
    rval <- do_gradient(h, x, y)
    rval$g <- sqrt(rval$gx^2 + rval$gy^2)
    rval
}


#' 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


#' Remove leading and trailing whitespace from strings (deprecated)
#'
#' This function will be removed from an upcoming version of oce, because
#' the base function [trimws()] does the same and more, without having
#' problems with encoding (see `https://github.com/dankelley/oce/issues/1993`).
#'
#' @param s vector of character strings
#'
#' @return a new vector formed by trimming leading and trailing whitespace
#' from the elements of `s`.
#'
#' @family deprecated functions
trimString <- function(s)
{
    .Deprecated("trimString", msg="Use trimws() instead, as of August 2022")
    gsub("^ *", "", gsub(" *$", "", s, perl=TRUE), perl=TRUE)
}

#' Perform 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
}

Try the oce package in your browser

Any scripts or data that you put into this service are public.

oce documentation built on July 9, 2023, 5:18 p.m.