#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Aquascat instruments (:mod:`hydrac.instruments.aquascat`)
=========================================================
.. autoclass:: Aquascat
:members:
:private-members:
"""
# from .instruments import instrument
from .acoustique import Acoustic
from hydrac.util import constants
import numpy as np
import struct as st
import glob
import os
import datetime
from tkinter import filedialog
from tkinter import Tk
from hydrac.util.parameters import Parameters
from hydrac.util.paramcontainer import ParamContainer
from hydrac.util.query import query_yes_no, query_number
[docs]class Aquascat(Acoustic):
""" Aquascat instrument class.
Base class : :mod:`hydrac.instruments.acoustique.Acoustic`
The Aquascat is a multifrequency profiler working at the kHz-MHz frequency
range developpe by Aquatec Group. Its full description is given in
http://www.aquatecgroup.com/11-products/78-aquascat. The example of a
Aquascat 1000R is shown below :
.. image:: ../_image/aqa.png
:width: 50%
:align: center
The Aquascat class reads the raw binary files from Aquascat instruments
and stores the valuable information into the ``param`` modified dictionnary
with the common shape handled by hydrac
(see :mod:`hydrac.instruments.acoustique`).
A typical example of how should the data considered is as follows ::
>>> A = Aquascat('Campaign_1')
The above line will prompt the user to select one or multiple files in a
directory, read and store each file data into ``param``
instanciated in :class: `hydrac.instruments.acoustique.Acoustic`, as
separate modified dictionnaries ``PX``, X being the file number. Ex: for
the first file loaded in ``param``, one can look at the differents
variables stored in ``param.P0`` ::
>>> A.param.P0.keys()
dict_keys(['BurstTime', 'PingRate', 'NumPingsTot',
'NumChannels', 'NumProfilesSamples', 'NumAverage'...])
One also gets the base Acoustic class instance ``instr_type``::
>>> A.instr_type
'acoustic'
Or methods::
>>> A.preproc_acoustic_data()
which will prompt the user for extra information necessary to the
processing of the hydroacoutic data such as water temperature. See
:func:`hydrac.instruments.acoustique.Acoustic.preproc_acoustic_data`
Parameters
----------
name : str, {'Campaign_1'}
Tag name, for personal usage, for instance to add a campaign name."""
def __init__(self, name):
Acoustic.__init__(self)
self.instr_name = 'aquascat'
self.name = name
self.tag_M = 'multi'
self.nearfield_min = 0.2
self.filepath, self.tempdir = self.file_select()
self.batch_read()
[docs] def file_select(self):
""" User input selection of the target files to read """
root = Tk()
filez = filedialog.askopenfilenames(parent=root, title='Choose a file')
root.withdraw() # use to hide tkinter window
filepath = list(filez)
filepath.sort()
filepath = tuple(filepath)
tempdir = os.path.dirname(os.path.abspath(filez[0]))
return filepath, tempdir
[docs] def param_shape_AQA(self, inp):
"""Reshapes the modified dictionnary ``param`` to fit the standard key
names and data shape used throughout hydrac"""
FNAME = ['BurstTime', 'PingRate', 'NumPings', 'NumAbsTimeSlots',
'AbsNumProfiles', 'AbsAverage', 'AbsProfileRate',
'AbsBinLengthMM', 'AbsStartBin', 'AbsNumBins', 'AbsRxChan',
'AbsTxChan', 'AbsTransducerName', 'AbsTransducerRadius',
'AbsTransducerBeamWidth', 'AbsTransducerKt', 'AbsTxFrequency',
'AbsRxFrequency', 'AbsTxPulseLength', 'AbsStartingGain',
'AbsTVG', 'AbsPowerLevel', 'AbsMean', 'AbsData',
'AbsBinRange', 'serial_num']
FNAME_out = ['BurstTime', 'PingRate', 'NumPingsTot', 'NumChannels',
'NumProfilesSamples', 'NumAverage', 'SampleRate',
'BinLengthMM', 'StartBin', 'NumBins', 'TransducerTag',
'TransducerTag2', 'TransducerName', 'TransducerRadius',
'TransducerBeamWidth', 'TransducerKt', 'Frequency',
'FrequencyRx', 'PulseLength', 'RxGain', 'TVG', 'TxGain',
'MeanProfile', 'RawIntensity', 'BinRange', 'SerialNum']
FNAME_aux = ['AuxData', 'AuxNumSamples', 'AuxChannelName',
'AuxChannelUnit', 'NumAuxChans', 'AuxSampleRate']
FNAME_aux_out = ['AuxData', 'AuxNumSamples', 'AuxChannelName',
'AuxChannelUnit', 'NumAuxChannels', 'AuxSampleRate']
inp2 = Parameters({'aquascat': {}})
[inp2.update({FNAME_out[tui]:inp[
FNAME[tui]]}) for tui in range(0, len(FNAME))]
[inp2.aquascat.update({FNAME_aux_out[tui]: inp[
FNAME_aux[tui]]}) for tui in range(0, len(FNAME_aux))]
inp2.time = np.arange(inp.BurstTime.timestamp(),
inp.BurstTime.timestamp() +
inp.NumPings / inp.AbsAverage[0] /
inp.AbsProfileRate[0],
1 / inp.AbsProfileRate[0])
# attention au cas où il n'y ait pas de AuxData => set tout à empty
if inp2.aquascat.NumAuxChannels == 0:
inp.AuxData = np.zeros((1, 8))
inp2.pressure = np.tile(inp.AuxData[:, 5],
(inp2.NumChannels, 1)).T
inp2.depth_raw = inp2.pressure * 0
inp2.temp = inp.AuxData[:, 4]
inp2.battery = inp.AuxData[:, -1]
inp2.BinRange = np.array(inp2.BinRange).T
inp2.StartBin = inp2.BinRange[np.int(inp2.StartBin[0]), :]
inp2.TransducerAngle = np.zeros(np.shape(inp2.StartBin))
inp2.sal = np.zeros(np.shape(inp2.battery))
inp2.TVG = np.zeros(np.shape(inp2.BinRange))
inp2.NumSamples = inp2.NumPingsTot / inp2.NumAverage
inp2.BurstTime = datetime.datetime.timestamp(inp2.BurstTime)
inp2.RawIntensity = inp2.RawIntensity**2
if inp2.aquascat.NumAuxChannels == 0:
inp2.depth = 0 * np.arange(0, inp2.NumSamples[0])
inp2.aquascat.time_aux = []
else:
inp2.depth = np.interp(np.arange(0, inp2.NumSamples[0]) /
inp2.NumSamples[0],
np.arange(0, len(inp2.pressure[:, 0])) /
len(inp2.pressure[:, 0]),
inp2.aquascat.AuxData[:, 5])
inp2.aquascat.time_aux = np.arange(inp.BurstTime.timestamp(),
inp.BurstTime.timestamp() +
inp2.aquascat.AuxNumSamples /
inp2.aquascat.AuxSampleRate,
1 / inp2.aquascat.AuxSampleRate)
self.paramvarnames = ['time', 'pressure', 'temp', 'battery', 'sal',
'time_aux', 'NumSamples']+FNAME_out+FNAME_aux_out
inp2.FileName = os.path.split(self.currfilepath)[-1]
inp2.calib = self._set_calib_aqa()
try:
TransducerKt = inp2.calib['SN' + str(inp2.SerialNum[0]) +
str(inp2.SerialNum[1])].Kt
inp2.TransducerKt = []
[inp2.TransducerKt.append([a, b]) for a, b in TransducerKt
if a in inp2.FrequencyRx]
except(AttributeError):
pass
ex = False
if 'SN' + str(inp2.SerialNum[0]) +\
str(inp2.SerialNum[1]) not in inp2.calib.keys():
inp2.calib.update({'SN' + str(inp2.SerialNum[0]) +
str(inp2.SerialNum[1]): Parameters({})})
ex = True
if query_yes_no('No Kt values for AQA SN' + str(
inp2.SerialNum[0]) +
str(inp2.SerialNum[1]) + ' - Do you want to' +
' enter Kt manually ? '):
Kt = list()
for ui in range(len(inp2.FrequencyRx)):
Kt.append((inp2.FrequencyRx[ui],
query_number('Kt channel F = ' +
str(inp2.FrequencyRx[ui]))))
inp2.calib['SN' + str(inp2.SerialNum[0]) +
str(inp2.SerialNum[1])].update({'Kt': Kt})
else:
Kt = []
for ui in range(len(inp2.FrequencyRx)):
Kt.append((inp2.FrequencyRx[ui], 1))
inp2.calib['SN' + str(inp2.SerialNum[0]) +
str(inp2.SerialNum[1])].update({'Kt': Kt})
if query_yes_no('No depth calibration values for AQA SN' +
str(inp2.SerialNum[0]) + str(inp2.SerialNum[1]) +
' - Do you want to enter depth calib manually ? '):
inp2.calib['SN' + str(inp2.SerialNum[0]) +
str(inp2.SerialNum[1])].update({'cal_depth1':
{'a':
query_number('a= '),
'b':
query_number('b= ')
}})
else:
inp2.calib['SN' + str(inp2.SerialNum[0]) +
str(inp2.SerialNum[1])].update({'cal_depth1':
{'a': 1, 'b': 0}})
TransducerKt = inp2.calib['SN' + str(inp2.SerialNum[0]) +
str(inp2.SerialNum[1])].Kt
inp2.TransducerKt = []
[inp2.TransducerKt.append([a, b]) for a, b in TransducerKt
if a in inp2.FrequencyRx]
hu=query_yes_no('Do you wosh to add other transducer Kts ?')
Fx=[]
while hu is True:
Fx.append((query_number('Choose a frequency ?'), query_number('Choose a Kt')))
hu=query_yes_no('Add another channel ?')
[inp2.calib['SN' + str(inp2.SerialNum[0]) +
str(inp2.SerialNum[1])].Kt.append(ui) for ui in Fx]
if ex:
if query_yes_no('Do you wish to save current' +
' instrument calibration parameters ?'):
self._update_calib_aqa('SN' + str(inp2.SerialNum[0]) +
str(inp2.SerialNum[1]),
inp2.calib['SN' +
str(inp2.SerialNum[0]) +
str(inp2.SerialNum[1])])
else:
print('je pas au bon endroit')
return inp2
[docs] def _set_calib_aqa(self):
"""Calibration coefficients for each transducer are unique to each
transducer. As such, each instrument comes with a different set of
calibration coefficients Kt. Those are stored in xml file generated
(if non-existant) or updated (with :func:`_update_calib_aqa`) with the
instrument serial number and the Kts. This xml is accessed/generated
through the class :class:`hydrac.util.paramcontainer.ParamContainer`.
"""
rr = glob.glob(constants.dir_aqa_calib + '*.xml')
rr.sort()
try:
aqa = Parameters(ParamContainer(path_file=rr[-1]))
except(IndexError, AttributeError):
aqa = Parameters({})
return aqa
[docs] def _update_calib_aqa(self, a, b=None, remove_elem=False):
"""
Update the calibration xml file either by adding a new instrument
or deleting the coefficients of a particular instrument. This xml is
accessed-generated through the class
:class:`hydrac.util.paramcontainer.ParamContainer`.
Parameters
----------
a : str
Instrument serial number read from binary files and stored
in ``param``.
b=None : dict
Contains the calibration coefficients, prompted by the user if the
instrument is used for the first time in hydrac.
remove_elem=False : bool
Removes the information relative to the serial number a in the xml.
"""
rr = glob.glob(constants.dir_aqa_calib + '*.xml')
rr.sort()
if remove_elem is False:
uiui = ParamContainer(path_file=rr[-1])
uiui._set_child(a, b)
uiui._save_as_xml(path_file=constants.dir_aqa_calib +
'calib_AQA.xml', find_new_name=False,
overwrite=True)
else:
uiui = Parameters(ParamContainer(path_file=rr[-1]))
try:
del uiui[a]
except(KeyError):
pass
yoyo = ParamContainer(tag='CalibAQA')
for a, b in uiui.items():
yoyo._set_child(a, b)
yoyo._save_as_xml(path_file=constants.dir_aqa_calib +
'calib_AQA.xml', find_new_name=False,
overwrite=True)
[docs] def batch_read(self):
"""Read a list of .aqa raw Aquascat file using the method
:func:`read_AQA`"""
for kk in range(len(self.filepath)):
self.currfilepath = self.filepath[kk]
tmp = self.param_shape_AQA(self.read_AQA(self.filepath[kk]))
self.param.update({"P"+str(kk): tmp})
del self.currfilepath
[docs] def psiCal(self, f, at, r):
"""Compute the near field correction for each transducer channel."""
PSI = np.zeros((len(r), len(f)))
lamb = 1500 / f
rn = np.pi * (at ** 2) / lamb
for k in range(0, len(f)):
for i in range(0, len(r)):
if r[i] > 2 * rn[k]:
PSI[i, k] = 1
else:
PSI[i, k] = (1/3) * (2+2 * rn[k] / r[i])
return PSI
[docs] def read_AQA(self, fname):
"""Read one .aqa raw Aquascat file"""
# typeNames = {
# 'int8' :'b',
# 'uint8' :'B',
# 'int16' :'h',
# 'uint16' :'H',
# 'int32' :'i',
# 'uint32' :'I',
# 'int64' :'q',
# 'uint64' :'Q',
# 'float' :'f',
# 'double' :'d',
# 'char' :'s'}
f = open(fname, 'rb')
f.seek(0, 2)
file_size = f.tell()
f.seek(0)
ExpType = 'SINGLE'
def Initialize_Arr(var, a, b, c=1):
if var:
pass
else:
Arr = np.zeros(a, b, c)
return Arr
def readNextAQUAscat1000Header(f, file_size):
PktType_ = f.read(1)
PktType, = st.unpack("B", PktType_)
PktVersion_ = f.read(1)
PktVersion, = st.unpack("B", PktVersion_)
PktSize_ = f.read(2)
PktSize, = st.unpack("H", PktSize_)
Pktchecksum_ = f.read(2)
Pktchecksum, = st.unpack("H", Pktchecksum_)
if f.tell() == file_size:
Status = 0
else:
Status = 1
return Status, PktType, PktSize
def findAQUAscat1000Packet(f, ktype, file_size):
f.seek(0)
while f.tell() != file_size:
Status, PktType, PktSize = readNextAQUAscat1000Header(
f, file_size)
if Status == 0:
break
elif PktType == ktype:
break
else:
f.seek(2 * PktSize, 1)
return Status, PktSize
def cleanString(strIn):
strOut = str(np.zeros(np.size(strIn)))
strTmp = strIn.strip()
strOut = strTmp
return strOut
# Read File Version from File
status, pktSize = findAQUAscat1000Packet(f, 19, file_size)
if status == 1:
sdata = [st.unpack("H", f.read(2))[0] for ui in range(
1, pktSize+1)]
FileVersionMajor = sdata[1]
FileVersionMinor = sdata[2]
else:
FileVersionMajor = 5
FileVersionMinor = 0
del sdata
# Read in the Burst Start Time Information
# Read File Version from File
status, pktSize = findAQUAscat1000Packet(f, 54, file_size)
if status == 1:
sdata = [st.unpack("H", f.read(2))[0] for ui in range(1, 6 + 1)]
WakeSource_ = f.read(2)
WakeSource, = st.unpack("H", WakeSource_)
BurstNumber_ = f.read(4)
BurstNumber, = st.unpack("I", BurstNumber_)
Junk_ = f.read(2)
Junk, = st.unpack("H", Junk_)
BurstTime = datetime.datetime(
sdata[0], sdata[1], sdata[2], sdata[3], sdata[4], sdata[5])
else:
BurstTime = datetime.datetime(0, 0, 0, 0, 0, 0)
BurstNumber = 0
WakeSource = 0
del sdata
# Deal with reading the personality
status, pktSize = findAQUAscat1000Packet(f, 53, file_size)
if status == 0:
print('Error')
pass
else:
PktStartPos = f.tell()
skip_ = f.read(2)
skip, = st.unpack("H", skip_) # Size of Packet
temp, = st.unpack("H", f.read(2))
BoardVersion = print('{}.{}'.format((temp >> 8) & 255,
(temp & 255)))
t1, = st.unpack("H", f.read(2))
t2, = st.unpack("H", f.read(2))
SerialNum = print('{}.{}'.format(t1, t2))
serial_num = [t1, t2]
MemorySize, = st.unpack("H", f.read(2))
LoggerType = [str(st.unpack("s", f.read(1)))[3]
for ui in range(1, 32+1)]
skip, = st.unpack("H", f.read(2))
# Number of ABS channels the system supports
NumAbsChannels, = st.unpack("H", f.read(2))
# Number of AUX channels the system supports (8)
NumAuxChannels, = st.unpack("H", f.read(2))
# Number of ABS Transducer Information that
# exist in personality table
NumAbsTransducers, = st.unpack("H", f.read(2))
# Not recorded at this time
BatteryCapacity, = st.unpack("f", f.read(4))
StandbyPower, = st.unpack("f", f.read(4))
ActivePower, = st.unpack("f", f.read(4))
f.seek(PktStartPos + 112, 0)
# These are the offsets into the Packet
PtrToAuxInfo, = st.unpack("H", f.read(2))
PtrToTransducerInfo, = st.unpack("H", f.read(2))
# Read in the important information for the Aux Channels
# First Need to assign the multiple dimension arrays
AuxChannelName = []
AuxChannelUnit = []
AuxFlags = []
AuxNumGain = []
AuxCalDate = []
AuxNumCoeff = []
AuxGainLabel = []
AuxGainCoeff = []
AuxGainMin = []
AuxGainMax = []
for i in range(1, NumAuxChannels+1):
PtrToThisAux = PktStartPos + PtrToAuxInfo * 2 + 400 * (i-1)
# Move to the start of the ABS information
f.seek(PtrToThisAux, 0)
skip = [st.unpack("H", f.read(2)) for ui in range(1, 2+1)]
TTmp = [str(
st.unpack("s", f.read(1)))[3] for ui in range(1, 16+1)]
sss = ''
TTmp2 = sss.join(TTmp)
TTmp2 = TTmp2.strip("\\")
AuxChannelName = np.append(AuxChannelName, [TTmp2], axis=0)
skip, = st.unpack("H", f.read(2))
TTmp = [str(
st.unpack("s", f.read(1)))[3] for ui in range(1, 8+1)]
sss = ''
TTmp2 = sss.join(TTmp)
TTmp2 = TTmp2.strip("\\")
AuxChannelUnit = np.append(AuxChannelUnit, [TTmp2], axis=0)
skip, = st.unpack("H", f.read(2))
AuxFlags = np.append(AuxFlags, st.unpack("H", f.read(2)))
skip = [st.unpack("H", f.read(2)) for ui in range(1, 2+1)]
AuxNumGain = np.append(
AuxNumGain, st.unpack("H", f.read(2)), axis=0)
skip, = st.unpack("H", f.read(2))
TTmp = [str(
st.unpack("s", f.read(1)))[3] for ui in range(1, 16+1)]
sss = ''
TTmp2 = sss.join(TTmp)
TTmp2 = TTmp2.strip("\\")
AuxCalDate = np.append(AuxCalDate, [TTmp2], axis=0)
skip = st.unpack("H", f.read(2))
AuxNumCoeff = np.append(
AuxNumCoeff, st.unpack("H", f.read(2)), axis=0)
skip = [st.unpack("H", f.read(2)) for ui in range(1, 5+1)]
f.seek(PtrToThisAux + 80, 0) # ensures aligned
AuxGainLabel.append([])
AuxGainCoeff.append([])
AuxGainMin.append([])
AuxGainMax.append([])
for j in range(1, int(AuxNumGain[i-1]+1)):
AuxGainLabel[i-1].append([])
AuxGainCoeff[i-1].append([])
AuxGainMin[i-1].append([])
AuxGainMax[i-1].append([])
TTmp = [str(
st.unpack("s", f.read(1)))[3] for ui in range(
1, 4+1)]
sss = ''
TTmp2 = sss.join(TTmp)
TTmp2 = TTmp2.strip("\\")
AuxGainLabel[i-1][j-1] = np.append(
AuxGainLabel[i-1][j-1], [TTmp2], axis=0)
skip = [st.unpack("H", f.read(2)) for ui in range(1, 4+1)]
TTmp = [st.unpack("f", f.read(4))[0] for ui in range(
1, 5+1)]
AuxGainCoeff[i-1][j-1] = np.append(
AuxGainCoeff[i-1][j-1], TTmp, axis=0)
# Minimum value used in calibration data
AuxGainMin[i-1][j-1] = st.unpack("f", f.read(4))[0]
# Maximum value used in calibration data
AuxGainMax[i-1][j-1] = st.unpack("f", f.read(4))
skip = [st.unpack("f", f.read(4)) for ui in range(1, 10+1)]
# Now Jump to the Transducer Info
TransducerSerialNum = []
TransducerFrequency = []
TransducerRadius = []
TransducerBeamWidth = []
TransducerKt = []
for i in range(1, NumAbsTransducers+1):
# Move to the start of the ABS information
f.seek(PktStartPos + PtrToTransducerInfo * 2 + 200 * (i-1), 0)
TTmp = [str(
st.unpack("s", f.read(1)))[3] for ui in range(1, 20+1)]
sss = ''
TTmp2 = sss.join(TTmp)
TTmp2 = TTmp2.strip("\\")
# GV-10
TransducerSerialNum = np.append(
TransducerSerialNum, [TTmp2], axis=0)
skip = st.unpack("H", f.read(2))
# In Hz
TransducerFrequency = np.append(
TransducerFrequency, st.unpack("f", f.read(4)))
# In meters
TransducerRadius = np.append(
TransducerRadius, st.unpack("f", f.read(4)))
# In Degrees (3dB beamdidth, derived from
# acoustic beam pattern)
TransducerBeamWidth = np.append(
TransducerBeamWidth, st.unpack("f", f.read(4))[0])
skip = [st.unpack("f", f.read(4)) for ui in range(1, 4+1)]
# This is only if set in the personality
TransducerKt = np.append(
TransducerKt, st.unpack("f", f.read(4)))
# -------------------------------------------------------------------------
# Read in the Regime (Logger Set-Up) Information
# -------------------------------------------------------------------------
# Read regime Details
status, pktSize = findAQUAscat1000Packet(f, 21, file_size)
if 0 == status:
print('Regime Packet not found abort')
else:
PktStartPos = f.tell()
# Session Information
# Not interested in session start time at the moment%
SessionControl = [st.unpack("H", f.read(2))
for ui in range(1, 11+1)]
TTmp = [str(st.unpack("s", f.read(1)))[3] for ui in range(1, 32+1)]
sss = ''
TTmp2 = sss.join(TTmp)
TTmp2 = TTmp2.strip("\\")
SessionTitle = TTmp2
# The Aux channels are sampled at the PingRate
# divided by the auxPingDiv
AuxPingDiv, = st.unpack("H", f.read(2))
# The serial pressure + Temperature sampling derived as above
SerialPingDiv, = st.unpack("H", f.read(2))
Junk = [st.unpack("H", f.read(2)) for ui in range(1, 2+1)]
# THIS IS INCORRECTLY SAVED BY THE SYSTEM, SHOULD BE 0 or 8
NumAuxChans, = st.unpack("H", f.read(2))
NumAuxSamples, = st.unpack("I", f.read(4))
# if NumAuxSamples!=0: # This is a temp trick to correct--|
# NumAuxChans=8 # for NumAuxChans always being 0--|
# # It should be eliminated when corrected by Aquatec---|
Junk = [st.unpack("H", f.read(2)) for ui in range(1, 4+1)]
# This is the base profile rate before averaging
PingRate, = st.unpack("H", f.read(2))
# This is number of profiles collected prior to averaging
NumPings, = st.unpack("I", f.read(4))
# Number of enabled ABS Channels
NumAbsTimeSlots, = st.unpack("H", f.read(2))
Junk = [st.unpack("H", f.read(2)) for ui in range(1, 3+1)]
# These are the offsets into the Packet
PtrToAuxInfo, = st.unpack("H", f.read(2))
PtrToAbsInfo, = st.unpack("H", f.read(2))
# Calculate Aux Specific Information
# if AuxPingDiv=0 then no aux channels are enabled
# and NumAuxChannels =0
if 0 == AuxPingDiv:
AuxSampleRate = 0
AuxNumSamples = 0
NumAuxChans = 0
else:
AuxSampleRate = PingRate / AuxPingDiv
# GV-10
AuxNumSamples = np.ceil(NumPings / AuxPingDiv)
if NumAuxChans == 0:
NumAuxChans = 8
# Serial Pressure + Temperature information
if 0 == SerialPingDiv:
NumSerialSamples = 0
SerialSampleRate = 0
else:
NumSerialSamples = np.ceil(NumPings / SerialPingDiv)
SerialSampleRate = PingRate / SerialPingDiv
# Nothing useful in the channel
# Now read in the ABS Channel information
# Move to the start of the ABS information
f.seek(PktStartPos + PtrToAbsInfo * 2, 0)
AbsComplex = []
AbsAverage = []
AbsDecimation = []
AbsBinLengthMM = []
AbsBinLength = []
AbsTransducerName = []
AbsTransducerRadius = []
AbsTransducerBeamWidth = []
AbsTransducerKt = []
AbsTxFrequency = []
AbsRxFrequency = []
AbsTxPulseLength = []
AbsStartingGain = []
AbsTVG = []
AbsPowerLevelpc = []
AbsPowerLevel = []
AbsStartBin = []
AbsRxChan = []
AbsTxChan = []
AbsNumBins = []
AbsNumProfiles = []
AbsProfileRate = []
AbsBinRange = []
for j in range(1, NumAbsTimeSlots+1):
AbsBinRange.append([])
# TransducerRadius = np.append(
# TransducerRadius,st.unpack("f",f.read(4))) # In meters
Tmp, = st.unpack("H", f.read(2))
# For magnitude=0,complex=2
AbsComplex = np.append(AbsComplex, Tmp & 2)
# No of bursts averaged before saving
AbsAverage = np.append(AbsAverage, st.unpack("H", f.read(2)))
# Raw sampling rate along a profile is 19.2MHz, AbsDecimation i
AbsDecimation = np.append(
AbsDecimation, st.unpack("H", f.read(2)))
# Converts to bin size in mm based on speed 1500ms-1
AbsBinLengthMM = np.append(
AbsBinLengthMM, 1.25 * 2 ** AbsDecimation[j-1])
# Stored as time in seconds
AbsBinLength = AbsBinLengthMM[j-1] / 1500
# Using the Trasnducer ID copy the relevent data across
# Used to look up transducer information from personality:
TransducerId, = st.unpack("H", f.read(2))
TransducerId = TransducerId+1
AbsTransducerName = np.append(
AbsTransducerName, TransducerSerialNum[TransducerId-1])
# in m
AbsTransducerRadius = np.append(
AbsTransducerRadius, TransducerRadius[TransducerId-1])
# in degs
AbsTransducerBeamWidth = np.append(
AbsTransducerBeamWidth, TransducerBeamWidth[
TransducerId-1])
AbsTransducerKt = np.append(
AbsTransducerKt, TransducerKt[TransducerId-1])
# In Hz
AbsTxFrequency = np.append(
AbsTxFrequency, st.unpack("f", f.read(4)))
# In Hz
AbsRxFrequency = np.append(
AbsRxFrequency, st.unpack("f", f.read(4)))
# In Hz
AbsTxPulseLength = np.append(
AbsTxPulseLength, st.unpack("f", f.read(4)))
Junk = st.unpack("f", f.read(4))
# In dB with reference to default (built-in) Gain of system
AbsStartingGain = np.append(
AbsStartingGain, st.unpack("f", f.read(4)))
# In dB / bin where first bin has StartingGain (not used, =0)
AbsTVG = np.append(AbsTVG, st.unpack("f", f.read(4)))
powerlevel, = st.unpack("H", f.read(2))
# Power Level in % of Maximum
AbsPowerLevelpc = np.append(
AbsPowerLevelpc, 100 / 2 ** (2 * powerlevel))
# Power Level in dB relative to Maximum Power
AbsPowerLevel = np.append(
AbsPowerLevel, -20 * np.log10(2 ** powerlevel))
# Number of Bins from Start of Tx pulse before recording
AbsStartBin = np.append(AbsStartBin, st.unpack("H", f.read(2)))
# Number of Bins recorded
AbsNumBins = np.append(AbsNumBins, st.unpack("H", f.read(2)))
# Indicates which channel on AQUAscat used for TX
AbsRxChan = np.append(AbsRxChan, st.unpack("B", f.read(1)))
# Indicates which channel on AQUAscat used for RX
AbsTxChan = np.append(AbsTxChan, st.unpack("B", f.read(1)))
Junk = [st.unpack("H", f.read(2)) for ui in range(1, 12+1)]
# Calculate the number of profiles that should
# be recorded for this channel
AbsNumProfiles = np.append(
AbsNumProfiles, NumPings / AbsAverage[j-1])
# Calculate the stored profile rate
AbsProfileRate = np.append(
AbsProfileRate, PingRate / AbsAverage[j-1])
# in m --
AbsBinRange[j-1] = np.append(AbsBinRange[j-1], np.array(
[np.arange((AbsStartBin[j-1]), (
AbsStartBin[j-1] + AbsNumBins[
j-1]))]) * AbsBinLengthMM[j-1] / 1000)
# --------------------------------------------------------------------
# Now deal with reading in the data
# --------------------------------------------------------------------
# Allocate Memory for the Data
AuxData = np.zeros((np.int(AuxNumSamples), np.int(NumAuxChans)))
AbsData = np.zeros((np.int(
AbsNumBins[0]), np.int(
AbsNumProfiles[0]), np.int(
NumAbsTimeSlots)), dtype=np.float, order='C')
PressTempData = np.zeros((NumSerialSamples, 2))
AuxIndex = 0
SerIndex = 0
AbsIndex = np.zeros((NumAbsTimeSlots, 1))
f.seek(0, 0)
# Now Read in all the Data
while(f.tell() != file_size): # 0 == feof(fid))
status, pktType, pktSize = readNextAQUAscat1000Header(f, file_size)
if 1 == status:
if pktType == 22:
chan = st.unpack("H", f.read(2))[0]+1
# Increase the Index
AbsIndex[chan-1] = AbsIndex[chan-1] + 1
Tmp = [st.unpack(
"H", f.read(2))[0] for ui in np.arange(
1, AbsNumBins[chan-1]+1)]
newList = [x / 65536 for x in Tmp]
AbsData[0:np.int(AbsNumBins[0]), np.int(
AbsIndex[chan-1]-1), np.int(
chan-1)] = np.array(newList[:])
# Case 41: Data were saved as Complex values
# (normalize by 32768)
# case (41)
elif pktType == 41:
chan = st.unpack("H", f.read(2))[0] + 1
# Increase the Index
AbsIndex[chan-1] = AbsIndex[chan-1] + 1
Tmp = [st.unpack(
"h", f.read(2))[0] for ui in np.arange(
1, 2 * AbsNumBins[chan-1]+1)]
sss = [x / 32768 for x in Tmp]
AbsData[:, np.int(
AbsIndex[chan-1]), np.int(chan-1)] = np.array(sss)
# Case 46: Auxiliary Channel Data
# case (46)
elif pktType == 46:
# Gain settings
temp = st.unpack("H", f.read(2))[0]
AuxGain = [(temp & 3) + 1, ((temp >> 2) & 3)+1, (
(temp >> 4) & 3)+1, (
(temp >> 6) & 3)+1, 1, 1, 1, 1]
# Flags
Junk = st.unpack("H", f.read(2))[0]
AuxIndex = AuxIndex + 1
Tmp = [st.unpack(
"H", f.read(2))[0] for ui in np.arange(
1, NumAuxChans+1)]
if AuxIndex == AuxNumSamples+1:
AuxIndex = AuxIndex-1
pass
else:
AuxData[AuxIndex-1, 0:NumAuxChans] = Tmp[:]
# AuxData[AuxIndex-1,0:NumAuxChans] = Tmp[:]
# Case 55: External Pressure Channel for upgraded
# ABS system - details to be provided
# case (55)
elif pktType == 55:
chan = st.unpack("H", f.read(2))[0]+1
# Increase the Index
SerIndex = SerIndex + 1
Tmp = [st.unpack(
"f", f.read(4))[0] for ui in np.arange(1, 2+1)]
PressTempData[SerIndex-1, :] = Tmp[:]
else: # otherwise
f.seek(pktSize * 2, 1)
print('Summary of Data Imported')
for j in range(1, NumAbsTimeSlots+1):
print('Timeslot {} Total Profiles = {}'.format(j-1, AbsIndex[j-1]))
print('Total Aux Samples = {}'.format(AuxIndex))
for i in range(0, NumAuxChans):
if not AuxGainCoeff[:][i]:
coeff = np.array(AuxGainCoeff[:][0])
AuxGainCoeff[:][i].append(coeff[0] * 0)
# coeff = AuxGainCoeff(:,AuxGain(i),i)
# Matlab orders polynomial coeff opposite to Aquatec
coeff = np.array(AuxGainCoeff[:][i])
coeff = coeff.reshape(1, 5)
coeff = np.fliplr(coeff)
coeff = coeff.reshape(5)
AuxData[:, i] = np.polyval(coeff, AuxData[:, i])
# Need to apply calibration coefficients for the Aux Channels, now that the
# gain is fixed, auto gain is not supported. Therefore will use last value
if NumAuxChans == 0:
AuxData = []
f.close # Close the *.aqa data file
AbsMean = np.zeros((np.int(
AbsNumBins[0]), np.int(
NumAbsTimeSlots)), dtype=np.float, order='C')
for i in range(0, NumAbsTimeSlots):
# This provides the correct mean for short files
AbsMean[:, i] = (np.sum(AbsData[:, :, i], axis=1)) / AbsIndex[i]
FNAME = ['SessionTitle', 'BurstTime', 'BurstNumber',
'WakeSource', 'PingRate', 'NumPings',
'NumAuxChans', 'NumAuxSamples', 'AuxSampleRate',
'AuxNumSamples', 'AuxChannelName', 'AuxChannelUnit',
'AuxData', 'NumAbsTimeSlots', 'AbsNumProfiles', 'AbsComplex',
'AbsAverage', 'AbsProfileRate', 'AbsDecimation',
'AbsBinLengthMM', 'AbsBinLength', 'AbsStartBin', 'AbsNumBins',
'AbsTransducerName', 'AbsTransducerRadius',
'AbsTransducerBeamWidth', 'AbsTransducerKt',
'AbsTxFrequency', 'AbsRxFrequency', 'AbsTxPulseLength',
'AbsStartingGain', 'AbsTVG', 'AbsPowerLevel', 'AbsRxChan',
'AbsTxChan', 'AbsMean', 'AbsData',
'AbsBinRange', 'serial_num']
fOut = {}
for tui in range(0, len(FNAME)):
fOut.update({FNAME[tui]: eval(FNAME[tui])})
f.close()
return Parameters(fOut)