Nothing
# vim: tw=80 shiftwidth=4 softtabstop=4 expandtab:
#' Read an adp File in Sontek Format
#'
#' Read a Sontek acoustic-Doppler profiler file (see reference 1).
#'
#' @param despike if `TRUE`, [despike()] will be used to clean
#' anomalous spikes in heading, etc.
#'
#' @param type A character string indicating the type of instrument.
#'
#' @template adpTemplate
#'
#' @template encodingIgnoredTemplate
#'
#' @references
#' 1. Information about Sontek profilers is available at https://www.sontek.com.
#'
#' @template adReadingMethodTemplate
#'
#' @author Dan Kelley and Clark Richards
#'
#' @family things related to adp data
#' @family functions that read adp data
read.adp.sontek <- function(
file, from = 1, to, by = 1, tz = getOption("oceTz"),
longitude = NA, latitude = NA, type = c("adp", "pcadp"), encoding = NA,
monitor = FALSE, despike = FALSE, processingLog,
debug = getOption("oceDebug"), ...) {
if (missing(file)) {
stop("must supply 'file'")
}
if (is.character(file)) {
if (!file.exists(file)) {
stop("cannot find file \"", file, "\"")
}
if (0L == file.info(file)$size) {
stop("empty file \"", file, "\"")
}
}
if (!interactive()) {
monitor <- FALSE
}
missing.to <- missing(to)
# In this function, comments in [] refer to logical page number of ADPManual_v710.pd; add 14 for file page number
profileStart <- NULL # prevent scope warning from rstudio; defined later anyway
bisectSontekAdp <- function(buf, t.find, add = 0, debug = 0) {
oceDebug(debug, "bisectSontekAdp(t.find=", format(t.find), ", add=", add, ", debug=", debug, ")\n")
len <- length(profileStart)
lower <- 1
upper <- len
passes <- floor(10 + log(len, 2)) # won't need this many; only do this to catch coding errors
for (pass in 1:passes) {
middle <- floor((upper + lower) / 2)
oceDebug(debug - 1, "pass=", pass, "middle=", middle, "\n")
oceDebug(debug - 1, "data=", buf[profileStart[middle] + seq(0, 100, 1)], "\n")
# oceDebug(debug-1, "profileStart=", profileStart[1:5], "...; profileStart2=", profileStart2[1:5], "...\n")
oceDebug(debug - 1, "profileStart=", profileStart[1:5], "\n")
# year <- readBin(buf[profileStart2[middle]+18], "integer", size=2, signed=FALSE, endian="little")
year <- readBin(buf[profileStart[middle] + 18:19], "integer", size = 2, signed = FALSE, endian = "little")
day <- as.integer(buf[profileStart[middle] + 20])
month <- as.integer(buf[profileStart[middle] + 21])
min <- as.integer(buf[profileStart[middle] + 22])
hour <- as.integer(buf[profileStart[middle] + 23])
# SIG p82 C code suggests sec100 comes before second.
sec100 <- as.integer(buf[profileStart[middle] + 24]) # FIXME: determine whether this is 1/100th second
sec <- as.integer(buf[profileStart[middle] + 25])
t <- ISOdatetime(year, month, day, hour, min, sec + sec100 / 100, tz = tz)
oceDebug(
debug, "t=", format(t),
" [year=", year, " month=", month, " day=", day, " hour=", hour, " sec=", sec, "sec100=", sec100, "]\n"
)
if (t.find < t) {
upper <- middle
} else {
lower <- middle
}
if (upper - lower < 2) {
break
}
oceDebug(debug, "middle=", middle, " lower=", lower, " upper=", upper, " pass=", pass, " of max=", passes, "\n")
}
middle <- middle + add # may use add to extend before and after window
if (middle < 1) {
middle <- 1
}
if (middle > len) {
middle <- len
}
t <- ISOdatetime(readBin(buf[profileStart[middle] + 18:19], "integer", size = 2, signed = FALSE, endian = "little"), # year
as.integer(buf[profileStart[middle] + 21]), # month
as.integer(buf[profileStart[middle] + 20]), # day
as.integer(buf[profileStart[middle] + 23]), # hour
as.integer(buf[profileStart[middle] + 22]), # min
as.integer(buf[profileStart[middle] + 25]), # sec
tz = tz
)
oceDebug(debug, "result: t=", format(t), " at profileStart[", middle, "]=", profileStart[middle], "\n")
return(list(index = middle, time = t)) # index is within vsd
}
oceDebug(debug, "read.adp.sontek(...,from=", from, ",to=", if (missing.to) "(missing)" else to, ",by=", by, "type=", type, "...)\n")
# parameters <- list(profile.byte1 = 0xa5, profile.byte2=0x10, profile.headerLength=80)
if (is.character(file)) {
if (0 == file.info(file)$size) {
stop("empty 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)) {
filename <- "(connection)"
open(file, "rb")
on.exit(close(file))
}
type <- match.arg(type)
res <- new("adp")
seek(file, 0, "end")
fileSize <- seek(file, 0, "start")
oceDebug(debug, "file", filename, "has", fileSize, "bytes\n")
buf <- readBin(file, what = "raw", n = fileSize, endian = "little")
# See if there is a header here. (May not be, since the file may have been chopped
# into parts by a deck unit.)
frequency <- NA
if (buf[1] == 0x10 && buf[2] == 0x02 && 96 == readBin(buf[3:4], "integer", n = 1, size = 2, signed = FALSE, endian = "little")) {
oceDebug(debug, "have a header, but ignoring it for now\n")
# Comments like [p83] refer to logical page number of ADPManual_v710.pd; add 14 for file page number
cpuSoftwareVerNum <- as.integer(buf[13]) / 10 # CPUSoftwareVerNum [p83]
dspSoftwareVerNum <- as.integer(buf[14]) / 10 # DSPSoftwareVerNum [p83]
boardRev <- readBin(buf[15], "character", n = 1, size = 1, signed = TRUE) # BoardRev [p83]
serialNumber <- readBin(buf[16:25], "character")
oceDebug(debug, "serialNumber=", serialNumber, "\n")
adp.type <- readBin(buf[26], what = "integer", n = 1, size = 1) # 0-3; 1=1.5; 2-750; 3-500 [p83]
oceDebug(debug, "adp.type=", adp.type, "\n")
frequency <- switch(adp.type + 1,
3000,
1500,
750,
500,
250
)
oceDebug(debug, "frequency=", frequency, "\n")
nbeams <- as.integer(buf[27])
oceDebug(debug, "nbeams=", nbeams, "\n")
beam.geometry <- as.integer(buf[28])
oceDebug(
debug, "beam.geometry=", beam.geometry,
"; 0 (2 beams); 1 (3 beams), 2 (4 beams with 1 vertical), 3 (4 beams, Janus)\n"
)
slant.angle <- readBin(buf[29:30], "integer", n = 1, size = 2, signed = FALSE) / 10
oceDebug(debug, "slant.angle=", slant.angle, "\n")
orientation <- switch(as.integer(buf[31]) + 1,
"down",
"up",
"side"
)
oceDebug(debug, "orientation=", orientation, "\n")
compass.installed <- switch(as.integer(buf[32]) + 1,
FALSE,
TRUE
) # nolint (variable not used)
recorder.installed <- switch(as.integer(buf[33]) + 1,
FALSE,
TRUE
) # nolint (variable not used)
temp.installed <- switch(as.integer(buf[34]) + 1,
FALSE,
TRUE
) # nolint (variable not used)
press.installed <- switch(as.integer(buf[35]) + 1,
FALSE,
TRUE
) # nolint (variable not used)
# 36 = spare
# 37 int[16], so I guess 2-byte ints, signed?
} else {
cpuSoftwareVerNum <- dspSoftwareVerNum <- boardRev <- "?"
adp.type <- nbeams <- slant.angle <- orientation <- "?"
compass.installed <- recorder.installed <- temp.installed <- press.installed <- "?"
}
# profileStart <- .Call("match2bytes", buf, parameters$profile.byte1, parameters$profile.byte2, FALSE)
# profileStart <- .Call("ldc_sontek_adp", buf, 0, 0, 0, 1, -1) # no ctd, no gps, no bottom-track; pcadp; all data
profileStart <- do_ldc_sontek_adp(buf, 0, 0, 0, 1, -1) # no ctd, no gps, no bottom-track; pcadp; all data
profileStart2 <- sort(c(profileStart, profileStart + 1)) # use this to subset for 2-byte reads
oceDebug(debug, "first 10 profileStart:", profileStart[1:10], "\n")
oceDebug(debug, "first 100 bytes of first profile:", paste(buf[profileStart[1]:(99 + profileStart[1])], collapse = " "), "\n")
# Examine the first profile to get numberOfBeams, etc.
headerLength <- 80 # FIXME: should this really be hard-wired??
s <- profileStart[1]
# Only read (important) things that don't change profile-by-profile
numberOfBeams <- as.integer(buf[s + 26])
oceDebug(debug, "numberOfBeams=", numberOfBeams, "\n")
if (numberOfBeams != 3) {
stop("there should be 3 beams, but the file indicates ", numberOfBeams)
}
orientation <- as.integer(buf[s + 27])
oceDebug(debug, "orientation=", orientation, "\n")
temp.mode <- as.integer(buf[s + 28])
oceDebug(debug, "temp.mode=", temp.mode, "\n")
originalCoordinate <- as.integer(buf[s + 29])
oceDebug(debug, "originalCoordinate=", originalCoordinate, "\n")
numberOfCells <- readBin(buf[s + 30:31], "integer", n = 1, size = 2, endian = "little", signed = FALSE)
oceDebug(debug, "numberOfCells=", numberOfCells, "\n")
cellSize <- readBin(buf[s + 32:33], "integer", n = 1, size = 2, endian = "little", signed = FALSE) / 100 # metres
oceDebug(debug, "cellSize=", cellSize, "m\n")
blankingDistance <- readBin(buf[s + 34:35], "integer", n = 1, size = 2, endian = "little", signed = FALSE) / 100 # metres
oceDebug(debug, "blankingDistance=", blankingDistance, "m\n")
soundSpeed <- readBin(buf[s + 60:61], "integer", n = 1, size = 2, endian = "little", signed = FALSE) / 10
oceDebug(debug, "soundSpeed=", soundSpeed, "m/s\n")
profilesInFile <- length(profileStart)
# id <- buf[profileStart]
bytesPerProfile <- diff(profileStart[1:2])
oceDebug(debug, "bytesPerProfile=", bytesPerProfile, "\n")
# File time range and deltat
measurementStart <- ISOdatetime(readBin(buf[profileStart[1] + 18:19], "integer", n = 1, size = 2, signed = FALSE, endian = "little"), # year
as.integer(buf[profileStart[1] + 21]), # month
as.integer(buf[profileStart[1] + 20]), # day
as.integer(buf[profileStart[1] + 23]), # hour
as.integer(buf[profileStart[1] + 22]), # min
as.integer(buf[profileStart[1] + 25]) + 0.01 * as.integer(buf[profileStart[1] + 24]), # sec (decimal)
tz = tz
)
oceDebug(debug, "measurementStart=", format(measurementStart), "\n")
oceDebug(debug, "length(profileStart)=", length(profileStart), " [FIXME: if to not given, use this??]\n")
measurementEnd <- ISOdatetime(readBin(buf[profileStart[profilesInFile] + 18:19], "integer", n = 1, size = 2, signed = FALSE, endian = "little"), # year
as.integer(buf[profileStart[profilesInFile] + 21]), # month
as.integer(buf[profileStart[profilesInFile] + 20]), # day
as.integer(buf[profileStart[profilesInFile] + 23]), # hour
as.integer(buf[profileStart[profilesInFile] + 22]), # min
as.integer(buf[profileStart[profilesInFile] + 25]) + 0.01 * as.integer(buf[profileStart[1] + 24]), # sec (decimal)
tz = tz
)
oceDebug(debug, "sampling.end=", format(measurementEnd), "\n")
measurementDeltat <- as.numeric(ISOdatetime(readBin(buf[profileStart[2] + 18:19], "integer", n = 1, size = 2, signed = FALSE, endian = "little"), # year
as.integer(buf[profileStart[2] + 21]), # month
as.integer(buf[profileStart[2] + 20]), # day
as.integer(buf[profileStart[2] + 23]), # hour
as.integer(buf[profileStart[2] + 22]), # min
as.integer(buf[profileStart[2] + 25]) + 0.01 * as.integer(buf[profileStart[1] + 24]), # sec
tz = tz
)) - as.numeric(measurementStart)
oceDebug(debug, "sampling.deltat=", format(measurementDeltat), "\n")
# Window data buffer, using bisection in case of a variable number of vd between sd pairs.
if (inherits(from, "POSIXt")) {
if (!inherits(to, "POSIXt")) {
stop("if 'from' is POSIXt, then 'to' must be, also")
}
fromPair <- bisectSontekAdp(buf, from, -1, debug - 1)
from <- fromIndex <- fromPair$index
toPair <- bisectSontekAdp(buf, to, 1, debug - 1)
to <- to.index <- toPair$index
oceDebug(
debug, " from=", format(fromPair$t), " yields profileStart[", fromIndex, "]\n",
" to =", format(toPair$t), " yields profileStart[", to.index, "]\n",
" by=", by, "s\n",
"profileStart[1:10]=", profileStart[1:10], "\n",
"profileStart[", fromPair$index, "]=", profileStart[fromPair$index], "at time", format(fromPair$t), "\n",
"profileStart[", toPair$index, "]=", profileStart[toPair$index], "at time", format(toPair$t), "\n"
)
# FIXME next line reads year incorrectly
two.times <- ISOdatetime(readBin(buf[profileStart[1:2] + 18:19], "integer", size = 2, signed = FALSE, endian = "little"), # year
as.integer(buf[profileStart[1:2] + 21]), # month
as.integer(buf[profileStart[1:2] + 20]), # day
as.integer(buf[profileStart[1:2] + 23]), # hour
as.integer(buf[profileStart[1:2] + 22]), # min
as.integer(buf[profileStart[1:2] + 25]) + 0.01 * as.integer(buf[profileStart[1] + 24]), # sec
tz = tz
)
dt <- as.numeric(difftime(two.times[2], two.times[1], units = "secs"))
oceDebug(debug, "dt=", dt, "s; at this stage, by=", by, "(not interpreted yet)\n")
profileStart <- profileStart[profileStart[fromIndex] < profileStart & profileStart < profileStart[to.index]]
if (is.character(by)) {
by <- floor(0.5 + ctimeToSeconds(by) / dt)
}
oceDebug(debug, "by=", by, "profiles (after change)\n")
profileStart <- profileStart[seq(1, length(profileStart), by = by)]
oceDebug(debug, "dt=", dt, "\n", "by=", by, "profileStart[1:10] after indexing:", profileStart[1:10], "\n")
} else {
fromIndex <- from
if (missing.to) {
to <- length(profileStart)
}
to.index <- to
if (to.index < 1 + fromIndex) {
stop("need more separation between from and to")
}
if (is.character(by)) {
stop("cannot have string for 'by' if 'from' and 'to' are integers")
}
profileStart <- profileStart[seq(from = from, to = to, by = by)]
oceDebug(debug, "profileStart[1:10] after indexing:", profileStart[1:10], "\n")
}
profilesToRead <- length(profileStart)
oceDebug(debug, "profilesInFile=", profilesInFile, "; profilesToRead=", profilesToRead, "\n")
profileStart2 <- sort(c(profileStart, profileStart + 1)) # use this to subset for 2-byte reads
year <- readBin(buf[profileStart2 + 18], "integer", n = profilesToRead, size = 2, endian = "little", signed = FALSE)
day <- as.integer(buf[profileStart + 20])
month <- as.integer(buf[profileStart + 21])
minute <- as.integer(buf[profileStart + 22])
hour <- as.integer(buf[profileStart + 23])
sec100 <- as.integer(buf[profileStart + 24])
second <- as.integer(buf[profileStart + 25])
time <- ISOdatetime(year, month, day, hour, minute, second + sec100 / 100, tz = tz)
rm(year, day, month, minute, hour, sec100, second)
temperature <- readBin(buf[profileStart2 + 46], "integer", n = profilesToRead, size = 2, endian = "little", signed = TRUE) / 100
oceDebug(debug, "temperature[1:10]=", temperature[1:10], "\n")
pressure <- readBin(buf[profileStart2 + 48], "integer", n = profilesToRead, size = 2, endian = "little", signed = FALSE) / 100
# FIXME: pressure (+else?) is wrong. Need to count bytes on p84 of ADPManual to figure out where to look [UGLY]
oceDebug(debug, "pressure[1:10]=", pressure[1:10], "\n")
heading <- readBin(buf[profileStart2 + 40], "integer", n = profilesToRead, size = 2, endian = "little", signed = TRUE) / 10
pitch <- readBin(buf[profileStart2 + 42], "integer", n = profilesToRead, size = 2, endian = "little", signed = TRUE) / 10
roll <- readBin(buf[profileStart2 + 44], "integer", n = profilesToRead, size = 2, endian = "little", signed = TRUE) / 10
oceDebug(debug, "time[1:10]=", format(time[1:10]), "\n")
v <- array(numeric(), dim = c(profilesToRead, numberOfCells, numberOfBeams))
a <- array(raw(), dim = c(profilesToRead, numberOfCells, numberOfBeams))
q <- array(raw(), dim = c(profilesToRead, numberOfCells, numberOfBeams))
nd <- numberOfCells * numberOfBeams
oceDebug(debug, "nd=", nd, "; headerLength=", headerLength, "\n")
if (type == "pcadp") {
nbeamMax <- 4 # Max numberOfBeams, not actual number
headerLength <- headerLength + 2 * (8 + nbeamMax) + 2 * nbeamMax + nbeamMax
# Below is C code from Sontek, for pulse-coherent adp (2-byte little-endian
# integers). FIXME: should perhaps read these things, but this is not a
# high priority, since in the data file for which the code was originally
# developed, all distances were set to 123 mm and all velocities to
# 9999 mm/s, suggestive of insignificant, place-holder values.
#
# typedef struct
# {
# unsigned int ResLag; /* in mm Used for single cell */
# unsigned int ResUa; /* in mm/s Ambiguity resolution */
# unsigned int ResStart; /* in mm Position of resolve */
# unsigned int ResLength; /* in mm cell */
# unsigned int PrfLag; /* in mm Used for full profile */
# unsigned int PrfUa; /* in mm/s */
# unsigned int PrfStart; /* in mm Position/Length of first*/
# unsigned int PrfLength; /* in mm cell in profile */
# unsigned int Range[MAX_BEAMS]; /* in mm Range to boundary */
# int Ures[MAX_BEAMS]; /* in mm/s Velocities from Resolve */
# /* lag */
# unsigned char Cres[MAX_BEAMS]; /* in % Correlations from */
# /* resolve lag */
# } PCrecordType;
}
velocityScale <- 1e-3
if (profilesToRead > 0) {
for (i in 1:profilesToRead) {
# nolint start object_usage_linter
v_ <- velocityScale * matrix(
readBin(buf[profileStart[i] + headerLength + seq(0, 2 * nd - 1)],
"integer",
n = nd, size = 2, signed = TRUE, endian = "little"
),
ncol = numberOfBeams, byrow = FALSE
)
a_ <- matrix(buf[profileStart[i] + headerLength + 2 * nd + seq(0, nd - 1)], ncol = numberOfBeams, byrow = FALSE)
q_ <- matrix(buf[profileStart[i] + headerLength + 3 * nd + seq(0, nd - 1)], ncol = numberOfBeams, byrow = FALSE)
# nolint end object_usage_linter
for (b in 1:numberOfBeams) {
# FIXME: probably could be speeded up
v[i, , b] <- v_[, b]
a[i, , b] <- a_[, b]
q[i, , b] <- q_[, b]
}
if (monitor) {
cat(".")
if (!(i %% 50)) {
cat(i, "\n")
}
}
}
rm(v_, a_, q_)
if (monitor) {
cat("\nRead", profilesToRead, "of the", profilesInFile, "profiles in", filename, "\n")
}
if (type == "pcadp") {
v <- v / 10 # it seems pcadp is in 0.1mm/s
}
} else {
stop("please request to read *some* profiles")
}
# interpolate headings (which may be less frequent than profiles ... FIXME: really???)
nheading <- length(heading)
nv <- dim(v)[1]
if (nheading != nv) {
warning("read.adp.sontek() interpolating ", nheading, " heading/pitch/roll values to the ", nv, " velocity profiles")
oceDebug(debug, "BEFORE: length(heading)=", nheading, ", nv=", nv, "\n")
# xout <- seq(1, nheading, length.out=nv)
heading <- approx(1:nheading, heading, seq(1, nheading, length.out = nv))$y
# print(data.frame(xout=xout, heading=heading))
pitch <- approx(1:nheading, pitch, seq(1, nheading, length.out = nv))$y
roll <- approx(1:nheading, roll, seq(1, nheading, length.out = nv))$y
oceDebug(debug, "AFTER: length(heading)=", length(heading), "\n")
}
res@data <- list(
v = v, a = a, q = q,
distance = seq(blankingDistance, by = cellSize, length.out = numberOfCells),
time = time,
temperature = temperature,
pressure = pressure,
heading = heading, pitch = pitch, roll = roll
)
oceDebug(debug, "slant.angle=", slant.angle, "; type=", type, "\n")
beamAngle <- if (slant.angle == "?") 25 else slant.angle
res@metadata$manufacturer <- "sontek"
res@metadata$type <- type
res@metadata$filename <- filename
res@metadata$serialNumber <- if (exists("serialNumber")) serialNumber else "?"
res@metadata$longitude <- longitude
res@metadata$latitude <- latitude
res@metadata$numberOfSamples <- dim(v)[1]
res@metadata$numberOfCells <- dim(v)[2]
res@metadata$numberOfBeams <- dim(v)[3]
res@metadata$velocityResolution <- velocityScale
res@metadata$velocityMaximum <- velocityScale * 2^15
res@metadata$measurementStart <- measurementStart
res@metadata$measurementEnd <- measurementEnd
res@metadata$measurementDeltat <- measurementDeltat
res@metadata$frequency <- frequency
res@metadata$cpuSoftwareVerNum <- cpuSoftwareVerNum
res@metadata$dspSoftwareVerNum <- dspSoftwareVerNum
res@metadata$boardRev <- boardRev
res@metadata$originalCoordinate <- c("beam", "xyz", "enu", "other")[originalCoordinate + 1]
res@metadata$oceCoordinate <- c("beam", "xyz", "enu", "other")[originalCoordinate + 1]
res@metadata$beamAngle <- beamAngle
res@metadata$oceBeamUnspreaded <- FALSE
res@metadata$orientation <- if (1 == orientation) "upward" else "downward"
if (numberOfBeams == 3) {
if (res@metadata$orientation == "upward") {
# S <- 1 / (3 * sin(25 * pi / 180)) # 0.7887339
# CS <- 1 / cos(30*pi/180) / sin(25*pi/180) / 2 # 1.366127 (30deg from 3-beam pattern)
# C <- 1 / (3 * cos(25 * pi / 180)) # 0.3677926
S <- 1 / (3 * sin(beamAngle * pi / 180)) # 0.7887339
CS <- 1 / cos(30 * pi / 180) / sin(beamAngle * pi / 180) / 2 # 1.366127 (30deg from 3-beam pattern)
C <- 1 / (3 * cos(beamAngle * pi / 180)) # 0.3677926
# FIXME: check up and down; also read it and check
res@metadata$transformationMatrix <- matrix(
c(
2 * S, -S, -S,
0, -CS, CS,
C, C, C
),
nrow = 3, byrow = TRUE
)
} else {
S <- 1 / (3 * sin(beamAngle * pi / 180)) # 0.7887339
CS <- 1 / cos(30 * pi / 180) / sin(beamAngle * pi / 180) / 2 # 1.366127 (30deg from 3-beam pattern)
C <- 1 / (3 * cos(beamAngle * pi / 180)) # 0.3677926
# warning("*****FIXME: check up and down; also read it and check*****")
res@metadata$transformationMatrix <- matrix(
c(
2 * S, -S, -S,
0, CS, -CS,
-C, -C, -C
),
nrow = 3, byrow = TRUE
)
}
# For later use, RC says that the PC-ADP uses
# T = 2.576 -1.288 -1.288
# 0.000 -2.230 2.230
# 0.345 0.345 0.345
# and these are by the same formulae, with 25 switched to 15 (different beamAngle)
} else {
stop("can only handle 3-beam devices")
}
res@metadata$units <- list(v = "m/s", distance = "m")
if (missing(processingLog)) {
pl <- processingLogItem(paste("read.adp.sontek(\"", filename, "\", from=", from, ", to=", to, ")", sep = ""))
} else {
pl <- processingLogAppend(
processingLogItem(processingLog),
paste("read.adp.sontek(\"", filename, "\", from=", from, ", to=", to, ")", sep = "")
)
}
res@processingLog <- pl
res
}
sontek.time <- function(t, tz = getOption("oceTz")) {
minute <- bcdToInteger(t[1])
second <- bcdToInteger(t[2])
day <- bcdToInteger(t[3])
hour <- bcdToInteger(t[4])
year <- bcdToInteger(t[5])
year <- year + if (year >= 90) 1900 else 2000 # page 51 of System Integrator Guide
month <- bcdToInteger(t[6])
milliseconds <- readBin(t[7:8], "integer", n = 1, size = 2, endian = "little", signed = FALSE)
ISOdatetime(year, month, day, hour, minute, second + milliseconds / 1000, tz = tz)
}
#' Read an adp File in Serial Sontek Format
#'
#' Read a Sontek acoustic-Doppler profiler file, in a serial form that
#' is possibly unique to Dalhousie University.
#'
#' @param beamAngle angle between instrument axis and beams, in degrees.
#'
#' @param type a character string indicating the type of instrument.
#'
#' @param orientation optional character string specifying the orientation of the
#' sensor, provided for those cases in which it cannot be inferred from the
#' data file. The valid choices are `"upward"`, `"downward"`, and
#' `"sideward"`.
#'
#' @template adpTemplate
#'
#' @template encodingIgnoredTemplate
#'
#' @template adReadingMethodTemplate
#'
#' @author Dan Kelley and Clark Richards
#'
#' @family things related to adp data
#' @family functions that read adp data
read.adp.sontek.serial <- function(
file,
from = 1, to, by = 1, tz = getOption("oceTz"), longitude = NA, latitude = NA, type = c("adp", "pcadp"),
beamAngle = 25, orientation, encoding = NA, monitor = FALSE, processingLog,
debug = getOption("oceDebug"), ...) {
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, "\"")
}
}
# Data format is described in
# SonTek/YSI
# ADPManual_v710.pdf
# A3. Profile Header/CTD/GPS/Bottom Track,/SonWave/Profile Data Structures
bisectAdpSontekSerial <- function(buf, t.find, add = 0, tz = "UTC", debug = 0) {
oceDebug(debug, "bisectAdpSontekSerial(t.find=", format(t.find), ", add=", add, "\n")
len <- length(p)
lower <- 1
upper <- len
passes <- floor(10 + log(len, 2)) # won't need this many; only do this to catch coding errors
for (pass in 1:passes) {
middle <- floor((upper + lower) / 2)
year <- readBin(buf[p[middle] + 18:19], what = "integer", n = 1, size = 2, signed = FALSE, endian = "little")
day <- readBin(buf[p[middle] + 20], what = "integer", n = 1, size = 1, signed = FALSE)
month <- readBin(buf[p[middle] + 21], what = "integer", n = 1, size = 1, signed = FALSE)
min <- readBin(buf[p[middle] + 22], what = "integer", n = 1, size = 1, signed = FALSE)
hour <- readBin(buf[p[middle] + 23], what = "integer", n = 1, size = 1, signed = FALSE)
sec100 <- readBin(buf[p[middle] + 24], what = "integer", n = 1, size = 1, signed = FALSE)
sec <- readBin(buf[p[middle] + 25], what = "integer", n = 1, size = 1, signed = FALSE)
t <- ISOdatetime(year = year, month = month, day = day, hour = hour, min = min, sec = sec + sec100 / 100, tz = tz)
oceDebug(debug, "t=", format(t), " y=", year, " m=", month, " d=", format(day, width = 2),
" h=", format(hour, width = 2),
" m=", format(min, width = 2),
" s=", format(sec + sec100 / 100, width = 4),
"| pass", format(pass, width = 2), "/", passes, "| middle=", middle,
"(", format(middle / upper * 100, digits = 4), "%)\n",
sep = ""
)
if (t.find < t) {
upper <- middle
} else {
lower <- middle
}
if (upper - lower < 2) {
break
}
}
middle <- middle + add # may use add to extend before and after window
if (middle < 1) {
middle <- 1
}
if (middle > len) {
middle <- len
}
t <- ISOdatetime(readBin(buf[p[middle] + 18:19], "integer", size = 2, signed = FALSE, endian = "little"),
as.integer(buf[p[middle] + 21]), # month
as.integer(buf[p[middle] + 20]), # day
as.integer(buf[p[middle] + 23]), # hour
as.integer(buf[p[middle] + 22]), # min
as.integer(buf[p[middle] + 25]) + 0.01 * as.integer(buf[p[middle] + 24]),
tz = tz
)
oceDebug(debug, "result: t=", format(t), " at d[", middle, "]=", p[middle], "\n")
return(list(index = middle, time = t))
}
oceDebug(debug, paste("read.adp.sontek.serial(file[1]=\"", file[1],
"\", from=", from,
if (missing(to)) "to," else sprintf(", to=%s, ", format(to)),
", by=", by,
", latitude=", latitude, ", longitude=", longitude,
", monitor=", monitor,
", processingLog=(not shown), debug=", debug, ") START\n",
sep = ""
), unindent = 1)
nfile <- length(file)
if (nfile > 1) {
# handle multiple files
oceDebug(debug, "handling multiple files\n")
buf <- NULL
for (i in 1:nfile) {
if (monitor) {
cat("\"", file[i], "\" ", sep = "")
}
thisFile <- file(file[i], "rb")
seek(thisFile, 0, "end", rw = "read")
fileSize <- seek(thisFile, 0, origin = "start", rw = "read")
if (monitor) {
cat(fileSize, "bytes\n")
}
buf <- c(buf, readBin(thisFile, what = "raw", n = fileSize, endian = "little"))
close(thisFile)
}
filename <- paste("(\"", file[i], "\", ...)", sep = "")
} else {
# handle single file (which might be a connection, etc)
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)) {
filename <- "(connection)"
open(file, "rb")
on.exit(close(file))
}
# read whole file into buffer
seek(file, 0, "end", rw = "read")
fileSize <- seek(file, 0, origin = "start", rw = "read")
oceDebug(debug, "filesize=", fileSize, "\n")
buf <- readBin(file, what = "raw", n = fileSize, endian = "little")
}
# p <- .Call("ldc_sontek_adp", buf, 0, 0, 0, 0, -1) # no ctd, no gps, no bottom-track; all data
p <- do_ldc_sontek_adp(buf, 0, 0, 0, 0, -1) # no ctd, no gps, no bottom-track; all data
# read some unchanging things from the first profile only
serialNumber <- paste(readBin(buf[p[1] + 4:13], "character", n = 10, size = 1), collapse = "")
numberOfBeams <- readBin(buf[p[1] + 26], "integer", n = 1, size = 1, signed = FALSE)
if (missing(orientation)) {
orientation <- readBin(buf[p[1] + 27], "integer", n = 1, size = 1, signed = FALSE)
if (orientation == 0) {
orientation <- "upward"
} else if (orientation == 1) {
orientation <- "downward"
} else if (orientation == 2) {
orientation <- "sideward"
} else {
stop("orientation=", orientation, "but must be 0 (upward), 1 (downward), or 2 (sideward)")
}
} else {
if (orientation != "upward" && orientation != "downward") {
stop("orientation \"", orientation, " is not allowed; try \"upward\" or \"downward\"")
}
}
# 28 is tempmode
# coord.system 0=beam 1=XYZ 2=ENU
originalCoordinate <- readBin(buf[p[1] + 29], "integer", n = 1, size = 1, signed = FALSE)
if (originalCoordinate == 0) {
originalCoordinate <- "beam"
} else if (originalCoordinate == 1) {
originalCoordinate <- "xyz"
} else if (originalCoordinate == 2) {
originalCoordinate <- "enu"
} else {
stop("originalCoordinate=", originalCoordinate, "but must be 0 (beam), 1 (xyz), or 2 (enu)")
}
numberOfCells <- readBin(buf[p[1] + 30:31], "integer", n = 1, size = 2, signed = FALSE, endian = "little")
cellSize <- 0.01 * readBin(buf[p[1] + 32:33], what = "integer", n = 1, size = 2, signed = FALSE, endian = "little")
blankingDistance <- 0.01 * readBin(buf[p[1] + 34:35], what = "integer", n = 1, size = 2, signed = FALSE)
distance <- blankingDistance + cellSize * seq(from = 0.5, by = cellSize, length.out = numberOfCells)
# trim, if from and to are integers
if (!missing(to)) {
if (is.numeric(from) && is.numeric(to) && is.numeric(by)) {
if (from < 1) {
stop("from=", from, " but must exceed 1")
}
if (to < 1) {
stop("to=", to, " but must exceed 1")
}
if (by < 1) {
stop("by=", by, " but must exceed 1")
}
if (to <= from) {
stop("from=", from, " must exceed to=", to)
}
p <- p[seq(from = from, to = to, by = by)]
} else {
if (inherits(from, "POSIXt")) {
if (!inherits(to, "POSIXt")) {
stop("if 'from' is POSIXt, then 'to' must be, also")
}
if (!is.numeric(by)) {
warning("'by' must be numeric")
by <- 1
}
fromPair <- bisectAdpSontekSerial(buf, from, add = -1, tz = tz, debug = debug - 1)
from <- fromIndex <- fromPair$index
toPair <- bisectAdpSontekSerial(buf, to, add = 1, tz = tz, debug = debug - 1)
to <- to.index <- toPair$index
oceDebug(
debug, "from=", format(fromPair$t), " yields p[", fromIndex, "]\n",
" to =", format(toPair$t), "yields p[", to.index, "]\n",
" by=", by, "(not yet decoded)\n",
vectorShow(p, "p:"),
" p[", fromPair$index, "]=", p[fromPair$index], "at time", format(fromPair$t), "\n",
" p[", toPair$index, "]=", p[toPair$index], "at time", format(toPair$t), "\n"
)
p <- p[seq(from = from, to = to, by = by)]
}
}
}
np <- length(p)
pp <- sort(c(p, p + 1)) # for 2-byte addressing ("int" in the Sontek docs)
# pppp <- sort(c(p, p+1, p+2, p+3)) # for 4-byte addressing ("long" in the Sontek docs)
# read profile-specific things profile by profile
# profile.number <- readBin(buf[pppp+14], "integer", n=np, size=4)
# FIXME: should check that profile number is monotonic ... it may
# help us with daily blank-outs, also!
year <- readBin(buf[pp + 18], "integer", n = np, size = 2, signed = FALSE)
day <- readBin(buf[p + 20], "integer", n = np, size = 1, signed = FALSE)
month <- readBin(buf[p + 21], "integer", n = np, size = 1, signed = FALSE)
min <- readBin(buf[p + 22], "integer", n = np, size = 1, signed = FALSE)
hour <- readBin(buf[p + 23], "integer", n = np, size = 1, signed = FALSE)
sec100 <- readBin(buf[p + 24], "integer", n = np, size = 1, signed = FALSE)
sec <- readBin(buf[p + 25], "integer", n = np, size = 1, signed = FALSE)
time <- ISOdatetime(year, month, day, hour, min, sec + 0.01 * sec100, tz = tz)
rm(year, day, month, min, hour, sec100, sec) # possibly this space will come in handy
heading <- 0.1 * readBin(buf[pp + 40], "integer", n = np, size = 2, signed = TRUE)
pitch <- 0.1 * readBin(buf[pp + 42], "integer", n = np, size = 2, signed = TRUE)
roll <- 0.1 * readBin(buf[pp + 44], "integer", n = np, size = 2, signed = TRUE)
temperature <- 0.01 * readBin(buf[pp + 46], "integer", n = np, size = 2, signed = TRUE)
v <- array(numeric(), dim = c(np, numberOfCells, numberOfBeams))
a <- array(raw(), dim = c(np, numberOfCells, numberOfBeams))
q <- array(raw(), dim = c(np, numberOfCells, numberOfBeams))
ndata <- numberOfCells * numberOfBeams
i1 <- seq(1, ndata)
i2 <- seq(1, 2 * ndata)
for (ip in 1:np) {
p0 <- p[ip] + 79
# nolint start object_usage_linter
v_ <- matrix(0.001 * readBin(buf[p0 + i2], "integer", endian = "little", n = ndata, size = 2, signed = TRUE),
ncol = numberOfBeams, byrow = FALSE
)
p0 <- p0 + 2 * ndata
# NOTE: q is std-dev; need to multiply by 0.001 to get in m/s
q_ <- matrix(buf[p0 + i1], ncol = numberOfBeams, byrow = FALSE)
p0 <- p0 + ndata
a_ <- matrix(buf[p0 + i1], ncol = numberOfBeams, byrow = FALSE)
# nolint end object_usage_linter
for (b in seq_len(numberOfBeams)) {
v[ip, , b] <- v[, b]
a[ip, , b] <- a[, b]
q[ip, , b] <- q[, b]
}
if (monitor) {
if (ip %% 50) {
cat(".")
} else {
cat(".", ip, "\n")
}
}
}
if (monitor) {
cat("\nRead", np, "of the", np, "profiles in", filename[1], "\n")
}
S <- sin(beamAngle * pi / 180)
C <- cos(beamAngle * pi / 180)
# FIXME: use the transformation.matrix, if it has been discovered in a header
if (orientation == "upward") {
# OAR explains the method of determining the matrix.
transformationMatrix <- rbind(
c(2 / 3 / S, -1 / 3 / S, -1 / 3 / S),
c(0, -1 / sqrt(3) / S, 1 / sqrt(3) / S),
c(1 / 3 / C, 1 / 3 / C, 1 / 3 / C)
)
} else {
transformationMatrix <- rbind(
c(2 / 3 / S, -1 / 3 / S, -1 / 3 / S),
c(0, 1 / sqrt(3) / S, -1 / sqrt(3) / S),
c(1 / 3 / C, -1 / 3 / C, -1 / 3 / C)
)
}
res <- new("adp")
res@metadata$manufacturer <- "sontek"
res@metadata$instrumentType <- "adp"
res@metadata$serialNumber <- serialNumber
res@metadata$filename <- filename
res@metadata$latitude <- latitude
res@metadata$longitude <- longitude
res@metadata$transformationMatrix <- transformationMatrix
res@metadata$measurementStart <- 0 # FIXME: should fill in
res@metadata$measurementEnd <- np # FIXME: should fill in
res@metadata$measurementDeltat <- mean(diff(as.numeric(time)))
res@metadata$subsampleStart <- 0 # FIXME: should fill in
res@metadata$subsampleEnd <- np
res@metadata$subsampleDeltat <- mean(diff(as.numeric(time)))
res@metadata$frequency <- NA # FIXME
res@metadata$numberOfSamples <- np
res@metadata$numberOfBeams <- numberOfBeams
res@metadata$originalCoordinate <- originalCoordinate
res@metadata$oceCoordinate <- originalCoordinate
res@metadata$beamAngle <- beamAngle
res@metadata$oceBeamUnspreaded <- FALSE
res@metadata$orientation <- orientation
res@metadata$units <- list(v = "m/s", distance = "m")
res@data <- list(
v = v, a = a, q = q,
distance = distance,
time = time,
heading = heading, pitch = pitch, roll = roll,
temperature = temperature,
pressure = rep(0, length(temperature)),
distance = distance
)
if (missing(processingLog)) {
processingLog <- paste(deparse(match.call()), sep = "", collapse = "")
}
hitem <- processingLogItem(processingLog)
res@processingLog <- hitem
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.