#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
EK60 instruments (:mod:`hydrac.instruments.ek60`)
=========================================================
.. autoclass:: EK60
:members:
:private-members:
"""
# from .instruments import instrument
from .acoustique import Acoustic
import numpy as np
import struct as st
import glob
import os
import datetime as dt
from tkinter import filedialog
from tkinter import Tk
from hydrac.util.parameters import Parameters
import sys
from scipy import interpolate
import scipy.io as sio
class InstrAmbiguity(Exception):
pass
[docs]class EK60(Acoustic):
""" EK60 instrument class.
Base class : :mod:`hydrac.instruments.acoustique.Acoustic`
The EK60 is a scientific mulitfrequency split beam echosounder developped
by SIMRAD (See https://www.simrad.com/ek60). It is a multifrequency
(frequency range from ~kHz) echosounder deployed from marine platforms such
as reasearch vessels. It is used for fish stock assessments and zooplankton
biomass estimations.
.. image:: ../_image/ek60.pdf
:width: 80%
:align: center
The EK60 class reads the raw binary files and matlab files generated by the
echopen software,and stores the valuable information into the ``param``
modified dictionnary with the common shape handled by hydrac
(see :mod:`hydrac.instruments.acoustique`).
The way the data are to be considered is similar to the example shown in
:mod:`hydrac.instruments.aquascat`.
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 = 'ek_60'
self.name = name
self.tag_M = 'multi'
self.nearfield_min = None
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()
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 .raw raw EK60 file using the methods:
- Raw files : :func:`read_EK60`
- Mat files : :func:`load_from_mat`
"""
for kk in range(len(self.filepath)):
self.currfilepath = self.filepath[kk]
if self.currfilepath[-3:] in ['raw', 'RAW']:
tmp = self.param_shape_EK(*self.read_EK60(self.currfilepath))
elif self.currfilepath[-3:] in ['mat', 'MAT']:
tmp = self.load_from_mat(self.currfilepath)
self.param.update({"P" + str(kk): tmp})
del self.currfilepath
[docs] def load_from_mat(self, fname):
"""Read one .mat EK60 file"""
mc2 = sio.loadmat(fname)
m_keys = []
[m_keys.append(i) for i in mc2.keys() if '__' not in i]
ptmp = self.param_def(mc2, m_keys)
ptmp2 = Parameters({})
pings = ptmp.data[0].pings
ptmp2.NumChannels = np.shape(ptmp.calParms)[0]
sampleinterval =\
np.array([np.squeeze(ptmp.calParms[i].sampleinterval[0])
for i in range(ptmp2.NumChannels)])
soundvelocity =\
np.array([np.squeeze(ptmp.calParms[i].soundvelocity[0])
for i in range(ptmp2.NumChannels)])
BinLengthMM =\
1000*np.array([sampleinterval[i]*soundvelocity[i]/2
for i in range(ptmp2.NumChannels)])
alpha = np.array([np.squeeze(ptmp.calParms[i].absorptioncoefficient[0])
for i in range(ptmp2.NumChannels)])
angleoffsetalong =\
np.array([np.squeeze(ptmp.calParms[i].anglesoffsetalongship[0])
for i in range(ptmp2.NumChannels)])
angleoffsetathwart =\
np.array([np.squeeze(ptmp.calParms[i].angleoffsetathwartship[0])
for i in range(ptmp2.NumChannels)])
psi = np.array([np.squeeze(ptmp.calParms[i].equivalentbeamangle[0])
for i in range(ptmp2.NumChannels)])
G = np.array([np.squeeze(ptmp.calParms[i].gain[0])
for i in range(ptmp2.NumChannels)])
pt = np.array([np.squeeze(ptmp.calParms[i].transmitpower[0])
for i in range(ptmp2.NumChannels)])
transducerdepth =\
np.array([np.squeeze(ptmp.calParms[i].transducerdepth[0])
for i in range(ptmp2.NumChannels)])
size_num = [pings[i].samplerange[1] for i in range(ptmp2.NumChannels)]
size_orig = np.array(size_num)
maxsize = max(size_num)
size_num.remove(maxsize)
pings = [self.scale_data(pings[i], size_num, maxsize)
for i in range(len(pings))]
dR = np.min([a * b for a, b in zip(soundvelocity, sampleinterval)])/2
r = np.arange(maxsize) * dR - 2 * dR
r[np.where(r <= 0)] = np.spacing(1)
ptmp2.BinRange = np.tile(r, (ptmp2.NumChannels, 1)).T
time = np.squeeze(pings[0].time)
ptmp2.BinLengthMM = 1000 *\
np.array([pings[i].sampleinterval[i]*pings[i].soundvelocity[i]/2
for i in range(len(pings))])
ptmp2.BurstTime =\
dt.datetime.timestamp(dt.datetime.fromordinal(int(time[0])-366) +
dt.timedelta(days=time[0] % 1))
ptmp2.time =\
np.array([dt.datetime.timestamp(dt.datetime.fromordinal(int(time[i]
)-366)
+ dt.timedelta(days=time[i] % 1))
for i in range(len(time))])
ptmp2.Frequency =\
np.array([np.squeeze(float(ptmp.calParms[i].frequency[0]))
for i in range(ptmp2.NumChannels)])
ptmp2.FrequencyRx =\
np.array([np.squeeze(float(ptmp.calParms[i].frequency[0]))
for i in range(ptmp2.NumChannels)])
ptmp2.NumAverage = np.ones(ptmp2.NumChannels)
ptmp2.NumBins = np.zeros(ptmp2.NumChannels) + maxsize
ptmp2.NumPingTot = maxsize
ptmp2.NumProfileSamples = np.zeros(ptmp2.NumChannels) + len(time)
ptmp2.NumSamples = np.zeros(ptmp2.NumChannels) + len(time)
ptmp2.PingRate = 1 / np.nanmean(np.diff(time))
ptmp2.PulseLength =\
np.array([np.squeeze(ptmp.calParms[i].pulselength[0])
for i in range(ptmp2.NumChannels)])
ptmp2.RxGain = np.zeros(ptmp2.NumChannels)
ptmp2.SampleRate = np.ones(ptmp2.NumChannels)
ptmp2.StartBin = np.ones(ptmp2.NumChannels) * r[0]
ptmp2.TVG = np.dot(np.atleast_2d(r).T, 2*np.atleast_2d(alpha)) +\
np.tile(np.atleast_2d(r).T, (1, ptmp2.NumChannels))
ptmp2.TransducerAngle = np.zeros(ptmp2.NumChannels)
ptmp2.TransducerAngle_ = [(a, b) for a, b
in zip(angleoffsetalong, angleoffsetathwart)]
ptmp2.TransducerBeamWidth =\
2 * np.arcsin(np.sqrt((10**(psi / 10))*2.5 / (4 * np.pi)
)) * 180 / np.pi
ptmp2.TransducerKt = []
[ptmp2.TransducerKt.append([a, b])
for a, b in zip(ptmp2.FrequencyRx, G)]
ptmp2.TransducerName = [str(d) + ' MHz' for d in ptmp2.Frequency]
ptmp2.TransducerRadius = np.zeros(np.shape(ptmp2.Frequency)) + np.NaN
ptmp2.TransducerTag = np.arange(0, ptmp2.NumChannels)
ptmp2.TransducerTag2 = np.arange(0, ptmp2.NumChannels)
ptmp2.TxGain = pt
ptmp2.battery = np.zeros(np.shape(ptmp2.Frequency)) + np.NaN
ptmp2.depth_raw = np.tile(transducerdepth, (len(time), 1))
ptmp2.depth = ptmp2.depth_raw[:, 0]
ptmp2.pressure = np.ones(len(ptmp2.depth))*np.nan
ptmp2.sal = np.ones(len(ptmp2.depth))*np.nan
temp = np.ones(len(ptmp2.depth)) * np.nan
ptmp2.Rawintensity = np.zeros((maxsize,
len(ptmp2.time),
len(ptmp2.Frequency)))
ptmp2.Intensity = np.zeros((maxsize, len(ptmp2.time),
len(ptmp2.Frequency)))
ptmp2.asv = np.zeros((maxsize,
len(ptmp2.time),
len(ptmp2.Frequency)))
ptmp2.FileName = os.path.split(self.currfilepath)[-1]
for i in range(len(pings)):
ptmp2.asv[:, :, i] = 10**(pings[i].Sv / 10)
ptmp2.Intensity[:, :, i] = 10**(pings[i].Sv / 10)
ptmp2.MeanProfile = np.nanmean(ptmp2.Intensity, 1)
ptmp2.calib = Parameters(ptmp.calParms)
ptmp2.ek_60 = Parameters({'AuxData': {}})
[ptmp2.ek_60.AuxData.update({k: Parameters(ptmp.data[0][k])})
for k in ptmp.data[0].keys()]
ptmp2.ek_60.update(Parameters({'size_orig': size_orig,
'fech_orig': 1/sampleinterval
}))
return ptmp2
[docs] def param_def(self, o, m_keys, p=None):
"""
Spans through the mat file to extract the basic and nested matlab
structures and read them in python language.
Parameters
----------
o : scipy.io.loadmat object
Raw matlab structure loaded
m_keys : str
list of top level structure names in the mat object.
"""
if p is None:
P = Parameters({})
else:
P = p
for u in m_keys:
try:
cdex = 1
siz1 = o[u].shape[0]
siz2 = o[u].shape[1]
except IndexError:
cdex = 0
siz1 = 1
siz2 = 1
if 'O' not in o[u].dtype.descr[0][1] or\
len(o[u].dtype.descr[0][0]) == 0:
siz1 = 1
siz2 = 1
cdex = 0
P[u] = Parameters([{}])
if siz1 * siz2 == siz1 or siz1 * siz2 == siz2:
if siz1 * siz2 == 1:
P[u] = Parameters([{}])
else:
P[u] = Parameters([{}] * max([siz1, siz2]))
else:
P[u] = Parameters([[{}] * siz1] * siz2)
for u1 in range(siz1):
for u2 in range(siz2):
if cdex == 0:
try:
tmp = np.atleast_1d(np.squeeze(o[u]))
except IndexError:
tmp = np.zeros(1)
else:
tmp = o[u][u1, u2]
tmp_descr = [tmp.dtype.descr[i][0]
for i in range(len(tmp.dtype.descr))]
if len(tmp_descr) > 1 and len(tmp_descr[0]) != 0:
self.param_def(tmp, tmp_descr, p=P[u][max([u1, u2])])
else:
P[u] = tmp
return P
[docs] def scale_data(self, pings, size_num, maxsize):
"""
As frequencies can range from a fex kHz to a few hundreds, the cell
size can greatly differ from one frequency to the other (the higher the
frequency, the finer the vertical resolution).
This method harmonizes the HAC profiles sizes from all transducer
channels according to the longest profile size by interpolation using a
nearest neighbour method.
Parameters
----------
pings : Paramcontainer object :class:`hydrac.util.parameters.AttrDict`
Dictionnary like object containing the intensity profiles at each
frequency
size_num : float
list of profiles sizes per frequency channel.
maxsize : float
Target size for each profile length
"""
any_in = lambda a, b: bool(set(a).intersection(b))
for k in pings.keys():
if any_in(np.shape(pings[k]), size_num):
u = list(np.shape(np.atleast_2d(pings[k])))
x = list(set(np.shape(pings[k])).intersection(set(size_num))
)[0]
u.remove(list(set(np.shape(pings[k])
).intersection(set(size_num)))[0])
y = u[0]
lp = np.zeros((maxsize, y))
ol = np.arange(x)
f1 = interpolate.interp1d(ol, pings[k].T, kind='nearest')
pings[k] = f1(np.linspace(ol[0], ol[-1], num=maxsize)).T
return pings
[docs] def param_shape_EK(self, inp0, inp):
"""Reshapes the modified dictionnary ``param`` to fit the standard key
names and data shape used throughout hydrac"""
calParms = self.readEKRaw_GetCalParms(inp0, inp)
inp2 = Parameters({})
inp2.BinLengthMM = 1000 *\
np.array([inp.pings[i].sampleinterval[i] *
inp.pings[i].soundvelocity[i]/2
for i in range(len(inp.pings))])
size_num = [inp.pings[i].samplerange[1] for i in range(len(inp.pings))]
size_orig = np.array(size_num)
maxsize = max(size_num)
size_num.remove(maxsize)
inp.pings = [self.scale_data(inp.pings[i], size_num, maxsize)
for i in range(len(inp.pings))]
t = np.array([calParms[n].sampleinterval
for n in range(len(inp.pings))])
c = np.array([calParms[n].soundvelocity
for n in range(len(inp.pings))])
c0 = c[0]
pt = np.array([calParms[n].transmitpower
for n in range(len(inp.pings))])
tau = np.array([calParms[n].pulselength
for n in range(len(inp.pings))])
G = np.array([calParms[n].gain
for n in range(len(inp.pings))])
psi = np.array([calParms[n].equivalentbeamangle
for n in range(len(inp.pings))])
alpha = np.array([calParms[n].absorptioncoefficient
for n in range(len(inp.pings))])
sac = 2*np.squeeze(np.array([[calParms[n].sacorrectiontable[j]
for j in range(len(calParms[n].pulselengthtable))
if calParms[n].pulselengthtable[j]==tau[n]]
for n in range(len(inp.pings))]))
angleoffsetalong = np.array([calParms[n].angleoffsetalongship
for n in range(len(inp.pings))])
angleoffsetathwart = np.array([calParms[n].angleoffsetathwartship
for n in range(len(inp.pings))])
transducerdepth = np.array([calParms[n].transducerdepth
for n in range(len(inp.pings))])
dR = np.min([a * b for a, b in zip(c, t)]) / 2
r = np.arange(maxsize) * dR - 2 * dR
r[np.where(r <= 0)] = np.spacing(1)
inp2.BinRange = np.tile(r, (len(inp.pings), 1)).T
inp2.BurstTime = inp0.time
inp2.time = inp.pings[0].time
inp2.FileName = os.path.split(self.currfilepath)[-1]
inp2.Frequency = np.array([calParms[n].frequency
for n in range(len(inp.pings))])
inp2.FrequencyRx = np.array([calParms[n].frequency
for n in range(len(inp.pings))])
csv = 10 * np.log10(pt * (c0 / inp2.Frequency)**2 * c0 * tau *
10**(psi / 10) / (32 * np.pi**2))
inp2.RawIntensity = np.transpose(np.array([inp.pings[n].power
for n in
range(len(inp.pings))]),
(1, 2, 0)) - np.tile(csv + sac,
(len(r),
len(inp2.time),
1))
inp2.RawIntensity = 10**(inp2.RawIntensity / 10)
inp2.MeanProfile = np.nanmean(inp2.RawIntensity, 1)
inp2.NumAverage = np.ones(len(inp.pings))
inp2.NumBins = np.ones(len(inp.pings)) * maxsize
inp2.NumChannels = len(inp.pings)
inp2.NumPingTot = maxsize
inp2.NumProfileSamples = np.ones(len(inp.pings)) *\
np.float(len(inp2.time))
inp2.NumSamples = np.ones(len(inp.pings)) * np.float(len(inp2.time))
inp2.PingRate = 1 / np.nanmean(np.diff(inp2.time))
inp2.PulseLength = tau
inp2.RxGain = np.zeros(len(inp.pings))
inp2.SampleRate = np.ones(len(inp.pings))
inp2.StartBin = np.ones(len(inp.pings)) * r[0]
inp2.TVG = np.dot(np.atleast_2d(r).T, 2 * np.atleast_2d(alpha)) +\
np.tile(np.atleast_2d(r).T, (1, len(inp.pings)))
inp2.TransducerAngle_ = [(a, b) for a, b in
zip(angleoffsetalong, angleoffsetathwart)]
inp2.TransducerAngle = np.zeros((len(inp2.Frequency)))
inp2.TransducerBeamWidth =\
2 * np.arcsin(np.sqrt((10**(psi / 10)) * 2.5 / (4 * np.pi))
) * 180/np.pi
inp2.TransducerKt = []
[inp2.TransducerKt.append([a, b]) for a, b in zip(inp2.FrequencyRx, G)]
inp2.TransducerName = [str(d) + ' MHz' for d in inp2.Frequency]
inp2.TransducerRadius = np.zeros(np.shape(inp2.Frequency)) + np.NaN
inp2.TransducerTag = np.arange(0, len(inp.pings))
inp2.TransducerTag2 = np.arange(0, len(inp.pings))
inp2.TxGain = pt
inp2.battery = np.zeros(np.shape(inp2.Frequency)) + np.NaN
inp2.depth_raw = np.tile(transducerdepth, (len(inp2.time), 1))
inp2.depth = inp2.depth_raw[:, 0]
inp2.pressure = inp2.depth_raw[:, 0] * np.nan
inp2.sal = inp2.depth_raw[:, 0] * np.nan
inp2.temp = np.ones(len(inp2.depth)) * np.nan
inp2.ek_60 = Parameters({'AuxData': {'gps': inp.gps,
'vlog': inp.vlog,
'vspeed': inp.vspeed},
'size_orig': size_orig,
'fech_orig': 1/t
})
inp2.calib = Parameters(calParms)
return inp2
def readEKRaw_ReadDgHeader(self, fid, timeOffset):
# read datagram type
dgType = [str(st.unpack("<s", fid.read(1)))[3] for ui in range(1, 4+1)]
# read datagram time (NT Time - number of 100-nanosecond
# intervals since January 1, 1601)
lowdatetime, = st.unpack("<I", fid.read(4))
highdatetime, = st.unpack("<I", fid.read(4))
# convert NT time to MATLAB serial time
ntSecs = (highdatetime * 2 ** 32 + lowdatetime) / 10000000
ntSecs_ent = int(np.fix(ntSecs))
ntSecs_micr = int(np.round(1e6*(ntSecs-ntSecs_ent)))
_dgTime = dt.datetime(1601, 1, 1) +\
dt.timedelta(seconds=ntSecs_ent, microseconds=ntSecs_micr)
dgTime = _dgTime.timestamp()
return dgType, dgTime
def readEKRaw_ReadConfigHeader(self, fid):
configheader = Parameters({})
configheader.surveyname = [str(st.unpack("<s", fid.read(1)))[3]
for ui in range(1, 128+1)]
configheader.transectname = [str(st.unpack("<s", fid.read(1)))[3]
for ui in range(1, 128+1)]
configheader.soundername = [str(st.unpack("<s", fid.read(1)))[3]
for ui in range(1, 128+1)]
configheader.spare = [str(st.unpack("<s", fid.read(1)))[3]
for ui in range(1, 128+1)]
configheader.transceivercount, = st.unpack("<i", fid.read(4))
return configheader
def readEKRaw_ReadTransceiverConfig(self, fid):
configXcvr = Parameters({})
configXcvr.channelid = [str(st.unpack("<s", fid.read(1)))[3]
for ui in range(1, 128+1)]
configXcvr.beamtype, = st.unpack("<i", fid.read(4))
configXcvr.frequency, = st.unpack("<f", fid.read(4))
configXcvr.gain, = st.unpack("<f", fid.read(4))
configXcvr.equivalentbeamangle, = st.unpack("<f", fid.read(4))
configXcvr.beamwidthalongship, = st.unpack("<f", fid.read(4))
configXcvr.beamwidthathwartship, = st.unpack("<f", fid.read(4))
configXcvr.anglesensitivityalongship, = st.unpack("<f", fid.read(4))
configXcvr.anglesensitivityathwartship, = st.unpack("<f", fid.read(4))
configXcvr.angleoffsetalongship, = st.unpack("<f", fid.read(4))
configXcvr.angleoffsetathwartship, = st.unpack("<f", fid.read(4))
configXcvr.posx, = st.unpack("<f", fid.read(4))
configXcvr.posy, = st.unpack("<f", fid.read(4))
configXcvr.posz, = st.unpack("<f", fid.read(4))
configXcvr.dirx, = st.unpack("<f", fid.read(4))
configXcvr.diry, = st.unpack("<f", fid.read(4))
configXcvr.dirz, = st.unpack("<f", fid.read(4))
configXcvr.pulselengthtable = [st.unpack("<f", fid.read(4))[0]
for ui in range(1, 5+1)]
configXcvr.spare2 = [str(st.unpack("<s", fid.read(1)))[3]
for ui in range(1, 8+1)]
configXcvr.gaintable = [st.unpack("<f", fid.read(4))[0]
for ui in range(1, 5+1)]
configXcvr.spare3 = [str(st.unpack("<s", fid.read(1)))[3]
for ui in range(1, 8+1)]
configXcvr.sacorrectiontable = [st.unpack("<f", fid.read(4))[0]
for ui in range(1, 5+1)]
configXcvr.spare4 = [str(st.unpack("<s", fid.read(1)))[3]
for ui in range(1, 52+1)]
return configXcvr
def readEKRaw_ReadHeader(self, fid):
tmp, = st.unpack("<i", fid.read(4))
dgType, dgTime = self.readEKRaw_ReadDgHeader(fid, 0)
configheader = self.readEKRaw_ReadConfigHeader(fid)
configheader.time = dgTime
# extract individual xcvr configurations and store list of frequencies
frequencies = np.zeros((1, configheader.transceivercount))
configXcvr = []
frequencies = []
for i in range(configheader.transceivercount):
configXcvr.append(self.readEKRaw_ReadTransceiverConfig(fid))
frequencies.append(configXcvr[i].frequency)
# create the configuration structure - store header and xcvr configs
header = Parameters({})
header.header = configheader
header.transceiver = configXcvr
st.unpack("<i", fid.read(4))
return header, frequencies
def readEKRaw_GetSampleCount(self, fid, CIDs, sampleRange, maxSampleRange,
ping):
# Raw datagram header size
HEADER_LEN = 12
# get position in file
fPosition = fid.tell()
# determine number of freqs and allocate return array
nXcvrs = len(CIDs)
nSamples = np.zeros((nXcvrs, 1))
# create array that tracks number of pings scanned
nPings = np.zeros((nXcvrs, 1))
# scan file and extract sample counts for specified freqs
while (nXcvrs > 0):
# read the next datagram
_len, = st.unpack("<i", fid.read(4))
# if (fid.tell() == self.file_size):
if (fid.tell() == self.file_size):
break
# scan file for RAW datagrams
[dgType, dgTime] = self.readEKRaw_ReadDgHeader(fid, 0)
dgType = ''.join(dgType)
if dgType == 'RAW0':
# read channel ID
channel, = st.unpack("<h", fid.read(2))
channel = channel-1
# skip to the samplecount record and read
fid.seek(66, 1)
sampleCount, = st.unpack("<i", fid.read(4))
# adjust sampleCount based on the sample range to be extracted
if sampleRange[0] <= sampleCount:
if sampleRange[1] == 65535:
endIdx = sampleCount
else:
endIdx = sampleRange[1]
if endIdx > maxSampleRange:
endIdx = maxSampleRange
sampleCount = (endIdx - sampleRange[0]) + 1
else:
sampleCount = 0
# store sample number if required/valid
# idx = find(CIDs == channel);
print('CID : {}; channel : {}'.format(CIDs, channel))
idx = np.where(CIDs[:] == channel)[0][0]
nPings[idx] = nPings[idx] + 1
if np.size(idx)>0 and (sampleCount > 0)\
and (nPings[idx] >= ping) and (nSamples[idx] == 0):
nSamples[idx] = sampleCount
nXcvrs = nXcvrs - 1
# skip to the next datagram
fid.seek(_len - 72 - HEADER_LEN, 1)
else:
# Skip to the next datagram
fid.seek(_len - HEADER_LEN, 1)
_tmp, = st.unpack("<i", fid.read(4))
# reset file pointer
fid.seek(fPosition, 0)
return nSamples
def readEKRaw_AllocateSampledata(self, length):
length = int(length)
sampledata = Parameters({})
sampledata.channel = 0
sampledata.mode = 0
sampledata.transducerdepth = 0
sampledata.frequency = 0
sampledata.transmitpower = 0
sampledata.pulselength = 0
sampledata.bandwidth = 0
sampledata.sampleinterval = 0
sampledata.soundvelocity = 0
sampledata.absorptioncoefficient = 0
sampledata.heave = 0
sampledata.roll = 0
sampledata.pitch = 0
sampledata.temperature = 0
sampledata.trawlupperdepthvalid = 0
sampledata.trawlopeningvalid = 0
sampledata.trawlupperdepth = 0
sampledata.trawlopening = 0
sampledata.offset = 0
sampledata.count = 0
sampledata.power = np.zeros((1, length)) - 999
sampledata.alongship = np.zeros((1, length))
sampledata.athwartship = np.zeros((1, length))
return sampledata
def readEKRaw_Allocatedata(self, xcvrCount, xcvrConfig, aSize, rData,
sampleRange, nSamples, userNMEA):
if (rData.skinny > 0):
bSize = 1
else:
bSize = aSize
data = Parameters({})
data.config = [[]] * np.size(xcvrConfig)
data.pings = []
[data.pings.append(Parameters({})) for i in range(np.size(xcvrConfig))]
# allocate transceiver based data arrays
for n in range(xcvrCount):
data.config[n] = xcvrConfig[n]
data.pings[n].number = np.zeros((1, aSize))
data.pings[n].time = np.zeros((1, aSize))
data.pings[n].mode = np.zeros((1, bSize))
data.pings[n].transducerdepth = np.zeros((1, bSize))
data.pings[n].frequency = np.zeros((1, bSize))
data.pings[n].transmitpower = np.zeros((1, bSize))
data.pings[n].pulselength = np.zeros((1, bSize))
data.pings[n].bandwidth = np.zeros((1, bSize))
data.pings[n].sampleinterval = np.zeros((1, bSize))
data.pings[n].soundvelocity = np.zeros((1, bSize))
data.pings[n].absorptioncoefficient = np.zeros((1, bSize))
if rData.heave:
data.pings[n].heave = np.zeros((1, aSize))
if (rData.roll):
data.pings[n].roll = np.zeros((1, aSize))
if (rData.pitch):
data.pings[n].pitch = np.zeros((1, aSize))
if (rData.temp):
data.pings[n].temperature = np.zeros((1, aSize))
if (rData.trawl):
data.pings[n].trawlopeningvalid = np.zeros((1, aSize))
data.pings[n].trawlupperdepthvalid = np.zeros((1, aSize))
data.pings[n].trawlupperdepth = np.zeros((1, aSize))
data.pings[n].trawlopening = np.zeros((1, aSize))
data.pings[n].offset = np.zeros((1, bSize))
data.pings[n].count = np.zeros((1, aSize))
if (rData.power):
data.pings[n].power = np.zeros((np.int(nSamples[n]),
aSize)) - 999
if (rData.angle):
data.pings[n].alongship_e = np.zeros((np.int(nSamples[n]),
aSize))
data.pings[n].athwartship_e = np.zeros((np.int(nSamples[n]),
aSize))
data.pings[n].samplerange = [sampleRange[0], nSamples[n]]
if (rData.gps):
data.pings[n].seg = np.zeros((1, aSize))
# allocate GPS data arrays
if (rData.gps > 0):
data.gps = Parameters({'time': np.zeros((aSize, 1)),
'lat': np.zeros((aSize, 1)),
'lon': np.zeros((aSize, 1)),
'seg': np.zeros((aSize, 1)),
'len': aSize})
# allocate raw NMEA data array
if (rData.rawNMEA > 0):
data.NMEA.time = np.zeros((aSize, 1))
data.NMEA.string = aSize*[[]]
data.NMEA._len = aSize
# allocate vessel log data arrays
if (rData.vlog):
data.vlog = Parameters({'time': np.zeros((aSize, 1)),
'vlog': np.zeros((aSize, 1)),
'len': aSize})
if (rData.gps):
data.vlog.seg = np.zeros((aSize, 1))
# allocate vessel speed data arrays
if (rData.vSpeed):
data.vspeed = Parameters({'time': np.zeros((aSize, 1)),
'speed': np.zeros((aSize, 1)),
'len': aSize})
if (rData.gps):
data.vspeed.seg = np.zeros((aSize, 1))
# allocate annotation data arrays
if (rData.annotations):
data.annotations = Parameters({'time': np.zeros((aSize, 1)),
'len': aSize})
data.annotations.text = aSize*[[]]
return data
[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)
def readEKRaw_ParseNMEAstring(self, nmeaString, rData):
try:
# extract the NMEA datagram type
idx = list(self.find_all(nmeaString, ','))
_type = nmeaString[1:idx[0]].upper()
nmeadata = nmeaString[idx[0] + 1:]
# remove checksum - add trailing comma for output lacking last
# field
idx = list(self.find_all(nmeadata, '*'))
if len(idx) == 0:
nmeadata = nmeadata + ','
else:
nmeadata = nmeadata[0: idx] + ','
# fill any empty fields (strread can't handle them)
nmeadata = nmeadata.replace(',,', ',0')
# process datagram by type - only process datagrams we're
# interested in
nmea = Parameters({'type': None,
'string': ''})
if _type[2:] == 'GGA':
if (rData.gps and _type == rData.gpsSource):
# parse the rest of the nmea text
form = '%2d %2d %f %2d %f %c %3d %f %c %d %d %f %f %c %f %c %f %d'
lll = nmeadata.split(',')
lll_ = [o.strip('\\') for o in lll]
# define GPS GPGGA datagram
nmea = Parameters({'type': _type,
'time': [np.double(lll_[0]),
np.double(lll_[1]),
float(lll_[2])],
'lat': np.double(lll_[3]) + float(lll_[4]) / 60,
'lat_hem': lll_[5],
'lon': np.double(lll_[6]) + float(lll_[7]) / 60,
'lon_hem': lll_[8],
'fix': np.double(lll_[9]),
'nsat': np.double(lll_[10]),
'precision': float(lll_[11]),
'msl_alt': float(lll_[12]),
'msl_unit': lll_[13],
'geoidal_alt': float(lll_[14]),
'geoidal_unit': lll_[15],
'dif_age': float(lll_[16]),
'dif_sta': np.double(lll_[17])})
elif _type[2:] == 'GLL':
if (rData.gps and _type == rData.gpsSource):
# parse the rest of the nmea text
lll = nmeadata.split(',')
lll_ = [o.strip('\\') for o in lll]
# convert status to GGA fix
fix = lll_[9] == 'A'
# define GPS geographic position datagram
nmea = Parameters({'type': _type,
'lat': np.double(lll_[0]) + float(lll_[1]) / 60,
'lat_hem': lll_[2],
'lon': np.double(lll_[3]) + float(lll_[4]) / 60,
'lon_hem': lll_[5],
'time': [np.double(lll_[6]),
np.double(lll_[7]),
float(lll_[8])],
'fix': fix})
elif _type[2:] == 'VTG':
if (rData.vSpeed):
# parse the rest of the nmea text
lll = nmeadata.split(',')
lll_ = [o.strip('\\') for o in lll]
# define Course Over Ground datagram
nmea = Parameters({'type': _type,
'true_cov': float(lll_[0]),
'tcov_label': lll_[1],
'mag_cov': float(lll_[2]),
'mcov_label': lll_[3],
'sog_knts': float(lll_[4]),
'sogn_unit': lll_[5],
'sog_kph': float(lll_[6]),
'sogk_unit': lll_[7]})
elif _type[2:] == 'VLW':
if (rData.vlog):
# parse the rest of the nmea text
lll = nmeadata.split(',')
lll_ = [o.strip('\\') for o in lll]
# define Distance Traveled through Water datagram
nmea = Parameters({'type': _type,
'total_cum_dist': float(lll_[0]),
'tcd_unit': lll_[1],
'dist_since_reset': float(lll_[2]),
'dsr_unit': lll_[3]})
else:
# unknown datagram type
nmea = Parameters({'type': _type, 'string': nmeadata})
except TypeError():
# malformed NMEA string
nmea = Parameters({'type': 'XXXXX', 'string': nmeaString})
return nmea
def readEKRaw_ReadSampledata(self, fid, sampledata):
sampledata.channel, = st.unpack("<h", fid.read(2))
mode_low, = st.unpack("<b", fid.read(1))
mode_high, = st.unpack("<b", fid.read(1))
sampledata.mode = 256 * mode_high + mode_low
sampledata.transducerdepth, = st.unpack("<f", fid.read(4))
sampledata.frequency, = st.unpack("<f", fid.read(4))
sampledata.transmitpower, = st.unpack("<f", fid.read(4))
sampledata.pulselength, = st.unpack("<f", fid.read(4))
sampledata.bandwidth, = st.unpack("<f", fid.read(4))
sampledata.sampleinterval, = st.unpack("<f", fid.read(4))
sampledata.soundvelocity, = st.unpack("<f", fid.read(4))
sampledata.absorptioncoefficient, = st.unpack("<f", fid.read(4))
sampledata.heave, = st.unpack("<f", fid.read(4))
sampledata.roll, = st.unpack("<f", fid.read(4))
sampledata.pitch, = st.unpack("<f", fid.read(4))
sampledata.temperature, = st.unpack("<f", fid.read(4))
sampledata.trawlupperdepthvalid, = st.unpack("<h", fid.read(2))
sampledata.trawlopeningvalid, = st.unpack("<h", fid.read(2))
sampledata.trawlupperdepth, = st.unpack("<f", fid.read(4))
sampledata.trawlopening, = st.unpack("<f", fid.read(4))
sampledata.offset, = st.unpack("<i", fid.read(4))
sampledata.count, = st.unpack("<i", fid.read(4))
if (sampledata.count > 0):
if (sampledata.mode != 2):
sampledata.power = np.zeros(sampledata.count, dtype='<h')
fid.readinto(sampledata.power)
sampledata.power = sampledata.power * 0.011758984205624
if (sampledata.mode > 1):
angle = np.zeros((sampledata.count, 2), dtype='<b')
fid.readinto(angle)
angle = angle.T
sampledata.athwartship = angle[0, :].T
sampledata.alongship = angle[1, :].T
return sampledata
[docs] def read_EK60(self, fname):
"""Read one .raw EK60 file"""
# typeNames = {
# 'int8' :'b',
# 'uint8' :'B',
# 'int16' :'h',
# 'uint16' :'H',
# 'int32' :'i',
# 'uint32' :'I',
# 'int64' :'q',
# 'uint64' :'Q',
# 'float' :'f',
# 'double' :'d',
# 'char' :'s'}
# define constants
# Bytes in datagram header
HEADER_LEN = 12
# size in array elements of allocation "chunks"
CHUNK_SIZE = 2500
# max distance (m) between rstat GPS fix and first read GPS fix
MAXRSTATGPSM = 500
rData = Parameters({})
# define default parameters
# Return all frequencies by default
returnFreqs = -1
# Assume logging PC time and timezone set correctly
timeOffset = 0
# Return alon and athw electrical angles by default
rData.angle = True
# Do not return annotations by default
rData.annotations = False
# default is non-continuation mode
rData.cont = False
# Return GPS data (even if none available)
rData.gps = True
# Return all data, not just GPS
rData.gpsOnly = False
# Return GPS as double
rData.gpsType = 'double'
# Use GGA datagrams as GPS data source
rData.gpsSource = 'GPGGA'
# Do not return heave by default
rData.heave = True
# Do not return pitch by default
rData.pitch = True
# Return power by default
rData.power = True
# Do not return raw NMEA data
rData.rawNMEA = False
# Do not return roll by default
rData.roll = True
# Return all ping-by-ping param data by default
rData.skinny = 0
# Do not return trawl data by default
rData.trawl = False
# Do not return temperature by default
rData.temp = False
# Do not return user defined NMEA data
rData.uNMEA = False
# Do not return vessel log data
rData.vlog = True
# Return vessel log as single
rData.vlogType = 'single'
# Use the SD talker ID for VLW datagrams
rData.vlwSource = 'SD'
# Do not return vessel speed data
rData.vSpeed = True
# Return vessel speed as single
rData.vSpeedType = 'single'
# Use the GP talker ID for VTG datagrams
rData.vtgSource = 'GP'
# Do not limit by ping and read all pings by default
pingRange = [1, 1, np.Inf]
# Do not limit by time by default
timeRange = [0, np.Inf]
# Do not limit by sample by default
sampleRange = [1, 2**16-1]
# Effectively no maximum limit
maxSampleRange = 30000
# Geographic region is undefined by default
geoRegion = -1
# No explicit polygon connectivity provided
geoConn = -1
# Don't limit data by geographic region (ROI)
rData.doGeo = False
# All points are in the ROI by default
isInRegion = True
# By default the reader parameters are unset.
rparms = -1
# By default, do not test for spatial gaps in data
maxSegGap = -1
# No progress callback function defined
progressCallback = []
# Minimum percent change in progress before callback is called
progressGranularity = 2
# By default the user defined NMEA array is undefined
userNMEA = {}
# Keep reading 5000 bytes past a bad datagram for a good one
maxBadBytes = 5000
# By default assume beam mode will not change.
allowBeamModeChange = False
# Initial beam mode is unset.
beamMode = -1
# open the raw file
fid = open(fname, 'rb')
fid.seek(0, 2)
self.file_size = fid.tell()
fid.seek(0)
# read configuration datagram 'CON0'
[config, frequencies] = self.readEKRaw_ReadHeader(fid)
# define channel ID array
CIDs = np.arange(config.header.transceivercount)
# no specific freqs set - return all frequencies in the file
returnFreqs = frequencies
header = config.header
# get initial sample counts for each freq - first try the 4th ping.
# certain files of mine had unusual sample numbers for the first
# ping so looking a bit into the file avoids this.
nSamples = self.readEKRaw_GetSampleCount(fid,
CIDs,
sampleRange,
maxSampleRange, 4)
# check for errant sampleRange parameters or very short files
badRanges = np.where(np.squeeze(nSamples) == 0)
if np.size(badRanges) > 0:
# try to get sample range from the first ping (very short file)
nSFP = self.readEKRaw_GetSampleCount(fid,
CIDs[badRanges],
sampleRange,
maxSampleRange,
1)
nSamples[badRanges] = nSFP
badRanges = np.where(np.squeeze(nSamples) == 0)[0]
if np.size(badRanges) > 0:
# still unable to get sample count, either the file contains
# no data or the sample count is out of range.
# Issue warning and remove channels
[config.transceiver.remove(config.transceiver[i])
for i in -np.sort(-badRanges)]
[CIDs.remove(CIDs[i]) for i in -np.sort(-badRanges)]
# update header data to reflect change
config.header.transceivercount =\
config.header.transceivercount - len(badRanges)
# return if we're left with no xcvrs
if (config.header.transceivercount <= 0):
print('The specified frequencies are not in the file.')
fid.close()
header = -1
data = -1
# allocate sampledata structure
sampledata = self.readEKRaw_AllocateSampledata(np.max(nSamples))
transceivercount = config.header.transceivercount
# allocate return structure
data = self.readEKRaw_Allocatedata(transceivercount,
config.transceiver,
CHUNK_SIZE,
rData,
sampleRange,
nSamples,
userNMEA)
# allocate array bookkeeping variables
nGPS = 1
nNMEA = 1
nVlog = 1
nVSpeed = 1
nTAG = 1
arraySize = np.zeros((transceivercount, 2), dtype=np.int)
arraySize[:, 0] = np.int(CHUNK_SIZE)
arraySize[:, 1] = nSamples.T
nPings = np.squeeze(np.tile(0, (1, transceivercount)))
nuNMEA = 0
# define reader state variables
# define state variables using defaults
gpsSeg = Parameters({'_in': 0, 'out': 0, 'inRegion': True})
pingCounter = np.squeeze(np.tile(0, (1, transceivercount)))
# progress callback - determine divisor
loopy = True
# read file, processing individual datagrams
oo = 0
while loopy is True:
# check if we're at the end of the file
oo += 1
b = ("Loading" + "..." + str(oo))
sys.stdout.write('\r'+b)
if (fid.tell() == self.file_size):
loopy = False
break
_len, = st.unpack("<i", fid.read(4))
# read datagram header
[dgType, dgTime] = self.readEKRaw_ReadDgHeader(fid, timeOffset)
# if reading subsets - check if we're done
if (dgTime > timeRange[1]):
# move file pointer back to beginning of this datagram
fid.seek(-(HEADER_LEN + 4), 1)
break
dgType = ''.join(dgType)
if dgType == 'NME0':
# Process NMEA datagram
# read datagram
text = [str(st.unpack("<s", fid.read(1)))[3]
for ui in range(1, _len - HEADER_LEN+1)]
text = ''.join(text)
# store the raw NMEA text
if (rData.rawNMEA):
# check len of arrays - extend if required
if (nNMEA > data.NMEA._len):
data.NMEA.time = np.concatenate((data.NMEA.time,
np.zeros((CHUNK_SIZE,
1))))
data.NMEA.string = data.NMEA.string + CHUNK_SIZE*[[]]
data.NMEA._len = data.NMEA._len + CHUNK_SIZE
# assign raw NMEA data
data.NMEA.time[nNMEA-1] = dgTime
data.NMEA.string[nNMEA-1] = text
nNMEA = nNMEA + 1
if rData.gps or rData.vlog or rData.vSpeed or rData.uNMEA:
# parse datagram
nmea = self.readEKRaw_ParseNMEAstring(text, rData)
# extract datagrams of interest, process and store
# GPS Fix data
if rData.gps and nmea.type == rData.gpsSource:
# check that GPS data has valid fix - drop if not
if (nmea.fix == 0):
st.unpack("<i", fid.read(4))
pass
# check len of arrays - extend if required
if (nGPS > data.gps.len):
data.gps.time =\
np.concatenate((data.gps.time,
np.zeros((CHUNK_SIZE, 1))))
data.gps.lat =\
np.concatenate((data.gps.lat,
np.zeros((CHUNK_SIZE, 1))))
data.gps.lon =\
np.concatenate((data.gps.lon,
np.zeros((CHUNK_SIZE, 1))))
data.gps.seg =\
np.concatenate((data.gps.seg,
np.zeros((CHUNK_SIZE, 1))))
data.gps.len = data.gps.len + CHUNK_SIZE
# set time
data.gps.time[nGPS-1] = dgTime
# set lat/lon signs and store values
if (nmea.lat_hem == 'S'):
data.gps.lat[nGPS-1] = -nmea.lat
else:
data.gps.lat[nGPS-1] = nmea.lat
if (nmea.lon_hem == 'W'):
data.gps.lon[nGPS-1] = -nmea.lon
else:
data.gps.lon[nGPS-1] = nmea.lon
# are we limiting data by geographic region?
if (rData.doGeo):
# yes - test if fix is in ROI
print('Cannot process geographic data yet...' +
'Switching to basic reading')
# if (geoConn[0] < 0):
# # no connectivity provided
# [isInRegion, bnd] = inpoly([data.gps.lat(nGPS),
# data.gps.lon(nGPS)
# ], geoRegion)
# else:
# # explicit connectivity provided
# [isInRegion, bnd] = inpoly([data.gps.lat(nGPS), data.gps.lon(nGPS)], geoRegion, geoConn)
else:
# no - always "in"
isInRegion = True
# set prior GPS fix
if (nGPS == 1):
print('no geographic processing yet')
# # no priors in data.gps
# if (isstruct(rparms)):
# # check if the rparms fix is close
# # this is a sanity check on the rstat GPS fix
# dist = vdist(data.gps.lat(nGPS), data.gps.lon(nGPS), rparms.lat, rparms.lon)
# if dist < MAXRSTATGPSM:
# # pick up last fix from rparms struct
# lastLat = rparms.lat
# lastLon = rparms.lon
# else:
# # too far apart - replicate first fix
# lastLat = data.gps.lat[nGPS-1]
# lastLon = data.gps.lon[nGPS-1]
# else:
# # no rparms - just replicate first fix
# lastLat = data.gps.lat[nGPS-1]
# lastLon = data.gps.lon[nGPS-1]
# # set initial segment value
# if (isInRegion):
# # only need to set in state...
# gpsSeg._in = 1
else:
# priors read - assign
lastLat = data.gps.lat[nGPS-2]
lastLon = data.gps.lon[nGPS-2]
# check if vessel moved and calculate delta
if (maxSegGap > 0) and\
((data.gps.lat[nGPS-1] - lastLat != 0) or
(data.gps.lon[nGPS-1] - lastLon != 0)):
# moved - calculate distance between fixes (in m)
dist = vdist(data.gps.lat(nGPS),
data.gps.lon(nGPS),
lastLat,
lastLon)
else:
# no movement
dist = 0
# determine segment value
if (maxSegGap > 0) and (dist > maxSegGap):
# distance exceeds threshold - increment segment
if (gpsSeg.inRegion):
# in ROI - increment "in" counter
gpsSeg._in = gpsSeg._in + int(1)
data.gps.seg[nGPS-1] = gpsSeg._in
else:
# out of ROI - increment "out" counter
gpsSeg.out = gpsSeg.out - int(1)
data.gps.seg[nGPS-1] = gpsSeg.out
elif (rData.doGeo):
if (isInRegion):
# in ROI
if gpsSeg.inRegion is False:
# transitioned - increment in counter
gpsSeg.inRegion = True
gpsSeg._in = gpsSeg._in + int(1)
data.gps.seg[nGPS-1] = gpsSeg._in
else:
# out of ROI
if (gpsSeg.inRegion):
# transitioned - decrement out counter
gpsSeg.inRegion = False
gpsSeg.out = gpsSeg.out - int(1)
data.gps.seg[nGPS-1] = gpsSeg.out
else:
# not using ROI's or segGap-always in one segment
data.gps.seg[nGPS-1] = gpsSeg._in
# increment GPS counter
nGPS += 1
# store vessel log data
elif rData.vlog and\
rData.vlwSource + 'VLW' == nmea.type and\
isInRegion:
# check len of arrays - extend if required
if (nVlog > data.vlog.len):
data.vlog.time =\
np.concatenate((data.vlog.time,
np.zeros((CHUNK_SIZE, 1))))
data.vlog.vlog =\
np.concatenate((data.vlog.vlog,
np.zeros((CHUNK_SIZE, 1))))
if (rData.gps):
data.vlog.seg =\
np.concatenate((data.vlog.seg,
np.zeros((CHUNK_SIZE, 1))))
data.vlog.len = data.vlog.len + CHUNK_SIZE
data.vlog.time[nVlog-1] = dgTime
data.vlog.vlog[nVlog-1] = nmea.total_cum_dist
if (rData.gps):
data.vlog.seg[nVlog-1] = gpsSeg._in
nVlog += 1
# store vessel speed data
elif rData.vSpeed and\
rData.vtgSource + 'VTG' == nmea.type and\
isInRegion:
# check len of arrays - extend if required
if (nVSpeed > data.vspeed.len):
data.vspeed.time =\
np.concatenate((data.vspeed.time,
np.zeros((CHUNK_SIZE, 1))))
data.vspeed.speed =\
np.concatenate((data.vspeed.speed,
np.zeros((CHUNK_SIZE, 1))))
if (rData.gps):
data.vspeed.seg =\
np.concatenate((data.vspeed.seg,
np.zeros((CHUNK_SIZE, 1))))
data.vspeed.len = data.vspeed.len + CHUNK_SIZE
data.vspeed.time[nVSpeed-1] = dgTime
data.vspeed.speed[nVSpeed-1] = nmea.sog_knts
data.vspeed.seg[nVSpeed-1] = gpsSeg._in
nVSpeed += 1
# store user defined NMEA strings
lastType = dgType
elif dgType == 'TAG0':
# Process annotation datagram
text = [str(st.unpack("<s", fid.read(1)))[3]
for ui in range(1, _len - HEADER_LEN+1)]
text = ''.join(text)
if (rData.annotations):
# check len of arrays - extend if required
if (nTAG > data.annotations.len):
data.annotations.time =\
np.concatenate((data.annotations.time,
np.zeros((CHUNK_SIZE, 1))))
data.annotations.text =\
data.annotations.text + CHUNK_SIZE * [[]]
data.annotations.len =\
data.annotations.len + CHUNK_SIZE
data.annotations.time[nTAG-1] = dgTime
data.annotations.text[nTAG-1] = text
nTAG += 1
lastType = dgType
elif dgType == 'RAW0':
# Process sample datagrams
if (CIDs[0] >= 0):
sampledata = self.readEKRaw_ReadSampledata(fid, sampledata)
idx2 = list(np.where(CIDs == sampledata.channel-1)[0])
idx = np.where(CIDs == sampledata.channel-1)[0][0]
else:
idx2 = []
idx = np.nf
text = [str(st.unpack("<s", fid.read(1)))[3]
for ui in range(1, _len - HEADER_LEN+1)]
text = ''.join(text)
# check if we're storing this frequency...
if (len(idx2) > 0):
# check if the beam mode has changed
if ~allowBeamModeChange:
if beamMode < 0:
# set the initial beam mode
beamMode = sampledata.mode
else:
# check if it has changed
if beamMode != sampledata.mode:
raise InstrAmbiguity('readEKRaw:ModeChange, ' +
'Your beam mode has ' +
'changed. Set the ' +
'AllowModeChange ' +
'parameter to true to ' +
'read this file.')
break
# increment total ping counter
pingCounter[idx] = pingCounter[idx] + 1
# if reading subsets - check if we're done
if (pingCounter[idx] > pingRange[2]):
# move file pointer back to beginning of this datagram
fid.seek(-(_len + 4), 1)
# tick back ping counter
pingCounter[idx] = pingCounter[idx] - 1
break
# check if we're storing this ping...
if (pingCounter[idx] >= pingRange[0])\
and (dgTime >= timeRange[0])\
and (isInRegion)\
and (np.mod(pingCounter[idx], pingRange[1]) == 0):
# check length of arrays - extend if required
if (nPings[idx] > arraySize[idx, 0]):
arraySize[idx, 0] = arraySize[idx, 0] + CHUNK_SIZE
data.pings[idx].number[0, nPings[idx]:arraySize[idx, 0]+1] = 0
data.pings[idx].time[0, nPings[idx]:arraySize[idx, 0]+1] = 0
if (rData.skinny)==0:
data.pings[idx].mode[0, nPings[idx]:arraySize[idx, 0]+1] = 0
data.pings[idx].transducerdepth[0, nPings[idx]:arraySize[idx, 0]+1] = 0
data.pings[idx].frequency[0, nPings[idx]:arraySize[idx, 0]+1] = 0
data.pings[idx].transmitpower[0, nPings[idx]:arraySize[idx, 0]+1] = 0
data.pings[idx].pulselength[0, nPings[idx]:arraySize[idx, 0]+1] = 0
data.pings[idx].bandwidth[0, nPings[idx]:arraySize[idx, 0]+1] = 0
data.pings[idx].sampleinterval[0, nPings[idx]:arraySize[idx, 0]+1] = 0
data.pings[idx].soundvelocity[0, nPings[idx]:arraySize[idx, 0]+1] = 0
data.pings[idx].absorptioncoefficient[0, nPings[idx]:arraySize[idx, 0]+1] = 0
data.pings[idx].offset[0, nPings[idx]:arraySize[idx, 0]+1] = 0
if (rData.heave):
data.pings[idx].heave[0, nPings[idx]:arraySize[idx, 0]+1] = 0
if (rData.roll):
data.pings[idx].roll[0, nPings[idx]:arraySize[idx, 0]+1] = 0
if (rData.pitch):
data.pings[idx].pitch[0, nPings[idx]:arraySize[idx, 0]+1] = 0
if (rData.temp):
data.pings[idx].temperature[0, nPings[idx]:arraySize[idx, 0]+1] = 0
if (rData.trawl):
data.pings[idx].trawlopeningvalid[0, nPings[idx]:arraySize[idx, 0]+1] = 0
data.pings[idx].trawlupperdepthvalid[0, nPings[idx]:arraySize[idx, 0]+1] = 0
data.pings[idx].trawlupperdepth[0, nPings[idx]:arraySize[idx, 0]+1] = 0
data.pings[idx].trawlopening[0, nPings[idx]:arraySize[idx, 0]+1] = 0
data.pings[idx].count[0, nPings[idx]:arraySize[idx, 0]+1] = 0
if (rData.power) and (beamMode != 2):
data.pings[idx].power[:, nPings[idx]:arraySize[idx, 0]+1] = -999
if (rData.angle) and ((beamMode > 1) or (beamMode < 0)):
data.pings[idx].alongship_e[:, nPings[idx]:arraySize[idx, 0]+1] = 0
data.pings[idx].athwartship_e[:, nPings[idx]:arraySize[idx, 0]+1] = 0
if (rData.gps):
data.pings[idx].seg[:, nPings[idx]:arraySize[idx, 0]+1] = 0
# copy this pings data into output arrays
data.pings[idx].number[0, nPings[idx]] =\
pingCounter[idx]
data.pings[idx].time[0, nPings[idx]] = dgTime
if (rData.skinny > 0):
if (rData.skinny == nPings[idx]):
# operating in "skinny" mode and this is
# the ping we want to read our parameters from
data.pings[idx].mode[0, 0] = sampledata.mode
data.pings[idx].transducerdepth[0, 0] =\
sampledata.transducerdepth
data.pings[idx].frequency[0, 0] =\
sampledata.frequency
data.pings[idx].transmitpower[0, 0] =\
sampledata.transmitpower
data.pings[idx].pulselength[0, 0] =\
sampledata.pulselength
data.pings[idx].bandwidth[0, 0] =\
sampledata.bandwidth
data.pings[idx].sampleinterval[0, 0] =\
sampledata.sampleinterval
data.pings[idx].soundvelocity[0, 0] =\
sampledata.soundvelocity
data.pings[idx].absorptioncoefficient[0, 0] =\
sampledata.absorptioncoefficient
data.pings[idx].offset[0, 0] =\
sampledata.offset
else:
# operating in "fat" mode. Store all parameters.
data.pings[idx].mode[0, nPings[idx]] =\
sampledata.mode
data.pings[idx].transducerdepth[0, nPings[idx]] =\
sampledata.transducerdepth
data.pings[idx].frequency[0, nPings[idx]] =\
sampledata.frequency
data.pings[idx].transmitpower[0, nPings[idx]] =\
sampledata.transmitpower
data.pings[idx].pulselength[0, nPings[idx]] =\
sampledata.pulselength
data.pings[idx].bandwidth[0, nPings[idx]] =\
sampledata.bandwidth
data.pings[idx].sampleinterval[0, nPings[idx]] =\
sampledata.sampleinterval
data.pings[idx].soundvelocity[0, nPings[idx]] =\
sampledata.soundvelocity
data.pings[idx].absorptioncoefficient[0, nPings[idx]] =\
sampledata.absorptioncoefficient
data.pings[idx].offset[0, nPings[idx]] =\
sampledata.offset
if (rData.heave):
data.pings[idx].heave[0, nPings[idx]] =\
sampledata.heave
if (rData.roll):
data.pings[idx].roll[0, nPings[idx]] =\
sampledata.roll
if (rData.pitch):
data.pings[idx].pitch[0, nPings[idx]] =\
sampledata.pitch
if (rData.temp):
data.pings[idx].temperature[0, nPings[idx]] =\
sampledata.temperature
if (rData.trawl):
data.pings[idx].trawlopeningvalid[0, nPings[idx]] =\
sampledata.trawlopeningvalid
data.pings[idx].trawlupperdepthvalid[0, nPings[idx]] =\
sampledata.trawlupperdepthvalid
data.pings[idx].trawlupperdepth[0, nPings[idx]] =\
sampledata.trawlupperdepth
data.pings[idx].trawlopening[0, nPings[idx]] =\
sampledata.trawlopening
# store sample data
if (sampledata.count > 0):
# handle "subsettable" data
if (sampleRange[0] <= sampledata.count):
# determine ending sample index
if (sampleRange[1] > sampledata.count):
endIdx = sampledata.count
# cap number of data samples stored
# to maxSampleRange
if endIdx > maxSampleRange:
endIdx = maxSampleRange
else:
endIdx = sampleRange[1]
# store actual sample count value
count = (endIdx - sampleRange[0])
data.pings[idx].count[0, nPings[idx]] = count
# check depth of arrays - extend if required
if (arraySize[idx, 1] < count):
nSampAdd = count - arraySize[idx, 1]
data.pings[idx].samplerange[1] = count
if (rData.power):
for i in idx:
data.pings[i].power =\
np.concatenate((data.pings[i].power,
-999+np.zeros((nSampAdd,
np.size(data.pings[i].power,
1)))))
if (rData.angle):
for i in idx:
data.pings[i].alongship_e =\
np.concatenate((data.pings[i].alongship_e,
np.zeros((nSampAdd,
np.size(data.pings[i].alongship_e,
1)))))
data.pings[i].athwartship_e =\
np.concatenate((data.pings[i].athwartship_e,
np.zeros((nSampAdd,
np.size(data.pings[i].athwartship_e,
1)))))
arraySize[idx, 1] = count
if (rData.power) and (sampledata.mode != 2):
# store power
data.pings[idx].power[0:count, nPings[idx]] =\
sampledata.power[sampleRange[0]:endIdx]
if (rData.angle) and (sampledata.mode > 1):
# store angles
data.pings[idx].alongship_e[0:count, nPings[idx]] =\
sampledata.alongship[sampleRange[0]:endIdx]
data.pings[idx].athwartship_e[0:count, nPings[idx]] =\
sampledata.athwartship[sampleRange[0]:endIdx]
else:
data.pings[idx].count[0, nPings[idx]] = 0
# store GPS line segment value
if (rData.gps):
data.pings[idx].seg[0, nPings[idx]] = gpsSeg._in
# increment stored ping counter
nPings[idx] = nPings[idx] + 1
lastIDX = idx
lastType = dgType
elif dgType == 'CON1':
# Read ME70 CON datagram
conTextRaw = [str(st.unpack("<s", fid.read(1)))[3]
for ui in range(1, _len - HEADER_LEN+1)]
conTextRaw = ''.join(text)
# return configuration data as a char array
data.conf.time = dgTime
data.conf.text = conTextRaw
lastType = dgType
elif dgType == 'SVP0':
# Read Sound Velocity Profile datagram
svpText = [str(st.unpack("<s", fid.read(1)))[3]
for ui in range(1, _len - HEADER_LEN+1)]
svpText = ''.join(text)
lastType = dgType
# Process unknown datagrams
else:
if (maxBadBytes > 0):
# Display warning - try to find next datagram
print('readEKRaw:Datagram, Invalid datagram at offset ' +
'{} (d) - Searching for next datagram...' +
''.format(fid.tell()))
# rewind to last good dgram header
fid.seek(-(lastLen + 8), 1)
ok = False
nBytes = 0
while ok is False:
try:
# read 4 bytes and look for known datagram
dgType = [str(st.unpack("<s", fid.read(1)))[3]
for ui in range(1, 4+1)]
dgType = ''.join(dgType)
except:
print('Unable to locate next valid datagram.')
break
nBytes = nBytes + 4
if (dgType == 'CON1')\
or (dgType == 'RAW0')\
or (dgType == 'TAG0')\
or (dgType == 'NME0')\
or (dgType == 'SVP0'):
if (nBytes < (lastLen - 8))\
and (lastType == 'RAW0'):
# only need to drop bad RAW0 datagrams
nPings[lastIDX-1] = nPings[lastIDX-1] - 1
# rewind to beginning of datagram
fid.seek(-12, 1)
ok = True
else:
# rewind 3 bytes
fid.seek(-3, 1)
nBytes = nBytes - 3
# check if we're giving up
if (nBytes > maxBadBytes):
print('Unable to locate valid datagram. ' +
'Giving up.')
break
if (ok is False):
break
else:
print('Invalid datagram at offset {} (d) - ' +
'Aborting...'.format(fid.tell()))
break
# datagram length is repeated...
lastLen, = st.unpack("<i", fid.read(4))
# close file.
fid.close()
# trim GPS/vlog/vspeed/annotation arrays
if (rData.gps):
data.gps.time = data.gps.time[0:nGPS]
data.gps.lat = data.gps.lat[0:nGPS]
data.gps.lon = data.gps.lon[0:nGPS]
data.gps.seg = data.gps.seg[0:nGPS]
if (rData.vlog):
data.vlog.time = data.vlog.time[0:nVlog]
data.vlog.vlog = data.vlog.vlog[0:nVlog]
if (rData.gps):
data.vlog.seg = data.vlog.seg[0:nVlog]
if (rData.vSpeed):
data.vspeed.time = data.vspeed.time[0:nVSpeed]
data.vspeed.speed = data.vspeed.speed[0:nVSpeed]
if (rData.gps):
data.vspeed.seg = data.vspeed.seg[0:nVSpeed]
if (rData.annotations):
data.annotations.time = data.annotations.time[0:nTAG]
data.annotations.text = data.annotations.text[0:nTAG]
if (rData.rawNMEA):
data.NMEA.time = data.NMEA.time[0:nTAG]
data.NMEA.string = data.NMEA.string[0:nTAG]
# for each transducer...
for idx in range(transceivercount):
if (nPings[idx] <= arraySize[idx, 0]):
# Trim output data arrays
data.pings[idx].number = data.pings[idx].number[0, 0:nPings[idx]]
data.pings[idx].time = data.pings[idx].time[0, 0:nPings[idx]]
if (rData.skinny == 0):
data.pings[idx].mode =\
data.pings[idx].mode[0, 0:nPings[idx]]
data.pings[idx].transducerdepth =\
data.pings[idx].transducerdepth[0, 0:nPings[idx]]
data.pings[idx].frequency =\
data.pings[idx].frequency[0, 0:nPings[idx]]
data.pings[idx].transmitpower =\
data.pings[idx].transmitpower[0, 0:nPings[idx]]
data.pings[idx].pulselength =\
data.pings[idx].pulselength[0, 0:nPings[idx]]
data.pings[idx].bandwidth =\
data.pings[idx].bandwidth[0, 0:nPings[idx]]
data.pings[idx].sampleinterval =\
data.pings[idx].sampleinterval[0, 0:nPings[idx]]
data.pings[idx].soundvelocity =\
data.pings[idx].soundvelocity[0, 0:nPings[idx]]
data.pings[idx].absorptioncoefficient =\
data.pings[idx].absorptioncoefficient[0, 0:nPings[idx]]
data.pings[idx].offset =\
data.pings[idx].offset[0, 0:nPings[idx]]
if (rData.heave):
data.pings[idx].heave =\
data.pings[idx].heave[0, 0:nPings[idx]]
if (rData.roll):
data.pings[idx].roll =\
data.pings[idx].roll[0, 0:nPings[idx]]
if (rData.pitch):
data.pings[idx].pitch =\
data.pings[idx].pitch[0, 0:nPings[idx]]
if (rData.temp):
data.pings[idx].temperature =\
data.pings[idx].temperature[0, 0:nPings[idx]]
if (rData.gps):
data.pings[idx].seg =\
data.pings[idx].seg[0, 0:nPings[idx]]
data.pings[idx].count =\
data.pings[idx].count[0, 0:nPings[idx]]
data.pings[idx].samplerange[1] =\
np.int(data.pings[idx].samplerange[1][0])
if (rData.power) and (beamMode != 2):
data.pings[idx].power =\
data.pings[idx].power[:, 0:nPings[idx]]
if (rData.angle) and ((beamMode > 1) or (beamMode < 0)):
data.pings[idx].alongship_e =\
data.pings[idx].alongship_e[:, 0:nPings[idx]]
data.pings[idx].athwartship_e =\
data.pings[idx].athwartship_e[:, 0:nPings[idx]]
# remove empty power or angle fields
if (beamMode == 2):
# this beam mode contains no power data - remove the power field
data.pings._delattr__('power')
if (beamMode == 1):
# this beam mode contains no power data - remove the power field
data.pings._delattr__('alongship_e')
data.pings._delattr__('athwartship_e')
return header, data
def readEKRaw_GetCalParms(self, header, data):
pingNo = 0
# determine the number of transceivers in the xcvr struct
nXcvrs = len(data.config)
# create the calibration parameters structure
calParms = []
d = {'soundername': '',
'frequency': 0,
'soundvelocity': 0,
'sampleinterval': 0,
'absorptioncoefficient': 0,
'gain': 0,
'equivalentbeamangle': 0,
'pulselengthtable': 0,
'gaintable': 0,
'sacorrectiontable': 0,
'transmitpower': 0,
'pulselength': 0,
'anglesensitivityalongship': 0,
'anglesensitivityathwartship': 0,
'angleoffsetalongship': 0,
'angleoffsetathwartship': 0,
'transducerdepth': 0}
[calParms.append(Parameters(d)) for i in range(nXcvrs)]
# loop thru available transceivers and extract calibration parameters
for n in range(nXcvrs):
calParms[n].soundername = header.soundername
calParms[n].frequency = data.config[n].frequency
calParms[n].soundvelocity = data.pings[n].soundvelocity[pingNo]
calParms[n].sampleinterval = data.pings[n].sampleinterval[pingNo]
calParms[n].absorptioncoefficient =\
data.pings[n].absorptioncoefficient[pingNo]
calParms[n].gain = data.config[n].gain
calParms[n].equivalentbeamangle =\
data.config[n].equivalentbeamangle
calParms[n].pulselengthtable = data.config[n].pulselengthtable
calParms[n].gaintable = data.config[n].gaintable
calParms[n].sacorrectiontable = data.config[n].sacorrectiontable
calParms[n].transmitpower = data.pings[n].transmitpower[pingNo]
calParms[n].pulselength = data.pings[n].pulselength[pingNo]
calParms[n].anglesensitivityalongship =\
data.config[n].anglesensitivityalongship
calParms[n].anglesensitivityathwartship =\
data.config[n].anglesensitivityathwartship
calParms[n].angleoffsetalongship =\
data.config[n].angleoffsetalongship
calParms[n].angleoffsetathwartship =\
data.config[n].angleoffsetathwartship
calParms[n].transducerdepth = data.pings[n].transducerdepth[pingNo]
return calParms