#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
UB-MES instruments (:mod:`hydrac.instruments.ub_mes`)
=========================================================
.. autoclass:: UB_MES
: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 math import fmod
from hydrac.util.parameters import Parameters
from hydrac.util.paramcontainer import ParamContainer
class SettingsAmbiguity(Exception):
pass
[docs]class UB_MES(Acoustic):
""" UB-MES instrument class.
Base class : :class:`hydrac.instruments.acoustique.Acoustic`
The UB_MES class reads the raw binary .udt files from UB-MES instruments
and stores the valuable information into the ``param`` modified dictionnary
with the common shape handled by hydrac
(see :mod:`hydrac.instruments.acoustique`).
.. image:: ../_image/ubmes.pdf
:width: 60%
:align: center
When loading the raw data files, the user is prompted with extra
informations. This instrument indeed works in bistatic configuration.
This configuration can be changed from one measurement to the other, while
keeping the same electronics. Each file thus must be completed with user
input distance (:math:`d_{0}`) from the emitter and reciever, as well as
bistatic angle between the principal axes of the emitter and receivers
(:math:`\gamma`) (see image below). This information will allow to
convert the travel times recorded by the instrument in terms of ranges
over the vertical.
.. image:: ../_image/bist.png
|
Bistatic instruments are more delicate to use than monostatic ones, as a
cell is not solely defined by half sound speed in water times the wave
travel time from the emitter to this cell (:math:`r = \\frac{ct}{2}`).
For bistatic instruments, a cell is defined by the region in space where
the travel time from the emitter to the cell and back to one receiver is
constant. It is an ellipsoïd surface centered on the emitter (E) - receiver
(R) (for one couple E-R) axis (:math:`ct_{j} = d_{j} + b_{j}`, see figure
above). The distance from the emitter to the cell j is thus given by :
.. math:: d_j = \\frac{d_0^2-(ct_j)^2}{2(d_0cos\\gamma-ct_j)}
The way the data are to be considered is similar to the example shown in
:mod:`hydrac.instruments.aquascat`, except the user will be prompted for
the geometric configuration of the sensors beforehand.
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 = 'ub_mes'
self.name = name
self.tag_M = 'mono'
self.nearfield_min = None
self.filepath, self.tempdir = self.file_select()
self.d0 = float(input('Distance between Emitter and receiver ' +
'(this value will be used for all selected' +
'files) :'))
self.gamma = float(input('Angle between Emitter/receiver axis and ' +
'vertical axis (this value will be used ' +
'for all selected files):'))
self.batch_read()
[docs] def _load_settings(self, currpath=None):
""" Data files from the UB_MES are two fold. One contains the data, the
other contains the settings of the instrument (ex. ping rate, pulse
length...). This function loads the settings from the setting file
associated to the udt file. The setting file is an xml file, that is
read here using the :class:`hydrac.util.paramcontainer.ParamContainer`
class. """
rr = glob.glob(self.tempdir + '/settings*.xml')
if (len(rr) != len(self.filepath) and len(rr) > 1) or rr == []:
raise SettingsAmbiguity('Number of settings files doesn t' +
'correspond to number of files')
elif (len(rr) != len(self.filepath) and len(rr) == 1):
import warnings
warnings.warn("Only one setting file is used for all the" +
"raw files")
pass
elif len(rr) == len(self.filepath):
_tag = [self.find_all(currpath, rr[i][-10:-4])
for i in range(0, len(rr))]
lp = [len(list(o)) for o in _tag]
from copy import copy
rr2 = copy(rr)
[rr2.remove(rr[i]) for i in range(0, len(lp)) if lp[i] == 0]
rr = rr2
try:
xml_set = Parameters(ParamContainer(path_file=rr[-1]))
xml_set.filepath = rr[-1]
except(IndexError, AttributeError):
xml_set = Parameters({})
return xml_set
[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()
filepath = list(filez)
filepath.sort()
filepath = tuple(filepath)
tempdir = os.path.dirname(os.path.abspath(filez[0]))
return filepath, tempdir
[docs] def batch_read(self):
"""Read a list of udt raw UB_MES file using the method :func:`read_UB`
"""
for kk in range(len(self.filepath)):
self.currfilepath = self.filepath[kk]
self.xml_set = self._load_settings(self.currfilepath)
self.xml_set.sequence.config.d0 = self.d0
self.xml_set.sequence.config.gamma = self.gamma
self.xml_set.sequence.config.gain_function.__rename__('a0', '_a0')
tmp = self.param_shape_UB(self.read_UB(self.filepath[kk]))
self.param.update({"P"+str(kk): tmp})
self.param['P' + str(kk)].ub_mes.update({'xml_set': self.xml_set})
self.param['P' + str(kk)].RawIntensity = 10 ** \
((70 - self.param['P' + str(kk)].ub_mes.xml_set.sequence.
config.gain_function._a0)/10) *\
self.param['P' + str(kk)].RawIntensity
# En bistatique, Kt est diff pour chaque porte. Implémenté ici
# directement sur les Raw int.
del self.currfilepath, self.xml_set
[docs] def find_all(self, a_str, sub):
start = 0
while True:
start = a_str.find(sub, start)
if start == -1:
return
yield start
start += len(sub) # use start += 1 to find overlapping matches
[docs] def _set_calib_ubmes(self):
"""see :func:`hydrac.instruments.aquascat.Aquascat._set_calib_aqa`.
In addition here, as working with a bistatic instrument, the Kt value
will change range as depending on the product of the transducers
beam patterns. The Kt is thus a refined array of Kts with respect to
the range from the emitter, proper to a bistatic configuration. This
vector can be interpolated at the ranges determined by the current
settings (that can evolve given the pulse length, speed of sound in
water...)
"""
rr = glob.glob(constants.dir_ubmes_calib + '*.xml')
rr.sort()
try:
ub = Parameters(ParamContainer(path_file=rr[-1]))
except(IndexError, AttributeError):
ub = Parameters({})
return ub
[docs] def _update_calib_ubmes(self, a, b=None, remove_elem=False):
""" see :func:`hydrac.instruments.aquascat.Aquascat._update_calib_aqa`
"""
rr = glob.glob(constants.dir_ubmes_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_ubmes_calib +
'calib_UB.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='CalibUB')
for a, b in uiui.items():
yoyo._set_child(a, b)
yoyo._save_as_xml(path_file=constants.dir_ubmes_calib +
'calib_UB.xml', find_new_name=False,
overwrite=True)
[docs] def param_shape_UB(self, inp):
"""Reshapes the modified dictionnary ``param`` to fit the standard key
names and data shape used throughout hydrac"""
tmp = self.xml_set.sequence.config
inp2 = Parameters({'ub_mes': {}})
inp2.BurstTime = inp.time[0]*tmp.n_ech
inp2.PingRate = tmp.prf
inp2.NumPingsTot = np.size(inp.u, 0) * tmp.n_ech
inp2.NumChannels = np.size(inp.RawIntensity, 2)
inp2.NumProfilesSamples = np.ones((inp2.NumChannels)) *\
np.size(inp.u, 0)
inp2.NumAverage = np.ones((inp2.NumChannels)) * tmp.n_ech
inp2.SampleRate = np.ones((inp2.NumChannels)) * tmp.prf / tmp.n_ech
inp2.BinLengthMM_bist = 1000 * \
np.tile(np.concatenate((np.diff(inp.r.T),
np.diff(inp.r[-2:].T)),
axis=1).T, (1, inp2.NumChannels))
inp2.BinLengthMM = 1000 * \
np.ones((inp2.NumChannels)) * \
np.nanmean(np.concatenate((np.diff(inp.r.T),
np.diff(inp.r[-2:].T)), axis=1).T)
inp2.StartBin = np.ones((inp2.NumChannels)) * tmp.r_vol1
inp2.NumBins = np.ones((inp2.NumChannels)) * np.size(inp.u, 1)
inp2.TransducerTag = np.array([1, 2])
inp2.TransducerRadius = np.array([0, 0])
inp2.TransducerBeamWidth = np.array([0, 0])
inp2.Frequency = np.ones((inp2.NumChannels)) * tmp.f0 * 1e6
inp2.FrequencyRx = np.ones((inp2.NumChannels)) * tmp.f0 * 1e6
inp2.TransducerKt = [[a, b] for a, b in zip(list(inp2.FrequencyRx.T),
list(np.sqrt(3 / (16 *
np.pi)) +
np.zeros((inp2.NumChannels
)).T))]
inp2.TransducerName = np.tile(np.array(str(inp2.Frequency) + ' MHz',
dtype=str),
(inp2.NumChannels))
inp2.PulseLength = np.ones((inp2.NumChannels)) * 2 / (tmp.f0 * 1e6)
inp2.RxGain = np.array([0, 0])
inp2.TxGain = tmp.gain_function._a0
inp2.TVG = np.tile(np.zeros((np.size(inp.u, 1), 1)),
(1, inp2.NumChannels))
inp2.MeanProfile = np.nanmean(inp.RawIntensity, 0)
inp2.RawIntensity = inp.RawIntensity
inp2.BinRange = np.tile(inp.r, (1, inp2.NumChannels))
# inp2.BinRange_bj = np.tile(inp.bj, (1, inp2.NumChannels))
inp2.SerialNum = ''
inp2.time = inp.time
inp2.pressure = np.tile(np.zeros((np.size(inp.u, 0), 1)),
(1, inp2.NumChannels))
inp2.depth_raw = np.tile(np.zeros((np.size(inp.u, 0), 1)),
(1, inp2.NumChannels))
inp2.depth = np.zeros((np.size(inp.u, 0)))
inp2.temp = np.zeros((np.size(inp.u, 0))) * np.nan
inp2.battery = np.zeros((np.size(inp.u, 0))) * np.nan
inp2.TransducerAngle = np.array([90-tmp.gamma * 180 /
np.pi, -(90 - tmp.gamma
* 180 / np.pi)])
inp2.sal = np.zeros((np.size(inp.u, 0))) * np.nan
inp2.NumSamples = np.ones((inp2.NumChannels)) * np.size(inp.u, 0)
inp2.FileName = os.path.split(self.currfilepath)[-1]
inp2.ub_mes = Parameters({'AuxData': {'u': inp.u, 'w': inp.w},
'alpha': inp.alpha,
'AuxNumSamples': inp2.NumSamples,
'AuxSampleRate': inp2.SampleRate,
'time_aux': inp2.time,
'BinRange_bj': np.tile(inp.bj,
(1, inp2.NumChannels))
})
_full_calib = self._set_calib_ubmes()
_full_calib = _full_calib.SN1
inp2.calib = Parameters({})
inp2.calib.kton, inp2.calib.ktoff = \
self.kt_decimation(_full_calib, inp2.BinRange[:, 0])
test=np.transpose(np.tile(np.atleast_3d(np.concatenate((np.atleast_2d(inp2.calib.kton),np.atleast_2d(inp2.calib.ktoff)),axis=0).T),(1,1,np.int(inp2.NumSamples[0]))),(0,2,1))
inp2.RawIntensity = inp2.RawIntensity/test
return inp2
[docs] def kt_decimation(self, a, b):
"""Interpolates the Kts relative to the current bistatic config to
the present range positions determined according to the settings
xml file."""
from scipy import interpolate
f1 = interpolate.interp1d(np.array(a.dj), np.array(a.kt2on))
f2 = interpolate.interp1d(np.array(a.dj), np.array(a.kt2off))
return np.abs(f1(b)), np.abs(f2(b))
[docs] def read_UB(self, fname):
"""Read one .udt 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')
header = [str(st.unpack("s", f.read(1)))[3] for ui in range(1, 6+1)]
header = ''.join(header)
print('header record file: {}'.format(header))
if header == 'UDT003' or header == 'UDT002':
data_format = 2
header2 = [str(st.unpack("s", f.read(1)))[3]
for ui in range(1, 36+1)]
header2 = ''.join(header2)
print('header2 : {}'.format(header2))
else:
data_format = 1
ind = 0
i_config = 0
i_prof = 0
loopy = 1
u_rad = []
intensity = []
t = []
t2 = []
while loopy == 1:
ind += 1
try:
profile = self.find_profile(f, 1, ind, i_prof, data_format)
u_rad.append(profile.velocity)
intensity.append(profile.amplitude ** 2)
t.append(profile.time_stamp)
t2.append(profile.t)
except EOFError:
print('supposed end of file')
loopy = 0
u1u2 = np.array(u_rad)
intensity = np.array(intensity)
u, w, r, alpha, bj = self.transform(u1u2)
time = np.array(t2)
BurstTime = datetime.datetime.fromtimestamp(time[0])
rawIntensity = np.transpose(np.array([intensity[:, 1::2],
intensity[:, 0::2]]), (2, 1, 0))
return Parameters({'BurstTime': BurstTime,
'RawIntensity': rawIntensity,
'u': u, 'w': w, 'time': time,
'r': r, 'alpha': alpha, 'bj': bj})
def find_profile(self, fd, config_num,
select_profile, i_prof, data_format):
loop = 1
total_size = 0
profile = []
tcp_flag = 20300
while loop:
i1 = fd.tell()
[flag, size, data, total_size] = self.read_chunck(fd, total_size)
if len(data) == 0:
loop = 0
i2 = fd.tell()
fd.seek(i1+8, 0)
if flag == tcp_flag:
[ref_config, num_config, profile] = \
self.extract_generic_profile(size, fd, data_format)
if num_config != config_num:
print('num_config = {}'.format(str(num_config)))
print('config_num = {}'.format(str(config_num)))
else:
i_prof += 1
loop = 0
fd.seek(i2, 0)
return profile
def _read_file(self, fd, _size, total_size):
_data = fd.read(_size)
if _data == b'':
print("{} byte read in the file".format(total_size))
raise EOFError
else:
if _size != len(_data):
print("WARNING, not all read")
total_size += _size
return _data
def read_chunck(self, fd, total_size):
flag = st.unpack('i', self._read_file(fd, st.calcsize('i'),
total_size))[0]
total_size += 4
size = st.unpack('i', self._read_file(fd, st.calcsize('i'),
total_size))[0]
total_size += 4
if size:
data = self._read_file(fd, size, total_size)
else:
print("chunck empty")
data = ''
crc = st.unpack('i', self._read_file(fd, st.calcsize('i'),
total_size))[0]
total_size += 4
return flag, size, data, total_size
def extract_generic_profile(self, size, fd, version):
if version == 2:
head_size = 20
ref, = st.unpack("i", fd.read(4))
size_data, = st.unpack("I", fd.read(4))
size_raw, = st.unpack("I", fd.read(4))
sec, = st.unpack("i", fd.read(4))
n_sec, = st.unpack("i", fd.read(4))
if size_data == 0 and size_raw == 0:
print('caca')
elif version == 1:
head_size = 24
n_vol, = st.unpack("i", fd.read(4))
ref, = st.unpack("i", fd.read(4))
sec, = st.unpack("i", fd.read(4))
n_sec, = st.unpack("i", fd.read(4))
origine, = st.unpack("f", fd.read(4))
interval, = st.unpack("f", fd.read(4))
if n_vol == 0:
print('caca')
else:
print('unknown data version')
ref_config = ref & 0xFFFFFF00 >> 8
num_config = ref & 0x000000FF
decimal1 = fmod(sec, 1.0)
integer1 = sec - decimal1
decimal1 = decimal1 * 1e9
decimal2 = fmod(n_sec + decimal1, 1.0)
integer2 = n_sec + decimal1 - decimal2
decimal2 = decimal2 * 1e9
t = sec + n_sec * 1e-9
sec = integer1 + integer2
if sec > 0 and decimal2 < 0:
sec = sec-1
decimal2 = decimal2 + 1000000000
if sec < 0 and decimal2 > 0:
sec = sec + 1
decimal2 = decimal2 - 1000000000
profile = Parameters({'t': t, 'time_stamp': sec, 'ref': ref,
'n_vol_total': 0, 'position': [], 'velocity': [],
'amplitude': [], 'snr': [], 'sat': [],
'ny_jp': [], 'sigma': [], 'n_ech': 0,
'I': [], 'Q': [], 'fft_doppler': [],
'freqs': [], 'raw_data': []})
if version == 2:
if size_data:
# floating + integer arrays
tab_size = size_data / 3
offset = head_size
velocity_b = []
[velocity_b.append(st.unpack("f", fd.read(4))[0])
for az in range(0, int(len(np.arange(offset, offset +
tab_size)) / 4))]
profile.velocity = velocity_b
offset = offset + tab_size
variance_b = []
[variance_b.append(st.unpack("f", fd.read(4))[0])
for az in range(0, int(len(np.arange(offset, offset +
tab_size)) / 4))]
profile.amplitude = np.sqrt(variance_b)
offset = offset + tab_size
facteur_qualite_b = []
[facteur_qualite_b.append(st.unpack("i", fd.read(4))[0])
for az in range(0, int(len(np.arange(offset, offset +
tab_size)) / 4))]
profile.facteur_qualite = facteur_qualite_b
offset = offset + tab_size
if profile.n_vol_total == 0:
profile.n_vol_total = len(profile.velocity)
elif profile.n_vol_total != len(profile.velocity):
print(['Error set_data_from_bin: Different length ' +
'between arrays, cant process !!! n_vol_tot' +
'= {}'.format(str(len(profile.velocity)))])
return
for i in range(profile.n_vol_total):
# 2013/03/04 Damien : plage SNR : [-10dB +10dB]
profile.snr.append((float(profile.facteur_qualite[i] &
0x000000FF) * 20.0 / 255.0)-10.0)
profile.sat.append(float((profile.facteur_qualite[i] >>
8) & 0x000000FF))
profile.ny_jp.append(float((profile.facteur_qualite[i] >>
16) & 0x000000FF))
profile.sigma.append(float((profile.facteur_qualite[i] >>
24) & 0x000000FF))
if size_raw + offset > size:
print('raw data size exceed size of data')
# raw can be anything, not always the IQ data,
# so not automatically extract IQ
if offset < size:
print(['Warning, extract generic data: still {} data' +
'to read !!'.format(str(size - offset))])
elif version == 1:
pass
return ref_config, num_config, profile