R/argo.R

In 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 = "", 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
        # message("argo[[..., i=\"", i, "\"]]\n")
        if (i %in% c(
            "CT", paste("Conservative", "Temperature"), "N2",
            "SA", paste("Absolute", "Salinity"), "sigmaTheta",
            "theta", "spice"
        )) {
            # message("about to get S, T, p")
            salinity <- x[["salinity", debug = debug - 1]]
            temperature <- x[["temperature", debug = debug - 1]]
            pressure <- x[["pressure", debug = debug - 1]]
            Sshape <- shape(salinity)
            # Actually, do not need longitude and latitude if eos="unesco"
            if (is.matrix(pressure)) {
                Nrow <- nrow(pressure)
                longitude <- rep(x@data$longitude, each = Nrow)
                latitude <- rep(x@data$latitude, each = Nrow)
                dim(longitude) <- Sshape
                dim(latitude) <- Sshape
            } else {
                Nrow <- length(pressure)
                longitude <- rep(x@data$longitude, each = Nrow)
                latitude <- rep(x@data$latitude, each = Nrow)
            }
            if (!identical(Sshape, shape(temperature))) stop("salinity and temperature have different lengths or dimensions")
            if (!identical(Sshape, shape(pressure))) stop("salinity and pressure have different lengths or dimensions")
            if (!identical(Sshape, shape(longitude))) stop("salinity and longitude have differing lengths or dimensions")
            if (!identical(Sshape, shape(latitude))) stop("salinity and latitude have differing lengths or dimensions")
            if (i %in% c("CT", "Conservative Temperature")) {
                res <- gsw::gsw_CT_from_t(x[["SA"]], temperature, pressure)
            } else if (i == "N2") {
                # nprofile <- dim[2]
                oceDebug(debug, "argo[[\"N2\"]] ...\n")
                res <- array(NA_real_, dim = Sshape)
                for (i in seq_len(Sshape[2])) {
                    # message("i=",i, ", nprofile=", nprofile)
                    # if (i == 14) browser()
                    if (sum(!is.na(pressure[, i])) > 2) {
                        oceDebug(debug - 1, "  OK i=", i, "\n")
                        ctd <- as.ctd(
                            salinity = salinity[, i],
                            temperature = temperature[, i],
                            pressure = pressure[, i],
                            longitude = longitude[1, i],
                            latitude = latitude[1, i]
                        )
                        res[, i] <- swN2(ctd)
                    } else {
                        oceDebug(debug - 1, "  BAD i=", i, "\n")
                        res[, i] <- rep(NA, length(salinity[, i]))
                    }
                }
            } else if (i %in% c("SA", "Absolute Salinity")) {
                # message("DAN 'SA' or 'Absolute Salinity'")
                # cat(vectorShow(salinity))
                # cat(vectorShow(temperature))
                # cat(vectorShow(pressure))
                # cat(vectorShow(longitude))
                # cat(vectorShow(latitude))
                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 == "spice") {
                if (missing(j)) {
                    res <- swSpice(x)
                } else {
                    # message("DAN 1 (spice about to check j)")
                    if (j != "gsw" && j != "unesco") {
                        stop("\"", j, "\" not allowed; use either \"gsw\" or \"unesco\"")
                    }
                    # message("DAN 2")
                    res <- swSpice(x, eos = j)
                }
            } else if (i == "sigmaTheta") {
                oceDebug(debug, "  argo[[\"sigmaTheta\"]] is about to call swSigmaTheta()\n")
                res <- swSigmaTheta(salinity,
                    temperature = temperature, pressure = pressure,
                    referencePressure = 0, longitude = longitude, latitude = latitude,
                    eos = getOption("oceEOS", default = "gsw"), debug = debug
                )
            } else if (i == "theta") {
                # message("argo.R: 281")
                res <- swTheta(salinity,
                    temperature = temperature, pressure = pressure,
                    referencePressure = 0, longitude = longitude, latitude = latitude,
                    eos = getOption("oceEOS", default = "gsw"),
                    debug = debug
                )
                # browser()
            } else {
                stop("argo[[ coding error: unknown item '", i, "'")
            }
            if (length(Sshape > 1)) {
                dim(res) <- Sshape
            }
        } else if (i == "z") {
            res <- swZ(x)
        } else if (i == "depth") {
            res <- swDepth(x)
        } 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 {
            oceDebug(debug, "calling lowest-level oce '[[' method (in AllClass.R)\n")
            res <- callNextMethod() # [[ defined in R/AllClass.R
        }
        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() START\n", sep = "", unindent = 1)
        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, "END subset,argo-method\n", sep = "", unindent = 1)
        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() START\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 `read.argo` returns 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")) {
            ncdf4::nc_close(file)
            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, "\", ...) START\n", sep = "", unindent = 1)
    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)) {
                # for the iconv() call and the gsub() test, see
                # https://github.com/dankelley/oce/issues/2206
                origValue <- value
                # cat("next is origValue:\n");print(origValue)
                value <- iconv(value, from = "latin1", to = "UTF-8")
                value <- try(
                    {
                        trimws(value)
                        # gsub("^[ \t\r\n](.*)[ \t\r\n]$", "\\1", value, useBytes=TRUE)
                    },
                    silent = TRUE
                )
                # cat("after trimws(), next is value:\n");print(value)
                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 = ""))
    ncdf4::nc_close(file)
    oceDebug(debug, "END read.argo()\n", sep = "", unindent = 1)
    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 spiciness 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),
            " ...) START\n",
            sep = "", unindent = 1
        )
        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, "END plot.argo()\n", unindent = 1)
        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
}