1105 lines
43 KiB
Python
1105 lines
43 KiB
Python
# -*- coding: utf-8 -*-
|
|
"""
|
|
Created on Wen Jun 08 2016
|
|
|
|
@author: Adrien Vergne
|
|
"""
|
|
|
|
###############################################################################
|
|
# #
|
|
# CLASSES and METHODS for loading AQUASCAT Data #
|
|
# #
|
|
###############################################################################
|
|
|
|
# <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|
# Adrien VERGNE - Irstea Lyon - 2016
|
|
# Program Python 3.5
|
|
# <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
|
|
|
|
|
###############################################################################
|
|
# LOADING LIBRARIES #
|
|
###############################################################################
|
|
|
|
import numpy as np # Matrix scientific computing
|
|
import os # Path / folder access management
|
|
import datetime # for time management
|
|
# Loading binary data
|
|
import struct as st
|
|
import glob
|
|
import pickle
|
|
import time
|
|
|
|
# ---------------------------------------------------------------------------#
|
|
# ------------------------------ CLASSES -----------------------------------#
|
|
# ---------------------------------------------------------------------------#
|
|
|
|
class RawAquascatData:
|
|
""" This class loads Aquascat data from Aquascat text file (aqa.txt) or
|
|
Aquascat raw file (.aqa) located in <path>
|
|
|
|
Data file info
|
|
@ivar path: string, full file location
|
|
@ivar file_name: string, file name
|
|
@ivar folder_name : string, full folder location
|
|
|
|
Aquascat data
|
|
@ivar NumProfiles: int, number of profiles
|
|
@ivar ProfileRate: int, profile recording rate in Hz
|
|
@ivar Freq: array of float, frequencies in MHz
|
|
@ivar freqText: list of strings, frequencies in text format in MHz
|
|
@ivar At: array of float, transducer radii
|
|
@ivar NumCells: int, number of cells along one profile
|
|
@ivar V: 3 dimensional array of floats, ABS data
|
|
- dim. 0 : r (distance from the transducer)
|
|
- dim. 1 : freq (transducer frequency)
|
|
- dim. 2 : profile (over all the profiles recorded in a burst)
|
|
@ivar r: array, vector column of int, range of the cells in m from
|
|
the transducer
|
|
@ivar cellSize : float, size of the cells in m
|
|
@ivar TxPulseLength: float, length of the pulse in second
|
|
@ivar RxGain: array of float, gain at reception in dB
|
|
@ivar TxGain: array of float, gain at emission in dB
|
|
@ivar Average: int, number of pings averaged in one profile
|
|
@ivar BeamWidth: array of float, width of the acoustic beam in degree
|
|
@ivar Kt: array of float, Kt constant
|
|
@ivar PingRate: int, number of pings per second, in Hz
|
|
@ivar Regime: string, regime name given when using the Aquascat
|
|
"""
|
|
|
|
def __init__(self, p):
|
|
"""
|
|
@param p: string, full file location path
|
|
"""
|
|
|
|
# ___________________ File info _______________________________________#
|
|
|
|
# Source file full location path
|
|
self.path = p
|
|
# Extracting file name
|
|
self.file_name = os.path.basename(self.path)
|
|
# Extracting folder location
|
|
self.folder_name = os.path.dirname(self.path)
|
|
# Extracting beginning time info
|
|
self.date = datetime.datetime(1, 1, 1)
|
|
|
|
# ________________________ Data _______________________________________#
|
|
|
|
# Number of profiles
|
|
self.NumProfiles = 0
|
|
# Profile rate (Hz)
|
|
self.ProfileRate = 0
|
|
# Transducer frequencies (float)
|
|
self.Freq = np.array([])
|
|
# Transducer frequencies (string text in MHz)
|
|
self.freqText = []
|
|
# Transducer radii
|
|
self.At = np.array([])
|
|
# Number of cells
|
|
self.NumCells = 0
|
|
# ABS data
|
|
self.V = np.array([])
|
|
# Cell range data
|
|
self.r = np.array([])
|
|
# Cell size
|
|
self.cellSize = 0
|
|
# Pulse length in second
|
|
self.TxPulseLength = 0
|
|
# Reception gain in dB
|
|
self.RxGain = np.array([])
|
|
# Gain at emission
|
|
self.TxGain = np.array([])
|
|
# Averaging
|
|
self.Average = 0
|
|
# Beam width
|
|
self.BeamWidth = np.array([])
|
|
# Kt constant
|
|
self.Kt = np.array([])
|
|
# Ping rate
|
|
self.PingRate = 0
|
|
# Regime name
|
|
self.Regime = ''
|
|
|
|
# _____________ Reading and loading Aquascat Data _____________________#
|
|
if p:
|
|
self.date = datetime.datetime(year=np.int(self.file_name[0:4]),
|
|
month=np.int(self.file_name[4:6]),
|
|
day=np.int(self.file_name[6:8]),
|
|
hour=np.int(self.file_name[8:10]),
|
|
minute=np.int(self.file_name[10:12]),
|
|
second=np.int(self.file_name[12:14]))
|
|
|
|
if self.file_name[-7:] == 'aqa.txt':
|
|
self.load_txt_file()
|
|
|
|
if self.file_name[-4:] == '.aqa':
|
|
self.load_aqa_file()
|
|
|
|
def load_txt_file(self):
|
|
# Opening the file ("r" for "read")
|
|
src = open(self.path, "r")
|
|
|
|
# Line 1 : skipped
|
|
src.readline().rstrip('\n')
|
|
|
|
# Lines 2-3 : number of profiles
|
|
src.readline().rstrip('\n') # .rstrip('\n') for removing the "end line"
|
|
# character
|
|
|
|
temp = np.fromstring(src.readline().rstrip('\n'), dtype=float, sep=',')[
|
|
0] # string is converted
|
|
|
|
self.NumProfiles = int(temp)
|
|
|
|
# Lines 4-5 : profile rate in Hz
|
|
src.readline().rstrip('\n')
|
|
temp = np.fromstring(src.readline().rstrip('\n'), dtype=float, sep=',')[
|
|
0]
|
|
self.ProfileRate = temp
|
|
|
|
# Lines 6-7 : transducer frequencies
|
|
src.readline().rstrip('\n')
|
|
temp = np.fromstring(src.readline().rstrip('\n'), dtype=float, sep=',')
|
|
self.Freq = temp
|
|
# Creating also an array of strings with the frequencies in MHz
|
|
for i in range(len(self.Freq)):
|
|
self.freqText.append(str(self.Freq[i] / 1e6) + ' MHz')
|
|
|
|
# Lines 8-9 : transducer radii
|
|
src.readline().rstrip('\n')
|
|
temp = np.fromstring(src.readline().rstrip('\n'), dtype=float, sep=',')
|
|
self.At = temp
|
|
|
|
# Lines 10-14 : ABS data
|
|
src.readline().rstrip('\n')
|
|
|
|
# Reading the lines one by one
|
|
for f in range(len(self.Freq)):
|
|
temp = np.fromstring(src.readline().rstrip('\n'), dtype=float,
|
|
sep=',')
|
|
# Computing the number of cells
|
|
self.NumCells = int(len(temp) / self.NumProfiles)
|
|
|
|
# Converting the 1D numpy array to 2D array
|
|
# dim. 0 : r ; dim.1 : profile
|
|
# Have to take the transpose cause numpy reshape function
|
|
# fills the lines one by one, so one profile will fill one
|
|
# line and not one column as we want
|
|
temp = temp.reshape(self.NumProfiles, self.NumCells).T
|
|
|
|
# stacking the frequency to the already existing array so we get
|
|
# a 3D array where dim. 2 (last dimension) is corresponds to the
|
|
# frequency
|
|
if f == 0:
|
|
self.V = temp
|
|
else:
|
|
self.V = np.dstack((self.V, temp))
|
|
|
|
# Finally we exchange axes 1 and 2 so to get dim.1 = frequency and
|
|
# dim. 2 = profile
|
|
self.V = self.V.swapaxes(1, 2)
|
|
|
|
# Lines 15-19 : Mean ABS over the burst. Not very useful, skipped
|
|
src.readline().rstrip('\n')
|
|
|
|
# Reading the lines one by one
|
|
for f in range(len(self.Freq)):
|
|
src.readline().rstrip('\n')
|
|
|
|
# Lines 20-24 : cell range
|
|
src.readline().rstrip('\n')
|
|
# Only reading the first line, they are all the same
|
|
temp = np.fromstring(src.readline().rstrip('\n'), dtype=float, sep=',')
|
|
# Reshaping the array in a column vector
|
|
temp = temp.reshape(self.NumCells, 1)
|
|
self.r = temp
|
|
|
|
# Computing the cell size
|
|
self.cellSize = self.r[1, 0]-self.r[0, 0]
|
|
|
|
# Skipping the other r lines
|
|
for f in range(1, len(self.Freq)):
|
|
src.readline()
|
|
|
|
# Lines 25-26 : pulse length
|
|
src.readline().rstrip('\n')
|
|
temp = np.fromstring(src.readline().rstrip('\n'), dtype=float, sep=',')[
|
|
0]
|
|
self.TxPulseLength = temp
|
|
|
|
# Lines 27-28 : Gain at reception
|
|
src.readline().rstrip('\n')
|
|
temp = np.fromstring(src.readline().rstrip('\n'), dtype=float, sep=',')
|
|
self.RxGain = temp
|
|
|
|
# Lines 29-30 : Gain at emission
|
|
src.readline().rstrip('\n')
|
|
temp = np.fromstring(src.readline().rstrip('\n'), dtype=float, sep=',')
|
|
self.TxGain = temp
|
|
|
|
# Lines 31-32 : Averaging
|
|
src.readline().rstrip('\n')
|
|
temp = np.fromstring(src.readline().rstrip('\n'), dtype=float, sep=',')[
|
|
0]
|
|
self.Average = int(temp)
|
|
|
|
# Lines 33-34 : Beam width
|
|
src.readline().rstrip('\n')
|
|
temp = np.fromstring(src.readline().rstrip('\n'), dtype=float, sep=',')
|
|
self.BeamWidth = temp
|
|
|
|
# Lines 35-36 : Kt
|
|
src.readline().rstrip('\n')
|
|
temp = np.fromstring(src.readline().rstrip('\n'), dtype=float, sep=',')
|
|
self.Kt = temp
|
|
|
|
# Lines xxx : looking for ping rate
|
|
line = src.readline().rstrip('\n')
|
|
i = 0
|
|
while line != 'PingRate' and i < 100:
|
|
line = src.readline().rstrip('\n')
|
|
i += 1
|
|
temp = eval(src.readline().rstrip('\n'))
|
|
self.PingRate = temp
|
|
|
|
# Lines xxxx : session title
|
|
src.readline().rstrip('\n')
|
|
temp = src.readline().rstrip('\n')
|
|
# Formating the string
|
|
temp = temp.replace(',', '')
|
|
temp = temp.rstrip('\0')
|
|
self.Regime = temp
|
|
|
|
# Closing the file
|
|
src.close()
|
|
|
|
def load_aqa_file(self):
|
|
# -------- From G. Fromant 2017 - LEGI, Grenoble (France) ---------- #
|
|
|
|
# __________________ sub-functions ________________________________ #
|
|
|
|
def read_next_aquascat1000_header(f, file_size):
|
|
# Reading packet type
|
|
pkt_type_ = f.read(1)
|
|
pkt_type, = st.unpack("B", pkt_type_)
|
|
# Packet version: skipped
|
|
pkt_version_ = f.read(1)
|
|
st.unpack("B", pkt_version_)
|
|
# Reading packet size
|
|
pkt_size_ = f.read(2)
|
|
pkt_size, = st.unpack("H", pkt_size_)
|
|
# Packet checksum: skipped
|
|
pkt_checksum_ = f.read(2)
|
|
st.unpack("H", pkt_checksum_)
|
|
|
|
if f.tell() == file_size:
|
|
status = 0
|
|
else:
|
|
status = 1
|
|
|
|
return status, pkt_type, pkt_size
|
|
|
|
def find_aquascat1000_packet(f, type, file_size):
|
|
f.seek(0)
|
|
while f.tell() != file_size:
|
|
status, pkt_type, pkt_size = read_next_aquascat1000_header(
|
|
f, file_size)
|
|
if status == 0:
|
|
break
|
|
elif pkt_type == type:
|
|
break
|
|
else:
|
|
f.seek(2 * pkt_size, 1)
|
|
return status, pkt_size
|
|
|
|
# Opening the file
|
|
f = open(self.path, 'rb')
|
|
f.seek(0, 2)
|
|
file_size = f.tell()
|
|
f.seek(0)
|
|
ExpType = 'SINGLE'
|
|
|
|
# Read File Version from File
|
|
status, pkt_size = find_aquascat1000_packet(f, 19, file_size)
|
|
|
|
if status == 1:
|
|
sdata = [st.unpack("H", f.read(2))[0] for ui in
|
|
range(1, pkt_size + 1)]
|
|
FileVersionMajor = sdata[1]
|
|
FileVersionMinor = sdata[2]
|
|
else:
|
|
FileVersionMajor = 5
|
|
FileVersionMinor = 0
|
|
|
|
del sdata
|
|
|
|
# Read in the Burst Start Time Information
|
|
status, pkt_size = find_aquascat1000_packet(f, 54, file_size)
|
|
if status == 1:
|
|
sdata = [st.unpack("H", f.read(2))[0] for ui in range(1, 6 + 1)]
|
|
wake_source_ = f.read(2)
|
|
WakeSource, = st.unpack("H", wake_source_)
|
|
burst_number_ = f.read(4)
|
|
BurstNumber, = st.unpack("I", burst_number_)
|
|
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, pkt_size = find_aquascat1000_packet(f, 53, file_size)
|
|
|
|
if status == 0:
|
|
print('Error')
|
|
pass
|
|
else:
|
|
pkt_start_pos = f.tell()
|
|
skip_ = f.read(2)
|
|
st.unpack("H", skip_) # Size of Packet
|
|
st.unpack("H", f.read(2))
|
|
t1, = st.unpack("H", f.read(2))
|
|
t2, = st.unpack("H", f.read(2))
|
|
st.unpack("H", f.read(2))
|
|
LoggerType = [str(st.unpack("s", f.read(1)))[3] for ui in
|
|
range(1, 32 + 1)]
|
|
st.unpack("H", f.read(2))
|
|
NumAbsChannels, = st.unpack("H", f.read(2))
|
|
# Number of ABS channels the system supports
|
|
NumAuxChannels, = st.unpack("H", f.read(2))
|
|
# Number of AUX channels the system supports (8)
|
|
NumAbsTransducers, = st.unpack("H", f.read(2))
|
|
# Number of ABS Transducer Information that exist in personality table
|
|
|
|
# Battery capacity (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(pkt_start_pos + 112, 0)
|
|
# 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 = pkt_start_pos + 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)
|
|
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)
|
|
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)
|
|
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)
|
|
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)
|
|
[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)
|
|
AuxGainMin[i - 1][j - 1] = st.unpack("f", f.read(4))[
|
|
0] # Minimum value used in calibration data
|
|
AuxGainMax[i - 1][j - 1] = st.unpack("f", f.read(
|
|
4)) # Maximum value used in calibration data
|
|
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):
|
|
f.seek(pkt_start_pos + PtrToTransducerInfo * 2 + 200 * (i - 1),
|
|
0) # Move to the start of the ABS information
|
|
TTmp = [str(st.unpack("s", f.read(1)))[3] for ui in
|
|
range(1, 20 + 1)]
|
|
sss = ''
|
|
TTmp2 = sss.join(TTmp)
|
|
TTmp2 = TTmp2.strip("\\")
|
|
TransducerSerialNum = np.append(TransducerSerialNum, [TTmp2],
|
|
axis=0) # GV-10
|
|
st.unpack("H", f.read(2))
|
|
|
|
# Transducer frequency in Hz
|
|
TransducerFrequency = np.append(TransducerFrequency,
|
|
st.unpack("f", f.read(4)))
|
|
# Transducer radius in meters
|
|
TransducerRadius = np.append(TransducerRadius,
|
|
st.unpack("f", f.read(4)))
|
|
|
|
# Transducers beam width 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)]
|
|
|
|
# Transducer Kt (only if set in the personality)
|
|
TransducerKt = np.append(TransducerKt,
|
|
st.unpack("f", f.read(4)))
|
|
|
|
# ----------------------------------------------------------------------
|
|
# Read in the Regime (Logger Set-Up) Information
|
|
# ----------------------------------------------------------------------
|
|
|
|
# regime details
|
|
status, pkt_size = find_aquascat1000_packet(f, 21, file_size)
|
|
if 0 == status:
|
|
print('Regime Packet not found abort')
|
|
|
|
else:
|
|
pkt_start_pos = 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))
|
|
|
|
# This is a temp trick to correct---------------------|
|
|
# if NumAuxSamples!=0:
|
|
# 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)]
|
|
# Ping rate (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
|
|
AuxSampleRate = 0
|
|
if 0 == AuxPingDiv:
|
|
AuxNumSamples = 0
|
|
NumAuxChans = 0
|
|
else:
|
|
AuxSampleRate = PingRate / AuxPingDiv
|
|
AuxNumSamples = np.ceil(NumPings / AuxPingDiv) # GV-10
|
|
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(pkt_start_pos + 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 relevant 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])
|
|
# Transducer radius in m
|
|
AbsTransducerRadius = np.append(
|
|
AbsTransducerRadius, TransducerRadius[TransducerId - 1])
|
|
# Transducer beam width in degs
|
|
AbsTransducerBeamWidth = np.append(
|
|
AbsTransducerBeamWidth,
|
|
TransducerBeamWidth[TransducerId - 1])
|
|
AbsTransducerKt = np.append(AbsTransducerKt,
|
|
TransducerKt[TransducerId - 1])
|
|
|
|
AbsTxFrequency = np.append(AbsTxFrequency,
|
|
st.unpack("f", f.read(4))) # In Hz
|
|
AbsRxFrequency = np.append(AbsRxFrequency,
|
|
st.unpack("f", f.read(4))) # In Hz
|
|
|
|
# Pulse length in seconds
|
|
AbsTxPulseLength = np.append(
|
|
AbsTxPulseLength, st.unpack("f", f.read(4)))
|
|
Junk = st.unpack("f", f.read(4))
|
|
# Gain 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))
|
|
AbsPowerLevel = np.append(AbsPowerLevel, -20 * np.log10(
|
|
2 ** powerlevel)) # Power Level in dB relative to Maximum Power
|
|
# 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])
|
|
|
|
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) # in m --
|
|
|
|
# ---------------------------------------------------------------------
|
|
# Now deal with reading in the data
|
|
# ---------------------------------------------------------------------
|
|
|
|
# Allocate Memory for the Data
|
|
AuxData = np.zeros((int(AuxNumSamples), NumAuxChans))
|
|
AbsData = np.zeros((int(AbsNumBins[0]), int(AbsNumProfiles[0]),
|
|
int(NumAbsTimeSlots)), dtype=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, pkt_size = read_next_aquascat1000_header(f,
|
|
file_size)
|
|
|
|
if 1 == status:
|
|
|
|
if pktType == 22:
|
|
|
|
# switch(pktType)
|
|
|
|
# Case 22: Data were saved as magnitude values
|
|
# (normalize by 65536)
|
|
|
|
# case (22)
|
|
chan = st.unpack("H", f.read(2))[0] + 1
|
|
AbsIndex[chan - 1] = AbsIndex[
|
|
chan - 1] + 1 # Increase the Index
|
|
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:int(AbsNumBins[0]),
|
|
int(AbsIndex[chan - 1] - 1),
|
|
int(chan - 1)] = np.array(newList[:])
|
|
|
|
# Case 41: Data were saved as Complex values
|
|
# (normalize by 32768)
|
|
|
|
elif pktType == 41: # case (41)
|
|
chan = st.unpack("H", f.read(2))[0] + 1
|
|
AbsIndex[chan - 1] = AbsIndex[
|
|
chan - 1] + 1 # Increase the Index
|
|
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]),
|
|
int(chan - 1)] = np.array(sss)
|
|
|
|
# Case 46: Auxiliary Channel Data
|
|
|
|
elif pktType == 46: # case (46)
|
|
temp = st.unpack("H", f.read(2))[0] # Gain settings
|
|
AuxGain = [(temp & 3) + 1, ((temp >> 2) & 3) + 1,
|
|
((temp >> 4) & 3) + 1, ((temp >> 6) & 3) + 1, 1,
|
|
1, 1, 1]
|
|
Junk = st.unpack("H", f.read(2))[0] # Flags
|
|
AuxIndex = AuxIndex + 1
|
|
Tmp = [st.unpack("H", f.read(2))[0] for ui in
|
|
np.arange(1, NumAuxChans + 1)]
|
|
try:
|
|
AuxData[AuxIndex - 1, 0:NumAuxChans] = Tmp[:]
|
|
except IndexError:
|
|
pass
|
|
|
|
# Case 55: External Pressure Channel for upgraded ABS system -
|
|
# details to be provided
|
|
|
|
elif pktType == 55: # case (55)
|
|
chan = st.unpack("H", f.read(2))[0] + 1
|
|
SerIndex = SerIndex + 1 # Increase the Index
|
|
Tmp = [st.unpack("f", f.read(4))[0] for ui in
|
|
np.arange(1, 2 + 1)]
|
|
PressTempData[SerIndex - 1, :] = Tmp[:]
|
|
|
|
else: # otherwise
|
|
f.seek(pkt_size * 2, 1)
|
|
|
|
|
|
for i in range(0, NumAuxChans):
|
|
if not AuxGainCoeff[:][i]:
|
|
coeff = np.array(AuxGainCoeff[:][0])
|
|
AuxGainCoeff[:][i].append(coeff[0] * 0)
|
|
|
|
# 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
|
|
|
|
####TO BE DONE USING ACTUAL AUXDATA
|
|
if NumAuxChans == 0:
|
|
AuxData = []
|
|
|
|
# ---------------- Add by Adrien Vergne: filling the class ---------- #
|
|
# Number of profiles
|
|
self.NumProfiles = int(AbsNumProfiles[0])
|
|
# Profile rate
|
|
self.ProfileRate = AbsProfileRate[0]
|
|
# Frequencies
|
|
self.Freq = AbsTxFrequency
|
|
# Frequencies in text format
|
|
for i in range(len(self.Freq)):
|
|
self.freqText.append(str(self.Freq[i] / 1e6) + ' MHz')
|
|
# Transducer radius
|
|
self.At = AbsTransducerRadius
|
|
# Backscatter data
|
|
self.V = AbsData
|
|
# Changing axis order: axis 0 range, axis1 frequency, axis 2 profile
|
|
self.V = self.V.swapaxes(2, 1)
|
|
# Range cells
|
|
self.r = AbsBinRange[0]
|
|
self.r = np.reshape(self.r, (np.size(self.r), 1))
|
|
# Number of cells
|
|
self.NumCells = np.size(self.r)
|
|
# Cell size
|
|
self.cellSize = self.r[1, 0] - self.r[0, 0]
|
|
# Pulse length
|
|
self.TxPulseLength = AbsTxPulseLength[0]
|
|
# Received gain
|
|
self.RxGain = AbsStartingGain
|
|
# Emitted gain
|
|
self.TxGain = np.round(AbsPowerLevel, 3)
|
|
# Averaging
|
|
self.Average = AbsAverage[0]
|
|
# Beam width
|
|
self.BeamWidth = AbsTransducerBeamWidth
|
|
# Kt
|
|
self.Kt = AbsTransducerKt
|
|
# Ping rate
|
|
self.PingRate = PingRate
|
|
# Session title
|
|
self.Regime = SessionTitle
|
|
|
|
f.close # Close the *.aqa data file
|
|
|
|
|
|
class MeanAquascatProfile:
|
|
""" This class creates an object that contains several data of the mean
|
|
profile computed for one burst
|
|
@ivar file_name: string, name of the raw data file
|
|
@ivar Freq: list of float, frequencies
|
|
@ivar freqText: list of strings, freq in text in MHz
|
|
@ivar NumProfiles: int, number of profiles averaged
|
|
@ivar r: array of float, range cells for the profile
|
|
@ivar rMin: int, minimum range
|
|
@ivar rMax: int, maximum range
|
|
@ivar V: 2D array of float, ABS mean data
|
|
@ivar Average_method: string, 'simple' or 'quadratic'
|
|
@ivar Rx: list of received gain (one for each frequency)
|
|
@ivar Tx: list of emission gain (one for each frequency)
|
|
@ivar average : int, number of pings averaged per profile in the
|
|
initial burst
|
|
@ivar pulseLength: array of floats, length of the pulse in second
|
|
"""
|
|
|
|
def __init__(self):
|
|
|
|
# File name of the raw data
|
|
self.file_name = ''
|
|
# Frequencies
|
|
self.Freq = []
|
|
# Frequencies in text format
|
|
self.freqText = []
|
|
# Number of profiles averaged
|
|
self.NumProfiles = 0
|
|
# Range cells
|
|
self.r = np.array([])
|
|
# Minimum range
|
|
self.rMin = 0
|
|
# Maximum range
|
|
self.rMax = 0
|
|
# ABS data
|
|
self.V = np.array([])
|
|
# Averaging method
|
|
self.Average_method = ''
|
|
# Reception gain
|
|
self.Rx = []
|
|
# Emission gain
|
|
self.Tx = []
|
|
# number of pings averaged per profile in the intial burst
|
|
self.average = 0
|
|
# pulse length (s)
|
|
self.pulseLength = 0
|
|
|
|
def compute_mean(self, abs_raw, r_min=0, r_max=-1, mean='quadratic'):
|
|
"""This function computes a mean backscatter profile from the
|
|
raw data abs_raw, between r_min and r_max, using average method
|
|
'quadratic' or 'simple'"""
|
|
# Filling the profile medata
|
|
self.file_name = abs_raw.file_name
|
|
self.Freq = abs_raw.Freq
|
|
self.freqText = abs_raw.freqText
|
|
self.NumProfiles = abs_raw.NumProfiles
|
|
self.Rx = abs_raw.RxGain
|
|
self.Tx = abs_raw.TxGain
|
|
self.average = abs_raw.Average
|
|
self.pulseLength = abs_raw.TxPulseLength
|
|
|
|
# Filling r values
|
|
if r_min == 0:
|
|
ind_rmin = 0
|
|
else:
|
|
ind_rmin = np.where(abs_raw.r <= r_min)[0][-1]
|
|
if r_max == -1:
|
|
ind_rmax = abs_raw.r.size
|
|
else:
|
|
ind_rmax = np.where(abs_raw.r >= r_max)[0][0]
|
|
self.r = abs_raw.r[ind_rmin:ind_rmax + 1]
|
|
self.rMin = self.r[0][0]
|
|
self.rMax = self.r[-1][0]
|
|
|
|
# Filling V values
|
|
if mean == 'quadratic':
|
|
self.Average_method = 'quadratic'
|
|
self.V = np.sqrt(
|
|
np.nanmean(abs_raw.V[ind_rmin:ind_rmax + 1, :, :] ** 2,
|
|
axis=2))
|
|
elif mean == 'simple':
|
|
self.Average_method = 'simple'
|
|
self.V = np.nanmean(abs_raw.V[ind_rmin:ind_rmax + 1, :, :], axis=2)
|
|
else:
|
|
raise ValueError("mean should be 'quadratic' or 'simple'")
|
|
|
|
def write_txt_file(self, path):
|
|
""" This method allows to write the mean profile data in a
|
|
text file"""
|
|
|
|
# Opening the text file to write in
|
|
file = open(path, "w")
|
|
|
|
# Writing the parameters
|
|
file.write("AQUAscat mean ABS profile\n\n")
|
|
|
|
# Burst name
|
|
file.write("Initial burst name\n")
|
|
file.write(self.file_name + '\n\n')
|
|
|
|
# Frequencies
|
|
file.write('Freq\n')
|
|
for i in self.Freq[0:-1]:
|
|
file.write('%s, ' % i)
|
|
file.write(str(self.Freq[-1])+'\n\n')
|
|
|
|
# Frequencies in text format
|
|
file.write('Freq in text format\n')
|
|
for i in self.freqText[0:-1]:
|
|
file.write(i + ', ')
|
|
file.write(self.freqText[-1]+'\n\n')
|
|
|
|
# Number of profiles averaged
|
|
file.write('Number of profiles averaged\n')
|
|
file.write('%s \n\n' % self.NumProfiles)
|
|
|
|
# Range cells
|
|
file.write('r\n')
|
|
for i in self.r[0:-1]:
|
|
file.write('%s, ' % i[0])
|
|
file.write('%s' % self.r[-1][0])
|
|
file.write('\n\n')
|
|
|
|
# Minimum range
|
|
file.write('rMin\n')
|
|
file.write(str(self.rMin) + '\n\n')
|
|
|
|
# Maximum range
|
|
file.write('rMax\n')
|
|
file.write(str(self.rMax) + '\n\n')
|
|
|
|
# ABS data
|
|
file.write('V\n')
|
|
for f in range(len(self.Freq)):
|
|
for i in self.V[0:-1, f]:
|
|
file.write('%s, ' % i)
|
|
file.write(str(self.V[-1, f])+'\n')
|
|
file.write('\n')
|
|
|
|
# Reception gain
|
|
file.write('Rx gain\n')
|
|
for i in self.Rx[0:-1]:
|
|
file.write('%s, ' % i)
|
|
file.write(str(self.Rx[-1])+'\n\n')
|
|
|
|
# Transmission gain
|
|
file.write('Tx gain\n')
|
|
for i in self.Tx[0:-1]:
|
|
file.write('%s, ' % i)
|
|
file.write(str(self.Tx[-1])+'\n\n')
|
|
|
|
# Number of pings average
|
|
file.write("Nb pings averaged in 1 profile\n")
|
|
file.write(str(self.average) + '\n\n')
|
|
|
|
# Pulse length
|
|
file.write("Pulse length (s)\n")
|
|
file.write(str(self.pulseLength) + '\n\n')
|
|
|
|
file.close()
|
|
|
|
def load_txt_file(self, path):
|
|
""" This method allows to load the mean profile data from a
|
|
text file"""
|
|
|
|
# Opening the text file to write in
|
|
file = open(path, "r")
|
|
|
|
# Line 1 and 2: skipped
|
|
file.readline().rstrip('\n')
|
|
file.readline().rstrip('\n')
|
|
|
|
# Lines 3-4 : initial burst name
|
|
file.readline().rstrip('\n') # .rstrip('\n') for removing the "end
|
|
# line" character
|
|
temp = file.readline().rstrip('\n')
|
|
self.file_name = temp
|
|
|
|
# Lines 6-7 : frequencies
|
|
file.readline().rstrip('\n')
|
|
file.readline().rstrip('\n')
|
|
temp = np.fromstring(file.readline().rstrip('\n'), dtype=float, sep=',')
|
|
self.Freq = temp
|
|
|
|
# Lines 9-10 : frequencies in text format
|
|
file.readline().rstrip('\n')
|
|
file.readline().rstrip('\n')
|
|
temp = file.readline().rstrip('\n')
|
|
self.freqText = temp.split(sep=', ')
|
|
|
|
# Line 12 - 13 : Number of profiles averaged
|
|
file.readline().rstrip('\n')
|
|
file.readline().rstrip('\n')
|
|
temp = eval(file.readline().rstrip('\n'))
|
|
self.NumProfiles = temp
|
|
|
|
# Lines 15-16 : range data
|
|
file.readline().rstrip('\n')
|
|
file.readline().rstrip('\n')
|
|
temp = np.fromstring(file.readline().rstrip('\n'), dtype=float, sep=',')
|
|
# Reshaping the array in a column vector
|
|
temp = temp.reshape(len(temp), 1)
|
|
self.r = temp
|
|
|
|
# Lines 18-19 : minimum range
|
|
file.readline().rstrip('\n')
|
|
file.readline().rstrip('\n')
|
|
temp = eval(file.readline().rstrip('\n'))
|
|
self.rMin = temp
|
|
|
|
# Lines 21-22 : maximum range
|
|
file.readline().rstrip('\n')
|
|
file.readline().rstrip('\n')
|
|
temp = eval(file.readline().rstrip('\n'))
|
|
self.rMax = temp
|
|
|
|
# Lines 24 --- ABS data
|
|
file.readline().rstrip('\n')
|
|
file.readline().rstrip('\n')
|
|
# Reading the lines one by one
|
|
for f in range(len(self.Freq)):
|
|
temp = np.fromstring(file.readline().rstrip('\n'), dtype=float,
|
|
sep=', ')
|
|
|
|
if f == 0:
|
|
self.V = temp
|
|
else:
|
|
# Adding a column
|
|
self.V = np.c_[self.V, temp]
|
|
|
|
# Rx gain
|
|
file.readline().rstrip('\n')
|
|
file.readline().rstrip('\n')
|
|
temp = np.fromstring(file.readline().rstrip('\n'), dtype=float, sep=',')
|
|
self.Rx = temp
|
|
|
|
# Tx gain
|
|
file.readline().rstrip('\n')
|
|
file.readline().rstrip('\n')
|
|
temp = np.fromstring(file.readline().rstrip('\n'), dtype=float, sep=',')
|
|
self.Tx = temp
|
|
|
|
# Number of pings averaged
|
|
file.readline().rstrip('\n')
|
|
file.readline().rstrip('\n')
|
|
temp = np.fromstring(file.readline().rstrip('\n'), dtype=float, sep=',')
|
|
self.average = temp
|
|
|
|
# Pulse length
|
|
file.readline().rstrip('\n')
|
|
file.readline().rstrip('\n')
|
|
temp = np.fromstring(file.readline().rstrip('\n'), dtype=float, sep=',')
|
|
self.pulseLength = temp
|
|
|
|
file.close()
|
|
|
|
# ------------------------- Test --------------------------------------#
|
|
start_time = time.time()
|
|
|
|
if __name__ == "__main__":
|
|
# path1 = r'C:\Users\vergne\Documents\Donnees_aquascat\2017_juillet - experience cuve sediments fins\2017_07_19 - mercredi eau claire\20170719114700.aqa'
|
|
path1 = r'//home/brahim.moudjed/Documents/3 Software_Project/river_inversion_project/Data/Aquascat data test/20171213135800.aqa'
|
|
data1 = RawAquascatData(path1)
|
|
|
|
# path2 = r'C:\Users\vergne\Documents\Donnees_aquascat\2017_juillet - experience cuve sediments fins\2017_07_19 - mercredi eau claire\20170719114700.aqa.txt'
|
|
path2 = r'//home/brahim.moudjed/Documents/3 Software_Project/river_inversion_project/Data/Aquascat data test/20171213135800.txt'
|
|
data2 = RawAquascatData(path2)
|
|
|
|
print(data1.PingRate)
|
|
print(data2.PingRate)
|
|
|
|
print("Computational time: %.2f min" %((time.time() - start_time)/60) )
|