acoused/Model/AquascatDataLoader.py

2230 lines
88 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) )
# **********************************************************************************************************************
# **********************************************************************************************************************
# **********************************************************************************************************************
# -*- 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=int(self.file_name[0:4]),
month=int(self.file_name[4:6]),
day=int(self.file_name[6:8]),
hour=int(self.file_name[8:10]),
minute=int(self.file_name[10:12]),
second=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 readNextAQUAscat1000Header(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 = readNextAQUAscat1000Header(
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), int(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)) - 1
f.seek(0, 0)
# Now Read in all the Data
while (f.tell() != file_size): # 0 == feof(fid))
status, pktType, pktSize = readNextAQUAscat1000Header(f, file_size)
if 1 == status:
if pktType == 22:
chan = st.unpack("H", f.read(2))[0] + 1
# Increase the Index
AbsIndex[chan - 1] = AbsIndex[chan - 1] + 1
Tmp = [st.unpack(
"H", f.read(2))[0] for ui in np.arange(
1, AbsNumBins[chan - 1] + 1)]
newList = [x / 65536 for x in Tmp]
AbsData[0: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)
# case (41)
elif pktType == 41:
chan = st.unpack("H", f.read(2))[0] + 1
# Increase the Index
AbsIndex[chan - 1] = AbsIndex[chan - 1] + 1
Tmp = [st.unpack(
"h", f.read(2))[0] for ui in np.arange(
1, 2 * AbsNumBins[chan - 1] + 1)]
sss = [x / 32768 for x in Tmp]
AbsData[:, int(
AbsIndex[chan - 1]), int(chan - 1)] = np.abs(
np.array(sss[0::2]) + complex(0, 1) * np.array(sss[1::2]))
# Case 46: Auxiliary Channel Data
# case (46)
elif pktType == 46:
# Gain settings
temp = st.unpack("H", f.read(2))[0]
AuxGain = [(temp & 3) + 1, ((temp >> 2) & 3) + 1, (
(temp >> 4) & 3) + 1, (
(temp >> 6) & 3) + 1, 1, 1, 1, 1]
# Flags
Junk = st.unpack("H", f.read(2))[0]
AuxIndex = AuxIndex + 1
Tmp = [st.unpack(
"H", f.read(2))[0] for ui in np.arange(
1, NumAuxChans + 1)]
if AuxIndex == AuxNumSamples + 1:
AuxIndex = AuxIndex - 1
pass
else:
AuxData[AuxIndex - 1, 0:NumAuxChans] = Tmp[:]
# AuxData[AuxIndex-1,0:NumAuxChans] = Tmp[:]
# Case 55: External Pressure Channel for upgraded
# ABS system - details to be provided
# case (55)
elif pktType == 55:
chan = st.unpack("H", f.read(2))[0] + 1
# Increase the Index
SerIndex = SerIndex + 1
Tmp = [st.unpack(
"f", f.read(4))[0] for ui in np.arange(1, 2 + 1)]
PressTempData[SerIndex - 1, :] = Tmp[:]
else: # otherwise
f.seek(pktSize * 2, 1)
print('Summary of Data Imported')
for j in range(1, NumAbsTimeSlots + 1):
print('Timeslot {} Total Profiles = {}'.format(j - 1, AbsIndex[j - 1]))
print('Total Aux Samples = {}'.format(AuxIndex))
for i in range(0, NumAuxChans):
if not AuxGainCoeff[:][i]:
coeff = np.array(AuxGainCoeff[:][0])
AuxGainCoeff[:][i].append(coeff[0] * 0)
# coeff = AuxGainCoeff(:,AuxGain(i),i)
# Matlab orders polynomial coeff opposite to Aquatec
coeff = np.array(AuxGainCoeff[:][i])
coeff = coeff.reshape(1, 5)
coeff = np.fliplr(coeff)
coeff = coeff.reshape(5)
AuxData[:, i] = np.polyval(coeff, AuxData[:, i])
# Need to apply calibration coefficients for the Aux Channels, now that the
# gain is fixed, auto gain is not supported. Therefore will use last value
if NumAuxChans == 0:
AuxData = []
# ---------------- 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))