Nothing
# vim:textwidth=80:expandtab:shiftwidth=4:softtabstop=4
testNewSerialNumberDecoder <- TRUE # see email exchange with CR dated 2023-10-05
# Data format overview
# hardware [a5 05 X1 X2] 48 bytes, 2*(short word made from X1 and X2)
# head [a5 04 X1 X2] 224 bytes, 2*(short word made from X1 and X2)
# user [a5 00 X1 X2] 512 bytes, 2*(short word made from X1 and X2)
# profiles, each starting with a5 2a [aquadoppHR] or a5 21 [aquadoopp profiler] [other]
# DOCUMENTATION BUGS
# 1. p38 System Integrator Guide says to offset 53 bytes for velocity, but I have to offset 54 to recover data
# that match the manufacturer's (*.v1, *.v2, *.v3) files.
# AQUADOPP
# notes for nortek:
# 1. "spare" at offset 74 (page 31) now seems to be salinity
# 2. extra byte
# 3. should state the order of headers at the start, not end
# 4. should state the algorithms to infer cellSize, blankingDistance, etc. from file
# 5. beamAngle should be in data file
# 6. generally, docs should indicate everything that is in the files, e.g. (prominently!)
# the beamAngles in the 'head' configuration section.
# 7. the C code suggests the velocityScale is in the second bit of conf.hMode
# but the docs suggest the fifth bit (page 31)
#' Decode a Nortek Header
#'
#' Decode data in a Nortek ADV or ADP header.
#'
#' Decodes the header in a binary-format Nortek ADV/ADP file. This function is
#' designed to be used by [read.adp()] and [read.adv()],
#' but can be used directly as well. The code is based on information in the
#' Nortek System Integrator Guide (2008) and on postings on the Nortek
#' ``knowledge center'' discussion board. One might assume that the latter is
#' less authoritative than the former. For example, the inference of cell size
#' follows advice found at
#' https://www.nortekusa.com/en/knowledge-center/forum/hr-profilers/736804717,
#' which contains a typo in an early posting that is
#' corrected later on.
#'
#' @param buf a ``raw'' buffer containing the header
#'
#' @param type type of device
#'
#' @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 \dots additional arguments, passed to called routines.
#'
#' @return A list containing elements `hardware`, `head`, `user`
#' and `offset`. The easiest way to find the contents of these is to run
#' this function with `debug=3`.
#'
#' @author Dan Kelley and Clark Richards
#'
#' @seealso Most users should employ the functions [read.adp()] and
#' [read.adv()] instead of this one.
#'
#' @references 1. Information on Nortek profilers (including the System
#' Integrator Guide, which explains the data format byte-by-byte) is available
#' at `https://www.nortekusa.com/usa?set_language=usa` after login.
#'
#' 2. The Nortek Knowledge Center
#' https://www.nortekusa.com/en/knowledge-center
#' may be of help if problems arise in dealing with data from Nortek instruments.
#'
#' 3. Nortek,
#' "Classic Integrators Guide: Aquadopp | Aquadopp DW | Aquadopp Profiler | HQ Aquadopp Profiler | Vector | AWAC."
#' Nortek AS, 2022.
decodeHeaderNortek <- function(
buf,
type = c("aquadoppHR", "aquadoppProfiler", "aquadopp", "aquadoppPlusMagnetometer", "vector"),
debug = getOption("oceDebug"),
...) {
type <- match.arg(type)
oceDebug(debug, "decodeHeaderNortek(buf, type=\"", type, "\", ...) {\n", sep = "", unindent = 1)
oceDebug(debug, "buf starts:", buf[1:20], "\n")
degToRad <- atan2(1, 1) / 45
syncCode <- as.raw(0xa5)
idHardwareConfiguration <- as.raw(0x05)
idHeadConfiguration <- as.raw(0x04)
idUserConfiguration <- as.raw(0x00)
headerLengthHardware <- 48
headerLengthHead <- 224
headerLengthUser <- 512
hardware <- head <- user <- list()
o <- 0 # offset
for (header in 1:3) {
# FIXME: code is needlessly written as if headers could be in different order
oceDebug(debug, "\n")
# oceDebug(debug, "examining buf[o+2]=", buf[o+2], "to see what type of header block is next...\n")
if (buf[o + 1] != syncCode) {
stop("expecting syncCode 0x", syncCode, " but got 0x", buf[o + 1], " instead (while reading header #", header, ")")
}
if (buf[o + 2] == idHardwareConfiguration) {
# see page 29 of System Integrator Guide
oceDebug(debug, "\nHARDWARE CONFIGURATION\n", unindent = 1)
hardware$size <- readBin(buf[o + 3:4], "integer", signed = FALSE, n = 1, size = 2, endian = "little")
if (hardware$size != 24) {
stop("size of hardware header expected to be 24 two-byte words, but got ", hardware$size)
}
if (2 * hardware$size != headerLengthHardware) {
stop("size of hardware header expected to be ", headerLengthHardware, "but got ", hardware$size)
}
oceDebug(debug, "hardware$size=", hardware$size, "\n")
# 2023-10-05: problem reading a CR file. The serial number has an
# 0x80 in it (which is a control character in some encodings).
# Let's just trim anything that is not in the ASCII table.
if (testNewSerialNumberDecoder) {
tmp <- buf[o + 5:18]
# Remove any non-printing char, and anything following it
if (any(tmp >= 0x80)) {
w <- which(tmp >= 0x80)[1]
tmp <- if (w > 1L) tmp[seq(1L, w - 1L)] else ""
warning("serial number contains control characters, so only preceeding characters retained\n")
}
hardware$serialNumber <- gsub(" *$", "", paste(readBin(tmp, "character", n = nchar(tmp), size = 1), collapse = ""))
} else {
hardware$serialNumber <- gsub(" *$", "", paste(readBin(buf[o + 5:18], "character", n = 14, size = 1), collapse = ""))
}
oceDebug(debug, "hardware$serialNumber: ", hardware$serialNumber, "\n")
hardware$config <- readBin(buf[o + 19:20], "raw", n = 2, size = 1)
oceDebug(debug, "hardware$config:", hardware$config, "\n")
hardware$frequency <- readBin(buf[o + 21:22], "integer", n = 1, size = 2, endian = "little", signed = FALSE) # not used
oceDebug(debug, "hardware$frequency:", hardware$frequency, "\n")
hardware$picVersion <- readBin(buf[o + 23:24], "integer", n = 1, size = 2, endian = "little")
oceDebug(debug, "hardware$picVersion=", hardware$picVersion, "\n")
hardware$hwRevision <- readBin(buf[o + 25:26], "integer", n = 1, size = 2, endian = "little")
oceDebug(debug, "hardware$hwRevision=", hardware$hwRevision, "\n")
hardware$recSize <- readBin(buf[o + 27:28], "integer", n = 1, size = 2, signed = FALSE, endian = "little")
oceDebug(debug, "hardware$recSize=", hardware$recSize, "\n")
hardware$velocityRange <- readBin(buf[o + 29:30], "integer", n = 1, size = 2, signed = FALSE, endian = "little")
oceDebug(debug, "hardware$velocityRange=", hardware$velocityRange, "\n")
hardware$fwVersion <- as.numeric(paste(readBin(buf[o + 43:46], "character", n = 4, size = 1), collapse = ""))
oceDebug(debug, "hardware$fw.version=", hardware$fw.version, "\n")
o <- o + 2 * hardware$size
} else if (buf[o + 2] == idHeadConfiguration) {
# see page 30 of System Integrator Guide
oceDebug(debug, "HEAD CONFIGURATION\n", unindent = 1)
# buf <- readBin(file, "raw", headerLengthHead)
head$size <- readBin(buf[o + 3:4], "integer", signed = FALSE, n = 1, size = 2)
if (2 * head$size != headerLengthHead) {
stop("size of head header expected to be ", headerLengthHead, "but got ", head$size)
}
oceDebug(debug, "head$size=", head$size, "\n")
# Nortek doc "system-integrator-manual_Mar2016.pdf" (page 23) says for the "head configuration":
# bit 0: Pressure sensor (0=no, 1=yes)
# bit 1: Magnetometer sensor (0=no, 1=yes)
# bit 2: Tilt sensor (0=no, 1=yes)
# bit 3: Tilt sensor mounting (0=up, 1=down)
# issue1220 debug <- debug + 10
# issue1220 head$config <- byteToBinary(buf[o+5:6], endian="little")
# issue1220 oceDebug(debug, "head$config=", head$config, "\n")
# issue1220 configTEST <- rawToBits(buf[o+5:6])
# issue1220 oceDebug(debug, "configTEST=", configTEST, "\n")
# issue1220 oceDebug(debug, "head$config=", head$config, "\n")
# issue1220 head$configPressureSensor <- substr(head$config[1], 1, 1) == "1"
# issue1220 oceDebug(debug, "head$configPressureSensor=", head$configPressureSensor, "\n")
# issue1220 head$configMagnetometerSensor <- substr(head$config[1], 2, 2) == "1"
# issue1220 oceDebug(debug, "head$configMagnetometerSensor=", head$configMagnetometerSensor, "\n")
# issue1220 head$configTiltSensor <- substr(head$config[1], 3, 3) == "1"
# issue1220 oceDebug(debug, "head$configTiltSensor=", head$configTiltSensor, "\n")
# issue1220 head$tiltSensorOrientation <- if (substr(head$config[1], 4, 4) == "1") "downward" else "upward"
# issue1220 oceDebug(debug, "head$tiltSensorOrientation=", head$tiltSensorOrientation, "\n")
tmpBits <- rawToBits(buf[o + 5])
head$configPressureSensor <- tmpBits[1] == as.raw(0x1)
head$configMagnetometerSensor <- tmpBits[2] == as.raw(0x1)
head$configTiltSensor <- tmpBits[3] == as.raw(0x1)
head$configTiltSensorOrientation <- ifelse(tmpBits[4] == as.raw(0x1), "downward", "upward")
oceDebug(debug, "head$configPressureSensor=", head$configPressureSensor, "\n")
oceDebug(debug, "head$configMagnetometerSensor=", head$configMagnetometerSensor, "\n")
oceDebug(debug, "head$configTiltSensor=", head$configTiltSensor, "\n")
oceDebug(debug, "head$configTiltSensorOrientation=", head$configTiltSensorOrientation, "\n")
head$frequency <- readBin(buf[o + 7:8], "integer", n = 1, size = 2, endian = "little", signed = FALSE)
oceDebug(debug, "head$frequency=", head$frequency, "kHz\n")
head$headType <- readBin(buf[o + 9:10], "integer", n = 1, size = 2, endian = "little")
oceDebug(debug, "head$headType=", head$headType, "\n")
head$headSerialNumber <- gsub(" *$", "", paste(readBin(buf[o + 11:22], "character", n = 12, size = 1), collapse = ""))
oceDebug(debug, "head$headSerialNumber=", head$headSerialNumber, "\n")
# NOTE: p30 of System Integrator Guide does not detail anything from offsets 23 to 119;
# the inference of beamAngles and transformationMatrix is drawn from other code.
# Since I don't trust any of this, I hard-wire beamAngle in at the end.
head$beamAngles <- readBin(buf[o + 23:30], "integer", n = 4, size = 2, endian = "little", signed = TRUE)
oceDebug(debug, "head$beamAngles=", paste(head$beamAngles, collapse = ", "), " deg\n")
# Transformation matrix (before division by 4096)
# FIXME: should we change the sign of rows 2 and 3 if pointed down??
head$transformationMatrix <- matrix(readBin(buf[o + 31:48], "integer", n = 9, size = 2, endian = "little"), nrow = 3, byrow = TRUE) / 4096
oceDebug(debug, "head$transformationMatrix\n")
oceDebug(debug, format(head$transformationMatrix[1, ], width = 15), "\n")
oceDebug(debug, format(head$transformationMatrix[2, ], width = 15), "\n")
oceDebug(debug, format(head$transformationMatrix[3, ], width = 15), "\n")
head$numberOfBeams <- readBin(buf[o + 221:222], "integer", n = 1, size = 2, endian = "little")
oceDebug(debug, "head$numberOfBeams=", head$numberOfBeams, "\n")
o <- o + 2 * head$size
} else if (buf[o + 2] == idUserConfiguration) {
# User Configuration [p30-32 of System Integrator Guide]
oceDebug(debug, "USER CONFIGURATION\n", unindent = 1)
user$size <- readBin(buf[o + 3:4], what = "integer", n = 1, size = 2, endian = "little")
if (2 * user$size != headerLengthUser) {
stop("expected header length is ", headerLengthUser, "but got ", user$size)
}
# buf <- readBin(file, "raw", headerLengthUser)
# Keep both the 'T' names and also oce-type names. The latter
# have names taken from reference 3, but these names are not
# literal, since we know from personal communications with Nortek
# that the physical quantities are computed via equations using T1,
# T2 etc. Reference 3 (page 29) describes the T variables as
# follows.
oceDebug(
debug, "T* variable meanings from Sec 5.1, page 29 of reference 3:\n",
" T1: transmit pulse length (counts)\n",
" T2: blanking distance (counts)\n",
" T3: receive length (counts)\n",
" T4: time between pings (counts)\n",
" T5: time between burst sequences (counts)\n"
)
user$T1 <- readBin(buf[o + 5:6], "integer", n = 1, size = 2, endian = "little", signed = FALSE)
oceDebug(debug, "user$T1=", user$T1, " (called System1 in .hdr files)\n")
user$T2 <- readBin(buf[o + 7:8], "integer", n = 1, size = 2, endian = "little", signed = FALSE)
oceDebug(debug, "user$T2=", user$T2, " (called System2 in .hdr files)\n")
user$T3 <- readBin(buf[o + 9:10], "integer", n = 1, size = 2, endian = "little", signed = FALSE)
oceDebug(debug, "user$T3=", user$T3, " (called System3 in .hdr files)\n")
user$T4 <- readBin(buf[o + 11:12], "integer", n = 1, size = 2, endian = "little", signed = FALSE)
oceDebug(debug, "user$T4=", user$T4, " (called System4 in .hdr files)\n")
user$T5 <- readBin(buf[o + 13:14], "integer", n = 1, size = 2, endian = "little", signed = FALSE)
oceDebug(debug, "user$T5=", user$T5, " (called System5 in .hdr files)\n")
user$transmitPulseLength <- user$T1
user$receiveLength <- user$T3
user$timeBetweenPings <- user$T4
oceDebug(debug, "user$timeBetweenPings=", user$timeBetweenPings, " counts\n")
# next is called NPings in Sec 5.1, page 29 of reference 3
user$numberOfBeamSequencesPerBurst <- readBin(buf[o + 15:16], "integer", n = 1, size = 2, endian = "little", signed = FALSE)
oceDebug(debug, "user$numberOfBeamSequencesPerBurst=", user$numberOfBeamSequencesPerBurst, " counts\n")
user$averagingInterval <- readBin(buf[o + 17:18], "integer", n = 1, size = 2, endian = "little")
oceDebug(debug, "user$averagingInterval=", user$averagingInterval, " seconds\n")
user$numberOfBeams <- readBin(buf[o + 19:20], "integer", n = 1, size = 2, endian = "little")
oceDebug(debug, "user$numberOfBeams=", user$numberOfBeams, "\n")
user$timCtrlReg <- readBin(buf[o + 21:22], "raw", n = 2)
oceDebug(debug, "user$timCtrlReg=", user$timCtrlReg, "\n")
# issue 2165 start new code
user$coordSystem <- readBin(buf[o + 33:34], "integer", signed = FALSE, n = 1, size = 2, endian = "little")
oceDebug(debug, "user$coordSystem=", user$coordSystem, " (issue 2165. 0=ENU, 1=XYZ, 2=BEAM?)\n")
user$nbins <- readBin(buf[o + 35:36], "integer", signed = FALSE, n = 1, size = 2, endian = "little")
oceDebug(debug, "user$nbins=", user$nbins, " (issue 2165: number of bins?)\n")
user$binLength <- readBin(buf[o + 37:38], "integer", signed = FALSE, n = 1, size = 2, endian = "little")
oceDebug(debug, "user$binLength=", user$binLength, " (issue 2165: cellSize??)\n")
# issue 2165 end new code
user$measurementInterval <- readBin(buf[o + 39:40], "integer", n = 1, size = 2, endian = "little")
oceDebug(debug, "user$measurementInterval=", user$measurementInterval, "\n")
user$deploymentName <- readBin(buf[o + 41:47], "character", n = 1, size = 6)
oceDebug(debug, "user$deploymentName=", user$deploymentName, "\n")
user$comments <- readBin(buf[o + 257 + 0:179], "character", n = 1, size = 180)
oceDebug(debug, "user$comments=", user$comments, "\n")
# Nortek doc "system-integrator-manual_Mar2016.pdf" (page 24-25) says for the "mode":
# bit 0: use user specified sound speed (0=no, 1=yes)
# bit 1: diagnostics/wave mode 0=disable, 1=enable)
# bit 2: analog output mode (0=disable, 1=enable)
# bit 3: output format (0=Vector, 1=ADV)
# bit 4: scaling (0=1 mm, 1=0.1 mm)
# bit 5: serial output (0=disable, 1=enable)
# bit 6: reserved EasyQ
# bit 7: stage (0=disable, 1=enable)
# bit 8: output power for analog input (0=disable, 1=enable)
# and we want bit 4. In R notation, that is rawToBits()[5]
#
# issue1220 user$mode <- byteToBinary(buf[o+59:60], endian="little")
# issue1220 user$velocityScale <- if (substr(user$mode[2], 4, 4) == "0") 0.001 else 0.0001
user$mode <- rawToBits(buf[o + 59:60])
oceDebug(debug, "user$mode: ", user$mode, "\n")
user$velocityScale <- ifelse(user$mode[5] == 0x00, 0.001, 0.0001)
oceDebug(debug, "user$velocityScale: ", user$velocityScale, "\n")
tmp.cs <- readBin(buf[o + 33:34], "integer", n = 1, size = 2, endian = "little")
if (tmp.cs == 0) {
user$originalCoordinate <- "enu"
} # page 31 of System Integrator Guide
else if (tmp.cs == 1) {
user$originalCoordinate <- "xyz"
} else if (tmp.cs == 2) {
user$originalCoordinate <- "beam"
} else {
stop("unknown originalCoordinate ", tmp.cs)
}
oceDebug(debug, "user$originalCoordinate: ", user$originalCoordinate, "\n")
user$numberOfCells <- readBin(buf[o + 35:36], "integer", n = 1, size = 2, endian = "little")
oceDebug(debug, "user$numberOfCells: ", user$numberOfCells, "\n")
user$hBinLength <- readBin(buf[o + 37:38], "integer", n = 1, size = 2, endian = "little", signed = FALSE)
oceDebug(debug, "user$hBinLength: ", user$hBinLength, " (p31 of System Integrator Guide)\n")
# The cell size is computed with different formulae for different devices, and
# different frequencies. See
# https://www.nortekusa.com/en/knowledge-center/forum/hr-profilers/736804717
# for the details (and note that there are some typos from the Nortek advisor, which
# get corrected later on that webpage).
if (type == "aquadoppHR") {
# (NOTE: 1674 is hBinLength)
# CS = (1674/256)*0.00675*cos(25) = 0.04 m, for a 2 MHz instrument
# For a 1 MHz instrument you must multiply by twice the sampling distance, i.e. 0.0135 m
if (isTRUE(all.equal.numeric(head$frequency, 2000))) {
oceDebug(debug, "computing user$cellSize and user$blankingDistance for type=\"aquadoppHR\" and frequency=2000\n")
user$cellSize <- user$hBinLength / 256 * 0.00675 * cos(25 * degToRad)
user$blankingDistance <- user$T2 * 0.00675 * cos(25 * degToRad) - user$cellSize
} else if (isTRUE(all.equal.numeric(head$frequency, 1000))) {
oceDebug(debug, "computing user$cellSize and user$blankingDistance for type=\"aquadoppHR\" and frequency=1000\n")
user$cellSize <- user$hBinLength / 256 * 0.01350 * cos(25 * degToRad)
user$blankingDistance <- user$T2 * 0.01350 * cos(25 * degToRad) - user$cellSize
} else {
oceDebug(debug, "computing user$cellSize and user$blankingDistance for type=\"aquadoppHR\" and frequency neither 1000 nor 2000\n")
warning("unknown frequency ", head$frequency, " (only understand 1MHz and 2MHz); using 2Mhz formula to calculate cell size")
user$cellSize <- user$hBinLength / 256 * 0.00675 * cos(25 * degToRad)
user$blankingDistance <- user$T2 * 0.00675 * cos(25 * degToRad) - user$cellSize
}
} else if (type == "aquadoppProfiler") {
# email from Evan Price (Nortek) dated 2023-10-30T12:05 (Halifax time)
# CS in mm:
# CS = (BinLength/256)*(750*17*15/(Frequency_in_kHz*4))*cosd(25) [mm]
# BD in mm:
# BD = T2*1000*(750/32768)*cosd(25) - CS[mm]
oceDebug(debug, "testing using Nortek formula from 2023-10-30T12:05 (Halifax time)\n")
oceDebug(debug, " user$binLength=", user$binLength, "\n")
oceDebug(debug, " user$hBinLength=", user$hBinLength, "\n")
cosAngle <- cos(25 * degToRad) # FIXME: is beam angle always 25deg for aquadopp?
CS <- (user$binLength / 256) * (750 * 17 * 15) / (head$frequency * 4) * cosAngle
oceDebug(debug, " CS = ", CS, "[mm], i.e. ", CS / 1000, "m\n")
BD <- user$T2 * 1000 * (750 / 32768) * cosAngle - CS
oceDebug(debug, " BD = ", BD, "[mm], i.e. ", BD / 1000, "m\n")
oceDebug(
debug, " NOTE: BD is 0.41 for a particular test file, but the .hdr file\n",
" for for a test file states 0.40. We learned from Nortek on 2023-11-07\n",
" that our computed value is correct, i.e. that the value in the .hdr\n",
" file is wrong\n"
)
user$cellSize <- CS / 1000
user$blankingDistance <- BD / 1000
#<issue 2165>if (isTRUE(all.equal.numeric(head$frequency, 2000))) {
#<issue 2165> oceDebug(debug, "computing user$cellSize and user$blankingDistance for type=\"aquadoppProfiler\" and frequency 2000\n")
#<issue 2165> user$cellSize <- user$hBinLength / 256 * 0.0239 * cos(25 * degToRad)
#<issue 2165>} else if (isTRUE(all.equal.numeric(head$frequency, 1000))) {
#<issue 2165> oceDebug(debug, "computing user$cellSize and user$blankingDistance for type=\"aquadoppProfiler\" and frequency 1000\n")
#<issue 2165> user$cellSize <- user$hBinLength / 256 * 0.0478 * cos(25 * degToRad)
#<issue 2165>} else if (isTRUE(all.equal.numeric(head$frequency, 600))) {
#<issue 2165> oceDebug(debug, "computing user$cellSize and user$blankingDistance for type=\"aquadoppProfiler\" and frequency 600\n")
#<issue 2165> user$cellSize <- user$hBinLength / 256 * 0.0797 * cos(25 * degToRad)
#<issue 2165>} else if (isTRUE(all.equal.numeric(head$frequency, 400))) {
#<issue 2165> oceDebug(debug, "computing user$cellSize and user$blankingDistance for type=\"aquadoppProfiler\" and frequency 400\n")
#<issue 2165> user$cellSize <- user$hBinLength / 256 * 0.1195 * cos(25 * degToRad)
#<issue 2165>} else {
#<issue 2165> warning("unknown frequency", head$frequency, "(only understand 1MHz and 2MHz); using 1Mhz formula to calculate cell size")
#<issue 2165> oceDebug(debug, "computing user$cellSize and user$blankingDistance for",
#<issue 2165> "type=\"aquadoppProfiler\" and frequency not 400, 600, 1000 or 2000\n")
#<issue 2165> user$cellSize <- user$hBinLength / 256 * 0.0478 * cos(25 * degToRad)
#<issue 2165> #user$cellSize <- user$hBinLength / 256 * 0.00675 * cos(25 * degToRad)
#<issue 2165>}
#<issue 2165>user$blankingDistance <- user$T2 * 0.0229 * cos(25 * degToRad) - user$cellSize
#<issue 2165>oceDebug(debug, "blankingDistance from old formula:", user$blankingDistance, " m\n")
# user$blankingDistanceTEST <- user$T2 * 100 * (750/32768) - user$T1 * 100 * ((750/(1000*head$frequency))*16)
# oceDebug(debug, "blankingDistance from new formula (not used):", user$blankingDistanceTEST, " m?\n")
} else if (type == "aquadopp" || type == "aquadoppPlusMagnetometer") {
oceDebug(debug, "computing user$cellSize and user$blankingDistance for type=\"aquadoppPlusMagnetometer\"\n")
# user$cellSize <- user$hBinLength / 256 * 0.00675 * cos(25 * degToRad)
# user$blankingDistance <- user$T2 * 0.00675 * cos(25 * degToRad) - user$cellSize
warning("using fixed cell size and blanking distance for Aquadopp, since cannot infer these from the file")
user$cellSize <- 0.75 # value for one particular test file
user$blankingDistance <- 0.37 # value for one particular test file
} else if (type == "vector") {
# Formula (revised) from Nortek https://www.nortekusa.com/en/knowledge-center/forum/hr-profilers/595666030
# Cell size (mm) = T3counts * 1000 * 750 / 480000
oceDebug(debug, "computing user$cellSize and user$blankingDistance for type=\"vector\"\n")
soundspeed <- 1500
user$cellSize <- user$receiveLength * (soundspeed / 2) / 480.0e3 # drop the 1000 to get m
user$blankingDistance <- 0 # ?
} else {
warning("unknown instrument type \"", type, "\", so calculating cell size as though it is a 2MHz AquadoppHR")
oceDebug(debug, "computing user$cellSize and user$blankingDistance for unknown type=\"", type, "\"\n", sep = "")
user$cellSize <- user$hBinLength / 256 * 0.00675 * cos(25 * degToRad)
user$blankingDistance <- user$T2 * 0.00675 * cos(25 * degToRad) - user$cellSize
}
# FIXME: eventually, incorporate the AWAC and Continental devices, using
# the above formulae, with coefficients as below (from the nortek forum)
#
# For the AWAC the table for cell size looks like this:
# 0.0478 for 1 MHz
# 0.0797 for 600 kHz
#
# And for the Continental:
# 0.0945 for 470 kHz
# 0.2221 for 190 kHz
oceDebug(debug, "cellSize=", user$cellSize, "m (FIXME: using formula from unknown source)\n")
user$measurementInterval <- readBin(buf[o + 39:40], "integer", n = 1, size = 2, endian = "little")
oceDebug(debug, "measurementInterval=", user$measurementInterval, "\n")
# bit 0: use user specified sound speed (0=no, 1=yes)
# bit 1: diagnostics/wave mode 0=disable, 1=enable)
# bit 2: analog output mode (0=disable, 1=enable)
# bit 3: output format (0=Vector, 1=ADV)
# bit 4: scaling (0=1 mm, 1=0.1 mm)
# bit 5: serial output (0=disable, 1=enable)
# bit 6: reserved EasyQ
# bit 7: stage (0=disable, 1=enable)
# bit 8: output power for analog input (0=disable, 1=enable)
# FIXME: Sample.cpp has 0.022888 for the factor on user$T2
# if (isTRUE(all.equal.numeric(head$frequency, 1000))) {
# #user$blankingDistance <- cos(25*degToRad) * (0.0135 * user$T2 - 12 * user$T1 / head$frequency)
# user$blankingDistance <- cos(25*degToRad) * (0.0135 * user$blankingDistance - 12 * user$transmitPulseLength / head$frequency)
# } else if (isTRUE(all.equal.numeric(head$frequency, 2000))) {
# #user$blankingDistance <- cos(25*degToRad) * (0.00675 * user$T2 - 12 * user$T1 / head$frequency)
# user$blankingDistance <- cos(25*degToRad) * (0.00675 * user$blankingDistance - 12 * user$transmitPulseLength / head$frequency)
# } else {
# user$blankingDistance <- 0
# }
# cat("adp.nortek.R:245 user$blankingDistance", user$blankingDistance, "\n")
oceDebug(debug, "blankingDistance=", user$blankingDistance, "; user$T1=", user$T1, "; user$T2=", user$T2, "\n")
user$swVersion <- readBin(buf[o + 73:74], "integer", n = 1, size = 2, endian = "little") / 10000
oceDebug(debug, "swVersion=", user$swVersion, "\n")
user$salinity <- readBin(buf[o + 75:76], "integer", n = 1, size = 2, endian = "little") * 0.1
oceDebug(debug, "salinity=", user$salinity, "\n")
o <- o + 2 * user$size
} else {
stop(
"cannot understand byte 0x", buf[o + 1], "; expecting one of the following: 0x",
idHardwareConfiguration, " [hardware configuration] 0x",
idHeadConfiguration, " [head configuration] or 0x", idUserConfiguration, " [user configuration]\n"
)
}
}
list(hardware = hardware, head = head, user = user, offset = o + 1)
}
#' Read an adp File in Nortek Aquadopp Format
#'
#' The R code is based on information in the Nortek System Integrator Guide
#' (2017), postings on the Nortek ``knowledge center'' discussion board, and
#' discussions with Nortek engineers (Dec. 2020).
#'
#' @param type Either "aquadopp" for a standard aquadopp file (the default), or
#' "aquadoppPlusMagnetometer" for a file which includes the raw magnetometer
#' data.
#'
#' @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"`.
#'
#' @param distance Optional vector holding the distances of bin centres from the
#' sensor. This argument is ignored except for Nortek profilers, and need not
#' be given if the function determines the distances correctly from the data.
#' The problem is that the distance is poorly documented in the Nortek System
#' Integrator Guide (2008 edition, page 31), so the function must rely on
#' word-of-mouth formulae that do not work in all cases.
#'
#' @param despike if `TRUE`, [despike()] will be used to clean
#' anomalous spikes in heading, etc.
#'
#' @references
#' 1. Information on Nortek profilers (including the System Integrator Guide,
#' which explains the data format byte-by-byte) is available at
#' `https://www.nortekusa.com/`. (One must join the site to see the
#' manuals.)
#'
#' 2. The Nortek Knowledge Center
#' `https://www.nortekusa.com/en/knowledge-center` may be of help if
#' problems arise in dealing with data from Nortek instruments.
#'
#' @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.aquadopp <- function(
file,
from = 1,
to,
by = 1,
tz = getOption("oceTz"),
longitude = NA,
latitude = NA,
type = "aquadopp",
orientation,
distance,
monitor = FALSE,
despike = FALSE,
encoding = NA,
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
}
return(read.adp.nortek(file,
from = from, to = to, by = by, tz = tz,
longitude = longitude, latitude = latitude,
type = type, orientation = orientation, distance = distance,
monitor = monitor, despike = despike, processingLog = processingLog,
debug = debug, ...
))
}
#' Read Nortek Aquadopp-HR File
#'
#' The R code is based on information in
#' the Nortek System Integrator Guide (2008) and on postings on the Nortek
#' ``knowledge center'' discussion board. One might assume that the latter is
#' less authoritative than the former. For example, the inference of cell size
#' follows advice found at
#' `https://www.nortekusa.com/en/knowledge-center/forum/hr-profilers/736804717`,
#' which contains a typo in an early posting that is
#' corrected later on.
#'
#' @param despike if `TRUE`, [despike()] will be used to clean
#' anomalous spikes in heading, etc.
#'
#' @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"`.
#'
#' @param distance Optional vector holding the distances of bin centres from the
#' sensor. This argument is ignored except for Nortek profilers, and need not
#' be given if the function determines the distances correctly from the data.
#' The problem is that the distance is poorly documented in the Nortek System
#' Integrator Guide (2008 edition, page 31), so the function must rely on
#' word-of-mouth formulae that do not work in all cases.
#'
#' @references
#' 1. Information on Nortek profilers (including the System Integrator Guide,
#' which explains the data format byte-by-byte) is available at
#' `https://www.nortekusa.com/`. (One must join the site to see the
#' manuals.)
#'
#' 2. The Nortek Knowledge Center
#' `https://www.nortekusa.com/en/knowledge-center` may be of help if
#' problems arise in dealing with data from Nortek instruments.
#'
#' @template adpTemplate
#'
#' @template encodingIgnoredTemplate
#'
#' @template adReadingMethodTemplate
#'
#' @author Dan Kelley
#'
#' @family things related to adp data
#' @family functions that read adp data
read.aquadoppHR <- function(
file,
from = 1,
to,
by = 1,
tz = getOption("oceTz"),
longitude = NA,
latitude = NA,
orientation = orientation,
distance,
monitor = FALSE,
despike = FALSE,
encoding = NA,
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
}
return(read.adp.nortek(file,
from = from, to = to, by = by, tz = tz,
longitude = longitude, latitude = latitude,
type = "aquadoppHR",
orientation = orientation, distance = distance,
monitor = monitor, despike = despike, processingLog = processingLog,
debug = debug, ...
))
}
#' Read an adp File in Nortek Aquadopp Format
#'
#' The R code is based on information in
#' the Nortek System Integrator Guide (2008) and on postings on the Nortek
#' ``knowledge center'' discussion board. One might assume that the latter is
#' less authoritative than the former. For example, the inference of cell size
#' follows advice found at
#' `https://www.nortekusa.com/en/knowledge-center/forum/hr-profilers/736804717`,
#' which contains a typo in an early posting that is
#' corrected later on.
#'
#' @param despike if `TRUE`, [despike()] will be used to clean
#' anomalous spikes in heading, etc.
#'
#' @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"`.
#'
#' @param distance Optional vector holding the distances of bin centres from the
#' sensor. This argument is ignored except for Nortek profilers, and need not
#' be given if the function determines the distances correctly from the data.
#' The problem is that the distance is poorly documented in the Nortek System
#' Integrator Guide (2008 edition, page 31), so the function must rely on
#' word-of-mouth formulae that do not work in all cases.
#'
#' @references
#' 1. Information on Nortek profilers (including the System Integrator Guide,
#' which explains the data format byte-by-byte) is available at
#' `https://www.nortekusa.com/`. (One must join the site to see the
#' manuals.)
#'
#' 2. The Nortek Knowledge Center
#' `https://www.nortekusa.com/en/knowledge-center` may be of help if
#' problems arise in dealing with data from Nortek instruments.
#'
#' @template adpTemplate
#'
#' @template encodingIgnoredTemplate
#'
#' @template adReadingMethodTemplate
#'
#' @author Dan Kelley
#'
#' @family things related to adp data
#' @family functions that read adp data
read.aquadoppProfiler <- function(
file,
from = 1,
to,
by = 1,
tz = getOption("oceTz"),
longitude = NA,
latitude = NA,
orientation,
distance,
monitor = FALSE,
despike = FALSE,
encoding = NA,
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, "\"")
}
}
return(read.adp.nortek(file,
from = from, to = to, by = by, tz = tz,
longitude = longitude, latitude = latitude,
type = "aquadoppProfiler",
orientation = orientation, distance = distance,
encoding = encoding,
monitor = monitor, despike = despike, processingLog = processingLog,
debug = getOption("oceDebug"), ...
))
}
#' Read an adp File in Nortek Format
#'
#' @param despike a logical value indicating whether to use [despike()] to remove
#' anomalous spikes in heading, etc.
#'
#' @param type a character string indicating the type of instrument.
#' If NULL (the default), then [oceMagic()] is used to infer the type. If
#' non-NULL, then the value must be one of: `"aquadoppHR"`,
#' `"aquadoppProfiler"`, `"aquadopp"`, or `"aquadoppPlusMagnetometer"`.
#'
#' @param orientation an 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"`.
#'
#' @param distance an optional vector holding the distances of bin centres from
#' the sensor. This argument is ignored except for Nortek profilers, and need
#' not be given if the function determines the distances correctly from the
#' data. The problem is that the distance is poorly documented in the Nortek
#' System Integrator Guide (2008 edition, page 31), so the function must rely
#' on word-of-mouth formulae that do not work in all cases.
#'
#' @references
#' 1. Information on Nortek profilers (including the System Integrator Guide,
#' which explains the data format byte-by-byte) is available at
#' `https://www.nortekusa.com/`. (One must join the site to see the
#' manuals.)
#'
#' 2. The Nortek Knowledge Center
#' `https://www.nortekusa.com/en/knowledge-center` may be of help if
#' problems arise in dealing with data from Nortek instruments.
#'
#' @template adpTemplate
#'
#' @template encodingIgnoredTemplate
#'
#' @template adReadingMethodTemplate
#'
#' @author Dan Kelley
#'
#' @family things related to adp data
#' @family functions that read adp data
read.adp.nortek <- function(
file, from = 1, to, by = 1, tz = getOption("oceTz"),
longitude = NA, latitude = NA,
type = NULL,
orientation, distance,
encoding = NA,
monitor = FALSE, despike = FALSE, processingLog,
debug = getOption("oceDebug"),
...) {
oceDebug(debug, "read.adp.nortek(...)\n", sep = "", unindent = 1)
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 (is.null(type)) {
oceDebug(debug, "inferring type using oceMagic()\n")
type <- gsub(".*/", "", oceMagic(file))
oceDebug(debug, "inferred type as \"", type, "\"\n")
}
typeAllowed <- c("aquadopp", "aquadoppHR", "aquadoppPlusMagnetometer", "aquadoppProfiler")
# message("typeAllowed= as follows\n");print(typeAllowed)
if (!type %in% typeAllowed) {
stop(
"type \"", type, "\" not permitted; try one of \"",
paste(typeAllowed, collapse = "\" \""), "\""
)
}
# message("DAN 1")
if (!interactive()) {
monitor <- FALSE
}
# degToRad <- atan2(1, 1) / 45
profileStart <- NULL # prevents scope warning in rstudio; defined later anyway
bisectAdpNortek <- function(buf, t.find, add = 0, debug = 0) {
oceDebug(debug, "bisectAdpNortek(t.find=", format(t.find), ", add=", add, "\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)
minute <- bcdToInteger(buf[profileStart[middle] + 4])
second <- bcdToInteger(buf[profileStart[middle] + 5])
day <- bcdToInteger(buf[profileStart[middle] + 6])
hour <- bcdToInteger(buf[profileStart[middle] + 7])
year <- bcdToInteger(buf[profileStart[middle] + 8])
year <- year + ifelse(year >= 90, 1900, 2000)
month <- bcdToInteger(buf[profileStart[middle] + 9])
sec1000 <- bcdToInteger(buf[profileStart[middle] + 10])
t <- ISOdatetime(year, month, day, hour, minute, second + sec1000 / 1000, tz = tz)
oceDebug(
debug, "t=", format(t), "| (from data",
sprintf("%4d-%02d-%02d", year, month, day),
sprintf("%02d:%02d:%02d.%03d", hour, minute, second, sec1000), ") | pass",
format(pass, width = 2), "/", passes, " | middle=", middle, "(", middle / upper * 100, "%)\n"
)
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
minute <- bcdToInteger(buf[profileStart[middle] + 4])
second <- bcdToInteger(buf[profileStart[middle] + 5])
day <- bcdToInteger(buf[profileStart[middle] + 6])
hour <- bcdToInteger(buf[profileStart[middle] + 7])
year <- bcdToInteger(buf[profileStart[middle] + 8])
year <- year + ifelse(year >= 90, 1900, 2000)
month <- bcdToInteger(buf[profileStart[middle] + 9])
sec1000 <- bcdToInteger(buf[profileStart[middle] + 10])
t <- ISOdatetime(year, month, day, hour, minute, second + sec1000 / 1000, tz = tz)
oceDebug(debug, "result: t=", format(t), " at vsd.start[", middle, "]=", profileStart[middle], "\n")
return(list(index = middle, time = t))
}
if (missing(to)) {
to <- NA
} # will catch this later
if (identical(type, "aquadoppPlusMagnetometer")) {
oceDebug(debug, "Reading an aquadopp file which includes raw magnetometer data\n")
}
res <- new("adp")
# fromKeep <- from
# toKeep <- to
# syncCode <- as.raw(0xa5)
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))
}
seek(file, 0, "start")
seek(file, 0, "start")
# go to the end, so the next seek (to get to the data) reveals file length
seek(file, where = 0, origin = "end")
fileSize <- seek(file, where = 0)
oceDebug(debug, "fileSize=", fileSize, "\n")
buf <- readBin(file, what = "raw", n = fileSize, size = 1)
header <- decodeHeaderNortek(buf, type = type, debug = debug - 1)
# averagingInterval <- header$user$averagingInterval
numberOfBeams <- header$numberOfBeams
numberOfCells <- header$numberOfCells
# bin1Distance <- header$bin1Distance
# xmitPulseLength <- header$xmitPulseLength
# cellSize <- header$cellSize
# profilesInFile <- readBin(buf[header$offset + 2:3], what="integer", n=1, size=2, endian="little")
oceDebug(debug, "profile data at buf[", header$offset, "] et seq.\n")
oceDebug(debug, "matching bytes: 0x", buf[header$offset], " 0x", buf[header$offset + 1], " 0x", buf[header$offset + 2], "\n", sep = "")
profileStart <- .Call("match3bytes", buf, buf[header$offset], buf[header$offset + 1], buf[header$offset + 2])
profilesInFile <- length(profileStart)
if (is.na(to)) {
to <- profilesInFile
}
oceDebug(debug, "profilesInFile=", profilesInFile, "\n")
measurementStart <- ISOdatetime(2000 + bcdToInteger(buf[profileStart[1] + 8]), # year FIXME: have to check if before 1990
bcdToInteger(buf[profileStart[1] + 9]), # month
bcdToInteger(buf[profileStart[1] + 6]), # day
bcdToInteger(buf[profileStart[1] + 7]), # hour
bcdToInteger(buf[profileStart[1] + 4]), # min
bcdToInteger(buf[profileStart[1] + 5]), # sec
tz = tz
)
measurementEnd <- ISOdatetime(2000 + bcdToInteger(buf[profileStart[profilesInFile] + 8]), # year FIXME: have to check if before 1990
bcdToInteger(buf[profileStart[profilesInFile] + 9]), # month
bcdToInteger(buf[profileStart[profilesInFile] + 6]), # day
bcdToInteger(buf[profileStart[profilesInFile] + 7]), # hour
bcdToInteger(buf[profileStart[profilesInFile] + 4]), # min
bcdToInteger(buf[profileStart[profilesInFile] + 5]), # sec
tz = tz
)
measurementDeltat <- as.numeric(ISOdatetime(2000 + bcdToInteger(buf[profileStart[2] + 8]), # year FIXME: have to check if before 1990
bcdToInteger(buf[profileStart[2] + 9]), # month
bcdToInteger(buf[profileStart[2] + 6]), # day
bcdToInteger(buf[profileStart[2] + 7]), # hour
bcdToInteger(buf[profileStart[2] + 4]), # min
bcdToInteger(buf[profileStart[2] + 5]), # sec
tz = tz
)) - as.numeric(measurementStart)
oceDebug(debug, "ORIG measurement.deltat=", measurementDeltat, "\n")
measurementDeltat <- (as.numeric(measurementEnd) - as.numeric(measurementStart)) / profilesInFile
oceDebug(debug, "NEW measurementDeltat=", measurementDeltat, "\n")
if (inherits(from, "POSIXt")) {
if (!inherits(to, "POSIXt")) {
stop("if 'from' is POSIXt, then 'to' must be, also")
}
fromPair <- bisectAdpNortek(buf, from, -1, debug - 1)
from <- fromIndex <- fromPair$index
toPair <- bisectAdpNortek(buf, to, 1, debug - 1)
to <- toIndex <- toPair$index
oceDebug(
debug, " from=", format(fromPair$t), " yields profileStart[", fromIndex, "]\n",
" to =", format(toPair$t), " yields profileStart[", toIndex, "]\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"
)
time1 <- ISOdatetime(2000 + bcdToInteger(buf[profileStart[1] + 8]), # year FIXME: have to check if before 1990
bcdToInteger(buf[profileStart[1] + 9]), # month
bcdToInteger(buf[profileStart[1] + 6]), # day
bcdToInteger(buf[profileStart[1] + 7]), # hour
bcdToInteger(buf[profileStart[1] + 4]), # min
bcdToInteger(buf[profileStart[1] + 5]), # sec
tz = tz
)
time2 <- ISOdatetime(2000 + bcdToInteger(buf[profileStart[2] + 8]), # year FIXME: have to check if before 1990
bcdToInteger(buf[profileStart[2] + 9]), # month
bcdToInteger(buf[profileStart[2] + 6]), # day
bcdToInteger(buf[profileStart[2] + 7]), # hour
bcdToInteger(buf[profileStart[2] + 4]), # min
bcdToInteger(buf[profileStart[2] + 5]), # sec
tz = tz
)
dt <- as.numeric(difftime(time2, time1, units = "secs"))
oceDebug(debug, "dt=", dt, "s; at this stage, by=", by, "(not interpreted yet)\n")
profileStart <- profileStart[profileStart[fromIndex] < profileStart & profileStart < profileStart[toIndex]]
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:",
paste(profileStart[1:10], " "), "\n"
)
} else {
fromIndex <- from
toIndex <- to
if (toIndex < 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:", paste(profileStart[1:10], " "), "\n")
}
profilesToRead <- length(profileStart)
oceDebug(debug, "profilesToRead=", profilesToRead, "\n")
profileStart2 <- sort(c(profileStart, profileStart + 1)) # use this to subset for 2-byte reads
numberOfCells <- header$user$numberOfCells
numberOfBeams <- header$head$numberOfBeams
oceDebug(debug, "numberOfCells=", numberOfCells, "\n")
oceDebug(debug, "numberOfBeams=", numberOfBeams, "\n")
items <- numberOfCells * numberOfBeams
time <- ISOdatetime(2000 + bcdToInteger(buf[profileStart + 8]), # year FIXME: have to check if before 1990
bcdToInteger(buf[profileStart + 9]), # month
bcdToInteger(buf[profileStart + 6]), # day
bcdToInteger(buf[profileStart + 7]), # hour
bcdToInteger(buf[profileStart + 4]), # min
bcdToInteger(buf[profileStart + 5]), # sec
tz = tz
)
class(time) <- c("POSIXt", "POSIXct") # FIXME do we need this?
attr(time, "tzone") <- getOption("oceTz") # Q: does file hold the zone?
# aquadopp error: see table 5.4 (p40) and table 5.10 (p53) of system-integrator-manual_jan2011.pdf
error <- readBin(buf[profileStart2 + 10], what = "integer", n = profilesToRead, size = 2, endian = "little", signed = FALSE)
heading <- 0.1 * readBin(buf[profileStart2 + 18], what = "integer", n = profilesToRead, size = 2, endian = "little", signed = TRUE)
pitch <- 0.1 * readBin(buf[profileStart2 + 20], what = "integer", n = profilesToRead, size = 2, endian = "little", signed = TRUE)
roll <- 0.1 * readBin(buf[profileStart2 + 22], what = "integer", n = profilesToRead, size = 2, endian = "little", signed = TRUE)
pressure.MSB <- readBin(buf[profileStart + 24], what = "integer", n = profilesToRead, size = 1, endian = "little", signed = FALSE)
pressure.LSW <- readBin(buf[profileStart2 + 26], what = "integer", n = profilesToRead, size = 2, endian = "little", signed = FALSE)
pressure <- (as.integer(pressure.MSB) * 65536 + pressure.LSW) * 0.001 # CHECK
temperature <- 0.01 * readBin(buf[profileStart2 + 28], what = "integer", n = profilesToRead, size = 2, endian = "little")
# if type=aquadoppPlusMagnetometer read the extra fields -- issue 1758, and SIG 2020:
if (type == "aquadoppPlusMagnetometer") {
soundSpeed <- 0.1 * readBin(buf[profileStart2 + 30], what = "integer", n = profilesToRead, size = 2, endian = "little")
ensCount <- readBin(buf[profileStart2 + 32], what = "integer", n = profilesToRead, size = 2, endian = "little")
compHx <- readBin(buf[profileStart2 + 34], what = "integer", n = profilesToRead, size = 2, endian = "little")
compHy <- readBin(buf[profileStart2 + 36], what = "integer", n = profilesToRead, size = 2, endian = "little")
compHz <- readBin(buf[profileStart2 + 38], what = "integer", n = profilesToRead, size = 2, endian = "little")
}
v <- array(double(), dim = c(profilesToRead, numberOfCells, numberOfBeams))
a <- array(raw(), dim = c(profilesToRead, numberOfCells, numberOfBeams)) # echo amplitude
q <- array(raw(), dim = c(profilesToRead, numberOfCells, numberOfBeams)) # correlation
velocityScale <- 1e-3 # FIXME: why not use the value in user$velocityScale?
# FIXME: why does 54 work, given 53 in docs? [see 38 of System Integrator Guide]
oShift <- switch(type,
aquadoppHR = 54,
aquadoppProfiler = 30,
aquadopp = 30,
aquadoppPlusMagnetometer = 40
)
for (i in 1:profilesToRead) {
o <- profileStart[i] + oShift
# oceDebug(debug, "getting data chunk",i," at file position',o,"\n")
v[i, , ] <- velocityScale * matrix(readBin(buf[o + seq(0, 2 * items - 1)], "integer",
n = items, size = 2,
endian = "little", signed = TRUE
), ncol = numberOfBeams, byrow = FALSE)
# BUG BUG BUG why difference in byrow for v vs a and q?
o <- o + items * 2
a[i, , ] <- matrix(buf[o + seq(0, items - 1)], ncol = items, byrow = TRUE)
o <- o + items
q[i, , ] <- matrix(buf[o + seq(0, items - 1)], ncol = items, byrow = TRUE) # FIXME: this is correlation, not quality
if (monitor) {
cat(".", ...)
if (!(i %% 50)) cat(i, "\n", ...)
}
}
if (monitor) cat("\nRead", profilesToRead, "of the", profilesInFile, "profiles in", filename, "\n", ...)
if (missing(distance)) {
distance <- seq(header$user$blankingDistance + header$user$cellSize, by = header$user$cellSize, length.out = header$user$numberOfCells)
} else {
if (length(distance) != dim(v)[2]) {
stop("the argument distance is of length ", length(distance), ",
which does not match the second dimension of the velocity matrix, ", dim(v)[2])
}
}
# get diagnostic data, if any, and trim them to same index range as conventional data
if (type == "aquadopp") {
diaStart <- .Call("match3bytes", buf, 0xa5, 0x80, 0x15)
oceDebug(debug, "diaStart range:", range(diaStart), "\n")
diaStart <- subset(diaStart, diaStart >= profileStart[fromIndex])
diaStart <- subset(diaStart, diaStart <= profileStart[toIndex])
oceDebug(debug, "LATER diaStart range:", range(diaStart), "\n")
diaToRead <- length(diaStart)
diaStart2 <- sort(c(diaStart, diaStart + 1))
timeDia <- ISOdatetime(2000 + bcdToInteger(buf[diaStart + 8]),
bcdToInteger(buf[diaStart + 9]), # month
bcdToInteger(buf[diaStart + 6]), # day
bcdToInteger(buf[diaStart + 7]), # hour
bcdToInteger(buf[diaStart + 4]), # min
bcdToInteger(buf[diaStart + 5]), # sec
tz = tz
)
# aquadopp error: see table 5.4 (p40) and table 5.10 (p53) of system-integrator-manual_jan2011.pdf
errorDia <- readBin(buf[diaStart2 + 10], what = "integer", n = diaToRead, size = 2, endian = "little", signed = FALSE)
headingDia <- 0.1 * readBin(buf[diaStart2 + 18], what = "integer", n = diaToRead, size = 2, endian = "little", signed = TRUE)
pitchDia <- 0.1 * readBin(buf[diaStart2 + 20], what = "integer", n = diaToRead, size = 2, endian = "little", signed = TRUE)
rollDia <- 0.1 * readBin(buf[diaStart2 + 22], what = "integer", n = diaToRead, size = 2, endian = "little", signed = TRUE)
pressureMSB <- readBin(buf[diaStart + 24], what = "integer", n = diaToRead, size = 1, endian = "little", signed = FALSE)
pressureLSW <- readBin(buf[diaStart2 + 26], what = "integer", n = diaToRead, size = 2, endian = "little", signed = FALSE)
pressureDia <- (as.integer(pressureMSB) * 65536 + pressureLSW) * 0.001
temperatureDia <- 0.01 * readBin(buf[diaStart2 + 28], what = "integer", n = diaToRead, size = 2, endian = "little")
vDia <- array(double(), dim = c(diaToRead, 1, 3))
vDia[, , 1] <- 0.001 * readBin(buf[diaStart2 + 30], what = "integer", n = diaToRead, size = 2, endian = "little", signed = TRUE)
vDia[, , 2] <- 0.001 * readBin(buf[diaStart2 + 32], what = "integer", n = diaToRead, size = 2, endian = "little", signed = TRUE)
vDia[, , 3] <- 0.001 * readBin(buf[diaStart2 + 34], what = "integer", n = diaToRead, size = 2, endian = "little", signed = TRUE)
aDia <- array(raw(), dim = c(diaToRead, 1, 3))
aDia[, , 1] <- buf[diaStart + 36]
aDia[, , 2] <- buf[diaStart + 37]
aDia[, , 3] <- buf[diaStart + 38]
} else if (type == "aquadoppPlusMagnetometer") {
diaStart <- .Call("match3bytes", buf, 0xa5, 0x80, 0x15)
oceDebug(debug, "diaStart range:", range(diaStart), "\n")
diaStart <- subset(diaStart, diaStart >= profileStart[fromIndex])
diaStart <- subset(diaStart, diaStart <= profileStart[toIndex])
oceDebug(debug, "LATER diaStart range:", range(diaStart), "\n")
diaToRead <- length(diaStart)
diaStart2 <- sort(c(diaStart, diaStart + 1))
timeDia <- ISOdatetime(2000 + bcdToInteger(buf[diaStart + 8]),
bcdToInteger(buf[diaStart + 9]), # month
bcdToInteger(buf[diaStart + 6]), # day
bcdToInteger(buf[diaStart + 7]), # hour
bcdToInteger(buf[diaStart + 4]), # min
bcdToInteger(buf[diaStart + 5]), # sec
tz = tz
)
# aquadopp error: see table 5.4 (p40) and table 5.10 (p53) of system-integrator-manual_jan2011.pdf
errorDia <- readBin(buf[diaStart2 + 10], what = "integer", n = diaToRead, size = 2, endian = "little", signed = FALSE)
headingDia <- 0.1 * readBin(buf[diaStart2 + 18], what = "integer", n = diaToRead, size = 2, endian = "little", signed = TRUE)
pitchDia <- 0.1 * readBin(buf[diaStart2 + 20], what = "integer", n = diaToRead, size = 2, endian = "little", signed = TRUE)
rollDia <- 0.1 * readBin(buf[diaStart2 + 22], what = "integer", n = diaToRead, size = 2, endian = "little", signed = TRUE)
pressureMSB <- readBin(buf[diaStart + 24], what = "integer", n = diaToRead, size = 1, endian = "little", signed = FALSE)
pressureLSW <- readBin(buf[diaStart2 + 26], what = "integer", n = diaToRead, size = 2, endian = "little", signed = FALSE)
pressureDia <- (as.integer(pressureMSB) * 65536 + pressureLSW) * 0.001
temperatureDia <- 0.01 * readBin(buf[diaStart2 + 28], what = "integer", n = diaToRead, size = 2, endian = "little")
# nolint start object_usage_linter
ensCountDia <- readBin(buf[profileStart2 + 32], what = "integer", n = profilesToRead, size = 2, endian = "little")
# nolint end object_usage_linter
soundSpeedDia <- 0.1 * readBin(buf[profileStart2 + 30], what = "integer", n = profilesToRead, size = 2, endian = "little")
compHxDia <- readBin(buf[profileStart2 + 34], what = "integer", n = profilesToRead, size = 2, endian = "little")
compHyDia <- readBin(buf[profileStart2 + 36], what = "integer", n = profilesToRead, size = 2, endian = "little")
compHzDia <- readBin(buf[profileStart2 + 38], what = "integer", n = profilesToRead, size = 2, endian = "little")
vDia <- array(double(), dim = c(diaToRead, 1, 3))
vDia[, , 1] <- 0.001 * readBin(buf[diaStart2 + 40], what = "integer", n = diaToRead, size = 2, endian = "little", signed = TRUE)
vDia[, , 2] <- 0.001 * readBin(buf[diaStart2 + 42], what = "integer", n = diaToRead, size = 2, endian = "little", signed = TRUE)
vDia[, , 3] <- 0.001 * readBin(buf[diaStart2 + 44], what = "integer", n = diaToRead, size = 2, endian = "little", signed = TRUE)
aDia <- array(raw(), dim = c(diaToRead, 1, 3))
aDia[, , 1] <- buf[diaStart + 46]
aDia[, , 2] <- buf[diaStart + 47]
aDia[, , 3] <- buf[diaStart + 48]
}
res@data <- list(
v = v,
a = a,
q = q,
distance = distance,
time = time,
pressure = pressure,
error = error,
temperature = temperature,
heading = heading,
pitch = pitch,
roll = roll
)
if (type == "aquadoppPlusMagnetometer") {
res@data <- c(
res@data,
list(
soundSpeed = soundSpeed,
ensCount = ensCount,
compHx = compHx,
compHy = compHy,
compHz = compHz
)
)
}
# Sometimes there can be an extra sample, with a time of
# NA. Because the times are continuous without the extra sample,
# for now I'll just remove the entire sample
tNA <- which(is.na(time))
if (length(tNA) > 0) {
for (field in names(res@data)) {
if (!(field %in% "distance")) {
if (field %in% c("v", "a", "q")) {
res@data[[field]] <- res@data[[field]][-tNA, , , drop = FALSE]
} else {
res@data[[field]] <- res@data[[field]][-tNA]
}
}
}
warning(paste("Found and removed", length(tNA), "NAs in the time vector."))
}
if (type == "aquadopp" || type == "aquadoppPlusMagnetometer" && diaToRead > 0) {
# FIXME: there may be other things here, e.g. does it try to measure salinity?
res@data$timeDia <- timeDia
res@data$errorDia <- errorDia
res@data$headingDia <- headingDia
res@data$pitchDia <- pitchDia
res@data$rollDia <- rollDia
res@data$pressureDia <- pressureDia
res@data$temperatureDia <- temperatureDia
res@data$vDia <- vDia
res@data$aDia <- aDia
if (type == "aquadoppPlusMagnetometer") {
res@data$soundSpeedDia <- soundSpeedDia
res@data$compHxDia <- compHxDia
res@data$compHyDia <- compHyDia
res@data$compHzDia <- compHzDia
}
}
if (missing(orientation)) {
orientation <- header$head$configTiltSensorOrientation # FIXME: is there a 'header$head'?
} else {
orientation <- match.arg(orientation, c("sideward", "upward", "downward"))
}
res@metadata$manufacturer <- "nortek"
res@metadata$fileType <- type
res@metadata$filename <- filename
res@metadata$latitude <- latitude
res@metadata$longitude <- longitude
res@metadata$numberOfSamples <- dim(v)[1]
res@metadata$numberOfCells <- dim(v)[2]
res@metadata$numberOfBeams <- dim(v)[3]
res@metadata$numberOfBeamSequencesPerBurst <- header$user$numberOfBeamSequencesPerBurst
res@metadata$measurementStart <- measurementStart
res@metadata$measurementEnd <- measurementEnd
res@metadata$measurementDeltat <- measurementDeltat
res@metadata$subsampleStart <- time[1]
res@metadata$subsampleEnd <- time[length(time)]
res@metadata$subsampleDeltat <- as.numeric(time[2]) - as.numeric(time[1])
res@metadata$size <- header$head$size
res@metadata$serialNumber <- header$hardware$serialNumber
res@metadata$internalCodeVersion <- header$hardware$picVersion
res@metadata$hardwareRevision <- header$hardware$hwRevision
res@metadata$recSize <- header$hardware$recSize
res@metadata$velocityRange <- header$hardware$velocityRange # FIXME: should check against velocityMaximum
res@metadata$firmwareVersion <- header$hardware$fwVersion
res@metadata$config <- header$hardware$config
res@metadata$configPressureSensor <- header$head$configPressureSensor
res@metadata$configMagnetometerSensor <- header$head$configMagnetometerSensor
res@metadata$configTiltSensor <- header$head$configTiltSensor
res@metadata$beamAngle <- 25 # FIXME: may change with new devices
res@metadata$tiltSensorOrientation <- header$head$tiltSensorOrientation
res@metadata$orientation <- orientation
res@metadata$frequency <- header$head$frequency
res@metadata$headSerialNumber <- header$head$headSerialNumber
res@metadata$bin1Distance <- header$user$blankingDistance # FIXME: is this right?
res@metadata$blankingDistance <- header$user$blankingDistance
res@metadata$measurementInteres <- header$user$measurementInteres
res@metadata$transformationMatrix <- header$head$transformationMatrix
res@metadata$deploymentName <- header$user$deploymentName
res@metadata$cellSize <- header$user$cellSize
res@metadata$velocityResolution <- velocityScale
res@metadata$velocityMaximum <- velocityScale * 2^15
res@metadata$originalCoordinate <- header$user$originalCoordinate
res@metadata$oceCoordinate <- header$user$originalCoordinate
res@metadata$oceBeamUnspreaded <- FALSE
res@metadata$units$v <- list(unit = expression(m / s), scale = "")
res@metadata$units$distance <- list(unit = expression(m), scale = "")
res@metadata$units$pressure <- list(unit = expression(dbar), scale = "")
res@metadata$units$salinity <- list(unit = expression(), scale = "PSS-78")
res@metadata$units$temperature <- list(unit = expression(degree * C), scale = "ITS-90")
res@metadata$units$soundSpeed <- list(unit = expression(m / s), scale = "")
res@metadata$units$heading <- list(unit = expression(degree), scale = "")
res@metadata$units$pitch <- list(unit = expression(degree), scale = "")
res@metadata$units$roll <- list(unit = expression(degree), scale = "")
res@metadata$units$headingStd <- list(unit = expression(degree), scale = "")
res@metadata$units$pitchStd <- list(unit = expression(degree), scale = "")
res@metadata$units$rollStd <- list(unit = expression(degree), scale = "")
res@metadata$units$attitude <- list(unit = expression(degree), scale = "")
if (missing(processingLog)) {
processingLog <- paste("read.adp.nortek(file=\"", filename, "\", from=", from, ", to=", to, ", by=", by, ")", sep = "")
}
res@processingLog <- processingLogItem(processingLog)
res
} # read.adp.nortek()
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.