Commit b48ddb72 authored by Michael Rudolf's avatar Michael Rudolf
Browse files

Merge branch 'devbranch_MR' into 'master'

Analysis and GUI Update

Closes #9 and #6

See merge request !1
parents a17373d6 0e881072
......@@ -69,14 +69,20 @@ For a more detailed description of the other components see the
[Documentation](https://gitext.gfz-potsdam.de/analab-code/rst-evaluation/blob/master/Documentation/)
and the readme files in the individual folders.
## Acknowledgements
The original scripts for version 0.0.1 (4260ec45) in this repository have been provided by M. Warsitzka @warsitzka_at_ig.cas.cz and contain scripts that have been developed by M. Ritter.
Warsitzka, Michael; Ge, Zhiyuan; Schönebeck, Jan-Michael; Pohlenz, Andre; Kukowski, Nina (2019): Ring-shear test data of foam glass beads used for analogue experiments in the Helmholtz Laboratory for Tectonic Modelling (HelTec) at the GFZ German Research Centre for Geosciences in Potsdam and the Institute of Geosciences, Friedrich Schiller University Jena. GFZ Data Services. http://doi.org/10.5880/GFZ.4.1.2019.002
This software reuses images for the user interface that are part of matplotlib (Hunter, 2007).
## Citation
Please cite this repository as:
Rudolf, Michael; Warsitzka, Michael (2021): RST Evaluation - Scripts for analysing shear experiments from the Schulze RST.pc01 ring shear tester. GFZ Data Services. https://doi.org/10.5880/GFZ.4.1.2021.001
## References
[J. D. Hunter, "Matplotlib: A 2D Graphics Environment", Computing in Science & Engineering, vol. 9, no. 3, pp. 90-95, 2007.](https://doi.org/10.1109/MCSE.2007.55)
\ No newline at end of file
......@@ -8,4 +8,5 @@ tqdm
traitlets
uncertainties
uncertainties.egg
pyqt5
\ No newline at end of file
pyqt5
terminal-animation
\ No newline at end of file
......@@ -12,25 +12,14 @@ GUI based approach for RSTpicking
"""
import os
import tkinter as tk
import matplotlib as mpl
from rstevaluation import mainGUI, utilities
from rstevaluation import mainGUI
def run():
app = mainGUI.RST_pick_GUI(None)
app.title('RSTpick GUI')
image_path = utilities.resource_path(
os.path.join('images', 'rst-evaluation_x256.ico')
)
try:
app.iconbitmap(image_path)
except tk.TclError:
pass
mpl.use('qt5agg')
app.mainloop()
......
""" Advanced data analysis functions for ringshear tester data """
import numpy as np
import scipy.optimize as spopt
import scipy.stats as stats
import scipy.signal as spsig
import scipy.stats as spstats
from uncertainties.core import ufloat
......@@ -25,8 +26,8 @@ def rst_analmut(x, y):
# M = M[np.nonzero(np.isfinite(M))]
# C = C[np.nonzero(np.isfinite(C))]
statsM = stats.norm.fit([d.n for d in M])
statsC = stats.norm.fit([d.n for d in C])
statsM = spstats.norm.fit([d.n for d in M])
statsC = spstats.norm.fit([d.n for d in C])
uncM = ufloat(statsM[0], 2*statsM[1])
uncC = ufloat(statsC[0], 2*statsC[1])
fric_mut = (uncM, uncC)
......@@ -80,7 +81,10 @@ def fitfric(x, a, b):
# %%=====================VELOCITY STEPPING ANALYSIS============================
def vst_analysis(data):
def vst_analysis(data, config):
"""
Simply fits the apparent friction with the current loading velocity.
"""
mu_app = data['shearstress']/data['normalstress'] # Apparent friction
logvel = np.log10(data['velocity'])
......@@ -91,11 +95,75 @@ def vst_analysis(data):
onlyfinite = np.isfinite(logvel)
mu_app = mu_app[onlyfinite]
logvel = logvel[onlyfinite]
# Remove values where logvel < -4
mu_app = mu_app[logvel >= -4]
logvel = logvel[logvel >= -4]
pfit, pcov, x, y = vst_fit(
fitfric, logvel, mu_app, config
)
perr = 2*np.sqrt(np.diag(pcov))
name_fit = (
r'y = ({:.4f}$\pm${:.5f})log10(x) +' +
'\n\t\t' +
r'({:.4f}$\pm${:.5f})').format(pfit[0],
perr[0],
pfit[1],
perr[1])
return pfit, perr, name_fit, logvel, mu_app, x, y
pfit, pcov = spopt.curve_fit(
fitfric,
logvel,
mu_app
# %%=====================VELOCITY STEPPING ANALYSIS============================
def vst_analysis_alternative(data, config):
"""
Alternative method for fitting apparent friction provided by E. Kosari when
the data shows strong stick-slip.
First the stick-slip events are detected with peak detection. Then the
peaks are fit with a log-linear equation.
"""
mu_app = data['shearstress']/data['normalstress'] # Apparent friction
logvel = np.log10(data['velocity'])
peaks = spsig.find_peaks(
mu_app,
prominence=config.getfloat('parameters', 'peak_prominence')
)[0].tolist()
troughs = spsig.find_peaks(
-mu_app,
prominence=config.getfloat('parameters', 'peak_prominence')
)[0].tolist()
# If the first trough is before a peak delete it
while troughs[0] < peaks[0]:
troughs.pop(0)
# Remove peaks if they are not followed by a trough:
while peaks[-1] > troughs[-1]:
peaks.pop(-1)
if len(peaks) != len(troughs):
troughs = [
troughs[int(np.argwhere(p < np.array(troughs))[0])]
for p in peaks
]
# Calculate means for fitting
mu_app_fit = np.array(
[(mu_app[p] + mu_app[t])/2 for p, t in zip(peaks, troughs)]
)
logvel_fit = logvel[peaks]
# Remove nonfinite values for fitting
onlyfinite = np.isfinite(mu_app_fit)
mu_app_fit = mu_app_fit[onlyfinite]
logvel_fit = logvel_fit[onlyfinite]
onlyfinite = np.isfinite(logvel_fit)
mu_app_fit = mu_app_fit[onlyfinite]
logvel_fit = logvel_fit[onlyfinite]
# Remove values where logvel < -4
mu_app_fit = mu_app_fit[logvel_fit >= -4]
logvel_fit = logvel_fit[logvel_fit >= -4]
# Do fitting
pfit, pcov, x, y = vst_fit(
fitfric, logvel_fit, mu_app_fit, config
)
perr = 2*np.sqrt(np.diag(pcov))
name_fit = (
......@@ -105,10 +173,36 @@ def vst_analysis(data):
perr[0],
pfit[1],
perr[1])
return pfit, perr, name_fit
return pfit, perr, name_fit, logvel_fit, mu_app_fit, x, y
def vst_fit(fitfunc, logvel, mu_app, config):
"""
Takes a special subsample per velocity step.
"""
vel_round = np.around(logvel, config.getint('parameters', 'vel_accuracy'))
vel_un = np.unique(vel_round)
# sample_counts = [len(np.argwhere(vel_round == v)) for v in vel_un]
# min_samples = np.min(sample_counts)
x = []
y = []
prct = np.linspace(1, 99, config.getint('parameters', 'fit_percentiles'))
for vu in vel_un:
# selection = np.argwhere(vel_round == vu)[:min_samples]
# x.extend(vel_round[selection])
# y.extend(mu_app[selection])
selection = mu_app[np.argwhere(vel_round == vu)]
for p in prct:
x.append(vu)
y.append(np.percentile(selection, p))
x = np.squeeze(np.array(x))
y = np.squeeze(np.array(y))
pfit, pcov = spopt.curve_fit(fitfunc, x, y)
return pfit, pcov, x, y
# %%======================BOOTSTRAP LINEAR REGRESSION==========================
def fit_bootstrap(datax, datay, function=None,
yerr_systematic=0.0, nsigma=2, nsets=100,
yerr=None):
......
""" Contains some custom classes for the use with the main GUI """
import tkinter as tk
from tkinter import ttk
class ToolTip(object):
"""
ToolTip object which is shown while hovering over an element.
Adapted from:
https://stackoverflow.com/questions/20399243/display-message-when-hovering-over-something-with-mouse-cursor-in-python
Possibly contains code by voidspace.org.uk (offline at the time of writing)
"""
def __init__(self, widget):
self.widget = widget
self.tipwindow = None
self.id = None
self.x = self.y = 0
def showtip(self, text):
"Display text in tooltip window"
self.text = text
if self.tipwindow or not self.text:
return
x, y, cx, cy = self.widget.bbox("insert")
x = x + self.widget.winfo_rootx() + 57
y = y + cy + self.widget.winfo_rooty() + 27
self.tipwindow = tw = tk.Toplevel(self.widget)
tw.wm_overrideredirect(1)
tw.wm_geometry("+%d+%d" % (x, y))
label = ttk.Label(
tw,
text=self.text, justify=tk.LEFT,
background="#ffffe0", relief=tk.SOLID, borderwidth=1,
font=("tahoma", "8", "normal")
)
label.pack(ipadx=1)
def hidetip(self):
tw = self.tipwindow
self.tipwindow = None
if tw:
tw.destroy()
def create_tooltip(widget, text):
toolTip = ToolTip(widget)
text = format_text(text)
def enter(event):
toolTip.showtip(text)
def leave(event):
toolTip.hidetip()
widget.bind('<Enter>', enter)
widget.bind('<Leave>', leave)
def format_text(text, num=80):
"""
Inserts new lines into text after 'num' characters when there is a space
character.
"""
split_text = text.split(' ')
new_text = ''
current_line = ''
for s in split_text:
if len(current_line) + len(s) < num:
current_line += s + ' '
else:
new_text += current_line + '\n'
current_line = s + ' '
if current_line:
new_text += current_line
return new_text
......@@ -9,17 +9,20 @@ import operator
import os
import shutil
import animation
import nptdms
import numpy as np
import uncertainties as unc
from scipy.signal import medfilt
from tqdm import tqdm
from uncertainties import unumpy as unp
from . import data as rstdat
from rstevaluation import data as rstdat
@animation.wait()
def convert(path, file_in, config):
"""Reads data from source files and returns a dictionary with data"""
if isinstance(config, configparser.ConfigParser):
var = {
k: json.loads(config['parameters'][k])
......@@ -66,7 +69,19 @@ def convert(path, file_in, config):
data['displacement'] = np.cumsum(data['time'][1]*data['velocity'])
data['name'] = file_in
data['shear_smth_nom'] = medfilt(
data['shearstress'],
config.getint('parameters', 'smoothwindow')
)
data['shear_smth'] = unp.uarray(
data['shear_smth_nom'],
np.full_like(
data['shear_smth_nom'],
# estimated accuracy of machine (+-10 Pa)
# and error of sensor (0.5%)
np.sqrt(10**2 + (data['shear_smth_nom']*0.005)**2)
)
)
print(file_in+' read')
return data
......
""" Provides PhotoImage objects for all images """
import tkinter as tk
from importlib import resources
def icon_rst():
""" Returns main application icon """
with resources.path('rstevaluation.icons', 'rst-evaluation.png') as fpath:
phim = tk.PhotoImage(
file=fpath
)
return phim
def icon_load():
""" Returns an icon for loading from a folder """
with resources.path('rstevaluation.icons', 'home.png') as fpath:
phim = tk.PhotoImage(
file=fpath
)
return phim
def icon_save():
""" Returns an icon for saving to a folder """
with resources.path('rstevaluation.icons', 'filesave.png') as fpath:
phim = tk.PhotoImage(
file=fpath
)
return phim
def icon_options():
""" Returns an icon for options """
with resources.path('rstevaluation.icons', 'subplots.png') as fpath:
phim = tk.PhotoImage(
file=fpath
)
return phim
def icon_start():
""" Returns an icon for processing """
with resources.path(
'rstevaluation.icons', 'qt4_editor_options.png'
) as fpath:
phim = tk.PhotoImage(
file=fpath
)
return phim
This diff is collapsed.
......@@ -10,8 +10,6 @@ import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import uncertainties as unc
from scipy.signal import medfilt
from uncertainties import unumpy as unp
# %%==================EVALUATION SHEAR CURVE DATA==============================
......@@ -73,21 +71,6 @@ def eval_shearstress(cur_dat, config, review=None):
def _auto_pick(cur_dat, config):
"""Helper function to automatically pick the points"""
# smoothing function: 1st#: window size, 2nd#: polynomial order
shear_smth = medfilt(
cur_dat['shearstress'],
config.getint('parameters', 'smoothwindow')
)
shear_smth = unp.uarray(
shear_smth,
np.full_like(
shear_smth,
# estimated accuracy of machine (+-10 Pa)
# and error of sensor (0.5%)
np.sqrt(10**2 + (shear_smth*0.005)**2)
)
)
# Searching for velocity changes to define data range
vel_above = np.nonzero(
cur_dat['velocity'] > config.getfloat('parameters', 'velnoise')
......@@ -98,16 +81,16 @@ def _auto_pick(cur_dat, config):
start_2 = end_1+np.max(np.diff(vel_above))+10
end_2 = vel_above[0][-1]
# Peak friction
i_peak = np.argmax(shear_smth[start_1:end_1])+start_1
tau_peak = shear_smth[i_peak]
i_peak = np.argmax(cur_dat['shear_smth'][start_1:end_1])+start_1
tau_peak = cur_dat['shear_smth'][i_peak]
# Dynamic friction
i_dyn = np.argmin(shear_smth[i_peak:end_1])+i_peak
tau_dyn = shear_smth[i_dyn]
i_dyn = np.argmin(cur_dat['shear_smth'][i_peak:end_1])+i_peak
tau_dyn = cur_dat['shear_smth'][i_dyn]
# Static friction (reactivation)
i_stat = np.argmax(shear_smth[start_2:end_2])+start_2
tau_stat = shear_smth[i_stat]
i_stat = np.argmax(cur_dat['shear_smth'][start_2:end_2])+start_2
tau_stat = cur_dat['shear_smth'][i_stat]
# Store picks in dictionary
picks = dict()
......@@ -153,7 +136,7 @@ def _manual_pick(cur_dat, config,
while not goodPick:
try:
(sel_dat, sel) = _select_data(cur_dat)
pdata = pick_func[pick](cur_dat['shearstress'][sel_dat])
pdata = pick_func[pick](cur_dat['shear_smth_nom'][sel_dat])
picked_val = unc.ufloat(
pdata, np.sqrt(10**2+pdata*.005**2)
)
......@@ -163,7 +146,7 @@ def _manual_pick(cur_dat, config,
ax.plot(
cur_dat['displacement'][sel_dat],
cur_dat['shearstress'][sel_dat],
cur_dat['shear_smth_nom'][sel_dat],
color=pick_col[pick]
)
plt.pause(0.001)
......@@ -261,12 +244,9 @@ def _pick_base_plot(cur_dat, config):
alpha=.5,
label='2$\sigma$ accuracy'
)
shear_smth = medfilt(
cur_dat['shearstress'], config.getint('parameters', 'smoothwindow')
)
ax.plot(
cur_dat['displacement'],
shear_smth,
cur_dat['shear_smth_nom'],
':',
label='filtered data',
color='C1'
......
""" Plotting functions for rst-evaluation """
import os
import animation
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np
......@@ -8,10 +9,11 @@ import scipy.stats as stats
import uncertainties as unc
from uncertainties.core import ufloat
from . import analysis as rst_analysis
from rstevaluation import analysis as rst_analysis
# ===================LINEAR REGRESSION=========================================
@animation.wait()
def plotstd(path, name, strength, fricdata):
""" Plots the results of the linear regression """
plt.rcParams['savefig.format'] = 'pdf'
......@@ -77,6 +79,7 @@ def plotstd(path, name, strength, fricdata):
# %%======================PLOT HISTOGRAMS======================================
@animation.wait()
def plothist(path, name, strength, data_mut):
""" Plots the result of the mutual linear regression as histograms """
title_mu = [
......@@ -185,6 +188,7 @@ def plothist(path, name, strength, data_mut):
# %%=======================PLOT TS============================================
@animation.wait()
def plotts(path, name, exp_data, normal_stress):
""" Plots all experiments into one plot """
plt.rcParams['figure.figsize'] = (10, 8)
......@@ -238,7 +242,8 @@ def plotts(path, name, exp_data, normal_stress):
# %%=================PLOT VST ANALYSIS=========================================
def plotVST(path, data, pfit, perr, name_fit):
@animation.wait()
def plotVST(path, data, pfit, perr, name_fit, x, y, xq, yq):
""" Plots VST analysis data """
plt.rcParams['savefig.format'] = 'pdf'
plt.rcParams['font.size'] = 10
......@@ -276,9 +281,10 @@ def plotVST(path, data, pfit, perr, name_fit):
ax2.set_xlim(0, np.nanmax(displ) + np.nanmax(displ)/100)
# shear velocity vs. fricton:
ax3.plot(vel, fric, 'k.', ms=2, rasterized=True)
ax3.plot(10**x, y, 'k.', ms=2, rasterized=True, label='original data')
ax3.plot(10**xq, yq, 'r.', ms=2, rasterized=True, label='fit data')
ax3.plot(
vel, rst_analysis.fitfric(np.log10(vel), *pfit),
10**x, rst_analysis.fitfric(x, *pfit),
'g-',
label='curve fit: ' + name_fit
)
......@@ -286,13 +292,13 @@ def plotVST(path, data, pfit, perr, name_fit):
ax3.set_xlabel(r'Shear velocity $v$ [$\mathregular{mm s^{-1}}$]')
ax3.set_ylabel('Friction')
ax3.set_xlim(
np.nanmin(vel)-(np.nanmin(vel)/5),
np.nanmax(vel)+np.nanmax(vel)/5
np.nanmin(10**x)-(np.nanmin(10**x)/5),
np.nanmax(10**x)+np.nanmax(10**x)/5
)
ax3.legend(fontsize=8,
facecolor='w',
edgecolor='k',
loc='upper right',
# loc='upper right',
framealpha=0.5)
fname = os.path.splitext(data['name'])[0]
fig.suptitle(fname, fontsize=14, y=0.95)
......
......@@ -5,49 +5,40 @@ import os
import subprocess
import sys
from . import analysis as rstanalysis
from . import files as rstfiles
from . import picking as rstpicking
from . import plots as rstplots
from tqdm import tqdm
from rstevaluation import analysis as rstanalysis
from rstevaluation import files as rstfiles
from rstevaluation import picking as rstpicking
from rstevaluation import plots as rstplots
def evaluate(config):
""" Evaluates the data coming from the config """
exp_data = get_data(config)
if config.getint('options', 'is_VST'):
if config.getint('options', 'rst'):
eval_data = [
rstpicking.eval_shearstress(
cur_dat,
config,
review='auto'
)
for cur_dat in exp_data
]
evaluate_data(eval_data, exp_data, config)
else:
print('Velocity Stepping Test')
for exp in exp_data:
pfit, perr, name_fit = rstanalysis.vst_analysis(exp)
for exp in tqdm(exp_data, desc='Analysing', leave=False):
if config.getint('options', 'stick_slip'):
pf, pe, nfit, x, y, xq, yq = rstanalysis.vst_analysis_alternative(
exp, config)
else: