R/cm.R

In dankelley/oce: Analysis of Oceanographic Data

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

#' Class to Store Current Meter Data
#'
#' This class stores current meter data, e.g. from an Interocean/S4 device
#' or an Aanderaa/RCM device.
#'
#' @templateVar class cm
#'
#' @templateVar dataExample The key items stored in this slot are `time`, `u` and `v`.
#'
#' @templateVar metadataExample {}
#'
#' @template slot_summary
#'
#' @template slot_put
#'
#' @template slot_get
#'
#' @family things related to cm data
#' @family classes provided by oce
#'
#' @author Dan Kelley
setClass("cm", contains = "oce")

#' Sample cm Data
#'
#' The result of using [read.cm()] on a current meter file holding measurements made with an
#' Interocean S4 device.  See [read.cm()] for some general cautionary notes on reading such
#' files. Note that the salinities in this sample dataset are known to be incorrect, perhaps
#' owing to a lack of calibration of an old instrument that had not been used in a long time.
#'
#' @name cm
#'
#' @docType data
#'
#' @usage data(cm)
#'
#' @examples
#' library(oce)
#' data(cm)
#' summary(cm)
#' plot(cm)
#'
#' @family datasets provided with oce
#' @family things related to cm data
NULL

#' Extract Something From a cm Object
#'
#' @param x a [cm-class] object.
#'
#' @section Details of the Specialized Method:
#'
#' * If `i` is `"?"`, then the return value is a list
#' containing four items, each of which is a character vector
#' holding the names of things that can be accessed with `[[`.
#' The `data` and `metadata` items hold the names of
#' entries in the object's data and metadata
#' slots, respectively. The `dataDerived`
#' and `metadataDerived` items are each NULL, because
#' no derived values are defined by [cm-class] objects.
#'
#' @template sub_subTemplate
#'
#' @author Dan Kelley
#'
#' @family things related to cm data
setMethod(
    f = "[[",
    signature(x = "cm", i = "ANY", j = "ANY"),
    definition = function(x, i, j, ...) {
        if (i == "?") {
            return(list(
                metadata = sort(names(x@metadata)),
                metadataDerived = NULL,
                data = sort(names(x@data)),
                dataDerived = NULL
            ))
        }
        callNextMethod() # [[
    }
)

#' Replace Parts of a cm Object
#'
#' @param x a [cm-class] object.
#'
#' @template sub_subsetTemplate
#'
#' @family things related to cm data
setMethod(
    f = "[[<-",
    signature(x = "cm", i = "ANY", j = "ANY"),
    definition = function(x, i, j, ..., value) {
        callNextMethod(x = x, i = i, j = j, ... = ..., value = value) # [[<-
    }
)

setMethod(
    f = "initialize",
    signature = "cm",
    definition = function(.Object, time = NULL, u = NULL, v = NULL, units, ...) {
        .Object <- callNextMethod(.Object, ...)
        # sample, direction, conductivity, salinity, temperature, pressure) {
        if (missing(units)) {
            .Object@metadata$units <- list(
                u = list(unit = expression(m / s), scale = ""),
                v = list(unit = expression(m / s), scale = "")
            )
        } else {
            .Object@metadata$units <- units # CAUTION: we are being quite trusting here
        }
        .Object@metadata$north <- "magnetic" # made "geographic" by applyMagneticDeclination()
        .Object@data$time <- time
        .Object@data$u <- u
        .Object@data$v <- v
        .Object@processingLog$time <- presentTime()
        .Object@processingLog$value <- "create 'cm' object"
        return(.Object)
    }
)

#' Summarize a cm Object
#'
#' Summarizes some of the data in a `cm` object, presenting such information
#' as the station name, sampling location, data ranges, etc.
#'
#' @param object A [cm-class] object.
#'
#' @param ... Further arguments passed to or from other methods.
#'
#' @examples
#' library(oce)
#' data(cm)
#' summary(cm)
#'
#' @seealso The documentation for the [cm-class] class explains the structure
#' of `cm` objects, and also outlines the other functions dealing with them.
#'
#' @family things related to cm data
#'
#' @author Dan Kelley
setMethod(
    f = "summary",
    signature = "cm",
    definition = function(object, ...) {
        cat("Cm summary\n----------\n\n", ...)
        showMetadataItem(object, "filename", "File source:   ", quote = TRUE)
        showMetadataItem(object, "type", "Instr. type:   ")
        showMetadataItem(object, "model", "Instr. model:  ")
        showMetadataItem(object, "serialNumber", "Serial Num.:   ")
        showMetadataItem(object, "version", "Version:       ")
        showMetadataItem(object, "north", "North:         ")
        showMetadataItem(object, "declination", "Declination:   ")
        invisible(callNextMethod()) # summary
    }
)


#' Subset a cm Object
#'
#' This function is somewhat analogous to [subset.data.frame()].
#'
#' @param x a [cm-class] object.
#'
#' @param subset a condition to be applied to the `data` portion of `x`.
#' See \dQuote{Details}.
#'
#' @param ... ignored.
#'
#' @return A new `cm` object.
#'
#' @examples
#' library(oce)
#' data(cm)
#' plot(cm)
#' plot(subset(cm, time < mean(range(cm[["time"]]))))
#'
#' @family things related to cm data
#' @family functions that subset oce objects
#'
#' @author Dan Kelley
setMethod(
    f = "subset",
    signature = "cm",
    definition = function(x, subset, ...) {
        subsetString <- paste(deparse(substitute(expr = subset, env = environment())), collapse = " ")
        res <- x
        dots <- list(...)
        debug <- getOption("oceDebug")
        if (length(dots) && ("debug" %in% names(dots))) {
            debug <- dots$debug
        }
        if (missing(subset)) {
            stop("must give 'subset'")
        }
        if (grepl("time", subsetString)) {
            oceDebug(debug, "subsetting a cm by time\n")
            keep <- eval(expr = substitute(expr = subset, env = environment()), envir = x@data, enclos = parent.frame(2))
            names <- names(x@data)
            oceDebug(debug, vectorShow(keep, "keeping bins:"))
            oceDebug(debug, "number of kept bins:", sum(keep), "\n")
            if (sum(keep) < 2) {
                stop("must keep at least 2 profiles")
            }
            res <- x
            # FIXME: are we handling slow timescale data?
            for (name in names(x@data)) {
                if (name == "time" || is.vector(x@data[[name]])) {
                    oceDebug(debug, "subsetting x@data$", name, ", which is a vector\n", sep = "")
                    res@data[[name]] <- x@data[[name]][keep] # FIXME: what about fast/slow
                } else if (is.matrix(x@data[[name]])) {
                    oceDebug(debug, "subsetting x@data$", name, ", which is a matrix\n", sep = "")
                    res@data[[name]] <- x@data[[name]][keep, ]
                } else if (is.array(x@data[[name]])) {
                    oceDebug(debug, "subsetting x@data$", name, ", which is an array\n", sep = "")
                    res@data[[name]] <- x@data[[name]][keep, , , drop = FALSE]
                }
            }
        }
        res@processingLog <- processingLogAppend(res@processingLog, paste("subset.cm(x, subset=", subsetString, ")", sep = ""))
        res
    }
)

#' Coerce Data Into a cm Object
#'
#' @param time A vector of times of observation, or an `oce` object from which time
#' and two velocity components can be inferred, e.g. an [adv-class] object, or
#' an [adp-class] object that has only one distance bin.  If `time` is an `oce` object,
#' then all of the following arguments are ignored.
#'
#' @param u,v optional numerical vectors containing the x and y components of velocity (m/s).
#'
#' @param pressure,conductivity,salinity,temperature optional numerical vectors
#' containing pressure (dbar), electrical conductivity, practical salinity,
#' and in-situ temperature (degree C).
#'
#' @param longitude,latitude optional position specified in degrees East and North.
#'
#' @param filename optional source file name.
#'
#' @template debugTemplate
#'
#' @examples
#' library(oce)
#' # Example 1: creation from scratch
#' t <- Sys.time() + 0:50
#' u <- sin(2 * pi * 0:50 / 5) + rnorm(51)
#' v <- cos(2 * pi * 0:50 / 5) + rnorm(51)
#' p <- 100 + rnorm(51)
#' summary(as.cm(t, u, v, p))
#'
#' # Example 2: creation from an adv object
#' data(adv)
#' summary(as.cm(adv))
#'
#' @family things related to cm data
as.cm <- function(
    time, u = NULL, v = NULL,
    pressure = NULL, conductivity = NULL, temperature = NULL, salinity = NULL,
    longitude = NA, latitude = NA, filename = "", debug = getOption("oceDebug")) {
    oceDebug(debug, "as.ctd() {\n", unindent = 1)
    rpd <- atan2(1, 1) / 45 # radians per degree (avoid 'pi')
    # Catch special cases
    firstArgIsOce <- inherits(time, "oce")
    if (firstArgIsOce) {
        x <- time
        mnames <- names(x@metadata)
        dnames <- names(x@data)
        if (!("time" %in% dnames)) {
            stop("first argument is an oce object, but it does not hold time")
        }
        time <- x@data$time
        if ("pressure" %in% dnames) {
            pressure <- x@data$pressure
        }
        if ("conductivity" %in% dnames) {
            conductivity <- x@data$conductivity
        }
        if ("temperature" %in% dnames) {
            temperature <- x@data$temperature
        }
        if ("salinity" %in% dnames) {
            salinity <- x@data$salinity
        }
        if ("longitude" %in% mnames) {
            longitude <- x@metadata$longitude
        }
        if ("latitude" %in% mnames) {
            latitude <- x@metadata$latitude
        }
        if ("filename" %in% mnames) {
            filename <- x@metadata$filename
        }
        if (inherits(x, "adp")) {
            if (1 == dim(x@data$v)[2]) {
                oceDebug(debug, "x is an adp object with just 1 cell\n")
                u <- as.vector(x@data$v[, 1, 1])
                v <- as.vector(x@data$v[, 1, 2])
            } else {
                stop("can't convert multi-cell adp; try as.cm(x[[\"time\"]],x[[\"v\"]][,1,1], x[[\"v\"]][,1,2])")
            }
        } else if (inherits(x, "adv")) {
            oceDebug(debug, "x is an adv object\n")
            u <- as.vector(x@data$v[, 1])
            v <- as.vector(x@data$v[, 2])
        } else {
            if ("u" %in% dnames && "v" %in% dnames) {
                u <- x@data$u
                v <- x@data$v
            } else if ("speedHorizontal" %in% dnames && "directionTrue" %in% dnames) {
                # NOTE: this can be generalized later to take e.g. 'speed', if some objects have that
                u <- x@data$speedHorizontal * cos(rpd * (90 - x@data$directionTrue))
                v <- x@data$speedHorizontal * sin(rpd * (90 - x@data$directionTrue))
            } else if ("speedHorizontal" %in% dnames && "direction" %in% dnames) {
                u <- x@data$speedHorizontal * cos(rpd * (90 - x@data$direction))
                v <- x@data$speedHorizontal * sin(rpd * (90 - x@data$direction))
            } else {
                stop("first argument must hold either 'u' plus 'v' or 'speed' plus 'directionTrue' or 'direction'")
            }
        }
    }
    direction <- 90 - atan2(v, u) / rpd
    direction <- ifelse(direction < 0, 360 + direction, direction) # put in range 0 to 360
    res <- new("cm", time = time, u = u, v = v)
    res@metadata$filename <- filename
    res@metadata$longitude <- longitude
    res@metadata$latitude <- latitude
    if (!is.null(conductivity)) {
        res@data$conductivity <- conductivity
        res@metadata$units$conductivity <- list(unit = expression(mS / cm), scale = "") # guessing on unit
    }
    if (!is.null(salinity)) {
        res@data$salinity <- salinity
        res@metadata$units$salinity <- list(unit = expression(), scale = "PSS-78")
    }
    if (!is.null(temperature)) {
        res@data$temperature <- temperature
        res@metadata$units$temperature <- list(unit = expression(degree * C), scale = "ITS-90")
    }
    if (!is.null(pressure)) {
        res@data$pressure <- pressure
        res@metadata$units$pressure <- list(unit = expression(dbar), scale = "")
    }
    if (firstArgIsOce) {
        # Copy some metadata that are used sometimes (esp. ODF files)
        if ("type" %in% mnames) {
            res@metadata$type <- x@metadata$type
        }
        if ("model" %in% mnames) {
            res@metadata$model <- x@metadata$model
        }
        if ("serialNumber" %in% mnames) {
            res@metadata$serialNumber <- x@metadata$serialNumber
        }
        if ("ship" %in% mnames) {
            res@metadata$ship <- x@metadata$ship
        }
        if ("scientist" %in% mnames) {
            res@metadata$scientist <- x@metadata$scientist
        }
        if ("cruise" %in% mnames) {
            res@metadata$cruise <- x@metadata$cruise
        }
        if ("station" %in% mnames) {
            res@metadata$station <- x@metadata$station
        }
        if ("countryInstituteCode" %in% mnames) {
            res@metadata$countryInstituteCode <- x@metadata$countryInstituteCode
        }
        if ("cruiseNumber" %in% mnames) {
            res@metadata$cruiseNumber <- x@metadata$cruiseNumber
        }
        if ("waterDepth" %in% mnames) {
            res@metadata$waterDepth <- x@metadata$waterDepth
        }
        # determine original names, where known
        if ("dataNamesOriginal" %in% mnames) {
            if (is.list(x@metadata$dataNamesOriginal)) {
                res@metadata$dataNamesOriginal <- x@metadata$dataNamesOriginal
            } else {
                nameMapping <- as.list(x@metadata$dataNamesOriginal)
                names(nameMapping) <- names(x@data)
                res@metadata$dataNamesOriginal <- lapply(
                    names(res@data),
                    function(name) {
                        if (is.null(nameMapping[[name]])) "-" else nameMapping[[name]]
                    }
                )
            }
        } else {
            res@metadata$dataNamesOriginal <- NULL
        }
        res@processingLog <- x@processingLog
    }
    res@processingLog <- processingLogAppend(
        res@processingLog,
        paste(deparse(match.call()), sep = "", collapse = "")
    )
    oceDebug(debug, "} # as.cm()\n", unindent = 1)
    res
}


#' Read a cm File
#'
#' Read a current-meter data file, producing a [cm-class] object.
#'
#' There function has been tested on only a single file, and the data-scanning
#' algorithm was based on visual inspection of that file.  Whether it will work
#' generally is an open question. It should be noted that the sample file had
#' several odd characteristics, some of which are listed below.
#'
#' * file contained two columns named `"Cond"`, which was guessed
#' to stand for conductivity. Since only the first contained data, the second was
#' ignored, but this may not be the case for all files.
#'
#' * The unit for `"Cond"` was stated in the file to be `"mS"`,
#' which makes no sense, so the unit was assumed to be mS/cm.
#'
#' * The file contained a column named `"T-Temp"`, which is not
#' something the author has seen in his career. It was assumed to stand for
#' in-situ temperature.
#'
#' * The file contained a column named `"Depth"`, which is not something
#' an instrument can measure. Presumably it was calculated from pressure (with
#' what atmospheric offset, though?) and so pressure was inferred from it using
#' [swPressure()].
#'
#' * The file contained several columns that lacked names. These were ignored.
#'
#' * The file contained several columns that seem to be derived from the
#' actual measured data, such as `"Speed"`, `"Dir"`, `"N-S Dist"`,
#' etc. These are ignored.
#'
#' * The file contained several columns that were basically a mystery to the
#' author, e.g. `"Hx"`, `"Hy"`, `"Vref"`, etc. These were ignored.
#'
#' Based on such considerations, [read.cm()] reads only the columns that
#' were reasonably well-understood based on the sample file. Users who need more
#' columns should contact the author. And a user who could produce a document
#' explaining the data format would be especially appreciated!
#'
#' @param file a connection or a character string giving the name of the file to
#' load.
#'
#' @param from index number of the first measurement to be read, or the time of
#' that measurement, as created with [as.POSIXct()] (hint: use
#' `tz="UTC"`).
#'
#' @param to indication of the last measurement to read, in a format matching that
#' of `from`.
#'
#' @param by an indication of the stride length to use while walking through the
#' file. If this is an integer, then `by-1` measurements are skipped between
#' each pair of profiles that is read. This may not make much sense, if the data
#' are not equi-spaced in time.  If `by` is a string representing a time
#' interval, in colon-separated format, then this interval is divided by the
#' sampling interval, to get the stride length. *BUG:* if the data are not
#' equi-spaced, then odd results will occur.
#'
#' @param longitude optional signed number indicating the longitude in degrees
#' East.
#'
#' @param latitude optional signed number indicating the latitude in degrees North.
#'
#' @param type character string indicating type of file (ignored at present).
#'
#' @param tz character string indicating time zone to be assumed in the data.
#'
#' @template encodingTemplate
#'
#' @param debug a flag that turns on debugging.  The value indicates the depth
#' within the call stack to which debugging applies.
#'
#' @param monitor ignored.
#'
#' @param processingLog if provided, the action item to be stored in the log.  This
#' parameter is typically only provided for internal calls; the default that it
#' provides is better for normal calls by a user.
#'
#' @return An [cm-class] object.
#'
#' The `data` slot will contain all the data in the file, with names
#' determined from the tokens in line 3 in that file, passed through
#' [make.names()], except that
#' `Vnorth` is renamed `v` (after conversion from cm/s to m/s),
#' `Veast` is renamed `u` (after conversion from cm/s to m/s),
#' `Cond` is renamed `conductivity`,
#' `T.Temp` is renamed `temperature`
#' and
#' `Sal` is renamed `salinity`, and a new
#' column named `time` (a POSIX time) is constructed
#' from the information in the file header, and another named
#' `pressure` is constructed from the column named `Depth`.
#' At least in the single file studied in the creation of this function,
#' there are some columns that are unnamed in line 3 of the header;
#' these yield data items with names `X`, `X.1`, etc.
#'
#' @section Historical note:
#' Prior to late July, 2016, the direction of current flow was stored in the
#' return value, but it is no longer stored, since it can be derived from the
#' `u` and `v` values.
#'
#' @section Changes:
#' * On 2023-02-09 an item named `north` was added to the `metadata` slot.  This
#' is initialized to `"magnetic"` by [read.cm()], but this is really just a
#' guess, and users ought to consider using [applyMagneticDeclination()] to take
#' magnetic declination into account.
#'
#' @section Sample of Usage:
#' \preformatted{
#' library(oce)
#' cm <- read.oce("cm_interocean_0811786.s4a.tab")
#' summary(cm)
#' plot(cm)
#' }
#'
#' @family things related to cm data
#'
#' @author Dan Kelley
read.cm <- function(
    file, from = 1, to, by = 1, tz = getOption("oceTz"), type = c("s4"),
    longitude = NA, latitude = NA, debug = getOption("oceDebug"),
    encoding = "latin1", monitor = FALSE, processingLog) {
    if (missing(file)) {
        stop("must supply 'file'")
    }
    if (is.character(file)) {
        if (!file.exists(file)) {
            stop("cannot find file \"", file, "\"")
        }
        if (0L == file.info(file)$size) {
            stop("empty file \"", file, "\"")
        }
    }
    oceDebug(debug, "read.cm(file=\"", file,
        "\", from=", format(from),
        ", to=", if (missing(to)) "(missing)" else format(to), ", by=", by, "type=", type, ", ...) {\n",
        sep = "", unindent = 1
    )
    type <- match.arg(type)
    if (type == "s4") {
        read.cm.s4(
            file = file, from = from, to = to, by = by, tz = tz,
            longitude = longitude, latitude = latitude,
            encoding = encoding, monitor = monitor, debug = debug - 1,
            processingLog = processingLog
        )
    } else {
        stop("unknown type of current meter")
    }
}

read.cm.s4 <- function(
    file, from = 1, to, by = 1, tz = getOption("oceTz"), longitude = NA, latitude = NA,
    monitor = FALSE, encoding = "latin1", debug = getOption("oceDebug"), processingLog) {
    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 (debug > 1) {
        debug <- 1L
    }
    oceDebug(debug, "read.cm.s4(file=\"", file,
        "\", from=", format(from),
        ", to=", if (missing(to)) "(missing)" else format(to), ", by=", by, ", ...) {\n",
        sep = "", unindent = 1
    )
    if (is.character(file)) {
        filename <- fullFilename(file)
        file <- file(file, "r", encoding = encoding)
        on.exit(close(file))
    }
    if (!inherits(file, "connection")) {
        stop("argument `file' must be a character string or connection")
    }
    if (!isOpen(file)) {
        filename <- "(connection)"
        open(file, "rb")
        on.exit(close(file))
    }
    # Examine the first line of the file, to get serial number, etc.
    items <- scan(file, "character", nlines = 1, sep = "\t", quiet = TRUE) # slow, but just one line
    oceDebug(debug, "line 1 contains: ", paste(items, collapse = " "), "\n")
    serialNumber <- "unknown"
    version <- "unknown"
    type <- "unknown"
    for (i in 1:(-1 + length(items))) {
        if (length(grep("Serial", items[i]))) {
            serialNumber <- items[i + 1]
        } else if (length(grep("Version", items[i]))) {
            version <- items[i + 1]
        } else if (length(grep("Type", items[i]))) {
            type <- items[i + 1]
        }
    }
    # Skip through the rest of the header, and start paying attention when
    # row number is 1, 2, and then 3.  These first rows give us the time
    # sequence.
    lines <- readLines(file, n = 20)
    for (i in 2:20) {
        items <- strsplit(lines[i], "\t")[[1]]
        oceDebug(debug, "line", i, "contains: ", paste(items, collapse = " "), "\n")
        if (items[1] == "Sample #") {
            # names <- sub("[ ]+$", "", sub("^[ ]+","", items))
            # names <- ifelse(0 == nchar(names), paste("column", seq_along(names), sep=""), names)
        } else if (items[1] == "1") {
            start.day <- items[2]
        } else if (items[1] == "2") {
            start.hms <- items[3]
        } else if (items[1] == "3") {
            t0 <- strptime(paste(start.day, start.hms), format = "%m/%d/%Y %H:%M:%S", tz = tz)
            t1 <- strptime(paste(start.day, items[3]), format = "%m/%d/%Y %H:%M:%S", tz = tz)
            deltat <- as.numeric(t1) - as.numeric(t0)
            break
        }
    }
    # Change names of known items to oce-style names, but leave others
    # as they are, in case the user knows what they mean
    dnames <- strsplit(lines[3], "\t")[[1]]
    dnames <- gsub(" *$", "", gsub("^ *", "", dnames)) # remove leading/trailing blanks
    namesOriginal <- dnames
    dnames <- make.names(dnames, unique = TRUE) # handle duplicated names
    dnames[dnames == "decS"] <- "decS"
    dnames[dnames == "Vnorth"] <- "v"
    dnames[dnames == "Veast"] <- "u"
    dnames[dnames == "Cond"] <- "conductivity"
    dnames[dnames == "T.Temp"] <- "temperature"
    dnames[dnames == "Sal"] <- "salinity"
    # Finally, read the data and chop last 2 lines (which contain footer info
    pushBack(lines, file)
    d <- read.delim(file, skip = 5, sep = "\t", col.names = dnames, stringsAsFactors = FALSE)
    d <- d[seq.int(1L, dim(d)[1] - 2), ]
    res <- new("cm") # will fill in later
    # namesOriginal[namesOriginal==""] <- "-"
    dataNamesOriginal <- as.list(namesOriginal)
    names(dataNamesOriginal) <- dnames
    res@metadata$dataNamesOriginal <- dataNamesOriginal
    d$u <- d$u / 100
    d$v <- d$v / 100
    d$pressure <- swPressure(d$Depth, eos = "gsw") # gsw is faster than unesco with essentially same results
    res@metadata$dataNamesOriginal$pressure <- "-"
    n <- length(d$u)
    time <- seq(t0, by = deltat, length.out = n)
    if (inherits(from, "POSIXt")) {
        if (!inherits(to, "POSIXt")) {
            stop("if \"from\" is POSIXt, then \"to\" must be, also")
        }
        if (!is.numeric(by) || by != 1) {
            stop("sorry, \"by\" must equal 1, in this version of read.cm.s4()")
        }
        # from.to.POSIX <- TRUE
        from.index <- which(time >= from)[1]
        if (is.na(from.index)) {
            from.index <- 1
        }
        to.index <- which(to <= time)[1]
        if (is.na(to.index)) {
            to.index <- n
        }
        oceDebug(
            debug, "Time-based trimming: from=", format(from),
            "to=", format(to), "yield from.index=", from.index, "and to.index=", to.index, "\n"
        )
        keep <- seq(from.index, to.index)
    } else {
        if (!is.numeric(from)) {
            stop("\"from\" must be either POSIXt or numeric")
        }
        to <- n
        if (!is.numeric(to)) {
            stop("\"to\" must be either POSIXt or numeric")
        }
        keep <- seq(from, to)
    }
    keep <- keep[1 <= keep]
    keep <- keep[keep <= n]
    d$time <- time
    d <- d[keep, ]
    d <- as.list(d)
    for (dname in names(d)) {
        if (dname != "Date" && dname != "Time" && dname != "time" && dname != "SRB.Time") {
            d[[dname]] <- as.numeric(d[[dname]])
        }
        if (length(grep("^X[.0-9]*$", dname))) {
            res@metadata$dataNamesOriginal[[dname]] <- "-"
        }
    }
    res@data <- d
    res@metadata$filename <- filename
    res@metadata$serialNumber <- serialNumber
    res@metadata$version <- version
    res@metadata$type <- type
    res@metadata$longitude <- longitude
    res@metadata$latitude <- latitude
    # res@metadata$units$u <- list(unit=expression(m/s), scale="")
    # res@metadata$units$v <- list(unit=expression(m/s), scale="")
    res@metadata$units$conductivity <- list(unit = expression(mS / cm), scale = "")
    res@metadata$units$salinity <- list(unit = expression(), scale = "PSS-78")
    res@metadata$units$temperature <- list(unit = expression(degree * C), scale = "ITS-90")
    res@metadata$units$pressure <- list(unit = expression(dbar), scale = "")
    if (missing(processingLog)) processingLog <- paste(deparse(match.call()), sep = "", collapse = "")
    res@processingLog <- processingLogAppend(res@processingLog, processingLog)
    warning("assuming the compass heading is magnetic; consider using applyMagneticDeclination()")
    oceDebug(debug, "} # read.cm()\n", unindent = 1)
    res
}


#' Plot a cm Object
#'
#' Creates a multi-panel summary plot of data measured by a current meter.
#'
#' The panels are controlled by the `which` argument, as follows.
#'
#' * `which=1` or `which="u"` for a time-series graph of eastward
#' velocity, `u`, as a function of time.
#'
#' * `which=2` or `which="v"` for a time-series graph of
#' northward velocity, `u`, as a function of time.
#'
#' * `which=3` or `"progressive vector"` for progressive-vector
#' plot
#'
#' * `which=4` or `"uv"` for a plot of `v` versus `u`.
#' (Dots are used for small datasets, and smoothScatter for large ones.)
#'
#' * `which=5` or `"uv+ellipse"` as the `"uv"` case, but
#' with an added indication of the tidal ellipse, calculated from the eigen
#' vectors of the covariance matrix.
#'
#' * `which=6` or `"uv+ellipse+arrow"` as the `"uv+ellipse"`
#' case, but with an added arrow indicating the mean current.
#'
#' * `which=7` or `"pressure"` for pressure
#'
#' * `which=8` or `"salinity"` for salinity
#'
#' * `which=9` or `"temperature"` for temperature
#'
#' * `which=10` or `"TS"` for a TS diagram
#'
#' * `which=11` or `"conductivity"` for conductivity
#'
#' * `which=20` or `"direction"` for the direction of flow
#'
#' @param x a [cm-class] object.
#'
#' @param which list of desired plot types.  These are graphed in panels running
#' down from the top of the page.  See \dQuote{Details} for the meanings of various
#' values of `which`.
#'
#' @param type type of plot, as for [plot()].
#'
#' @param xlim,ylim optional limit to the x and y axes, passed to [oce.plot.ts()] for
#' time-series plots.
#'
#' @param xaxs,yaxs optional controls over the limits of the x and y axes,
#' passed to [oce.plot.ts()] for time-series plots.  These values default
#' to `"r"`, meaning to use the regular method of extend the plot past
#' its normal limits.  It is common to use `"i"` to make the graph extend
#' to the panel limits.
#'
#' @param drawTimeRange boolean that applies to panels with time as the horizontal
#' axis, indicating whether to draw the time range in the top-left margin of the
#' plot.
#'
#' @param drawZeroLine boolean that indicates whether to draw zero lines on
#' velocities.
#'
#' @param mgp 3-element numerical vector to use for `par(mgp)`, and also for
#' `par(mar)`, computed from this.  The default is tighter than the R default,
#' in order to use more space for the data and less for the axes.
#'
#' @param mar value to be used with [`par`]`("mar")`.
#'
#' @param small an integer indicating the size of data set to be considered
#' "small", to be plotted with points or lines using the standard
#' [plot()] function.  Data sets with more than `small` points will
#' be plotted with [smoothScatter()] instead.
#'
#' @param main main title for plot, used just on the top panel, if there are
#' several panels.
#'
#' @param tformat optional argument passed to [oce.plot.ts()], for plot
#' types that call that function.  (See [strptime()] for the format
#' used.)
#'
#' @param debug a flag that turns on debugging.  Set to 1 to get a moderate amount
#' of debugging information, or to 2 to get more.
#'
#' @param ... Optional arguments passed to plotting functions.
#'
#' @examples
#' library(oce)
#' data(cm)
#' summary(cm)
#' plot(cm)
#'
#' @family functions that plot oce data
#' @family things related to cm data
#'
#' @aliases plot.cm
#'
#' @author Dan Kelley
setMethod(
    f = "plot",
    signature = signature("cm"),
    definition = function(x,
                          which = c(1:2),
                          type = "l",
                          xlim, ylim, xaxs = "r", yaxs = "r",
                          drawTimeRange = getOption("oceDrawTimeRange"),
                          drawZeroLine = FALSE,
                          mgp = getOption("oceMgp"),
                          mar = c(mgp[1] + 1.5, mgp[1] + 1.5, 1.5, 1.5),
                          small = 2000,
                          main = "",
                          tformat,
                          debug = getOption("oceDebug"),
                          ...) {
        oceDebug(debug, "plot.cm() {\n", unindent = 1)
        oceDebug(debug, "  par(mar)=", paste(par("mar"), collapse = " "), "\n")
        oceDebug(debug, "  par(mai)=", paste(par("mai"), collapse = " "), "\n")
        if (3 != sum(c("time", "u", "v") %in% names(x@data))) {
            warning("In plot,cm-method() :\n  cannot plot a cm object unless its 'data' slot contains 'time', 'u' and 'v'", call. = FALSE)
            return(invisible(NULL))
        }
        opar <- par(no.readonly = TRUE)
        lw <- length(which)
        oceDebug(debug, "length(which) =", lw, "\n")
        if (lw > 1) {
            on.exit(par(opar))
        }
        par(mgp = mgp, mar = mar)
        dots <- list(...)
        oceDebug(debug, "later on in plot,cm-method:\n")
        oceDebug(debug, "  par(mar)=", paste(par("mar"), collapse = " "), "\n")
        oceDebug(debug, "  par(mai)=", paste(par("mai"), collapse = " "), "\n")
        oceDebug(debug, "which:", which, "\n")
        which <- oce.pmatch(
            which,
            list(
                u = 1, v = 2, "progressive vector" = 3,
                "uv" = 4, "uv+ellipse" = 5, "uv+ellipse+arrow" = 6,
                pressure = 7, salinity = 8, temperature = 9, TS = 10, conductivity = 11,
                direction = 20
            )
        )
        oceDebug(debug, "which:", which, "\n")
        tt <- x@data$time
        class(tt) <- "POSIXct" # otherwise image() gives warnings
        if (lw > 1) {
            par(mfrow = c(lw, 1))
            oceDebug(debug, "calling par(mfrow=c(", lw, ", 1)\n")
        }
        len <- length(x@data$u)
        for (w in 1:lw) {
            oceDebug(debug, "which[", w, "]=", which[w], "; drawTimeRange=", drawTimeRange, "\n")
            if (which[w] == 1) {
                oce.plot.ts(x@data$time, x@data$u,
                    type = type,
                    xlim = xlim, ylim = ylim, xaxs = xaxs, yaxs = yaxs,
                    xlab = "", ylab = resizableLabel("u"),
                    main = main, mgp = mgp, mar = c(mgp[1], mgp[1] + 1.5, 1.5, 1.5),
                    tformat = tformat
                )
            } else if (which[w] == 2) {
                oce.plot.ts(x@data$time, x@data$v,
                    type = type,
                    xlim = xlim, ylim = ylim, xaxs = xaxs, yaxs = yaxs,
                    xlab = "", ylab = resizableLabel("v"),
                    main = main, mgp = mgp, mar = c(mgp[1], mgp[1] + 1.5, 1.5, 1.5),
                    tformat = tformat
                )
            } else if (which[w] == 3) {
                # or "progressive vector"
                oceDebug(debug, "progressive vector plot\n")
                dt <- as.numeric(difftime(x@data$time[2], x@data$time[1], units = "sec")) # FIXME: assuming equal dt
                mPerKm <- 1000
                u <- x@data$u
                v <- x@data$v
                u[is.na(u)] <- 0 # zero out missing
                v[is.na(v)] <- 0
                x.dist <- cumsum(u) * dt / mPerKm
                y.dist <- cumsum(v) * dt / mPerKm
                plot(x.dist, y.dist,
                    xlab = resizableLabel("distance km"), ylab = resizableLabel("distance km"),
                    type = "l", asp = 1, ...
                )
            } else if (which[w] %in% 4:6) {
                # "uv" (if 4), "uv+ellipse" (if 5), or "uv+ellipse+arrow" (if 6)
                oceDebug(debug, "\"uv\", \"uv+ellipse\", or \"uv+ellipse+arrow\" plot\n")
                if (len <= small) {
                    plot(x@data$u, x@data$v,
                        type = type,
                        xlab = resizableLabel("u"), ylab = resizableLabel("v"),
                        asp = 1, ...
                    )
                } else {
                    smoothScatter(x@data$u, x@data$v,
                        xlab = resizableLabel("u"), ylab = resizableLabel("v"),
                        asp = 1, ...
                    )
                }
                if (which[w] >= 5) {
                    oceDebug(debug, "\"uv+ellipse\", or \"uv+ellipse+arrow\" plot\n")
                    ok <- !is.na(x@data$u) & !is.na(x@data$v)
                    e <- eigen(cov(data.frame(u = x@data$u[ok], v = x@data$v[ok])))
                    major <- sqrt(e$values[1])
                    minor <- sqrt(e$values[2])
                    theta <- seq(0, 2 * pi, length.out = 360 / 5)
                    xx <- major * cos(theta)
                    yy <- minor * sin(theta)
                    theta0 <- atan2(e$vectors[2, 1], e$vectors[1, 1])
                    rotate <- matrix(c(cos(theta0), -sin(theta0), sin(theta0), cos(theta0)), nrow = 2, byrow = TRUE)
                    xxyy <- rotate %*% rbind(xx, yy)
                    col <- if ("col" %in% names(dots)) col else "darkblue"
                    lines(xxyy[1, ], xxyy[2, ], lwd = 5, col = "yellow")
                    lines(xxyy[1, ], xxyy[2, ], lwd = 2, col = col)
                    if (which[w] >= 6) {
                        oceDebug(debug, "\"uv+ellipse+arrow\" plot\n")
                        umean <- mean(x@data$u, na.rm = TRUE)
                        vmean <- mean(x@data$v, na.rm = TRUE)
                        arrows(0, 0, umean, vmean, lwd = 5, length = 1 / 10, col = "yellow")
                        arrows(0, 0, umean, vmean, lwd = 2, length = 1 / 10, col = col)
                    }
                }
            } else if (which[w] == 7) {
                # an older version had depth stored
                p <- if ("pressure" %in% names(x@data)) {
                    x@data$pressure
                } else {
                    swPressure(x@data$depth, eos = "gsw")
                }
                oce.plot.ts(x@data$time, p,
                    type = type,
                    xlim = xlim, ylim = ylim, xaxs = xaxs, yaxs = yaxs,
                    xlab = "", ylab = resizableLabel("p", "y"),
                    main = main, mgp = mgp, mar = c(mgp[1], mgp[1] + 1.5, 1.5, 1.5),
                    tformat = tformat
                )
            } else if (which[w] == 8) {
                oce.plot.ts(x@data$time, x[["salinity"]],
                    type = type,
                    xlim = xlim, ylim = ylim, xaxs = xaxs, yaxs = yaxs,
                    xlab = "", ylab = resizableLabel("S", "y"),
                    main = main, mgp = mgp, mar = c(mgp[1], mgp[1] + 1.5, 1.5, 1.5),
                    tformat = tformat
                )
            } else if (which[w] == 9) {
                oce.plot.ts(x@data$time, x[["temperature"]],
                    type = type,
                    xlim = xlim, ylim = ylim, xaxs = xaxs, yaxs = yaxs,
                    xlab = "", ylab = resizableLabel("T", "y"),
                    main = main, mgp = mgp, mar = c(mgp[1], mgp[1] + 1.5, 1.5, 1.5),
                    tformat = tformat
                )
            } else if (which[w] == 10) {
                plotTS(as.ctd(x[["salinity"]], x[["temperature"]], x[["pressure"]]), main = main, ...)
            } else if (which[w] == 11) {
                cu <- x[["conductivityUnit"]]
                if (is.list(cu)) {
                    cu <- as.character(cu$unit)
                }
                oce.plot.ts(x@data$time, x@data$conductivity,
                    type = type,
                    xlim = xlim, ylim = ylim, xaxs = xaxs, yaxs = yaxs,
                    xlab = "",
                    ylab = if (0 == length(cu)) {
                        resizableLabel("C", "y")
                    } else if (cu == "mS/cm") {
                        resizableLabel("conductivity mS/cm", "y")
                    } else if (cu == "S/m") {
                        resizableLabel("conductivity S/m", "y")
                    } else {
                        "conductivity (unknown unit"
                    },
                    main = main, mgp = mgp, mar = c(mgp[1], mgp[1] + 1.5, 1.5, 1.5),
                    tformat = tformat
                )
            } else if (which[w] == 20) {
                oce.plot.ts(x@data$time, x@data$direction,
                    type = type,
                    xlim = xlim, ylim = ylim, xaxs = xaxs, yaxs = yaxs,
                    xlab = "", ylab = resizableLabel("direction"),
                    main = main, mgp = mgp, mar = c(mgp[1], mgp[1] + 1.5, 1.5, 1.5),
                    tformat = tformat
                )
            } else {
                stop("unknown value of which (", which[w], ")")
            }
        }
        oceDebug(debug, "} # plot.cm()\n", unindent = 1)
        invisible(NULL)
    }
)

#' Alter a cm Object to Account for Magnetic Declination
#'
#' Current-meter (`cm`) instruments determine directions from onboard compasses,
#' so interpreting velocity components in geographical coordinates requires that
#' magnetic declination be taken into account.  This is what the present
#' function does (see \dQuote{Details}).
#'
#' @template declinationTemplate
#'
#' @param object a [cm-class] object.
#'
#' @param declination numeric value holding magnetic declination in degrees,
#' positive for clockwise from north.
#'
#' @template debugTemplate
#'
#' @return A [cm-class] object, adjusted as outlined in \dQuote{Details}.
#'
#' @seealso Use [magneticField()] to determine the declination,
#' inclination and intensity at a given spot on the world, at a given time.
#'
#' @author Dan Kelley
#'
#' @family things related to magnetism
#' @family things related to cm data
setMethod(
    f = "applyMagneticDeclination",
    signature(object = "cm", declination = "ANY", debug = "ANY"),
    definition = function(object, declination = 0.0, debug = getOption("oceDebug")) {
        oceDebug(debug, "applyMagneticDeclination,cm-method(object, declination=", declination, ") {\n", sep = "", unindent = 1)
        if (length(declination) != 1L) {
            stop("length of 'declination' must equal 1")
        }
        # Note that, unlike the adp and adv methods, we do not check on
        # metadata$coordinate, because it does not exist for cm-class objects.
        if (identical(object@metadata$north, "geographic")) {
            warning("a declination has already been applied, so expect odd results")
        }
        res <- object
        S <- sin(-declination * pi / 180)
        C <- cos(-declination * pi / 180)
        r <- matrix(c(C, S, -S, C), nrow = 2)
        uvr <- r %*% rbind(object@data$u, object@data$v)
        res@data$u <- uvr[1, ]
        res@data$v <- uvr[2, ]
        oceDebug(debug, "originally, u[1:3] =", object@data$u[1:3], "\n")
        oceDebug(debug, "    set u[1:3] to", res@data$u[1:3], "\n")
        oceDebug(debug, "originally, v[1:3] =", object@data$v[1:3], "\n")
        oceDebug(debug, "    set v[1:3] to", res@data$v[1:3], "\n")
        dataNames <- names(object@data)
        for (headingName in c("Hdg", "Hdg.1")) {
            oceDebug(debug, "rotating ", headingName, "\n")
            if (headingName %in% dataNames) {
                oceDebug(debug, "originally, ", headingName, "[1:3] =", object@data[headingName][1:3], "\n", sep = "")
                res@data[[headingName]] <- object@data[[headingName]] + declination
                oceDebug(debug, "    set ", headingName, "[1:3] to ", object@data[headingName][1:3], "\n", sep = "")
            }
        }
        res@metadata$north <- "geographic"
        res@metadata$declination <- declination[1]
        res@processingLog <- processingLogAppend(
            res@processingLog,
            paste0("applyMagneticDeclinationCm(x, declination=", declination[1], ")")
        )
        oceDebug(debug, "} # applyMagneticDeclinationCm\n", unindent = 1)
        res
    }
)