Nothing
# 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
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.