Implemented the computation of the VBI on a common grid of 'r' values for the two selected frequencies. Exporting to a .csv file and updating the sand SSC plots still need work.

Model/TableModel.py

@@ -1,5 +1,5 @@
 from PyQt5.QtCore import Qt, QAbstractTableModel
-
+import numpy as np
 
 class TableModel(QAbstractTableModel):
     def __init__(self, data):
@@ -9,20 +9,17 @@ class TableModel(QAbstractTableModel):
     def data(self, index, role):
         if role == Qt.DisplayRole:
             value = self._data.iloc[index.row(), index.column()]
-            # if role == Qt.TextAlignmentRole:
-            #     value = self._data.iloc[index.row(), index.column()]
-            #     if isinstance(value, int) or isinstance(value, float):
-            #         return Qt.AlignVCenter + Qt.AlignRight
 
             if isinstance(value, float):
+                
+                if np.isnan(value):
+                    return "NaN"
+                
                 # Render float to 2 dp
-                if len(str(value).split(".")[1]) <= 3:
+                elif len(str(value).split(".")[1]) <= 3:
                     return "%.2f" % value
                 else:
                     return "%.2e" % value
-            # if isinstance(value, str):
-            #     # Render strings with quotes
-            #     return '"%s"' % value
 
             return value
 
@@ -41,4 +38,4 @@ class TableModel(QAbstractTableModel):
             if orientation == Qt.Horizontal:
                 return str(self._data.columns[section])
             if orientation == Qt.Vertical:
-                return str(self._data.index[section])
+                return str(self._data.index[section])

Model/acoustic_data_loader_UBSediFlow.py

@@ -206,8 +206,8 @@ class AcousticDataLoaderUBSediFlow:
                 self._nb_profiles[k] += 1
 
 
-            # Finally, form the 2D arrays for each frequency and store them in 3D NumPy arrays:
+            # Finally, store the 2D arrays for each frequency in 3D NumPy arrays:
-
+            
             self._BS_raw_data[k,:,:] = bs_list[k]
             self._time[k,:] = time_list[k]
             self._r[k,:] = r_list[k]
@@ -215,10 +215,48 @@ class AcousticDataLoaderUBSediFlow:
         
 
     def reshape_BS_raw_data(self):
-        BS_raw_cross_section = np.reshape(self._BS_raw_data,
+        
-                                          (self._r.shape[1]*self._time.shape[1], len(self._freq)),
+        '''
-                                          order="F")
+        Form and return a table from the raw BS (backscatter) data loaded into AcouSed. That table
+        is the one to be displayed in the "Table of values" section of the "Acoustic data" tab.
+
+        The table is a L*K NumPy array, where 'K' stands for the number of frequencies in the
+        loaded dataset, and L = R*T with 'R' the number of vertical coordinates and 'T' the number
+        of time stamps.
+
+        Finally, the BS array is flattened (i.e., converted to a 1D array) for each frequency so that
+        the index spanning the vertical coordinates runs the fastest i.e.
+
+        Line 0  : BS data at (t, r) = (t_0, r_0)
+        Line 1  : BS data at (t, r) = (t_0, r_1)
+        Line 2  : BS data at (t, r) = (t_0, r_2)
+        ...
+        Line R  : BS data at (t, r) = (t_0, r_R)
+        Line R+1: BS data at (t, r) = (t_1, r_0)
+        Line R+2: BS data at (t, r) = (t_1, r_1)
+        ...
+
+        '''
+
+
+        # Create and initialise the returned 2D array. Note:
+        # R = self._r.shape[1] = self._BS_raw_data.shape[1]
+        # T = self._time.shape[1] = self._BS_raw_data.shape[2]
+        # K = len(self._freq)
+        
+        BS_raw_cross_section = np.zeros( ( self._r.shape[1] * self._time.shape[1], len(self._freq) ) )
+
+        
+        # Fill the k-th column of 'BS_raw_cross_section' with the BS data recorded with the k-th frequency:
+        
+        for k in range( len(self._freq) ):
+
+            BS_raw_cross_section[:, k] = np.reshape( self._BS_raw_data[k,:,:],
+                                                     self._BS_raw_data.shape[1] * self._BS_raw_data.shape[2],
+                                                     order = "F" )
+            
         return BS_raw_cross_section
+    
 
     def reshape_r(self):
         r = np.zeros((self._r.shape[1]*self._time.shape[1], len(self._freq)))

View/acoustic_data_tab.py

@@ -1718,6 +1718,9 @@ class AcousticDataTab(QWidget):
         stg.kt2D.append(np.array([]))
         stg.kt3D.append(np.array([]))
         stg.J_cross_section.append([np.array([]), np.array([])])
+        stg.J_cross_section_interp.append( [ np.array([]), np.array([]) ] )
+        stg.r2D_interp.append(np.array([]))
+        stg.t2D_interp.append(np.array([]))
         stg.VBI_cross_section.append(np.array([]))
         stg.SSC_fine.append(np.array([]))
         stg.SSC_sand.append(np.array([]))

View/acoustic_inversion_tab.py

@@ -25,6 +25,7 @@ import gc
 import time
 import logging
 import numpy as np
+import scipy as sp
 import pandas as pd
 from math import isinf
 
@@ -529,7 +530,75 @@ class AcousticInversionTab(QWidget):
             self.plot_measured_vs_inverted_SSC_sand()
 
     def compute_VBI(self):
+
+        '''
+        Compute the VBI (Volume Backscatter Index) for the selected dataset in the combobox.
+        If the data has been recorded by an Ubertone UB SediFlow, the data for the two
+        calibration frequencies is first interpolated on the same spatial and temporal
+        grids.
+        
+        '''
+
+        # Get the dataset ID (integer relating the dataset selected in the combobox
+        # to the corresponding dataset loaded into AcouSed):
+        
         data_id = self.combobox_acoustic_data_choice.currentIndex()
+        
+
+        # Get the spatial and temporal grids for the two selected frequencies, and create common
+        # grids to be used for the VBI and SSC computations (reminder: the grids have
+        # the same number of points, but the intervals are different). Here, we assume that the
+        # time delay between the channels (i.e., frequencies) is negligible, so that the i-th
+        # time stamp is taken as the average between the i-th time stamp of the first frequency
+        # and the i-th time stamp of the second frequency:
+
+        rGrid_freq1 = stg.depth_2D[data_id][ stg.frequencies_for_calibration[0][1] ]
+        rGrid_freq2 = stg.depth_2D[data_id][ stg.frequencies_for_calibration[1][1] ]
+
+        rVect_freq1 = rGrid_freq1[:,0] #Vector of vertical coordinates (m) for the first frequency
+        rVect_freq2 = rGrid_freq2[:,0] #Vector of vertical coordinates (m) for the second frequency
+
+        tVect_freq1 = stg.time[data_id][ stg.frequencies_for_calibration[0][1] ]
+        tVect_freq2 = stg.time[data_id][ stg.frequencies_for_calibration[1][1] ]
+
+        rMin = np.min( np.array( [ np.min(rVect_freq1), np.min(rVect_freq2) ] ) ) #Minimum vertical position (m)
+        rMax = np.max( np.array( [ np.max(rVect_freq1), np.max(rVect_freq2) ] ) ) #Maximum vertical position (m)
+
+        rVect_interp = np.linspace( rMin, rMax, rVect_freq1.shape[0] ) #Vector of vertical coordinates (m) used for interpolation
+        tVect_interp = ( tVect_freq1 + tVect_freq2 ) / 2 #Vector of time stamps (s) used for interpolation
+        
+        rGrid_interp = np.zeros( ( rVect_interp.shape[0], rGrid_freq1.shape[1] ) )
+        tGrid_interp = np.zeros( rGrid_interp.shape  )
+
+        for j in range( rGrid_interp.shape[1] ):
+
+            rGrid_interp[:, j] = rVect_interp
+            tGrid_interp[:, j] = tVect_interp[j] * np.ones( rGrid_interp.shape[0] )
+
+
+        stg.r2D_interp[data_id] = rGrid_interp
+        stg.t2D_interp[data_id] = tGrid_interp
+        
+
+
+        # Interpolate on the same grid the quantity 'J' that has been previously
+        # computed for both frequencies:
+
+        j_interpolator_freq1 = sp.interpolate.RegularGridInterpolator( ( rVect_freq1, tVect_interp ),
+                                                                       stg.J_cross_section[data_id][0],
+                                                                       method="linear", fill_value=np.nan )
+
+        j_interpolator_freq2 = sp.interpolate.RegularGridInterpolator( ( rVect_freq2, tVect_interp ),
+                                                                       stg.J_cross_section[data_id][1],
+                                                                       method="linear", fill_value=np.nan )
+
+        j_interp_freq1 = j_interpolator_freq1( ( stg.r2D_interp[data_id], stg.t2D_interp[data_id] ) )
+        j_interp_freq2 = j_interpolator_freq2( ( stg.r2D_interp[data_id], stg.t2D_interp[data_id] ) )
+        
+        stg.J_cross_section_interp[data_id] = [ j_interp_freq1, j_interp_freq2 ]
+        
+        
+        # Compute the Volume Backscatter Index (VBI) on the common grid:
 
         stg.VBI_cross_section[data_id] = np.array([])
 
@@ -537,11 +606,9 @@ class AcousticInversionTab(QWidget):
             freq1=stg.frequencies_for_calibration[0][0],
             freq2=stg.frequencies_for_calibration[1][0],
             zeta_freq1=stg.zeta[0], zeta_freq2=stg.zeta[1],
-            j_cross_section_freq1=stg.J_cross_section[data_id][0],
+            j_cross_section_freq1=stg.J_cross_section_interp[data_id][0],
-            j_cross_section_freq2=stg.J_cross_section[data_id][1],
+            j_cross_section_freq2=stg.J_cross_section_interp[data_id][1],
-            r2D=stg.depth_2D[data_id][
+            r2D=stg.r2D_interp[data_id],
-                stg.frequency_for_inversion[1]
-            ],
             water_attenuation_freq1=stg.alpha_s[0],
             water_attenuation_freq2=stg.alpha_s[1],
             X=stg.X_exponent[0]
@@ -555,39 +622,24 @@ class AcousticInversionTab(QWidget):
         '''
         stg.SSC_fine[data_id] = self.inv_hc.SSC_fine(
             zeta=stg.zeta[0],
-            r2D=stg.depth_2D[data_id][stg.frequencies_for_calibration[0][1]],
+            r2D=stg.r2D_interp[data_id],
             VBI=stg.VBI_cross_section[data_id],
             freq=stg.frequencies_for_calibration[0][0],
             X=stg.X_exponent[0],
             j_cross_section=stg.J_cross_section[data_id][0],
-            alpha_w=np.full(
+            alpha_w=stg.water_attenuation[data_id][ stg.frequencies_for_calibration[0][1] ]
-                shape=stg.depth_2D[data_id][
-                    stg.frequencies_for_calibration[0][1]
-                ].shape,
-                fill_value=stg.water_attenuation[data_id][
-                    stg.frequencies_for_calibration[0][1]
-                ]
-            )
         ) #Inversion using the first frequency
         '''
 
-
         
         stg.SSC_fine[data_id] = self.inv_hc.SSC_fine(
             zeta=stg.zeta[1],
-            r2D=stg.depth_2D[data_id][stg.frequency_for_inversion[1]],
+            r2D=stg.r2D_interp[data_id],
             VBI=stg.VBI_cross_section[data_id],
             freq=stg.frequencies_for_calibration[1][0],
             X=stg.X_exponent[0],
-            j_cross_section=stg.J_cross_section[data_id][1],
+            j_cross_section=stg.J_cross_section_interp[data_id][1],
-            alpha_w=np.full(
+            alpha_w=stg.water_attenuation[data_id][ stg.frequency_for_inversion[1] ]
-                shape=stg.depth_2D[data_id][
-                    stg.frequency_for_inversion[1]
-                ].shape,
-                fill_value=stg.water_attenuation[data_id][
-                    stg.frequency_for_inversion[1]
-                ]
-            )
         ) #Inversion using the second frequency
         
 
@@ -658,9 +710,9 @@ class AcousticInversionTab(QWidget):
                 time_data = stg.time
                 depth_data = stg.depth
 
+
             pcm_SSC_fine = self.axis_SSC_fine.pcolormesh(
-                time_data[data_id][stg.frequency_for_inversion[1]],
+                stg.t2D_interp[data_id], -stg.r2D_interp[data_id],
-                -depth_data[data_id][stg.frequency_for_inversion[1]],
                 stg.SSC_fine[data_id],
                 cmap='rainbow', norm=LogNorm(vmin=1e0, vmax=15),
                 shading='gouraud'
@@ -668,21 +720,15 @@ class AcousticInversionTab(QWidget):
 
             if stg.depth_bottom[data_id].shape != (0,):
                 self.axis_SSC_fine.plot(
-                    time_data[data_id][stg.frequency_for_inversion[1]],
+                    stg.t2D_interp[data_id][0,:],
                     -stg.depth_bottom[data_id],
                     color='black', linewidth=1, linestyle="solid"
                 )
 
             self.pcm_SSC_fine_vertical_line,  = self.axis_SSC_fine.plot(
-                time_data[data_id][stg.frequency_for_inversion[1],
+                stg.t2D_interp[data_id][0, self.slider_fine.value() - 1]
-                                   self.slider_fine.value() - 1]
+                * np.ones( ( stg.r2D_interp[data_id].shape[0], 1 ) ),
-                * np.ones(
+                - stg.r2D_interp[data_id][:,0], linestyle="solid", color='r', linewidth=2
-                    depth_data[data_id][
-                        stg.frequency_for_inversion[1]
-                    ].shape
-                ),
-                -depth_data[data_id][stg.frequency_for_inversion[1]],
-                linestyle="solid", color='r', linewidth=2
             )
 
             # --- Plot samples of fine sediments ---
@@ -789,36 +835,19 @@ class AcousticInversionTab(QWidget):
 
             self.plot_fine , = self.axis_vertical_profile_SSC_fine.plot(
                 stg.SSC_fine[data_id][:, self.slider_fine.value() -1],
-                -depth_data[data_id][stg.frequency_for_inversion[1]],
+                -stg.r2D_interp[data_id][:,0],
                 linestyle="solid", linewidth=1, color="k"
             )
 
             self.pcm_SSC_fine_vertical_line.set_data(
-                time_data[data_id][
+                stg.t2D_interp[data_id][0, self.slider_fine.value() - 1]
-                    stg.frequency_for_inversion[1],
+                * np.ones( ( stg.r2D_interp[data_id].shape[0], 1 ) ),
-                    self.slider_fine.value() - 1
+                -stg.r2D_interp[data_id][:,0] )
-                ]
+
-                * np.ones(
-                    depth_data[data_id][
-                        stg.frequency_for_inversion[1]
-                    ].shape
-                ),
-                -depth_data[data_id][
-                    stg.frequency_for_inversion[1]
-                ]
-            )
             self.figure_SSC_fine.canvas.draw_idle()
 
             self.axis_vertical_profile_SSC_fine.set_ylim(
-                [
+                [ - np.nanmax( stg.r2D_interp[data_id][:,0] ), -np.nanmin( stg.r2D_interp[data_id][:,0] ) ] )
-                    -np.nanmax(
-                        depth_data[data_id][stg.frequency_for_inversion[1]]
-                    ),
-                    -np.nanmin(
-                        depth_data[data_id][stg.frequency_for_inversion[1]]
-                    )
-                ]
-            )
 
             self.axis_vertical_profile_SSC_fine.set_xlim(0, 10)
             self.axis_vertical_profile_SSC_fine.set_xlabel("Inverted Fine SSC (g/L)")
@@ -845,30 +874,19 @@ class AcousticInversionTab(QWidget):
 
         self.axis_vertical_profile_SSC_fine.plot(
             stg.SSC_fine[data_id][:, self.slider_fine.value() - 1],
-            -depth_data[data_id][stg.frequency_for_inversion[1]],
+            -stg.r2D_interp[data_id][:,0],
             linestyle="solid", linewidth=1, color="k"
         )
 
         self.pcm_SSC_fine_vertical_line.set_data(
-            time_data[data_id][stg.frequency_for_inversion[1],
+            stg.t2D_interp[data_id][0, self.slider_fine.value() - 1]
-                               self.slider_fine.value() - 1]
+            * np.ones( ( stg.r2D_interp[data_id].shape[0], 1 ) ),
-            * np.ones(depth_data[data_id][
+            -stg.r2D_interp[data_id][:,0]
-                stg.frequency_for_inversion[1]
-            ].shape),
-            -depth_data[data_id][stg.frequency_for_inversion[1]]
         )
         self.figure_SSC_fine.canvas.draw_idle()
 
         self.axis_vertical_profile_SSC_fine.set_ylim(
-            [
+                [ - np.nanmax( stg.r2D_interp[data_id][:,0] ), -np.nanmin( stg.r2D_interp[data_id][:,0] ) ] )
-                -np.nanmax(
-                    depth_data[data_id][stg.frequency_for_inversion[1]]
-                ),
-                -np.nanmin(
-                    depth_data[data_id][stg.frequency_for_inversion[1]]
-                )
-            ]
-        )
 
         for f in stg.fine_sample_profile:
             time_fine_calibration = (
@@ -1218,8 +1236,7 @@ class AcousticInversionTab(QWidget):
                 depth_data = stg.depth
 
             pcm_SSC_sand = self.axis_SSC_sand.pcolormesh(
-                time_data[data_id][stg.frequency_for_inversion[1]],
+                stg.t2D_interp[data_id], -stg.r2D_interp[data_id],
-                -depth_data[data_id][stg.frequency_for_inversion[1]],
                 stg.SSC_sand[data_id],
                 cmap='rainbow', norm=LogNorm(vmin=1e0, vmax=10),
                 shading='gouraud'
@@ -1227,18 +1244,15 @@ class AcousticInversionTab(QWidget):
 
             if stg.depth_bottom[data_id].shape != (0,):
                 self.axis_SSC_sand.plot(
-                    time_data[data_id][stg.frequency_for_inversion[1]],
+                    stg.t2D_interp[data_id][0,:],
                     -stg.depth_bottom[data_id],
                     color='black', linewidth=1, linestyle="solid"
                 )
 
             self.pcm_SSC_sand_vertical_line, = self.axis_SSC_sand.plot(
-                time_data[data_id][stg.frequency_for_inversion[1],
+                stg.t2D_interp[data_id][0, self.slider_sand.value() - 1]
-                                   self.slider_sand.value() - 1]
+                * np.ones( ( stg.r2D_interp[data_id].shape[0], 1 ) ),
-                * np.ones(
+                -stg.r2D_interp[data_id][:,0],
-                    depth_data[data_id][stg.frequency_for_inversion[1]].shape
-                ),
-                -depth_data[data_id][stg.frequency_for_inversion[1]],
                 linestyle="solid", color='r', linewidth=2
             )
 

settings.py

@@ -260,6 +260,9 @@ lin_reg = []                    # (lin_reg_compute.slope, lin_reg_compute.interc
 
 # Variables names               # Description                                                               # Type
 
+r2D_interp = []                 # Spatial grid on which the BS (through 'J') is interpolated                # List of 2D arrays
+t2D_interp = []                 # Temporal grid on which the BS (through 'J') is inteprolated               # List of 2D arrays
+J_cross_section_interp = []     # List of 'J_cross_section' interpolated on 'r2D_interp'                    # List of list of arrays
 VBI_cross_section = []          # Volume Backscattering Index                                               # List of 2D array
 fine_sample_position = []       # Fine samples indexes position for time and depth [(time_index, depth_index)]  # List of tuples
 sand_sample_position = []       # Sand samples indexes position for time and depth [(time_index, depth_index)]  # List of tuples