R/argo.R

In dankelley/oce: Analysis of Oceanographic Data

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

#' Class to Store Argo Profiler Data
#'
#' This class stores data from Argo floats.
#'
#' An `argo` object may be read with [read.argo()] or
#' created with [as.argo()].  Argo data can be gridded to constant
#' pressures with [argoGrid()] or subsetted with
#' [subset,argo-method()].  Plots can be made with
#' [plot,argo-method()], while [summary,argo-method()]
#' produces statistical summaries and `show` produces overviews.
#'
#' @templateVar class argo
#'
# nolint start line_length_inter
#' @templateVar dataExample The key items stored in this slot include  equal-length vectors `time`, `longitude`, `latitude` and equal-dimension matrices `pressure`, `salinity`, and `temperature`.
#'
#' @templateVar metadataExample Examples that are of common interest include `id`, a vector of ID codes for the profiles, and `dataMode`, a vector of strings indicating whether the profile is in archived mode (`"A"`), realtime mode (`"R"`), or delayed mode (`"D"`).
# nolint end line_length_inter
#'
#' @template slot_summary
#'
#' @template slot_put
#'
#' @template slot_get
#'
#' @author Dan Kelley and Clark Richards
#'
#' @family classes provided by oce
#' @family things related to argo data
setClass("argo", contains = "oce")

#' Sample argo Data
#'
#' This holds data from ARGO 6900388 in the North Atlantic.
#'
#' Below is the official citation (note that this DOI has web links for
#' downloads):
#'
#' Argo (2017). Argo float data and metadata from Global Data Assembly Centre
#' (Argo GDAC) - Snapshot of Argo GDAC of July, 8st 2017. SEANOE.
#' DOI:\code{10.17882/42182#50865}
#'
#' @name argo
#' @docType data
#'
#' @examples
#' library(oce)
#' data(argo)
#' summary(argo)
#' data(coastlineWorld)
#' plot(argo, which = "trajectory")
#'
#' @source The netcdf file used by [read.argo()] to create this [argo-class]
#' object was downloaded using FTP to
#' \code{ftp.ifremer.fr/ifremer/argo/dac/bodc/6900388/6900388_prof.nc}
#' on 2020 June 24.
#'
#' @family datasets provided with oce
#' @family things related to argo data
NULL


#' Extract Something From an argo Object
#'
#' @param x an [argo-class] object.
#'
#' @templateVar class argo
#'
#' @section Details of the Specialized Method:
#'
#' Note that [argo-class] data may contain both unadjusted data and adjusted
#' data.  By default, this extraction function refers to the former, but a
#' preference for the latter may be set with [preferAdjusted()], the
#' documentation of which explains (fairly complex) details.
#'
#' The results from `argo[[i]]` or `argo[[i,j]]` depend on the
#' nature of `i` and (if provided) `j`. The details are as follows.
#'
#' * 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 hold the names of things
#' that can be inferred from the object's contents, e.g.
#' `"SA"` is named in `dataDerived`, indicating that
#' `argo[["SA"]]` is permitted (to compute Absolute Salinity).
#'
#' * If `i` is `"profile"` and `j` is an integer vector,
#' then an argo object is returned, as specified by `j`. For example,
#' `argo[["profile", 2:5]]` is equivalent to
#' `subset(argo, profile %in% 2:5)`.
#'
#' * If `i` is `"CT"`, then
#' Conservative Temperature is returned, as computed with
#' [`gsw::gsw_CT_from_t`]`(SA,t,p)`, where
#' first `SA` is computed as explained
#' in the next item, `t` is in-situ temperature,
#' and `p` is pressure.
#'
#' * If `i` is `"N2"`, then the square of buoyancy is returned,
#' as computed with [swN2()].
#'
#' * If `i` is `"SA"`, then
#' Absolute Salinity is returned, as computed with
#' [gsw::gsw_SA_from_SP()].
#'
#' * If `i` is `"sigmaTheta"`, then
#' potential density anomaly (referenced to zero
#' pressure) is computed, with [swSigmaTheta()], where the
#' equation of state is taken to be
#' [`getOption`]`("oceEOS",default="gsw")`.
#'
#' * If `i` is `"sigma0"`, `"sigma1"`, `"sigma2"`, `"sigma3"` or `"sigma4"`,
#' then the associated function in the \CRANpkg{gsw} package.
#' For example, `"sigma0"` uses [gsw::gsw_sigma0()], which returns
#' potential density anomaly referenced to 0 dbar,
#' according to the gsw equation of state.
#'
#' * If `i` is `"theta"`, then
#' potential temperature (referenced to zero
#' pressure) is computed, with [swTheta()], where the
#' equation of state is taken to be
#' [`getOption`]`("oceEOS",default="gsw")`.
#'
#' * If `i` is `"depth"`, then
#' a matrix of depths is returned.
#'
#' * If `i` is `"id"` or `"ID"`, then the `id` element within
#' the `metadata` slot is returned.
#'
#' * If `i` is in the `data` slot of `x`,
#' then it is returned, otherwise if it is in the `metadata` slot,
#' then that is returned, otherwise `NULL` is returned.
#'
#' @template sub_subTemplate
#'
#' @examples
#' data(argo)
#' # 1. show that dataset has 223 profiles, each with 56 levels
#' dim(argo[["temperature"]])
#'
#' # 2. show importance of focussing on data flagged 'good'
#' fivenum(argo[["salinity"]], na.rm = TRUE)
#' fivenum(argo[["salinity"]][argo[["salinityFlag"]] == 1], na.rm = TRUE)
#'
#' @author Dan Kelley
#'
#' @family things related to argo data
setMethod(
    f = "[[",
    signature(x = "argo", i = "ANY", j = "ANY"),
    definition = function(x, i, j, ...) {
        res <- NULL
        dots <- list(...)
        debug <- if ("debug" %in% names(dots)) dots$debug else 0
        oceDebug(debug, "[[,argo-method(\"", i, "\") {\n", sep = "", style = "bold", unindent = 1)
        metadataDerived <- c("ID", "cycle", "*Flag", "*Unit")
        dataDerived <- c(
            "profile", "CT", "N2", "SA", "sigmaTheta",
            "theta",
            "z", "depth",
            paste("Absolute", "Salinity"),
            paste("Conservative", "Temperature"),
            paste("sigma", 0:4, sep = ""),
            "spice",
            paste("spiciness", 0:2, sep = "")
        )
        if (i == "?") {
            return(list(
                metadata = sort(names(x@metadata)),
                metadataDerived = sort(metadataDerived),
                data = sort(names(x@data)),
                dataDerived = sort(dataDerived)
            ))
        }
        if (i == "profile") {
            # This assignment to profile is merely to prevent a NOTE from
            # the syntax checker. It is needed because of issues with non-standard
            # evaluation in subset() calls. This is a problem that many
            # package authors have encountered; see e.g.
            # nolint start line_length_linter
            # https://stackoverflow.com/questions/23475309/in-r-is-it-possible-to-suppress-note-no-visible-binding-for-global-variable?noredirect=1&lq=1
            # https://stackoverflow.com/questions/9439256/how-can-i-handle-r-cmd-check-no-visible-binding-for-global-variable-notes-when
            # nolint end line_length_linter

            profile <- NULL # does *not* affect the subset() call that follows
            if (missing(j)) {
                stop("must provide an integer vector to access, e.g. argo[[\"profile\", 1:3]]")
            }
            return(subset(x, profile %in% j))
        }
        namesData <- names(x@data)
        # handle some computed items
        if (i %in% c(
            "CT", paste("conservative", "temperature"), "N2",
            "SA", paste("Absolute", "Salinity"), "sigmaTheta",
            "theta", "spice"
        )) {
            salinity <- x[["salinity", debug = debug - 1]]
            pressure <- x[["pressure", debug = debug - 1]]
            temperature <- x[["temperature", debug = debug - 1]]
            dim <- dim(salinity)
            # Do not need longitude and latitude if eos="unesco", but retain for code clarity
            longitude <- rep(x@data$longitude, each = dim[1])
            latitude <- rep(x@data$latitude, each = dim[1])
            if (i %in% c("CT", "Conservative Temperature")) {
                res <- gsw::gsw_CT_from_t(x[["SA"]], temperature, pressure)
            } else if (i == "N2") {
                # nprofile <- dim[2]
                res <- array(NA_real_, dim = dim)
                for (i in seq_len(dim[2])) {
                    # message("i=",i, ", nprofile=", nprofile)
                    # if (i == 14) browser()
                    if (sum(!is.na(pressure[, i])) > 2) {
                        ctd <- as.ctd(
                            salinity = salinity[, i],
                            temperature = temperature[, i],
                            pressure = pressure[, i],
                            longitude = x@data$longitude[i],
                            latitude = x@data$latitude[i]
                        )
                        res[, i] <- swN2(ctd)
                    } else {
                        res[, i] <- rep(NA, length(salinity[, i]))
                    }
                }
            } else if (i %in% c("SA", "Absolute Salinity")) {
                res <- gsw::gsw_SA_from_SP(salinity, pressure, longitude = longitude, latitude = latitude)
            } else if (i %in% paste("sigma", 0:4, sep = "")) {
                SA <- gsw::gsw_SA_from_SP(salinity, pressure, longitude = longitude, latitude = latitude)
                CT <- gsw::gsw_CT_from_t(SA, temperature, pressure)
                res <- switch(i,
                    "sigma0" = gsw::gsw_sigma0(SA, CT),
                    "sigma1" = gsw::gsw_sigma1(SA, CT),
                    "sigma2" = gsw::gsw_sigma2(SA, CT),
                    "sigma3" = gsw::gsw_sigma3(SA, CT),
                    "sigma4" = gsw::gsw_sigma4(SA, CT)
                )
            } else if (i %in% "spice") {
                if (missing(j)) {
                    res <- swSpice(x)
                } else {
                    if (!j %in% c("gsw", "unesco")) {
                        stop("\"", j, "\" not allowed; use either \"gsw\" or \"unesco\"")
                    }
                    res <- swSpice(x, eos = j)
                }
            } else if (i == "sigmaTheta") {
                res <- swSigmaTheta(salinity,
                    temperature = temperature, pressure = pressure,
                    referencePressure = 0, longitude = longitude, latitude = latitude,
                    eos = getOption("oceEOS", default = "gsw")
                )
            } else if (i == "theta") {
                res <- swTheta(salinity,
                    temperature = temperature, pressure = pressure,
                    referencePressure = 0, longitude = longitude, latitude = latitude,
                    eos = getOption("oceEOS", default = "gsw")
                )
            } else {
                stop("argo[[ coding error: unknown item '", i, "'")
            }
            dim(res) <- dim
        } else if (i == "z") {
            # See note at "depth", below.
            if (is.matrix(x@data$pressure)) {
                n <- dim(x@data$pressure)[1]
                latitude <- matrix(rep(x@data$latitude, each = n), nrow = n, byrow = TRUE)
                res <- -swDepth(x@data$pressure, latitude)
            } else {
                res <- -swDepth(x@data$pressure, x@data$latitude)
            }
        } else if (i == "depth") {
            # This accessor added for issue 1333. Note that the
            # fix for that issue was sometimes calling with
            # vector-form argo object. I don't know how that vector
            # form is arising, but it is likely an index without
            # a drop=FALSE condition ... if I find it, I'll fix it,
            # but the following works fine, so I don't really care too
            # much.
            if (is.matrix(x@data$pressure)) {
                n <- dim(x@data$pressure)[1]
                latitude <- matrix(rep(x@data$latitude, each = n), nrow = n, byrow = TRUE)
                res <- swDepth(x@data$pressure, latitude)
            } else {
                res <- swDepth(x@data$pressure, x@data$latitude)
            }
        } else if (i == "ID" || i == "id") {
            res <- x@metadata$id
        } else if (i == "cycleNumber" || i == "cycle") {
            res <- x@metadata$cycle
        } else if (i == "latitude") {
            res <- x@data$latitude
        } else if (i == "longitude") {
            res <- x@data$longitude
        } else if (i %in% c(namesData, paste0(namesData, "Flag"), paste0(namesData, "Unit"))) {
            # Select adjusted or unadusted variants of things stored in the data slot, or
            # their cousins stored in metadata$units and metadata$flags.
            wantFlag <- grepl("Flag", i)
            wantUnit <- grepl("Unit", i)
            which <- x@metadata$adjustedWhich
            fallback <- x@metadata$adjustedFallback
            iBase <- gsub("Unit$", "", gsub("Flag$", "", i)) # base name for i
            iBaseAdjusted <- paste0(iBase, "Adjusted")
            oceDebug(debug, "i=\"", i,
                "\", iBase=\"", iBase, "\", iBaseAdjusted=\"", iBaseAdjusted,
                "\", wantFlag=", wantFlag, ", wantUnit=", wantUnit, "\n",
                sep = ""
            )
            unadjusted <- if (wantUnit) {
                x@metadata$units[[iBase]]
            } else if (wantFlag) {
                x@metadata$flags[[iBase]]
            } else {
                x@data[[iBase]]
            }
            adjusted <- if (wantUnit) {
                x@metadata$units[[iBaseAdjusted]]
            } else if (wantFlag) {
                x@metadata$flags[[iBaseAdjusted]]
            } else {
                x@data[[iBaseAdjusted]]
            }
            if (!is.null(which) && !is.null(fallback)) {
                # The handling of adjusted/unadjusted preference is carried out in cases;
                # the debugging statements explain the logic flow here in the code, and
                # also expose it to the user (in hopes that users may notice if there are
                # errors with respect to the documented behaviour).
                if (which == "all" || i %in% which) {
                    if (is.null(adjusted)) {
                        oceDebug(debug, "Case 1: returning unadjusted item, since the adjusted item does not exist.\n")
                        res <- unadjusted
                    } else {
                        if (any(is.finite(adjusted))) {
                            oceDebug(debug, "Case 2: returning adjusted data.\n")
                            res <- adjusted
                        } else {
                            if (fallback) {
                                oceDebug(
                                    debug, "Case 3: returning unadjusted data, since all adjusted ",
                                    "values are NA and metadata$adjustedFallback=", fallback, "\n"
                                )
                                res <- unadjusted
                            } else {
                                oceDebug(debug, "Case 4: returning adjusted data, even though all are are NA, ",
                                    "because metadata$adjustedFallback=", fallback, "\n",
                                    sep = ""
                                )
                                res <- adjusted
                            }
                        }
                    }
                } else {
                    oceDebug(debug, "Case 5: returning unadjusted data, because \"", i, "\" is not in metadata$adjustedWhich.\n", sep = "")
                    res <- unadjusted
                }
            } else {
                oceDebug(debug, "Case 6: returning unadjusted data, since metadata slot does not contain adjustedWhich and adjustedFallback\n")
                res <- unadjusted
            }
        } else {
            res <- callNextMethod() # [[ defined in R/AllClass.R
        }
        oceDebug(debug, "} # [[,argo-method\n", sep = "", style = "bold", unindent = 1)
        res
    }
) # [[


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

setMethod(
    f = "initialize",
    signature = "argo",
    definition = function(.Object, time, id, longitude, latitude, salinity, temperature, pressure, filename, dataMode, ...) {
        .Object <- callNextMethod(.Object, ...)
        if (!missing(time)) .Object@data$time <- time
        if (!missing(id)) .Object@metadata$id <- id
        if (!missing(longitude)) .Object@data$longitude <- longitude
        if (!missing(latitude)) .Object@data$latitude <- latitude
        if (!missing(salinity)) .Object@data$salinity <- salinity
        if (!missing(temperature)) .Object@data$temperature <- temperature
        if (!missing(pressure)) .Object@data$pressure <- pressure
        .Object@metadata$filename <- if (missing(filename)) "" else filename
        .Object@metadata$dataMode <- if (missing(dataMode)) "" else dataMode
        .Object@processingLog$time <- presentTime()
        .Object@processingLog$value <- "create 'argo' object"
        .Object <- initializeFlagScheme(.Object, "argo")
        return(.Object)
    }
) # initialize

# a local function -- no need to pollute namespace with it
maybeLC <- function(s, lower) {
    if (lower) tolower(s) else s
}

# a local function -- no need to pollute namespace with it
getData <- function(file, name, quiet = FALSE) {
    res <- NA
    # wrap in capture.output to silence what seems like direct printing to stdout()
    # or stderr() by ncvar_get().
    capture.output({
        res <- try(ncdf4::ncvar_get(file, name), silent = TRUE)
    })
    if (inherits(res, "try-error")) {
        if (!quiet) {
            warning(file$filename, " has no variable named '", name, "'\n", sep = "")
        }
        res <- NULL
    }
    # FIXME: the logic of the next linechanged 2023-01-03 because it looked wrong to me then
    if (is.array(res) && 1 == length(dim(res))) {
        res <- matrix(res)
    } # FIXME: is this right?
    res
}

#' Convert Argo Data Name to Oce Name
#'
#' This function is used internally by [read.argo()] to convert Argo-convention
#' data names to oce-convention names. Users should not call this directly, since
#' its return value may be changed at any moment (e.g. to include units as well
#' as names).
#'
#' Initially, Feb 2016, the inference of names was initially done
#' by an inspection of some data files, based on reference 1. Later, in June
#' 2023, broader inspection of more files and documents yielded about ten
#' additions, and a single correction: `VRSpH` was renamed
#' `phSensorVoltageDifference`, to match related names that had been added.
#'
#' It should be noted that the data files examined contain some names that are not
#' documented in reference 1, and others that are listed only in its changelog,
#' with no actual definitions being given. For example, the files had six distinct
#' variable names that seem to relate to phase in the oxygen sensor, but
#' these are not translated by the present function because these
#' variable names are not defined in reference 1, or not defined uniquely
#' in reference 2.
#'
#' The names are converted with
#' [gsub()], using the `ignore.case` argument of the present
#' function.
#' The procedure
#' is to first handle the items listed in the following table, with string
#' searches anchored to the start of the string. After that,
#' the qualifiers
#' `_ADJUSTED`, `_ERROR` and `_QC`,
#' are translated to `Adjusted`, `Error`, and `QC`, respectively.
#'
#' An incomplete list of name translations is as follows, where `~`
#' represents digit sequences in some instances and letters
#' in others.  Note that until June 2023, `pHSensorVoltageDifference`
#' was called `VRSpH`.
#'
#' \tabular{ll}{
#' **Argo name** \tab **oce name**\cr
#' `BBP` \tab `bbp`\cr
#' `BETA_BACKSCATTERING` \tab `betaBackscattering`\cr
#' `BPHASE_OXY` \tab `bphaseOxygen`\cr
#' `C~PHASE_DOXY` \tab `C~phaseOxygen`\cr
#' `CDOM` \tab `CDOM`\cr
#' `CNDC` \tab `conductivity`\cr
#' `CHLA` \tab `chlorophyllA`\cr
#' `CP` \tab `beamAttenuation`\cr
#' `CYCLE_NUMBER` \tab `cycleNumber` (both this and `cycle` are handled by the [[ operator)\cr
#' `DATA_CENTRE` \tab `dataCentre`\cr
#' `DATA_MODE` \tab `dataMode`\cr
#' `DATA_STATE_INDICATOR` \tab `dataStateIndicator`\cr
#' `DC_REFERENCE` \tab `DCReference`\cr
#' `DIRECTION` \tab `direction`\cr
#' `DOWN_IRRADIANCE` \tab `downwellingIrradiance`\cr
#' `DOWNWELLING_PAR` \tab `downwellingPAR`\cr
#' `FIRMWARE_VERSION` \tab `firmwareVersion`\cr
#' `FIT_ERROR_NITRATE` \tab `fitErrorNitrate`\cr
#' `FLUORESCENCE_CDOM` \tab `fluorescenceCDOM`\cr
#' `FLUORESCENCE_CHLA` \tab `fluorescenceChlorophyllA`\cr
#' `IB_PH` \tab `pHBaseCurrent`\cr
#' `IK_PH` \tab `pHCounterCurrent`\cr
#' `INST_REFERENCE` \tab `instReference`\cr
#' `JULD` \tab `juld` (and used to compute `time`)\cr
#' `JULD_QC_LOCATION` \tab `juldQCLocation`\cr
#' `LATITUDE` \tab `latitude`\cr
#' `LONGITUDE` \tab `longitude`\cr
#' `MOLAR_DOXY` \tab `oxygenUncompensated`\cr
#' `MTIME` \tab `mtime`\cr
#' `NB_SAMPLE_CTD` \tab `nbSampleCtd`\cr
#' `PH_IN_SITU_FREE` \tab `pHFree`\cr
#' `PH_IN_SITU_TOTAL` \tab `pH`\cr
#' `PI_NAME` \tab `PIName`\cr
#' `PLATFORM_NUMBER` \tab `id`\cr
#' `POSITION_ACCURACY` \tab `positionAccuracy`\cr
#' `POSITIONING_SYSTEM` \tab `positioningSystem`\cr
#' `PROFILE` \tab `profile`\cr
#' `PROJECT_NAME` \tab `projectName`\cr
#' `RAW_DOWNWELLING_IRRADIANCE` \tab `rawDownwellingIrradiance`\cr
#' `RAW_DOWNWELLING_PAR` \tab `rawDownwellingPAR`\cr
#' `RAW_UPWELLING_RADIANCE` \tab `rawUpwellingRadiance`\cr
#' `STATION_PARAMETERS` \tab `stationParameters`\cr
#' `TEMP` \tab `temperature`\cr
#' `TEMP_CPU_CHLA` \tab `temperatureCPUChlorophyllA`\cr
#' `TEMP_DOXY` \tab `temperatureOxygen`\cr
#' `TEMP_NITRATE` \tab `temperatureNitrate`\cr
#' `TEMP_PH` \tab `temperaturePH`\cr
#' `TEMP_SPECTROPHOTOMETER_NITRATE` \tab `temperatureSpectrophotometerNitrate`\cr
#' `TILT` \tab `tilt`\cr
#' `TPHASE_DOXY` \tab `tphaseOxygen`\cr
#' `TURBIDITY` \tab `turbidity`\cr
#' `UP_RADIANCE` \tab `upwellingRadiance`\cr
#' `UV_INTENSITY` \tab `UVIntensity`\cr
#' `UV_INTENSITY_DARK_NITRATE` \tab `UVIntensityDarkNitrate`\cr
#' `UV_INTENSITY_NITRATE` \tab `UVIntensityNitrate`\cr
#' `VRS_PH` \tab `pHSensorVoltageDifference`\cr
#' `WMO_INST_TYPE` \tab `WMOInstType`\cr
#' }
#'
#' @param names vector of character strings containing names in the Argo convention.
#'
#' @param ignore.case a logical value passed to [gsub()], indicating whether to
#' ignore the case of input strings. The default is set to `TRUE` because some data
#' files use lower-case names, despite the fact that the Argo documentation specifies
#' upper-case.
#'
#' @return A character vector of the same length as `names`, but with
#' replacements having been made for all known quantities.
#'
#' @references
#' 1. Argo User's Manual Version 3.3, Nov 89th, 2019, available at
#' `https://archimer.ifremer.fr/doc/00187/29825/` online.
#'
#' 2. Argo list of parameters in an excel spreadsheet, available at
#' `http://www.argodatamgt.org/content/download/27444/187206/file/argo-parameters-list-core-and-b.xlsx`
#'
#' @author Dan Kelley, with help from Anna Victor
#'
#' @family things related to argo data
argoNames2oceNames <- function(names, ignore.case = TRUE) {
    # do NOT change the order below, because we are working with partial strings.
    names <- gsub("^BBP([0-9_]*)", "BBP\\1", names, ignore.case = ignore.case)
    names <- gsub("^BETA_BACKSCATTERING([0-9_]*)", "betaBackscattering\\1", names, ignore.case = ignore.case)
    names <- gsub("^BPHASE_DOXY", "bphaseOxygen", names, ignore.case = ignore.case)
    names <- gsub("^C([1-9])PHASE_DOXY", "c\\1phaseOxygen", names, ignore.case = ignore.case)
    names <- gsub("^CHLA", "chlorophyllA", names, ignore.case = ignore.case)
    names <- gsub("^CDOM", "CDOM", names, ignore.case = ignore.case)
    names <- gsub("^CNDC([0-9_]*)", "conductivity\\1", names, ignore.case = ignore.case)
    names <- gsub("^CP([0-9_]*)", "beamAttenuation\\1", names, ignore.case = ignore.case)
    names <- gsub("^CYCLE_NUMBER", "cycleNumber", names, ignore.case = ignore.case)
    names <- gsub("^DOWN_IRRADIANCE", "downwellingIrradiance", names, ignore.case = ignore.case)
    names <- gsub("^DOWNWELLING_PAR", "downwellingPAR", names, ignore.case = ignore.case)
    names <- gsub("^DOXY", "oxygen", names, ignore.case = ignore.case)
    names <- gsub("^FIT_ERROR_NITRATE", "fitErrorNitrate", names, ignore.case = ignore.case)
    names <- gsub("^FLUORESCENCE_CDOM", "fluorescenceCDOM", names, ignore.case = ignore.case)
    names <- gsub("^FLUORESCENCE_CHLA", "fluorescenceChlorophyllA", names, ignore.case = ignore.case)
    names <- gsub("^IB_PH", "pHSensorBaseCurrent", names, ignore.case = ignore.case)
    names <- gsub("^IK_PH", "pHSensorCounterCurrent", names, ignore.case = ignore.case)
    names <- gsub("^LATITUDE", "latitude", names, ignore.case = ignore.case)
    names <- gsub("^LONGITUDE", "longitude", names, ignore.case = ignore.case)
    names <- gsub("^MOLAR_DOXY", "oxygenUncompensated", names, ignore.case = ignore.case)
    names <- gsub("^MTIME", "mtime", names, ignore.case = ignore.case)
    names <- gsub("^NB_SAMPLE_CTD", "nbSampleCTD", names, ignore.case = ignore.case)
    names <- gsub("^NB_SAMPLE_(.*)", "nbSample\\1", names, ignore.case = ignore.case)
    names <- gsub("^PH_IN_SITU_FREE", "pHFree", names, ignore.case = ignore.case)
    names <- gsub("^PH_IN_SITU_TOTAL", "pH", names, ignore.case = ignore.case)
    names <- gsub("^NITRATE", "nitrate", names, ignore.case = ignore.case)
    names <- gsub("^PHASE_DELAY_DOXY", "phaseDelayOxygen", names, ignore.case = ignore.case)
    names <- gsub("^POSITION_ACCURACY", "positionAccuracy", names, ignore.case = ignore.case)
    names <- gsub("^PRES", "pressure", names, ignore.case = ignore.case)
    names <- gsub("^PSAL", "salinity", names, ignore.case = ignore.case)
    names <- gsub("^RAW_DOWNWELLING_IRRADIANCE", "rawDownwellingIrradiance", names, ignore.case = ignore.case)
    names <- gsub("^RAW_DOWNWELLING_PAR", "rawDownwellingPAR", names, ignore.case = ignore.case)
    names <- gsub("^RAW_UPWELLING_RADIANCE", "rawUpwellingRadiance", names, ignore.case = ignore.case)
    names <- gsub("^TEMP_CNDC", "temperatureConductivityCell", names, ignore.case = ignore.case)
    names <- gsub("^TEMP_CPU_CHLA", "temperatureCPUChlA", names, ignore.case = ignore.case)
    names <- gsub("^TEMP_DOXY", "temperatureOxygen", names, ignore.case = ignore.case)
    names <- gsub("^TEMP_NITRATE", "temperatureNitrate", names, ignore.case = ignore.case)
    names <- gsub("^TEMP_PH", "temperaturePH", names, ignore.case = ignore.case)
    names <- gsub("^TEMP_SPECTROPHOTOMETER_NITRATE", "temperatureSpectrophotometerNitrate", names, ignore.case = ignore.case)
    names <- gsub("^TEMP_VOLTAGE_DOXY", "temperatureVoltageOxygen", names, ignore.case = ignore.case)
    names <- gsub("^TEMP_", "temperature_", names, ignore.case = ignore.case)
    names <- gsub("^TEMP([0-9_]*)$", "temperature\\1", names, ignore.case = ignore.case)
    names <- gsub("^TILT([0-9_]*)$", "tilt\\1", names, ignore.case = ignore.case)
    names <- gsub("^TPHASE_DOXY", "tphaseOxygen", names, ignore.case = ignore.case)
    names <- gsub("^TRANSMITTANCE_PARTICLE_BEAM_ATTENUATION([0-9_]*)$", "transmittanceParticleBeamAttenuation\\1", names, ignore.case = ignore.case)
    names <- gsub("^TURBIDITY([0-9_]*)$", "turbidity\\1", names, ignore.case = ignore.case)
    names <- gsub("^UP_RADIANCE", "upwellingRadiance", names, ignore.case = ignore.case)
    names <- gsub("^UV_INTENSITY_DARK_NITRATE", "UVIntensityDarkNitrate", names, ignore.case = ignore.case)
    names <- gsub("^UV_INTENSITY_NITRATE", "UVIntensityNitrate", names, ignore.case = ignore.case)
    names <- gsub("^UV_INTENSITY", "UVIntensity", names, ignore.case = ignore.case)
    names <- gsub("^VRS_PH", "pHSensorVoltageDifference", names, ignore.case = ignore.case) # 2023-06-03: was
    names <- gsub("^VK_PH", "pHSensorCounterVoltage", names, ignore.case = ignore.case)
    names <- gsub("_ADJUSTED", "Adjusted", names, ignore.case = ignore.case)
    names <- gsub("_ERROR", "Error", names, ignore.case = ignore.case)
    names <- gsub("_MED", "Med", names, ignore.case = ignore.case)
    names <- gsub("_QC", "QC", names, ignore.case = ignore.case)
    names <- gsub("_STD", "Std", names, ignore.case = ignore.case)
    names
}

#' Subset an argo Object
#'
#' Subset an argo object, either by selecting just the "adjusted" data
#' or by subsetting by pressure or other variables.
#'
#' @details
#' If `subset` is the string `"adjusted"`, then `subset`
#' replaces the station variables with their adjusted counterparts. In
#' the argo notation, e.g. `PSAL` is replaced with `PSAL_ADJUSTED`;
#' in the present notation, this means that `salinity` in the `data`
#' slot is replaced with `salinityAdjusted`, and the latter is deleted.
#' Similar replacements are also done with the flags stored in the `metadata`
#' slot.
#'
#' If `subset` is an expression, then the action is somewhat similar
#' to other `subset` functions, but with the restriction that
#' only one independent variable may be
#' used in in any call to the function, so that
#' repeated calls will be necessary to subset based on more than one
#' independent variable.  Subsetting may be done by anything
#' stored in the data, e.g. `time`,
#' `latitude`, `longitude`, `profile`, `dataMode`,
#' or `pressure` or by `profile` (a made-up variable),
#' `id` (from the `metadata` slot) or `ID` (a synonym for `id`).
#' Note that subsetting by `pressure`
#' preserves matrix shape, by setting discarded values to `NA`, as opposed
#' to dropping data (as is the case with `time`, for example).
#'
#' @param x an [argo-class] object.
#'
#' @param subset An expression indicating how to subset `x`.
#'
#' @param ... optional arguments, of which only the first is examined. The
#' only possibility is `within`, a polygon enclosing data to be
#' retained. This must be either a list or data frame, containing items
#' named either `x` and `y` or `longitude` and
#' `latitude`; see Example 4.  If `within` is given,
#' then `subset` is ignored.
#'
#' @return An [argo-class] object.
#'
#' @examples
#' library(oce)
#' data(argo)
#'
#' # Example 1: subset by time, longitude, and pressure
#' par(mfrow = c(2, 2))
#' plot(argo)
#' plot(subset(argo, time > mean(time)))
#' plot(subset(argo, longitude > mean(longitude)))
#' plot(subset(argoGrid(argo), pressure > 500 & pressure < 1000), which = 5)
#'
#' @section Sample of Usage:
#' \preformatted{
#' # Example 2: restrict attention to delayed-mode profiles.
#' par(mfrow=c(1, 1))
#' plot(subset(argo, dataMode == "D"))
#'
#' # Example 3: contrast adjusted and unadjusted data
#' par(mfrow=c(1, 2))
#' plotTS(argo)
#' plotTS(subset(argo, "adjusted"))
#'
#' # Example 2. Subset by a polygon determined with locator()
#' par(mfrow=c(1, 2))
#' plot(argo, which="map")
#' # Can get a boundary with e.g. locator(4)
#' boundary <- list(x=c(-65, -40, -40, -65), y=c(65, 65, 45, 45))
#' argoSubset <- subset(argo, within=boundary)
#' plot(argoSubset, which="map")
#' }
#
#' @author Dan Kelley
#'
#' @family things related to argo data
#' @family functions that subset oce objects
#' @aliases subset.argo
setMethod(
    f = "subset",
    signature = "argo",
    definition = function(x, subset, ...) {
        subsetString <- paste(deparse(substitute(expr = subset, env = environment())), collapse = " ")
        res <- x
        dots <- list(...)
        dotsNames <- names(dots)
        withinGiven <- length(dots) && ("within" %in% dotsNames)
        debug <- if (length(dots) && ("debug" %in% names(dots))) dots$debug else getOption("oceDebug")
        oceDebug(debug, "subset,argo-method() {\n", sep = "", unindent = 1, style = "bold")
        if (withinGiven) {
            oceDebug(debug, "subsetting with 'within' method\n")
            polygon <- dots$within
            polygonNames <- names(polygon)
            latp <- if ("y" %in% polygonNames) {
                polygon$y
            } else if ("latitude" %in% polygonNames) {
                polygon$latitude
            } else {
                stop("'within' must contain either 'y' or 'latitude'")
            }
            lonp <- if ("x" %in% polygonNames) {
                polygon$x
            } else if ("longitude" %in% polygonNames) {
                polygon$longitude
            } else {
                stop("'within' must contain either 'x' or 'longitude'")
            }
            lon <- x[["longitude", "byStation"]]
            lat <- x[["latitude", "byStation"]]
            poly <- sf::st_polygon(list(outer = cbind(c(lonp, lonp[1]), c(latp, latp[1]))))
            points <- sf::st_multipoint(cbind(lon, lat))
            inside <- sf::st_intersection(points, poly)
            keep <- matrix(points %in% inside, ncol = 2)[, 1]
            #<OLD>if (!all.equal(keepNew, keep)) {
            #<OLD>    warning("subset,argo-method disagreement between old 'sp' method and new 'sf' method\n")
            #<OLD>} else {
            #<OLD>    oceDebug(debug, "subset,argo-method: old 'sp' method and new 'sf' method gave identical results\n")
            #<OLD>}
            # }NEW
            # Metadata
            for (name in names(x@metadata)) {
                oceDebug(debug, "subsetting metadata item named '", name, "'.\n", sep = "")
                # Pass some things through directly.
                # 20200831 if (name %in% c("units", "flags", "filename", "flagScheme", "dataNamesOriginal")) {
                if (name %in% c("units", "filename", "flagScheme", "dataNamesOriginal")) {
                    next
                }
                item <- x@metadata[[name]]
                # Handle things that are encoded as characters in a string, namely 'direction', 'juldQC', 'positionQC',
                # and some other 'QC` things that are found by grepping.
                if (name == "direction" || grepl("QC$", name)) {
                    oceDebug(debug, "  \"", name, "\" is a special string (\"direction\" or \"*QC\"), being subsetted by character number\n", sep = "")
                    res@metadata[[name]] <- paste(strsplit(item, "")[[1]][keep], collapse = "")
                } else if (is.list(item)) {
                    oceDebug(debug, "  \"", name, "\" is a list, with each element being subsetted\n", sep = "")
                    for (l in seq_along(item)) {
                        res@metadata[[name]][[l]] <- item[[l]][, keep, drop = FALSE]
                    }
                } else if (is.vector(name)) {
                    res@metadata[[name]] <- item[keep]
                } else if (is.matrix(name)) {
                    res@metadata[[name]] <- item[, keep, drop = FALSE]
                } else if (is.array(name)) {
                    oceDebug(debug, "name=", name, " has dim ", paste(dim(res@metadata[[name]]), collapse = " "), "\n")
                    if (length(dim(res@metadata[[name]])) <= 3) {
                        res@metadata[[name]] <- item[, , keep, drop = FALSE]
                    } else {
                        warning("not subsetting \"", name, "\" in metadata, because it is an array of rank > 3")
                    }
                } else {
                    stop("cannot subset metadata item named '", name, "' because it is not a length-one string, a vector, or a matrix")
                }
            }
            # Data
            for (name in names(x@data)) {
                oceDebug(debug, "subsetting data item named '", name, "'\n", sep = "")
                item <- x@data[[name]]
                if ("time" == name) {
                    res@data$time <- item[keep]
                } else if (is.vector(item)) {
                    res@data[[name]] <- item[keep]
                } else if (is.matrix(item)) {
                    res@data[[name]] <- item[, keep]
                } else {
                    stop("argo object has data item '", name, "' that is neither a vector nor a matrix, so we cannot subset it")
                }
            }
            res@processingLog <- processingLogAppend(
                res@processingLog,
                paste("subset(x, within) kept ", sum(keep), " of ",
                    length(keep), " stations",
                    sep = ""
                )
            )
        } else {
            if (is.character(substitute(expr = subset, env = environment()))) {
                if (subset != "adjusted") {
                    stop("if subset is a string, it must be \"adjusted\"")
                }
                dataNames <- names(x@data)
                # Seek 'Adjusted' names
                adjustedIndices <- grep(".*Adjusted$", dataNames)
                for (i in adjustedIndices) {
                    adjusted <- dataNames[i]
                    base <- gsub("Adjusted$", "", adjusted)
                    adjustedError <- paste(adjusted, "Error", sep = "")
                    res@data[[base]] <- res@data[[adjusted]]
                    res@data[[adjusted]] <- NULL
                    res@data[[adjustedError]] <- NULL
                }
                flagNames <- names(x@metadata$flags)
                adjustedIndices <- grep("Adjusted", flagNames)
                for (i in adjustedIndices) {
                    adjusted <- flagNames[i]
                    base <- gsub("Adjusted", "", adjusted)
                    adjustedError <- paste(adjusted, "Error", sep = "")
                    res@metadata$flags[[base]] <- res@metadata$flags[[adjusted]]
                    res@metadata$flags[[adjusted]] <- NULL
                    res@metadata$flags[[adjustedError]] <- NULL
                }
                res@processingLog <- processingLogAppend(
                    res@processingLog,
                    paste("subset.argo(x, subset=\"", as.character(subset), "\")", sep = "")
                )
            } else {
                subsetString <- paste(deparse(substitute(expr = subset, env = environment())), collapse = " ")
                res <- x
                if (length(grep("time", subsetString)) ||
                    length(grep("longitude", subsetString)) || length(grep("latitude", subsetString))) {
                    keep <- eval(expr = substitute(expr = subset, env = environment()), envir = x@data, enclos = parent.frame(2))
                } else if (length(grep("id", subsetString, ignore.case = TRUE))) {
                    # add id into the data, then do as usual
                    tmp <- x@data
                    tmp$id <- x@metadata$id
                    keep <- eval(expr = substitute(expr = subset, env = environment()), envir = tmp, enclos = parent.frame(2))
                    rm(tmp)
                } else if (length(grep("profile", subsetString))) {
                    # add profile into the data, then do as usual
                    tmp <- x@data
                    tmp$profile <- seq_along(x@data$time)
                    keep <- eval(expr = substitute(expr = subset, env = environment()), envir = tmp, enclos = parent.frame(2))
                    rm(tmp)
                } else if (length(grep("pressure", subsetString))) {
                    # issue1628 # check that it is a "gridded" argo
                    # issue1628 gridded <- ifelse(all(apply(x@data$pressure, 1, diff) == 0, na.rm=TRUE), TRUE, FALSE)
                    # issue1628 if (gridded) {
                    # issue1628     x@data$pressure <- x@data$pressure[, 1] # FIXME: have to convert pressure to vector
                    # issue1628     keep <- eval(substitute(subset), x@data, parent.frame(2))
                    # issue1628     x@data$pressure <- res@data$pressure # FIXME: convert back to original for subsetting below
                    # issue1628 } else {
                    # issue1628     stop("cannot subset ungridded argo by pressure -- use argoGrid() first", call.=FALSE)
                    # issue1628 }
                    keep <- eval(expr = substitute(expr = subset, env = environment()), envir = x@data, enclos = parent.frame(2))
                } else if (length(grep("dataMode", subsetString))) {
                    keep <- eval(expr = substitute(expr = subset, env = environment()), envir = x@metadata, enclos = parent.frame(2))
                } else {
                    stop("can only subset by time, longitude, latitude, pressure, dataMode, and not by combinations", call. = FALSE)
                }
                if (length(grep("pressure", subsetString))) {
                    # Now do the subset. Note that we preserve matrix dimensions, by setting
                    # discarded values to NA.
                    fieldname <- names(x@data)
                    for (field in fieldname) {
                        if (field != "time" & field != "longitude" & field != "latitude") { # DEBUG: see issue 1327
                            ifield <- which(field == fieldname)
                            # debug message("ifield=", ifield, ", field=", field,
                            # debug        "\n\tlength(keep)=", length(keep),
                            # debug        "\n\tsum(keep)=", sum(keep))
                            if (is.matrix(res@data[[ifield]])) {
                                res@data[[ifield]][!keep] <- NA
                            } else {
                                res@data[[ifield]][!keep] <- NA
                            }
                        }
                    }
                    fieldname <- names(x@metadata$flags)
                    for (field in fieldname) {
                        ifield <- which(field == fieldname)
                        res@metadata$flags[[ifield]][!keep] <- NA
                    }
                    # res@data$salinity <- x@data$salinity[keep, ]
                    # res@data$temperature <- x@data$temperature[keep, ]
                    # res@data$pressure <- x@data$pressure[keep, ]
                    res@processingLog <- processingLogAppend(res@processingLog, paste("subset.argo(x, subset=", subsetString, ")", sep = ""))
                } else {
                    res@data$time <- x@data$time[keep]
                    res@data$longitude <- x@data$longitude[keep]
                    res@data$latitude <- x@data$latitude[keep]
                    res@data$profile <- x@data$profile[keep]
                    res@metadata$dataMode <- x@metadata$dataMode[keep]
                    fieldname <- names(x@data)
                    for (field in fieldname) {
                        if (field != "time" && field != "longitude" && field != "latitude" && field != "profile") {
                            ifield <- which(field == fieldname)
                            res@data[[ifield]] <- if (is.matrix(x@data[[ifield]])) {
                                x@data[[ifield]][, keep, drop = FALSE]
                            } else {
                                x@data[[ifield]][keep]
                            }
                        }
                    }
                    fieldname <- names(x@metadata$flags)
                    for (field in fieldname) {
                        ifield <- which(field == fieldname)
                        res@metadata$flags[[ifield]] <- res@metadata$flags[[ifield]][, keep]
                    }
                    # if (sum(keep) < 1) warning("In subset.argo() :\n  removed all profiles", call.=FALSE)
                    # res@data$salinity <- x@data$salinity[, keep]
                    # res@data$temperature <- x@data$temperature[, keep]
                    # res@data$pressure <- x@data$pressure[, keep]
                }
                res@processingLog <- processingLogAppend(res@processingLog, paste("subset.argo(x, subset=", subsetString, ")", sep = ""))
            }
        }
        oceDebug(debug, "} # subset,argo-method\n", sep = "", unindent = 1, style = "bold")
        res
    }
)


#' Summarize an argo Object
#'
#' @description Summarizes some of the data in an `argo` object.
#'
#' @details Pertinent summary information is presented.
#' @param object  an object of class `"argo"`, usually, a result of a call to [read.argo()].
#' @param ... Further arguments passed to or from other methods.
#'
#' @return A matrix containing statistics of the elements of the `data` slot.
#' @examples
#' library(oce)
#' data(argo)
#' summary(argo)
#'
#' @family things related to argo data
#'
#' @aliases summary.argo
#'
#' @author Dan Kelley
setMethod(
    f = "summary",
    signature = "argo",
    definition = function(object, ...) {
        cat("Argo Summary\n------------\n\n")
        showMetadataItem(object, "filename", "Source:              ", quote = TRUE)
        nid <- length(unique(object@metadata$id))
        if (1 == nid) {
            cat("* ID:                  \"", object@metadata$id[1], "\"\n", sep = "")
        } else {
            cat("* ID list:             \"", object@metadata$id[1], "\", \"", object@metadata$id[2], "\", ...\n", sep = "")
        }
        if ("featureType" %in% names(object@metadata)) {
            cat("* Feature type:        \"", object@metadata$featureType, "\"\n", sep = "")
        }
        nD <- sum(object@metadata$dataMode == "D")
        nA <- sum(object@metadata$dataMode == "A")
        nR <- sum(object@metadata$dataMode == "R")
        cat("* Profiles:            ", nD, " delayed; ", nA, " adjusted; ", nR, " realtime", "\n", sep = "")
        if ("time" %in% names(object@metadata)) {
            cat("* Time:                ", format(object@metadata$time), "\n", sep = "")
        }
        invisible(callNextMethod()) # summary
    }
)

ncdfFixMatrix <- function(x) {
    if (length(dim(x)) == 1L) as.vector(x) else x
}

#' Grid Argo Float Data
#'
#' Grid an Argo float, by interpolating to fixed pressure levels.
#' The gridding is done with [approx()].  If there is
#' sufficient user demand, other methods may be added, by analogy to
#' [sectionGrid()].
#'
#' @template flagDeletionTemplate
#'
#' @param argo A `argo` object to be gridded.
#'
#' @param p Optional indication of the pressure levels to which interpolation
#' should be done.  If this is not supplied, the pressure levels will be
#' calculated based on the existing values, using medians. If `p="levitus"`,
#' then pressures will be set to be those of the Levitus atlas, given by
#' [standardDepths()], trimmed to the maximum pressure in `argo`.
#' If `p` is a single numerical value, it is taken as the number of
#' subdivisions to use in a call to [seq()] that has range from 0 to the
#' maximum pressure in `argo`.  Finally, if a vector numerical values is
#' provided, then it is used as is.
#'
#' @param debug A flag that turns on debugging.  Higher values provide deeper
#' debugging.
#'
#' @param ... Optional arguments to [approx()], which is used to do the
#' gridding.
#'
#' @return x an [argo-class] object.
#'
#' @examples
#' library(oce)
#' data(argo)
#' g <- argoGrid(argo, p = seq(0, 100, 1))
#' par(mfrow = c(2, 1))
#' t <- g[["time"]]
#' z <- -g[["pressure"]][, 1]
#' # Set zlim because of spurious temperatures.
#' imagep(t, z, t(g[["temperature"]]), ylim = c(-100, 0), zlim = c(0, 20))
#' imagep(t, z, t(g[["salinity"]]), ylim = c(-100, 0))
#'
#' @family things related to argo data
#' @author Dan Kelley and Clark Richards
argoGrid <- function(argo, p, debug = getOption("oceDebug"), ...) {
    oceDebug(debug, "argoGrid() {\n", sep = "", unindent = 1)
    warningMessages <- NULL
    dim <- dim(argo@data$pressure)
    # ndepth <- dim[1]
    nprofile <- dim[2]
    # FIXME: modify sal, temp, and pre.  In the end, pre constant along first index
    res <- argo
    res[["flags"]] <- NULL
    warningMessages <- c(
        warningMessages,
        "Data flags are omitted from the gridded argo object. Use handleFlags() first to remove bad data."
    )
    pressure <- argo[["pressure"]]
    if (missing(p)) {
        pt <- apply(pressure, 1, median, na.rm = TRUE)
    } else if (length(p) == 1 && p == "levitus") {
        pt <- standardDepths()
        pt <- pt[pt < max(pressure, na.rm = TRUE)]
    } else if (is.numeric(p)) {
        if (length(p) == 1) {
            if (p < 1) {
                stop("\"p\" must exceed 1")
            }
            pt <- seq(0, max(pressure, na.rm = TRUE), length.out = p)
        } else {
            pt <- p
        }
    } else {
        stop("\"p\" must be numeric, or \"levitus\"")
    }
    npt <- length(pt)
    res@data$pressure <- matrix(NA, ncol = nprofile, nrow = npt)
    for (field in names(res@data)) {
        if (!(field %in% c("time", "longitude", "latitude"))) {
            res@data[[field]] <- matrix(NA, ncol = nprofile, nrow = npt)
            for (profile in 1:nprofile) {
                ndata <- sum(!is.na(argo@data[[field]][, profile]))
                if (ndata > 2 && sum(is.finite(diff(pressure[, profile]))) &&
                    0 < max(abs(diff(pressure[, profile])), na.rm = TRUE)) {
                    res@data[[field]][, profile] <- approx(pressure[, profile], argo@data[[field]][, profile], pt, ...)$y
                } else {
                    res@data[[field]][, profile] <- rep(NA, npt)
                }
                res@data$pressure[, profile] <- pt
            }
        }
    }
    res@processingLog <- processingLogAppend(res@processingLog, paste("Grid to regular pressures with: ", deparse(match.call()), sep = "", collapse = ""))
    for (w in warningMessages) {
        res@processingLog <- processingLogAppend(res@processingLog, w)
    }
    res
}

argoDecodeFlags <- function(f) # local function
{
    res <- unlist(lapply(seq_along(f), function(i) strsplit(f[i], split = "")))
    dim(res) <- c(length(res) / length(f), length(f))
    mode(res) <- "numeric"
    res
}


#' Read an Argo Data File
#'
#' `read.argo` is used to read an Argo file, producing an [argo-class] object.
#' The file must be in the ARGO-style NetCDF format described
#' in the Argo documentation (see references 2 and 3).
#'
#' @details
#'
#' See the Argo documentation (see references 2 and 3) for some details on what files contain.
#' Many items listed in section 2.2.3 of reference 3 are read from the
#' file and stored in the `metadata` slot, with the exception of
#' `longitude` and `latitude`, which are stored in the
#' `data` slot, alongside hydrographic information.  The details of storage
#' in the return value are somewhat complex, although the following notes
#' might be helpful to readers seeking to learn more.
#'
#' *1. Variable renaming.*
#'
#' The names of several data parameters stored within the netCDF file
#' are altered to fit the oce context. For example, `PRES` becomes `pressure`,
#' matching the name of this variable in other oce data types.
#' The original names are reported by `summary,argo-method`, and
#' data may be extracted with `[[,argo-method` using those names, so
#' the renaming should not be too inconvenient to Argo experts who
#' are new to oce.
#'
#' Argo netcdf files employ a `"SNAKE_CASE"` naming scheme (sometimes
#' using lower case) that is inconsistent with the `"camelCase"` scheme
#' used in oce. Since argo objects are just a small part of oce, a decision
#' was made to rename argo items. For example, `"CYCLE_NUMBER"` in the netcdf file
#' becomes `"cycleNumber"` in the oce object returned by `read.argo`.
#' (Note that `[[,argo-method` also accepts `"cycle"` for this item.)
#' The conversion for objects in the `data` slot often also involves
#' expanding on argo abbreviations, e.g. `"PSAL"` becomes `"salinity"`.
#' The renaming work is carried out with [argoNames2oceNames()] for
#' handles both name expansion for several dozen special cases,
#' and with [snakeToCamel()] with the `specialCase` argument
#' set to `"QC"`. While this results in variable names that should make
#' sense in the general oce context (where, for example, salinity is expected
#' to be stored in a variable named `"salinity"`), it may be confusing
#' to argo experts who are just starting to use oce.  Such people might
#' find it helpful to use e.g. `sort(names(x[["metadata"]]))` to get a list
#' of all items in the `metadata` slot (or similar with `"data"`), since working
#' in reverse may be easier than simply guessing at what names oce has chosen.
#' (Note that prior to 2020 June 24, some metadata items were stored in
#' `"SNAKE_CASE"`.)
#'
#' *2. Metadata.*
#'
#' Several of the netCDF global attributes are also renamed before
#' placement in the `metadata` slot of the return value.  These include
#' `conventions`, `featureType`, `history`, `institution`,
#' `nParameters`, `nProfiles`,  `references`, `source`, `title`,
#' and `userManualVersion`.
#' These names are derived from those in the netcdf
#' file, and mainly follow the pattern explained in the
#' \dQuote{Variable renaming convention} section.
#'
#' For profile data (as indicated by the NetCDF global attribute
#' named `"featureType"` being equal to `"trajectoryProfile"`),
#' the NetCDF item named `"STATION_PARAMETERS"` controls
#' whether variables in the source file will be stored in the
#' `metadata` or `data` slot of the returned object.
#' If `STATION_PARAMETERS` is not present, as is the case for
#' trajectory files (which are detected by `featureType` being
#' `"trajectory"`), some guesses are made as to what goes in
#' `data` and `metadata` slots.
#'
#' *3. Data variants.*
#'
#' Each data item can have variants, as
#' described in Sections 2.3.4 of reference 3.
#' For example, if `"PRES"` is found in `STATION_PARAMETERS`,
#' then `PRES` (pressure) data are sought in the file, along with
#' `PRES_QC`, `PRES_ADJUSTED`, `PRES_ADJUSTED_QC`, and
#' `PRES_ERROR`. The same pattern works for other profile data. The variables
#' are stored with names created as explained in the
#' \dQuote{Variable renaming convention} section below. Note that
#' flags, which are stored variables ending in `"_QC"` in the netcdf
#' file, are stored in the `flags` item within the `metadata` slot
#' of the returned object; thus, for example,
#' `PRES_QC` is stored as `pressure` in `flags`.
#'
#' *4. How time is handled.*
#'
#' The netcdf files for profile data store time in an item named `juld`,
#' which holds the overall profile time, in what the Argo documentation
#' calls Julian days, with respect to a reference time that is also stored
#' in the file.  Based on this information, a [POSIXct] value named `time`
#' is stored in the `metadata` slot of the returned value, and this
#' may be found with e.g. `a[["time"]]`, where `a` is that returned value.
#' Importantly, this value matches the time listed in profile index files.
#' In addition, some profile data files contain a field called `MTIME`,
#' which holds the offset (in days) between the time of individual measurements and the
#' overall profile time. For such files, the measurement times may be
#' computed with `a[["time"]]+86400*a[["mtime"]]`. (This formula is used by
#' [as.ctd()], if its first argument is an [argo-class] object created
#' by supplying [read.argo()] with such a data file.)
#'
#' *5. Data sources.*
#'
#' Argo data are made available at several websites. A bit of detective
#' work can be required to track down the data.
#'
#' Some servers provide data for floats that surfaced in a given ocean
#' on a given day, the anonymous FTP server
#' \code{usgodae.org/pub/outgoing/argo/geo/} being an example.
#'
#' Other servers provide data on a per-float basis. A complicating
#' factor is that these data tend to be categorized by "dac" (data
#' archiving centre), which makes it difficult to find a particular
#' float. For example,
#' \code{https://www.usgodae.org/ftp/outgoing/argo/} is the top level of
#' a such a repository. If the ID of a float is known but not the
#' "dac", then a first step is to download the text file
#' \code{https://www.usgodae.org/ftp/outgoing/argo/ar_index_global_meta.txt}
#' and search for the ID. The first few lines of that file are header,
#' and after that the format is simple, with columns separated by slash
#' (`/`). The dac is in the first such column and the float ID in the
#' second. A simple search will reveal the dac.
#' For example `data(argo)` is based on float 6900388, and the line
#' containing that token is
#' `bodc/6900388/6900388_meta.nc,846,BO,20120225005617`, from
#' which the dac is seen to be the British Oceanographic Data Centre
#' (`bodc`). Armed with that information, visit
#' \code{https://www.usgodae.org/ftp/outgoing/argo/dac/bodc/6900388}
#' and see a directory called `profiles` that contains a NetCDF
#' file for each profile the float made. These can be read with
#' `read.argo`. It is also possible, and probably more common,
#' to read a NetCDF file containing all the profiles together and for
#' that purpose the file
#' \code{https://www.usgodae.org/ftp/outgoing/argo/dac/bodc/6900388/6900388_prof.nc}
#' should be downloaded and provided as the `file` argument to
#' `read.argo`.  This can be automated as in Example 2,
#' although readers are cautioned that URL structures tend to change
#' over time.
#'
#' Similar steps can be followed on other servers.
#'
#'
#' @param file A character string giving the name of the file to load.
#'
#' @template encodingIgnoredTemplate
#'
#' @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 processingLog If provided, the action item to be stored in the log.
#' (Typically only provided for internal calls; the default that it provides is
#' better for normal calls by a user.)
#'
#' @param ... additional arguments, passed to called routines.
#'
#' @return
#' An [argo-class] object.
#'
#' @section Sample of Usage:
#' \preformatted{
#' # Example 1: read from a local file
#' library(oce)
#' d <- read.argo("~/data/OAR/6900388_prof.nc")
#' summary(d)
#' plot(d)
#'
#' # Example 2: construct URL for download (brittle)
#' id <- "6900388"
#' url <- "https://www.usgodae.org/ftp/outgoing/argo"
#' if (!length(list.files(pattern="argo_index.txt")))
#'     download.file(paste(url, "ar_index_global_meta.txt", sep="/"), "argo_index.txt")
#' index <- readLines("argo_index.txt")
#' line <- grep(id, index)
#' if (0 == length(line))
#'     stop("id ", id, " not found")
#' if (1 < length(line))
#'     stop("id ", id, " found multiple times")
#' dac <- strsplit(index[line], "/")[[1]][1]
#' profile <- paste(id, "_prof.nc", sep="")
#' float <- paste(url, "dac", dac, id, profile, sep="/")
#' download.file(float, profile)
#' argo <- read.argo(profile)
#' summary(argo)
#' }
#'
#'
#' @seealso
#' The documentation for the [argo-class] class explains the structure of argo
#' objects, and also outlines the other functions dealing with them.
#'
#' @references
#' 1. `https://argo.ucsd.edu`
#'
#' 2. Argo User's Manual Version 3.2, Dec 29th, 2015, available at
#' `https://archimer.ifremer.fr/doc/00187/29825/` online.
#'
#' 3. User's Manual (ar-um-02-01) 13 July 2010, available at
#' `http://www.argodatamgt.org/content/download/4729/34634/file/argo-dm-user-manual-version-2.3.pdf`,
#' which is the main document describing argo data.
#'
#'
#' @author Dan Kelley
#' @family things related to argo data
read.argo <- function(file, encoding = NA, 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, "\"")
        }
    }
    debug <- max(0, min(2, floor(as.numeric(debug))))
    if (!requireNamespace("ncdf4", quietly = TRUE)) {
        stop("must install.packages(\"ncdf4\") to read argo data")
    }
    if (missing(processingLog)) processingLog <- paste(deparse(match.call()), sep = "", collapse = "")
    # ofile <- file
    filename <- ""
    # NOTE: need to name ncdf4 package because otherwise R checks give warnings.
    if (is.character(file)) {
        filename <- fullFilename(file)
        file <- ncdf4::nc_open(file)
        on.exit(ncdf4::nc_close(file))
    } else {
        if (!inherits(file, "connection")) {
            stop("argument `file' must be a character string or connection")
        }
        if (!isOpen(file)) {
            file <- ncdf4::nc_open(file)
            on.exit(ncdf4::nc_close(file))
        }
    }
    oceDebug(debug, "read.argo(file=\"", filename, "\", ...) {\n", sep = "", unindent = 1, style = "bold")
    varNames <- names(file$var)
    # 'lc' will be TRUE if the data names are in lower case
    lc <- "data_type" %in% varNames
    oceDebug(debug, "File convention inferred to be ", if (lc) "lower-case" else "upper-case", ".\n", sep = "")
    res <- new("argo")
    # columnNames <- gsub(" *$", "", gsub("^ *", "", unique(as.vector(ncvar_get(f, maybeLC("STATION_PARAMETERS", lc))))))
    # QCNames <- paste(columnNames, "_QC",  sep="")
    # global attributes (see https://github.com/dankelley/oce/issues/1528)
    getGlobalAttribute <- function(file, attname) {
        a <- ncdf4::ncatt_get(nc = file, varid = 0, attname = attname)
        res <- if (a$hasatt) a$value else NULL
        # message("\"", attname, "\" value=\"", res, "\"")
        res
    }
    res@metadata$title <- getGlobalAttribute(file, "title")
    res@metadata$institution <- getGlobalAttribute(file, "institution")
    res@metadata$source <- getGlobalAttribute(file, "source")
    res@metadata$history <- getGlobalAttribute(file, "history")
    res@metadata$references <- getGlobalAttribute(file, "references")
    res@metadata$userManualVersion <- getGlobalAttribute(file, "user_manual_version")
    res@metadata$conventions <- getGlobalAttribute(file, "Conventions")
    res@metadata$featureType <- getGlobalAttribute(file, "featureType")
    varNamesOmit <- function(v, o) {
        where <- which(tolower(o) == tolower(v))
        if (length(where)) {
            v <- v[-where[1]]
        }
        v
    }
    oceDebug(debug - 1, "At processing step  1, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    res@metadata$id <- if (maybeLC("PLATFORM_NUMBER", lc) %in% varNames) {
        as.vector(trimws(ncdf4::ncvar_get(file, maybeLC("PLATFORM_NUMBER", lc))))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "PLATFORM_NUMBER")
    res@metadata$projectName <- if (maybeLC("PROJECT_NAME", lc) %in% varNames) {
        as.vector(trimws(ncdf4::ncvar_get(file, maybeLC("PROJECT_NAME", lc))))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "PROJECT_NAME")
    oceDebug(debug - 1, "Extracting PROJECT_NAME\n")
    oceDebug(debug - 1, "At processing step  2, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    res@metadata$PIName <- if (maybeLC("PI_NAME", lc) %in% varNames) {
        as.vector(trimws(ncdf4::ncvar_get(file, maybeLC("PI_NAME", lc))))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "PI_NAME")
    oceDebug(debug - 1, "Extracting PI_NAME\n")
    oceDebug(debug - 1, "At processing step  3, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    res@metadata$stationParameters <- NULL
    res@metadata$nParameters <- file$dim$N_PARAM$len
    res@metadata$nProfiles <- file$dim$N_PROF$len
    if (maybeLC("STATION_PARAMETERS", lc) %in% varNames) {
        res@metadata$stationParameters <- trimws(ncdf4::ncvar_get(file, maybeLC("STATION_PARAMETERS", lc)))
        if (is.null(res@metadata$stationParameters)) {
            warning(
                "This file has nothing listed in its STATION_PARAMETERS item, ",
                "so pressure, salinity, temperature, etc. are being stored in the metadata slot",
                "instead of the data slot. This will cause problems in further processing."
            )
        }
    } else {
        # print(sort(names(file$var)))
        if (res@metadata$featureType != "trajectory") {
            warning(
                "This 'profile'-type file lacks a STATION_PARAMETERS item, so guesses ",
                "were on whether to store items in the 'metadata' or 'data' slot. This may lead cause problems."
            )
        }
        for (want in c("PSAL", "TEMP", "PRES")) {
            if (want %in% toupper(varNames)) {
                res@metadata$stationParameters <- c(res@metadata$stationParameters, want)
                oceDebug(debug, "Will try to extract \"", want, "\" as a special case, because STATION_PARAMETERS is missing.\n", sep = "")
            }
        }
    }
    varNames <- varNamesOmit(varNames, "STATION_PARAMETERS")
    oceDebug(debug - 1, "Extracting STATION_PARAMETERS (if it exists)\n")
    oceDebug(debug - 1, "At processing step  4, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    oceDebug(debug - 1, "Extracting CYCLE_NUMBER\n")
    res@metadata$cycleNumber <- if (maybeLC("CYCLE_NUMBER", lc) %in% varNames) {
        as.vector(ncdf4::ncvar_get(file, maybeLC("CYCLE_NUMBER", lc)))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "CYCLE_NUMBER")
    oceDebug(debug - 1, "Extracting CYCLE_NUMBER_INDEX\n")
    res@metadata$cycleNumberIndex <- if (maybeLC("CYCLE_NUMBER_INDEX", lc) %in% varNames) {
        as.vector(ncdf4::ncvar_get(file, maybeLC("CYCLE_NUMBER_INDEX", lc)))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "CYCLE_NUMBER_INDEX")
    oceDebug(debug - 1, "Extracting CYCLE_NUMBER_ADJUSTED\n")
    res@metadata$cycleNumberAdjusted <- if (maybeLC("CYCLE_NUMBER_ADJUSTED", lc) %in% varNames) {
        as.vector(ncdf4::ncvar_get(file, maybeLC("CYCLE_NUMBER_ADJUSTED", lc)))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "CYCLE_NUMBER_ADJUSTED")
    oceDebug(debug - 1, "Extracting CYCLE_NUMBER_ADJUSTED_INDEX\n")
    res@metadata$cycleNumberAdjustedIndex <- if (maybeLC("CYCLE_NUMBER_ADJUSTED_INDEX", lc) %in% varNames) {
        as.vector(ncdf4::ncvar_get(file, maybeLC("CYCLE_NUMBER_ADJUSTED_INDEX", lc)))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "CYCLE_NUMBER_ADJUSTED_INDEX")
    oceDebug(debug - 1, "At processing step  5, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    res@metadata$direction <- if (maybeLC("DIRECTION", lc) %in% varNames) {
        as.vector(ncdf4::ncvar_get(file, maybeLC("DIRECTION", lc)))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "DIRECTION")
    oceDebug(debug - 1, "Extracting DIRECTION\n")
    oceDebug(debug - 1, "At processing step  6, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    res@metadata$dataCentre <- if (maybeLC("DATA_CENTRE", lc) %in% varNames) {
        as.vector(ncdf4::ncvar_get(file, maybeLC("DATA_CENTRE", lc)))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "DATA_CENTRE")
    oceDebug(debug - 1, "Extracting DATA_CENTRE\n")
    oceDebug(debug - 1, "At processing step  7, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    res@metadata$DCReference <- if (maybeLC("DC_REFERENCE", lc) %in% varNames) {
        as.vector(trimws(ncdf4::ncvar_get(file, maybeLC("DC_REFERENCE", lc))))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "DC_REFERENCE")
    oceDebug(debug - 1, "Extracting DC_REFERENCE\n")
    oceDebug(debug - 1, "At processing step  8, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    res@metadata$dataStateIndicator <- if (maybeLC("DATA_STATE_INDICATOR", lc) %in% varNames) {
        as.vector(trimws(ncdf4::ncvar_get(file, maybeLC("DATA_STATE_INDICATOR", lc))))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "DATA_STATE_INDICATOR")
    oceDebug(debug - 1, "Extracting DATA_STATE_INDICATOR\n")
    oceDebug(debug - 1, "At processing step  9, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    res@metadata$dataMode <- if (maybeLC("DATA_MODE", lc) %in% varNames) {
        strsplit(ncdf4::ncvar_get(file, maybeLC("DATA_MODE", lc)), "")[[1]]
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "DATA_MODE")
    oceDebug(debug - 1, "Extracting DATA_MODE\n")
    oceDebug(debug - 1, "At processing step 10, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    res@metadata$instReference <- if (maybeLC("INST_REFERENCE", lc) %in% varNames) {
        as.vector(trimws(ncdf4::ncvar_get(file, maybeLC("INST_REFERENCE", lc))))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "INST_REFERENCE")
    oceDebug(debug - 1, "Extracting INST_REFERENCE\n")
    oceDebug(debug - 1, "At processing step 11, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    res@metadata$firmwareVersion <- if (maybeLC("FIRMWARE_VERSION", lc) %in% varNames) {
        as.vector(trimws(ncdf4::ncvar_get(file, maybeLC("FIRMWARE_VERSION", lc))))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "FIRMWARE_VERSION")
    oceDebug(debug - 1, "Extracting FIRMWARE_REFERENCE\n")
    oceDebug(debug - 1, "At processing step 12, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    res@metadata$WMOInstType <- if (maybeLC("WMO_INST_TYPE", lc) %in% varNames) {
        as.vector(trimws(ncdf4::ncvar_get(file, maybeLC("WMO_INST_TYPE", lc))))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "WMO_INST_TYPE")
    oceDebug(debug - 1, "Extracting WMO_INST_TYPE\n")
    oceDebug(debug - 1, "At processing step 13, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    res@metadata$juld <- if (maybeLC("JULD", lc) %in% varNames) {
        as.vector(ncdf4::ncvar_get(file, maybeLC("JULD", lc)))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "JULD")
    oceDebug(debug - 1, "Extracting JULD\n")
    oceDebug(debug - 1, "At processing step 14, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    # set up 'time' also
    t0s <- as.vector(ncdf4::ncvar_get(file, maybeLC("REFERENCE_DATE_TIME", lc)))
    varNames <- varNamesOmit(varNames, "REFERENCE_DATE_TIME")
    oceDebug(debug - 1, "Extracting REFERENCE_DATE_TIME\n")
    oceDebug(debug - 1, "At processing step 15, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    t0 <- strptime(t0s, "%Y%m%d%M%H%S", tz = "UTC")
    oceDebug(debug, vectorShow(t0))
    oceDebug(debug, vectorShow(res@metadata$juld))
    # julianDayTime <- as.vector(ncdf4::ncvar_get(file, maybeLC("JULD", lc)))
    res@metadata$time <- t0 + res@metadata$juld * 86400
    oceDebug(debug, vectorShow(res@metadata$time))
    rm(list = c("t0s", "t0")) # no longer needed
    res@metadata$juldQC <- if (maybeLC("JULD_QC", lc) %in% varNames) {
        as.vector(ncdf4::ncvar_get(file, maybeLC("JULD_QC", lc)))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "JULD_QC")
    oceDebug(debug - 1, "Extracting JULD_QC\n")
    oceDebug(debug - 1, "At processing step 16, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    res@metadata$juldLocation <- if (maybeLC("JULD_LOCATION", lc) %in% varNames) {
        as.vector(ncdf4::ncvar_get(file, maybeLC("JULD_LOCATION", lc)))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "JULD_LOCATION")
    oceDebug(debug - 1, "Extracting JULD_LOCATION\n")
    oceDebug(debug - 1, "At processing step 17, the ", length(varNames),
        " varnames are: c(\"", paste(sort(varNames), collapse = "\",\""), "\")\n",
        sep = ""
    )
    # Now for the data.
    res@metadata$dataNamesOriginal <- list() # NB. will store upper-case names
    if (maybeLC("LATITUDE", lc) %in% varNames) {
        res@data$latitude <- as.vector(ncdf4::ncvar_get(file, maybeLC("LATITUDE", lc)))
        res@metadata$dataNamesOriginal$latitude <- "LATITUDE"
        varNames <- varNamesOmit(varNames, "LATITUDE")
        latitudeNA <- ncdf4::ncatt_get(file, maybeLC("LATITUDE", lc), "_FillValue")$value
        res@data$latitude[res@data$latitude == latitudeNA] <- NA
        rm(list = "latitudeNA") # no longer needed
        res@metadata$units$latitude <- if (1 == grepl("north", ncdf4::ncatt_get(file, maybeLC("LATITUDE", lc), "units")$value, ignore.case = TRUE)) {
            list(unit = expression(degree * N), scale = "")
        } else {
            list(unit = expression(degree * S), scale = "")
        }
    }
    if (maybeLC("LONGITUDE", lc) %in% varNames) {
        res@data$longitude <- as.vector(ncdf4::ncvar_get(file, maybeLC("LONGITUDE", lc)))
        res@metadata$dataNamesOriginal$longitude <- "LONGITUDE"
        varNames <- varNamesOmit(varNames, "LONGITUDE")
        longitudeNA <- ncdf4::ncatt_get(file, maybeLC("LONGITUDE", lc), "_FillValue")$value
        res@data$longitude[res@data$longitude == longitudeNA] <- NA
        rm(list = "longitudeNA") # no longer needed
        res@metadata$units$longitude <- if (1 == length(grep("east",
            ncdf4::ncatt_get(file, maybeLC("LONGITUDE", lc), "units")$value,
            ignore.case = TRUE
        ))) {
            list(unit = expression(degree * E), scale = "")
        } else {
            list(unit = expression(degree * W), scale = "")
        }
    }
    if (maybeLC("MTIME", lc) %in% varNames) {
        res@data$mtime <- ncdf4::ncvar_get(file, maybeLC("MTIME", lc))
        res@metadata$dataNamesOriginal$mtime <- "MTIME"
        varNames <- varNamesOmit(varNames, "MTIME")
        res@metadata$units$mtime <- list(unit = expression(day), scale = "")
    }
    res@metadata$positionQC <- if (maybeLC("POSITION_QC", lc) %in% varNames) {
        as.vector(ncdf4::ncvar_get(file, maybeLC("POSITION_QC", lc)))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "POSITION_QC")
    oceDebug(debug - 1, "POSITION_QC\n")
    oceDebug(debug - 1, "varNames=", paste(varNames, collapse = ","), "\n")
    res@metadata$positioningSystem <- if (maybeLC("POSITIONING_SYSTEM", lc) %in% varNames) {
        as.vector(trimws(ncdf4::ncvar_get(file, maybeLC("POSITIONING_SYSTEM", lc))))
    } else {
        NULL
    }
    varNames <- varNamesOmit(varNames, "POSITIONING_SYSTEM")
    oceDebug(debug - 1, "POSITIONING_SYSTEM\n")
    oceDebug(debug - 1, "varNames=", paste(varNames, collapse = ","), "\n")
    stationParameters <- unique(as.vector(res@metadata$stationParameters)) # will be PRES, TEMP etc
    # Handle units before getting the data. (We must do it here, because when
    # we read the data, we remove entries from varNames ... that is how
    # we know what is left over at the end, to be shoved into metadata.)
    # FIXME: it would be nice to automate this more, to handle cases not handled
    # in the hardwired code below.
    # Oxygen: DOXY DOXY_ADJUSTED DOXY_ADJUSTED_ERROR
    if (maybeLC("DOXY", lc) %in% varNames) {
        attTMP <- ncdf4::ncatt_get(file, maybeLC("DOXY", lc), "units")
        if (attTMP$hasatt) {
            if (attTMP$value == "micromole/kg") {
                res@metadata$units$oxygen <- list(unit = expression(mu * mol / kg), scale = "")
            } else {
                warning("skipping oxygen unit '", attTMP$value, "' because only understood unit is 'micromole/kg'", sep = "")
            }
        }
    }
    if (maybeLC("DOXY_ADJUSTED", lc) %in% varNames) {
        # print(ncdf4::ncatt_get(file, maybeLC("DOXY", lc), "long_name"))
        attTMP <- ncdf4::ncatt_get(file, maybeLC("DOXY_ADJUSTED", lc), "units")
        if (attTMP$hasatt) {
            if (attTMP$value == "micromole/kg") {
                res@metadata$units$oxygenAdjusted <- list(unit = expression(mu * mol / kg), scale = "")
            } else {
                warning("skipping oxygenAdjusted unit '", attTMP$value, "' because only understood unit is 'micromole/kg'", sep = "")
            }
        }
    }
    if (maybeLC("DOXY_ADJUSTED_ERROR", lc) %in% varNames) {
        # print(ncdf4::ncatt_get(file, maybeLC("DOXY", lc), "long_name"))
        attTMP <- ncdf4::ncatt_get(file, maybeLC("DOXY_ADJUSTED_ERROR", lc), "units")
        if (attTMP$hasatt) {
            if (attTMP$value == "micromole/kg") {
                res@metadata$units$oxygenAdjustedError <- list(unit = expression(mu * mol / kg), scale = "")
            } else {
                warning("skipping oxygenAdjustedError unit '", attTMP$value,
                    "' because only understood unit is 'micromole/kg'",
                    sep = ""
                )
            }
        }
    }
    # Temperature: TEMP, TEMP_ADJUSTED, TEMP_ADJUSTED_ERROR
    if (maybeLC("TEMP", lc) %in% varNames) {
        if (1 == length(grep("ITS-90", ncdf4::ncatt_get(file, maybeLC("TEMP", lc), "long_name")$value, ignore.case = TRUE))) {
            res@metadata$units$temperature <- list(unit = expression(degree * C), scale = "ITS-90")
        } else {
            res@metadata$units$temperature <- list(unit = expression(degree * C), scale = "ITS-90")
        }
    }
    if (maybeLC("TEMP_ADJUSTED", lc) %in% varNames) {
        if (1 == length(grep("ITS-90", ncdf4::ncatt_get(file, maybeLC("TEMP_ADJUSTED", lc), "long_name")$value, ignore.case = TRUE))) {
            res@metadata$units$temperatureAdjusted <- list(unit = expression(degree * C), scale = "ITS-90")
        } else {
            res@metadata$units$temperatureAdjusted <- list(unit = expression(degree * C), scale = "ITS-90")
        }
    }
    if (maybeLC("TEMP_ADJUSTED_ERROR", lc) %in% varNames) {
        if (1 == length(grep("ITS-90", ncdf4::ncatt_get(file, maybeLC("TEMP_ADJUSTED_ERROR", lc), "long_name")$value, ignore.case = TRUE))) {
            res@metadata$units$temperatureAdjustedError <- list(unit = expression(degree * C), scale = "ITS-90")
        } else {
            res@metadata$units$temperatureAdjustedError <- list(unit = expression(degree * C), scale = "ITS-90")
        }
    }
    # salinity: PSAL, PSAL_ADJUSTED, PSAL_ADJUSTED_ERROR
    if (maybeLC("PSAL", lc) %in% varNames) {
        if (1 == length(grep("PRACTICAL", ncdf4::ncatt_get(file, maybeLC("PSAL", lc), "long_name")$value, ignore.case = TRUE))) {
            res@metadata$units$salinity <- list(unit = expression(), scale = "PSS-78")
        } else {
            res@metadata$units$salinity <- list(unit = expression(), scale = "PSS-78")
        }
    }
    if (maybeLC("PSAL_ADJUSTED", lc) %in% varNames) {
        if (1 == length(grep("PRACTICAL", ncdf4::ncatt_get(file, maybeLC("PSAL_ADJUSTED", lc), "long_name")$value, ignore.case = TRUE))) {
            res@metadata$units$salinityAdjusted <- list(unit = expression(), scale = "PSS-78")
        } else {
            res@metadata$units$salinityAdjusted <- list(unit = expression(), scale = "PSS-78")
        }
    }
    if (maybeLC("PSAL_ADJUSTED_ERROR", lc) %in% varNames) {
        if (1 == length(grep(maybeLC("PRACTICAL", lc),
            ncdf4::ncatt_get(file, maybeLC("PSAL_ADJUSTED_ERROR", lc), "long_name")$value,
            ignore.case = TRUE
        ))) {
            res@metadata$units$salinityAdjustedError <- list(unit = expression(), scale = "PSS-78")
        } else {
            res@metadata$units$salinityAdjustedError <- list(unit = expression(), scale = "PSS-78")
        }
    }
    # pressure: PRES, PRES_ADJUSTED, PRES_ADJUSTED_ERROR
    if (maybeLC("PRES", lc) %in% varNames) {
        if (1 == length(grep("decibar", ncdf4::ncatt_get(file, maybeLC("PRES", lc), "units")$value, ignore.case = TRUE))) {
            res@metadata$units$pressure <- list(unit = expression(dbar), scale = "")
        } else {
            res@metadata$units$pressure <- list(unit = expression(dbar), scale = "")
        }
    }
    if (maybeLC("PRES_ADJUSTED", lc) %in% varNames) {
        if (1 == length(grep("decibar", ncdf4::ncatt_get(file, maybeLC("PRES_ADJUSTED", lc), "units")$value, ignore.case = TRUE))) {
            res@metadata$units$pressureAdjusted <- list(unit = expression(dbar), scale = "")
        } else {
            res@metadata$units$pressureAdjusted <- list(unit = expression(dbar), scale = "")
        }
    }
    if (maybeLC("PRES_ADJUSTED_ERROR", lc) %in% varNames) {
        if (1 == length(grep("decibar", ncdf4::ncatt_get(file, maybeLC("PRES_ADJUSTED_ERROR", lc), "units")$value, ignore.case = TRUE))) {
            res@metadata$units$pressureAdjustedError <- list(unit = expression(dbar), scale = "")
        } else {
            res@metadata$units$pressureAdjustedError <- list(unit = expression(dbar), scale = "")
        }
    }
    # Fix up names of flags. This became required with changes made to argoNames2oceNames() in Dec 17-18, 2016. Arguably, I
    # should find out why the change occurred, but fixing the names now is just as easy, and might be clearer to the reader.
    names(res@metadata$flags) <- gsub("QC$", "", names(res@metadata$flags))
    # Now, work through the columnar data. Note how we remove entries from varNames again,
    # after having interrupted that practice whilst finding units and flags.
    # Skip MTIME, if it exists, since it is handled separately.
    if ("MTIME" %in% stationParameters) {
        stationParameters <- stationParameters[stationParameters != "MTIME"]
    }
    oceDebug(debug, "About to process stationParameters: c(\"",
        paste(stationParameters, collapse = "\",\""), "\")\n",
        sep = ""
    )
    for (item in stationParameters) {
        # some files have unnamed variables, so we skip them
        if (!nchar(item)) {
            next
        }
        n <- item
        d <- getData(file, maybeLC(n, lc))
        varNames <- varNamesOmit(varNames, n)
        if (!is.null(d)) {
            oceDebug(debug, "  Storing \"", n, "\" as \"", argoNames2oceNames(n), "\" in the data slot.\n", sep = "")
            res@data[[argoNames2oceNames(n)]] <- d
            res@metadata$dataNamesOriginal[[argoNames2oceNames(n)]] <- n
        } else {
            oceDebug(debug, "  Set item = \"", n, "\" in data slot to NULL, since the data file contains no data for this.\n", sep = "")
            res@data[[argoNames2oceNames(n)]] <- NULL
        }
        oceDebug(debug - 1, "  Remaining ", length(varNames), "are: =", paste(varNames, collapse = " "), "\n")
        n <- paste(item, maybeLC("_QC", lc), sep = "")
        oceDebug(debug - 2, "  about to try to get '", n, "' from netcdf file\n", sep = "")
        d <- getData(file, maybeLC(n, lc), quiet = TRUE)
        oceDebug(debug - 2, "  ... got it\n", sep = "")
        varNames <- varNamesOmit(varNames, n)
        oceDebug(debug - 1, n, "\n")
        oceDebug(debug - 1, "  B varNames=", paste(sort(varNames), collapse = ","), "\n")
        if (!is.null(d)) res@metadata$flags[[argoNames2oceNames(n)]] <- argoDecodeFlags(d)
        n <- paste(item, maybeLC("_ADJUSTED", lc), sep = "")
        if (n %in% varNames) {
            d <- getData(file, maybeLC(n, lc))
            varNames <- varNamesOmit(varNames, n)
            oceDebug(debug - 1, n, "\n")
            oceDebug(debug - 1, "C varNames=", paste(sort(varNames), collapse = ","), "\n")
            if (!is.null(d)) {
                res@data[[argoNames2oceNames(n)]] <- d
                res@metadata$dataNamesOriginal[[argoNames2oceNames(n)]] <- n
            } else {
                res@data[[argoNames2oceNames(n)]] <- NULL
            }
        }
        n <- paste(item, maybeLC("_ADJUSTED_QC", lc), sep = "")
        if (n %in% varNames) {
            d <- getData(file, maybeLC(n, lc))
            varNames <- varNamesOmit(varNames, n)
            oceDebug(debug - 1, n, "\n")
            oceDebug(debug - 1, "D varNames=", paste(sort(varNames), collapse = ","), "\n")
            if (!is.null(d)) res@metadata$flags[[argoNames2oceNames(n)]] <- argoDecodeFlags(d)
        }
        n <- paste(item, maybeLC("_ADJUSTED_ERROR", lc), sep = "")
        if (n %in% varNames) {
            d <- getData(file, maybeLC(n, lc))
            varNames <- varNamesOmit(varNames, n)
            oceDebug(debug - 1, n, "\n")
            oceDebug(debug - 1, "E varNames=", paste(sort(varNames), collapse = ","), "\n")
            if (!is.null(d)) {
                res@data[[argoNames2oceNames(n)]] <- d
                res@metadata$dataNamesOriginal[[argoNames2oceNames(n)]] <- n
            } else {
                res@data[[argoNames2oceNames(n)]] <- NULL
            }
        }
    }
    oceDebug(debug, "After processing stationParameters, flag names are: c(\"",
        paste(names(res@metadata$flags), collapse = "\",\""), "\").\n",
        sep = ""
    )
    if (length(res@metadata$flags)) {
        names(res@metadata$flags) <- gsub("QC$", "", names(res@metadata$flags))
    }
    oceDebug(debug, "After trimming QC, flag names are: c(\"",
        paste(names(res@metadata$flags), collapse = "\",\""), "\")\n",
        sep = ""
    )
    res@metadata$filename <- filename
    # cat("@L1391 varnames: ", paste(sort(varNames), sep=" "), "\n")
    # Now, insert any unprocessed items from varNames into metadata. We
    # need to check for access failures because we get the error
    #     > Error in R_nc4_get_vara_text: NetCDF: Index exceeds dimension bound
    #     > Var: HISTORY_PARAMETER  Ndims: 3   Start: 0,0,0 Count: 0,223,16
    # for "HISTORY_INSTITUTION", "HISTORY_STEP", "HISTORY_SOFTWARE",
    # "HISTORY_SOFTWARE_RELEASE", "HISTORY_REFERENCE",
    # "HISTORY_DATE", "HISTORY_ACTION", "HISTORY_PARAMETER",
    # "HISTORY_START_PRES", "HISTORY_STOP_PRES",
    # "HISTORY_PREVIOUS_VALUE", "HISTORY_QCTEST"
    for (name in varNames) {
        ocename <- snakeToCamel(name, specialCases = c("QC"))
        oceDebug(debug, "Inserting \"", name, "\" as \"", ocename, "\" in the metadata slot.\n", sep = "")
        value <- NA
        o <- capture.output({
            value <- try(ncdf4::ncvar_get(file, name), silent = TRUE)
        })
        if (inherits(value, "try-error")) {
            # see https://github.com/dankelley/oce/issues/1522 for a discussion of the fact
            # that the file used for data(argo) has zero-length HISTORY_* items.
            if (length(grep("Index exceeds dimension", o))) {
                # FIXME: this code is brittle, being dependent on the layout of
                # FIXME: the output from nc_open(), which might be subject to change.
                # FIXME: The code worked on 2019-04-13.
                oceDebug(
                    debug, "ncdf4::ncvar_get() error diagnosis ... name=", name, " len=",
                    file$var[[name]]$dim[[file$var[[name]]$ndim]]$len,
                    " (if this is 0, will not save '", name, "' to metadata)\n"
                )
                if (0 != file$var[[name]]$dim[[file$var[[name]]$ndim]]$len) {
                    oceDebug(debug, "ncvar_get() failed for \"", name, "\" (Index exceeds dimension), so it isn't stored in metadata\n")
                }
            } else {
                oceDebug(debug, "ncvar_get() failed for \"", name, "\", so it isn't stored in metadata\n")
            }
        } else {
            # Make a vector, if it is a single-column matrix
            if (1 == length(dim(value))) {
                value <- as.vector(value)
            }
            # Trim leading/trailing whitespace, if it is a string
            if (is.character(value)) {
                origValue <- value
                value <- try(
                    {
                        trimws(value)
                    },
                    silent = TRUE
                )
                if (inherits(value, "try-error")) {
                    warning("cannot trim leading/trailing whitespace in metadata$", ocename, "")
                    value <- origValue
                }
            }
            res@metadata[[ocename]] <- value
        }
    }
    # Record a log item
    res@processingLog <- processingLogAppend(res@processingLog, paste("read.argo(file=\"", filename, "\")", sep = ""))
    oceDebug(debug, "} # read.argo()\n", sep = "", unindent = 1, style = "bold")
    res
}

#' Coerce Data Into an argo Object
#'
#' Coerce a dataset into an argo dataset. This is not the right way to
#' read official argo datasets, which are provided in NetCDF format and may
#' be read with [read.argo()].
#'
#' @param time a vector of POSIXct times.
#' @param longitude a vector of longitudes.
#' @param latitude a vector of latitudes.
#' @param salinity a vector of salinities.
#' @param temperature a vector of temperatures.
#' @param pressure a vector of pressures.
#' @param units an optional list containing units. If `NULL`, the default,
#' then `"degree east"` is used for `longitude`,
#' `"degree north"` for `latitude`,
#' `""` for `salinity`,
#' `"ITS-90"` for `temperature`, and
#' `"dbar"` for `pressure`.
#' @param id an identifier for the argo float, typically a number, but stored within
#' the object in a character form. (For example, the dataset retrieved with `data(argo)`
#' has an `id` of `"6900388"`.)
#' @param filename a source filename, which defaults to an empty string.
#' @param missingValue an optional missing value, indicating data values that should be
#' taken as `NA`.
#'
#' @return
#' An [argo-class] object.
#'
#' @seealso
#' The documentation for the [argo-class] class explains the structure of argo
#' objects, and also outlines the other functions dealing with them.
#'
#' @author Dan Kelley
#' @family things related to argo data
as.argo <- function(time, longitude, latitude, salinity, temperature, pressure, units = NULL, id, filename = "", missingValue) {
    if (inherits(class, "data.frame")) {
        df <- time
        names <- names(df)
        time <- if ("time" %in% names) df$time else NULL
        salinity <- if ("salinity" %in% names) df$salinity else NULL
        temperature <- if ("temperature" %in% names) df$temperature else NULL
        pressure <- if ("pressure" %in% names) df$pressure else NULL
        longitude <- if ("longitude" %in% names) df$longitude else NULL
        latitude <- if ("latitude" %in% names) df$latitude else NULL
        id <- if ("id" %in% names) df$id else NULL
    } else {
        if (missing(time)) {
            stop("must give time")
        }
        if (missing(longitude)) {
            stop("must give longitude")
        }
        if (missing(latitude)) {
            stop("must give latitude")
        }
        if (missing(temperature)) {
            stop("must give temperature")
        }
        if (missing(salinity)) {
            stop("must give salinity")
        }
        if (missing(pressure)) {
            stop("must give pressure")
        }
        if (missing(id)) {
            stop("must give id")
        }
    }
    res <- new("argo",
        time = time, id = id,
        longitude = longitude, latitude = latitude, salinity = salinity,
        temperature = temperature, pressure = pressure, filename = filename
    )
    res@metadata$units <- if (!is.null(units)) {
        units
    } else {
        list(
            longitude = list(expression(degree * E), scale = ""),
            latitude = list(expression(degree * N), scale = ""),
            salinity = list(unit = expression(), scale = "PSS-78"), # assuming a particular scale
            temperature = list(unit = expression(degree * C), scale = "ITS-90"), # assuming a particular scale
            pressure = list(unit = expression(dbar), scale = "")
        ) # assuming a particular unit
    }
    res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep = "", collapse = ""))
    res
}


#' Plot an argo Object
#'
#' Plot a summary diagram for argo data.
#'
#' @param x an [argo-class] object.
#'
#' @param which list of desired plot types, one of the following. Note
#' that [oce.pmatch()] is used to try to complete partial
#' character matches, and that an error will occur if the match is
#' not complete (e.g. `"salinity"` matches to both
#' `"salinity ts"` and `"salinity profile"`.).
#'
#' * `which=1`, `which="trajectory"` or `which="map"` gives a
#' plot of the argo trajectory, with the coastline, if one is provided.
#'
#' * `which=2` or `"salinity ts"` gives a time series of
#' salinity at the indicated level(s)
#'
#' * `which=3` or `"temperature ts"` gives a time series
#' of temperature at the indicated level(s)
#'
#' * `which=4` or `"TS"` gives a TS diagram at the
#' indicated level(s)
#'
#' * `which=5` or `"salinity profile"` gives a salinity profile
#'
#' * `which=6` or `"temperature profile"` gives a temperature profile
#'
#' * `which=7` or `"sigma0 profile"` gives a sigma0 profile
#'
#' * `which=8` or `"spice profile"` gives a spiciness
#' profile, referenced to the surface.  (This is the
#' same as using `which=9`.)
#'
#' * `which=9` or `"spiciness0 profile"` gives
#' a profile of spicininess referenced to a pressure
#' of 0 dbar, i.e. the surface.  (This is the
#' same as using `which=8`.)
#'
#' * `which=10` or `"spiciness1 profile"` gives
#' a profile of spiciness referenced to a pressure
#' of 1000 dbar.
#'
#' * `which=11` or `"spiciness2 profile"` gives
#' a profile of spiciness referenced to a pressure
#' of 2000 dbar.
#'
#' @param level depth pseudo-level to plot, for `which=2` and higher.
#' May be an integer, in which case it refers to an index of depth (1
#' being the top) or it may be the string "all" which means to plot
#' all data.
#'
#' @param coastline character string giving the coastline to be used
#' in an Argo-location map, or `"best"` to pick the one with highest
#' resolution, or `"none"` to avoid drawing the coastline.
#'
#' @param cex size of plotting symbols to be used if `type="p"`.
#'
#' @param pch type of plotting symbols to be used if `type="p"`.
#'
#' @param type plot type, either `"l"` or `"p"`.
#'
#' @param col optional list of colors for plotting.
#'
#' @param fill either a logical, indicating whether to fill the land with
#' light-gray, or a color name.  Owing to problems with some projections, the
#' default is not to fill.
#'
#' @param mgp a 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 projection character value indicating the projection to be used
#' in trajectory maps. If this is `NULL`, no projection is used, although
#' the plot aspect ratio will be set to yield zero shape distortion at the
#' mean float latitude.  If `projection="automatic"`, then one
#' of two projections is used: stereopolar (i.e. `"+proj=stere +lon_0=X"`
#' where `X` is the mean longitude), or Mercator (i.e. `"+proj=merc"`)
#' otherwise.  Otherwise, `projection` must be a character string specifying
#' a projection in the notation used by [oceProject()] and [mapPlot()].
#'
#' @param mar value to be used with `par("mar")`.
#'
#' @param tformat optional argument passed to [oce.plot.ts()], for plot
#' types that call that function.  (See [strptime()] for the format
#' used.)
#'
#' @param debug debugging flag.
#'
#' @param ... optional arguments passed to plotting functions.
#'
#' @return None.
#'
#' @examples
#' library(oce)
#' data(argo)
#' tc <- cut(argo[["time"]], "year")
#' # Example 1: plot map, which reveals float trajectory.
#' plot(argo, pch = as.integer(tc))
#' year <- substr(levels(tc), 1, 4)
#' data(topoWorld)
#' contour(topoWorld[["longitude"]], topoWorld[["latitude"]],
#'     topoWorld[["z"]],
#'     add = TRUE
#' )
#' legend("bottomleft", pch = seq_along(year), legend = year, bg = "white", cex = 3 / 4)
#'
#' # Example 2: plot map, TS, T(z) and S(z). Note the use
#' # of handleFlags(), to skip over questionable data.
#' plot(handleFlags(argo), which = c(1, 4, 6, 5))
#'
#' @family things related to argo data
#' @family functions that plot oce data
#'
#' @aliases plot.argo
#'
#' @author Dan Kelley
setMethod(
    f = "plot",
    signature = signature("argo"),
    definition = function(x, which = 1, level,
                          coastline = c("best", "coastlineWorld", "coastlineWorldMedium", "coastlineWorldFine", "none"),
                          cex = 1, pch = 1, type = "p", col = 1, fill = FALSE,
                          projection = NULL, mgp = getOption("oceMgp"),
                          mar = c(mgp[1] + 1.5, mgp[1] + 1.5, 1.5, 1.5),
                          tformat, debug = getOption("oceDebug"), ...) {
        debug <- min(3L, max(0L, as.integer(debug)))
        if (!inherits(x, "argo")) {
            stop("method is only for objects of class '", "argo", "'")
        }
        oceDebug(debug, "plot.argo(x, which=c(", paste(which, collapse = ","), "),",
            argShow(mgp),
            argShow(mar),
            argShow(cex),
            " ...) {\n",
            sep = "", unindent = 1, style = "bold"
        )
        coastline <- match.arg(coastline)
        nw <- length(which)
        if (nw > 1) {
            par(mfcol = c(1, nw))
        }
        par(mgp = mgp, mar = mar)
        if (missing(level) || level == "all") {
            level <- seq(1L, dim(x@data$temperature)[1])
        }
        longitude <- x[["longitude"]]
        latitude <- x[["latitude"]]
        dim <- dim(x@data$salinity)
        if (length(longitude) < prod(dim)) {
            # Copy across depths. This is inside a conditional because possibly
            # argo[["longitude"]] should mimic section[["longitude"]], in doing
            # the lengthing by itself unless the second argument is "byStation"
            # (issue 1273 ... under consideration 2017jul12)
            longitude <- rep(x[["longitude"]], each = dim[1])
            latitude <- rep(x[["latitude"]], each = dim[1])
        }
        ctd <- as.ctd(x@data$salinity, x@data$temperature, x@data$pressure,
            longitude = longitude, latitude = latitude,
            units = list(
                temperature = list(unit = expression(degree * C), scale = "ITS-90"),
                conductivity = list(list = expression(), scale = "")
            )
        ) # guess on units
        whichOrig <- which
        which <- oce.pmatch(
            which,
            list(
                "trajectory" = 1,
                "map" = 1,
                "salinity ts" = 2,
                "temperature ts" = 3,
                "TS" = 4,
                "salinity profile" = 5,
                "temperature profile" = 6,
                "sigma0 profile" = 7,
                "spice profile" = 8,
                "spiciness0 profile" = 9,
                "spiciness1 profile" = 10,
                "spiciness2 profile" = 11
            )
        )
        # if (any(is.na(which)))
        #    stop("In plot,argo-method() :\n  unrecognized value(s) of which: ", paste(whichOrig[is.na(which)], collapse=", "), call.=FALSE)
        lw <- length(which)
        par(mgp = mgp)
        if (lw == 2) {
            par(mfcol = c(2, 1))
        } else if (lw == 3) {
            par(mfcol = c(3, 1))
        } else if (lw == 4) {
            par(mfrow = c(2, 2))
        } else if (lw != 1) {
            nnn <- floor(sqrt(lw))
            par(mfcol = c(nnn, ceiling(lw / nnn)))
            rm(nnn)
        }
        for (w in 1:nw) {
            oceDebug(debug, "handling which[", w, "]=\"", which[w], "\"\n", sep = "")
            if (is.na(which[w])) {
                oceDebug(debug, "not a special case, so passing 'which' to plot,ctd-method\n")
                plot(ctd, which = whichOrig[w], debug = debug - 1, ...)
            } else if (which[w] == 1) {
                oceDebug(debug, "note: par(\"mfrow\") = ", paste(par("mfrow"), collapse = ","), "\n")
                # map
                # FIXME: coastline selection should be DRY
                haveCoastline <- FALSE
                haveOcedata <- requireNamespace("ocedata", quietly = TRUE)
                lonr <- range(x[["longitude"]], na.rm = TRUE)
                latr <- range(x[["latitude"]], na.rm = TRUE)
                if (coastline == "best") {
                    if (haveOcedata) {
                        bestcoastline <- coastlineBest(lonRange = lonr, latRange = latr)
                        oceDebug(debug, " 'best' coastline is: \"", bestcoastline, '\"\n', sep = "")
                        if (bestcoastline == "coastlineWorld") {
                            data(list = bestcoastline, package = "oce", envir = environment())
                        } else {
                            data(list = bestcoastline, package = "ocedata", envir = environment())
                        }
                        coastline <- get(bestcoastline)
                    } else {
                        bestcoastline <- coastlineBest(lonRange = lonr, latRange = latr)
                        oceDebug(debug, " using \"coastlineWorld\" because ocedata package not installed\n")
                        data("coastlineWorld", package = "oce", envir = environment())
                        coastline <- get("coastlineWorld")
                    }
                    haveCoastline <- TRUE
                } else {
                    if (coastline != "none") {
                        if (coastline == "coastlineWorld") {
                            data("coastlineWorld", package = "oce", envir = environment())
                            coastline <- get("coastlineWorld")
                        } else if (haveOcedata && coastline == "coastlineWorldFine") {
                            data("coastlineWorldFine", package = "ocedata", envir = environment())
                            coastline <- get("coastlineWorldFine")
                        } else if (haveOcedata && coastline == "coastlineWorldMedium") {
                            data("coastlineWorldMedium", package = "ocedata", envir = environment())
                            coastline <- get("coastlineWorldMedium")
                        } else {
                            stop("there is no built-in coastline file of name \"", coastline, "\"")
                        }
                        haveCoastline <- TRUE
                    }
                }
                if (!is.null(projection)) {
                    oceDebug(debug, "drawing an argo map with a projection\n")
                    meanlat <- mean(x[["latitude"]], na.rm = TRUE)
                    meanlon <- mean(x[["longitude"]], na.rm = TRUE)
                    # id <- pmatch(projection, "automatic")
                    if (!is.na(pmatch(projection, "automatic"))) {
                        projection <- if (meanlat > 70) {
                            paste("+proj=stere +lon_0=", meanlon, sep = "")
                        } else {
                            "+proj=merc"
                        }
                        oceDebug(debug, "using", projection, "projection (chosen automatically)\n")
                    } else {
                        oceDebug(debug, "using", projection, "projection (specified)\n")
                    }
                    mapPlot(x[["longitude"]], x[["latitude"]],
                        projection = projection, type = "p", cex = cex, pch = pch, col = col, debug = debug - 1
                    )
                    if (is.logical(fill) && fill) {
                        mapPolygon(coastline[["longitude"]], coastline[["latitude"]], col = "lightgray")
                    } else {
                        if (is.character(fill)) {
                            mapPolygon(coastline[["longitude"]], coastline[["latitude"]], col = fill)
                        } else {
                            mapPolygon(coastline[["longitude"]], coastline[["latitude"]])
                        }
                    }
                } else {
                    oceDebug(debug, "drawing an argo map without a projection\n")
                    asp <- 1 / cos(mean(range(x@data$latitude, na.rm = TRUE)) * atan2(1, 1) / 45)
                    plot(x@data$longitude, x@data$latitude,
                        asp = asp,
                        type = type, cex = cex, pch = pch, col = col,
                        xlab = resizableLabel("longitude"), ylab = resizableLabel("latitude"), ...
                    )
                    if (haveCoastline) {
                        if (!is.null(coastline@metadata$fillable) && coastline@metadata$fillable) {
                            polygon(coastline[["longitude"]], coastline[["latitude"]], col = "lightgray", lwd = 3 / 4)
                            polygon(coastline[["longitude"]] + 360, coastline[["latitude"]], col = "lightgray", lwd = 3 / 4)
                        } else {
                            lines(coastline[["longitude"]], coastline[["latitude"]], col = "darkgray")
                            lines(coastline[["longitude"]] + 360, coastline[["latitude"]], col = "darkgray")
                        }
                    }
                    if (!missing(coastline)) {
                        polygon(coastline[["longitude"]], coastline[["latitude"]], col = "lightgray")
                        if (type[w] == "l") {
                            lines(x@data$longitude, x@data$latitude)
                        } else {
                            points(x@data$longitude, x@data$latitude, cex = cex, pch = pch, col = col)
                        }
                    }
                }
                par(mar = mar)
            } else if (which[w] == 2) {
                # salinity timeseries
                if (0 != sum(!is.na(x@data$salinity))) {
                    nlevels <- dim(x@data$salinity)[1]
                    t <- if (length(level) > 1) {
                        numberAsPOSIXct(t(matrix(rep(x@data$time, nlevels), byrow = FALSE, ncol = nlevels)))
                    } else {
                        x@data$time
                    }
                    oce.plot.ts(t, as.vector(x@data$salinity[level, ]),
                        ylab = resizableLabel("S", "y"),
                        cex = cex, pch = pch, col = col, type = type, tformat = tformat
                    )
                } else {
                    warning("no non-missing salinity data")
                }
            } else if (which[w] == 3) {
                # temperature timeseries
                if (0 != sum(!is.na(x@data$temperature))) {
                    nlevels <- dim(x@data$temperature)[1]
                    t <- if (length(level) > 1) {
                        numberAsPOSIXct(t(matrix(rep(x@data$time, nlevels), byrow = FALSE, ncol = nlevels)))
                    } else {
                        x@data$time
                    }
                    oce.plot.ts(t, x@data$temperature[level, ],
                        ylab = resizableLabel("T", "y"),
                        cex = cex, pch = pch, col = col, type = type, tformat = tformat
                    )
                } else {
                    warning("no non-missing temperature data")
                }
            } else if (which[w] == 4) {
                # TS
                if (0 != sum(!is.na(x@data$temperature)) && 0 != sum(!is.na(x@data$salinity))) {
                    plotTS(ctd, cex = cex, pch = pch, col = col, type = type, debug = debug - 1)
                } else {
                    warning("no non-missing salinity data")
                }
            } else if (which[w] == 5) { # "salinity profile"
                plotProfile(ctd,
                    xtype = "salinity",
                    ylim = rev(range(x@data$pressure, na.rm=TRUE)),
                    cex = cex, pch = pch, col = col, type = type
                )
            } else if (which[w] == 6) { # "temperature profile"
                plotProfile(ctd,
                    xtype = "temperature",
                    ylim = rev(range(x@data$pressure, na.rm=TRUE)),
                    cex = cex, pch = pch, col = col, type = type
                )
            } else if (which[w] == 7) { "sigma0 profile"
                # sigma profile
                sigma0lim <- quantile(ctd[["sigma0"]], c(0.01, 0.99), na.rm = TRUE)
                plotProfile(ctd,
                    xtype = "sigma0",
                    ylim = rev(range(x@data$pressure, na.rm=TRUE)),
                    cex = cex, pch = pch, col = col, type = type
                )
            } else if (which[w] == 8) { "spice profile"
                plotProfile(ctd,
                    xtype = "spice",
                    ylim = rev(range(x@data$pressure, na.rm=TRUE)),
                    xlab = resizableLabel("spice"),
                    cex = cex, pch = pch, col = col, type = type
                )
            } else if (which[w] == 9) { "spiciness0 profile"
                plotProfile(ctd,
                    xtype = ctd[["spiciness0"]],
                    ylim = rev(range(x@data$pressure, na.rm=TRUE)),
                    xlab = resizableLabel("spiciness0"),
                    cex = cex, pch = pch, col = col, type = type
                )
            } else if (which[w] == 10) { "spiciness1 profile"
                plotProfile(ctd,
                    xtype = ctd[["spiciness1"]],
                    ylim = rev(range(x@data$pressure, na.rm=TRUE)),
                    xlab = resizableLabel("spiciness1"),
                    cex = cex, pch = pch, col = col, type = type
                )
            } else if (which[w] == 11) { "spiciness2 profile"
                plotProfile(ctd,
                    xtype = ctd[["spiciness2"]],
                    ylim = rev(range(x@data$pressure, na.rm=TRUE)),
                    xlab = resizableLabel("spiciness2"),
                    cex = cex, pch = pch, col = col, type = type
                )
            } else {
                stop("Unknown value of which=", which[w], "\n", call. = FALSE)
            }
        }
        oceDebug(debug, "} # plot.argo()\n", unindent = 1, style = "bold")
        invisible(NULL)
    }
)

# DEVELOPERS: please pattern functions and documentation on the 'ctd' code, for uniformity.
# DEVELOPERS: You will need to change the docs, and the 3 spots in the code
# DEVELOPERS: marked '# DEVELOPER 1:', etc.
#' @title Handle Flags in argo Objects
#'
#' @param object an [argo-class] object.
#'
#' @template handleFlagsTemplate
#'
#' @references
#' 1. Wong, Annie, Robert Keeley, Thierry Carval, and Argo Data Management Team.
#' "Argo Quality Control Manual for CTD and Trajectory Data," January 1, 2020.
#' `https://archimer.ifremer.fr/doc/00228/33951/`.
#'
#' @examples
#' library(oce)
#' data(argo)
#' argoNew <- handleFlags(argo)
#' # Demonstrate replacement, looking at the second profile
#' f <- argo[["salinityFlag"]][, 2]
#' df <- data.frame(flag = f, orig = argo[["salinity"]][, 2], new = argoNew[["salinity"]][, 2])
#' df[11:15, ] # notice line 13
#'
#' @author Dan Kelley
#'
#' @family things related to argo data
#' @aliases handleFlags.argo
setMethod("handleFlags",
    signature = c(object = "argo", flags = "ANY", actions = "ANY", where = "ANY", debug = "ANY"),
    definition = function(object, flags = NULL, actions = NULL, where = NULL, debug = getOption("oceDebug")) {
        # DEVELOPER 1: alter the next comment to explain your setup
        if (is.null(flags)) {
            flags <- defaultFlags(object)
            if (is.null(flags)) {
                stop("must supply 'flags', or use initializeFlagScheme() on the argo object first")
            }
        }
        if (is.null(actions)) {
            actions <- list("NA") # DEVELOPER 3: alter this line to suit a new data class
            names(actions) <- names(flags)
        }
        if (any(names(actions) != names(flags))) {
            stop("names of flags and actions must match")
        }
        handleFlagsInternal(object = object, flags = flags, actions = actions, where = where, debug = debug)
    }
)


#' Set Preference for Adjusted Values
#'
#' [argo-class] data can contain "adjusted" forms of data items,
#' which may be more trustworthy than the original
#' data, and `preferAdjusted` lets the user express a preference
#' for such adjusted data.  This means that using
#' \code{\link{[[,argo-method}} on the results returned by `preferAdjusted`
#' will (if possible) return adjusted data, and also use those adjusted
#' data in computations of derived quantities such as Absolute Salinity.
#' The preference applies also to units and to data-quality flags,
#' both of which can be returned by \code{\link{[[,argo-method}}, as
#' discussed in \dQuote{Details}.
#'
#' `preferAdjusted()` merely sets two items in the `metadata` slot of the
#' returned [argo-class] object. The real action is carried out by
#' \code{\link{[[,argo-method}} but, for convenience, the details are explained here.
#'
#' Consider salinity, for example.
#' If `which` equals `"all"`, or if it is a character
#' vector containing `"salinity"`, then using
#' \code{\link{[[,argo-method}} on the returned object
#' will yield the adjusted forms of the salinity data,
#' its associated flags, or its units.  Thus, in the salinity
#' case,
#' * `argo[["salinity"]]` will attempt to return `argo@data$salinityAdjusted`
#' instead of returning `argo@data$salinity`, although if the adjusted values
#' are all `NA` then, depending on the value of `fallback`, the
#' unadjusted values may be returned; similarly
#' * `argo[["salinityFlags"]]` will attempt to return
#' `argo@metadata$flags$salinityAdjusted`
#' instead of `argo@metadata$flags$salinity`, and
#' * `argo[["salinityUnits"]]` will attempt to return
#' `argo@metadata$units$salinityAdjusted`
#' instead of `argo@metadata$units$salinity`.
#'
#' The default value, `which="all"`, indicates that this
#' preference for adjusted values will apply to all the
#' elements of the `data` slot of the returned object, along
#' with associated flags and units. This can be handy for quick
#' work, but analysts may also choose to restrict their use of
#' adjusted values to a subset of variables, based on their own
#' decisions about data quality or accuracy.
#'
#' The default value `fallback=TRUE` indicates that later calls to
#' \code{\link{[[,argo-method}} should return unadjusted values for any
#' data items that have `NA` for all the adjusted values.  This
#' condition is rare for core variables (salinity, temperature and
#' pressure) but is annoyingly common for biogeochemical variables; see
#' e.g. Section 2.2.5 of Reference 1 for a discussion of
#' the conditions under which Argo netcdf files contain
#' adjusted values. Setting `fallback=FALSE` means that adjusted
#' values (if they exist) will always be returned, even if they
#' are a useless collection of `NA` values.
#'
#' Error fields, such as `salinityAdjustedError`, are returned
#' as-is by \code{\link{[[,argo-method}}, regardless of whether
#' the object was created by `preferAdjusted`.
#'
#' It should be noted that, regardless of whether `preferAdjusted`
#' has been used, the analyst can always access either unadjusted
#' or adjusted data directly, using the original variable names stored
#' in the source netcdf file.  For example, `argo[["PSAL"]]`
#' yields unadjusted salinity values, and
#' `argo[["PSAL_ADJUSTED"]]` yields adjusted values (if they exist, or
#' `NULL` if they do not).
#' Similarly, adjusted value can always be obtained by using a form
#' like `argo[["salinityAdjusted"]]`.
#'
#' @param argo An [argo-class] object.
#'
#' @param which A character vector naming the items for which
#' (depending also on the value of `fallback`) adjusted values
#' are to be sought by future calls to \code{\link{[[,argo-method}}.
#' The short names are used, e.g. `which="oxygen"` means that
#' adjusted oxygen is to be returned in future calls
#' such as `argo[["oxygen"]]`.  The default,
#' `"all"`, means to  use adjusted values for any item in `argo`
#' that has adjusted values.
#'
#' @param fallback A logical value indicating whether to fall back
#' to unadjusted values for any data field in which the
#' adjusted values are all `NA`.  The default value, `TRUE`,
#' avoids a problem with biogeochemical fields, where adjustment
#' of any one field may lead to insertion of "adjusted" values for
#' other fields that consist of nothing more than `NA`s.
#'
#' @return An [argo-class] object its `metadata` slot altered
#' (in its `adjustedWhich` and `adjustedFallback` elements)
#' as a signal for how \code{\link{[[,argo-method}} should
#' function on the object.
#'
#' @examples
#' library(oce)
#' data(argo)
#' argoAdjusted <- preferAdjusted(argo)
#' all.equal(argo[["salinityAdjusted"]], argoAdjusted[["salinity"]])
#' all.equal(argo[["salinityFlagsAdjusted"]], argoAdjusted[["salinityFlags"]])
#' all.equal(argo[["salinityUnitsAdjusted"]], argoAdjusted[["salinityUnits"]])
#'
#' @references
#' 1. Argo Data Management Team. "Argo User's Manual V3.3." Ifremer,
#' November 28, 2019.
#' \doi{10.13155/29825}
#'
#' @author Dan Kelley, based on discussions with Jaimie Harbin (with
#' respect to the \code{\link{[[,argo-method}} interface) and Clark Richards
#' (with respect to storing the preference in the `metadata` slot).
preferAdjusted <- function(argo, which = "all", fallback = TRUE) {
    if (!inherits(argo, "argo")) {
        stop("'argo' must be an oce 'argo' object")
    }
    if (!is.logical(fallback)) {
        stop("fallback must be a logical value")
    }
    argo@metadata$adjustedFallback <- fallback
    argo@metadata$adjustedWhich <- which
    argo
}

#' Convert Time to Argo Julian Day (juld)
#'
#' This converts a POSIXct time into an Argo julian day, with the convention
#' that juld=0 at 1950-01-01.
#'
#' @param t A POSIXct time or a string that can be converted to a POSIXct time
#'
#' @examples
#' timeToArgoJuld("2020-07-01")
#'
#' @author Jaimie Harbin and Dan Kelley
timeToArgoJuld <- function(t) {
    oce::julianDay(as.POSIXct(t, tz = "UTC")) - oce::julianDay(as.POSIXct("1950-01-01", tz = "UTC"))
}

#' Convert Argo Julian Day to R Time
#'
#' @param jday A numerical value indicating the julian day in the Argo convention,
#' with day=0 at 1950-01-01.
#'
#' @examples
#' argoJuldToTime(25749)
#'
#' @author Jaimie Harbin and Dan Kelley
argoJuldToTime <- function(jday) {
    as.POSIXct("1950-01-01", tz = "UTC") + jday * 86400
}