Acoustic inversion method for high concentration is implemented. It is tested with Aquascat data acquired during Adrien PhD (Confluence Rhône-Isère, 2018.01.07)
brahim
parent
3dd6a4e310
commit
cb54fdd7ac
|
|
@ -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 <nb_modes> 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()
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -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__":
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
||||
|
|
|
|||
3
main.py
3
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())
|
||||
|
||||
|
|
|
|||
|
|
@ -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([])
|
||||
|
|
|
|||
Loading…
Reference in New Issue