diff --git a/Model/GrainSizeTools.py b/Model/GrainSizeTools.py new file mode 100644 index 0000000..0783f82 --- /dev/null +++ b/Model/GrainSizeTools.py @@ -0,0 +1,269 @@ +# -*- coding: utf-8 -*- +""" +Created on Thur Jan 10 2019 + +@author: Adrien Vergne +""" + +############################################################################### +# # +# GRAIN SIZE TOOLS # +# # +############################################################################### + +# <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +# Adrien VERGNE - Irstea Lyon - 2019 +# Program Python 3.5 +# <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> + +# This script is designed for loading PSD data from Malvern Mastersizer +# Excek file output, as well as for fitting Gausian mixture models to +# the measured PSD + +########################### Loading libraries ######################### +import numpy as np # Numpy (array and matrix computing) +import pandas as pd # DataFrames +from math import sqrt, pi +from scipy.interpolate import Akima1DInterpolator # an interpolator +from sklearn.mixture import GaussianMixture # Gaussian mixture model (demodulation) +import matplotlib.pyplot as plt # matplotlib scientific plotting + + + +########################### Loading Malvern PSD data ######################### +class load_malvern: + """This class loads Malvern particle size data from an Excel file""" + def __init__(self, folder, file): + # Excel file + self.file_folder = folder + self.file_name = file + self.file_path = folder + '/' + file + + # Data + self.raw_data = None + self.size_inf = None + self.size_sup = None + self.size_center = None + self.size_proba = None + self.D10 = 0 + self.D50 = 0 + self.D90 = 0 + + # Loading data + self.load() + + def load(self): + self.raw_data = pd.read_excel(self.file_path, 0) + self.size_inf = np.array(self.raw_data['Size x'][1:]) * 1e-6 / 2 + self.size_sup = np.array(self.raw_data['Size y'][1:]) * 1e-6 / 2 + self.size_center = (self.size_inf + self.size_sup)/2 + self.size_proba = np.array(self.raw_data['% In Size'][1:])/100 + self.D10 = self.raw_data['µm Diameter'][2] + self.D50 = self.raw_data['µm Diameter'][5] + self.D90 = self.raw_data['µm Diameter'][8] + + +########################### Fitting GMM ######################### +def mix_gaussian_model(x, m_list, s_list, w_list): + """This function computes a sum of gaussians""" + res = np.zeros(np.shape(x)) + for i in range(len(m_list)): + res += w_list[i] * (1 / sqrt(2 * pi * s_list[i]**2) * np.exp(-(x - m_list[i])**2 / (2 * s_list[i]**2))) + return res + +def simple_gaussian(x, m, s, w): + """This function computes a simple guaussian""" + + res = w * (1 / sqrt(2 * pi * s**2) * np.exp(-(x - m)**2 / (2 * s**2))) + return res + +class demodul_granulo: + """This class demodulates a grain size distribution + INPUTS: + @size_classes: center of each size class (radius in meters) + @proba_cumul_vol: cumulative volume probability of each size class + @interpolate_min_size: start (i.e. minimum size) of cumulative volume probability interpolation (radius in meters) + @interpolate_max_size: end (i.e. maximum size) of cumulative volume probability interpolation (radius in meters) + """ + + def __init__(self, size_classes, proba_cumul_vol, interpolate_min_size, interpolate_max_size): + ###### Data containers ###### + # ~~~~~~~~ Inputs ~~~~~~~~~~# + # size classes (radius in meters) + self.size_classes = size_classes + # volume cumulative probability (from 0 to 1) + self.proba_cumul_vol = proba_cumul_vol + # start (i.e. minimum size) of cumulative volume probability interpolation + self.interpolate_min_size = interpolate_min_size + # end (i.e. maximum size) of cumulative volume probability interpolation + self.interpolate_max_size = interpolate_max_size + + # ~~~~~~~~ Outputs ~~~~~~~~~~# + # Resampled size classes in logscale + self.sizes_log_resampled = np.array([]) + # Interpolated cumulative probability + self.cumul_interpolated = np.array([]) + # Size class volume probability (no cumul) from resampled size classes + self.resampled_proba_vol = np.array([]) + # Size class number probability (no cumul) from resampled size classes + self.resampled_proba_num = np.array([]) + # List of Gaussian mixture models (from 1 to 10 modes) + self.gaussian_mixture_model_list = [] + # List of outputs data for each model (from 1 to 10 modes) --> see class "model_XX_modes_values" + self.demodul_data_list = [] + + # Computing + self.interpolate_data() + self.demodulate() + + + class model_XX_modes_values: + """This sub-class is a data container for the outputs of a Gaussian Mixture Models of modes""" + def __init__(self, x, nb_modes, mu_list, sigma_list, w_list): + # Number of modes + self.nb_modes = nb_modes + # List of mode centers (mu) + self.mu_list = mu_list + # List of mode scalling parameters (sigma) + self.sigma_list = sigma_list + # List of mode weights (w) + self.w_list = w_list + + # Gaussian Mixture model values + self.model_total = mix_gaussian_model(x, self.mu_list, self.sigma_list, self.w_list) / \ + np.sum(mix_gaussian_model(x, self.mu_list, self.sigma_list, self.w_list)) + # list of single modes values + self.modes_list = [] + for i in range(self.nb_modes): + # Computing the mode values + self.modes_list.append(simple_gaussian(x, self.mu_list[i], self.sigma_list[i], self.w_list[i]) / + np.sum(mix_gaussian_model(x, self.mu_list, self.sigma_list, self.w_list))) + + + def interpolate_data(self): + """This method interpolates the cumulative probability in logscale, between the specified min and max size""" + # Phi steps + log_scale = 0.005 + # Computing maximum phi (minimum diameter) + log_min = np.log(self.interpolate_min_size) + # Computing minimum phi (maximum diameter) + log_max = np.log(self.interpolate_max_size) + # Creating an array of equally spaced phi + self.sizes_log_resampled = np.arange(log_min, log_max, log_scale) + # Akima interpolation + f_logsize = Akima1DInterpolator(np.log(self.size_classes), self.proba_cumul_vol) + self.cumul_interpolated = f_logsize(self.sizes_log_resampled) + + # Computing volume and number probability from interpolated data + self.resampled_proba_vol = np.diff(self.cumul_interpolated) + # Computing number probability, + ss = np.sum(self.resampled_proba_vol /(np.exp(self.sizes_log_resampled[:-1]))**3) + self.resampled_proba_num = (self.resampled_proba_vol / (np.exp(self.sizes_log_resampled[:-1]))**3) / ss + + + def demodulate(self): + """This function demodulates the interpolated probability distribution""" + # Sampling the size classes. Need to normalize the probability to get sum = 1.0 + ech_demod = np.random.choice(self.sizes_log_resampled[:-1], (10000,1), + p=(self.resampled_proba_vol / np.sum(self.resampled_proba_vol))) + # fitting models from 1 to 10 components + N = np.arange(1, 11) + self.gaussian_mixture_model_list = [None for i in range(len(N))] + + for i in range(len(N)): + self.gaussian_mixture_model_list[i] = GaussianMixture(n_components=N[i], + covariance_type='full').fit(ech_demod) + mu_list = [m[0] for m in self.gaussian_mixture_model_list[i].means_] + s_list = [sqrt(s[0][0]) for s in self.gaussian_mixture_model_list[i].covariances_] + w_list = [w for w in self.gaussian_mixture_model_list[i].weights_] + + self.demodul_data_list.append(self.model_XX_modes_values(self.sizes_log_resampled, i + 1, mu_list, s_list, w_list)) + + def plot_interpolation(self): + """This method plots the cumulative probability interpolation""" + # ------- Plotting interpolation----------- # + fig,(ax)=plt.subplots(1, 2, figsize = (11,4), dpi=100) + + # Data + ax[0].plot(np.log(self.size_classes), self.proba_cumul_vol, linestyle='', marker='+') + ax[1].plot(np.log(self.size_classes), self.proba_cumul_vol, linestyle='', marker='+') + + # Interpolation + ax[0].plot(self.sizes_log_resampled, self.cumul_interpolated, color='red') + ax[1].plot(self.sizes_log_resampled, self.cumul_interpolated, color='red') + + ax[0].set_xlabel("Log(a)") + ax[0].set_ylabel("Cumulative fraction") + + ax[1].set_yscale('log') + ax[1].set_ylim((1e-6, 2)) + ax[1].set_xlabel("Log(a)") + ax[1].set_ylabel("Cumulative fraction") + ax[1].set_title('Double log scale') + + plt.show() + + + def print_mode_data(self, n): + """This method shows the parameters of the Gaussian Mixture Model of n modes""" + if (n <= 0) or (n >10): + print('Error: number of modes should be between 1 and 10') + else: + print('##### Number of modes: %d #####' %n) + model = self.gaussian_mixture_model_list[n - 1] + mode_center_txt = '' + sigma_txt = '' + weights_txt = '' + sorted_indices = [] + for i in range(n): + sorted_indices.append(np.where(model.means_ == np.sort(model.means_, axis=0)[i])[0][0]) + mode_center_txt = '%s %2.2f um' %(mode_center_txt, np.exp(model.means_[sorted_indices[i]])*1e6) + sigma_txt = '%s %2.2f' %(sigma_txt, sqrt(model.covariances_[sorted_indices[i]][0][0])) + weights_txt = '%s %2.1f %%' %(weights_txt, model.weights_[sorted_indices[i]]*100) + print('Centers (radius):', mode_center_txt) + print('Sigmas:', sigma_txt) + print('Weigths:', weights_txt) + + def plot_modes(self, n): + """This method plots the Gaussian Mixture Model of n components""" + if (n <= 0) or (n >10): + print('Error: number of modes should be between 1 and 10') + else: + # ---------------- Plotting ----------- # + fig,(ax)=plt.subplots(2, 1, figsize = (9,6), dpi=100) + # Data + ax[0].plot(np.exp(self.sizes_log_resampled[0:-1])*1e6, self.resampled_proba_vol, + label='Data - (proba vol.)', ls='--', linewidth=2, color='black') + ax[1].plot(np.exp(self.sizes_log_resampled[0:-1])*1e6, self.resampled_proba_vol, + label='Data - (proba vol.)', ls='--', linewidth=2, color='black') + + # Gaussian mixture model + ax[0].plot(np.exp(self.sizes_log_resampled)*1e6, + self.demodul_data_list[n - 1].model_total, color='red', label='Model (%d modes)' %(n)) + ax[1].plot(np.exp(self.sizes_log_resampled)*1e6, + self.demodul_data_list[n - 1].model_total, color='red', label='Model (%d modes)' %(n)) + + # Modes + for i in range(n): + sorted_indice = np.where(self.gaussian_mixture_model_list[n - 1].means_ == + np.sort(self.gaussian_mixture_model_list[n - 1].means_, axis=0)[i])[0][0] + ax[0].plot(np.exp(self.sizes_log_resampled)*1e6, + self.demodul_data_list[n - 1].modes_list[sorted_indice]/np.sum(self.demodul_data_list[n - 1].model_total), + label='Mode %d' %(i+1), linestyle=':') + + ax[1].plot(np.exp(self.sizes_log_resampled)*1e6, + self.demodul_data_list[n - 1].modes_list[sorted_indice]/np.sum(self.demodul_data_list[n - 1].model_total), + label='Mode %d' %(i+1), linestyle=':') + # Axes properties + ax[0].legend() + ax[0].set_xlabel('Radius (um)', fontsize=13, weight='bold') + ax[0].set_ylabel('Volume fraction', fontsize=13, weight='bold') + ax[0].set_title('Gaussian Mixture Model -- %d modes' %(n), fontsize=14, weight='bold') + + ax[1].legend(loc=0) + ax[1].set_xlabel('Radius (um)', fontsize=13, weight='bold') + ax[1].set_ylabel('Volume fraction', fontsize=13, weight='bold') + ax[1].set_xscale('log') + fig.tight_layout() + + plt.show() \ No newline at end of file diff --git a/Model/acoustic_data_loader_UBSediFlow.py b/Model/acoustic_data_loader_UBSediFlow.py index 95f38c9..b5c3487 100644 --- a/Model/acoustic_data_loader_UBSediFlow.py +++ b/Model/acoustic_data_loader_UBSediFlow.py @@ -29,6 +29,9 @@ from Model.udt_extract.raw_extract import raw_extract # ("/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/Data/APAVER_2021/" # "Rhone_20210519/Rhone_20210519/record/") +# path_BS_raw_data0 = ("/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/Data/Raw_data_udt/") +# filename0 = "raw_20210519_130000.udt" + # filename = "raw_20210519_115128.udt" # "raw_20210526_153310.udt" @@ -174,6 +177,13 @@ class AcousticDataLoaderUBSediFlow: if config == 1: BS_data = np.array([data_us_dicts[config][channel]['echo_avg_profile']['data'][0]]) # print("BS_data shape ", BS_data.shape) + print("******************************") + date_list = [np.abs(datetime.datetime(2021, 5, 19, 13, 10, 00).timestamp() + - date.timestamp()) for date in data_us_dicts[config][channel]['echo_avg_profile']['time']] + print(np.where(date_list == np.min(date_list))) + # == datetime.datetime(2021, 5, 19, 14, 10, 2, 644000)) + print("******************************") + for i in range(self._time.shape[1]): BS_data = np.append(BS_data, np.array([data_us_dicts[config][channel]['echo_avg_profile']['data'][i]]), @@ -335,8 +345,8 @@ class AcousticDataLoaderUBSediFlow: # --- Plot Backscatter US data --- # fig, ax = plt.subplots(nrows=1, ncols=1, layout="constrained") - # pcm = ax.pcolormesh(self._time[0, :], -self._r[0, :], np.log(self._BS_raw_data[0, :, :]), - # cmap='Blues')#, shading='gouraud') + # pcm = ax.pcolormesh(self._time[0, :], -self._r[0, 2:], np.log(self._BS_raw_data[0, 2:, :]), + # cmap='plasma')#, shading='gouraud') # # pcm = ax.pcolormesh(list(range(self._BS_raw_data.shape[2])), list(range(self._BS_raw_data.shape[1])), # # np.log(self._BS_raw_data[0, :, :]), # # cmap='Blues') # , shading='gouraud') @@ -370,7 +380,13 @@ class AcousticDataLoaderUBSediFlow: # fig1, ax1 = plt.subplots(nrows=1, ncols=1, layout="constrained") # pcm1 = ax1.pcolormesh(self._time[0, :], -self._r[0, :], np.log(BS_smooth[:, :]), cmap='Blues') # fig1.colorbar(pcm1, ax=ax1, shrink=1, location='right') + + # print("find value in depth", np.where(np.abs(self._r - 3.3) == np.min(np.abs(self._r - 3.3)))) # + # fig, ax = plt.subplots(nrows=1, ncols=1, layout="constrained") + # ax.plot(self._r[0, :], self._BS_raw_data[0, :, 766]) + # plt.show() + # # --- Plot vertical profile for bottom detection --- # n = 60 # t0 = 200 diff --git a/Model/acoustic_inversion_method_high_concentration.py b/Model/acoustic_inversion_method_high_concentration.py index dc5868b..cefb418 100644 --- a/Model/acoustic_inversion_method_high_concentration.py +++ b/Model/acoustic_inversion_method_high_concentration.py @@ -1,5 +1,7 @@ +import matplotlib.pyplot as plt import numpy as np import settings as stg +from Model.GrainSizeTools import demodul_granulo, mix_gaussian_model class AcousticInversionMethodHighConcentration(): @@ -23,14 +25,14 @@ class AcousticInversionMethodHighConcentration(): return C # -------- Computing water attenuation coefficient ----------- # - def water_attenuation(self, freq1, freq2, T): + def water_attenuation(self, freq, T): """Computing attenuation from François and Garrison 1982""" if T > 20: alpha = (3.964 * 1e-4 - 1.146 * 1e-5 * T + 1.45 * 1e-7 * T ** 2 - 6.5 * 1e-10 * T ** 3) * 1e-3 * \ - (np.log(10) / 20) * (np.array([freq1, freq2]) * 1e-3) ** 2 + (np.log(10) / 20) * (freq * 1e-3) ** 2 else: alpha = (4.937 * 1e-4 - 2.59 * 1e-5 * T + 9.11 * 1e-7 * T ** 2 - 1.5 * 1e-8 * T ** 3) * 1e-3 * \ - (np.log(10) / 20) * (np.array([freq1, freq2]) * 1e-3) ** 2 + (np.log(10) / 20) * (freq * 1e-3) ** 2 return alpha # ---------- Conmpute FBC ---------- @@ -44,22 +46,90 @@ class AcousticInversionMethodHighConcentration(): # np.log(2 * self.alpha_w_function(self.Freq[f], self.temperature) * self.r_3D[:, f, :]) # return FCB - # ------------- Computing ks ------------- # + # --- Gaussian mixture --- + def compute_particle_size_distribution_in_number_of_particles(self, num_sample, r_grain, frac_vol_cumul): + min_demodul = 1e-6 + max_demodul = 500e-6 - def ks(self, ind): - ks0 = np.array([0.0446681, 0.11490388, 0.35937832, 2.5025668]) - ks = ks0[ind] + sample_demodul = demodul_granulo(r_grain, frac_vol_cumul[num_sample, :], min_demodul, max_demodul) + + resampled_log_array = np.log(np.logspace(-10, -2, 3000)) + proba_vol_demodul = mix_gaussian_model(resampled_log_array, + sample_demodul.demodul_data_list[2].mu_list, + sample_demodul.demodul_data_list[2].sigma_list, + sample_demodul.demodul_data_list[2].w_list) + + # N_modes = 3 + # sample_demodul.print_mode_data(N_modes) + # sample_demodul.plot_interpolation() + # sample_demodul.plot_modes(N_modes) + + print(f"mu_list : {sample_demodul.demodul_data_list[3 - 1].mu_list}") + print(f"sigma_list : {sample_demodul.demodul_data_list[3 - 1].sigma_list}") + print(f"w_list : {sample_demodul.demodul_data_list[3 - 1].w_list}") + + proba_vol_demodul = proba_vol_demodul / np.sum(proba_vol_demodul) + ss = np.sum(proba_vol_demodul / np.exp(resampled_log_array) ** 3) + proba_num = proba_vol_demodul / np.exp(resampled_log_array) ** 3 / ss + + return proba_num + + # ------------- Computing ks ------------- # + def form_factor_function_MoateThorne2012(self, a, freq, C=1500): + """This function computes the form factor based on the equation of + Moate and Thorne (2012)""" + # computing the wave number + k = 2 * np.pi * freq / C + x = k * a + f = (x ** 2 * (1 - 0.25 * np.exp(-((x - 1.5) / 0.35) ** 2)) * (1 + 0.6 * np.exp(-((x - 2.9) / 1.15) ** 2))) / ( + 42 + 28 * x ** 2) + print(f"form factor = {f}") + return f + + def ks(self, freq, pdf): + # --- Calcul de la fonction de form --- + # form_factor = self.form_factor_function_MoateThorne2012(a, freq) + # print(f"form_factor shape = {form_factor}") + # print(f"form_factor = {form_factor}") + + #--- Particle size distribution --- + proba_num = ( + self.compute_particle_size_distribution_in_number_of_particles( + num_sample=2, r_grain=stg.radius_grain, frac_vol_cumul=stg.frac_vol_sand_cumul)) + + print(f"proba_num : {proba_num}") + + # --- Compute k_s by dividing two integrals --- + resampled_log_array = np.log(np.logspace(-10, -2, 3000)) + a2f2pdf = 0 + a3pdf = 0 + for i in range(len(resampled_log_array)): + a = np.exp(resampled_log_array)[i] + a2f2pdf += a**2 * self.form_factor_function_MoateThorne2012(a, freq)**2 * proba_num[i] + a3pdf += a**3 * proba_num[i] + + print(f"a2f2pdf = {a2f2pdf}") + print(f"a3pdf = {a3pdf}") + + ks = np.sqrt(a2f2pdf / a3pdf) + + # ks = np.array([0.04452077, 0.11415143, 0.35533713, 2.47960051]) + # ks = ks0[ind] return ks # ------------- Computing sv ------------- # - - def sv(self): - pass + def sv(self, ks, M_sand): + # print(f"ks = {ks}") + # print(f"M_sand = {M_sand}") + sv = (3 / (16 * np.pi)) * (ks ** 2) * M_sand + # sv = np.full((stg.r.shape[1], stg.t.shape[1]), sv0) + return sv # ------------- Computing X ------------- # - def X_exponent(self, ind): - X0 = [3.688664080521131, 3.451418735488537, 0, 3.276577078823426, 0, 0] - X = X0[ind] + def X_exponent(self, freq1, freq2, sv_freq1, sv_freq2): + # X0 = [3.450428714146802, 3.276478927777019, 3.6864638665972893, 0] + # X = X0[ind] + X = np.log(sv_freq1 / sv_freq2) / np.log(freq1 / freq2) return X # -------------- Computing Kt -------------- # @@ -78,34 +148,160 @@ class AcousticInversionMethodHighConcentration(): # Gain gain = 10 ** ((RxGain + TxGain) / 20) # Computing Kt - kt0 = kt_ref * gain * np.sqrt(tau * cel / (tau_ref * c_ref)) # 1D numpy array - kt = np.reshape(kt0, (1, 2)) # convert to 2d numpy array to compute J_cross_section + kt = kt_ref * gain * np.sqrt(tau * cel / (tau_ref * c_ref)) # 1D numpy array + # kt = np.reshape(kt0, (1, 2)) # convert to 2d numpy array to compute J_cross_section + # print(f"kt = {kt}") + # kt_2D = np.repeat(np.array([kt]), stg.r.shape[1], axis=0) + # print("kt 2D ", kt_2D) + # print("kt 2D shape ", kt_2D.shape) + # # kt_3D = np.zeros((kt_2D.shape[1], kt_2D.shape[0], stg.t.shape[1])) + # # for k in range(kt_2D.shape[1]): + # # kt_3D[k, :, :] = np.repeat(kt_2D, stg.t.shape[1], axis=1)[:, k * stg.t.shape[1]:(k + 1) * stg.t.shape[1]] + # kt_3D = np.repeat(kt_2D.transpose()[:, :, np.newaxis], stg.t.shape[1], axis=2) + # # print("kt 3D ", kt_3D) + # print("kt 3D shape ", kt_3D.shape) + return kt # ------------- Computing J_cross_section ------------- # def j_cross_section(self, BS, r2D, kt): - J_cross_section = np.zeros((2, BS.shape[0], BS.shape[2])) # 2 because it's a pair of frequencies - for k in range(2): - J_cross_section[k, :, :] = (3 / (16 * np.pi)) * ((BS[:, k, :]**2 * r2D**2) / kt[k, :, :]**2) - # J_cross_section[J_cross_section == 0] = np.nan - print("compute j_cross_section finished") + # J_cross_section = np.zeros((1, BS.shape[1], BS.shape[2])) # 2 because it's a pair of frequencies + # print("BS.shape", BS.shape) + # print("r2D.shape", r2D.shape) + # print("kt.shape", kt.shape) + + # if stg.ABS_name == "Aquascat 1000R": + # print("--------------------------------") + # print("BS : ", BS) + # print("BS min : ", np.nanmin(BS)) + # print("BS max : ", np.nanmax(BS)) + # print("r2D : ", r2D) + print("kt shape : ", kt.shape) + print("kt : ", kt) + # print("--------------------------------") + # for k in range(1): + # J_cross_section[k, :, :] = (3 / (16 * np.pi)) * ((BS[k, :, :]**2 * r2D[k, :, :]**2) / kt[k, :, :]**2) + J_cross_section = (3 / (16 * np.pi)) * ((BS**2 * r2D**2) / kt**2) + # J_cross_section[J_cross_section == 0] = np.nan + # print("J_cross_section.shape", J_cross_section.shape) + # elif stg.ABS_name == "UB-SediFlow": + # for k in range(1): + # J_cross_section[k, :, :] = (3 / (16 * np.pi)) * ((BS[k, :, :]**2 * r2D[0, :, :]**2) / kt[k, :, :]**2) + # print("compute j_cross_section finished") return J_cross_section # ------------- Computing alpha_s ------------- # + def alpha_s(self, sv, j_cross_section, depth, alpha_w): + alpha_s = (np.log(sv / j_cross_section) / (4 * depth)) - alpha_w + print("----------------------------") + print(f"sv = {sv}") + print(f"j_cross_section = {j_cross_section}") + print(f"depth = {depth}") + print(f"alpha_w = {alpha_w}") + print(f"(np.log(sv / j_cross_section) / (4 * depth)) = {(np.log(sv / j_cross_section) / (4 * depth))}") + print(f"alpha_s {alpha_s}") - def alpha_s(self): - pass + return alpha_s # ------------- Computing interpolation of fine SSC data obtained from water sampling ------------- # ------------- collected at various depth in the vertical sample ------------- - def M_profile_SCC_fine_interpolated(self): - pass + def M_profile_SCC_fine_interpolated(self, sample_depth, M_profile, range_cells, r_bottom): + res = np.zeros((len(range_cells),)) * np.nan + for i in range(len(M_profile) - 1): + # print(f"i = {i}") + r_ini = sample_depth[i] + # print(f"r_ini = {r_ini}") + c_ini = M_profile[i] + # print(f"c_ini = {c_ini}") + r_end = sample_depth[i + 1] + # print(f"r_end = {r_end}") + c_end = M_profile[i + 1] + # print(f"c_end = {c_end}") + + # Computing the linear equation + a = (c_end - c_ini) / (r_end - r_ini) + # print(f"a = {a}") + b = c_ini - a * r_ini + # print(f"b = {b}") + + # Finding the indices of r_ini and r_end in the interpolated array + # print(f"range_cells = {range_cells}") + loc = (range_cells >= r_ini) * (range_cells < r_end) + # print(f"loc = {loc}") + # Filling the array with interpolation values + res[loc] = range_cells[loc] * a + b + # print(f"res = {res}") + + # Filling first and last values + i = 0 + while np.isnan(res[i]): + res[i] = M_profile[0] + i += 1 + + # Filling the last values + i = -1 + while np.isnan(res[i]): + res[i] = M_profile[-1] + i += -1 + + if stg.r_bottom.size != 0: + res[np.where(range_cells > r_bottom)] = np.nan + + # print(f"res.shape = {res.shape}") + # print(f"res = {res}") + + loc_point_lin_interp0 = range_cells[np.where((range_cells > sample_depth[0]) & (range_cells < sample_depth[-1]))] + # print(f"range_cells : {range_cells}") + # print(f"loc_point_lin_interp0 : {loc_point_lin_interp0}") + + loc_point_lin_interp = loc_point_lin_interp0[np.where(loc_point_lin_interp0 > range_cells[0])] + # print(f"loc_point_lin_interp : {loc_point_lin_interp}") + + # fig, ax = plt.subplots(nrows=1, ncols=1) + # ax.plot(loc_point_lin_interp, res[:len(loc_point_lin_interp)], c="blue") + # ax.plot(sample_depth, M_profile, color="k", marker="o") + # plt.show() + + return (loc_point_lin_interp, res) # ------------- Computing zeta ------------- # - def zeta(self, ind1, ind2): - # zeta = alpha_s / ((1/r) * (M_profile_SSC_fine)) - zeta0 = np.array([0.04341525, 0.04832906, 0.0847188, np.nan]) - zeta = zeta0[[ind1, ind2]] + def zeta(self, alpha_s, r, M_profile_fine): + + delta_r = r[1] - r[0] + zeta = alpha_s / (np.sum(M_profile_fine*delta_r)) + + # zeta0 = np.array([0.021, 0.035, 0.057, 0.229]) + # zeta = zeta0[ind] + # zeta0 = np.array([0.04341525, 0.04832906, 0.0847188, np.nan]) + # zeta = zeta0[[ind1, ind2]] + + # for k in range(3): + # for p in range(3): + # if np.isnan(ind_X_min_around_sample[p, k]): + # zeta_list_exp.append(np.nan) + # else: + # ind_X_min = int(ind_X_min_around_sample[p, k]) + # ind_X_max = int(ind_X_max_around_sample[p, k]) + # ind_r_min = int(ind_r_min_around_sample[p, k]) + # ind_r_max = int(ind_r_max_around_sample[p, k]) + # + # R_temp = R_cross_section[ind_r_min:ind_r_max, :, ind_X_min:ind_X_max] + # J_temp = J_cross_section[ind_r_min:ind_r_max, :, ind_X_min:ind_X_max] + # aw_temp = aw_cross_section[ind_r_min:ind_r_max, :, ind_X_min:ind_X_max] + # sv_temp_1 = np.repeat([sv_list_temp[3 * k + p]], np.shape(R_temp)[0], axis=0) + # sv_temp = np.swapaxes(np.swapaxes(np.repeat([sv_temp_1], np.shape(R_temp)[2], axis=0), 1, 0), 2, 1) + # ind_depth = np.where(R_cross_section[:, 0, 0] >= M_list_temp[k][0, p + 1])[0][0] + # # Using concentration profile + # zeta_temp = alpha_s / ((1 / M_list_temp[k][0, p + 1]) * (R_cross_section[0, 0, 0] * M_list_temp[k][1, 0] + + # delta_r * np.sum(M_interpolate_list[k][:ind_depth]))) + # zeta_temp = (1 / (4 * R_temp) * + # np.log(sv_temp / J_temp) - aw_temp) / ((1 / M_list_temp[k][0, p + 1]) * + # (R_cross_section[0, 0, 0] * M_list_temp[k][ + # 1, 0] + + # delta_r * np.sum( + # M_interpolate_list[k][:ind_depth]))) + # zeta_list_exp.append(np.mean(np.mean(zeta_temp, axis=0), axis=1)) + return zeta # ------------- Computing VBI ------------- # @@ -122,6 +318,14 @@ class AcousticInversionMethodHighConcentration(): # print('self.water_attenuation[ind_i]', self.water_attenuation[ind_i]) # print('self.x_exp[0.3-1 MHz]', self.x_exp['0.3-1 MHz'].values[0]) # print("start computing VBI") + # print("================================") + # print(f"zeta_freq2 : {zeta_freq2}") + # print(f"j_cross_section_freq1 : {j_cross_section_freq1.shape}") + # print(f"r2D : {r2D.shape}") + # print(f"water_attenuation_freq1 : {water_attenuation_freq1}") + # print(f"freq1 : {freq1}") + # print(f"X : {X}") + # print("================================") logVBI = ((zeta_freq2 * np.log(j_cross_section_freq1 * np.exp(4 * r2D * water_attenuation_freq1) / @@ -131,6 +335,14 @@ class AcousticInversionMethodHighConcentration(): (freq2 ** X))) / (zeta_freq2 - zeta_freq1)) + # logVBI = (( np.full((stg.r.shape[1], stg.t.shape[1]), zeta_freq2) * + # np.log(j_cross_section_freq1 * np.exp(4 * r2D * np.full((stg.r.shape[1], stg.t.shape[1]), water_attenuation_freq1)) / + # (freq1 ** X)) - + # np.full((stg.r.shape[1], stg.t.shape[1]), zeta_freq1) * + # np.log(j_cross_section_freq2 * np.exp(4 * r2D * np.full((stg.r.shape[1], stg.t.shape[1]), water_attenuation_freq2)) / + # (freq2 ** X))) / + # (zeta_freq2 - zeta_freq1)) + print("compute VBI finished") return np.exp(logVBI) diff --git a/Model/granulo_loader.py b/Model/granulo_loader.py index a0d8d24..c67b4ca 100644 --- a/Model/granulo_loader.py +++ b/Model/granulo_loader.py @@ -1,6 +1,7 @@ +import matplotlib.pyplot as plt import numpy as np import pandas as pd - +from Model.GrainSizeTools import demodul_granulo, mix_gaussian_model class GranuloLoader: @@ -15,11 +16,11 @@ class GranuloLoader: self._y = np.array(self._data.iloc[:, 1]) # distance from left bank (m) self._z = np.array(self._data.iloc[:, 2]) # depth (m) - self._r_grain = np.array(self._data.columns.values)[5:] # grain radius (um) + self._r_grain = 1e-6*np.array(self._data.columns.values)[5:].astype(float) / 2 # grain radius (m) self._Ctot = np.array(self._data.iloc[:, 3]) # Total concentration (g/L) self._D50 = np.array(self._data.iloc[:, 4]) # median diameter (um) - self._frac_vol = np.array(self._data.iloc[:, 5:]) # Volume fraction (%) + self._frac_vol = np.array(self._data.iloc[:, 5:])/100 # Volume fraction (%) self._frac_vol_cumul = np.cumsum(self._frac_vol, axis=1) # Cumulated volume fraction (%) @@ -38,11 +39,125 @@ class GranuloLoader: # self._Ctot_fine_per_cent = 100 * self._Ctot_fine / (self._Ctot_fine + self._Ctot_sand) # self._Ctot_sand_per_cent = 100 * self._Ctot_sand / (self._Ctot_fine + self._Ctot_sand) - # print(self._time) + # ============================================================================================================== + # ============================================================================================================== + +# print(self._r_grain.shape) +# N_sample = 2 +# +# fig, ax = plt.subplots(1, 2) +# ax[0].plot(self._r_grain, self._frac_vol[N_sample, :], color="k", marker='.') +# ax[0].set_xscale('log') +# ax[0].set_xlabel('Radius ($\mu m$)') +# ax[0].set_ylabel('Class size volume fraction') +# +# ax[1].plot([self._r_grain[i+1]-self._r_grain[i] for i in range(self._r_grain.shape[0]-1)], list(range(self._r_grain.shape[0]-1)), color="k", marker="x") +# ax[1].set_xlabel('Ecart inter-class') +# ax[1].set_ylabel('n° échantillon') +# +# plt.show() +# +# print(f"self._r_grain : {self._r_grain}") +# print(f"self._frac_vol_cumul[N_sample, :] : {self._frac_vol_cumul[N_sample, :]}") +# +# min_demodul = 1e-6 +# max_demodul = 500e-6 +# sample_demodul = demodul_granulo(self._r_grain, +# self._frac_vol_cumul[N_sample, :], +# min_demodul, max_demodul) +# +# print(f"sample_demodul : {sample_demodul.demodul_data_list}") +# +# N_modes = 3 +# sample_demodul.print_mode_data(N_modes) +# sample_demodul.plot_interpolation() +# sample_demodul.plot_modes(N_modes) +# +# print(f"mu_list : {sample_demodul.demodul_data_list[3 - 1].mu_list}") +# print(f"sigma_list : {sample_demodul.demodul_data_list[3 - 1].sigma_list}") +# print(f"w_list : {sample_demodul.demodul_data_list[3 - 1].w_list}") +# +# resampled_log_array = np.log(np.logspace(-10, -2, 3000)) +# proba_vol_demodul = mix_gaussian_model(resampled_log_array, +# sample_demodul.demodul_data_list[2].mu_list, +# sample_demodul.demodul_data_list[2].sigma_list, +# sample_demodul.demodul_data_list[2].w_list) +# +# +# +# proba_vol_demodul = proba_vol_demodul / np.sum(proba_vol_demodul) +# ss = np.sum(proba_vol_demodul / np.exp(resampled_log_array) ** 3) +# proba_num = proba_vol_demodul / np.exp(resampled_log_array) ** 3 / ss +# +# print(f"proba_num : {proba_num}") +# freq = 5e6 +# +# a2f2pdf = 0 +# a3pdf = 0 +# for i in range(len(resampled_log_array)): +# a = np.exp(resampled_log_array)[i] +# a2f2pdf += a**2 * form_factor_function_MoateThorne2012(a, freq)**2 * proba_num[i] +# a3pdf += a**3 * proba_num[i] +# +# print(f"a2f2pdf = {a2f2pdf}") +# print(f"a3pdf = {a3pdf}") +# +# ks = (a2f2pdf / a3pdf) +# +# print(f"ks = {ks}") +# +# +# def form_factor_function_MoateThorne2012(a_s, freq, C=1500): +# """This function computes the form factor based on the equation of +# Moate and Thorne (2012)""" +# # computing the wave number +# k = 2 * np.pi * freq / C +# x = k * a_s +# f = (x ** 2 * (1 - 0.25 * np.exp(-((x - 1.5) / 0.35) ** 2)) * ( +# 1 + 0.6 * np.exp(-((x - 2.9) / 1.15) ** 2))) / ( +# 42 + 28 * x ** 2) +# return f + # if __name__ == "__main__": +# # GranuloLoader("/home/bmoudjed/Documents/3 SSC acoustic meas project/Graphical interface project/Data/Granulo_data/" +# # "fine_sample_file.ods") # 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/" # "sand_sample_file.ods") + + + + + + + + + +# # form_factor = form_factor_function_MoateThorne2012(GranuloLoader.) +# +# def ks(a_s, rho_s, freq, pdf): +# # --- Calcul de la fonction de form --- +# form_factor = form_factor_function_MoateThorne2012(a_s, freq) +# print(f"form_factor shape = {len(form_factor)}") +# # print(f"form_factor = {form_factor}") +# +# # --- Gaussian mixture --- +# # sample_demodul = demodul_granulo(a_s.astype(float), pdf, 0.17e-6, 200e-6) +# # sample_demodul.plot_interpolation() +# +# # ss = np.sum(pdf / a_s ** 3) +# # proba_num = (pdf / a_s ** 3) / ss +# +# # --- Compute k_s by dividing two integrals --- +# a2f2pdf = 0 +# a3pdf = 0 +# for i in range(len(pdf)): +# a2f2pdf += a_s[i] ** 2 * form_factor[i] * proba_num[i] +# a3pdf += a_s[i] ** 3 * proba_num[i] +# +# ks = np.sqrt(a2f2pdf / (rho_s * a3pdf)) +# +# # ks = np.array([0.04452077, 0.11415143, 0.35533713, 2.47960051]) +# # ks = ks0[ind] +# return ks \ No newline at end of file diff --git a/View/acoustic_inversion_tab.py b/View/acoustic_inversion_tab.py index 42aa969..3305ac5 100644 --- a/View/acoustic_inversion_tab.py +++ b/View/acoustic_inversion_tab.py @@ -93,7 +93,7 @@ class AcousticInversionTab(QWidget): self.gridLayout_groupbox_acoustic_inversion_options.addWidget(self.label_sample_choice, 1, 0, 1, 1) self.label_sample_choice.setText("Calibration samples : ") - self.combobox_calibration_samples = CheckableComboBox() + self.combobox_calibration_samples = QComboBox() self.gridLayout_groupbox_acoustic_inversion_options.addWidget(self.combobox_calibration_samples, 1, 1, 1, 1) self.combobox_calibration_samples.currentIndexChanged.connect(self.sample_choice) @@ -119,7 +119,7 @@ class AcousticInversionTab(QWidget): self.gridLayout_groupbox_acoustic_inversion_settings_parameter.addWidget( self.label_frequencies_pairs_to_compute_VBI, 1, 0, 1, 1) - self.combobox_frequencies_VBI = CheckableComboBox() + self.combobox_frequencies_VBI = QComboBox() self.gridLayout_groupbox_acoustic_inversion_settings_parameter.addWidget( self.combobox_frequencies_VBI, 1, 1, 1, 1) self.combobox_frequencies_VBI.currentIndexChanged.connect(self.frequencies_pair_choice_to_compute_VBI) @@ -129,7 +129,7 @@ class AcousticInversionTab(QWidget): self.gridLayout_groupbox_acoustic_inversion_settings_parameter.addWidget( self.label_frequency_to_compute_SSC, 2, 0, 1, 1) - self.combobox_frequency_SSC = CheckableComboBox() + self.combobox_frequency_SSC = QComboBox() self.gridLayout_groupbox_acoustic_inversion_settings_parameter.addWidget( self.combobox_frequency_SSC, 2, 1, 1, 1) self.combobox_frequency_SSC.currentIndexChanged.connect(self.frequency_choice_to_compute_SSC) @@ -220,166 +220,687 @@ class AcousticInversionTab(QWidget): # --- add items in combobox of samples to calibrate acoustic inversion method --- samples_vertical_line = np.split(stg.samples, np.where(np.diff(stg.sample_time) != 0)[0]+1) + self.combobox_calibration_samples.addItem(" ") for s in samples_vertical_line: self.combobox_calibration_samples.addItem(" - ".join([i for i in s])) - for i in range(len(samples_vertical_line)): - self.combobox_calibration_samples.setItemChecked(i, False) + # for i in range(len(samples_vertical_line)): + # self.combobox_calibration_samples.setItemChecked(i, False) # --- add items in combobox of frequencies for VBI computation --- + self.combobox_frequencies_VBI.addItem(" ") for k in combinations(stg.freq_text, 2): self.combobox_frequencies_VBI.addItem(k[0] + " - " + k[1]) # print(k) - for i in range(len(list(combinations(stg.freq_text, 2)))): - self.combobox_frequencies_VBI.setItemChecked(i, False) + # for i in range(len(list(combinations(stg.freq_text, 2)))): + # self.combobox_frequencies_VBI.setItemChecked(i, False) + + print(f"stg.fine_sediment_columns length : {len(stg.fine_sediment_columns)}") + print(f"stg.fine_sediment_columns : {stg.fine_sediment_columns}") + print(f"stg.frac_vol_fine.shape : {stg.frac_vol_fine.shape}") + print(f"stg.frac_vol_fine : {stg.frac_vol_fine}") def sample_choice(self): sample_position = [] - for i in range(self.combobox_calibration_samples.count()): - if self.combobox_calibration_samples.itemChecked(i): - sample_position.append(i) - elif (i in sample_position) and (not self.combobox_calibration_samples.itemChecked(i)): - sample_position.remove(i) + # for i in range(self.combobox_calibration_samples.count()): + # if self.combobox_calibration_samples.itemChecked(i): + sample_position.append(self.combobox_calibration_samples.currentIndex()) + # elif (i in sample_position) and (not self.combobox_calibration_samples.itemChecked(i)): + # sample_position.remove(i) + print(f"sample position : {sample_position}") def frequencies_pair_choice_to_compute_VBI(self): freq_combination = list(combinations(stg.freq, 2)) frequencies_position = [] - for i in range(self.combobox_frequencies_VBI.count()): - if self.combobox_frequencies_VBI.itemChecked(i): - frequencies_position.append(i) - elif (i in frequencies_position) and (not self.combobox_frequencies_VBI.itemChecked(i)): - frequencies_position.remove(i) + # for i in range(self.combobox_frequencies_VBI.count()): + # if self.combobox_frequencies_VBI.itemChecked(i): + # frequencies_position.append(i) + # elif (i in frequencies_position) and (not self.combobox_frequencies_VBI.itemChecked(i)): + # frequencies_position.remove(i) + frequencies_position.append(self.combobox_frequencies_VBI.currentIndex()) + print(f"frequencies_position : {frequencies_position}") if len(frequencies_position) != 0: # print(frequencies_position) # print(freq_combination[frequencies_position[0]][0], freq_combination[frequencies_position[0]][1]) stg.frequencies_to_compute_VBI = ( - np.array([[int(np.where(stg.freq == freq_combination[frequencies_position[0]][0])[0][0]), - freq_combination[frequencies_position[0]][0]], - [int(np.where(stg.freq == freq_combination[frequencies_position[0]][1])[0][0]), - freq_combination[frequencies_position[0]][1]]])) + np.array([[int(np.where(stg.freq == freq_combination[frequencies_position[0]-1][0])[0][0]), + freq_combination[frequencies_position[0]-1][0]], + [int(np.where(stg.freq == freq_combination[frequencies_position[0]-1][1])[0][0]), + freq_combination[frequencies_position[0]-1][1]]])) + print(f"stg.frequencies_to_compute_VBI : {stg.frequencies_to_compute_VBI}") # --- add items in combobox of frequency for SSC computation --- - for k in range(stg.frequencies_to_compute_VBI.shape[0]): - self.combobox_frequency_SSC.addItem(str(1e-6*stg.frequencies_to_compute_VBI[k, 1]) + " MHz") - for i in range(stg.frequencies_to_compute_VBI.shape[0]): - self.combobox_frequency_SSC.setItemChecked(i, False) + # for k in range(stg.frequencies_to_compute_VBI.shape[0]+1): + # self.combobox_frequency_SSC.removeItem(k) + self.combobox_frequency_SSC.clear() + for k in range(stg.frequencies_to_compute_VBI.shape[0]+1): + if k == 0: + self.combobox_frequency_SSC.addItem(" ") + else: + self.combobox_frequency_SSC.addItem(str(1e-6*stg.frequencies_to_compute_VBI[k-1, 1]) + " MHz") + # for i in range(stg.frequencies_to_compute_VBI.shape[0]): + # self.combobox_frequency_SSC.setItemChecked(i, False) print("frequencies to compute VBI", stg.frequencies_to_compute_VBI) - def frequency_choice_to_compute_SSC(self, index): - # print(self.combobox_frequency_SSC.currentText()) - # print(self.combobox_frequency_SSC.currentIndex()) + def frequency_choice_to_compute_SSC(self): + print(self.combobox_frequency_SSC.currentText()) + print(self.combobox_frequency_SSC.currentIndex()) # print(self.combobox_frequency_SSC.itemChecked(index)) - if self.combobox_frequency_SSC.itemChecked(index): + # if self.combobox_frequency_SSC.itemChecked(index): # # itemChecked(index)): # currentIndex() == 0) or (self.combobox_frequency_SSC.currentIndex() == 1): # print(self.combobox_frequency_SSC.currentText()) # print(stg.freq_text) # print(np.where(np.array(stg.freq_text) == self.combobox_frequency_SSC.currentText())) - stg.frequency_to_compute_SSC \ - = np.array([int(np.where(np.array(stg.freq_text) == self.combobox_frequency_SSC.currentText())[0][0]), - stg.freq[int(np.where(np.array(stg.freq_text) == self.combobox_frequency_SSC.currentText())[0][0])]]) - print("stg.frequency_to_compute_SSC ", stg.frequency_to_compute_SSC) + # stg.frequency_to_compute_SSC \ + # = np.array([int(np.where(np.array(stg.freq_text) == self.combobox_frequency_SSC.currentText())[0][0]), + # stg.freq[int(np.where(np.array(stg.freq_text) == self.combobox_frequency_SSC.currentText())[0][0])]]) + # + # print("stg.frequency_to_compute_SSC ", stg.frequency_to_compute_SSC) def temperature_value(self): stg.temperature = self.spinbox_temperature.value() # print(stg.temperature) - def compute_acoustic_inversion_method_high_concentration(self): - stg.water_attenuation = self.inv_hc.water_attenuation(stg.frequencies_to_compute_VBI[0, 1], stg.frequencies_to_compute_VBI[1, 1], stg.temperature) - # print("water attenuation ", stg.water_attenuation) - + def compute_sound_velocity(self): stg.water_velocity = self.inv_hc.water_velocity(stg.temperature) - # print("water velocity ", stg.water_velocity) + print("water velocity ", stg.water_velocity) - stg.kt_corrected = self.inv_hc.kt_corrected(stg.r, stg.water_velocity, - stg.gain_rx[[int(stg.frequencies_to_compute_VBI[0, 0]), int(stg.frequencies_to_compute_VBI[1, 0])]], - stg.gain_tx[[int(stg.frequencies_to_compute_VBI[0, 0]), int(stg.frequencies_to_compute_VBI[1, 0])]], - stg.kt[[int(stg.frequencies_to_compute_VBI[0, 0]), int(stg.frequencies_to_compute_VBI[1, 0])]]) - # print("kt ", stg.kt_corrected) - # print("kt shape ", stg.kt_corrected.shape) + # def compute_kt(self, freq_ind): + # + # # stg.kt_corrected = self.inv_hc.kt_corrected(stg.r[int(stg.frequency_to_compute_SSC[0]), :], stg.water_velocity, + # # stg.gain_rx[[int(stg.frequencies_to_compute_VBI[0, 0]), int(stg.frequencies_to_compute_VBI[1, 0])]], + # # stg.gain_tx[[int(stg.frequencies_to_compute_VBI[0, 0]), int(stg.frequencies_to_compute_VBI[1, 0])]], + # # stg.kt[[int(stg.frequencies_to_compute_VBI[0, 0]), int(stg.frequencies_to_compute_VBI[1, 0])]]) + # + # if stg.ABS_name == "Aquascat 1000R": + # kt_corrected = self.inv_hc.kt_corrected(stg.r[0, :], stg.water_velocity, + # stg.gain_rx[freq_ind], stg.gain_tx[freq_ind], stg.kt[0, freq_ind]) + # kt_corrected_2D = np.repeat(kt_corrected, stg.r.shape[1], axis=0) + # print("kt 2D ", kt_corrected_2D) + # print("kt 2D shape ", kt_corrected_2D.shape) + # kt_corrected_3D = np.zeros((kt_corrected_2D.shape[1], kt_corrected_2D.shape[0], stg.t.shape[1])) + # for k in range(kt_corrected_2D.shape[1]): + # kt_corrected_3D[k, :, :] = np.repeat(kt_corrected_2D, stg.t.shape[1], axis=1)[:, + # k * stg.t.shape[1]:(k + 1) * stg.t.shape[1]] + # print("kt 3D ", kt_corrected_3D) + # print("kt 3D shape ", kt_corrected_3D.shape) + # elif stg.ABS_name == "UB-SediFlow": + # stg.kt_corrected = np.array([[0.003, 0.006]]) + # print("kt ", stg.kt_corrected) + # print("kt shape ", stg.kt_corrected.shape) + # stg.kt_corrected_2D = np.array(np.repeat(stg.kt_corrected, stg.r.shape[1], axis=0)) + # print("kt 2D shape ", stg.kt_corrected_2D.shape) + # stg.kt_corrected_3D = np.zeros((stg.kt_corrected_2D.shape[1], stg.kt_corrected_2D.shape[0], stg.t.shape[1])) + # for k in range(stg.kt_corrected_2D.shape[1]): + # stg.kt_corrected_3D[k, :, :] = np.repeat(stg.kt_corrected_2D, stg.t.shape[1], axis=1)[:, + # k * stg.t.shape[1]:(k + 1) * stg.t.shape[1]] + # print("kt corrected 3D zeros shape ", stg.kt_corrected_3D.shape) + # + # print("kt ", stg.kt_corrected) + # # print("kt shape ", stg.kt_corrected.shape) + # + # # stg.kt_corrected_2D = np.repeat(stg.kt_corrected, stg.r.shape[1], axis=0) + # + # # stg.kt_corrected_3D = np.zeros((stg.kt_corrected_2D.shape[1], stg.kt_corrected_2D.shape[0], stg.t.shape[1])) + # # print("stg.t.shape ", stg.t.shape) + # # print("kt corrected 3D zeros shape ", stg.kt_corrected_3D.shape) + # # for k in range(stg.kt_corrected_2D.shape[1]): + # # stg.kt_corrected_3D[k, :, :] = np.repeat(stg.kt_corrected_2D, stg.t.shape[1], axis=1)[:, k*stg.t.shape[1]:(k+1)*stg.t.shape[1]] + # + # # print("kt 2D", np.repeat(stg.kt_corrected[:, :, np.newaxis], stg.t.shape[0], axis=2)) - stg.kt_corrected_2D = np.repeat(stg.kt_corrected, stg.r.shape[0], axis=0) - # print("kt 2D ", stg.kt_corrected_2D) - # print("kt 2D shape ", stg.kt_corrected_2D.shape) - stg.kt_corrected_3D = np.zeros((stg.kt_corrected_2D.shape[1], stg.kt_corrected_2D.shape[0], stg.t.shape[0])) - # print("stg.t.shape ", stg.t.shape) - # print("kt corrected 3D zeros shape ", stg.kt_corrected_3D.shape) - for k in range(stg.kt_corrected_2D.shape[1]): - stg.kt_corrected_3D[k, :, :] = np.repeat(stg.kt_corrected_2D, stg.t.shape[0], axis=1)[:, k*stg.t.shape[0]:(k+1)*stg.t.shape[0]] - print("kt 3D ", stg.kt_corrected_3D) - print("kt 3D shape ", stg.kt_corrected_3D.shape) - # print("kt 2D", np.repeat(stg.kt_corrected[:, :, np.newaxis], stg.t.shape[0], axis=2)) + def compute_J(self, freq_ind, kt): + if stg.ABS_name == "Aquascat 1000R": + print(f"stg.BS_stream_bed.shape : {stg.BS_stream_bed.shape}") + if stg.BS_stream_bed.size == 0: - if stg.BS_data_section.size == 0: - stg.J_cross_section = ( - self.inv_hc.j_cross_section( - stg.BS_data[:, [int(stg.frequencies_to_compute_VBI[0, 0]), - int(stg.frequencies_to_compute_VBI[1, 0])], :], - stg.r_2D, stg.kt_corrected_3D)) - elif (stg.BS_data_section_averaged.size == 0) and (stg.BS_data_section_SNR_filter.size == 0): - stg.J_cross_section = ( - self.inv_hc.j_cross_section( - stg.BS_data_section[:, [int(stg.frequencies_to_compute_VBI[0, 0]), - int(stg.frequencies_to_compute_VBI[1, 0])], :], - stg.r_2D, stg.kt_corrected_3D)) - elif (stg.BS_data_section_averaged.size != 0) and (stg.BS_data_section_SNR_filter.size == 0): - stg.J_cross_section = ( - self.inv_hc.j_cross_section( - stg.BS_data_section_averaged[:, [int(stg.frequencies_to_compute_VBI[0, 0]), - int(stg.frequencies_to_compute_VBI[1, 0])], :], - stg.r_2D, stg.kt_corrected_3D)) - else: - print("stg.r_2D.shape ", stg.r_2D.shape) - print("stg.kt_corrected_3D.shape ", stg.kt_corrected_3D.shape) - stg.J_cross_section = ( - self.inv_hc.j_cross_section( - stg.BS_data_section_SNR_filter[:, [int(stg.frequencies_to_compute_VBI[0, 0]), - int(stg.frequencies_to_compute_VBI[1, 0])], :], - stg.r_2D, stg.kt_corrected_3D)) + if stg.BS_cross_section_averaged.size != 0: + print("1/ stg.BS_cross_section_averaged") + stg.J_cross_section = ( + self.inv_hc.j_cross_section(stg.BS_cross_section_averaged[freq_ind, :, :], + stg.r_2D[0, :, :stg.BS_cross_section_averaged.shape[2]], + kt[freq_ind, :, :])) + # stg.J_cross_section = ( + # self.inv_hc.j_cross_section( + # stg.BS_cross_section_averaged[ + # [int(stg.frequencies_to_compute_VBI[0, 0]), int(stg.frequencies_to_compute_VBI[1, 0])], + # :, :], + # stg.r_2D[int(stg.frequency_to_compute_SSC[0]), :, :stg.BS_cross_section_averaged.shape[2]], + # stg.kt_corrected_3D)) - print("J ", stg.J_cross_section) - print("J sahpe ", stg.J_cross_section.shape) + elif stg.BS_cross_section_SNR_filter.size != 0: + print("2/ stg.BS_cross_section_SNR_filter") + stg.J_cross_section = ( + self.inv_hc.j_cross_section( + stg.BS_cross_section_SNR_filter[freq_ind, :, :], + stg.r_2D[freq_ind, :, :stg.BS_cross_section_SNR_filter.shape[2]], + kt[freq_ind, :, :])) - stg.X_exponent = self.inv_hc.X_exponent(self.combobox_frequencies_VBI.currentIndex()) - print("X ", stg.X_exponent) + # stg.J_cross_section = ( + # self.inv_hc.j_cross_section( + # stg.BS_cross_section_SNR_filter[ + # [int(stg.frequencies_to_compute_VBI[0, 0]), int(stg.frequencies_to_compute_VBI[1, 0])], + # :, :], + # stg.r_2D[int(stg.frequency_to_compute_SSC[0]), :, + # :stg.BS_cross_section_SNR_filter.shape[2]], + # stg.kt_corrected_3D)) - stg.zeta = self.inv_hc.zeta(int(stg.frequencies_to_compute_VBI[0, 0]), int(stg.frequencies_to_compute_VBI[1, 0])) - print("zeta ", stg.zeta) + else: + print("3/ stg.BS_cross_section") + stg.J_cross_section = ( + self.inv_hc.j_cross_section( + stg.BS_cross_section[freq_ind, :, :], + stg.r_2D[freq_ind, :, :stg.BS_cross_section.shape[2]], + kt[freq_ind, :, :])) - stg.ks = self.inv_hc.ks(int(stg.frequency_to_compute_SSC[0])) - print("ks ", stg.ks) + # stg.J_cross_section = ( + # self.inv_hc.j_cross_section( + # stg.BS_cross_section[ + # [int(stg.frequencies_to_compute_VBI[0, 0]), int(stg.frequencies_to_compute_VBI[1, 0])], + # :, :], + # stg.r_2D[int(stg.frequency_to_compute_SSC[0]), :, + # :stg.BS_cross_section.shape[2]], + # stg.kt_corrected_3D)) - stg.VBI_cross_section = self.inv_hc.VBI_cross_section(stg.frequencies_to_compute_VBI[0, 1], stg.frequencies_to_compute_VBI[1, 1], - stg.zeta[0], # zeta is already limited to the frequencies pairs so that we just need to select indices 0 and 1 - stg.zeta[1], - stg.J_cross_section[0, :, :], - stg.J_cross_section[1, :, :], - stg.r_2D, - stg.water_attenuation[0], - stg.water_attenuation[1], - stg.X_exponent) + else: + + if stg.BS_stream_bed_averaged.size != 0: + print("1/ stg.BS_stream_bed_averaged") + stg.J_cross_section = ( + self.inv_hc.j_cross_section( + stg.BS_stream_bed_averaged[freq_ind, :, :], + stg.r_2D[freq_ind, :, :stg.BS_stream_bed_averaged.shape[2]], + kt[freq_ind, :, :])) + + # stg.J_stream_bed = ( + # self.inv_hc.j_cross_section( + # stg.BS_stream_bed_averaged[ + # [int(stg.frequencies_to_compute_VBI[0, 0]), int(stg.frequencies_to_compute_VBI[1, 0])], + # :, :], + # stg.r_2D[int(stg.frequency_to_compute_SSC[0]), :, :stg.BS_stream_bed_averaged.shape[2]], + # stg.kt_corrected_3D)) + + elif stg.BS_stream_bed_SNR_filter.size != 0: + print("2/ stg.BS_stream_bed_SNR_filter") + stg.J_cross_section = ( + self.inv_hc.j_cross_section( + stg.BS_stream_bed_SNR_filter[freq_ind, :, :], + stg.r_2D[freq_ind, :, :stg.BS_stream_bed_SNR_filter.shape[2]], + kt[freq_ind, :, :])) + + # stg.J_stream_bed = ( + # self.inv_hc.j_cross_section( + # stg.BS_stream_bed_SNR_filter[ + # [int(stg.frequencies_to_compute_VBI[0, 0]), int(stg.frequencies_to_compute_VBI[1, 0])], + # :, :], + # stg.r_2D[int(stg.frequency_to_compute_SSC[0]), :, :stg.BS_stream_bed_SNR_filter.shape[2]], + # stg.kt_corrected_3D)) + + else: + print("3/ stg.BS_stream_bed") + stg.J_cross_section = ( + self.inv_hc.j_cross_section( + stg.BS_stream_bed[freq_ind, :, :], + stg.r_2D[freq_ind, :, :stg.BS_stream_bed.shape[2]], + kt[freq_ind, :, :])) + + # stg.J_stream_bed = ( + # self.inv_hc.j_cross_section( + # stg.BS_stream_bed[ + # [int(stg.frequencies_to_compute_VBI[0, 0]), int(stg.frequencies_to_compute_VBI[1, 0])], + # :, :], + # stg.r_2D[int(stg.frequency_to_compute_SSC[0]), :, + # :stg.BS_stream_bed.shape[2]], + # stg.kt_corrected_3D)) + + # elif (stg.BS_data_section_averaged.size == 0) and (stg.BS_data_section_SNR_filter.size == 0): + # stg.J_cross_section = ( + # self.inv_hc.j_cross_section( + # stg.BS_data_section[[int(stg.frequencies_to_compute_VBI[0, 0]), + # int(stg.frequencies_to_compute_VBI[1, 0])], :, :], + # stg.r_2D[int(stg.frequency_to_compute_SSC[0]), :, :stg.BS_data_section.shape[2]], stg.kt_corrected_3D)) + # elif (stg.BS_data_section_averaged.size != 0) and (stg.BS_data_section_SNR_filter.size == 0): + # print("Je suis dans ce J") + # stg.J_cross_section = ( + # self.inv_hc.j_cross_section( + # stg.BS_data_section_averaged[[0, 2], :, :], + # stg.r_2D[[0, 2], :, :stg.BS_data_section_averaged.shape[2]], + # stg.kt_corrected_3D[[0, 1], :, :])) + # stg.J_cross_section = ( + # self.inv_hc.j_cross_section( + # stg.BS_data_section_averaged[[int(stg.frequencies_to_compute_VBI[0, 0]), + # int(stg.frequencies_to_compute_VBI[1, 0])], :, :], + # stg.r_2D[int(stg.frequency_to_compute_SSC[0]), :, :stg.BS_data_section_averaged.shape[2]], stg.kt_corrected_3D)) + # else: + # print("stg.r_2D.shape ", stg.r_2D.shape) + # print("stg.kt_corrected_3D.shape ", stg.kt_corrected_3D.shape) + # stg.J_cross_section = ( + # self.inv_hc.j_cross_section( + # stg.BS_data_section_SNR_filter[[int(stg.frequencies_to_compute_VBI[0, 0]), + # int(stg.frequencies_to_compute_VBI[1, 0])], :, :], + # stg.r_2D[int(stg.frequency_to_compute_SSC[0]), :, :stg.BS_data_section_SNR_filter.shape[2]], stg.kt_corrected_3D)) + # + # print("J ", stg.J_cross_section) + # print("J sahpe ", stg.J_cross_section.shape) + + elif stg.ABS_name == "UB-SediFlow": + + if stg.BS_stream_bed.size == 0: + + if stg.BS_cross_section_averaged.size != 0: + print("1/ stg.BS_cross_section_averaged") + stg.J_cross_section = ( + self.inv_hc.j_cross_section(stg.BS_cross_section_averaged[freq_ind, :, :], + stg.r_2D[0, :, :stg.BS_cross_section_averaged.shape[2]], + kt[freq_ind, :, :])) + + elif stg.BS_cross_section_SNR_filter.size != 0: + print("2/ stg.BS_cross_section_SNR_filter") + stg.J_cross_section = ( + self.inv_hc.j_cross_section( + stg.BS_cross_section_SNR_filter[freq_ind, :, :], + stg.r_2D[freq_ind, :, :stg.BS_cross_section_SNR_filter.shape[2]], + kt[freq_ind, :, :])) + + else: + print("3/ stg.BS_cross_section") + stg.J_cross_section = ( + self.inv_hc.j_cross_section( + stg.BS_cross_section[freq_ind, :, :], + stg.r_2D[freq_ind, :, :stg.BS_cross_section.shape[2]], + kt[freq_ind, :, :])) + + else: + + if stg.BS_stream_bed_averaged.size != 0: + print("1/ stg.BS_stream_bed_averaged") + stg.J_cross_section = ( + self.inv_hc.j_cross_section( + stg.BS_stream_bed_averaged[freq_ind, :, :], + stg.r_2D[freq_ind, :, :stg.BS_stream_bed_averaged.shape[2]], + kt[freq_ind, :, :])) + + elif stg.BS_stream_bed_SNR_filter.size != 0: + print("2/ stg.BS_stream_bed_SNR_filter") + stg.J_cross_section = ( + self.inv_hc.j_cross_section( + stg.BS_stream_bed_SNR_filter[freq_ind, :, :], + stg.r_2D[freq_ind, :, :stg.BS_stream_bed_SNR_filter.shape[2]], + kt[freq_ind, :, :])) + + else: + print("3/ stg.BS_stream_bed") + stg.J_cross_section = ( + self.inv_hc.j_cross_section( + stg.BS_stream_bed[freq_ind, :, :], + stg.r_2D[freq_ind, :, :stg.BS_stream_bed.shape[2]], + kt[freq_ind, :, :])) + + + # if stg.BS_data_section.size == 0: + # stg.J_cross_section = ( + # self.inv_hc.j_cross_section( + # stg.BS_data[ + # [int(stg.frequencies_to_compute_VBI[0, 0]), int(stg.frequencies_to_compute_VBI[1, 0])], + # :, :], stg.r_2D[int(stg.frequency_to_compute_SSC[0]), :, :stg.BS_data.shape[2]], + # stg.kt_corrected_3D)) + # elif (stg.BS_data_section_averaged.size == 0) and (stg.BS_data_section_SNR_filter.size == 0): + # stg.J_cross_section = ( + # self.inv_hc.j_cross_section( + # stg.BS_data_section[[int(stg.frequencies_to_compute_VBI[0, 0]), + # int(stg.frequencies_to_compute_VBI[1, 0])], :, :], + # np.array([stg.r_2D[0, :, :stg.BS_data_section.shape[2]]]), + # stg.kt_corrected_3D)) + # elif (stg.BS_data_section_averaged.size != 0) and (stg.BS_data_section_SNR_filter.size == 0): + # print("Je suis dans ce J") + # stg.J_cross_section = ( + # self.inv_hc.j_cross_section( + # stg.BS_data_section_averaged[[0, 2], :, :], + # stg.r_2D[[0, 2], :, :stg.BS_data_section_averaged.shape[2]], + # stg.kt_corrected_3D[[0, 1], :, :])) + # # stg.J_cross_section = ( + # # self.inv_hc.j_cross_section( + # # stg.BS_data_section_averaged[[int(stg.frequencies_to_compute_VBI[0, 0]), + # # int(stg.frequencies_to_compute_VBI[1, 0])], :, :], + # # stg.r_2D[int(stg.frequency_to_compute_SSC[0]), :, :stg.BS_data_section_averaged.shape[2]], stg.kt_corrected_3D)) + # else: + # print("stg.r_2D.shape ", stg.r_2D.shape) + # print("stg.kt_corrected_3D.shape ", stg.kt_corrected_3D.shape) + # stg.J_cross_section = ( + # self.inv_hc.j_cross_section( + # stg.BS_data_section_SNR_filter[[int(stg.frequencies_to_compute_VBI[0, 0]), + # int(stg.frequencies_to_compute_VBI[1, 0])], :, :], + # stg.r_2D[int(stg.frequency_to_compute_SSC[0]), :, :stg.BS_data_section_SNR_filter.shape[2]], + # stg.kt_corrected_3D)) + # + # print("J ", stg.J_cross_section) + # print("J sahpe ", stg.J_cross_section.shape) + + return stg.J_cross_section + + # def compute_alpha_s(self, freq, freq_ind): + # stg.sv = self.compute_sv(freq) + # j_cross_section = self.compute_J(freq_ind) + # stg.alpha_s = np.log(stg.sv / j_cross_section) / (4 * stg.sample_depth[10]) - stg.water_attenuation + # return + + # def compute_zeta(self): + # stg.zeta_freq1 = self.inv_hc.zeta(0, 2) + # stg.zeta_freq2 = self.inv_hc.zeta() + # # stg.zeta = self.inv_hc.zeta(int(stg.frequencies_to_compute_VBI[0, 0]), int(stg.frequencies_to_compute_VBI[1, 0])) + # print("zeta ", stg.zeta) + + # def compute_X(self): + # stg.X_exponent = self.inv_hc.X_exponent(2) #self.inv_hc.X_exponent(self.combobox_frequencies_VBI.currentIndex()) + # print("X ", stg.X_exponent) + + # def compute_VBI(self): + # stg.VBI_cross_section = self.inv_hc.VBI_cross_section(3e5, + # 1e6, + # stg.zeta[0], + # # zeta is already limited to the frequencies pairs so that we just need to select indices 0 and 1 + # stg.zeta[1], + # stg.J_cross_section[0, :, :], + # stg.J_cross_section[1, :, :], + # stg.r_2D[0, :, :stg.t.shape[1]], + # stg.water_attenuation[0], + # stg.water_attenuation[1], + # stg.X_exponent) + + # stg.VBI_cross_section = self.inv_hc.VBI_cross_section(stg.frequencies_to_compute_VBI[0, 1], stg.frequencies_to_compute_VBI[1, 1], + # stg.zeta[0], # zeta is already limited to the frequencies pairs so that we just need to select indices 0 and 1 + # stg.zeta[1], + # stg.J_cross_section[0, :, :], + # stg.J_cross_section[1, :, :], + # stg.r_2D[int(stg.frequency_to_compute_SSC[0]), :, :stg.t.shape[1]], + # stg.water_attenuation[0], + # stg.water_attenuation[1], + # stg.X_exponent) # print("VBI shape ", stg.VBI_cross_section.shape) # print(int(self.combobox_frequency_SSC.currentIndex())) # print(stg.zeta[int(self.combobox_frequency_SSC.currentIndex())]) - stg.SSC_fine = self.inv_hc.SSC_fine(stg.zeta[int(self.combobox_frequency_SSC.currentIndex())], - stg.r_2D, - stg.VBI_cross_section, - stg.frequency_to_compute_SSC[1], - stg.X_exponent, - stg.J_cross_section[self.combobox_frequency_SSC.currentIndex(), :, :]) - # print("SSC fine shape ", stg.SSC_fine.shape) + def compute_acoustic_inversion_method_high_concentration(self): - stg.SSC_sand = self.inv_hc.SSC_sand(stg.VBI_cross_section, - stg.frequency_to_compute_SSC[1], - stg.X_exponent, - stg.ks) + self.temperature_value() - # print("SSC sand shape ", stg.SSC_sand.shape) + if stg.ABS_name == "Aquascat 1000R": + + freq1 = 0 # 0 = 300kHz et 1 = 500kHz + freq2 = 2 # 2 = 1MHz + + stg.water_velocity = self.inv_hc.water_velocity(stg.temperature) + + alpha_w_freq1, alpha_w_freq2 = (self.inv_hc.water_attenuation(freq=stg.freq[freq1], T=stg.temperature), + self.inv_hc.water_attenuation(freq=stg.freq[freq2], T=stg.temperature)) + print(f"alpha_w_freq1 = {alpha_w_freq1}, alpha_w_freq2 = {alpha_w_freq2}") + + self.compute_sound_velocity() + + # print("ks value : ", self.inv_hc.ks(a_s=stg.sand_sediment_columns[5:], rho_s=2500, freq=stg.freq[0], pdf=stg.frac_vol_sand[2, :])) + + ks_freq1, ks_freq2 = ( + self.inv_hc.ks(freq=stg.freq[freq1], pdf=stg.frac_vol_sand_cumul[2, :]), + self.inv_hc.ks(freq=stg.freq[freq2], pdf=stg.frac_vol_sand_cumul[2, :])) + # 6 = V11 , 10 = V12 sur le notebook d'Adrien + # ks_freq1, ks_freq2 = self.inv_hc.ks()[0], self.inv_hc.ks()[2] + print(f"ks_freq1 = {ks_freq1}, ks_freq2 = {ks_freq2}") + + # f = self.inv_hc.form_factor_function_MoateThorne2012(np.log(np.logspace(-10, -2, 3000)), stg.freq[2]) + # fig, ax = plt.subplots(nrows=1, ncols=1) + # ax.plot(list(range(len(np.log(np.logspace(-10, -2, 3000))))), f, color="k", ls="solid") + # plt.show() + + sv_freq1, sv_freq2 = ( + self.inv_hc.sv(ks=ks_freq1, M_sand=stg.Ctot_sand[2]), self.inv_hc.sv(ks=ks_freq2, M_sand=stg.Ctot_sand[2])) + print(f"sv_freq1 = {sv_freq1}, sv_freq2 = {sv_freq2}") + + X_exponent = self.inv_hc.X_exponent(freq1=stg.freq[freq1], freq2=stg.freq[freq2], sv_freq1=sv_freq1, sv_freq2=sv_freq2) + # X_exponent = self.inv_hc.X_exponent(ind=2) + print(f"X_exponent = {X_exponent}") + + stg.kt = np.array([[0.04], [0.01838], [0.03267], [0.01376]]) + kt = np.array([]) + for i, v in enumerate(stg.kt): + kt = np.append(kt, self.inv_hc.kt_corrected(r=stg.r[i, :], + water_velocity=stg.water_velocity, + RxGain=stg.gain_rx[i], + TxGain=stg.gain_tx[i], + kt_ref=stg.kt[i])) + kt = np.reshape(kt, (len(kt), 1)) + print(f"kt = {kt}") + kt2D = np.repeat(np.array(kt), stg.r.shape[1], axis=1) + print(f"kt2D.shape = {kt2D.shape}") + # print(f"kt2D = {kt2D}") + kt3D = np.repeat(kt2D[:, :, np.newaxis], stg.t.shape[1], axis=2) + print(f"kt3D.shape = {kt3D.shape}") + # print(f"kt3D = {kt3D}") + + # kt_freq1, kt_freq2 = ( + # self.inv_hc.kt_corrected(r=stg.r[0, :], water_velocity=stg.water_velocity, RxGain=stg.gain_rx[0], + # TxGain=stg.gain_tx[0], kt_ref=stg.kt[0]), + # self.inv_hc.kt_corrected(r=stg.r[2, :], water_velocity=stg.water_velocity, RxGain=stg.gain_rx[2], + # TxGain=stg.gain_tx[2], kt_ref=stg.kt[2])) + # print(f"kt_freq1 = {kt_freq1}, kt_freq2 = {kt_freq2}") + # kt2D_freq1 = np.full((stg.r.shape[1], stg.t.shape[1]), kt_freq1) + # kt2D_freq2 = np.full((stg.r.shape[1], stg.t.shape[1]), kt_freq2) + # print(f"kt2D_freq1 = {kt2D_freq1.shape}") + # print(f"kt2D_freq2 = {kt2D_freq2.shape}") + # kt3D = np.repeat(kt2D[:, :, np.newaxis], stg.t.shape[1], axis=2) + # print(f"kt_3D = {kt3D.shape}") + # print(f"kt_3D = {kt3D}") + + # kt_freq1, kt_freq2 = ( + # np.full((stg.r.shape[1], stg.t.shape[1]), stg.kt[0]), + # np.full((stg.r.shape[1], stg.t.shape[1]), stg.kt[2])) + # print(f"kt_freq1 = {kt_freq1.shape}, kt_freq2 = {kt_freq2.shape}") + + J_freq1, J_freq2 = self.compute_J(freq_ind=freq1, kt=kt3D), self.compute_J(freq_ind=freq2, kt=kt3D) + print(f"J_freq1 = {J_freq1.shape}, J_freq2 = {J_freq2.shape}") + + # fig, ax = plt.subplots(nrows=1, ncols=1) + # pcm = ax.pcolormesh(stg.t[0, :], stg.r[0, :], J_freq2, cmap='rainbow', shading='gouraud') + # fig.colorbar(pcm, ax=ax, shrink=1, location='right') + # plt.show() + + # print("***************************") + # print(sv_freq1.shape) + # print(J_freq1.shape) + # print(np.repeat(stg.r[0, :][:, np.newaxis], stg.t.shape[1], axis=1).shape) + # print(np.full((stg.r.shape[1], stg.t.shape[1]), alpha_w_freq1).shape) + # print("***************************") + + r_sample_ind = [] + t_sample_ind = [] + for i in range(len(stg.sample_depth[:3])): + # print(-stg.sample_depth[i]) + # print(-stg.sample_time[i]) + r_sample_ind.append(np.where(np.abs(stg.r[0, :] - (-stg.sample_depth[i])) == + np.nanmin(np.abs(stg.r[0, :] - (-stg.sample_depth[i]))))) + # print(np.abs(stg.t[0, :] - (-stg.sample_time[i]))) + t_sample_ind.append(np.where(np.abs(stg.t[0, :] - (stg.sample_time[i])) == + np.nanmin(np.abs(stg.t[0, :] - (stg.sample_time[i]))))) + print(f"r_sample_ind = {r_sample_ind}") + print(f"t_sample_ind = {t_sample_ind}") + # print("J_freq1[r_sample_ind[1][0], t_sample_ind[1][0]] ", J_freq1[r_sample_ind[1][0], t_sample_ind[1][0]]) + + # ind_r_min = int(ind_X_min_around_sample[p, k]) + # print(f"ind_r_min = {ind_r_min}") + # ind_r_max = int(ind_X_max_around_sample[p, k]) + # print(f"ind_r_max = {ind_r_max}") + # ind_t_min = int(ind_r_min_around_sample[p, k]) + # print(f"ind_t_min = {ind_t_min}") + # ind_t_max = int(ind_r_max_around_sample[p, k]) + # print(f"ind_t_max = {ind_t_max}") + + + print(f"J_freq1[r_sample_ind[-1], t_sample_ind[-1]] {J_freq1[r_sample_ind[-1], t_sample_ind[-1]]}") + print(f"J_freq2[r_sample_ind[-1], t_sample_ind[-1]] {J_freq2[r_sample_ind[-1], t_sample_ind[-1]]}") + + print(f"stg.r[0, r_sample_ind[-1]] {stg.r[2, r_sample_ind[-1]]}") + + + alpha_s_freq1, alpha_s_freq2 = ( + self.inv_hc.alpha_s(sv=sv_freq1, j_cross_section=J_freq1[r_sample_ind[-1], t_sample_ind[-1]], depth=stg.r[freq1, r_sample_ind[-1]], alpha_w=alpha_w_freq1), + self.inv_hc.alpha_s(sv=sv_freq2, j_cross_section=J_freq2[r_sample_ind[-1], t_sample_ind[-1]], depth=stg.r[freq2, r_sample_ind[-1]], alpha_w=alpha_w_freq2)) + print(f"alpha_s_freq1 = {alpha_s_freq1}, alpha_s_freq2 = {alpha_s_freq2}") + + range_lin_interp, M_profile_fine = self.inv_hc.M_profile_SCC_fine_interpolated(sample_depth=-stg.sample_depth[:3], + M_profile=stg.Ctot_fine[:3], + range_cells=stg.r[0, :], + r_bottom=stg.r_bottom[t_sample_ind[0][0]]) + M_profile_fine = M_profile_fine[:len(range_lin_interp)] + + print(f"range_lin_interp : {range_lin_interp}") + print(f"M_profile_fine : {M_profile_fine}") + + # print("----------------------") + # print(f"r_sample_ind[-1][0] = {r_sample_ind[-1][0][0]}") + # print(f"M_profile_fine[:r_sample_ind[-1][0]] = ", M_profile_fine[:r_sample_ind[-1][0][0]]) + # print(f"stg.r[0, r_sample_ind[-1][0]] = ", stg.r[0, r_sample_ind[-1][0][0]]) + # print("----------------------") + + zeta_freq1, zeta_freq2 = (self.inv_hc.zeta(alpha_s=alpha_s_freq1, r=stg.r[freq1, :], M_profile_fine=M_profile_fine), + self.inv_hc.zeta(alpha_s=alpha_s_freq2, r=stg.r[freq2, :], M_profile_fine=M_profile_fine)) + # zeta_freq1, zeta_freq2 = ( + # self.inv_hc.zeta(r=stg.r[0, :r_sample_ind[-1][0][0]], M_profile_fine=M_profile_fine[:r_sample_ind[-1][0][0]]), + # self.inv_hc.zeta(r=stg.r[0, :r_sample_ind[-1][0][0]], M_profile_fine=M_profile_fine[:r_sample_ind[-1][0][0]])) + print(f"zeta_freq1 = {zeta_freq1}, zeta_freq2 = {zeta_freq2}") + # plt.figure() + # plt.plot(stg.r[0, :], Mprofile, 'b.', -stg.sample_depth[:3], stg.Ctot_fine[:3], 'ko') + # plt.show() + + stg.VBI_cross_section = self.inv_hc.VBI_cross_section(stg.freq[freq1], stg.freq[freq2], + zeta_freq1, zeta_freq2, + J_freq1, J_freq2, + stg.r_2D[freq1, :, :stg.t.shape[1]], + alpha_w_freq1, alpha_w_freq2, + X_exponent) + + stg.SSC_fine = self.inv_hc.SSC_fine(zeta_freq2, stg.r_2D[freq2, :, :stg.t.shape[1]], + stg.VBI_cross_section, stg.freq[freq2], X_exponent, J_freq2) + + stg.SSC_sand = self.inv_hc.SSC_sand(stg.VBI_cross_section, stg.freq[freq2], X_exponent, ks_freq2) + + elif stg.ABS_name == "UB-SediFlow": + + freq1 = 0 # 0 = 500kHz + freq2 = 1 # 1 = 1MHz + + stg.water_velocity = self.inv_hc.water_velocity(stg.temperature) + + alpha_w_freq1, alpha_w_freq2 = (self.inv_hc.water_attenuation(freq=stg.freq[freq1], T=stg.temperature), + self.inv_hc.water_attenuation(freq=stg.freq[freq2], T=stg.temperature)) + print(f"alpha_w_freq1 = {alpha_w_freq1}, alpha_w_freq2 = {alpha_w_freq2}") + + self.compute_sound_velocity() + + # ks_freq1, ks_freq2 = 0.11373812635175432, 0.35705575378038723 + ks_freq1, ks_freq2 = (self.inv_hc.form_factor_function_MoateThorne2012(a=100e-6, freq=stg.freq[freq1])/np.sqrt(2500*100e-6), + self.inv_hc.form_factor_function_MoateThorne2012(a=100e-6, freq=stg.freq[freq2])/np.sqrt(2500*100e-6)) + print(f"ks_freq1 = {ks_freq1}, ks_freq2 = {ks_freq2}") + + sv_freq1, sv_freq2 = ( + self.inv_hc.sv(ks=ks_freq1, M_sand=stg.Ctot_sand[-1]), + self.inv_hc.sv(ks=ks_freq2, M_sand=stg.Ctot_sand[-1])) + print(f"sv_freq1 = {sv_freq1}, sv_freq2 = {sv_freq2}") + + X_exponent = self.inv_hc.X_exponent(freq1=stg.freq[freq1], freq2=stg.freq[freq2], sv_freq1=sv_freq1, + sv_freq2=sv_freq2) + # X_exponent = self.inv_hc.X_exponent(ind=2) + print(f"X_exponent = {X_exponent}") + + stg.kt = np.array([[1.38e-3], [6.02e-4]]) # Values of kt for 500kHz and 1MHz + kt = np.array([]) + for i, v in enumerate(stg.kt): + kt = np.append(kt, self.inv_hc.kt_corrected(r=stg.r[i, :], + water_velocity=stg.water_velocity, + RxGain=0, + TxGain=0, + kt_ref=stg.kt[i])) + kt = np.reshape(kt, (len(kt), 1)) + print(f"kt = {kt}, kt.shape = {kt.shape}") + kt2D = np.repeat(np.array(kt), stg.r.shape[1], axis=1) + print(f"kt2D.shape = {kt2D.shape}") + # print(f"kt2D = {kt2D}") + kt3D = np.repeat(kt2D[:, :, np.newaxis], stg.t.shape[1], axis=2) + print(f"kt3D.shape = {kt3D.shape}") + + J_freq1, J_freq2 = self.compute_J(freq_ind=freq1, kt=kt3D), self.compute_J(freq_ind=freq2, kt=kt3D) + print(f"J_freq1 = {J_freq1.shape}, J_freq2 = {J_freq2.shape}") + + r_sample_ind = [] + t_sample_ind = [] + for i in range(len(stg.sample_depth[:])): + # print(-stg.sample_depth[i]) + # print(-stg.sample_time[i]) + r_sample_ind.append(np.where(np.abs(stg.r[freq1, :] - (-stg.sample_depth[i])) == + np.nanmin(np.abs(stg.r[freq1, :] - (-stg.sample_depth[i]))))) + # print(np.abs(stg.t[0, :] - (-stg.sample_time[i]))) + t_sample_ind.append(np.where(np.abs(stg.t[freq1, :] - (stg.sample_time[i])) == + np.nanmin(np.abs(stg.t[freq1, :] - (stg.sample_time[i]))))) + print(f"r_sample_ind = {r_sample_ind}") + print(f"t_sample_ind = {t_sample_ind}") + + alpha_s_freq1, alpha_s_freq2 = ( + self.inv_hc.alpha_s(sv=sv_freq1, j_cross_section=J_freq1[r_sample_ind[-1], t_sample_ind[-1]], + depth=stg.r[freq1, r_sample_ind[-1]], alpha_w=alpha_w_freq1), + self.inv_hc.alpha_s(sv=sv_freq2, j_cross_section=J_freq2[r_sample_ind[-1], t_sample_ind[-1]], + depth=stg.r[freq2, r_sample_ind[-1]], alpha_w=alpha_w_freq2)) + print(f"alpha_s_freq1 = {alpha_s_freq1}, alpha_s_freq2 = {alpha_s_freq2}") + + # range_lin_interp, M_profile_fine = self.inv_hc.M_profile_SCC_fine_interpolated( + # sample_depth=-stg.sample_depth[:], + # M_profile=stg.Ctot_sand[:], + # range_cells=stg.r[0, :], + # r_bottom=np.array([])) # stg.r_bottom[t_sample_ind[0][0]] is empty + # M_profile_fine = M_profile_fine[:len(range_lin_interp)] + range_lin_interp, M_profile_fine = 3.2, 1 + + print(f"range_lin_interp : {range_lin_interp}") + print(f"M_profile_fine : {M_profile_fine}") + + zeta_freq1, zeta_freq2 = (alpha_s_freq1 / (M_profile_fine*(3.19127224 - 0.52102404)), + alpha_s_freq2 / (M_profile_fine*(3.19127224 - 0.52102404))) + # zeta_freq1, zeta_freq2 = ( + # self.inv_hc.zeta(alpha_s=alpha_s_freq1, r=stg.r[freq1, :], M_profile_fine=M_profile_fine), + # self.inv_hc.zeta(alpha_s=alpha_s_freq2, r=stg.r[freq2, :], M_profile_fine=M_profile_fine)) + + print(f"zeta_freq1 = {zeta_freq1}, zeta_freq2 = {zeta_freq2}") + + # zeta_freq1, zeta_freq2 = 1e-10*self.inv_hc.zeta(ind=1), 1e1*self.inv_hc.zeta(ind=2) + # print(f"zeta_freq1 = {zeta_freq1}, zeta_freq2 = {zeta_freq2}") + + + stg.VBI_cross_section = self.inv_hc.VBI_cross_section(stg.freq[freq1], stg.freq[freq2], + zeta_freq1, zeta_freq2, + J_freq1, J_freq2, + stg.r_2D[freq2, :, :stg.t.shape[1]], + alpha_w_freq1, alpha_w_freq2, + X_exponent) + + stg.SSC_fine = self.inv_hc.SSC_fine(zeta_freq2, stg.r_2D[freq2, :, :stg.t.shape[1]], + stg.VBI_cross_section, stg.freq[freq2], X_exponent, J_freq2) + # stg.SSC_fine = self.inv_hc.SSC_fine(stg.zeta[int(self.combobox_frequency_SSC.currentIndex())], + # stg.r_2D[int(stg.frequency_to_compute_SSC[0]), :, :stg.t.shape[1]], + # stg.VBI_cross_section, + # stg.frequency_to_compute_SSC[1], + # stg.X_exponent, + # stg.J_cross_section[self.combobox_frequency_SSC.currentIndex(), :, :]) + # print("SSC fine shape ", stg.SSC_fine.shape) + + stg.SSC_sand = self.inv_hc.SSC_sand(stg.VBI_cross_section, stg.freq[freq2], X_exponent, ks_freq2) + # stg.SSC_sand = self.inv_hc.SSC_sand(stg.VBI_cross_section, + # stg.frequency_to_compute_SSC[1], + # stg.X_exponent, + # stg.ks) + + # print("SSC sand shape ", stg.SSC_sand.shape) def plot_SSC_2D_fields(self): @@ -394,23 +915,56 @@ class AcousticInversionTab(QWidget): def plot_SSC_fine(self): - val_min = 1e-2 - val_max = 15 - # val_min = np.nanmin(stg.SSC_fine) - # val_max = np.nanmax(stg.SSC_fine) - print('val_min fine = ', np.nanmin(stg.SSC_fine)) - print('val_max fine =', np.nanmax(stg.SSC_fine)) + # val_min = 1e-2 + # val_max = 15 + val_min = np.nanmin(stg.SSC_fine) + val_max = np.nanmax(stg.SSC_fine) + print('val_min fine = ', val_min) + print('val_max fine =', val_max) + # print('val_min fine = ', np.nanmin(stg.VBI_cross_section)) + # print('val_max fine =', np.nanmax(stg.VBI_cross_section)) # if val_min == 0: # val_min = 0.5 # print('val_min update =', val_min) - pcm_SSC_fine = self.axis_SSC_2D_field[0].pcolormesh(stg.t, -stg.r, stg.SSC_fine, - cmap='rainbow', - norm=LogNorm(vmin=val_min, vmax=val_max), - shading='gouraud') + # pcm_SSC_fine = self.axis_SSC_2D_field[0].pcolormesh(stg.t[0, :], + # -stg.r[0, :], + # stg.SSC_fine, + # # np.log(stg.VBI_cross_section/1.1932980954310682e-27), + # # cmap='jet', + # # vmin=val_min, vmax=val_max) + # cmap = 'rainbow', + # norm=LogNorm(vmin=1e-1, vmax=15)) + # # shading='gouraud') + if stg.ABS_name == "Aquascat 1000R": + pcm_SSC_fine = self.axis_SSC_2D_field[0].pcolormesh(stg.t[0, :], + -stg.r[0, :], + stg.SSC_fine, + cmap='rainbow', + norm=LogNorm(vmin=1e-2, vmax=15), + shading='gouraud') + # pcm_SSC_fine = self.axis_SSC_2D_field[0].pcolormesh(stg.t[int(stg.frequency_to_compute_SSC[0]), :], + # -stg.r[int(stg.frequency_to_compute_SSC[0]), :], + # stg.SSC_fine, + # cmap='rainbow', + # norm=LogNorm(vmin=1e-2, vmax=15), + # shading='gouraud') + elif stg.ABS_name == "UB-SediFlow": + pcm_SSC_fine = self.axis_SSC_2D_field[0].pcolormesh(stg.t[0, :], + -stg.r[0, :], + stg.SSC_fine, + cmap='rainbow', + norm=LogNorm(vmin=1e-2, vmax=1), + shading='gouraud') + if stg.r_bottom.size != 0: - self.axis_SSC_2D_field[0].plot(stg.t, -stg.r_bottom, color='black', linewidth=1, linestyle="solid") + self.axis_SSC_2D_field[0].plot(stg.t[0, :], + -stg.r_bottom, + color='black', linewidth=1, linestyle="solid") + # self.axis_SSC_2D_field[0].plot(stg.t[int(stg.frequency_to_compute_SSC[0]), :], + # -stg.r_bottom, + # color='black', linewidth=1, linestyle="solid") self.figure_SSC_2D_field.supxlabel("Time (sec)", fontsize=10) self.figure_SSC_2D_field.supylabel("Depth (m)", fontsize=10) @@ -424,22 +978,57 @@ class AcousticInversionTab(QWidget): val_min = 1e-2 val_max = 2 - # val_min = np.nanmin(stg.SSC_sand) - # val_max = np.nanmax(stg.SSC_sand) - print('val_min sand = ', np.nanmin(stg.SSC_sand)) - print('val_max sand = ', np.nanmax(stg.SSC_sand)) + val_min = np.nanmin(stg.SSC_sand) + val_max = np.nanmax(stg.SSC_sand) + # val_min = np.nanmin(np.log(stg.VBI_cross_section)) + # val_max = np.nanmax(np.log(stg.VBI_cross_section)) + print('val_min sand = ', val_min) + print('val_max sand =', val_max) + # print('val_min sand = ', np.nanmin(stg.VBI_cross_section)) + # print('val_max sand = ', np.nanmax(stg.VBI_cross_section)) # if val_min == 0: # val_min = 0.5 # print('val_min update =', val_min) - pcm_SSC_sand = self.axis_SSC_2D_field[1].pcolormesh(stg.t, -stg.r, stg.SSC_sand, - cmap='rainbow', - # vmin=val_min, vmax=val_max, - norm=LogNorm(vmin=val_min, vmax=val_max), - shading='gouraud') + # print(stg.SSC_sand) + + if stg.ABS_name == "Aquascat 1000R": + pcm_SSC_sand = self.axis_SSC_2D_field[1].pcolormesh(stg.t[0, :], + -stg.r[0, :], + (stg.SSC_sand), + cmap='rainbow', + # vmin=-60, vmax=-45) + # vmin=1e-2, vmax=10) + # vmin=val_min, vmax=val_max, + norm=LogNorm(vmin=1e-2, vmax=2), + shading='gouraud') + + elif stg.ABS_name == "UB-SediFlow": + pcm_SSC_sand = self.axis_SSC_2D_field[1].pcolormesh(stg.t[0, :], + -stg.r[0, :], + (stg.SSC_sand), + cmap='rainbow', + # vmin=-60, vmax=-45) + # vmin=1e-2, vmax=10) + # vmin=val_min, vmax=val_max, + norm=LogNorm(vmin=1e-2, vmax=1), + shading='gouraud') + + # pcm_SSC_sand = self.axis_SSC_2D_field[1].pcolormesh(stg.t[int(stg.frequency_to_compute_SSC[0]), :], + # -stg.r[int(stg.frequency_to_compute_SSC[0]), :], + # stg.SSC_sand, + # cmap='rainbow', + # # vmin=val_min, vmax=val_max, + # norm=LogNorm(vmin=1e-2, vmax=2), + # shading='gouraud') if stg.r_bottom.size: - self.axis_SSC_2D_field[1].plot(stg.t, -stg.r_bottom, color='black', linewidth=1, linestyle="solid") + self.axis_SSC_2D_field[1].plot(stg.t[0, :], + -stg.r_bottom, + color='black', linewidth=1, linestyle="solid") + # self.axis_SSC_2D_field[1].plot(stg.t[int(stg.frequency_to_compute_SSC[0]), :], + # -stg.r_bottom, + # color='black', linewidth=1, linestyle="solid") self.figure_SSC_2D_field.supxlabel("Time (sec)", fontsize=10) self.figure_SSC_2D_field.supylabel("Depth (m)", fontsize=10) @@ -453,32 +1042,62 @@ class AcousticInversionTab(QWidget): sample_depth_position = [] for i in range(stg.sample_depth.shape[0]): sample_depth_position.append( - np.where(np.abs(stg.r + stg.sample_depth[i]) == np.min(np.abs(stg.r + stg.sample_depth[i])))[0][0]) + np.where(np.abs(stg.r[0, :] + stg.sample_depth[i]) == + np.min(np.abs(stg.r[0, :] + stg.sample_depth[i])))[0][0]) sample_time_position = [] for j in range(stg.sample_time.shape[0]): sample_time_position.append( - np.where(np.abs(stg.t - stg.sample_time[j]) == np.min(np.abs(stg.t - stg.sample_time[j])))[0][0]) + np.where(np.abs(stg.t[0, :] - stg.sample_time[j]) == + np.min(np.abs(stg.t[0, :] - stg.sample_time[j])))[0][0]) if self.canvas_SSC_sample_vs_inversion == None: print("Ctot fine : ", stg.Ctot_fine) print("SCC fine : ", stg.SSC_fine[sample_depth_position, sample_time_position]) - print("Ctot fine : ", stg.Ctot_sand) - print("SCC fine : ", stg.SSC_sand[sample_depth_position, sample_time_position]) + print("Ctot sand : ", stg.Ctot_sand) + print("SCC sand : ", stg.SSC_sand[sample_depth_position, sample_time_position]) self.figure_SSC_sample_vs_inversion, self.axis_SSC_sample_vs_inversion = plt.subplots(nrows=1, ncols=1, layout="constrained") self.canvas_SSC_sample_vs_inversion = FigureCanvas(self.figure_SSC_sample_vs_inversion) self.verticalLayout_groupbox_SSC_sample_vs_inversion.addWidget(self.canvas_SSC_sample_vs_inversion) - self.axis_SSC_sample_vs_inversion.plot(stg.Ctot_fine, stg.SSC_fine[sample_depth_position, sample_time_position], ls=" ", marker='v', color='black', label='Fine SSC') - self.axis_SSC_sample_vs_inversion.plot(stg.Ctot_sand, stg.SSC_sand[sample_depth_position, sample_time_position], ls=" ", marker='x', color='black', label='Sand SSC') - self.axis_SSC_sample_vs_inversion.set_xscale('log') - self.axis_SSC_sample_vs_inversion.set_yscale('log') - # self.axis_SSC_sample_vs_inversion.plot([0, 100], [0, 100], color='black', lw=1) - self.axis_SSC_sample_vs_inversion.set_xlabel('Measured SSC (g/l)', weight='bold') - self.axis_SSC_sample_vs_inversion.set_ylabel('Inverse SSC (g/l)', weight='bold') - self.axis_SSC_sample_vs_inversion.legend() + if stg.ABS_name == "Aquascat 1000R": + self.axis_SSC_sample_vs_inversion.plot(stg.Ctot_fine, stg.SSC_fine[sample_depth_position, sample_time_position], ls=" ", marker='v', color='black', label='Fine SSC') + self.axis_SSC_sample_vs_inversion.plot(stg.Ctot_sand, + stg.SSC_sand[sample_depth_position, sample_time_position], + ls=" ", marker='x', color='black', label='Sand SSC') + + self.axis_SSC_sample_vs_inversion.set_xscale('log') + self.axis_SSC_sample_vs_inversion.set_yscale('log') + self.axis_SSC_sample_vs_inversion.plot([0, 10], [0, 10], color='black', lw=1) + self.axis_SSC_sample_vs_inversion.set_xlabel('Measured SSC (g/l)', weight='bold') + self.axis_SSC_sample_vs_inversion.set_ylabel('Inverse SSC (g/l)', weight='bold') + self.axis_SSC_sample_vs_inversion.legend() + + elif stg.ABS_name == "UB-SediFlow": + # self.axis_SSC_sample_vs_inversion.plot(stg.Ctot_sand, + # stg.SSC_sand[sample_depth_position, sample_time_position], + # ls=" ", marker='x', color='black', label='Sand SSC') + + # fig, ax = plt.subplots(nrows=2, ncols=10) + # for j in range(2): + # for i in range(10): + # ax[j, i].plot(stg.r[1, :145], stg.SSC_sand[:145, 750+j*10+i]) + # plt.show() + # fig, ax = plt.subplots(nrows=1, ncols=1) + # for i in range(20): + # ax.plot(stg.r[1, :145], stg.SSC_sand[:145, 750+i]) + # plt.show() + + # self.axis_SSC_sample_vs_inversion.plot(stg.r[1, :145], stg.SSC_sand[:145, 40], color='blue', ls='solid') + self.axis_SSC_sample_vs_inversion.plot(stg.r[1, 2:145], stg.SSC_sand[2:145, 762], + color='red', ls='solid', marker='o', mec='k', mfc='k', ms=3) + self.axis_SSC_sample_vs_inversion.set_xlabel('Depth (m)', weight='bold') + self.axis_SSC_sample_vs_inversion.set_ylabel('Sand concentration (g/l)', weight='bold') + + + diff --git a/View/note_tab.py b/View/note_tab.py index 71de6da..08b7bf2 100644 --- a/View/note_tab.py +++ b/View/note_tab.py @@ -200,12 +200,12 @@ class NoteTab(QWidget): f" Fine sediments data file: {stg.fine_sediment_path}/{stg.fine_sediment_filename} \n" f" Sand sediments data file: {stg.sand_sediment_path}/{stg.sand_sediment_filename} \n" "\n\n" - "------------------------------------------------------------------------- \n\n\n" + "------------------------------------------------------------------------- \n\n\n") - "Acoustic Inversion parameters: \n" - f" frequencies to compute VBI: {stg.freq_text[int(stg.frequencies_to_compute_VBI[0, 0])]}, " - f"{stg.freq_text[int(stg.frequencies_to_compute_VBI[1, 0])]} \n" - f" frequency to compute SSC: {stg.freq_text[int(stg.frequency_to_compute_SSC[0])]}") + # "Acoustic Inversion parameters: \n" + # f" frequencies to compute VBI: {stg.freq_text[int(stg.frequencies_to_compute_VBI[0, 0])]}, " + # f"{stg.freq_text[int(stg.frequencies_to_compute_VBI[1, 0])]} \n" + # f" frequency to compute SSC: {stg.freq_text[int(stg.frequency_to_compute_SSC[0])]}") # if __name__ == "__main__": diff --git a/View/sample_data_tab.py b/View/sample_data_tab.py index c97b4a2..263c0db 100644 --- a/View/sample_data_tab.py +++ b/View/sample_data_tab.py @@ -773,7 +773,7 @@ class SampleDataTab(QWidget): def plot_sample_position_on_transect(self): - if stg.BS_data.size == 0: + if stg.BS_cross_section.size == 0: msgBox = QMessageBox() msgBox.setWindowTitle("Plot transect Error") msgBox.setIcon(QMessageBox.Warning) @@ -968,7 +968,8 @@ class SampleDataTab(QWidget): else: indices = [] for i in position_list: - indices.append(np.where(stg.t == stg.sample_time[i])[0][0]) + # indices.append(np.where(stg.t == stg.sample_time[i])[0][0]) + indices.append(np.where(np.abs(stg.t - stg.sample_time[i]) == np.nanmin(np.abs(stg.t - stg.sample_time[i])))[0][0]) self.axis_total_concentration.scatter(stg.Ctot_fine[position_list], -stg.sample_depth[position_list] / stg.r_bottom[indices], # np.max(stg.sample_depth[position_list]), @@ -1002,7 +1003,10 @@ class SampleDataTab(QWidget): else: indices = [] for i in position_list: - indices.append(np.where(stg.t == stg.sample_time[i])[0][0]) + # indices.append(np.where(stg.t == stg.sample_time[i])[0][0]) + indices.append(np.where( + np.abs(stg.t - stg.sample_time[i]) == np.nanmin(np.abs(stg.t - stg.sample_time[i])))[0][ + 0]) self.axis_total_concentration.scatter(stg.Ctot_fine_per_cent[position_list], -stg.sample_depth[position_list] / stg.r_bottom[indices], # np.max(stg.sample_depth[position_list]), @@ -1050,7 +1054,10 @@ class SampleDataTab(QWidget): else: indices = [] for i in position_list: - indices.append(np.where(stg.t == stg.sample_time[i])[0][0]) + # indices.append(np.where(stg.t == stg.sample_time[i])[0][0]) + indices.append(np.where( + np.abs(stg.t - stg.sample_time[i]) == np.nanmin(np.abs(stg.t - stg.sample_time[i])))[0][ + 0]) self.axis_total_concentration.scatter(stg.Ctot_sand[position_list], -stg.sample_depth[position_list] / stg.r_bottom[indices], # np.max(stg.sample_depth[position_list]), @@ -1084,7 +1091,10 @@ class SampleDataTab(QWidget): else: indices = [] for i in position_list: - indices.append(np.where(stg.t == stg.sample_time[i])[0][0]) + # indices.append(np.where(stg.t == stg.sample_time[i])[0][0]) + indices.append(np.where( + np.abs(stg.t - stg.sample_time[i]) == np.nanmin(np.abs(stg.t - stg.sample_time[i])))[0][ + 0]) self.axis_total_concentration.scatter(stg.Ctot_sand_per_cent[position_list], -stg.sample_depth[position_list] / stg.r_bottom[indices], # np.max(stg.sample_depth[position_list]), @@ -1138,7 +1148,10 @@ class SampleDataTab(QWidget): else: indices = [] for i in position_list: - indices.append(np.where(stg.t == stg.sample_time[i])[0][0]) + # indices.append(np.where(stg.t == stg.sample_time[i])[0][0]) + indices.append(np.where( + np.abs(stg.t - stg.sample_time[i]) == np.nanmin(np.abs(stg.t - stg.sample_time[i])))[0][ + 0]) self.axis_total_concentration.scatter(stg.Ctot_fine[position_list], -stg.sample_depth[position_list] / stg.r_bottom[indices], #np.max(stg.sample_depth[position_list]), s=100, facecolor=color_list, edgecolor="None", alpha=0.5) @@ -1180,7 +1193,10 @@ class SampleDataTab(QWidget): else: indices = [] for i in position_list: - indices.append(np.where(stg.t == stg.sample_time[i])[0][0]) + # indices.append(np.where(stg.t == stg.sample_time[i])[0][0]) + indices.append(np.where( + np.abs(stg.t - stg.sample_time[i]) == np.nanmin(np.abs(stg.t - stg.sample_time[i])))[0][ + 0]) self.axis_total_concentration.scatter(stg.Ctot_fine_per_cent[position_list], -stg.sample_depth[position_list] / stg.r_bottom[indices], #np.max(stg.sample_depth[position_list]), s=100, facecolor=color_list, edgecolor="None", alpha=0.5) diff --git a/View/signal_processing_tab.py b/View/signal_processing_tab.py index da596af..7b0a48c 100644 --- a/View/signal_processing_tab.py +++ b/View/signal_processing_tab.py @@ -633,9 +633,9 @@ class SignalProcessingTab(QWidget): msgBox.setStandardButtons(QMessageBox.Ok) msgBox.exec() else: - if stg.r_bottom.size == 0: - stg.BS_stream_bed = deepcopy(stg.BS_cross_section) - elif stg.r_bottom.size != 0: + # if stg.r_bottom.size == 0: + # stg.BS_stream_bed = deepcopy(stg.BS_cross_section) + if stg.r_bottom.size != 0: stg.BS_stream_bed = deepcopy(stg.BS_cross_section) for f, _ in enumerate(stg.freq): for k, _ in enumerate(stg.r_bottom): @@ -757,50 +757,52 @@ class SignalProcessingTab(QWidget): else: filter_convolve = np.ones(self.spinbox_average_horizontal.value()) - # print(filter_convolve) + print(filter_convolve) if stg.BS_stream_bed_SNR_filter.size != 0: # stg.BS_stream_bed_averaged = np.zeros((stg.r.shape[0], stg.freq.shape[0], stg.t.shape[0])) # stg.BS_stream_bed_averaged = np.zeros((stg.freq.shape[0], stg.r.shape[1], stg.t.shape[1])) stg.BS_stream_bed_averaged = deepcopy(stg.BS_stream_bed_SNR_filter) - # print(f"stg.BS_cross_section_averaged : {stg.BS_cross_section_averaged.shape}") + print(f"1/ stg.BS_stream_bed_averaged : {stg.BS_stream_bed_averaged}") for f, _ in enumerate(stg.freq): for i in range(stg.r.shape[1]): stg.BS_stream_bed_averaged[f, i, :] \ - = convolve1d(stg.BS_stream_bed_SNR_filter[f, i, :], weights=filter_convolve)# / filter_convolve.shape[0] + = convolve1d(stg.BS_stream_bed_SNR_filter[f, i, :], weights=filter_convolve) / len(filter_convolve) # stg.BS_stream_bed_averaged[i, f, :] \ # = convolve1d(stg.BS_cross_section[i, f, :], weights=filter_convolve) / filter_convolve.shape[0] elif stg.BS_cross_section_SNR_filter.size != 0: stg.BS_cross_section_averaged = deepcopy(stg.BS_cross_section_SNR_filter) - print(f"stg.BS_cross_section_averaged : {stg.BS_cross_section_averaged.shape}") + print(f"2/ stg.BS_cross_section_averaged : {stg.BS_cross_section_averaged}") for f, _ in enumerate(stg.freq): for i in range(stg.r.shape[1]): stg.BS_stream_bed_averaged[f, i, :] \ = convolve1d(stg.BS_cross_section_SNR_filter[f, i, :], - weights=filter_convolve) # / filter_convolve.shape[0] + weights=filter_convolve) / len(filter_convolve) elif stg.BS_stream_bed.size != 0: stg.BS_stream_bed_averaged = deepcopy(stg.BS_stream_bed) - print(f"stg.BS_cross_section_averaged : {stg.BS_cross_section_averaged.shape}") + print(f"3/ stg.BS_stream_bed_averaged : {stg.BS_stream_bed_averaged.shape}") + print(stg.BS_stream_bed_averaged[0, 10, :50]) for f, _ in enumerate(stg.freq): for i in range(stg.r.shape[1]): stg.BS_stream_bed_averaged[f, i, :] \ = convolve1d(stg.BS_stream_bed[f, i, :], - weights=filter_convolve) # / filter_convolve.shape[0] + weights=filter_convolve) / len(filter_convolve) + print(stg.BS_stream_bed_averaged[0, 10, :50]) elif stg.BS_cross_section.size != 0: stg.BS_cross_section_averaged = deepcopy(stg.BS_cross_section) - print(f"stg.BS_cross_section_averaged : {stg.BS_cross_section_averaged.shape}") + print(f"4/ stg.BS_cross_section_averaged : {stg.BS_cross_section_averaged.shape}") for f, _ in enumerate(stg.freq): for i in range(stg.r.shape[1]): - stg.BS_stream_bed_averaged[f, i, :] \ + stg.BS_cross_section_averaged[f, i, :] \ = convolve1d(stg.BS_cross_section[f, i, :], - weights=filter_convolve) # / filter_convolve.shape[0] + weights=filter_convolve) / len(filter_convolve) if stg.ABS_name == "Aquascat 1000R": self.label_cells_horizontal.clear() @@ -983,32 +985,42 @@ class SignalProcessingTab(QWidget): self.verticalLayout_groupbox_display_profile_position.addWidget(self.canvas_plot_profile_position_on_transect) # --- Plot transect with profile position --- - val_min = np.nanmin(stg.BS_stream_bed[stg.freq_bottom_detection, :, :]) - val_max = np.nanmax(stg.BS_stream_bed[stg.freq_bottom_detection, :, :]) - if val_min == 0: - val_min = 1e-5 if stg.ABS_name == "Aquascat 1000R": if stg.BS_stream_bed.size != 0: + + val_min = np.nanmin(stg.BS_stream_bed[stg.freq_bottom_detection, :, :]) + val_max = np.nanmax(stg.BS_stream_bed[stg.freq_bottom_detection, :, :]) + if val_min == 0: + val_min = 1e-5 + self.axis_plot_profile_position_on_transect.pcolormesh( stg.t[stg.freq_bottom_detection, :], -stg.r[stg.freq_bottom_detection, :], stg.BS_stream_bed[stg.freq_bottom_detection, :, :], cmap='viridis', norm=LogNorm(vmin=val_min, vmax=val_max)) else: + + val_min = np.nanmin(stg.BS_cross_section[stg.freq_bottom_detection, :, :]) + val_max = np.nanmax(stg.BS_cross_section[stg.freq_bottom_detection, :, :]) + if val_min == 0: + val_min = 1e-5 + self.axis_plot_profile_position_on_transect.pcolormesh( stg.t[stg.freq_bottom_detection, :], -stg.r[stg.freq_bottom_detection, :], stg.BS_cross_section[stg.freq_bottom_detection, :, :], cmap='viridis', norm=LogNorm(vmin=val_min, vmax=val_max)) + elif stg.ABS_name == "UB-SediFlow": + if stg.BS_stream_bed.size != 0: self.axis_plot_profile_position_on_transect.pcolormesh( stg.t[stg.freq_bottom_detection, :], -stg.r[stg.freq_bottom_detection, :], - np.log(stg.BS_cross_section[stg.freq_bottom_detection, :, :]), + np.log(stg.BS_stream_bed[stg.freq_bottom_detection, :, :]), cmap='Blues') else: self.axis_plot_profile_position_on_transect.pcolormesh( stg.t[stg.freq_bottom_detection, :], -stg.r[stg.freq_bottom_detection, :], - np.log(stg.BS_stream_bed[stg.freq_bottom_detection, :, :]), + np.log(stg.BS_cross_section[stg.freq_bottom_detection, :, :]), cmap='Blues') if stg.r_bottom.size != 0: @@ -1116,8 +1128,8 @@ class SignalProcessingTab(QWidget): self.axis_plot_profile_position_on_transect.cla() - val_min = np.nanmin(stg.BS_stream_bed[self.combobox_frequency.currentIndex(), :, :]) - val_max = np.nanmax(stg.BS_stream_bed[self.combobox_frequency.currentIndex(), :, :]) + val_min = np.nanmin(stg.BS_cross_section[self.combobox_frequency.currentIndex(), :, :]) + val_max = np.nanmax(stg.BS_cross_section[self.combobox_frequency.currentIndex(), :, :]) if val_min == 0: val_min = 1e-5 diff --git a/main.py b/main.py index 8026364..e4bec1f 100644 --- a/main.py +++ b/main.py @@ -12,6 +12,9 @@ from View.acoustic_inversion_tab import AcousticInversionTab from View.note_tab import NoteTab from View.user_manual_tab import UserManualTab +import matplotlib.pyplot as plt +plt.close("all") + # Check encoding used # print(sys.getdefaultencoding()) diff --git a/settings.py b/settings.py index e5ba377..bbd5c45 100644 --- a/settings.py +++ b/settings.py @@ -134,6 +134,8 @@ VBI_stream_bed = np.array([[[]]]) frequency_to_compute_SSC = 0 ks = 0 +sv = 0 +alpha_s = 0 SSC_fine = np.array(()) SSC_sand = np.array([])