diff --git a/Model/AquascatDataLoader.py b/Model/AquascatDataLoader.py new file mode 100644 index 0000000..0598855 --- /dev/null +++ b/Model/AquascatDataLoader.py @@ -0,0 +1,1104 @@ +# -*- 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 + + 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((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) ) diff --git a/Model/TableModel.py b/Model/TableModel.py new file mode 100644 index 0000000..54ac9b7 --- /dev/null +++ b/Model/TableModel.py @@ -0,0 +1,39 @@ +from PyQt5.QtCore import Qt, QAbstractTableModel + + +class TableModel(QAbstractTableModel): + def __init__(self, data): + super(TableModel, self).__init__() + self._data = data + + def data(self, index, role): + if role == Qt.DisplayRole: + value = self._data.iloc[index.row(), index.column()] + # if role == Qt.TextAlignmentRole: + # value = self._data.iloc[index.row(), index.column()] + # if isinstance(value, int) or isinstance(value, float): + # return Qt.AlignVCenter + Qt.AlignRight + return str(value) + # if isinstance(value, float): + # # Render float to 2 dp + # return "%.2f" % value + # if isinstance(value, str): + # # Render strings with quotes + # return '"%s"' % value + + def rowCount(self, index): + # The length of the outer list. + return self._data.shape[0] + + def columnCount(self, index): + # The following takes the first sub-list, and returns + # the length (only works if all rows are an equal length) + return self._data.shape[1] + + def headerData(self, section, orientation, role): + # section is the index of the column/row. + if role == Qt.DisplayRole: + if orientation == Qt.Horizontal: + return str(self._data.columns[section]) + if orientation == Qt.Vertical: + return str(self._data.index[section]) diff --git a/Model/__init__.py b/Model/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/Model/acoustic_data_loader.py b/Model/acoustic_data_loader.py new file mode 100644 index 0000000..0259b97 --- /dev/null +++ b/Model/acoustic_data_loader.py @@ -0,0 +1,94 @@ +from Model.AquascatDataLoader import RawAquascatData + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from matplotlib.colors import LogNorm + +path_BS_raw_data = "/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/" \ + "Data/Acoustic_data/20180107123500.aqa" +path_noise_data = "/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/" \ + "Data/AcousticNoise_data/20180107121600.aqa" + +class AcousticDataLoader(): + + def __init__(self, path_BS_raw_data: str): + + self.path_BS_raw_data = path_BS_raw_data + + # --- Load Backscatter acoustic raw data with RawAquascatData class --- + + self._data_BS = RawAquascatData(self.path_BS_raw_data) + + self._BS_raw_data = self._data_BS.V + self._r = self._data_BS.r + self._freq = self._data_BS.Freq + self._freq_text = self._data_BS.freqText + self._time = np.array([t / self._data_BS.PingRate for t in range(self._data_BS.NumProfiles)]) + + self._date = self._data_BS.date.date() + self._hour = self._data_BS.date.time() + self._nb_profiles = self._data_BS.NumProfiles + self._nb_profiles_per_sec = self._data_BS.ProfileRate + self._nb_cells = self._data_BS.NumCells + self._cell_size = self._data_BS.cellSize + self._pulse_length = self._data_BS.TxPulseLength + self._nb_pings_per_sec = self._data_BS.PingRate + self._nb_pings_averaged_per_profile = self._data_BS.Average + self._kt = self._data_BS.Kt + self._gain_rx = self._data_BS.RxGain + self._gain_tx = self._data_BS.TxGain + + self._snr = np.array([]) + self._snr_reshape = np.array([]) + self._time_snr = np.array([]) + + # print(["BS - " + f for f in self._freq_text]) + # print(self._time.shape[0]*self._r.shape[0]*4) + + # fig, ax = plt.subplots(nrows=1, ncols=1) + # ax.pcolormesh(self._time, self._r, np.flipud(self._BS_raw_data[:, 1, :]), + # cmap='viridis', + # norm=LogNorm(vmin=1e-5, vmax=np.max(self._BS_raw_data[:, 0, :]))) # , shading='gouraud') + # ax.pcolormesh(range(self._BS_raw_data.shape[2]), range(self._BS_raw_data.shape[0]), self._BS_raw_data[:, 1, :], cmap='viridis', + # norm=LogNorm(vmin=1e-5, vmax=np.max(self._BS_raw_data[:, 0, :]))) # , shading='gouraud') + # plt.show() + + + # print(self.reshape_BS_raw_cross_section()[0, 0]) + # self.reshape_r() + # self.reshape_t() + + # dataframe = self.concatenate_data() + # print(dataframe) + + def reshape_BS_raw_cross_section(self): + BS_raw_cross_section = np.reshape(self._BS_raw_data, + (self._r.shape[0]*len(self._time), self._freq.shape[0]), + order="F") + return BS_raw_cross_section + + def reshape_r(self): + r = np.reshape(np.repeat(self._r, self._time.shape[0], axis=1), + self._r.shape[0]*self._time.shape[0], + order="F") + return r + + def reshape_t(self): + t = np.reshape(np.repeat(self._time, self._r.shape[0]), (self._time.shape[0]*self._r.shape[0], 1)) + return t + + def concatenate_data(self): + self.reshape_t() + self.reshape_BS_raw_cross_section() + # print(self.reshape_t().shape) + # print(se.lf.reshape_BS_raw_cross_section().shape) + df = pd.DataFrame(np.concatenate((self.reshape_t(), self.reshape_BS_raw_cross_section()), axis=1), + columns=["time"] + self._freq_text) + return df + + +# if __name__ == "__main__": +# AcousticDataLoader(path_BS_raw_data) + + diff --git a/Model/granulo_loader.py b/Model/granulo_loader.py new file mode 100644 index 0000000..80d51ef --- /dev/null +++ b/Model/granulo_loader.py @@ -0,0 +1,41 @@ +import numpy as np +import pandas as pd + + +class GranuloLoader: + + """ This class allows to load granulo data file """ + + def __init__(self, path_fine: str, path_sand: str): + + # --- Load fine sediments data file --- + self.path_fine = path_fine + self._data_fine = pd.read_excel(self.path_fine, engine="odf", header=0) + + self._y = np.array(self._data_fine.iloc[:, 0]) # distance from left bank (m) + self._z = np.array(self._data_fine.iloc[:, 1]) # depth (m) + self._r_grain = np.array(self._data_fine.columns.values)[4:] # grain radius (um) + + self._Ctot_fine = np.array(self._data_fine.iloc[:, 2]) # Total concentration (g/L) + self._D50_fine = np.array(self._data_fine.iloc[:, 3]) # median diameter (um) + self._frac_vol_fine = np.array(self._data_fine.iloc[:, 4:]) # Volume fraction (%) + + self._frac_vol_fine_cumul = np.cumsum(self._frac_vol_fine, axis=1) # Cumulated volume fraction (%) + + # --- Load sand sediments data file --- + self.path_sand = path_sand + self._data_sand = pd.read_excel(self.path_sand, engine="odf", header=0) + + self._Ctot_sand = np.array(self._data_sand.iloc[:, 2]) # Total concentration (g/L) + self._D50_sand = np.array(self._data_sand.iloc[:, 3]) # median diameter (um) + self._frac_vol_sand = np.array(self._data_sand.iloc[:, 4:]) # Volume fraction (%) + + self._frac_vol_sand_cumul = np.cumsum(self._frac_vol_sand, axis=1) # Cumulated volume fraction (%) + + +if __name__ == "__main__": + GranuloLoader("/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/Data/Granulo_data/" + "fine_sample_file.ods", + "/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/Data/Granulo_data/" + "fine_sample_file.ods") +