R/echosounder.R

In dankelley/oce: Analysis of Oceanographic Data

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

# REFERENCES:
#   [1] "DT4 Data File Format Specification" [July, 2010] DT4_format_2010.pdf


#' Class to Store Echosounder Data
#'
#' This class stores echosounder data. Echosounder objects may be
#' read with [read.echosounder()],
#' summarized with [summary,echosounder-method()],
#' and plotted with [plot,echosounder-method()].
#' The [findBottom()]
#' function infers the ocean bottom from tracing the strongest reflector from
#' ping to ping.
#'
#' @templateVar class echosounder
#'
#' @templateVar dataExample {}
#'
#' @templateVar metadataExample {}
#'
#' @template slot_summary
#'
#' @template slot_put
#'
#' @template slot_get
#'
#' @details
#'
#' * An infrequently updated record of the instrument position, in
#' `timeSlow`, `longitudeSlow` and `latitudeSlow`.  These are
#' used in plotting maps with [plot,echosounder-method()].
#'
#' * An interpolated record of the instrument position, in `time`,
#' `longitude`, and `latitude`.  Linear interpolation is used to
#' infer the longitude and latitude from the variables listed above.
#'
#' * `depth`, vector of depths of echo samples (measured positive
#' downwards in the water column).  This is calculated from the inter-sample
#' time interval and the sound speed provided as the `soundSpeed` argument
#' to [read.echosounder()], so altering the value of the latter will
#' alter the echosounder plots provided by [plot,echosounder-method()].
#'
#' * The echosounder signal amplitude `a`, a matrix whose number of
#' rows matches the length of `time`, etc., and number of columns equal to
#' the length of `depth`.  Thus, for example, `a[100,]` represents
#' the depth-dependent amplitude at the time of the 100th ping.
#'
#' * A matrix named `b` exists for dual-beam and split-beam cases.
#' For dual-beam data, this is the wide-beam data, whereas `a` is the
#' narrow-beam data.  For split-beam data, this is the x-angle data.
#'
#' * A matrix named `c` exists for split-beam data, containing the
#' y-angle data.
#'
#' * In addition to these matrices, ad-hoc calculated matrices named
#' `Sv` and `TS` may be accessed as explained in the next section.
#'
#' @name echosounder-class
#'
#' @docType class
#'
#' @author Dan Kelley
#'
#' @family things related to echosounder data
setClass("echosounder", contains = "oce")


#' Sample echosounder Data
#'
#' This is degraded subsample of measurements that were made with a Biosonics
#' scientific echosounder, as part of the St Lawrence Internal Wave Experiment
#' (SLEIWEX).
#'
#' @name echosounder
#'
#' @docType data
#'
#' @author Dan Kelley
#'
#' @source This file came from the SLEIWEX-2008 experiment, and was decimated
#' using [decimate()] with `by=c()`.
#' @family datasets provided with oce
#' @family things related to echosounder data
NULL

setMethod(
    f = "initialize",
    signature = "echosounder",
    definition = function(.Object, filename = "", ...) {
        .Object <- callNextMethod(.Object, ...)
        .Object@metadata$filename <- filename
        .Object@processingLog$time <- presentTime()
        .Object@processingLog$value <- "create 'echosounder' object"
        return(.Object)
    }
)


#' Summarize an echosounder Object
#'
#' Summarizes some of the data in an [echosounder-class] object.
#'
#' @param object an object of class `"echosounder"`, usually, a result of
#' a call to [read.echosounder()], [read.oce()], or
#' [as.echosounder()].
#'
#' @param \dots further arguments passed to or from other methods.
#'
#' @author Dan Kelley
#'
#' @family things related to echosounder data
setMethod(
    f = "summary",
    signature = "echosounder",
    definition = function(object, ...) {
        cat("Echosounder Summary\n-------------------\n\n")
        showMetadataItem(object, "filename", "File source:         ", quote = TRUE)
        showMetadataItem(object, "transducerSerialNumber", "Transducer serial #: ", quote = FALSE)
        metadataNames <- names(object@metadata)
        cat(sprintf("* File type:           %s\n", object[["fileType"]]))
        if ("beamType" %in% metadataNames) {
            cat(sprintf("* Beam type:           %s\n", object[["beamType"]]))
        }
        time <- object[["time"]]
        tz <- attr(time[1], "tzone")
        nt <- length(time)
        cat(sprintf("* Channel:             %d\n", object[["channel"]]))
        cat(sprintf("* Measurements:        %s %s to %s %s\n", format(time[1]), tz, format(time[nt]), tz))
        cat(sprintf("* Sound speed:         %.2f m/s\n", object[["soundSpeed"]]))
        # cat(sprintf("* Time between pings:  %.2e s\n", object[["samplingDeltat"]]))
        if ("pulseDuration" %in% metadataNames) {
            cat(sprintf("* Pulse duration:      %g s\n", object[["pulseDuration"]] / 1e6))
        }
        cat(sprintf("* Frequency:           %f\n", object[["frequency"]]))
        cat(sprintf("* Blanked samples:     %d\n", object[["blankedSamples"]]))
        cat(sprintf("* Pings in file:       %d\n", object[["pingsInFile"]]))
        cat(sprintf("* Samples per ping:    %d\n", object[["samplesPerPing"]]))
        invisible(callNextMethod()) # summary
    }
)


#' Extract Something From an echosounder Object
#'
#' @param x an [echosounder-class] object.
#'
#' @templateVar class echosounder
#'
#' @section Details of the Specialized Method:
#'
#' * If `i` is `"?"`, then the return value is a list
#' containing four items, each of which is a character vector
#' holding the names of things that can be accessed with `[[`.
#' The `data` and `metadata` items hold the names of
#' entries in the object's data and metadata
#' slots, respectively. The `metadataDerived` item
#' is NULL, while the `dataDerived` item holds
#' `"Sv"` and `"TS"` (see next).
#'
#' * If `i` is `"Sv"`, then the following is returned.
#' \preformatted{
#' 20*log10(a) -
#'   (x@@metadata$sourceLevel+x@@metadata$receiverSensitivity+x@@metadata$transmitPower) +
#'   20*log10(r) +
#'   2*absorption*r -
#'   x@@metadata$correction +
#'   10*log10(soundSpeed*x@@metadata$pulseDuration/1e6*psi/2)
#' }
#'
#' * If `i` is `"TS"`, then the following is returned.
#' \preformatted{
#' 20*log10(a) -
#'   (x@@metadata$sourceLevel+x@@metadata$receiverSensitivity+x@@metadata$transmitPower) +
#'   40*log10(r) +
#'   2*absorption*r +
#'   x@@metadata$correction
#' }
#'
#' @template sub_subTemplate
#'
#' @author Dan Kelley
#'
#' @family things related to echosounder data
setMethod(
    f = "[[",
    signature(x = "echosounder", i = "ANY", j = "ANY"),
    definition = function(x, i, j, ...) {
        if (i == "?") {
            return(list(
                metadata = sort(names(x@metadata)), metadataDerived = NULL,
                data = sort(names(x@data)), dataDerived = c("Sv", "TS")
            ))
        }
        if (i %in% c("Sv", "TS")) {
            range <- rev(x@data$depth)
            a <- x@data$a
            # biosonics has /20 because they have bwx in 0.1deg
            psi <- x@metadata$beamwidthX / 2 * x@metadata$beamwidthY / 2 / 10^3.16
            r <- matrix(rev(range), nrow = nrow(a), ncol = length(range), byrow = TRUE)
            absorption <- swSoundAbsorption(x@metadata$frequency, 35, 10, mean(range))
            soundSpeed <- x@metadata$soundSpeed
            if (i == "Sv") {
                Sv <- 20 * log10(a) -
                    (x@metadata$sourceLevel + x@metadata$receiverSensitivity + x@metadata$transmitPower) +
                    20 * log10(r) +
                    2 * absorption * r -
                    x@metadata$correction +
                    10 * log10(soundSpeed * x@metadata$pulseDuration / 1e6 * psi / 2)
                Sv[!is.finite(Sv)] <- NA
                Sv
            } else if (i == "TS") {
                TS <- 20 * log10(a) -
                    (x@metadata$sourceLevel + x@metadata$receiverSensitivity + x@metadata$transmitPower) +
                    40 * log10(r) +
                    2 * absorption * r +
                    x@metadata$correction
                TS[!is.finite(TS)] <- NA
                TS
            }
        } else {
            callNextMethod() # [[
        }
    }
)

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


#' Subset an echosounder Object
#'
#' This function is somewhat analogous to [subset.data.frame()].
#' Subsetting can be by `time` or `depth`, but these may not be
#' combined; use a sequence of calls to subset by both.
#'
#' @param x an [echosounder-class] object.
#'
#' @param subset a condition to be applied to the `data` portion of
#' `x`.  See \dQuote{Details}.
#'
#' @param \dots ignored.
#'
#' @return An [echosounder-class] object.
#'
#' @author Dan Kelley
#'
#' @examples
#' library(oce)
#' data(echosounder)
#' plot(echosounder)
#' plot(subset(echosounder, depth < 10))
#' plot(subset(echosounder, time < mean(range(echosounder[["time"]]))))
#'
#' @family things related to echosounder data
#' @family functions that subset oce objects
setMethod(
    f = "subset",
    signature = "echosounder",
    definition = function(x, subset, ...) {
        subsetString <- paste(deparse(substitute(expr = subset, env = environment())), collapse = " ")
        res <- x
        dots <- list(...)
        debug <- if (length(dots) && ("debug" %in% names(dots))) dots$debug else getOption("oceDebug")
        if (missing(subset)) {
            stop("must give 'subset'")
        }
        if (length(grep("time", subsetString))) {
            oceDebug(debug, "subsetting an echosounder object by time\n")
            keep <- eval(substitute(expr = subset, env = environment()), envir = x@data, enclos = parent.frame(2))
            oceDebug(debug, "keeping", 100 * sum(keep) / length(keep), "% of the fast-sampled data\n")
            res <- x
            # trim fast variables, handling matrix 'a' differently, and skipping 'distance'
            dataNames <- names(res@data)
            res@data$a <- x@data$a[keep, ]
            if ("b" %in% dataNames) {
                res@data$b <- x@data$b[keep, ]
            }
            if ("c" %in% dataNames) {
                res@data$c <- x@data$c[keep, ]
            }
            # lots of debugging in here, in case other data types have other variable names
            oceDebug(debug, "dataNames (orig):", dataNames, "\n")
            if (length(grep("^a$", dataNames))) {
                dataNames <- dataNames[-grep("^a$", dataNames)]
            }
            if (length(grep("^b$", dataNames))) {
                dataNames <- dataNames[-grep("^b$", dataNames)]
            }
            if (length(grep("^c$", dataNames))) {
                dataNames <- dataNames[-grep("^c$", dataNames)]
            }
            oceDebug(debug, "dataNames (step 2):", dataNames, "\n")
            if (length(grep("^depth$", dataNames))) {
                dataNames <- dataNames[-grep("^depth$", dataNames)]
            }
            oceDebug(debug, "dataNames (step 3):", dataNames, "\n")
            if (length(grep("Slow", dataNames))) {
                dataNames <- dataNames[-grep("Slow", dataNames)]
            }
            oceDebug(debug, "dataNames (final), i.e. fast dataNames to be trimmed by time:", dataNames, "\n")
            for (dataName in dataNames) {
                oceDebug(debug, "fast variable:", dataName, "orig length", length(x@data[[dataName]]), "\n")
                res@data[[dataName]] <- x@data[[dataName]][keep]
                oceDebug(debug, "fast variable:", dataName, "new length", length(res@data[[dataName]]), "\n")
            }
            # trim slow variables
            subsetStringSlow <- gsub("time", "timeSlow", subsetString)
            oceDebug(debug, "subsetting slow variables with string:", subsetStringSlow, "\n")
            keepSlow <- eval(parse(text = subsetStringSlow), x@data, parent.frame(2))
            oceDebug(debug, "keeping", 100 * sum(keepSlow) / length(keepSlow), "% of the slow-sampled data\n")
            for (slowName in names(x@data)[grep("Slow", names(x@data))]) {
                oceDebug(debug, "slow variable:", slowName, "orig length", length(x@data[[slowName]]), "\n")
                res@data[[slowName]] <- x@data[[slowName]][keepSlow]
                oceDebug(debug, "slow variable:", slowName, "new length", length(res@data[[slowName]]), "\n")
            }
        } else if (length(grep("depth", subsetString))) {
            oceDebug(debug, "subsetting an echosounder object by depth\n")
            keep <- eval(expr = substitute(expr = subset, env = environment()), envir = x@data, enclos = parent.frame(2))
            res <- x
            res[["depth"]] <- res[["depth"]][keep]
            dataNames <- names(res@data)
            res[["a"]] <- res[["a"]][, keep]
            if ("b" %in% dataNames) {
                res@data$b <- x@data$b[, keep]
            }
            if ("c" %in% dataNames) {
                res@data$c <- x@data$c[, keep]
            }
        } else {
            stop("can only subset an echosounder object by 'time' or 'depth'")
        }
        res@processingLog <- processingLogAppend(res@processingLog, paste("subset.echosounder(x, subset=", subsetString, ")", sep = ""))
        res
    }
) # subset,echosounder-method


#' Coerce Data Into an echosounder Object
#'
#' Coerces a dataset into a echosounder dataset.
#'
#' Creates an echosounder file.  The defaults for e.g.  `transmitPower`
#' are taken from the `echosounder` dataset, and they are unlikely to make
#' sense generally.  The first three parameters must be supplied, and the
#' dimension of `a` must align with the lengths of `time` and `depths`. The
#' other parameters have defaults that are unlikely to be correct for
#' arbitrary application, but they are not essential for basic plots,
#' etc.
#'
#' Those who use the \CRANpkg{readHAC} package to read echosounder
#' data should note that it stores the `a` matrix in a flipped and
#' transposed format. See that package's demo code for a function
#' named `flip()` that transforms the matrix as required by
#' [as.echosounder()]. Indeed, users working with HAC
#' data ought to study the whole of the \CRANpkg{readHAC}
#' documentation, to learn what data are stored, so that
#' [oceSetMetadata()] and [oceSetData()] can be used as needed
#' to flesh out the contents returned by [as.echosounder()].
#'
#' @param time times of pings.
#'
#' @param depth depths of samples within pings.
#'
#' @param a matrix of echo amplitudes, as would be stored with
#' [read.echosounder()].
#'
#' @param src optional string indicating data source.
#'
#' @param sourceLevel source level, in dB (uPa at 1m), denoted `sl` in
#' reference 1 p15, where it is in units 0.1dB (uPa at 1m).
#'
#' @param receiverSensitivity receiver sensitivity of the main element, in
#' dB(counts/uPa), denoted `rs` in reference 1 p15, where it is in units of
#' 0.1dB(counts/uPa)
#'
#' @param transmitPower transmit power reduction factor, in dB, denoted
#' `tpow` in reference 1 p10, where it is in units 0.1 dB.
#'
#' @param pulseDuration duration of transmitted pulse in us
#'
#' @param beamwidthX x-axis -3dB one-way beamwidth in deg, denoted `bwx`
#' in reference 1 p16, where the unit is 0.2 deg
#'
#' @param beamwidthY y-axis -3dB one-way beamwidth in deg, denoted `bwx`
#' in reference 1 p16, where the unit is 0.2 deg
#'
#' @param frequency transducer frequency in Hz, denoted `fq` in reference 1 p16
#'
#' @param correction user-defined calibration correction in dB, denoted
#' `corr` in reference 1 p14, where the unit is 0.01dB.
#'
#' @return An [echosounder-class] object.
#'
#' @author Dan Kelley
#'
#' @family things related to echosounder data
as.echosounder <- function(
    time, depth, a, src = "",
    sourceLevel = 220, receiverSensitivity = -55.4, transmitPower = 0, pulseDuration = 400,
    beamwidthX = 6.5, beamwidthY = 6.5, frequency = 41800, correction = 0) {
    res <- new("echosounder", filename = src)
    res@metadata$channel <- 1
    res@metadata$soundSpeed <- swSoundSpeed(35, 10, 1, eos = "unesco") # don't bother with gsw
    res@metadata$samplingDeltat <- as.numeric(time[2]) - as.numeric(time[1])
    dim <- dim(a)
    res@metadata$pingsInFile <- dim[1]
    res@metadata$samplesPerPing <- dim[2]
    # args
    res@metadata$sourceLevel <- sourceLevel
    res@metadata$receiverSensitivity <- receiverSensitivity
    res@metadata$transmitPower <- transmitPower
    res@metadata$pulseDuration <- pulseDuration
    res@metadata$beamwidthX <- beamwidthX
    res@metadata$beamwidthY <- beamwidthY
    res@metadata$frequency <- frequency
    res@metadata$correction <- correction
    # FIXME: what about timeLocation, latitude, and longitude?
    res@data$time <- time
    res@data$depth <- depth
    res@data$a <- a
    names <- names(res@data)
    if ("latitude" %in% names) res@metadata$units$latitude <- list(unit = expression(degree * N), scale = "")
    if ("longitude" %in% names) res@metadata$units$longitude <- list(unit = expression(degree * E), scale = "")
    if ("latitudeSlow" %in% names) res@metadata$units$latitudeSlow <- list(unit = expression(degree * N), scale = "")
    if ("longitudeSlow" %in% names) res@metadata$units$longitudeSlow <- list(unit = expression(degree * E), scale = "")
    if ("depth" %in% names) {
        res@metadata$units$depth <- list(unit = expression(m), scale = "")
    }
    res@processingLog <- processingLogAppend(res@processingLog, paste(deparse(match.call()), sep = "", collapse = ""))
    res
}



#' Find the Ocean Bottom in an Echosounder Object
#'
#' Finds the depth in a Biosonics echosounder file, by finding the strongest
#' reflector and smoothing its trace.
#'
#' @param x an [echosounder-class] object.
#'
#' @param ignore number of metres of data to ignore, near the surface.
#'
#' @param clean a function to clean the inferred depth of spikes.
#'
#' @return A list with elements: the `time` of a ping, the `depth` of
#' the inferred depth in metres, and the `index` of the inferred bottom
#' location, referenced to the object's `depth` vector.
#'
#' @author Dan Kelley
#'
#' @seealso The documentation for [echosounder-class] explains the
#' structure of `echosounder` objects, and also outlines the other
#' functions dealing with them.
#'
#' @family things related to echosounder data
findBottom <- function(x, ignore = 5, clean = despike) {
    a <- x[["a"]]
    keep <- x[["depth"]] >= ignore
    wm <- clean(apply(a[, keep], 1, which.max))
    depth <- x[["depth"]][wm]
    list(time = x[["time"]], depth = depth, index = wm)
}


#' Plot an echosounder Object
#'
#' Plot echosounder data.
#' Simple linear approximation is used when a `newx` value is specified
#' with the `which=2` method, but arguably a gridding method should be
#' used, and this may be added in the future.
#'
#' @param x an [echosounder-class] object.
#'
#' @param which list of desired plot types: `which=1` or \code{which="zt
#' image"} gives a z-time image, `which=2` or `which="zx image"`
#' gives a z-distance image, and `which=3` or `which="map"` gives a
#' map showing the cruise track.  In the image plots, the display is of
#' [log10()] of amplitude, trimmed to zero for any amplitude values
#' less than 1 (including missing values, which equal 0).  Add 10 to the
#' numeric codes to get the secondary data (non-existent for single-beam files,
#'
#' @param beam a more detailed specification of the data to be plotted.  For
#' single-beam data, this may only be `"a"`.  For dual-beam data, this may
#' be `"a"` for the narrow-beam signal, or `"b"` for the wide-beam
#' signal.  For split-beam data, this may be `"a"` for amplitude,
#' `"b"` for x-angle data, or `"c"` for y-angle data.
#'
#' @param newx optional vector of values to appear on the horizontal axis if
#' `which=1`, instead of time.  This must be of the same length as the
#' time vector, because the image is remapped from time to `newx` using
#' [approx()].
#'
#' @param xlab,ylab optional labels for the horizontal and vertical axes; if
#' not provided, the labels depend on the value of `which`.
#'
#' @param xlim optional range for x axis.
#'
#' @param ylim optional range for y axis.
#'
#' @param zlim optional range for color scale.
#'
#' @param type type of graph, `"l"` for line, `"p"` for points, or
#' `"b"` for both.
#'
#' @param col a function providing the color scale for image plots.
#' This value is passed to [imagep()], which draws
#' the images.  Since [imagep()] defaults `col` to
#' [oceColorsViridis()], that is effectively also the default
#' for the present function. (Prior to 2023-03-18, the present
#' function defaulted `col` to [oceColorsJet()].)
#'
#' @param lwd line width (ignored if `type="p"`).
#'
#' @param atTop optional vector of time values, for labels at the top of the
#' plot produced with `which=2`.  If `labelsTop` is provided, then it
#' will hold the labels.  If `labelsTop` is not provided, the labels will
#' be constructed with the [format()] function, and these may be
#' customized by supplying a `format` in the `...` arguments.
#'
#' @param labelsTop optional vector of character strings to be plotted above
#' the `atTop` times.  Ignored unless `atTop` was provided.
#'
#' @param tformat optional argument passed to [imagep()], for plot
#' types that call that function.  (See [strptime()] for the format
#' used.)
#'
#' @param despike remove vertical banding by using [smooth()] to
#' smooth across image columns, row by row.
#'
#' @param drawBottom optional flag used for section images.  If `TRUE`,
#' then the bottom is inferred as a smoothed version of the ridge of highest
#' image value, and data below that are grayed out after the image is drawn.
#' If `drawBottom` is a color, then that color is used, instead of
#' white.  The bottom is detected with [findBottom()], using the
#' `ignore` value described next.
#'
#' @param ignore optional flag specifying the thickness in metres of a surface
#' region to be ignored during the bottom-detection process.  This is ignored
#' unless `drawBottom=TRUE`.
#'
#' @param drawTimeRange if `TRUE`, the time range will be drawn at the
#' top.  Ignored except for `which=2`, i.e. distance-depth plots.
#'
#' @param drawPalette if `TRUE`, the palette will be drawn.
#'
#' @param radius radius to use for maps; ignored unless `which=3` or
#' `which="map"`.
#'
#' @param coastline coastline to use for maps; ignored unless `which=3` or
#' `which="map"`.
#'
#' @param mgp 3-element numerical vector to use for [`par`]`("mgp")`, and also
#' for [`par`]`("mar")`, computed from this.  The default is tighter than the R
#' default, in order to use more space for the data and less for the axes.
#'
#' @param mar value to be used with [`par`]`("mar")`.
#'
#' @param debug set to an integer exceeding zero, to get debugging information
#' during processing.
#'
#' @param \dots optional arguments passed to plotting functions.  For example,
#' for maps, it is possible to specify the radius of the view in kilometres,
#' with `radius`.
#'
#' @return A list is silently returned, containing `xat` and `yat`,
#' values that can be used by [oce.grid()] to add a grid to the plot.
#'
#' @author Dan Kelley, with extensive help from Clark Richards
#'
#' @examples
#' library(oce)
#' data(echosounder)
#' plot(echosounder, drawBottom = TRUE)
#'
#' @family things related to echosounder data
#' @aliases plot.echosounder
setMethod(
    f = "plot",
    signature = signature("echosounder"),
    definition = function(x, which = 1, # 1=z-t section 2=dist-t section 3=map
                          beam = "a",
                          newx,
                          xlab, ylab,
                          xlim, ylim, zlim,
                          type = "l", col, lwd = 2,
                          despike = FALSE,
                          drawBottom, ignore = 5,
                          drawTimeRange = FALSE, drawPalette = TRUE,
                          radius, coastline,
                          mgp = getOption("oceMgp"),
                          mar = c(mgp[1], mgp[1] + 1.5, mgp[2] + 1 / 2, 1 / 2),
                          atTop, labelsTop,
                          tformat,
                          debug = getOption("oceDebug"),
                          ...) {
        dots <- list(...)
        res <- list(xat = NULL, yat = NULL)
        dotsNames <- names(dots)
        oceDebug(debug, "plot() { # for echosounder\n", unindent = 1, style = "bold")
        lw <- length(which)
        if (length(beam) < lw) {
            beam <- rep(beam, lw)
        }
        # opar <- par(no.readonly = TRUE)
        par(mgp = mgp, mar = mar)
        if (lw > 1) {
            # on.exit(par(opar))
            if (lw > 2) {
                layout(matrix(1:4, nrow = 2, byrow = TRUE))
            } else {
                layout(matrix(1:2, nrow = 2, byrow = TRUE))
            }
        }
        oceDebug(debug, "which before converting to numbers: c(", paste(which, collapse = ","), ")\n")
        whichOrig <- which
        which <- oce.pmatch(which, list("zt image" = 1, "zx image" = 2, map = 3))
        oceDebug(debug, "which after converting to numbers: c(", paste(which, collapse = ","), ")\n")
        for (w in seq_along(which)) {
            oceDebug(debug, "this which:", which[w], "\n")
            if (which[w] == 1) {
                oceDebug(debug, "handling which[", w, "]==1\n", sep = "")
                time <- x[["time"]]
                xInImage <- time
                if (!length(time)) {
                    stop("plot.echosounder() cannot plot, because @data$time has zero length")
                }
                signal <- x[[beam[w]]]
                newxGiven <- !missing(newx)
                if (newxGiven) {
                    t <- as.numeric(time)
                    if (length(newx) != length(t)) {
                        stop("length of 'newx' must match that of time within the object")
                    }
                    xInImage <- newx
                }
                if (despike) {
                    signal <- apply(signal, 2, smooth)
                }
                if (beam[w] == "Sv" || beam[w] == "TS") {
                    oceDebug(debug, "using signal for z")
                    z <- signal
                } else {
                    oceDebug(debug, "using log10(signal) for z\n")
                    z <- log10(signal)
                }
                z[!is.finite(z)] <- NA # prevent problem in computing range
                if (!missing(drawBottom)) {
                    oceDebug(debug, "drawBottom is TRUE\n")
                    if (is.logical(drawBottom) && drawBottom) {
                        drawBottom <- "lightgray"
                    }
                    waterDepth <- findBottom(x, ignore = ignore)$depth
                    axisBottom <- par("usr")[3]
                    deepestWater <- max(abs(waterDepth))
                    ats <- imagep(xInImage,
                        y = -x[["depth"]], z = z,
                        xlab = if (missing(xlab)) "" else xlab, # time
                        ylab = if (missing(ylab)) "z [m]" else ylab, # depth
                        xlim = xlim,
                        ylim = if (missing(ylim)) c(-deepestWater, 0) else ylim,
                        zlim = if (missing(zlim)) c(if (beam[w] %in% c("Sv", "TS")) min(z, na.rm = TRUE) else 0, max(z, na.rm = TRUE)) else zlim,
                        col = col,
                        mgp = mgp, mar = mar,
                        tformat = tformat,
                        drawPalette = drawPalette,
                        debug = debug - 1,
                        zlab = beam[w],
                        ...
                    )
                    axisBottom <- par("usr")[3]
                    waterDepth <- c(axisBottom, -waterDepth, axisBottom)
                    time <- x[["time"]]
                    if (newxGiven) {
                        newx2 <- approx(as.numeric(time), newx, as.numeric(time))$y
                        newx2 <- c(newx2[1], newx2, newx2[length(newx2)])
                        polygon(newx2, waterDepth, col = drawBottom)
                    } else {
                        time2 <- c(time[1], time, time[length(time)])
                        polygon(time2, waterDepth, col = drawBottom)
                    }
                } else {
                    oceDebug(debug, "drawBottom is FALSE\n")
                    ats <- imagep(xInImage,
                        y = -x[["depth"]], z = z,
                        xlab = if (missing(xlab)) "" else xlab, # time
                        ylab = if (missing(ylab)) "z [m]" else ylab, # depth
                        xlim = xlim,
                        ylim = if (missing(ylim)) c(-max(abs(x[["depth"]])), 0) else ylim,
                        zlim = if (missing(zlim)) c(if (beam[w] %in% c("Sv", "TS")) min(z, na.rm = TRUE) else 0, max(z, na.rm = TRUE)) else zlim,
                        col = col,
                        mgp = mgp, mar = mar,
                        tformat = tformat,
                        drawPalette = drawPalette,
                        debug = debug - 1,
                        zlab = beam[w],
                        ...
                    )
                }
                res$xat <- ats$xat
                res$yat <- ats$yat
                if (newxGiven) {
                    if (!missing(atTop)) {
                        at <- approx(as.numeric(x[["time"]]), newx, as.numeric(atTop))$y
                        if (missing(labelsTop)) {
                            labelsTop <- format(atTop, format = if ("format" %in% dotsNames) dots$format else "%H:%M:%S")
                        }
                        axis(side = 3, at = at, labels = labelsTop, cex.axis = par("cex"))
                    } else {
                        pretty <- pretty(time)
                        labels <- format(pretty, format = "%H:%M:%S")
                        at <- approx(as.numeric(time), newx, as.numeric(pretty))$y
                        axis(3, at = at, labels = labels, cex.axis = par("cex"))
                    }
                }
            } else if (which[w] == 2) {
                oceDebug(debug, "handling which[", w, "]==2\n", sep = "")
                latitude <- x[["latitude"]]
                longitude <- x[["longitude"]]
                jitter <- rnorm(length(latitude), sd = 1e-8) # jitter lat by equiv 1mm to avoid colocation
                distance <- geodDist(longitude, jitter + latitude, alongPath = TRUE) # FIXME: jitter should be in imagep
                depth <- x[["depth"]]
                a <- x[["a"]]
                if (despike) {
                    a <- apply(a, 2, smooth)
                }
                z <- log10(ifelse(a > 1, a, 1)) # FIXME: make an argument for this 1
                deepestWater <- max(abs(depth))
                if (!missing(drawBottom)) {
                    if (is.logical(drawBottom) && drawBottom) {
                        drawBottom <- "lightgray"
                    }
                    waterDepth <- findBottom(x, ignore = ignore)
                    axisBottom <- par("usr")[3]
                    deepestWater <- max(abs(waterDepth$depth))
                }
                ats <- imagep(distance, -depth, z,
                    xlab = if (missing(xlab)) "Distance [km]" else xlab,
                    ylab = if (missing(ylab)) "z [m]" else ylab,
                    ylim = if (missing(ylim)) c(-deepestWater, 0) else ylim,
                    zlim = if (missing(zlim)) c(if (beam[w] %in% c("Sv", "TS")) min(z, na.rm = TRUE) else 0, max(z, na.rm = TRUE)) else zlim,
                    mgp = mgp, mar = mar,
                    tformat = tformat,
                    col = col,
                    drawPalette = drawPalette,
                    debug = debug - 1,
                    zlab = beam[w]
                )
                if (!missing(drawBottom)) {
                    if (is.logical(drawBottom) && drawBottom) {
                        drawBottom <- "white"
                    }
                    ndistance <- length(distance)
                    distance2 <- c(distance[1], distance, distance[ndistance])
                    axisBottom <- par("usr")[3]
                    depth2 <- c(axisBottom, -depth[waterDepth$index], axisBottom)
                    polygon(distance2, depth2, col = drawBottom)
                }
                if (!missing(atTop)) {
                    at <- approx(as.numeric(x[["time"]]), distance, as.numeric(atTop))$y
                    if (missing(labelsTop)) {
                        labelsTop <- format(atTop, format = if ("format" %in% dotsNames) dots$format else "%H:%M:%S")
                    }
                    axis(side = 3, at = at, labels = labelsTop, cex.axis = par("cex"))
                }
                if (drawTimeRange) {
                    timeRange <- range(x[["time"]])
                    label <- paste(timeRange[1], timeRange[2], sep = " to ")
                    mtext(label, side = 3, cex = 0.9 * par("cex"), adj = 0)
                }
                res$xat <- ats$xat
                res$yat <- ats$yat
            } else if (which[w] == 3) {
                oceDebug(debug, "handling which[", w, "]==3\n", sep = "")
                lat <- x[["latitude"]]
                lon <- x[["longitude"]]
                asp <- 1 / cos(mean(range(lat, na.rm = TRUE)) * pi / 180)
                latm <- mean(lat, na.rm = TRUE)
                lonm <- mean(lon, na.rm = TRUE)
                if (missing(radius)) {
                    radius <- max(geodDist(lonm, latm, lon, lat))
                } else {
                    radius <- max(radius, geodDist(lonm, latm, lon, lat))
                }
                km_per_lat_deg <- geodDist(lonm, latm, lonm, latm + 1)
                km_per_lon_deg <- geodDist(lonm, latm, lonm + 1, latm)
                lonr <- lonm + radius / km_per_lon_deg * c(-2, 2)
                latr <- latm + radius / km_per_lat_deg * c(-2, 2)
                plot(lonr, latr,
                    asp = asp, type = "n",
                    xlab = if (missing(xlab)) "Longitude" else xlab,
                    ylab = if (missing(ylab)) "Latitude" else ylab
                )
                xaxp <- par("xaxp")
                xat <- seq(xaxp[1], xaxp[2], length.out = 1 + xaxp[3])
                yaxp <- par("yaxp")
                yat <- seq(yaxp[1], yaxp[2], length.out = 1 + yaxp[3])
                ats <- list(xat = xat, yat = yat)

                if (!missing(coastline)) {
                    coastline <- coastline
                    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)
                        box()
                    } else {
                        lines(coastline[["longitude"]], coastline[["latitude"]], col = "darkgray")
                        lines(coastline[["longitude"]] + 360, coastline[["latitude"]], col = "darkgray")
                    }
                }
                lines(lon, lat, col = if (!is.function(col)) col else "black", lwd = lwd)
            } else {
                stop("In plot,echosounder-method : unknown value of which (",
                    if (is.character(whichOrig[w])) {
                        paste("\"", whichOrig[w], "\"", sep = "")
                    } else {
                        whichOrig[w]
                    }, "); try 1, 2 or 3, or, equivalently, \"zt image\", \"zx image\" or \"map\"",
                    call. = FALSE
                )
            }
        }
        oceDebug(debug, "} # plot.echosounder()\n", unindent = 1, style = "bold")
        invisible(res)
    }
)


#' Read an echosounder File
#'
#' Reads a biosonics echosounder file.  This function was written for and
#' tested with single-beam, dual-beam, and split-beam Biosonics files of type
#' V3, and may not work properly with other file formats.
#'
#' @param file a connection or a character string giving the name of the file
#' to load.
#'
#' @param channel sequence number of channel to extract, for multi-channel
#' files.
#'
#' @param soundSpeed sound speed, in m/s. If not provided, this is calculated
#' using [`swSoundSpeed`]`(35,15,30,eos="unesco")`.  (In theory,
#' it could be calculated using the temperature and salinity that are stored
#' in the data file, but these will just be nominal values, anyway.
#'
#' @param tz character string indicating time zone to be assumed in the data.
#'
#' @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.
#'
#' @return An [echosounder-class] object.
#'
#' @section Bugs: Only the amplitude information (in counts) is determined.  A
#' future version of this function may provide conversion to dB, etc.  The
#' handling of dual-beam and split-beam files is limited.  In the dual-beam
#' cse, only the wide beam signal is processed (I think ... it could be the
#' narrow beam, actually, given the confusing endian tricks being played).  In
#' the split-beam case, only amplitude is read, with the x-axis and y-axis
#' angle data being ignored.
#'
#' @author Dan Kelley, with help from Clark Richards
#'
#' @seealso The documentation for [echosounder-class] explains the
#' structure of `ctd` objects, and also outlines the other functions
#' dealing with them.
#'
#' @references Various echosounder instruments provided by BioSonics are
#' described at the company website, `https://www.biosonicsinc.com/`.  The
#' document listed as reference 1 below was provided to the author of this function in
#' November 2011, which suggests that the data format was not changed since
#' July 2010.
#'
#' 1. Biosonics, 2010.  DT4 Data File Format Specification.  BioSonics
#' Advanced Digital Hydroacoustics. July, 2010.  SOFTWARE AND ENGINEERING
#' LIBRARY REPORT BS&E-2004-07-0009-2.0.
#' @family things related to echosounder data
read.echosounder <- function(
    file,
    channel = 1,
    soundSpeed,
    tz = getOption("oceTz"),
    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, "\"")
        }
    }
    oceDebug(debug, "read.echosounder(file=\"", file, "\", tz=\"", tz, "\", debug=", debug, ") {\n", sep = "", unindent = 1, style = "bold")
    filename <- NULL
    if (is.character(file)) {
        filename <- fullFilename(file)
        file <- file(file, "rb")
        on.exit(close(file))
    }
    if (!inherits(file, "connection")) {
        stop("argument `file' must be a character string or connection")
    }
    if (!isOpen(file)) {
        open(file, "rb")
        on.exit(close(file))
    }
    res <- new("echosounder", filename = filename)
    seek(file, 0, "start")
    seek(file, where = 0, origin = "end")
    fileSize <- seek(file, where = 0)
    oceDebug(debug, "fileSize=", fileSize, "\n")
    buf <- readBin(file, what = "raw", n = fileSize, endian = "little")

    # [1 p9]: "After this comes the main data section of the file.
    # This will typically contain some or all of the following tuples:
    # ping tuples (type 0x0015, 0x001C or 0x001D),
    # time tuples (type 0x000F or 0x0020),
    # position tuples (type 0x000E),
    # navigation string tuples (type 0x0011 or 0x0030),
    # transducer orientation tuples (type 0x0031),
    # bottom pick tuples (type 0x0032),
    # single echoes tuples (type 0x0033), and
    # comment tuples (type 0x0034).
    #
    # [1 sec 3.3 ] describes the file format.  NB: files are little endian.
    #
    # Data are organized as a sequence of tuples, in the following format:
    #   N = 2-byte unsigned int that indicates the size of the tuple.
    #   code = 2-byte code (see table below)
    #   data = N bytes (depends on code)
    #   N6 = 2 bytes that must equal N+6, or the data are corrupted
    #
    # The codes, from the table in [1 sec 3.5] are as follows.
    # The first tuple in a file must have code 0xFFFF, and the
    # second must have code 001E, 0018, or 0001.
    #
    #   0xFFFF Signature (for start of file)
    #   0x001E V3 File Header
    #   0x0018 V2 File Header
    #   0x0001 V1 File Header
    #   0x0012 Channel Descriptor
    #   0x0036 Extended Channel Descriptor
    #   0x0015 Single-Beam Ping
    #   0x001C Dual-Beam Ping
    #   0x001D Split-Beam Ping
    #   0x000F or 0x0020 Time
    #   0x000E Position
    #   0x0011 Navigation String
    #   0x0030 Timestamped Navigation String
    #   0x0031 Transducer Orientation
    #   0x0032 Bottom Pick
    #   0x0033 Single Echoes
    #   0x0034 Comment
    #   0xFFFE End of File
    tuple <- 1
    offset <- 0
    timeSlow <- latitudeSlow <- longitudeSlow <- NULL # accumulate using c() because length unknown
    timeLast <- 0
    # first <- TRUE
    scan <- 1
    # intensity <- list()
    time <- list()
    # samplingDeltat <- 2.4e-05 # a guess, to avoid being unknown if the header cannot be read
    channelNumber <- NULL
    # channelID <- NULL
    channelDeltat <- NULL
    blankedSamples <- 0
    fileType <- "unknown"
    range <- NULL
    beamType <- "unknown"
    # The next three lines are just to prevent code-diagnostic warnings;
    # These matrices are redefined later, when we know the geometry
    a <- matrix(NA_real_, nrow = 1, ncol = 1)
    b <- matrix(NA_real_, nrow = 1, ncol = 1)
    c <- matrix(NA_real_, nrow = 1, ncol = 1)
    # FIXME: find out whether samplesPerPing is always defined prior to use in the code1==0x15 blocks.
    # The Rstudio code-diagnostic complains that this variable is used before being defined,
    # but when I run test code there is no problem, because the variable has been defined. What I
    # do *not* know is whether files will always have these byte groupsing in this order, but at
    # least setting to a zero value is likely to cause an error, if that ever occurs. (I may just
    # need to reorder some code, if problems arise.)
    samplesPerPing <- 0 # overriddent later; here just to prevent code-diagnostic warning
    while (offset < fileSize) {
        # print <- debug && tuple < 200
        N <- .C("uint16_le", buf[offset + 1:2], 1L, res = integer(1), NAOK = TRUE, PACKAGE = "oce")$res
        code1 <- buf[offset + 3]
        code2 <- buf[offset + 4]
        # Next line said 'endian="small"' prior to 2019-12-18, an error kindly pointed out
        # to Dan Kelley in a private email. However, this was not used here, except for an
        # erroneous use (where the test should have been on code2). But note that we do not
        # actually *use* the 'code' defined previously
        # See https://github.com/dankelley/oce/issues/1634
        # code <- readBin(buf[offset+3:4], "integer", size=2, n=1, endian="little", signed=FALSE)
        if (debug > 3) cat("buf[", 3 + offset, "] = code1 = 0x", code1, sep = "")
        if (debug > 3) cat("buf[", 4 + offset, "] = code2 = 0x", code2, sep = "")
        # Interpret code1 and code2, which signal beam type.  These are listed in [1 table 3.5],
        # as follows (note that the table lists in little-endian ordering, so the first two
        # digits are code2 here, and the second two digits are code1 here.
        # * 0x0015 (code1=0x00 code2=0x15): single-beam ping
        # * 0x001c (code1=0x00 code2=0x1c): dual-bam ping
        # * 0x001d (code1=0x00 code2=0x1d): split-beam ping
        #
        # The ordering of the code1 tests is not too systematic here; frequently-encountered
        # codes are placed first, but then it's a bit random.
        if (code2 == 0x00 && (code1 == 0x15 || code1 == 0x1c || code1 == 0x1d)) {
            # single-beam, dual-beam, or split-beam tuple
            thisChannel <- .C("uint16_le", buf[offset + 4 + 1:2], 1L, res = integer(1), NAOK = TRUE, PACKAGE = "oce")$res
            pingNumber <- readBin(buf[offset + 6 + 1:4], "integer", size = 4L, n = 1L, endian = "little")
            pingElapsedTime <- 0.001 * readBin(buf[offset + 10 + 1:4], "integer", size = 4L, n = 1L, endian = "little")
            # message("samplersPerPing=", samplesPerPing)
            ns <- .C("uint16_le", buf[offset + 14 + 1:2], 1L, res = integer(1), NAOK = TRUE, PACKAGE = "oce")$res # number of samples
            if (thisChannel == channelNumber[channel]) {
                if (debug > 3) {
                    cat("buf[", 1 + offset, ", ...] (0x", code1, " single-beam ping)",
                        " scan=", scan,
                        " ping=", pingNumber,
                        " ns=", ns,
                        " channel=", thisChannel,
                        " elapsedTime=", pingElapsedTime,
                        "\n",
                        sep = ""
                    )
                }
                # Note the time reversal in the assignment to the data matrix 'a'
                # FIXME: is it faster to flip the data matrix later?
                if (code1 == 0x15) {
                    # In [1 table 3.5] this case is listed as "0x0015 Single-Beam Ping"
                    tmp <- do_biosonics_ping(buf[offset + 16 + 1:(2 * ns)], samplesPerPing, ns, 0)
                    beamType <- "single-beam"
                } else if (code1 == 0x1c) {
                    # In [1 table 3.5] this case is listed as "0x001C Dual-Beam Ping"
                    tmp <- do_biosonics_ping(buf[offset + 16 + 1:(4 * ns)], samplesPerPing, ns, 1)
                    beamType <- "dual-beam"
                } else if (code1 == 0x1d) {
                    # In [1 table 3.5] this case is listed as "0x001D Split-Beam Ping"
                    # e.g. 01-Fish.dt4 sample file from Biosonics
                    tmp <- do_biosonics_ping(buf[offset + 16 + 1:(4 * ns)], samplesPerPing, ns, 2)
                    beamType <- "split-beam"
                } else {
                    stop("unknown 'tuple' 0x", code1, sep = "")
                }
                a[scan, ] <- rev(tmp$a)
                b[scan, ] <- rev(tmp$b)
                c[scan, ] <- rev(tmp$c)
                time[[scan]] <- timeLast # FIXME many pings between times, so this is wrong
                scan <- scan + 1
                if (debug > 3) cat("channel:", thisChannel, "ping:", pingNumber, "pingElapsedTime:", pingElapsedTime, "\n")
            } else {
                if (debug > 0) {
                    cat("buf[", 1 + offset, ", ...] = 0x", code1,
                        " ping=", pingNumber, " ns=", ns, " channel=", thisChannel, " IGNORED since wrong channel)\n",
                        sep = ""
                    )
                }
            }
        } else if (code2 == 0x00 && (code1 == 0x0f || code1 == 0x20)) {
            # In [1 table 3.5] these cases are listed as "0x000F or 0x0020 Time"
            # time
            timeSec <- readBin(buf[offset + 4 + 1:4], what = "integer", endian = "little", size = 4, n = 1)
            timeSubSec <- .C("biosonics_ss", buf[offset + 10], res = numeric(1), NAOK = TRUE, PACKAGE = "oce")$res
            timeFull <- timeSec + timeSubSec
            timeElapsedSec <- readBin(buf[offset + 10 + 1:4], what = "integer", endian = "little", size = 4, n = 1) / 1e3
            # centisecond = ss & 0x7F [1 sec 4.7]
            timeLast <- timeSec + timeSubSec # used for positioning
            if (debug > 3) cat(sprintf(" time calendar: %s   elapsed %.2f\n", timeFull + as.POSIXct("1970-01-01 00:00:00", tz = "UTC"), timeElapsedSec))
        } else if (code2 == 0x00 && code1 == 0x0e) {
            # In [1 table 3.5] this case is listed as "0x000E Position" position
            lat <- readBin(buf[offset + 4 + 1:4], "integer", endian = "little", size = 4, n = 1) / 6e6
            lon <- readBin(buf[offset + 8 + 1:4], "integer", endian = "little", size = 4, n = 1) / 6e6
            latitudeSlow <- c(latitudeSlow, lat)
            longitudeSlow <- c(longitudeSlow, lon)
            timeSlow <- c(timeSlow, timeLast)
            if (debug > 3) cat(" position", lat, "N", lon, "E\n")
        } else if (code2 == 0xff) {
            if (tuple == 1) {
                if (debug > 1) cat(" file start code\n")
            } else {
                break
            }
        } else if (code1 == 0x11) {
            if (debug > 3) cat(" navigation string IGNORED\n")
        } else if (code1 == 0x12) {
            channelNumber <- c(channelNumber, .C("uint16_le", buf[offset + 4 + 1:2], 1L, res = integer(1), NAOK = TRUE, PACKAGE = "oce")$res)
            ib <- .C("uint16_le", buf[offset + 20 + 1:2], 1L, res = integer(1), NAOK = TRUE, PACKAGE = "oce")$res
            blankedSamples <- ib
            channelDeltat <- c(channelDeltat, 1e-9 * .C("uint16_le", buf[offset + 12 + 1:2], 1L,
                res = integer(1), NAOK = TRUE,
                PACKAGE = "oce"
            )$res)
            pingsInFile <- readBin(buf[offset + 6 + 1:4], "integer", n = 1, size = 4, endian = "little")
            samplesPerPing <- .C("uint16_le", buf[offset + 10 + 1:2], 1L, res = integer(1), NAOK = TRUE, PACKAGE = "oce")$res # ssp [p13 1]
            sp <- 1e-9 * readBin(buf[offset + 12 + 1:2], "integer", n = 1L, size = 2L, endian = "little") # [p13 1] time between samples (ns)
            if (debug > 1) cat("sp: ", sp, " sample period in ns\n", sep = "")
            pud <- readBin(buf[offset + 16 + 1:2], "integer", n = 1L, size = 2L, endian = "little")
            if (debug > 1) cat("pud: ", pud, " pulse duration in us (expect 400 for 01-Fish.dt4)\n", sep = "")
            pp <- .C("uint16_le", buf[offset + 18 + 1:2], 1L, res = integer(1), NAOK = TRUE, PACKAGE = "oce")$res # [p13 1] ping period (ms)
            if (debug > 1) cat("pp: ", pp, "\n", sep = "")
            ib <- .C("uint16_le", buf[offset + 20 + 1:2], 1L, res = integer(1), NAOK = TRUE, PACKAGE = "oce")$res # [p13 1]
            if (debug > 1) cat("ib: ", ib, "\n", sep = "")
            # next 2 bytes contain u2, ignored
            th <- 0.01 * readBin(buf[offset + 24 + 1:2], "integer", n = 1L, size = 2L, endian = "little") # [p13 1]
            if (debug > 1) cat("th: ", th, " data collection threshold in dB (expect -65 for 01-Fish.dt4)\n", sep = "")
            rxee <- buf[offset + 26 + 1:128] # [1 p13] # receiver EEPROM image (FIXME: why is serialnum out by 1? INDEX question)
            # corr <- .C("int16_le", buf[offset+282+1:2], 1L, res=integer(1), NAOK=TRUE, PACKAGE="oce")$res # [p13 1]
            corr <- 0.01 * readBin(buf[offset + 282 + 1:2], "integer", n = 1L, size = 2L, endian = "little") # [p13 1]
            if (debug > 1) cat("corr: ", corr, " user-defined calibration correction in dB (expect 0 for 01-Fish.dt4)\n", sep = "")
            if (1 == length(channelNumber)) {
                # get space
                a <- matrix(NA_real_, nrow = pingsInFile, ncol = samplesPerPing)
                b <- matrix(NA_real_, nrow = pingsInFile, ncol = samplesPerPing)
                c <- matrix(NA_real_, nrow = pingsInFile, ncol = samplesPerPing)
            }
            if (debug > 3) {
                cat(
                    " channel descriptor ",
                    " number=", tail(channelNumber, 1),
                    " blankedSamples=", blankedSamples,
                    " dt=", tail(channelDeltat, 1),
                    " pingsInFile=", pingsInFile,
                    " samplesPerPing=", samplesPerPing,
                    "\n"
                )
            }
        } else if (code1 == 0x30) {
            if (debug > 3) cat(" time-stamped navigation string IGNORED\n")
        } else if (code1 == 0xff && tuple > 1) {
            if (debug > 3) cat("  EOF\n")
        } else if (code1 == 0x1e) {
            if (debug > 1) cat(" V3 file header\n")
            fileType <- if (buf[offset + 1] == 0x10 && buf[offset + 2] == 0x00) "DT4 v2.3" else "DT4 pre v2.3"
            res@metadata$temperature <- 0.01 * .C("uint16_le", buf[offset + 8 + 1:2], 1L, res = integer(1), NAOK = TRUE, PACKAGE = "oce")$res
            if (debug > 1) cat("temperature=", res@metadata$temperature, "degC (expect 14 for 1-Fish.dt4)\n")
            res@metadata$salinity <- 0.01 * .C("uint16_le", buf[offset + 10 + 1:2], 1L, res = integer(1), NAOK = TRUE, PACKAGE = "oce")$res
            if (debug > 1) cat("salinity=", res@metadata$salinity, "PSU (expect 30 for 1-Fish.dt4)\n")
            # res@metadata$soundSpeed <- swSoundSpeed(res@metadata$salinity, res@metadata$temperature, 30,
            #                                        eos="unesco")
            res@metadata$transmitPower <- 0.01 * .C("uint16_le", buf[offset + 12 + 1:2], 1L, res = integer(1), NAOK = TRUE, PACKAGE = "oce")$res
            if (debug > 1) cat("transmitPower=", res@metadata$transmitPower, "(expect 0 for 1-Fish.dt4; called mTransmitPower)\n")
            res@metadata$tz <- readBin(buf[offset + 16 + 1:2], "integer", size = 2)
            if (debug > 1) cat("tz=", res@metadata$tz, "(expect -420 for 1-Fish.dt4)\n")
            dst <- readBin(buf[offset + 18 + 1:2], "integer", size = 2)
            res@metadata$dst <- dst != 0
            if (debug > 1) cat("dst=", res@metadata$dst, "(daylight saings time) expect TRUE for 1-Fish.dt4)\n")
        } else if (code1 == 0x18) {
            warning("Biosonics file of type 'V2' detected ... errors may crop up")
            fileType <- "V2"
        } else if (code1 == 0x01) {
            warning("Biosonics file of type 'V1' detected ... errors may crop up")
            fileType <- "V1"
        } else if (code1 == 0x32) {
            if (debug > 3) cat(" bottom pick tuple [1 sec 4.12] ")
            # thisChannel <- .C("uint16_le", buf[offset+4+1:2], 1L, res=integer(1))$res
            # thisPing <- .C("uint16_le", buf[offset+4+1:2], 1L, res=integer(1))$res
            foundBottom <- .C("uint16_le", buf[offset + 14 + 1:2], 1L, res = integer(1), NAOK = TRUE, PACKAGE = "oce")$res
            if (foundBottom) {
                thisRange <- readBin(buf[offset + 20 + 1:4], "numeric", size = 4, n = 1, signed = TRUE, endian = "little")
            } else {
                thisRange <- NA
            }
            range <- c(range, thisRange)
            if (debug > 3) cat(" thisRange:", thisRange)
        } else if (code1 == 0x36) {
            if (debug > 3) cat(" extended channel descriptor IGNORED\n")
        } else if (code1 == 0x33) {
            if (debug > 3) cat(" single echo tuple IGNORED ... see p26 of DT4_format_2010.pdf\n")
        } else if (code1 == 0x34) {
            if (debug > 1) cat(" comment tuple [1 sec 4.14 p28]\n")
            # FIXME: other info could be gleaned from the comment, if needed
            numbytes <- .C("uint16_le", buf[offset + 34:35], 1L, res = integer(1), NAOK = TRUE, PACKAGE = "oce")$res
            if (debug > 1) cat("numbytes:", numbytes, " ... NOTHING ELSE DECODED in this version of oce.\n")
        } else {
            if (debug > 3) cat(" unknown code IGNORED\n")
        }
        N6 <- .C("uint16_le", buf[offset + N + 5:6], 1L, res = integer(1), NAOK = TRUE, PACKAGE = "oce")$res
        if (N6 != N + 6) {
            stop("error reading tuple number ", tuple, " (mismatch in redundant header-length flags)")
        }
        offset <- offset + N + 6
        tuple <- tuple + 1
    }
    res@metadata$beamType <- beamType
    res@metadata$channel <- channel
    res@metadata$fileType <- fileType
    res@metadata$blankedSamples <- blankedSamples
    if (missing(soundSpeed)) {
        res@metadata$soundSpeed <- swSoundSpeed(35, 10, 30, eos = "unesco")
    } else {
        res@metadata$soundSpeed <- soundSpeed
    }
    res@metadata$soundSpeed <- if (missing(soundSpeed)) swSoundSpeed(35, 10, 30, eos = "unesco") else soundSpeed
    res@metadata$samplingDeltat <- channelDeltat[1] # nanoseconds
    res@metadata$pingsInFile <- pingsInFile
    res@metadata$samplesPerPing <- samplesPerPing
    # channel info, with names matching [1 p13]
    # res@metadata$samplingDeltat <- sp
    res@metadata$pulseDuration <- pud
    res@metadata$pp <- pp
    res@metadata$ib <- ib
    res@metadata$th <- th
    res@metadata$rxee <- rxee
    res@metadata$correction <- corr
    depth <- rev(blankedSamples + seq(1, dim(a)[2])) * res@metadata$soundSpeed * res@metadata$samplingDeltat / 2
    # test: for 01-Fish.dt4, have as follows:
    #     <mPingBeginRange_m>0.988429069519043</mPingBeginRange_m>
    #     <mPingEndRange_m>64.787033081054688</mPingEndRange_m>
    # but we get as follows:
    #     range(d[["depth"]])
    #     [1]  1.001714 64.485323
    # i.e. a 12cm error at top and a 20cm error at bottom.  Note that
    #    > diff(d[["depth"]][2:1])
    #    [1] 0.01788775
    # FIXME: check depth mismatch relates to (a) sound speed or (b) geometry.  (Small error; low priority.)
    time <- as.numeric(time)
    # interpolate to "fast" latitude and longitude, after extending to ensure spans
    # enclose the ping times.
    n <- length(latitudeSlow)
    # t <- c(2*timeSlow[1]-timeSlow[2], timeSlow, 2*timeSlow[n] - timeSlow[n-1])
    approx2 <- function(x, y, xout) {
        nxout <- length(xout)
        before <- y[1] + (y[2] - y[1]) * (xout[1] - x[1]) / (x[2] - x[1])
        after <- y[n - 1] + (y[n] - y[n - 1]) * (xout[nxout] - x[n - 1]) / (x[n] - x[n - 1])
        approx(c(xout[1], x, xout[nxout]), c(before, y, after), xout)$y
    }
    latitude <- approx2(timeSlow, latitudeSlow, time)
    longitude <- approx2(timeSlow, longitudeSlow, time)
    res@metadata$transducerSerialNumber <- readBin(rxee[2 + 1:8], "character") # [1 p16] offset=2 length 8
    if (debug > 1) cat("transducerSerialNumber '", res@metadata$transducerSerialNumber, "' (expect DT600085 for 01-Fish.dt4)\n", sep = "")
    res@metadata$calibrationTime <- numberAsPOSIXct(readBin(rxee[36 + 1:4], "integer"), tz = "UTC") # [1 p16] offset=36
    if (debug > 1) cat("calibrationTime: ", format(res@metadata$calibrationTime), " (expect Thu Apr 13 02:36:38 2000 for 01-Fish.dt4)\n", sep = "")
    # res@metadata$sl <- 0.1 * .C("int16_le", rxee[58+1:2], 1L, res=integer(1), NAOK=TRUE, PACKAGE="oce")$res # [p16 1]
    res@metadata$sourceLevel <- 0.1 * readBin(rxee[58 + 1:2], "integer", n = 1L, size = 2L, endian = "little") # [p16 1]
    if (debug > 1) cat("sl=", res@metadata$sourceLevel, " (expect 220 for 01-Fish.dt4 source level SourceLevel_dBuPa)\n", sep = "")
    # res@metadata$rs <- 0.1 * .C("int16_le", rxee[1+64+1:2], 1L, res=integer(1), NAOK=TRUE, PACKAGE="oce")$res # [p16 1]
    res@metadata$receiverSensitivity <- 0.1 * readBin(rxee[64 + 1:2], "integer", n = 1L, size = 2) # [p16 1]
    if (debug > 1) cat("rs=", res@metadata$receiverSensitivity, " receiver sensitivity in dB (counts/micro Pa) (expect -58.8 for 01-Fish.dt4)\n")
    res@metadata$trtype <- .C("uint16_le", rxee[84 + 1:2], 1L, res = integer(1), NAOK = TRUE, PACKAGE = "oce")$res # [p16 1]
    if (debug > 1) {
        cat(
            "trtype=", res@metadata$trtype, " hardware [not deployment] transducer type ",
            "(0 single, 3 dual, 4 split) expect 1 for 01-Fish.dt4\n"
        )
    }
    res@metadata$frequency <- readBin(rxee[86 + 1:4], "integer", n = 1L, size = 4) # [p16 1]
    if (debug > 1) cat("fq=", res@metadata$frequency, " transducer frequency, Hz (expect 208000 for 01-Fish.dt4)\n", sep = "")
    res@metadata$beamwidthX <- 0.1 * as.numeric(rxee[1 + 100]) # [1 p16] offset=100
    res@metadata$beamwidthY <- 0.1 * as.numeric(rxee[1 + 101]) # [1 p16] offset=101
    if (debug > 1) cat("beamwidthX=", res@metadata$beamwidthX, " (expect 6.5 for 01-Fish.dt4; called BeamWidthX_deg)\n", sep = "")
    if (debug > 1) cat("beamwidthY=", res@metadata$beamwidthY, " (expect 6.5 for 01-Fish.dt4; called BeamWidthY_deg)\n", sep = "")
    res@data <- list(
        timeSlow = timeSlow + as.POSIXct("1970-01-01 00:00:00", tz = "UTC"),
        latitudeSlow = latitudeSlow,
        longitudeSlow = longitudeSlow,
        depth = depth,
        # range=range,
        time = time + as.POSIXct("1970-01-01 00:00:00", tz = "UTC"),
        latitude = latitude,
        longitude = longitude,
        a = a, b = b, c = c
    )
    if (res@metadata$beamType == "single-beam") {
        res@data$b <- NULL
        res@data$c <- NULL
    }
    names <- names(res@data)
    if ("latitude" %in% names) res@metadata$units$latitude <- list(unit = expression(degree * N), scale = "")
    if ("longitude" %in% names) res@metadata$units$longitude <- list(unit = expression(degree * E), scale = "")
    if ("latitudeSlow" %in% names) res@metadata$units$latitudeSlow <- list(unit = expression(degree * N), scale = "")
    if ("longitudeSlow" %in% names) res@metadata$units$longitudeSlow <- list(unit = expression(degree * E), scale = "")
    if ("depth" %in% names) res@metadata$units$depth <- list(unit = expression(m), scale = "")
    if (!missing(processingLog)) {
        res@processingLog <- processingLogItem(processingLog)
    }
    res@processingLog <- processingLogAppend(
        res@processingLog,
        paste("read.echosounder(\"", filename, "\", channel=", channel, ", soundSpeed=",
            if (missing(soundSpeed)) "(missing)" else soundSpeed, ", tz=\"", tz, "\", debug=", debug, ", processingLog)",
            sep = ""
        )
    )
    .C("biosonics_free_storage", PACKAGE = "oce") # clear temporary storage space
    oceDebug(debug, "} read.echosounder()\n", sep = "", unindent = 1, style = "bold")
    res
}