Formatting and Dilation Analysis

- Added the possibility to get dilation data
- Replaced Savgol filter with median filter because of pick distortion
- Some minor formatting changes

.gitignore

@@ -7,7 +7,6 @@ Examples/
 # ignore data files
 *.tdms
 *.h5
-*.txt
 *.asc
 
 # ignore pictures and pdfs

RSTpicking/RST_Func.py

@@ -15,26 +15,28 @@ import csv
 import os
 import nptdms
 from scipy import stats
-from scipy.signal import savgol_filter
 from scipy import optimize
+from scipy.signal import medfilt
 import time
 import sys
+import collections
+from operator import itemgetter
 
 
 # %%==============CONVERSION===================================================
-def convert(path, file, var):
+def convert(path, file_in, var):
     """Reads data from source files and returns a dictionary with data"""
     # possible inputs and linked helper functions
     file_convert = {'.asc': _readasc,
                     '.dat': _readdat,
                     '.tdms': _readtdms}
 
-    (root, ext) = os.path.splitext(path+file)
-    data = file_convert[ext](path, file)  # data contains basic content
+    (root, ext) = os.path.splitext(path+file_in)
+    data = file_convert[ext](path, file_in)  # data contains basic content
 
     # Corrections
     data['time'] = data['time']-data['time'][0]
-    data['liddispl'] = -(data['liddispl']-data['liddispl'][0])
+    # data['liddispl'] = -(data['liddispl']-data['liddispl'][0])
     data['velocity'][data['velocity'] < 0] = 0
 
     # Additional data
@@ -42,7 +44,7 @@ def convert(path, file, var):
     data['normalstress'] = data['normalforce']/var['A']
     data['displacement'] = np.cumsum(data['time'][1]*data['velocity'])
 
-    print(file+' read')
+    print(file_in+' read')
     return data
 
 
@@ -54,7 +56,7 @@ def eval_shearstress(cur_dat, var, review=None):
     been chosen by the automatic picking algorithm, or skip autopicking.
 
     Possible review options
-    -----------------------
+    ---
     None (default):
         automatic picking without revision
     'auto':
@@ -68,10 +70,7 @@ def eval_shearstress(cur_dat, var, review=None):
     def _auto_pick(cur_dat, var):
         """Helper function to automatically pick the points"""
         # smoothing function: 1st#: window size, 2nd#: polynomial order
-        shear_smth = savgol_filter(
-            cur_dat['shearstress'],
-            int(float(var['smoothwindow'])),
-            3)
+        shear_smth = medfilt(cur_dat['shearstress'], 11)
         # Processing
         velnoise, stressnoise = var['velnoise'], var['stressnoise']
         # Searching for velocity changes to define data range
@@ -157,6 +156,7 @@ def eval_shearstress(cur_dat, var, review=None):
         return picks
 
     def _evaluate_picks(cur_dat, picks):
+        """ Postprocessing of picked data """
         # Calculate values
         peak = [cur_dat['displacement'][picks['peak'][1]],
                 picks['peak'][0]]
@@ -165,14 +165,25 @@ def eval_shearstress(cur_dat, var, review=None):
         stat = [cur_dat['displacement'][picks['static'][1]],
                 picks['static'][0]]
 
+        # Store pick indices for later analysis
+        p_ind = dict()
+        p_ind['peak'] = picks['peak'][1]
+        p_ind['dynamic'] = picks['dynamic'][1]
+        p_ind['static'] = picks['static'][1]
+
         # Calculate dilation for peak
         i_d1 = np.argmin(cur_dat['liddispl'][:picks['peak'][1]])
         d1 = cur_dat['liddispl'][i_d1]
-        i_d2 = np.argmax(
-            cur_dat['liddispl'][picks['peak'][1]:picks['dynamic'][1]])
+        try:  # Sometimes autodetect finds two peaks at the same spot
+            i_d2 = np.argmax(
+                cur_dat['liddispl'][picks['peak'][1]:picks['dynamic'][1]])
+        except ValueError as e:  # Fallback = peak pick
+            print('Detected two peaks at same position!')
+            print('Please do a manual revision!!!')
+            i_d2 = picks['peak'][1]
         d2 = cur_dat['liddispl'][i_d2]
         deltad = d2-d1
-        return (peak, dyn, stat, deltad)
+        return (peak, dyn, stat, deltad, p_ind)
 
     def _select_data(cur_dat):
         """Creates a slice of +-20 points around a selected point"""
@@ -197,10 +208,7 @@ def eval_shearstress(cur_dat, var, review=None):
         ax.plot(cur_dat['displacement'],
                 cur_dat['shearstress'],
                 label='original data')
-        shear_smth = savgol_filter(
-            cur_dat['shearstress'],
-            int(float(var['smoothwindow'])),
-            3)
+        shear_smth = medfilt(cur_dat['shearstress'], 11)
         ax.plot(cur_dat['displacement'],
                 shear_smth,
                 ':',
@@ -222,7 +230,7 @@ def eval_shearstress(cur_dat, var, review=None):
 
         base_col_list = ['r', 'y', 'b']
         col_list = [base_col_list[i] for i in range(len(picks))]
-        (peak, dyn, stat, deltad) = _evaluate_picks(cur_dat, picks)
+        (peak, dyn, stat, deltad, p_ind) = _evaluate_picks(cur_dat, picks)
         # fig.show()
         pnt_peak = ax.plot(peak[0], peak[1], 's',
                            color='r',
@@ -265,7 +273,7 @@ def eval_shearstress(cur_dat, var, review=None):
         def _cb_peak(event):
             nonlocal pnt_peak
             btn['peak']()
-            (peak, dyn, stat, deltad) = _evaluate_picks(cur_dat, picks)
+            (peak, dyn, stat, deltad, p_ind) = _evaluate_picks(cur_dat, picks)
             for p in pnt_peak:
                 p.set_xdata(peak[0])
                 p.set_ydata(peak[1])
@@ -273,7 +281,7 @@ def eval_shearstress(cur_dat, var, review=None):
         def _cb_dyn(event):
             nonlocal pnt_dyn
             btn['dyn']()
-            (peak, dyn, stat, deltad) = _evaluate_picks(cur_dat, picks)
+            (peak, dyn, stat, deltad, p_ind) = _evaluate_picks(cur_dat, picks)
             for p in pnt_dyn:
                 p.set_xdata(dyn[0])
                 p.set_ydata(dyn[1])
@@ -281,7 +289,7 @@ def eval_shearstress(cur_dat, var, review=None):
         def _cb_stat(event):
             nonlocal pnt_stat
             btn['stat']()
-            (peak, dyn, stat, deltad) = _evaluate_picks(cur_dat, picks)
+            (peak, dyn, stat, deltad, p_ind) = _evaluate_picks(cur_dat, picks)
             for p in pnt_stat:
                 p.set_xdata(stat[0])
                 p.set_ydata(stat[1])
@@ -328,21 +336,23 @@ def eval_shearstress(cur_dat, var, review=None):
     else:
         print('No revision set.')
         pass
-    (peak, dyn, stat, deltad) = _evaluate_picks(cur_dat, picks)
+    (peak, dyn, stat, deltad, p_ind) = _evaluate_picks(cur_dat, picks)
     # Further calculations for output
     normal = int(np.mean(cur_dat['normalstress'])/var['prec'])*var['prec']
     # Calculate relative weakening (after Ritter et al., 2016)
     weak = 1-(dyn[1]/peak[1])
     weak_p = (peak[1]/dyn[1])-1
 
-    return (normal, peak, dyn, stat, deltad, weak, weak_p)
+    return (normal, peak, dyn, stat, deltad, weak, weak_p, p_ind)
 
 
 # %%=======================ANALYSIS OF DATA====================================
 # =======================MUTUAL LINEAR REGRESSION==============================
 def rst_analmut(x, y):
-    # calculation of friction coefficients (slope, M) and y-axis intercept
-    # (cohesion, C) by mutual two point linear regressions:
+    """
+    calculation of friction coefficients (slope, M) and y-axis intercept
+    (cohesion, C) by mutual two point linear regressions:
+    """
     n = len(x)
     ys = [y[k] for k in range(len(x)) for j in range(len(x)-k)]
     xs = [x[k] for k in range(len(x)) for j in range(len(x)-k)]
@@ -363,7 +373,7 @@ def rst_analmut(x, y):
 
 # %%=======================STANDARD LINEAR REGRESSION==========================
 def rst_analstd(x, y):
-    # Correlation of normal and shear stress:  y = slope * x+intercept
+    """ Correlation of normal and shear stress:  y = slope * x+intercept """
     x = np.array(x)
     y = np.array(y)
     n = len(x)
@@ -440,10 +450,11 @@ def fit_bootstrap(p0, datax, datay, function=None,
 # %%========================PLOT===============================================
 # ===================LINEAR REGRESSION=========================================
 def plotstd(path, name, strength, fricdata):
+    """ Plots the results of the linear regression """
     plt.rcParams['savefig.format'] = 'pdf'
     plt.rcParams['font.size'] = 10
     plt.rcParams['font.family'] = 'Arial Unicode MS'
-    plt.rcParams['savefig.dpi'] = 1000
+    plt.rcParams['savefig.dpi'] = 300
     plt.rcParams['figure.figsize'] = (8, 8)
     fig1 = plt.figure()
     linecolor = ['royalblue', 'r', 'g']
@@ -493,6 +504,7 @@ def plotstd(path, name, strength, fricdata):
 
 # %%======================PLOT HISTOGRAMS======================================
 def plothist(path, name, strength, data_mut):
+    """ Plots the result of the mutual linear regression as histograms """
     title_mu = ['Peak friction coefficient $\\mu_P$',
                 'Dynamic friction coefficient $\\mu_D$',
                 'Reactivation friction coefficient $\\mu_R$']
@@ -582,11 +594,20 @@ def plothist(path, name, strength, data_mut):
 
 # %%=======================PLOT  TS============================================
 def plotts(path, name, exp_data, normal_stress):
+    """ Plots all experiments into one plot """
     plt.rcParams['figure.figsize'] = (10, 8)
 
     # =======================PLOT TIME SERIES==================================
     uq_normal = np.sort(np.unique(normal_stress))
-    linecolor = ['red', 'orange', 'gold', 'green', 'royalblue', 'k']
+    linecolor = ['red',
+                 'orange',
+                 'gold',
+                 'green',
+                 'royalblue',
+                 'k',
+                 'yellow',
+                 'teal',
+                 'aqua']
     # Plot everything
     fig3, ax3 = plt.subplots()
     for cur_set, cur_norm in zip(exp_data, normal_stress):
@@ -625,6 +646,7 @@ def plotts(path, name, exp_data, normal_stress):
 
 # %%========================SAVE===============================================
 def saveTS(path, name, exp_data):
+    """ Saves time series in a txt file """
     # create header
     header = ['Time [s]']
     n_stress = ['%.0f' % np.mean(m['normalstress']) for m in exp_data]
@@ -746,28 +768,125 @@ def saveFric(path, name, fricmut, fricstd):
     f.closed
 
 
-# ---- Helper function to read *.asc file ----
-def _readasc(path, file):
-    with codecs.open(path+file, encoding='utf-8-sig') as f:
+def savelidpos(path, name, p_ind, exp_data):
+    """ Saves lid position at peak location to estimate densities """
+    lid_pos = collections.defaultdict(list)
+    for p, ex in zip(p_ind, exp_data):
+        for k in p.keys():
+            lid_pos[k].append(ex['liddispl'][p[k]])
+    names = [ex['name'] for ex in exp_data]
+    header = [k for k in lid_pos.keys()]
+    header.insert(0, 'experiment')
+    data = [lid_pos[k] for k in lid_pos.keys()]
+    rows = [[names[i], data[0][i], data[1][i], data[2][i]]
+            for (i, _) in enumerate(data[0])]
+    rows.sort(key=itemgetter(0))
+    with open(path+name+'_lidpos.txt', 'w+') as f:
+        w = csv.writer(f,
+                       delimiter=',')
+        w.writerow(header)
+        w.writerows(rows)
+    return rows
+
+
+def dilation(path, name, lid_pos, exp_data, prop_file, cfg):
+    """ Analysis of dilation and density during experiment """
+
+    def _getexpnum(exp):
+        prefix = exp['name'].split(' ')[0]
+        if '_' in prefix:
+            number = int(prefix.split('_')[-1])
+        else:
+            number = int(prefix.split('-')[-1])
+        return number
+
+    def _getposnum(lpos):
+        prefix = lpos[0].split(' ')[0]
+        if '_' in prefix:
+            number = int(prefix.split('_')[-1])
+        else:
+            number = int(prefix.split('-')[-1])
+        return number
+
+    props = _readpropfile(prop_file)
+    num_props = np.array([int(d[0].split('-')[-1]) for d in props['data']])
+    num_exp = np.array([_getexpnum(exp) for exp in exp_data])
+    num_pos = np.array([_getposnum(lpos) for lpos in lid_pos])
+
+    ref_mass = cfg.getfloat('parameters', 'cellweight')
+    masses = np.array([float(d[2]) for d in props['data']]) - ref_mass
+    normloads = np.array([float(d[1]) for d in props['data']])
+    tardens = np.array([float(d[3]) for d in props['data']])
+    area = cfg.getfloat('parameters', 'A')
+
+    linecolor = ['royalblue', 'r', 'g', 'k']
+    marker = ['o', '^', 's', 'h']
+    label2 = ['Peak', 'Dynamic', 'Reactivation', 'Initial']
+
+    try:
+        compact_density = [float(d[3]) for d in props['data']]
+    except expression as identifier:
+        pass
+
+    fig, (ax1, ax2) = plt.subplots(ncols=2,
+                                   sharey=True,
+                                   figsize=(10, 5))
+    densities = []
+    for (i, (mass, norm, tden)) in enumerate(zip(masses, normloads, tardens)):
+        n_p = num_props[i]
+        n_x = np.nonzero(num_exp == n_p)[0][0]
+        n_l = np.nonzero(num_pos == n_p)[0][0]
+        exp = exp_data[n_x]
+        ini_lid = exp['liddispl'][0]
+        lid = lid_pos[n_l][1:]
+        lid.append(ini_lid)
+        for (j, l) in enumerate(lid):
+            density = (mass/1000)/(area*(0.04-l/1000))
+            densities.append(density)
+            ax1.plot(norm, density,
+                     marker=marker[j],
+                     color=linecolor[j],
+                     label=label2[j] if i == 0 else '')
+            ax2.plot(tden, density,
+                     marker=marker[j],
+                     color=linecolor[j])
+    all_dens = np.concatenate((densities, tardens))
+    fig.legend()
+    ax1.set_ylabel('Density during shearing (kg/m$^3$)')
+    ax1.set_xlabel('Normal stress (Pa)')
+    ax2.set_xlabel('Target density (kg/m$^3$')
+    # ax2.set_xlim(
+    #     np.min(all_dens),
+    #     np.max(all_dens)
+    # )
+    # ax2.set_ylim(ax2.get_xlim())
+    fig.savefig(path+name+'_sheardensity.png')
+    plt.close()
+
+
+def _readasc(path, file_in):
+    """ Helper function to read *.asc file """
+    with codecs.open(path+file_in, encoding='utf-8-sig') as f:
         data_load = np.loadtxt(f)  # load file for further calculations
 
     # extract data and convert to SI-units
     time = data_load[:, 0]  # in s
-    shearforce = data_load[:, 1] * 9.81  # kg -> N
-    velocity = data_load[:, 2] / 60  # mm/min -> mm/s
+    shearforce = data_load[:, 2] * 9.81  # kg -> N
+    velocity = data_load[:, 4] / 60  # mm/min -> mm/s
     liddispl = data_load[:, 3]  # in mm
-    normalforce = data_load[:, 4] * 9.81  # kg -> N
+    normalforce = data_load[:, 1] * 9.81  # kg -> N
     data = {'time': time,
             'velocity': velocity,
             'normalforce': normalforce,
             'shearforce': shearforce,
-            'liddispl': liddispl}
+            'liddispl': liddispl,
+            'name': file_in.split('.')[0]}
     return data
 
 
-# ---- Helper function to read *.dat file ----
-def _readdat(path, file):
-    with codecs.open(path+file, encoding='utf-8-sig') as f:
+def _readdat(path, file_in):
+    """ Helper function to read *.dat file  """
+    with codecs.open(path+file_in, encoding='utf-8-sig') as f:
         data_load = np.loadtxt(f, delimiter=';', skiprows=3)
     time = data_load[:, 0]  # in s
     normalforce = data_load[:, 1]  # in N
@@ -778,12 +897,13 @@ def _readdat(path, file):
             'velocity': velocity,
             'normalforce': normalforce,
             'shearforce': shearforce,
-            'liddispl': liddispl}
+            'liddispl': liddispl,
+            'name': file_in.split('.')[0]}
     return data
 
 
-# ---- Helper function to read *.tdms file ----
 def _readtdms(path, file_in):
+    """ Helper function to read *.tdms file """
     f = nptdms.TdmsFile(path+file_in)
     df_raw = dict()
     time = f.group_channels('Untitled')[0].time_track()  # in [s]
@@ -798,5 +918,17 @@ def _readtdms(path, file_in):
             'velocity': velocity,
             'normalforce': normalforce,
             'shearforce': shearforce,
-            'liddispl': liddispl}
+            'liddispl': liddispl,
+            'name': file_in.split('.')[0]}
     return data
+
+
+def _readpropfile(prop_file):
+    """ Reads the properties of the experiment from the csv file """
+    props = dict()
+    with open(prop_file, 'r') as f:
+        csv_reader = csv.reader(f)
+        props['names'] = next(csv_reader)
+        props['units'] = next(csv_reader)
+        props['data'] = [row for row in csv_reader]
+    return props

RSTpicking/RST_pick_GUI.py

@@ -216,6 +216,10 @@ class RST_pick_GUI(tk.Tk):
         if event.widget._name == '!button':  # Input path
             self.entry_input.delete(0, 'end')
             self.entry_input.insert(0, f_path + '/')
+            self.entry_output.delete(0, 'end')
+            self.entry_output.insert(0, f_path + '/')
+            self.entry_projname.delete(0, 'end')
+            self.entry_projname.insert(0, f_path.split('/')[-1])
         elif event.widget._name == '!button2':  # Output path
             self.entry_output.delete(0, 'end')
             self.entry_output.insert(0, f_path + '/')
@@ -266,6 +270,7 @@ class RST_pick_GUI(tk.Tk):
         dilation = [ev[4] for ev in eval_data]
         weak = [ev[5] for ev in eval_data]
         weak_p = [ev[6] for ev in eval_data]
+        p_ind = [ev[7] for ev in eval_data]
 
         # ==========ANALYSIS OF FRICTION COEFFICIENTS AND COHESIONS=======
         # ===========Analysis 1: Mutual linear regression=================
@@ -301,12 +306,22 @@ class RST_pick_GUI(tk.Tk):
         if self.rst.get() == 1:
             rfnc.saveStrength(path_out, projectname, strength)
             rfnc.saveFric(path_out, projectname, fric_mut, fric_std)
+            lid_pos = rfnc.savelidpos(path_out, projectname, p_ind, exp_data)
             print('>>> Friction data saved')
+            prop_file = [path_in+f for f in os.listdir(path_in)
+                         if f.endswith('csv')]
+            if prop_file:
+                print('>>> Prop file detected! Starting dilation analysis...')
+                rfnc.dilation(path_out, projectname, lid_pos, exp_data,
+                              prop_file[0], self.cfg)
+            else:
+                print('>>> There is no prop file! Going on as usual.')
 
         if self.save_ts.get() == 1:
             rfnc.saveTS(path_out, projectname, exp_data)
             print('>>> Time series data saved')
 
+        # ======= Exit if 'exit on processing' is selected =====
         if self.quit_onprc.get() == 1:
             sys.exit(0)
 
@@ -390,7 +405,7 @@ class RST_pick_GUI(tk.Tk):
             'prec': 500,        # precision for normal stress in filenames
             'velnoise': 1e-2,
             'stressnoise': 0.025,
-            'smoothwindow': 51  # define window for smooth shear stress curve
+            'cellweight': 2185.3  # weight of shear cell
         }
 
         config['units'] = {
@@ -401,7 +416,7 @@ class RST_pick_GUI(tk.Tk):
             'prec': 'Pa',
             'velnoise': 'mm/min',
             'stressnoise': 'Pa',
-            'smoothwindow': 'samples'
+            'cellweight': 'g'
         }
         self.path_cfg.set('default.ini')
         with open(self.path_cfg.get(), 'w') as f:

requirements.txt

+scipy==1.3.1
+npTDMS==0.15.1
+numpy==1.17.1
+pandas==0.25.1
+matplotlib==3.1.1
+h5py==2.9.0